From 14c1be09f0639c83910d089247e4ffef714da6c6 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Wed, 6 Mar 2024 14:54:20 +0100 Subject: [PATCH 01/15] Create a DPD integrator --- src/Integrator/DPD.cuh | 74 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 src/Integrator/DPD.cuh diff --git a/src/Integrator/DPD.cuh b/src/Integrator/DPD.cuh new file mode 100644 index 00000000..2d152bad --- /dev/null +++ b/src/Integrator/DPD.cuh @@ -0,0 +1,74 @@ +#pragma once + +#include "Integrator/VerletNVE.cuh" +#include "Interactor/PairForces.cuh" +#include "Interactor/Potential/DPD.cuh" + +namespace uammd{ + class DPDIntegrator: public Integrator{ + std::shared_ptr verlet; + int steps; + public: + struct Parameters{ + //VerletNVE + real energy = 0; //Target energy, ignored if initVelocities is false + real dt = 0; + bool is2D = false; + bool initVelocities = false; + + //DPD + real cutOff; + real gamma; + real A; + real temperature; + real3 L; + }; + + DPDIntegrator(shared_ptr pg, Parameters par): + Integrator(pg, "DPDIntgrator"), + steps(0){ + //Initialize verletNVE + VerletNVE::Parameters verletpar; + verletpar.energy = par.energy; + verletpar.dt = par.dt; + verletpar.is2D = par.is2D; + verletpar.initVelocities = par.initVelocities; //=false? + + verlet = std::make_shared(pg, verletpar); + + //Initialize DPD and add to interactor list. + Potential::DPD::Parameters dpdPar; + dpdPar.temperature = par.temperature; + dpdPar.cutOff = par.cutOff; + dpdPar.gamma = par.gamma; + dpdPar.A = par.A; + dpdPar.dt = par.dt; + auto dpd = std::make_shared(dpdPar); + using PF = PairForces; + typename PF::Parameters pfpar; + pfpar.box = Box(par.L); + //From the example in PairForces + auto dpd_interactor = std::make_shared(pg, pfpar, dpd); + verlet->addInteractor(dpd_interactor); + } + + DPDIntegrator(shared_ptr pd, Parameters par): + DPDIntegrator(std::make_shared(pd, "All"), par){} + + ~DPDIntegrator(){} + + virtual void forwardTime() override { + steps++; + if (steps == 1) + for(auto forceComp: interactors){ + verlet->addInteractor(forceComp); + } + verlet->forwardTime(); + } + + virtual real sumEnergy() override { + return verlet->sumEnergy(); + } + + }; +} From 2de0490e480a22f609778e63218a3059a72db162 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Wed, 6 Mar 2024 14:54:51 +0100 Subject: [PATCH 02/15] Replace ForceTransverser by Transverser --- src/Interactor/Potential/DPD.cuh | 52 +++++++++++++++++--------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/src/Interactor/Potential/DPD.cuh b/src/Interactor/Potential/DPD.cuh index 6a8f544c..90af42ee 100644 --- a/src/Interactor/Potential/DPD.cuh +++ b/src/Interactor/Potential/DPD.cuh @@ -84,7 +84,7 @@ namespace uammd{ } - struct ForceTransverser{ + struct Transverser{ private: real4* pos; real3 * vel; @@ -100,19 +100,19 @@ namespace uammd{ real A; public: - ForceTransverser(real4* pos, real3 *vel, real4* force, - ullint seed, ullint step, - Box box, - int N, - real rcut, - DissipativeStrength gamma, real sigma, real A): + Transverser(real4* pos, real3 *vel, real4* force, + ullint seed, ullint step, + Box box, + int N, + real rcut, + DissipativeStrength gamma, real sigma, real A): pos(pos), vel(vel), force(force), step(step), seed(seed), box(box), N(N), invrcut(1.0/rcut), gamma(gamma), sigma(sigma), A(A){} - using returnInfo = real3; + using returnInfo = ForceEnergyVirial; struct Info{ real3 vel; @@ -130,41 +130,43 @@ namespace uammd{ Saru rng(ij, seed, step); const real rmod = sqrt(dot(rij,rij)); //There is an indetermination at r=0 - if(rmod == real(0)) return make_real3(0); + if(rmod == real(0)) return {}; const real invrmod = real(1.0)/rmod; //The force is 0 beyond rcut - if(invrmod<=invrcut) return make_real3(0); + if(invrmod<=invrcut) return {}; const real wr = real(1.0) - rmod*invrcut; //This weight function is arbitrary as long as wd = wr*wr const real Fc = A*wr*invrmod; const real wd = wr*wr; //Wd must be such as wd = wr^2 to ensure fluctuation dissipation balance const real g = gamma.dissipativeStrength(i, j, pi, pj, infoi.vel, infoj.vel); const real Fd = -g*wd*invrmod*invrmod*dot(rij, vij); const real Fr = rng.gf(real(0.0), sigma*sqrt(g)*wr*invrmod).x; - return (Fc+Fd+Fr)*rij; + returnInfo fev; + fev.force = (Fc+Fd+Fr)*rij; + return fev; } inline __device__ Info getInfo(int pi){return {vel[pi], pi};} - inline __device__ void set(uint pi, const returnInfo &total){ force[pi] += make_real4(total);} + inline __device__ void set(uint pi, const returnInfo &total){ force[pi] += make_real4(total.force);} }; - ForceTransverser getForceTransverser(Box box, shared_ptr pd){ - auto pos = pd->getPos(access::location::gpu, access::mode::read); - auto vel = pd->getVel(access::location::gpu, access::mode::read); - auto force = pd->getForce(access::location::gpu, access::mode::readwrite); - static auto seed = pd->getSystem()->rng().next(); - step++; - int N = pd->getNumParticles(); - return ForceTransverser(pos.raw(), vel.raw(), force.raw(), seed, step, box, N, rcut, gamma, sigma, A); - } - //Notice that no getEnergyTransverser is present, this is not a problem as modules using this potential will fall back to a BasicNullTransverser when the method getEnergyTransverser is not found and the energy will not be computed altogether. - - BasicNullTransverser getForceEnergyTransverser(Box box, shared_ptr pd){ + Transverser getTransverser(Interactor::Computables comp, Box box, std::shared_ptr pd){ + if (comp.energy){ System::log("[DPD] No way of measuring energy in DPD"); - return BasicNullTransverser(); } + if (comp.virial){ + System::log("[DPD] No way of measuring virial in DPD"); + } + auto pos = pd->getPos(access::location::gpu, access::mode::read).raw(); + auto vel = pd->getVel(access::location::gpu, access::mode::read).raw(); + auto force = pd->getForce(access::location::gpu, access::mode::readwrite).raw(); + static auto seed = pd->getSystem()->rng().next(); + step++; + int N = pd->getNumParticles(); + return Transverser(pos, vel, force, seed, step, box, N, rcut, gamma, sigma, A); + } }; template<> From d073192d7b39eaf7e52a973bd449eeec7a602e3b Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Wed, 6 Mar 2024 14:57:01 +0100 Subject: [PATCH 03/15] Add DPD integrator --- src/Integrator/DPD.cuh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Integrator/DPD.cuh b/src/Integrator/DPD.cuh index 2d152bad..ff814cb2 100644 --- a/src/Integrator/DPD.cuh +++ b/src/Integrator/DPD.cuh @@ -59,10 +59,11 @@ namespace uammd{ virtual void forwardTime() override { steps++; - if (steps == 1) + if (steps == 1){ for(auto forceComp: interactors){ verlet->addInteractor(forceComp); } + } verlet->forwardTime(); } From 69ebe009e03eeb437814913fcc0f19cef7113241 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Wed, 6 Mar 2024 18:25:06 +0100 Subject: [PATCH 04/15] Add test for DPDinteractor --- test/DPD/testDPDIntegrator.cu | 158 ++++++++++++++++++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 test/DPD/testDPDIntegrator.cu diff --git a/test/DPD/testDPDIntegrator.cu b/test/DPD/testDPDIntegrator.cu new file mode 100644 index 00000000..26e0674c --- /dev/null +++ b/test/DPD/testDPDIntegrator.cu @@ -0,0 +1,158 @@ +/*Raul P. Pelaez 2022. Tests for the Lanczos algorithm. + Tests the result of sqrt(M)*v for increasingly complex matrices and vectors of several sizes. + */ +#include +#include +#include "gmock/gmock.h" +#include "utils/ParticleSorter.cuh" +#include "utils/container.h" +#include "utils/execution_policy.cuh" +#include +#include +#include "uammd.cuh" +#include "Integrator/DPD.cuh" +#include "Interactor/Potential/Potential.cuh" +#include "utils/InitialConditions.cuh" +#include "Interactor/NeighbourList/VerletList.cuh" +#include "Interactor/NeighbourList/CellList.cuh" + + +using namespace uammd; +using std::make_shared; + +struct Parameters{ + real3 boxSize; + real dt; + int numberParticles; + int numberSteps, printSteps; + real temperature, friction; + real strength; + real cutOff; + real radius; +}; + +//Let us group the UAMMD simulation in this struct that we can pass around +struct UAMMD{ + std::shared_ptr sys; + std::shared_ptr pd; + Parameters par; +}; + +//Creates the DPD integrator +auto createIntegratorDPD(UAMMD sim){ + DPDIntegrator::Parameters par; + par.temperature = sim.par.temperature; + par.cutOff = sim.par.cutOff; + par.gamma = sim.par.friction; + par.A = sim.par.strength; + par.dt = sim.par.dt; + par.initVelocities=false; + par.L = sim.par.boxSize; + auto dpd = make_shared(sim.pd, par); + return dpd; +} + +//Creates a VerletNVE integrator and adds the DPD potential +auto createIntegratorVerletNVEWithDPD(UAMMD sim){ + Potential::DPD::Parameters par; + par.temperature = sim.par.temperature; + par.cutOff = sim.par.cutOff; + par.gamma = sim.par.friction; + par.A = sim.par.strength; + par.dt = sim.par.dt; + auto dpd = make_shared(par); + + using PFDPD = PairForces; + typename PFDPD::Parameters paramsdpd; + paramsdpd.box = Box(sim.par.boxSize); + auto pairforces = make_shared(sim.pd, paramsdpd, dpd); + + using NVE = VerletNVE; + NVE::Parameters params; + params.dt = par.dt; + params.initVelocities=false; + auto verlet = make_shared(sim.pd, params); + verlet->addInteractor(pairforces); + return verlet; +} + +//Set the initial positons for the particles in a fcc lattice +void initializeParticles(UAMMD sim){ + auto pos = sim.pd->getPos(access::cpu, access::write); + auto initial = initLattice(sim.par.boxSize, sim.par.numberParticles, fcc); + std::copy(initial.begin(), initial.end(), pos.begin()); +} + +//Creates the LJ interactor +auto createLJInteractor(UAMMD sim){ + auto pot = make_shared(); + Potential::LJ::InputPairParameters par; + par.epsilon = 1.0; + par.shift = false; + par.sigma = 2*sim.par.radius; + par.cutOff = 2.5*par.sigma; + pot->setPotParameters(0, 0, par); + + using PairForces = PairForces; + PairForces::Parameters params; + params.box = Box(sim.par.boxSize); //Box to work on + auto pairforceslj = make_shared(sim.pd, params, pot); + return pairforceslj; +} + +bool compareReal3(real3 v1, real3 v2){ + return v1.x == v2.x and v1.y == v2.y and v1.z == v2.z; +} + +bool checkSamePositions(UAMMD sim1, UAMMD sim2){ + auto pos1 = sim1.pd->getPos(access::cpu, access::read); + auto pos2 = sim2.pd->getPos(access::cpu, access::read); + bool samePositions = true; + fori(0, pos1.size()){ + samePositions = compareReal3(make_real3(pos1[i]), make_real3(pos2[i])); + } + return samePositions; +} + +TEST(DPD, PositionsConsistentAcrossMethods){ + auto sys = make_shared(); + Parameters par; + par.temperature = 1.25; + par.radius = 1; + par.cutOff = 2*2.5*par.radius; + par.friction = 1.3; + par.strength = 1.5; + par.dt = 0.001; + par.boxSize = make_real3(25,25,25); + par.numberParticles = 1000; + par.numberSteps = 100; + + + //Prepare the simulation that runs the VerletNVE integrator with DPD potential + UAMMD sim1; + sim1.sys = sys; + sim1.par = par; + sim1.pd = make_shared(par.numberParticles, sim1.sys); + initializeParticles(sim1); + auto integrator1 = createIntegratorVerletNVEWithDPD(sim1); + auto lj1 = createLJInteractor(sim1); + integrator1->addInteractor(lj1); + + //Prepare the simulation that runs the DPD integrator + UAMMD sim2; + sim2.sys = sys; + sim2.par = par; + sim2.pd = make_shared(sim2.par.numberParticles, sim2.sys); + initializeParticles(sim2); + auto integrator2 = createIntegratorDPD(sim2); + auto lj2 = createLJInteractor(sim2); + integrator2->addInteractor(lj2); + + //Run both simulations and check they move the particles in the same way + fori(0, par.numberSteps){ + integrator1->forwardTime(); + integrator2->forwardTime(); + ASSERT_TRUE(checkSamePositions(sim1, sim2)); + } + sys->finish(); +} From 2a7082395a4a7bec14c06a5e0f370e754c700fae Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Wed, 6 Mar 2024 18:25:43 +0100 Subject: [PATCH 05/15] Update CMakeLists.txt --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f61d5d26..b8ce90a5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -72,6 +72,7 @@ add_subdirectory(utils) add_subdirectory(misc/chebyshev) add_subdirectory(misc/ibm) add_subdirectory(misc/lanczos) +add_subdirectory(DPD) IF (CMAKE_BUILD_TYPE MATCHES "Debug") if(CMAKE_CUDA_COMPILER_ID STREQUAL "NVIDIA") From d427159e7f29c12acacdb96ce2bb68851f56749d Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Wed, 6 Mar 2024 18:26:31 +0100 Subject: [PATCH 06/15] Add authorship and brief description to Integrator/DPD.cuh --- src/Integrator/DPD.cuh | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Integrator/DPD.cuh b/src/Integrator/DPD.cuh index ff814cb2..7534980e 100644 --- a/src/Integrator/DPD.cuh +++ b/src/Integrator/DPD.cuh @@ -1,3 +1,14 @@ +/*Pablo Palacios-Alonso 2024. DPD integrator module + + This module simulates particle dynamics using a Dissipative Particle Dynamics (DPD) simulation. + To achieve this, it employs the Verlet NVE algorithm to update the positions of the particles + in conjunction with an interactor that treats the typical DPD forces [1] as a UAMMD potential. + + + References: + [1] On the numerical treatment of dissipative particle dynamics and related systems. Leimkuhler and Shang 2015. https://doi.org/10.1016/j.jcp.2014.09.008 +*/ + #pragma once #include "Integrator/VerletNVE.cuh" From f394f0426f33564315c9586dadfa2c0b3c5ca0d6 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Thu, 7 Mar 2024 10:47:28 +0100 Subject: [PATCH 07/15] Ad CMakeLists.txt to DPD --- test/DPD/CMakeLists.txt | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 test/DPD/CMakeLists.txt diff --git a/test/DPD/CMakeLists.txt b/test/DPD/CMakeLists.txt new file mode 100644 index 00000000..5038a965 --- /dev/null +++ b/test/DPD/CMakeLists.txt @@ -0,0 +1,11 @@ +add_executable(testDPDIntegrator testDPDIntegrator.cu) +target_include_directories(testDPDIntegrator PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(testDPDIntegrator cufft cublas) +set_target_properties(testDPDIntegrator PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin") +target_link_libraries( + testDPDIntegrator + GTest::gtest_main + GTest::gmock_main +) +include(GoogleTest) +gtest_discover_tests(testDPDIntegrator) From ed45f8bfaac042c6de6c8fefdc0146f0d9720acc Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Thu, 7 Mar 2024 11:20:19 +0100 Subject: [PATCH 08/15] Update examples to use the DPDIntegrator module --- examples/generic_md/generic_simulation.cu | 29 +++++++------------ examples/integration_schemes/integrators.cu | 32 +++++++-------------- 2 files changed, 20 insertions(+), 41 deletions(-) diff --git a/examples/generic_md/generic_simulation.cu b/examples/generic_md/generic_simulation.cu index b7bfb7df..fc78babf 100644 --- a/examples/generic_md/generic_simulation.cu +++ b/examples/generic_md/generic_simulation.cu @@ -70,7 +70,7 @@ Which has a lot of information. From basic functionality to descriptions and ref #include "Integrator/Integrator.cuh" #include "Integrator/VerletNVE.cuh" #include "Integrator/VerletNVT.cuh" -#include "Interactor/Potential/DPD.cuh" +#include "Integrator/DPD.cuh" #include"Interactor/SPH.cuh" #include"utils/InputFile.h" #include @@ -238,26 +238,17 @@ auto createIntegratorVerletNVE(UAMMD sim){ } //Dissipative Particle Dynamics -//DPD is handled by UAMMD as a VerletNVE integrator with a special short range interaction auto createIntegratorDPD(UAMMD sim){ - using NVE = VerletNVE; - NVE::Parameters par; + DPDIntegrator::Parameters par; par.dt = sim.par.dt; - par.initVelocities = false; - auto verlet = std::make_shared(sim.pd, par); - using DPD = PairForces; - Potential::DPD::Parameters dpd_params; - dpd_params.cutOff = sim.par.cutOff_dpd; - dpd_params.temperature = sim.par.temperature; - dpd_params.gamma = sim.par.gamma_dpd; - dpd_params.A = sim.par.A_dpd; - dpd_params.dt = par.dt; - auto pot = std::make_shared(dpd_params); - DPD::Parameters params; - params.box = Box(sim.par.L); - auto pairforces = std::make_shared(sim.pd, params, pot); - verlet->addInteractor(pairforces); - return verlet; + par.cutOff = sim.par.cutOff_dpd; + par.temperature = sim.par.temperature; + par.gamma = sim.par.gamma_dpd; + par.A = sim.par.A_dpd; + par.dt = par.dt; + par.box = Box(sim.par.L); + auto dpd = std::make_shared(sim.pd, par); + return dpd; } //Smoothed Particle Hydrodynamics diff --git a/examples/integration_schemes/integrators.cu b/examples/integration_schemes/integrators.cu index 1a479fa7..9d646cd9 100644 --- a/examples/integration_schemes/integrators.cu +++ b/examples/integration_schemes/integrators.cu @@ -78,30 +78,18 @@ std::shared_ptr createIntegratorVerletNVE(UAMMD sim){ return std::make_shared(sim.pd, par); } -#include "Integrator/VerletNVE.cuh" -#include "Interactor/Potential/DPD.cuh" -#include"Interactor/PairForces.cuh" -//DPD is handled by UAMMD as a VerletNVE integrator with a special short range interaction +#include "Integrator/DPD.cuh" std::shared_ptr createIntegratorDPD(UAMMD sim){ - using NVE = VerletNVE; - NVE::Parameters par; - par.dt = 1.0; - par.initVelocities = false; - auto verlet = std::make_shared(sim.pd, par); - using DPD = PairForces; - Potential::DPD::Parameters dpd_params; - dpd_params.cutOff = 1.0; - dpd_params.temperature = 1.0; - dpd_params.gamma = 1.0; - dpd_params.A = 1.0; - dpd_params.dt = 0.1; - auto pot = std::make_shared(dpd_params); - DPD::Parameters params; + DPDIntegrator::Parameters par; + par.cutOff = 1.0; + par.temperature = 1.0; + par.gamma = 1.0; + par.A = 1.0; + par.dt = 0.1; real3 L = make_real3(32,32,32); - params.box = Box(L); - auto pairforces = std::make_shared(sim.pd, params, pot); - verlet->addInteractor(pairforces); - return verlet; + par.box = Box(L); + auto dpd = std::make_shared(sim.pd, par); + return dpd; } #include "Integrator/VerletNVE.cuh" From 4959123592733de351e3dadce9a2f749caf18f31 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Thu, 7 Mar 2024 11:21:00 +0100 Subject: [PATCH 09/15] Use box as parameter instead of boxSize --- src/Integrator/DPD.cuh | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/Integrator/DPD.cuh b/src/Integrator/DPD.cuh index 7534980e..eafcab9c 100644 --- a/src/Integrator/DPD.cuh +++ b/src/Integrator/DPD.cuh @@ -22,17 +22,15 @@ namespace uammd{ public: struct Parameters{ //VerletNVE - real energy = 0; //Target energy, ignored if initVelocities is false real dt = 0; bool is2D = false; - bool initVelocities = false; - + //DPD real cutOff; real gamma; real A; real temperature; - real3 L; + Box box; }; DPDIntegrator(shared_ptr pg, Parameters par): @@ -40,10 +38,9 @@ namespace uammd{ steps(0){ //Initialize verletNVE VerletNVE::Parameters verletpar; - verletpar.energy = par.energy; verletpar.dt = par.dt; verletpar.is2D = par.is2D; - verletpar.initVelocities = par.initVelocities; //=false? + verletpar.initVelocities = false; verlet = std::make_shared(pg, verletpar); @@ -57,7 +54,7 @@ namespace uammd{ auto dpd = std::make_shared(dpdPar); using PF = PairForces; typename PF::Parameters pfpar; - pfpar.box = Box(par.L); + pfpar.box = par.box; //From the example in PairForces auto dpd_interactor = std::make_shared(pg, pfpar, dpd); verlet->addInteractor(dpd_interactor); From 46e26aee6f9d00b4171a7975ec711c2270e2f1a6 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Thu, 7 Mar 2024 11:27:40 +0100 Subject: [PATCH 10/15] Small change in test --- test/DPD/testDPDIntegrator.cu | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/DPD/testDPDIntegrator.cu b/test/DPD/testDPDIntegrator.cu index 26e0674c..def7fdbb 100644 --- a/test/DPD/testDPDIntegrator.cu +++ b/test/DPD/testDPDIntegrator.cu @@ -1,5 +1,7 @@ -/*Raul P. Pelaez 2022. Tests for the Lanczos algorithm. - Tests the result of sqrt(M)*v for increasingly complex matrices and vectors of several sizes. +/* P. Palacios-Alonso 2024 + Test of the implementation of DPDIntegrator, it will simulate a LJ gas + using the new DPDIntegrator and a VerletNVE integrator with the DPD potential and check + that the reuslts are the same */ #include #include @@ -46,8 +48,7 @@ auto createIntegratorDPD(UAMMD sim){ par.gamma = sim.par.friction; par.A = sim.par.strength; par.dt = sim.par.dt; - par.initVelocities=false; - par.L = sim.par.boxSize; + par.box = Box(sim.par.boxSize); auto dpd = make_shared(sim.pd, par); return dpd; } From 09957f053138323fd13da80a7d8fdaf911a99593 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Mon, 11 Mar 2024 12:47:33 +0100 Subject: [PATCH 11/15] Test ci.yml --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e4ec6d11..071cf63a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -111,6 +111,7 @@ jobs: make install echo "#include" > test.cu echo "int main(){return 0;}" >> test.cu + echo "${CONDA_PREFIX}" nvcc -std=c++14 -I"${CONDA_PREFIX}/include/uammd" -I"${CONDA_PREFIX}/include/uammd/third_party" test.cu -o /tmp/test From e5cde4b5802bb4546051de99531fb1123ce6e913 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Mon, 11 Mar 2024 12:51:18 +0100 Subject: [PATCH 12/15] Test ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 071cf63a..83726369 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -111,7 +111,7 @@ jobs: make install echo "#include" > test.cu echo "int main(){return 0;}" >> test.cu - echo "${CONDA_PREFIX}" + run: echo "${CONDA_PREFIX}" nvcc -std=c++14 -I"${CONDA_PREFIX}/include/uammd" -I"${CONDA_PREFIX}/include/uammd/third_party" test.cu -o /tmp/test From 694c136842fed40e11e3eb270e756d6e9296a7a7 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Mon, 11 Mar 2024 12:53:29 +0100 Subject: [PATCH 13/15] Test ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 83726369..071cf63a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -111,7 +111,7 @@ jobs: make install echo "#include" > test.cu echo "int main(){return 0;}" >> test.cu - run: echo "${CONDA_PREFIX}" + echo "${CONDA_PREFIX}" nvcc -std=c++14 -I"${CONDA_PREFIX}/include/uammd" -I"${CONDA_PREFIX}/include/uammd/third_party" test.cu -o /tmp/test From 8b3ffe9f43cfd54c5ea6d03a43cfa649ca8ca90a Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Mon, 11 Mar 2024 12:55:29 +0100 Subject: [PATCH 14/15] Test ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 071cf63a..694cc301 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -111,7 +111,7 @@ jobs: make install echo "#include" > test.cu echo "int main(){return 0;}" >> test.cu - echo "${CONDA_PREFIX}" + echo "${CONDA_PREFIX}" nvcc -std=c++14 -I"${CONDA_PREFIX}/include/uammd" -I"${CONDA_PREFIX}/include/uammd/third_party" test.cu -o /tmp/test From e10d165a48a8ef017ef47135e7ce2b8cd4760a94 Mon Sep 17 00:00:00 2001 From: PabloPalaciosAlonso Date: Mon, 11 Mar 2024 15:46:30 +0100 Subject: [PATCH 15/15] Create a SPHIntegrator that combines the VerletNVE integrator with the SPH interactor --- docs/Integrator/LangevinDynamics.rst | 33 ++-- examples/generic_md/generic_simulation.cu | 37 ++-- examples/integration_schemes/integrators.cu | 27 ++- .../integration_schemes/others/SPH_test.cu | 32 ++-- src/Integrator/SPH.cuh | 84 +++++++++ test/CMakeLists.txt | 3 +- test/{ => Languevin}/DPD/CMakeLists.txt | 0 test/{ => Languevin}/DPD/testDPDIntegrator.cu | 0 test/Languevin/SPH/CMakeLists.txt | 11 ++ test/Languevin/SPH/testSPHIntegrator.cu | 159 ++++++++++++++++++ 10 files changed, 315 insertions(+), 71 deletions(-) create mode 100644 src/Integrator/SPH.cuh rename test/{ => Languevin}/DPD/CMakeLists.txt (100%) rename test/{ => Languevin}/DPD/testDPDIntegrator.cu (100%) create mode 100644 test/Languevin/SPH/CMakeLists.txt create mode 100644 test/Languevin/SPH/testSPHIntegrator.cu diff --git a/docs/Integrator/LangevinDynamics.rst b/docs/Integrator/LangevinDynamics.rst index 2da20bf7..b5f8930a 100644 --- a/docs/Integrator/LangevinDynamics.rst +++ b/docs/Integrator/LangevinDynamics.rst @@ -284,40 +284,41 @@ where :math:`\nu` is a provided viscosity, and :math:`\vec{v}_{ij} = \vec{v}_j - Usage ------------- -We couple an SPH :ref:`Interactor` with a :ref:`VerletNVE` :ref:`Integrator`. +For simplicity, an :ref:`Integrator` is defined that couples the SPH :ref:Interactor with a :ref:`VerletNVE` :ref:`Integrator`. -The following parameters are available for the SPH :ref:`Interactor`: +The following parameters are available for the SPH :ref:`Integrator`, they are a combination of the parameters of the SPH :ref:`Interactor` and VerletNVE :ref:`Integrator`: * :cpp:`real support`. The length scale :math:`h`. * :cpp:`real viscosity`. The prefactor for the artificial viscosity, :math:`\mu`. * :cpp:`real gasStiffness`. The prefactor for the ideal gas equation of state, :math:`K`. * :cpp:`real restDensity`. The rest density in the equation of state, :math:`\rho_0`. * :cpp:`Box box`. A :cpp:any:`Box` with the simulation domain information. - + * :cpp:`real dt` The time step. + * :cpp:`bool initVelocities=true` Modify starting velocities to ensure the target temperature from the start. When :cpp:`false` the velocities of the particles are left untouched at initialization. The default is true and sets particle velocities following the botzmann distribution. + * :code:`bool is2D = false` Set to true if the system is 2D + * :cpp:`real energy` Target energy per particle, can be omitted if :cpp:`initVelocities=false`. + * :cpp:`real mass = -1` Mass of all the particles. If >0 all particles will have this mass, otherwise the mass for each particle in :ref:`ParticleData` will be used. If masses have not been set in :ref:`ParticleData` the default mass is 1 for all particles. .. code:: c++ #include "Integrator/VerletNVE.cuh" #include "Interactor/SPH.cuh" using namespace uammd; - //SPH is handled by UAMMD as a VerletNVE integrator with a special interaction + //SPHIntegrator combines a VerletNVE integrator with a special interaction encoding SPH std::shared_ptr createIntegratorSPH(std::shared_ptr pd){ - using NVE = VerletNVE; - NVE::Parameters par; + SPHIntegrator::Parameters par; par.dt = 0.1; par.initVelocities = false; - auto verlet = std::make_shared(pd, par); - SPH::Parameters params; + real3 L = make_real3(32,32,32); - params.box = Box(L); + par.box = Box(L); //Pressure for a given particle "i" in SPH will be computed as gasStiffness·(density_i - restDensity) //Where density is computed as a function of the masses of the surroinding particles //Particle mass starts as 1, but you can change this in customizations.cuh - params.support = 2.4; //Cut off distance for the SPH kernel - params.viscosity = 1.0; //Environment viscosity - params.gasStiffness = 1.0; - params.restDensity = 1.0; - auto sph = std::make_shared(pd, params); - verlet->addInteractor(sph); - return verlet; + par.support = 2.4; //Cut off distance for the SPH kernel + par.viscosity = 1.0; //Environment viscosity + par.gasStiffness = 1.0; + par.restDensity = 1.0; + auto sph_verlet = std::make_shared(pd, par); + return sph_verlet; } diff --git a/examples/generic_md/generic_simulation.cu b/examples/generic_md/generic_simulation.cu index fc78babf..fe1df4c9 100644 --- a/examples/generic_md/generic_simulation.cu +++ b/examples/generic_md/generic_simulation.cu @@ -61,17 +61,17 @@ Which has a lot of information. From basic functionality to descriptions and ref #include"Interactor/AngularBondedForces.cuh" #include"Interactor/TorsionalBondedForces.cuh" #include"Integrator/BrownianDynamics.cuh" -#include "Integrator/BDHI/FIB.cuh" -#include "Integrator/Hydro/ICM.cuh" -#include "Integrator/BDHI/BDHI_EulerMaruyama.cuh" +#include"Integrator/BDHI/FIB.cuh" +#include"Integrator/Hydro/ICM.cuh" +#include"Integrator/BDHI/BDHI_EulerMaruyama.cuh" #include"Integrator/BDHI/BDHI_PSE.cuh" #include"Integrator/BDHI/BDHI_FCM.cuh" #include"Interactor/SpectralEwaldPoisson.cuh" -#include "Integrator/Integrator.cuh" -#include "Integrator/VerletNVE.cuh" -#include "Integrator/VerletNVT.cuh" -#include "Integrator/DPD.cuh" -#include"Interactor/SPH.cuh" +#include"Integrator/Integrator.cuh" +#include"Integrator/VerletNVE.cuh" +#include"Integrator/VerletNVT.cuh" +#include"Integrator/DPD.cuh" +#include"Integrator/SPH.cuh" #include"utils/InputFile.h" #include #include @@ -253,23 +253,20 @@ auto createIntegratorDPD(UAMMD sim){ //Smoothed Particle Hydrodynamics auto createIntegratorSPH(UAMMD sim){ - using NVE = VerletNVE; - NVE::Parameters par; + SPHIntegrator::Parameters par; par.dt = sim.par.dt; + par.box = Box(sim.par.L); par.initVelocities = false; - auto verlet = std::make_shared(sim.pd, par); - SPH::Parameters params; - params.box = Box(sim.par.L); + //Pressure for a given particle "i" in SPH will be computed as gasStiffness·(density_i - restDensity) //Where density is computed as a function of the masses of the surroinding particles //Particle mass starts as 1, but you can change this in customizations.cuh - params.support = sim.par.support_sph; //Cut off distance for the SPH kernel - params.viscosity = sim.par.viscosity; //Environment viscosity - params.gasStiffness = sim.par.gasStiffness_sph; - params.restDensity = sim.par.restDensity_sph; - auto sph = std::make_shared(sim.pd, params); - verlet->addInteractor(sph); - return verlet; + par.support = sim.par.support_sph; //Cut off distance for the SPH kernel + par.viscosity = sim.par.viscosity; //Environment viscosity + par.gasStiffness = sim.par.gasStiffness_sph; + par.restDensity = sim.par.restDensity_sph; + auto sph = std::make_shared(sim.pd, par); + return sph; } //Creates a triply periodic Brownian Dynamics with Hydrodynamic Interactions integration module diff --git a/examples/integration_schemes/integrators.cu b/examples/integration_schemes/integrators.cu index 9d646cd9..1490df61 100644 --- a/examples/integration_schemes/integrators.cu +++ b/examples/integration_schemes/integrators.cu @@ -92,28 +92,21 @@ std::shared_ptr createIntegratorDPD(UAMMD sim){ return dpd; } -#include "Integrator/VerletNVE.cuh" -#include "Interactor/SPH.cuh" -//SPH is handled by UAMMD as a VerletNVE integrator with a special interaction +#include "Integrator/SPH.cuh" std::shared_ptr createIntegratorSPH(UAMMD sim){ - using NVE = VerletNVE; - NVE::Parameters par; - par.dt = 0.1; - par.initVelocities = false; - auto verlet = std::make_shared(sim.pd, par); - SPH::Parameters params; + SPHIntegrator::Parameters par; + par.dt = 0.1; real3 L = make_real3(32,32,32); - params.box = Box(L); + par.box = Box(L); //Pressure for a given particle "i" in SPH will be computed as gasStiffness·(density_i - restDensity) //Where density is computed as a function of the masses of the surroinding particles //Particle mass starts as 1, but you can change this in customizations.cuh - params.support = 2.4; //Cut off distance for the SPH kernel - params.viscosity = 1.0; //Environment viscosity - params.gasStiffness = 1.0; - params.restDensity = 1.0; - auto sph = std::make_shared(sim.pd, params); - verlet->addInteractor(sph); - return verlet; + par.support = 2.4; //Cut off distance for the SPH kernel + par.viscosity = 1.0; //Environment viscosity + par.gasStiffness = 1.0; + par.restDensity = 1.0; + auto sph = std::make_shared(sim.pd, par); + return sph; } #include "Integrator/Hydro/ICM.cuh" diff --git a/examples/integration_schemes/others/SPH_test.cu b/examples/integration_schemes/others/SPH_test.cu index b847656f..b40b1a5b 100644 --- a/examples/integration_schemes/others/SPH_test.cu +++ b/examples/integration_schemes/others/SPH_test.cu @@ -14,7 +14,7 @@ You can visualize the reuslts with superpunto #include"uammd.cuh" //The rest can be included depending on the used modules #include"Integrator/VerletNVE.cuh" -#include"Interactor/SPH.cuh" +#include"Integrator/SPH.cuh" #include"Interactor/ExternalForces.cuh" #include "utils/ForceEnergyVirial.cuh" #include"utils/InitialConditions.cuh" @@ -95,30 +95,28 @@ int main(int argc, char *argv[]){ auto vel = pd->getVel(access::location::gpu, access::mode::write); thrust::fill(thrust::cuda::par, vel.begin(), vel.end(), real3()); } - VerletNVE::Parameters par; + //VerletNVE with SPH interactor + SPHIntegrator::Parameters par; par.dt = dt; //If set to true (default), VerletNVE will compute Energy at step 0 and modify the velocities accordingly par.initVelocities = false; - auto verlet = make_shared(pg, par); + par.box = box; //Box to work on + //These are the default parameters, + //if any parameter is not present, it will revert to the default in the .cuh + par.support = support; + par.viscosity = viscosity; + par.gasStiffness = gasStiffness; + par.restDensity = restDensity; + auto sph_verlet = make_shared(pg, par); { //Harmonic walls acting on different particle groups //This two interactors will cause particles in group pg2 to stick to a wall in -Lz/4 //And the ones in pg3 to +Lz/4 auto extForces = make_shared>(pg, make_shared(box.boxSize)); - //Add interactors to integrator. - verlet->addInteractor(extForces); + //Add interactors to integrator (any other interactor can be added). + sph_verlet->addInteractor(extForces); } - SPH::Parameters params; - params.box = box; //Box to work on - //These are the default parameters, - //if any parameter is not present, it will revert to the default in the .cuh - params.support = support; - params.viscosity = viscosity; - params.gasStiffness = gasStiffness; - params.restDensity = restDensity; - auto sph = make_shared(pg, params); - //You can add as many modules as necessary - verlet->addInteractor(sph); + //You can issue a logging event like this, a wide variety of log levels exists (see System.cuh). //A maximum log level is set in System.cuh, every logging event with a level superior to the max will result in // absolutely no overhead, so dont be afraid to write System::DEBUGX log calls. @@ -136,7 +134,7 @@ int main(int argc, char *argv[]){ forj(0,nsteps){ //This will instruct the integrator to take the simulation to the next time step, //whatever that may mean for the particular integrator (i.e compute forces and update positions once) - verlet->forwardTime(); + sph_verlet->forwardTime(); //Write results if(printSteps and j%printSteps==0) { diff --git a/src/Integrator/SPH.cuh b/src/Integrator/SPH.cuh new file mode 100644 index 00000000..250064f0 --- /dev/null +++ b/src/Integrator/SPH.cuh @@ -0,0 +1,84 @@ +/*Pablo Palacios-Alonso 2024. SPH integrator module + + This module simulates particle dynamics using a Dissipative Particle Dynamics (DPD) simulation. + To achieve this, it employs the Verlet NVE algorithm to update the positions of the particles + in conjunction with an interactor that treats the typical DPD forces [1] as a UAMMD potential. + + + References: + [1] On the numerical treatment of dissipative particle dynamics and related systems. Leimkuhler and Shang 2015. https://doi.org/10.1016/j.jcp.2014.09.008 +*/ + +#pragma once + +#include "Integrator/VerletNVE.cuh" +#include "Interactor/SPH.cuh" + +namespace uammd{ + class SPHIntegrator: public Integrator{ + std::shared_ptr verlet; + int steps; + public: + struct Parameters{ + using NeighbourList = VerletList; + //VerletNVE + real dt = 0; + bool is2D = false; + bool initVelocities = true; //Modify initial particle velocities to enforce the provided energy + real energy = 0; //Target energy, ignored if initVelocities is false + real mass = -1; + //SPH + Box box; + real support = 1.0; + real viscosity = 50.0; + real gasStiffness = 100.0; + real restDensity = 0.4; + std::shared_ptr nl = nullptr; + }; + + SPHIntegrator(shared_ptr pg, Parameters par): + Integrator(pg, "SPHIntgrator"), + steps(0){ + //Initialize verletNVE + VerletNVE::Parameters verletpar; + verletpar.dt = par.dt; + verletpar.is2D = par.is2D; + verletpar.initVelocities = par.initVelocities; + verletpar.energy = par.energy; + + + verlet = std::make_shared(pg, verletpar); + + //Initialize SPH and add to interactor list. + SPH::Parameters sphPar; + sphPar.support = par.support; + sphPar.gasStiffness = par.gasStiffness; + sphPar.restDensity = par.restDensity; + sphPar.viscosity = par.viscosity; + sphPar.box = par.box; + + auto sph = std::make_shared(pd, sphPar); + verlet->addInteractor(sph); + } + + SPHIntegrator(shared_ptr pd, Parameters par): + SPHIntegrator(std::make_shared(pd, "All"), par){} + + ~SPHIntegrator(){} + + virtual void forwardTime() override { + steps++; + if (steps == 1){ + for(auto forceComp: interactors){ + verlet->addInteractor(forceComp); + } + } + verlet->forwardTime(); + } + + virtual real sumEnergy() override { + return verlet->sumEnergy(); + } + + }; +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b8ce90a5..1cc770b9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -72,7 +72,8 @@ add_subdirectory(utils) add_subdirectory(misc/chebyshev) add_subdirectory(misc/ibm) add_subdirectory(misc/lanczos) -add_subdirectory(DPD) +add_subdirectory(Languevin/DPD) +add_subdirectory(Languevin/SPH) IF (CMAKE_BUILD_TYPE MATCHES "Debug") if(CMAKE_CUDA_COMPILER_ID STREQUAL "NVIDIA") diff --git a/test/DPD/CMakeLists.txt b/test/Languevin/DPD/CMakeLists.txt similarity index 100% rename from test/DPD/CMakeLists.txt rename to test/Languevin/DPD/CMakeLists.txt diff --git a/test/DPD/testDPDIntegrator.cu b/test/Languevin/DPD/testDPDIntegrator.cu similarity index 100% rename from test/DPD/testDPDIntegrator.cu rename to test/Languevin/DPD/testDPDIntegrator.cu diff --git a/test/Languevin/SPH/CMakeLists.txt b/test/Languevin/SPH/CMakeLists.txt new file mode 100644 index 00000000..fad7e16b --- /dev/null +++ b/test/Languevin/SPH/CMakeLists.txt @@ -0,0 +1,11 @@ +add_executable(testSPHIntegrator testSPHIntegrator.cu) +target_include_directories(testSPHIntegrator PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(testSPHIntegrator cufft cublas) +set_target_properties(testSPHIntegrator PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin") +target_link_libraries( + testSPHIntegrator + GTest::gtest_main + GTest::gmock_main +) +include(GoogleTest) +gtest_discover_tests(testSPHIntegrator) diff --git a/test/Languevin/SPH/testSPHIntegrator.cu b/test/Languevin/SPH/testSPHIntegrator.cu new file mode 100644 index 00000000..83405a45 --- /dev/null +++ b/test/Languevin/SPH/testSPHIntegrator.cu @@ -0,0 +1,159 @@ +/* P. Palacios-Alonso 2024 + Test of the implementation of SPHIntegrator, it will simulate a LJ gas + using the new SPHIntegrator and a VerletNVE integrator with the SPH potential and check + that the reuslts are the same + */ +#include +#include +#include "gmock/gmock.h" +#include "utils/ParticleSorter.cuh" +#include "utils/container.h" +#include "utils/execution_policy.cuh" +#include +#include +#include "uammd.cuh" +#include "Integrator/SPH.cuh" +#include "Interactor/SPH.cuh" +#include "Interactor/Potential/Potential.cuh" +#include "Interactor/PairForces.cuh" +#include "utils/InitialConditions.cuh" +#include "Interactor/NeighbourList/VerletList.cuh" +#include "Interactor/NeighbourList/CellList.cuh" + + +using namespace uammd; +using std::make_shared; + +struct Parameters{ + real3 boxSize; + real dt; + int numberParticles; + int numberSteps, printSteps; + real restDensity; + real support; + real viscosity; + real gasStiffness; + real radius; + +}; + +//Let us group the UAMMD simulation in this struct that we can pass around +struct UAMMD{ + std::shared_ptr sys; + std::shared_ptr pd; + Parameters par; +}; + +//Creates the SPH integrator +auto createIntegratorSPH(UAMMD sim){ + SPHIntegrator::Parameters par; + par.support = sim.par.support; + par.viscosity = sim.par.viscosity; + par.gasStiffness = sim.par.gasStiffness; + par.restDensity = sim.par.restDensity; + par.dt = sim.par.dt; + par.box = Box(sim.par.boxSize); + par.initVelocities = false; + auto sph = make_shared(sim.pd, par); + return sph; +} + +//Creates a VerletNVE integrator and adds the SPH potential +auto createIntegratorVerletNVEWithSPH(UAMMD sim){ + SPH::Parameters par; + par.support = sim.par.support; + par.viscosity = sim.par.viscosity; + par.gasStiffness = sim.par.gasStiffness; + par.restDensity = sim.par.restDensity; + par.box = Box(sim.par.boxSize); + auto sph = make_shared(sim.pd, par); + + using NVE = VerletNVE; + NVE::Parameters params; + params.dt = sim.par.dt; + params.initVelocities = false; + auto verlet = make_shared(sim.pd, params); + verlet->addInteractor(sph); + return verlet; +} + +//Set the initial positons for the particles in a fcc lattice +void initializeParticles(UAMMD sim){ + auto pos = sim.pd->getPos(access::cpu, access::write); + auto initial = initLattice(sim.par.boxSize, sim.par.numberParticles, fcc); + std::copy(initial.begin(), initial.end(), pos.begin()); +} + +//Creates the LJ interactor +auto createLJInteractor(UAMMD sim){ + auto pot = make_shared(); + Potential::LJ::InputPairParameters par; + par.epsilon = 1.0; + par.shift = false; + par.sigma = 2*sim.par.radius; + par.cutOff = 2.5*par.sigma; + pot->setPotParameters(0, 0, par); + + using PairForces = PairForces; + PairForces::Parameters params; + params.box = Box(sim.par.boxSize); //Box to work on + auto pairforceslj = make_shared(sim.pd, params, pot); + return pairforceslj; +} + +bool compareReal3(real3 v1, real3 v2){ + return v1.x == v2.x and v1.y == v2.y and v1.z == v2.z; +} + +bool checkSamePositions(UAMMD sim1, UAMMD sim2){ + auto pos1 = sim1.pd->getPos(access::cpu, access::read); + auto pos2 = sim2.pd->getPos(access::cpu, access::read); + bool samePositions = true; + fori(0, pos1.size()){ + samePositions = compareReal3(make_real3(pos1[i]), make_real3(pos2[i])); + } + return samePositions; +} + +TEST(SPH, PositionsConsistentAcrossMethods){ + auto sys = make_shared(); + Parameters par; + par.support = 2.4; //Cut off distance for the SPH kernel + par.viscosity = 1.0; //Environment viscosity + par.gasStiffness = 1.0; + par.restDensity = 1.0; + par.radius = 1.0; + par.boxSize = make_real3(32,32,32); + par.dt = 0.001; + par.numberParticles = 1000; + par.numberSteps = 10000; + par.printSteps = 100; + //Prepare the simulation that runs the VerletNVE integrator with SPH potential + UAMMD sim1; + sim1.sys = sys; + sim1.par = par; + sim1.pd = make_shared(par.numberParticles, sim1.sys); + initializeParticles(sim1); + auto integrator1 = createIntegratorVerletNVEWithSPH(sim1); + auto lj1 = createLJInteractor(sim1); + integrator1->addInteractor(lj1); + + //Prepare the simulation that runs the SPH integrator + UAMMD sim2; + sim2.sys = sys; + sim2.par = par; + sim2.pd = make_shared(sim2.par.numberParticles, sim2.sys); + initializeParticles(sim2); + auto integrator2 = createIntegratorSPH(sim2); + auto lj2 = createLJInteractor(sim2); + integrator2->addInteractor(lj2); + //Run both simulations and check they move the particles in the same way + fori(0, par.numberSteps){ + integrator1->forwardTime(); + integrator2->forwardTime(); + if(i%par.printSteps == 0){ + ASSERT_TRUE(checkSamePositions(sim1, sim2)); + } + } + sys->finish(); +}