src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #include "ParaMEDMEMTest.hxx"
  20 #include <cppunit/TestAssert.h>
  21
  22 #include "CommInterface.hxx"
  23 #include "ProcessorGroup.hxx"
  24 #include "MPIProcessorGroup.hxx"
  25 #include "Topology.hxx"
  26 #include "DEC.hxx"
  27 #include "MxN_Mapping.hxx"
  28 #include "IntersectionDEC.hxx"
  29 #include "ParaMESH.hxx"
  30 #include "ParaFIELD.hxx"
  31 #include "ICoCoMEDField.hxx"
  32 #include "MEDLoader.hxx"
  33
  34 #include <string>
  35
  36 // use this define to enable lines, execution of which leads to Segmentation Fault
  37 #define ENABLE_FAULTS
  38
  39 // use this define to enable CPPUNIT asserts and fails, showing bugs
  40 #define ENABLE_FORCED_FAILURES
  41
  42
  43 using namespace std;
  44 using namespace ParaMEDMEM;
  45
  46 void ParaMEDMEMTest::testIntersectionDEC_2D()
  47 {
  48   testIntersectionDEC_2D_("P0","P0");
  49 }
  50
  51 void ParaMEDMEMTest::testIntersectionDEC_2DP0P1()
  52 {
  53   //testIntersectionDEC_2D_("P0","P1");
  54 }
  55
  56 /*
  57  * Check methods defined in IntersectionDEC.hxx
  58  *
  59  IntersectionDEC();
  60  IntersectionDEC(ProcessorGroup& local_group, ProcessorGroup& distant_group);
  61  virtual ~IntersectionDEC();
  62  void synchronize();
  63  void recvData();
  64  void sendData();
  65 */
  66
  67 void ParaMEDMEMTest::testIntersectionDEC_2D_(const char *srcMeth, const char *targetMeth)
  68 {
  69   std::string srcM(srcMeth);
  70   std::string targetM(targetMeth);
  71   int size;
  72   int rank;
  73   MPI_Comm_size(MPI_COMM_WORLD,&size);
  74   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  75
  76   //the test is meant to run on five processors
  77   if (size !=5) return ;
  78
  79   int nproc_source = 3;
  80   set<int> self_procs;
  81   set<int> procs_source;
  82   set<int> procs_target;
  83
  84   for (int i=0; i<nproc_source; i++)
  85     procs_source.insert(i);
  86   for (int i=nproc_source; i<size; i++)
  87     procs_target.insert(i);
  88   self_procs.insert(rank);
  89
  90   ParaMEDMEM::CommInterface interface;
  91
  92   ParaMEDMEM::ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
  93   ParaMEDMEM::ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
  94   ParaMEDMEM::ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
  95
  96   //loading the geometry for the source group
  97
  98   ParaMEDMEM::IntersectionDEC dec (*source_group,*target_group);
  99
 100   ParaMEDMEM::MEDCouplingUMesh* mesh;
 101   ParaMEDMEM::ParaMESH* paramesh;
 102   ParaMEDMEM::ParaFIELD* parafield;
 103   ICoCo::Field* icocofield ;
 104
 105   string data_dir                   = getenv("MED_ROOT_DIR");
 106   string tmp_dir                    = getenv("TMP");
 107   if (tmp_dir == "")
 108     tmp_dir = "/tmp";
 109   string filename_xml1              = data_dir + "/share/salome/resources/med/square1_split";
 110   string filename_xml2              = data_dir + "/share/salome/resources/med/square2_split";
 111   string filename_seq_wr            = tmp_dir + "/";
 112   string filename_seq_med           = tmp_dir + "/myWrField_seq_pointe221.med";
 113
 114   // To remove tmp files from disk
 115   ParaMEDMEMTest_TmpFilesRemover aRemover;
 116
 117   MPI_Barrier(MPI_COMM_WORLD);
 118   if (source_group->containsMyRank())
 119     {
 120       string master = filename_xml1;
 121
 122       ostringstream strstream;
 123       strstream <<master<<rank+1<<".med";
 124       ostringstream meshname ;
 125       meshname<< "Mesh_2_"<< rank+1;
 126
 127       mesh=MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0);
 128
 129
 130       paramesh=new ParaMESH (mesh,*source_group,"source mesh");
 131
 132       //      ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
 133       ParaMEDMEM::ComponentTopology comptopo;
 134       if(srcM=="P0")
 135         parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
 136       else
 137         parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
 138       int nb_local;
 139       if(srcM=="P0")
 140         nb_local=mesh->getNumberOfCells();
 141       else
 142         nb_local=mesh->getNumberOfNodes();
 143       //      double * value= new double[nb_local];
 144       double *value=parafield->getField()->getArray()->getPointer();
 145       for(int ielem=0; ielem<nb_local;ielem++)
 146         value[ielem]=1.0;
 147
 148       //      ICoCo::Field* icocofield=new ICoCo::MEDField(paramesh,parafield);
 149       icocofield=new ICoCo::MEDField(paramesh,parafield);
 150       dec.setMethod(srcMeth);
 151       dec.attachLocalField(icocofield);
 152     }
 153
 154   //loading the geometry for the target group
 155   if (target_group->containsMyRank())
 156     {
 157       string master= filename_xml2;
 158       ostringstream strstream;
 159       strstream << master<<(rank-nproc_source+1)<<".med";
 160       ostringstream meshname ;
 161       meshname<< "Mesh_3_"<<rank-nproc_source+1;
 162       mesh = MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0);
 163
 164       paramesh=new ParaMESH (mesh,*target_group,"target mesh");
 165       //      ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
 166       ParaMEDMEM::ComponentTopology comptopo;
 167       if(targetM=="P0")
 168         parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
 169       else
 170         parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
 171       int nb_local;
 172       if(targetM=="P0")
 173         nb_local=mesh->getNumberOfCells();
 174       else
 175         nb_local=mesh->getNumberOfNodes();
 176       //      double * value= new double[nb_local];
 177       double *value=parafield->getField()->getArray()->getPointer();
 178       for(int ielem=0; ielem<nb_local;ielem++)
 179         value[ielem]=0.0;
 180       //      ICoCo::Field* icocofield=new ICoCo::MEDField(paramesh,parafield);
 181       icocofield=new ICoCo::MEDField(paramesh,parafield);
 182       dec.setMethod(targetMeth);
 183       dec.attachLocalField(icocofield);
 184     }
 185
 186
 187   //attaching a DEC to the source group
 188   double field_before_int;
 189   double field_after_int;
 190
 191   if (source_group->containsMyRank())
 192     {
 193       field_before_int = parafield->getVolumeIntegral(0);
 194       dec.synchronize();
 195       cout<<"DEC usage"<<endl;
 196       dec.setForcedRenormalization(false);
 197
 198       dec.sendData();
 199       MEDLoader::writeParaMesh("./sourcesquareb",paramesh);
 200       if (source_group->myRank()==0)
 201         aRemover.Register("./sourcesquareb");
 202       ostringstream filename;
 203       filename<<"./sourcesquareb_"<<source_group->myRank()+1;
 204       aRemover.Register(filename.str().c_str());
 205       MEDLoader::writeParaField("./sourcesquareb","boundary",parafield);
 206
 207       dec.recvData();
 208       cout <<"writing"<<endl;
 209       MEDLoader::writeParaMesh("./sourcesquare",paramesh);
 210       if (source_group->myRank()==0)
 211         aRemover.Register("./sourcesquare");
 212       MEDLoader::writeParaField("./sourcesquare","boundary",parafield);
 213
 214
 215       filename<<"./sourcesquare_"<<source_group->myRank()+1;
 216       aRemover.Register(filename.str().c_str());
 217       field_after_int = parafield->getVolumeIntegral(0);
 218
 219
 220       //      MPI_Bcast(&field_before_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
 221       //       MPI_Bcast(&field_after_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
 222
 223       CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6);
 224
 225     }
 226
 227   //attaching a DEC to the target group
 228   if (target_group->containsMyRank())
 229     {
 230       dec.synchronize();
 231       dec.setForcedRenormalization(false);
 232
 233       dec.recvData();
 234       MEDLoader::writeParaMesh("./targetsquareb",paramesh);
 235       MEDLoader::writeParaField("./targetsquareb", "boundary",parafield);
 236       if (target_group->myRank()==0)
 237         aRemover.Register("./targetsquareb");
 238       ostringstream filename;
 239       filename<<"./targetsquareb_"<<target_group->myRank()+1;
 240       aRemover.Register(filename.str().c_str());
 241       dec.sendData();
 242       MEDLoader::writeParaMesh("./targetsquare",paramesh);
 243       MEDLoader::writeParaField("./targetsquare", "boundary",parafield);
 244
 245       if (target_group->myRank()==0)
 246         aRemover.Register("./targetsquareb");
 247
 248       filename<<"./targetsquareb_"<<target_group->myRank()+1;
 249       aRemover.Register(filename.str().c_str());
 250       //    double field_before_int, field_after_int;
 251       //       MPI_Bcast(&field_before_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
 252       //       MPI_Bcast(&field_after_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
 253
 254       //      CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6);
 255
 256     }
 257
 258   delete source_group;
 259   delete target_group;
 260   delete self_group;
 261   delete parafield;
 262   delete paramesh;
 263   mesh->decrRef();
 264
 265   delete icocofield;
 266
 267   MPI_Barrier(MPI_COMM_WORLD);
 268   cout << "end of IntersectionDEC_2D test"<<endl;
 269 }
 270
 271 //Synchronous tests without interpolation with native mode (AllToAll(v) from lam/MPI:
 272 void ParaMEDMEMTest::testSynchronousEqualIntersectionWithoutInterpNativeDEC_2D()
 273 {
 274   testAsynchronousIntersectionDEC_2D(0.1,1,0.1,1,false,false,false,"P0","P0");
 275 }
 276
 277 //Synchronous tests without interpolation :
 278 void ParaMEDMEMTest::testSynchronousEqualIntersectionWithoutInterpDEC_2D()
 279 {
 280   testAsynchronousIntersectionDEC_2D(0.1,1,0.1,1,true,false,false,"P0","P0");
 281 }
 282
 283 //Synchronous tests with interpolation :
 284 void ParaMEDMEMTest::testSynchronousEqualIntersectionDEC_2D()
 285 {
 286   testAsynchronousIntersectionDEC_2D(0.1,1,0.1,1,true,false,true,"P0","P0");
 287 }
 288 void ParaMEDMEMTest::testSynchronousFasterSourceIntersectionDEC_2D()
 289 {
 290   testAsynchronousIntersectionDEC_2D(0.09,1,0.1,1,true,false,true,"P0","P0");
 291 }
 292 void ParaMEDMEMTest::testSynchronousSlowerSourceIntersectionDEC_2D()
 293 {
 294   testAsynchronousIntersectionDEC_2D(0.11,1,0.1,1,true,false,true,"P0","P0");
 295 }
 296 void ParaMEDMEMTest::testSynchronousSlowSourceIntersectionDEC_2D()
 297 {
 298   testAsynchronousIntersectionDEC_2D(0.11,1,0.01,1,true,false,true,"P0","P0");
 299 }
 300 void ParaMEDMEMTest::testSynchronousFastSourceIntersectionDEC_2D()
 301 {
 302   testAsynchronousIntersectionDEC_2D(0.01,1,0.11,1,true,false,true,"P0","P0");
 303 }
 304
 305 //Asynchronous tests with interpolation :
 306 void ParaMEDMEMTest::testAsynchronousEqualIntersectionDEC_2D()
 307 {
 308   testAsynchronousIntersectionDEC_2D(0.1,1,0.1,1,true,true,true,"P0","P0");
 309 }
 310 void ParaMEDMEMTest::testAsynchronousFasterSourceIntersectionDEC_2D()
 311 {
 312   testAsynchronousIntersectionDEC_2D(0.09,1,0.1,1,true,true,true,"P0","P0");
 313 }
 314 void ParaMEDMEMTest::testAsynchronousSlowerSourceIntersectionDEC_2D()
 315 {
 316   testAsynchronousIntersectionDEC_2D(0.11,1,0.1,1,true,true,true,"P0","P0");
 317 }
 318 void ParaMEDMEMTest::testAsynchronousSlowSourceIntersectionDEC_2D()
 319 {
 320   testAsynchronousIntersectionDEC_2D(0.11,1,0.01,1,true,true,true,"P0","P0");
 321 }
 322 void ParaMEDMEMTest::testAsynchronousFastSourceIntersectionDEC_2D()
 323 {
 324   testAsynchronousIntersectionDEC_2D(0.01,1,0.11,1,true,true,true,"P0","P0");
 325 }
 326
 327 /*!
 328  * Tests an asynchronous exchange between two codes
 329  * one sends data with dtA as an interval, the max time being tmaxA
 330  * the other one receives with dtB as an interval, the max time being tmaxB
 331  */
 332 void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA,
 333                                                         double dtB, double tmaxB, bool WithPointToPoint, bool Asynchronous,
 334                                                         bool WithInterp, const char *srcMeth, const char *targetMeth)
 335 {
 336   std::string srcM(srcMeth);
 337   std::string targetM(targetMeth);
 338   int size;
 339   int rank;
 340   MPI_Comm_size(MPI_COMM_WORLD,&size);
 341   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 342
 343   //the test is meant to run on five processors
 344   if (size !=5) return ;
 345
 346   int nproc_source = 3;
 347   set<int> self_procs;
 348   set<int> procs_source;
 349   set<int> procs_target;
 350
 351   for (int i=0; i<nproc_source; i++)
 352     procs_source.insert(i);
 353   for (int i=nproc_source; i<size; i++)
 354     procs_target.insert(i);
 355   self_procs.insert(rank);
 356
 357   ParaMEDMEM::CommInterface interface;
 358
 359   ParaMEDMEM::ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
 360   ParaMEDMEM::ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
 361   ParaMEDMEM::ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
 362
 363   //loading the geometry for the source group
 364
 365   ParaMEDMEM::IntersectionDEC dec (*source_group,*target_group);
 366
 367   ParaMEDMEM::MEDCouplingUMesh* mesh;
 368   ParaMEDMEM::ParaMESH* paramesh;
 369   ParaMEDMEM::ParaFIELD* parafield;
 370
 371   double * value ;
 372   ICoCo::Field* icocofield ;
 373
 374   string data_dir                   = getenv("MED_ROOT_DIR");
 375   string tmp_dir                    = getenv("TMP");
 376   if (tmp_dir == "")
 377     tmp_dir = "/tmp";
 378   string filename_xml1              = data_dir + "/share/salome/resources/med/square1_split";
 379   string filename_xml2              = data_dir + "/share/salome/resources/med/square2_split";
 380   string filename_seq_wr            = tmp_dir + "/";
 381   string filename_seq_med           = tmp_dir + "/myWrField_seq_pointe221.med";
 382
 383   // To remove tmp files from disk
 384   ParaMEDMEMTest_TmpFilesRemover aRemover;
 385
 386   MPI_Barrier(MPI_COMM_WORLD);
 387
 388   if (source_group->containsMyRank())
 389     {
 390       string master = filename_xml1;
 391
 392       ostringstream strstream;
 393       strstream <<master<<rank+1<<".med";
 394       ostringstream meshname ;
 395       meshname<< "Mesh_2_"<< rank+1;
 396
 397       mesh=MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0);
 398
 399       paramesh=new ParaMESH (mesh,*source_group,"source mesh");
 400
 401       //      ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
 402       ParaMEDMEM::ComponentTopology comptopo;
 403       if(srcM=="P0")
 404         parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
 405       else
 406         parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
 407
 408       int nb_local;
 409       if(srcM=="P0")
 410         nb_local=mesh->getNumberOfCells();
 411       else
 412         nb_local=mesh->getNumberOfNodes();
 413       //      double * value= new double[nb_local];
 414       double *value=parafield->getField()->getArray()->getPointer();
 415       for(int ielem=0; ielem<nb_local;ielem++)
 416         value[ielem]=0.0;
 417
 418       //      ICoCo::Field* icocofield=new ICoCo::MEDField(paramesh,parafield);
 419       icocofield=new ICoCo::MEDField(paramesh,parafield);
 420
 421       dec.attachLocalField(icocofield);
 422
 423
 424     }
 425
 426   //loading the geometry for the target group
 427   if (target_group->containsMyRank())
 428     {
 429       string master= filename_xml2;
 430       ostringstream strstream;
 431       strstream << master<<(rank-nproc_source+1)<<".med";
 432       ostringstream meshname ;
 433       meshname<< "Mesh_3_"<<rank-nproc_source+1;
 434
 435       mesh = MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0);
 436
 437       paramesh=new ParaMESH (mesh,*target_group,"target mesh");
 438       //      ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
 439       ParaMEDMEM::ComponentTopology comptopo;
 440       if(targetM=="P0")
 441         parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
 442       else
 443         parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
 444
 445       int nb_local;
 446       if(targetM=="P0")
 447         nb_local=mesh->getNumberOfCells();
 448       else
 449         nb_local=mesh->getNumberOfNodes();
 450
 451       double *value=parafield->getField()->getArray()->getPointer();
 452       for(int ielem=0; ielem<nb_local;ielem++)
 453         value[ielem]=0.0;
 454       //      ICoCo::Field* icocofield=new ICoCo::MEDField(paramesh,parafield);
 455       icocofield=new ICoCo::MEDField(paramesh,parafield);
 456
 457       dec.attachLocalField(icocofield);
 458     }
 459
 460
 461   //attaching a DEC to the source group
 462
 463   if (source_group->containsMyRank())
 464     {
 465       cout<<"DEC usage"<<endl;
 466       dec.setAsynchronous(Asynchronous);
 467       if ( WithInterp ) {
 468         dec.setTimeInterpolationMethod(LinearTimeInterp);
 469       }
 470       if ( WithPointToPoint ) {
 471         dec.setAllToAllMethod(PointToPoint);
 472       }
 473       else {
 474         dec.setAllToAllMethod(Native);
 475       }
 476       dec.synchronize();
 477       dec.setForcedRenormalization(false);
 478       for (double time=0; time<tmaxA+1e-10; time+=dtA)
 479         {
 480           cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time
 481                << " dtA " << dtA << " tmaxA " << tmaxA << endl ;
 482           if ( time+dtA < tmaxA+1e-7 ) {
 483             dec.sendData( time , dtA );
 484           }
 485           else {
 486             dec.sendData( time , 0 );
 487           }
 488           double* value = parafield->getField()->getArray()->getPointer();
 489           int nb_local=parafield->getField()->getMesh()->getNumberOfCells();
 490           for (int i=0; i<nb_local;i++)
 491             value[i]= time+dtA;
 492
 493
 494         }
 495     }
 496
 497   //attaching a DEC to the target group
 498   if (target_group->containsMyRank())
 499     {
 500       cout<<"DEC usage"<<endl;
 501       dec.setAsynchronous(Asynchronous);
 502       if ( WithInterp ) {
 503         dec.setTimeInterpolationMethod(LinearTimeInterp);
 504       }
 505       if ( WithPointToPoint ) {
 506         dec.setAllToAllMethod(PointToPoint);
 507       }
 508       else {
 509         dec.setAllToAllMethod(Native);
 510       }
 511       dec.synchronize();
 512       dec.setForcedRenormalization(false);
 513       vector<double> times;
 514       for (double time=0; time<tmaxB+1e-10; time+=dtB)
 515         {
 516           cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time
 517                << " dtB " << dtB << " tmaxB " << tmaxB << endl ;
 518           dec.recvData( time );
 519           double vi = parafield->getVolumeIntegral(0);
 520           cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time
 521                << " VolumeIntegral " << vi
 522                << " time*10000 " << time*10000 << endl ;
 523
 524           CPPUNIT_ASSERT_DOUBLES_EQUAL(vi,time*10000,0.001);
 525         }
 526
 527     }
 528
 529   delete source_group;
 530   delete target_group;
 531   delete self_group;
 532   delete parafield ;
 533   delete paramesh ;
 534   mesh->decrRef() ;
 535   delete icocofield ;
 536
 537   cout << "testAsynchronousIntersectionDEC_2D" << rank << " MPI_Barrier " << endl ;
 538
 539   if (Asynchronous) MPI_Barrier(MPI_COMM_WORLD);
 540   cout << "end of IntersectionDEC_2D test"<<endl;
 541 }