test/unit/test_nvt_integrator_mpi.cc

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2009-2014 The Regents of
   4 the University of Michigan All rights reserved.
   5
   6 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   7 copyright is held, by various Contributors who have granted The Regents of the
   8 University of Michigan the right to modify and/or distribute such Contributions.
   9
  10 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  11 and binary forms, provided you abide by the following conditions:
  12
  13 * Redistributions of source code must retain the above copyright notice, this
  14 list of conditions, and the following disclaimer both in the code and
  15 prominently in any materials provided with the distribution.
  16
  17 * Redistributions in binary form must reproduce the above copyright notice, this
  18 list of conditions, and the following disclaimer in the documentation and/or
  19 other materials provided with the distribution.
  20
  21 * All publications and presentations based on HOOMD-blue, including any reports
  22 or published results obtained, in whole or in part, with HOOMD-blue, will
  23 acknowledge its use according to the terms posted at the time of submission on:
  24 http://codeblue.umich.edu/hoomd-blue/citations.html
  25
  26 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  27 http://codeblue.umich.edu/hoomd-blue/
  28
  29 * Apart from the above required attributions, neither the name of the copyright
  30 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  31 promote products derived from this software without specific prior written
  32 permission.
  33
  34 Disclaimer
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  37 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  38 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  39 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  40
  41 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  42 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  43 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  44 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  46 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  47 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48 */
  49
  50 //! name the boost unit test module
  51 #define BOOST_TEST_MODULE NVTUpdaterTestsMPI
  52 #include "boost_utf_configure.h"
  53
  54 #include "HOOMDMath.h"
  55 #include "ExecutionConfiguration.h"
  56 #include "System.h"
  57 #include "TwoStepNVT.h"
  58 #include "IntegratorTwoStep.h"
  59 #include "AllPairPotentials.h"
  60 #include "NeighborListBinned.h"
  61 #include "RandomGenerator.h"
  62
  63 #include <boost/python.hpp>
  64 #include <boost/mpi.hpp>
  65 #include <boost/shared_ptr.hpp>
  66
  67 #include <math.h>
  68
  69 #include "Communicator.h"
  70 #include "DomainDecomposition.h"
  71
  72 #ifdef ENABLE_CUDA
  73 #include "CellListGPU.h"
  74 #include "NeighborListGPUBinned.h"
  75 #include "TwoStepNVTGPU.h"
  76 #include "ComputeThermoGPU.h"
  77 #include "CommunicatorGPU.h"
  78 #endif
  79
  80 using namespace boost;
  81
  82 void test_nvt_integrator_mpi(boost::shared_ptr<ExecutionConfiguration> exec_conf)
  83 {
  84     // initialize random particle system
  85     Scalar phi_p = 0.2;
  86     unsigned int N = 20000;
  87     Scalar L = pow(M_PI/6.0/phi_p*Scalar(N),1.0/3.0);
  88     BoxDim box_g(L);
  89     RandomGenerator rand_init(exec_conf, box_g, 12345, 3);
  90     std::vector<string> types;
  91     types.push_back("A");
  92     std::vector<unsigned int> bonds;
  93     std::vector<string> bond_types;
  94     rand_init.addGenerator((int)N, boost::shared_ptr<PolymerParticleGenerator>(new PolymerParticleGenerator(exec_conf, 1.0, types, bonds, bonds, bond_types, 100, 3)));
  95     rand_init.setSeparationRadius("A", .4);
  96
  97     rand_init.generate();
  98
  99     boost::shared_ptr<SnapshotSystemData> snap;
 100     snap = rand_init.getSnapshot();
 101
 102     boost::shared_ptr<DomainDecomposition> decomposition(new DomainDecomposition(exec_conf,snap->global_box.getL(), 0));
 103
 104     boost::shared_ptr<SystemDefinition> sysdef_1(new SystemDefinition(snap, exec_conf,decomposition));
 105
 106     // initialize a second system (single proc) on rank zero
 107     boost::shared_ptr<SystemDefinition> sysdef_2;
 108     if (exec_conf->getRank() == 0)
 109         sysdef_2 = boost::shared_ptr<SystemDefinition>(new SystemDefinition(snap, exec_conf));
 110
 111     boost::shared_ptr<ParticleData> pdata_1 = sysdef_1->getParticleData();
 112
 113     boost::shared_ptr<ParticleData> pdata_2;
 114     if (exec_conf->getRank() == 0)
 115         pdata_2 = sysdef_2->getParticleData();
 116
 117     boost::shared_ptr<Communicator> comm;
 118 #ifdef ENABLE_CUDA
 119     if (exec_conf->isCUDAEnabled())
 120         comm = boost::shared_ptr<Communicator>(new CommunicatorGPU(sysdef_1, decomposition));
 121     else
 122 #endif
 123         comm = boost::shared_ptr<Communicator>(new Communicator(sysdef_1,decomposition));
 124
 125     boost::shared_ptr<ParticleSelector> selector_all_1(new ParticleSelectorTag(sysdef_1, 0, pdata_1->getNGlobal()-1));
 126     boost::shared_ptr<ParticleGroup> group_all_1(new ParticleGroup(sysdef_1, selector_all_1));
 127
 128     boost::shared_ptr<ParticleSelector> selector_all_2;
 129     boost::shared_ptr<ParticleGroup> group_all_2;
 130     if (exec_conf->getRank() ==0)
 131         {
 132         selector_all_2 = boost::shared_ptr<ParticleSelector>(new ParticleSelectorTag(sysdef_2, 0, pdata_2->getNGlobal()-1));
 133         group_all_2 = boost::shared_ptr<ParticleGroup>(new ParticleGroup(sysdef_2, selector_all_2));
 134         }
 135
 136     Scalar r_cut = Scalar(3.0);
 137     Scalar r_buff = Scalar(0.8);
 138     boost::shared_ptr<NeighborList> nlist_1(new NeighborListBinned(sysdef_1, r_cut, r_buff));
 139
 140     nlist_1->setStorageMode(NeighborList::full);
 141     nlist_1->setCommunicator(comm);
 142     boost::shared_ptr<PotentialPairLJ> fc_1 = boost::shared_ptr<PotentialPairLJ>(new PotentialPairLJ(sysdef_1, nlist_1));
 143
 144     fc_1->setRcut(0, 0, r_cut);
 145
 146     // setup some values for alpha and sigma
 147     Scalar epsilon = Scalar(1.0);
 148     Scalar sigma = Scalar(1.0);
 149     Scalar alpha = Scalar(1.0);
 150     Scalar lj1 = Scalar(4.0) * epsilon * pow(sigma,Scalar(12.0));
 151     Scalar lj2 = alpha * Scalar(4.0) * epsilon * pow(sigma,Scalar(6.0));
 152     fc_1->setParams(0,0,make_scalar2(lj1,lj2));
 153
 154     boost::shared_ptr<NeighborList> nlist_2;
 155     boost::shared_ptr<PotentialPairLJ> fc_2;
 156     if (exec_conf->getRank() == 0)
 157         {
 158         nlist_2 = boost::shared_ptr<NeighborList>(new NeighborListBinned(sysdef_2, r_cut, r_buff));
 159         nlist_2->setStorageMode(NeighborList::full);
 160         fc_2 = boost::shared_ptr<PotentialPairLJ>(new PotentialPairLJ(sysdef_2, nlist_2));
 161         fc_2->setRcut(0, 0, r_cut);
 162         fc_2->setParams(0,0,make_scalar2(lj1,lj2));
 163         }
 164
 165     Scalar deltaT = Scalar(0.005);
 166     Scalar Q = Scalar(2.0);
 167     Scalar T = Scalar(1.5/3.0);
 168     Scalar tau = sqrt(Q / (Scalar(3.0) * T));
 169     boost::shared_ptr<VariantConst> T_variant_1(new VariantConst(T));
 170     boost::shared_ptr<IntegratorTwoStep> nvt_1(new IntegratorTwoStep(sysdef_1, deltaT));
 171     boost::shared_ptr<ComputeThermo> thermo_1 = boost::shared_ptr<ComputeThermo>(new ComputeThermo(sysdef_1,group_all_1));
 172     thermo_1->setCommunicator(comm);
 173
 174     boost::shared_ptr<VariantConst> T_variant_2;
 175     boost::shared_ptr<IntegratorTwoStep> nvt_2;
 176     boost::shared_ptr<ComputeThermo> thermo_2;
 177
 178     if (exec_conf->getRank()==0)
 179         {
 180         T_variant_2 = boost::shared_ptr<VariantConst>(new VariantConst(T));
 181         nvt_2  = boost::shared_ptr<IntegratorTwoStep>(new IntegratorTwoStep(sysdef_2, deltaT));
 182         thermo_2 = boost::shared_ptr<ComputeThermo>(new ComputeThermo(sysdef_2,group_all_2));
 183         }
 184
 185     boost::shared_ptr<TwoStepNVT> two_step_nvt_1;
 186     boost::shared_ptr<TwoStepNVT> two_step_nvt_2;
 187 #ifdef ENABLE_CUDA
 188     if (exec_conf->isCUDAEnabled())
 189         {
 190         two_step_nvt_1 = boost::shared_ptr<TwoStepNVT>(new TwoStepNVTGPU(sysdef_1, group_all_1, thermo_1, tau, T_variant_1));
 191         if (exec_conf->getRank() == 0)
 192             two_step_nvt_2 = boost::shared_ptr<TwoStepNVT>(new TwoStepNVTGPU(sysdef_2, group_all_2, thermo_2, tau, T_variant_2));
 193         }
 194     else
 195 #endif
 196         {
 197         two_step_nvt_1 = boost::shared_ptr<TwoStepNVT>(new TwoStepNVT(sysdef_1, group_all_1, thermo_1, tau, T_variant_1));
 198         if (exec_conf->getRank() == 0)
 199             two_step_nvt_2 = boost::shared_ptr<TwoStepNVT>(new TwoStepNVT(sysdef_2, group_all_2, thermo_2, tau, T_variant_2));
 200         }
 201
 202     nvt_1->addIntegrationMethod(two_step_nvt_1);
 203     nvt_1->addForceCompute(fc_1);
 204     nvt_1->setCommunicator(comm);
 205
 206     if (exec_conf->getRank() == 0)
 207         {
 208         nvt_2->addIntegrationMethod(two_step_nvt_2);
 209         nvt_2->addForceCompute(fc_2);
 210         }
 211
 212     unsigned int ndof = nvt_1->getNDOF(group_all_1);
 213     if (exec_conf->getRank() == 0)
 214         BOOST_CHECK_EQUAL(ndof, nvt_2->getNDOF(group_all_2));
 215
 216     thermo_1->setNDOF(ndof);
 217
 218     if (exec_conf->getRank() == 0)
 219         thermo_2->setNDOF(ndof);
 220
 221     nvt_1->prepRun(0);
 222
 223     if (exec_conf->getRank() == 0)
 224         nvt_2->prepRun(0);
 225
 226     SnapshotParticleData snap_1(N);
 227     SnapshotParticleData snap_2(N);
 228     for (int i=0; i< 100; i++)
 229         {
 230         // compare temperatures
 231         if (exec_conf->getRank() == 0)
 232             BOOST_CHECK_CLOSE(thermo_1->getTemperature(), thermo_2->getTemperature(), tol_small);
 233
 234 //       if (world->rank() ==0)
 235 //           std::cout << "step " << i << std::endl;
 236         Scalar rough_tol = 15.0;
 237         Scalar abs_tol = 1e-5;
 238
 239         // in the first five steps, compare all accelerations and velocities
 240         // beyond this number of steps, trajectories will generally diverge, since they are chaotic
 241         // compare the snapshot of the parallel simulation
 242         if (i < 5)
 243             {
 244             pdata_1->takeSnapshot(snap_1);
 245             // ... against the serial simulation
 246
 247             if (exec_conf->getRank() == 0)
 248                 {
 249                 pdata_2->takeSnapshot(snap_2);
 250                 // check position, velocity and acceleration
 251                 for (unsigned int j = 0; j < N; j++)
 252                     {
 253                     // we do not check positions (or we would need to pull back vectors over the boundaries)
 254                     //MY_BOOST_CHECK_CLOSE(snap_1.pos[j].x, snap_2.pos[j].x, rough_tol);
 255                     //MY_BOOST_CHECK_CLOSE(snap_1.pos[j].y, snap_2.pos[j].y, rough_tol);
 256                     //MY_BOOST_CHECK_CLOSE(snap_1.pos[j].z, snap_2.pos[j].z, rough_tol);
 257
 258                     if (fabsf(snap_1.vel[j].x) < abs_tol)
 259                         BOOST_CHECK_SMALL(snap_2.vel[j].x, 2*abs_tol);
 260                     else
 261                         MY_BOOST_CHECK_CLOSE(snap_1.vel[j].x, snap_2.vel[j].x, rough_tol);
 262
 263                     if (fabsf(snap_1.vel[j].y) < abs_tol)
 264                         BOOST_CHECK_SMALL(snap_2.vel[j].y, 2*abs_tol);
 265                     else
 266                         MY_BOOST_CHECK_CLOSE(snap_1.vel[j].y, snap_2.vel[j].y, rough_tol);
 267
 268                     if (fabsf(snap_1.vel[j].z) < abs_tol)
 269                         BOOST_CHECK_SMALL(snap_2.vel[j].z, 2*abs_tol);
 270                     else
 271                         MY_BOOST_CHECK_CLOSE(snap_1.vel[j].z, snap_2.vel[j].z, rough_tol);
 272
 273                     if (fabsf(snap_1.accel[j].x) < abs_tol)
 274                         BOOST_CHECK_SMALL(snap_2.accel[j].x, 2*abs_tol);
 275                     else
 276                         MY_BOOST_CHECK_CLOSE(snap_1.accel[j].x, snap_2.accel[j].x, rough_tol);
 277
 278                     if (fabsf(snap_1.accel[j].y) < abs_tol)
 279                         BOOST_CHECK_SMALL(snap_2.accel[j].y, 2*abs_tol);
 280                     else
 281                         MY_BOOST_CHECK_CLOSE(snap_1.accel[j].y, snap_2.accel[j].y, rough_tol);
 282
 283                     if (fabsf(snap_1.accel[j].z) < abs_tol)
 284                         BOOST_CHECK_SMALL(snap_2.accel[j].z, 2*abs_tol);
 285                     else
 286                         MY_BOOST_CHECK_CLOSE(snap_1.accel[j].z, snap_2.accel[j].z, rough_tol);
 287                     }
 288                 }
 289             }
 290         nvt_1->update(i);
 291         if (exec_conf->getRank() == 0)
 292             nvt_2->update(i);
 293         }
 294
 295 }
 296
 297 //! Tests MPI domain decomposition with NVT integrator
 298 BOOST_AUTO_TEST_CASE( DomainDecomposition_NVT_test )
 299     {
 300     test_nvt_integrator_mpi(boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 301     }
 302
 303 #ifdef ENABLE_CUDA
 304 //! Tests MPI domain decomposition with NVT integrator on the GPU
 305 BOOST_AUTO_TEST_CASE( DomainDecomposition_NVT_test_GPU )
 306     {
 307     test_nvt_integrator_mpi(boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 308     }
 309 #endif