src/ParaMEDMEMTest/ParaMEDMEMTest_ICocoTrio.cxx

   1 // Copyright (C) 2007-2013  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
  20 #include "ParaMEDMEMTest.hxx"
  21 #include <string>
  22 #include "CommInterface.hxx"
  23 #include "ProcessorGroup.hxx"
  24 #include "MPIProcessorGroup.hxx"
  25 #include "DEC.hxx"
  26 #include "InterpKernelDEC.hxx"
  27 #include <set>
  28 #include <time.h>
  29 #include "ICoCoTrioField.hxx"
  30 #include <iostream>
  31 #include <assert.h>
  32
  33 using namespace std;
  34 using namespace ParaMEDMEM;
  35 using namespace ICoCo;
  36
  37 typedef enum {sync_and,sync_or} synctype;
  38 void synchronize_bool(bool& stop, synctype s)
  39 {
  40   int my_stop;
  41   int my_stop_temp = stop?1:0;
  42   if (s==sync_and)
  43     MPI_Allreduce(&my_stop_temp,&my_stop,1,MPI_INTEGER,MPI_MIN,MPI_COMM_WORLD);
  44   else if (s==sync_or)
  45     MPI_Allreduce(&my_stop_temp,&my_stop,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD);
  46   stop =(my_stop==1);
  47 }
  48
  49 void synchronize_dt(double& dt)
  50 {
  51   double dttemp=dt;
  52   MPI_Allreduce(&dttemp,&dt,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
  53 }
  54
  55
  56 void affiche( const TrioField&   field)
  57 {
  58   cout <<field.getName()<<endl;
  59   for (int ele=0;ele<field._nb_elems;ele++)
  60     cout <<ele <<": "<<field._field[ele]<<endl;;
  61
  62 }
  63
  64 void remplit_coord(double* coords)
  65 {
  66   coords[0*3+0]=0.;
  67   coords[0*3+1]=0.;
  68   coords[0*3+2]=0.;
  69
  70   coords[1*3+0]=1.;
  71   coords[1*3+1]=0.;
  72   coords[1*3+2]=0.;
  73
  74
  75   coords[2*3+0]=0.;
  76   coords[2*3+1]=0.;
  77   coords[2*3+2]=1.;
  78
  79   coords[3*3+0]=1.;
  80   coords[3*3+1]=0.;
  81   coords[3*3+2]=1.;
  82
  83   for (int i=4;i<8;i++)
  84     {
  85       for (int d=0;d<3;d++)
  86         coords[i*3+d]=coords[(i-4)*3+d];
  87       coords[i*3+1]+=1e-5;
  88     }
  89
  90 }
  91
  92 void init_quad(TrioField& champ_quad)
  93 {
  94
  95   champ_quad.setName("champ_quad");
  96   champ_quad._space_dim=3;
  97   champ_quad._mesh_dim=2;
  98   champ_quad._nbnodes=8;
  99   champ_quad._nodes_per_elem=4;
 100   champ_quad._nb_elems=2;
 101   champ_quad._itnumber=0;
 102   champ_quad._time1=0;
 103   champ_quad._time2=1;
 104   champ_quad._nb_field_components=1;
 105
 106   champ_quad._coords=new double[champ_quad._nbnodes*champ_quad._space_dim];
 107   //memcpy(afield._coords,sommets.addr(),champ_quad._nbnodes*champ_quad._space_dim*sizeof(double));
 108
 109   remplit_coord(champ_quad._coords);
 110
 111
 112   champ_quad._connectivity=new int[champ_quad._nb_elems*champ_quad._nodes_per_elem];
 113   champ_quad._connectivity[0*champ_quad._nodes_per_elem+0]=0;
 114   champ_quad._connectivity[0*champ_quad._nodes_per_elem+1]=1;
 115   champ_quad._connectivity[0*champ_quad._nodes_per_elem+2]=3;
 116   champ_quad._connectivity[0*champ_quad._nodes_per_elem+3]=2;
 117   champ_quad._connectivity[1*champ_quad._nodes_per_elem+0]=4;
 118   champ_quad._connectivity[1*champ_quad._nodes_per_elem+1]=5;
 119   champ_quad._connectivity[1*champ_quad._nodes_per_elem+2]=7;
 120   champ_quad._connectivity[1*champ_quad._nodes_per_elem+3]=6;
 121
 122
 123   champ_quad._has_field_ownership=false;
 124   champ_quad._field=0;
 125   //champ_quad._field=new double[champ_quad._nb_elems];
 126   //  assert(champ_quad._nb_field_components==1);
 127 }
 128 void init_triangle(TrioField& champ_triangle)
 129 {
 130
 131   champ_triangle.setName("champ_triangle");
 132   champ_triangle._space_dim=3;
 133   champ_triangle._mesh_dim=2;
 134   champ_triangle._nbnodes=8;
 135   champ_triangle._nodes_per_elem=3;
 136   champ_triangle._nb_elems=4;
 137   champ_triangle._itnumber=0;
 138   champ_triangle._time1=0;
 139   champ_triangle._time2=1;
 140   champ_triangle._nb_field_components=1;
 141
 142   champ_triangle._coords=new double[champ_triangle._nbnodes*champ_triangle._space_dim];
 143   //memcpy(afield._coords,sommets.addr(),champ_triangle._nbnodes*champ_triangle._space_dim*sizeof(double));
 144   remplit_coord(champ_triangle._coords);
 145
 146   champ_triangle._connectivity=new int[champ_triangle._nb_elems*champ_triangle._nodes_per_elem];
 147   champ_triangle._connectivity[0*champ_triangle._nodes_per_elem+0]=0;
 148   champ_triangle._connectivity[0*champ_triangle._nodes_per_elem+1]=1;
 149   champ_triangle._connectivity[0*champ_triangle._nodes_per_elem+2]=2;
 150   champ_triangle._connectivity[1*champ_triangle._nodes_per_elem+0]=1;
 151   champ_triangle._connectivity[1*champ_triangle._nodes_per_elem+1]=3;
 152   champ_triangle._connectivity[1*champ_triangle._nodes_per_elem+2]=2;
 153
 154   champ_triangle._connectivity[2*champ_triangle._nodes_per_elem+0]=4;
 155   champ_triangle._connectivity[2*champ_triangle._nodes_per_elem+1]=5;
 156   champ_triangle._connectivity[2*champ_triangle._nodes_per_elem+2]=7;
 157   champ_triangle._connectivity[3*champ_triangle._nodes_per_elem+0]=4;
 158   champ_triangle._connectivity[3*champ_triangle._nodes_per_elem+1]=7;
 159   champ_triangle._connectivity[3*champ_triangle._nodes_per_elem+2]=6;
 160
 161   champ_triangle._has_field_ownership=false;
 162   // champ_triangle._field=new double[champ_triangle._nb_elems];
 163   champ_triangle._field=0;
 164 }
 165
 166 void ParaMEDMEMTest::testICocoTrio1()
 167 {
 168   int size;
 169   int rank;
 170   MPI_Comm_size(MPI_COMM_WORLD,&size);
 171   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 172
 173   //the test is meant to run on five processors
 174   if (size !=2) return ;
 175
 176   CommInterface comm;
 177   set<int> emetteur_ids;
 178   set<int> recepteur_ids;
 179   emetteur_ids.insert(0);
 180   recepteur_ids.insert(1);
 181
 182   MPIProcessorGroup recepteur_group(comm,recepteur_ids);
 183   MPIProcessorGroup emetteur_group(comm,emetteur_ids);
 184
 185
 186   string cas;
 187   if (recepteur_group.containsMyRank())
 188     {
 189       cas="recepteur";
 190
 191     }
 192   else
 193     cas="emetteur";
 194
 195   InterpKernelDEC dec_emetteur(emetteur_group, recepteur_group);
 196
 197   TrioField champ_emetteur, champ_recepteur;
 198
 199   init_triangle(champ_emetteur);
 200   //init_triangle(champ_emetteur);
 201   init_quad(champ_recepteur);
 202   //init_emetteur(champ_recepteur);
 203
 204   if (cas=="emetteur")
 205     {
 206       champ_emetteur._field=new double[champ_emetteur._nb_elems];
 207       for (int ele=0;ele<champ_emetteur._nb_elems;ele++)
 208         champ_emetteur._field[ele]=1;
 209
 210       champ_emetteur._has_field_ownership=true;
 211     }
 212
 213
 214   MPI_Barrier(MPI_COMM_WORLD);
 215
 216   clock_t clock0= clock ();
 217   int compti=0;
 218
 219   bool init=true; // first time step ??
 220   bool stop=false;
 221   //boucle sur les pas de quads
 222   while (!stop) {
 223
 224     compti++;
 225     clock_t clocki= clock ();
 226     cout << compti << " CLOCK " << (clocki-clock0)*1.e-6 << endl;
 227     for (int non_unif=0;non_unif<2;non_unif++)
 228       {
 229         // if (champ_recepteur._field)
 230         //   delete [] champ_recepteur._field;
 231         champ_recepteur._field=0;
 232         // champ_recepteur._has_field_ownership=false;
 233
 234
 235
 236         if (cas=="emetteur")
 237           if (non_unif)
 238             champ_emetteur._field[0]=40;
 239         //bool ok=false; // Is the time interval successfully solved ?
 240
 241         // Loop on the time interval tries
 242         if(1)
 243           {
 244             if (cas=="emetteur")
 245               dec_emetteur.attachLocalField((ICoCo::Field*) &champ_emetteur);
 246             else
 247               dec_emetteur.attachLocalField((ICoCo::Field*) &champ_recepteur);
 248
 249             dec_emetteur.setNature(ConservativeVolumic);
 250
 251             if(init)
 252               dec_emetteur.synchronize();
 253             init=false;
 254
 255             if (cas=="emetteur")
 256               {
 257                 dec_emetteur.sendData();
 258                 affiche(champ_emetteur);
 259               }
 260             else if (cas=="recepteur")
 261               {
 262                 dec_emetteur.recvData();
 263                 affiche(champ_recepteur);
 264               }
 265             else
 266               throw 0;
 267           }
 268         stop=true;
 269       }
 270   }
 271 }