Ignore:
Timestamp:
Apr 10, 2012, 1:55:49 PM (14 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
a40eea
Parents:
d24c2f
Message:

Merge recent changes of the vmg library into ScaFaCos.

Includes a fix for the communication problems on Jugene.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1666 5161e1c8-67bf-11de-9fd5-51895aff932f

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/comm/comm_serial.cpp

    rd24c2f rac6d04  
    2828#include "base/vector.hpp"
    2929#include "comm/comm_serial.hpp"
    30 #include "comm/comm.hpp"
    3130#include "grid/multigrid.hpp"
    3231#include "grid/tempgrid.hpp"
    … …  
    4342}
    4443
    45 void CommSerial::CommSubgrid(Grid& grid_old, Grid& grid_new)
     44void CommSerial::CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
    4645{
    4746  Grid::iterator iter;
    … …  
    5453}
    5554
    56 void CommSerial::CommAddSubgrid(Grid& grid_old, Grid& grid_new)
     55void CommSerial::CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
    5756{
    5857  Grid::iterator iter;
    … …  
    7574}
    7675
    77 void CommSerial::CommFromGhosts(Grid& mesh)
     76void CommSerial::CommFromGhosts(Grid& grid)
    7877{
    7978  Grid::iterator iter;
    8079
    81   const GridIteratorSuite& iterators = mesh.Iterators();
    82   const Index& size = mesh.Local().Size();
     80  const GridIteratorSuite& iterators = grid.Iterators();
     81  const Index& size = grid.Local().Size();
    8382
    8483  if (BoundaryConditions().X() == Periodic) {
    8584
    8685    for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
    87       mesh((*iter).X()+size.X(), (*iter).Y(), (*iter).Z()) += mesh.GetVal(*iter);
     86      grid((*iter).X()+size.X(), (*iter).Y(), (*iter).Z()) += grid.GetVal(*iter);
    8887
    8988    for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
    90       mesh((*iter).X()-size.X(), (*iter).Y(), (*iter).Z()) += mesh.GetVal(*iter);
     89      grid((*iter).X()-size.X(), (*iter).Y(), (*iter).Z()) += grid.GetVal(*iter);
    9190  }
    9291
    … …  
    9493
    9594    for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
    96       mesh((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z()) += mesh.GetVal(*iter);
     95      grid((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z()) += grid.GetVal(*iter);
    9796
    9897    for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
    99       mesh((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z()) += mesh.GetVal(*iter);
     98      grid((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z()) += grid.GetVal(*iter);
    10099
    101100  }
    … …  
    104103
    105104    for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
    106       mesh((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z()) += mesh.GetVal(*iter);
     105      grid((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z()) += grid.GetVal(*iter);
    107106
    108107    for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
    109       mesh((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z()) += mesh.GetVal(*iter);
     108      grid((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z()) += grid.GetVal(*iter);
    110109
    111110  }
    112111
    113112#ifdef DEBUG_MATRIX_CHECKS
    114   mesh.IsConsistent();
     113  grid.IsConsistent();
    115114#endif
    116115}
    117116
    118 void CommSerial::CommToGhosts(Grid& mesh)
     117void CommSerial::CommToGhosts(Grid& grid)
    119118{
    120119  Grid::iterator iter;
    121120
    122   const GridIteratorSuite& iterators = mesh.Iterators();
    123   const Index& size = mesh.Local().Size();
     121  const GridIteratorSuite& iterators = grid.Iterators();
     122  const Index& size = grid.Local().Size();
    124123
    125124  if (BoundaryConditions().Z() == Periodic) {
    126125
    127126    for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
    128       mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z());
     127      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z());
    129128
    130129    for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
    131       mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z());
     130      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z());
    132131
    133132  }
    … …  
    136135
    137136    for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
    138       mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z());
     137      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z());
    139138
    140139    for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
    141       mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z());
     140      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z());
    142141
    143142  }
    … …  
    146145
    147146    for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
    148       mesh(*iter) = mesh.GetVal((*iter).X()+size.X(), (*iter).Y(), (*iter).Z());
     147      grid(*iter) = grid.GetVal((*iter).X()+size.X(), (*iter).Y(), (*iter).Z());
    149148
    150149    for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
    151       mesh(*iter) = mesh.GetVal((*iter).X()-size.X(), (*iter).Y(), (*iter).Z());
     150      grid(*iter) = grid.GetVal((*iter).X()-size.X(), (*iter).Y(), (*iter).Z());
    152151
    153152  }
    154153
    155154#ifdef DEBUG_MATRIX_CHECKS
    156   mesh.IsConsistent();
     155  grid.IsConsistent();
    157156#endif
    158157}
    159158
    160 void CommSerial::CommToGhostsRB(Grid& grid, const int& offset)
     159void CommSerial::CommToGhostsAsyncStart(Grid& grid)
    161160{
    162161  CommToGhosts(grid);
    163162}
    164163
    165 void CommSerial::CommToGhostsAsyncStart(Grid& grid)
    166 {
    167   CommToGhosts(grid);
    168 }
    169 
    170 void CommSerial::CommToGhostsAsyncFinish()
     164void CommSerial::CommToGhostsAsyncFinish(Grid& grid)
    171165{
    172166}
    … …  
    179173void CommSerial::CommFromGhostsAsyncFinish(Grid& grid)
    180174{
    181 }
    182 
    183 Grid& CommSerial::GetGlobalCoarseGrid(Multigrid& multigrid)
    184 {
    185   return multigrid(multigrid.MinLevel());
    186175}
    187176
    … …  
    199188}
    200189
    201 void CommSerial::CommLCListGhosts(const Grid& grid, Particle::LinkedCellList& lc)
    202 {
    203   std::list<Particle::Particle>::iterator p_iter;
    204   Grid::iterator g_iter;
    205 
    206   for (int i=2; i>=0; --i) {
    207 
    208     if (BoundaryConditions()[i] == Periodic) {
    209 
    210       for (g_iter=lc.Iterators().Boundary1()[i].Begin(); g_iter!=lc.Iterators().Boundary1()[i].End(); ++g_iter)
    211         lc(*g_iter).clear();
    212 
    213       for (g_iter=lc.Iterators().NearBoundary1()[i].Begin(); g_iter!=lc.Iterators().NearBoundary1()[i].End(); ++g_iter)
    214         for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
    215           lc.AddParticleToComplete(p_iter->Pos() + grid.Extent().Size(), p_iter->Charge(),
    216                                    p_iter->Pot(), p_iter->Force(),
    217                                    p_iter->Rank(), p_iter->Index());
    218 
    219       for (g_iter=lc.Iterators().Boundary2()[i].Begin(); g_iter!=lc.Iterators().Boundary2()[i].End(); ++g_iter)
    220         lc(*g_iter).clear();
    221 
    222       for (g_iter=lc.Iterators().NearBoundary2()[i].Begin(); g_iter!=lc.Iterators().NearBoundary2()[i].End(); ++g_iter)
    223         for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
    224           lc.AddParticleToComplete(p_iter->Pos() - grid.Extent().Size(), p_iter->Charge(),
    225                                    p_iter->Pot(), p_iter->Force(),
    226                                    p_iter->Rank(), p_iter->Index());
    227 
    228     }
    229   }
    230 }
    231 
    232 void CommSerial::CommAddPotential(std::list<Particle::Particle>& particles)
     190void CommSerial::CommParticlesBack(std::list<Particle::Particle>& particles)
    233191{
    234192  std::list<Particle::Particle>::iterator iter;
    235193  vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
    236 
    237194  for (iter=particles.begin(); iter!=particles.end(); ++iter)
    238     p[iter->Index()] += iter->Pot();
    239 }
    240 
    241 void CommSerial::CommAddPotential(Particle::LinkedCellList& lc)
    242 {
    243   Grid::iterator lc_iter;
    244   std::list<Particle::Particle>::iterator p_iter;
    245   vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
    246 
    247   for (lc_iter=lc.Iterators().Local().Begin(); lc_iter!=lc.Iterators().Local().End(); ++lc_iter)
    248     for (p_iter=lc(*lc_iter).begin(); p_iter!=lc(*lc_iter).end(); ++p_iter)
    249       p[p_iter->Index()] += p_iter->Pot();
    250 }
     195    p[iter->Index()] = iter->Pot();
     196}
     197
     198void CommSerial::CommLCListToGhosts(const Grid& grid, Particle::LinkedCellList& lc)
     199{
     200  // std::list<Particle::Particle*>::iterator p_iter;
     201  // Grid::iterator g_iter;
     202
     203  // for (int i=2; i>=0; --i) {
     204
     205  //   if (BoundaryConditions()[i] == Periodic) {
     206
     207  //     for (g_iter=lc.Iterators().Halo1()[i].Begin(); g_iter!=lc.Iterators().Halo1()[i].End(); ++g_iter)
     208  //    lc(*g_iter).clear();
     209
     210  //     for (g_iter=lc.Iterators().NearBoundary1()[i].Begin(); g_iter!=lc.Iterators().NearBoundary1()[i].End(); ++g_iter)
     211  //    for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
     212  //      lc.AddParticleToComplete((*p_iter)->Pos() + grid.Extent().Size(), (*p_iter)->Charge(),
     213  //                               (*p_iter)->Pot(), (*p_iter)->Force(),
     214  //                               (*p_iter)->Rank(), (*p_iter)->Index());
     215
     216  //     for (g_iter=lc.Iterators().Halo2()[i].Begin(); g_iter!=lc.Iterators().Halo2()[i].End(); ++g_iter)
     217  //    lc(*g_iter).clear();
     218
     219  //     for (g_iter=lc.Iterators().NearBoundary2()[i].Begin(); g_iter!=lc.Iterators().NearBoundary2()[i].End(); ++g_iter)
     220  //    for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
     221  //      lc.AddParticleToComplete((*p_iter)->Pos() - grid.Extent().Size(), (*p_iter)->Charge(),
     222  //                               (*p_iter)->Pot(), (*p_iter)->Force(),
     223  //                               (*p_iter)->Rank(), (*p_iter)->Index());
     224
     225  //   }
     226  // }
     227}
     228
     229void CommSerial::CommLCListFromGhosts(const Grid& grid, Particle::LinkedCellList& lc)
     230{
     231
     232}
     233
     234// void CommSerial::CommAddPotential(std::list<Particle::Particle>& particles)
     235// {
     236//   std::list<Particle::Particle>::iterator iter;
     237//   vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
     238
     239//   for (iter=particles.begin(); iter!=particles.end(); ++iter)
     240//     p[iter->Index()] += iter->Pot();
     241// }
    251242
    252243void CommSerial::PrintString(const char* format, ...)
    … …  
    286277}
    287278
    288 void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& mesh, const char* information, bool inner)
     279void CommSerial::PrintGrid(Grid& grid, const char* information)
     280{
     281  Index i;
     282  std::stringstream out;
     283  std::ofstream out_file;
     284
     285  OpenFileAndPrintHeader(out_file, grid, information);
     286
     287  for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z())
     288    for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
     289      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
     290        out << std::scientific << grid.GetVal(i) << std::endl;
     291
     292  out_file << out.str();
     293  out_file.close();
     294}
     295
     296void CommSerial::PrintDefect(Grid& sol, Grid& rhs, const char* information)
     297{
     298  Grid::iterator iter;
     299  std::ofstream out_file;
     300  std::stringstream out;
     301
     302  TempGrid *temp = MG::GetTempGrid();
     303  temp->SetProperties(sol);
     304  temp->ImportFromResidual(sol, rhs);
     305
     306  OpenFileAndPrintHeader(out_file, *temp, information);
     307
     308  for (iter=temp->Iterators().Local().Begin(); iter!=temp->Iterators().Local().End(); ++iter)
     309    out << std::scientific << temp->GetVal(*iter) << std::endl;
     310
     311  out_file << out.str();
     312  out_file.close();
     313}
     314
     315void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& grid, const char* information)
    289316{
    290317  char path_str[129];
    … …  
    294321
    295322  out.open(path_str, std::ios::trunc);
    296 
    297   if (!out.good()) {
    298 #ifdef DEBUG_OUTPUT
    299     printf("Multigrid: File %s not accessible.\n", path_str);
    300 #endif /* DEBUG_OUTPUT */
    301     return;
    302   }
    303323
    304324  out << "# vtk DataFile Version 2.0" << std::endl
    305325      << count << ": " << information << std::endl
    306326      << "ASCII" << std::endl
    307       << "DATASET STRUCTURED_POINTS" << std::endl;
    308 
    309   if (inner) {
    310     out << "DIMENSIONS "
    311         << mesh.Local().Size().X() << " "
    312         << mesh.Local().Size().Y() << " "
    313         << mesh.Local().Size().Z() << std::endl;
    314   }else {
    315     out << "DIMENSIONS "
    316         << mesh.Local().SizeTotal().X() << " "
    317         << mesh.Local().SizeTotal().Y() << " "
    318         << mesh.Local().SizeTotal().Z() << std::endl;
    319   }
    320 
    321   out   << "ORIGIN 0 0 0" << std::endl
    322         << "SPACING " << mesh.MeshWidth() << " "
    323         << mesh.MeshWidth() << " "
    324         << mesh.MeshWidth() << std::endl;
    325 
    326   if (inner)
    327     out << "POINT_DATA " << mesh.Local().Size().Product() << std::endl;
    328   else
    329     out << "POINT_DATA " << mesh.Local().SizeTotal().Product() << std::endl;
    330 
    331   out << "SCALARS residual double 1" << std::endl
     327      << "DATASET STRUCTURED_POINTS" << std::endl
     328      << grid.Local().Size().X() << " "
     329      << grid.Local().Size().Y() << " "
     330      << grid.Local().Size().Z() << std::endl
     331      << "ORIGIN 0 0 0" << std::endl
     332      << "SPACING " << grid.Extent().MeshWidth().X() << " "
     333      << grid.Extent().MeshWidth() << " "
     334      << grid.Extent().MeshWidth() << std::endl
     335      << "POINT_DATA " << grid.Local().Size().Product() << std::endl
     336      << "SCALARS residual double 1" << std::endl
    332337      << "LOOKUP_TABLE default" << std::endl;
    333 }
    334 
    335 void CommSerial::PrintDefect(Grid& sol, Grid& rhs, const char* information)
    336 {
    337   std::ofstream out_file;
    338   std::stringstream out;
    339 
    340   TempGrid *temp = MG::GetTempGrid();
    341   temp->SetProperties(sol);
    342   temp->ImportFromResidual(sol, rhs);
    343 
    344 #ifdef DEBUG_MATRIX_CHECKS
    345   temp->IsConsistent();
    346 #endif
    347 
    348   OpenFileAndPrintHeader(out_file, *temp, information, true);
    349 
    350   assert(out_file.good());
    351 
    352   if (!out_file.good())
    353     return;
    354 
    355   for (int k=temp->Local().Begin().Z(); k<temp->Local().End().Z(); k++)
    356     for (int j=temp->Local().Begin().Y(); j<temp->Local().End().Y(); j++)
    357       for (int i=temp->Local().Begin().X(); i<temp->Local().End().X(); i++)
    358         out << std::scientific << temp->GetVal(i,j,k) << std::endl;
    359 
    360   out_file << out.str();
    361 
    362   out_file.close();
    363 }
    364 
    365 void CommSerial::PrintGrid(Grid& mesh, const char* information)
    366 {
    367   std::stringstream out;
    368   std::ofstream out_file;
    369 
    370   CommToGhosts(mesh);
    371 
    372   OpenFileAndPrintHeader(out_file, mesh, information, false);
    373 
    374   for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    375     for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
    376       for (int i=0; i<mesh.Local().SizeTotal().X(); i++)
    377         out << std::scientific << mesh.GetVal(i,j,k) << std::endl;
    378 
    379   out_file << out.str();
    380   out_file.close();
    381 }
    382 
    383 void CommSerial::PrintCorrectedGrid(Grid& mesh, const char* information)
    384 {
    385   std::ofstream out_file;
    386   std::stringstream out;
    387 
    388   CommToGhosts(mesh);
    389 
    390   OpenFileAndPrintHeader(out_file, mesh, information, false);
    391 
    392   for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    393     for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
    394       for (int i=0; i<mesh.Local().SizeTotal().X(); i++)
    395         out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;
    396 
    397   out_file << out.str();
    398 
    399   out_file.close();
    400 }
    401 
    402 void CommSerial::PrintInnerGrid(Grid& mesh, const char* information)
    403 {
    404   std::ofstream out_file;
    405   std::stringstream out;
    406 
    407   CommToGhosts(mesh);
    408 
    409   OpenFileAndPrintHeader(out_file, mesh, information, true);
    410 
    411   for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    412     for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
    413       for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    414         out << std::scientific << mesh.GetVal(i,j,k) << std::endl;
    415 
    416   out_file << out.str();
    417 
    418   out_file.close();
    419 }
    420 
    421 void CommSerial::PrintInnerCorrectedGrid(Grid& mesh, const char* information)
    422 {
    423   std::ofstream out_file;
    424   std::stringstream out;
    425 
    426   CommToGhosts(mesh);
    427 
    428   OpenFileAndPrintHeader(out_file, mesh, information, true);
    429 
    430   for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    431     for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
    432       for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    433         out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;
    434 
    435   out_file << out.str();
    436 
    437   out_file.close();
    438338}
    439339
    … …  
    450350  for (int i=mg->MinLevel(); i<=mg->MaxLevel(); ++i) {
    451351
    452     out << "GRID LEVEL " << i << std::endl
    453 
    454         << "GLOBAL" << std::endl
    455         << "BEGINFINEST: " << (*mg)(i).Global().BeginFinest() << std::endl
    456         << "ENDFINEST: " << (*mg)(i).Global().EndFinest() << std::endl
    457         << "SIZEFINEST: " << (*mg)(i).Global().SizeFinest() << std::endl
    458         << "BEGINLOCAL: " << (*mg)(i).Global().BeginLocal() << std::endl
    459         << "ENDLOCAL: " << (*mg)(i).Global().EndLocal() << std::endl
    460         << "SIZELOCAL: " << (*mg)(i).Global().SizeLocal() << std::endl
    461         << "SIZEGLOBAL: " << (*mg)(i).Global().SizeGlobal() << std::endl
    462         << "BOUNDARYTYPE: " << (*mg)(i).Global().BoundaryType() << std::endl
    463 
    464         << "LOCAL" << std::endl
    465         << "SIZE: " << (*mg)(i).Local().Size() << std::endl
    466         << "SIZETOTAL: " << (*mg)(i).Local().SizeTotal() << std::endl
    467         << "BEGIN: " << (*mg)(i).Local().Begin() << std::endl
    468         << "END: " << (*mg)(i).Local().End() << std::endl
    469         << "HALOBEGIN1: " << (*mg)(i).Local().HaloBegin1() << std::endl
    470         << "HALOEND1: " << (*mg)(i).Local().HaloEnd1() << std::endl
    471         << "HALOBEGIN2: " << (*mg)(i).Local().HaloBegin2() << std::endl
    472         << "HALOEND2: " << (*mg)(i).Local().HaloEnd2() << std::endl
    473         << "BOUNDARYBEGIN1: " << (*mg)(i).Local().BoundaryBegin1() << std::endl
    474         << "BOUNDARYEND1: " << (*mg)(i).Local().BoundaryEnd1() << std::endl
    475         << "BOUNDARYBEGIN2: " << (*mg)(i).Local().BoundaryBegin2() << std::endl
    476         << "BOUNDARYEND2: " << (*mg)(i).Local().BoundaryEnd2() << std::endl
    477         << "ALIGNMENTBEGIN: " << (*mg)(i).Local().AlignmentBegin() << std::endl
    478         << "ALIGNMENTEND: " << (*mg)(i).Local().AlignmentEnd() << std::endl
    479 
    480         << "EXTENT" << std::endl
    481         << "SIZE: " << (*mg)(i).Extent().Size() << std::endl
    482         << "SIZEFACTOR: " << (*mg)(i).Extent().SizeFactor() << std::endl
    483         << "BEGIN: " << (*mg)(i).Extent().Begin() << std::endl
    484         << "END: " << (*mg)(i).Extent().End() << std::endl
    485         << "MESHWIDTH: " << (*mg)(i).Extent().MeshWidth() << std::endl;
     352    out << "###########################################################" << std::endl
     353        << "LEVEL:                " << i << std::endl
     354        << "GLOBAL:" << std::endl
     355        << "  LOCAL_BEGIN:        " << (*mg)(i).Global().LocalBegin() << std::endl
     356        << "  LOCAL_END:          " << (*mg)(i).Global().LocalEnd() << std::endl
     357        << "  LOCAL_SIZE:         " << (*mg)(i).Global().LocalSize() << std::endl
     358        << "  GLOBAL_FINER_BEGIN: " << (*mg)(i).Global().GlobalFinerBegin() << std::endl
     359        << "  GLOBAL_FINER_END:   " << (*mg)(i).Global().GlobalFinerEnd() << std::endl
     360        << "  GLOBAL_FINER_SIZE:  " << (*mg)(i).Global().GlobalFinerSize() << std::endl
     361        << "  LOCAL_FINER_BEGIN:  " << (*mg)(i).Global().LocalFinerBegin() << std::endl
     362        << "  LOCAL_FINER_END:    " << (*mg)(i).Global().LocalFinerEnd() << std::endl
     363        << "  LOCAL_FINER_SIZE:   " << (*mg)(i).Global().LocalFinerSize() << std::endl
     364        << "  FINEST_ABS_BEGIN:   " << (*mg)(i).Global().FinestAbsBegin() << std::endl
     365        << "  FINEST_ABS_END:     " << (*mg)(i).Global().FinestAbsEnd() << std::endl
     366        << "  FINEST_ABS_SIZE:    " << (*mg)(i).Global().FinestAbsSize() << std::endl
     367        << "  GLOBAL_SIZE:        " << (*mg)(i).Global().GlobalSize() << std::endl
     368        << "  BOUNDARY_TYPE:      " << (*mg)(i).Global().BoundaryType() << std::endl
     369        << "LOCAL:" << std::endl
     370        << "  BEGIN:              " << (*mg)(i).Local().Begin() << std::endl
     371        << "  END:                " << (*mg)(i).Local().End() << std::endl
     372        << "  SIZE:               " << (*mg)(i).Local().Size() << std::endl
     373        << "  SIZE_TOTAL:         " << (*mg)(i).Local().SizeTotal() << std::endl
     374        << "  HALO_BEGIN_1:       " << (*mg)(i).Local().HaloBegin1() << std::endl
     375        << "  HALO_END_1:         " << (*mg)(i).Local().HaloEnd1() << std::endl
     376        << "  HALO_SIZE_1:        " << (*mg)(i).Local().HaloSize1() << std::endl
     377        << "  HALO_BEGIN_2:       " << (*mg)(i).Local().HaloBegin2() << std::endl
     378        << "  HALO_END_2:         " << (*mg)(i).Local().HaloEnd2() << std::endl
     379        << "  HALO_SIZE_2:        " << (*mg)(i).Local().HaloSize2() << std::endl
     380        << "  BOUNDARY_BEGIN_1:   " << (*mg)(i).Local().BoundaryBegin1() << std::endl
     381        << "  BOUNDARY_END_1:     " << (*mg)(i).Local().BoundaryEnd1() << std::endl
     382        << "  BOUNDARY_SIZE_1:    " << (*mg)(i).Local().BoundarySize1() << std::endl
     383        << "  BOUNDARY_BEGIN_2:   " << (*mg)(i).Local().BoundaryBegin2() << std::endl
     384        << "  BOUNDARY_END_2:     " << (*mg)(i).Local().BoundaryEnd2() << std::endl
     385        << "  BOUNDARY_SIZE_2:    " << (*mg)(i).Local().BoundarySize2() << std::endl
     386        << "  FINER_BEGIN:        " << (*mg)(i).Local().FinerBeginFoo() << std::endl
     387        << "  FINER_END:          " << (*mg)(i).Local().FinerEndFoo() << std::endl
     388        << "  FINER_SIZE:         " << (*mg)(i).Local().FinerSizeFoo() << std::endl
     389        << "EXTENT:" << std::endl
     390        << "  BEGIN:              " << (*mg)(i).Extent().Begin() << std::endl
     391        << "  END:                " << (*mg)(i).Extent().End() << std::endl
     392        << "  SIZE:               " << (*mg)(i).Extent().Size() << std::endl
     393        << "  MESH_WIDTH:         " << (*mg)(i).Extent().MeshWidth() << std::endl
     394        << "###########################################################" << std::endl;
    486395
    487396  }
    … …  
    489398  assert(out.good());
    490399  out.close();
    491 }
    492 
    493 void CommSerial::DebugPrintError(const Grid& sol, const char* information)
    494 {
    495   Index i;
    496   std::ofstream out_file;
    497   std::stringstream out;
    498 
    499   OpenFileAndPrintHeader(out_file, sol, information, false);
    500 
    501   for (i.Z()=0; i.Z()<sol.Local().SizeTotal().Z(); ++i.Z())
    502     for (i.Y()=0; i.Y()<sol.Local().SizeTotal().Y(); ++i.Y())
    503       for (i.X()=0; i.X()<sol.Local().SizeTotal().X(); ++i.X())
    504         out << std::scientific << sol.DebugKnownSolution(i) - sol.GetVal(i) << std::endl;
    505 
    506   out_file << out.str();
    507 
    508   out_file.close();
    509 }
    510 
    511 void CommSerial::DebugPrintErrorNorm(Grid& sol)
    512 {
    513 #ifdef HAVE_LIBGSL
    514   char path_str[129];
    515   std::ofstream out;
    516   const vmg_float sp = sol.MeshWidth();
    517   const vmg_float sp_3_inv = 1.0/(sp*sp*sp);
    518   const int numsamples = 1e6;
    519   const vmg_float numsamples_inv = 1.0 / numsamples;
    520   vmg_float val1, val2, diff;
    521 
    522   vmg_float norm_L1 = 0.0;
    523   vmg_float norm_L2 = 0.0;
    524   vmg_float norm_Linfty = 0.0;
    525   vmg_float norm_u_L1 = 0.0;
    526   vmg_float norm_u_L2 = 0.0;
    527   vmg_float norm_u_Linfty = 0.0;
    528 
    529   Vector x, x_le, x_ri;
    530   Index ind;
    531 
    532   gsl_qrng *q = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 3);
    533 
    534   sprintf(path_str, "%serror.txt", defects_path_str.c_str());
    535 
    536   out.open(path_str, std::ios::app);
    537 
    538   assert(out.good());
    539 
    540   CommToGhosts(sol);
    541 
    542   for (int i=0; i<numsamples; i++) {
    543 
    544     gsl_qrng_get(q, x.vec());
    545 
    546     ind = x / sp;
    547 
    548     x_le = x - sp * static_cast<Vector>(ind);
    549     x_ri = sp - x_le;
    550 
    551     for (int j=0; j<3; ++j)
    552       if (this->BoundaryConditions()[j] == Periodic)
    553         ind[j] += sol.Local().Begin()[j];
    554 
    555     val1 = ( sol.GetVal(ind) * x_ri[0] * x_ri[1] * x_ri[2] +
    556              sol.GetVal(ind[0]+1, ind[1]  , ind[2]  ) * x_le[0] * x_ri[1] * x_ri[2] +
    557              sol.GetVal(ind[0]  , ind[1]+1, ind[2]  ) * x_ri[0] * x_le[1] * x_ri[2] +
    558              sol.GetVal(ind[0]  , ind[1]  , ind[2]+1) * x_ri[0] * x_ri[1] * x_le[2] +
    559              sol.GetVal(ind[0]+1, ind[1]+1, ind[2]  ) * x_le[0] * x_le[1] * x_ri[2] +
    560              sol.GetVal(ind[0]+1, ind[1]  , ind[2]+1) * x_le[0] * x_ri[1] * x_le[2] +
    561              sol.GetVal(ind[0]  , ind[1]+1, ind[2]+1) * x_ri[0] * x_le[1] * x_le[2] +
    562              sol.GetVal(ind[0]+1, ind[1]+1, ind[2]+1) * x_le[0] * x_le[1] * x_le[2] ) * sp_3_inv;
    563 
    564     val2 = sol.DebugKnownSolution(x);
    565 
    566     diff = fabs(val1 - val2);
    567 
    568     norm_L1 += diff * numsamples_inv;
    569     norm_L2 += diff*diff * numsamples_inv;
    570     norm_Linfty = std::max(norm_Linfty, diff);
    571 
    572     norm_u_L1 += fabs(val2) * numsamples_inv;
    573     norm_u_L2 += val2*val2 * numsamples_inv;
    574     norm_u_Linfty = std::max(norm_u_Linfty, fabs(val2));
    575   }
    576 
    577   norm_L2 = sqrt(norm_L2);
    578   norm_u_L2 = sqrt(norm_u_L2);
    579 
    580   std::cout << std::scientific << std::endl
    581             << "Error L1-norm:              " << norm_L1 << std::endl
    582             << "Error L2-norm:              " << norm_L2 << std::endl
    583             << "Error Linfty-norm:          " << norm_Linfty << std::endl
    584             << "Relative error L1-norm:     " << norm_L1 / norm_u_L1 << std::endl
    585             << "Relative error L2-norm:     " << norm_L2 / norm_u_L2 << std::endl
    586             << "Relative error Linfty-norm: " << norm_Linfty / norm_u_Linfty << std::endl;
    587 
    588   gsl_qrng_free(q);
    589 
    590   out << std::scientific << error_norm_count++ << " "
    591       << norm_L1 << " "
    592       << norm_L2 << " "
    593       << norm_Linfty << " "
    594       << norm_L1/norm_u_L1 << " "
    595       << norm_L2/norm_u_L2 << " "
    596       << norm_Linfty/norm_u_Linfty << std::endl;
    597 
    598   out.close();
    599 
    600 #endif
    601400}
    602401
    … …  
    611410void CommSerial::PrintGridStructureLevel(Grid& grid, std::ofstream& out)
    612411{
    613   const vmg_float sp = grid.MeshWidth();
     412  const Vector& sp = grid.Extent().MeshWidth();
    614413  int numLines;
    615414
    … …  
    644443      for (int j=0; j<grid.Local().Size().Y(); j++)
    645444        for (int k=0; k<grid.Local().Size().Z(); k++)
    646           out << grid.Extent().Begin().X() + i * sp << " "
    647               << grid.Extent().Begin().Y() + j * sp << " "
    648               << grid.Extent().Begin().Z() + k * sp << " ";
     445          out << grid.Extent().Begin().X() + i * sp.X() << " "
     446              << grid.Extent().Begin().Y() + j * sp.Y() << " "
     447              << grid.Extent().Begin().Z() + k * sp.Z() << " ";
    649448  }else {
    650449    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
    651450      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
    652451        for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
    653           out << grid.Extent().Begin().X() + i * sp << " "
    654               << grid.Extent().Begin().Y() + j * sp << " "
    655               << grid.Extent().Begin().Z() + k * sp << " ";
     452          out << grid.Extent().Begin().X() + i * sp.X() << " "
     453              << grid.Extent().Begin().Y() + j * sp.Y() << " "
     454              << grid.Extent().Begin().Z() + k * sp.Z() << " ";
    656455  }
    657456
Note: See TracChangeset for help on using the changeset viewer.