Ignore:
Timestamp:
Nov 22, 2011, 9:22:10 PM (14 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
facba0
Parents:
66f24d
Message:

Major vmg update.

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/comm/comm_serial.cpp

    r66f24d rdfed1c  
    1717#endif
    1818
    19 #include <iostream>
    20 #include <fstream>
    21 #include <cassert>
    22 #include <cmath>
     19#include <cstdarg>
     20#include <cstdio>
     21#include <sstream>
    2322
    2423#include "base/discretization.hpp"
     24#include "base/helper.hpp"
     25#include "base/index.hpp"
     26#include "base/linked_cell_list.hpp"
    2527#include "base/stencil.hpp"
    2628#include "base/vector.hpp"
    … …  
    3335using namespace VMG;
    3436
    35 inline vmg_float pow_3(vmg_float val)
    36 {
    37   return val*val*val;
    38 }
    39 
    40 void CommSerial::Init()
    41 {
    42   output_count = 0;
     37static char print_buffer[512];
     38
     39void CommSerial::InitCommSerial()
     40{
    4341  error_norm_count = 0;
    4442  residual_count = 0;
    45 
    46   /* Create directory to output defects */
    47 #ifdef HAVE_BOOST_FILESYSTEM
    48   char buffer[129];
    49   time_t rawtime;
    50   struct tm *timeinfo;
    51 
    52   time(&rawtime);
    53   timeinfo = localtime(&rawtime);
    54   strftime(buffer, 128, "output/%Y_%m_%d_%H_%M_%S/", timeinfo);
    55   defects_path_str = buffer;
    56 #else
    57   defects_path_str = "";
    58 #endif
    59 }
    60 
    61 void CommSerial::CreateOutputDirectory()
    62 {
    63 #ifdef HAVE_BOOST_FILESYSTEM
    64   if (!fs::exists("output"))
    65     fs::create_directory("output");
    66 
    67   if (!fs::exists(defects_path_str.c_str()))
    68     fs::create_directory(defects_path_str.c_str());
    69 #endif
     43}
     44
     45void CommSerial::CommSubgrid(Grid& grid_old, Grid& grid_new)
     46{
     47  Grid::iterator iter;
     48
     49  if (&grid_old == &grid_new)
     50    return;
     51
     52  for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
     53    grid_new(*iter) = grid_old.GetVal(*iter);
     54}
     55
     56void CommSerial::CommAddSubgrid(Grid& grid_old, Grid& grid_new)
     57{
     58  Grid::iterator iter;
     59
     60  if (&grid_old == &grid_new)
     61    return;
     62
     63  for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
     64    grid_new(*iter) += grid_old.GetVal(*iter);
    7065}
    7166
    … …  
    8075}
    8176
    82 vmg_float CommSerial::ComputeResidualNorm(Multigrid& sol, Multigrid& rhs)
    83 {
    84 #ifdef DEBUG
    85   sol().IsCompatible(rhs());
    86   sol().IsConsistent();
    87   rhs().IsConsistent();
    88 #endif
    89 
    90   Stencil::iterator iter;
    91   Index i;
    92   vmg_float val;
    93   vmg_float norm = 0.0;
    94 
    95   const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol());
    96   const Stencil& A = MG::GetDiscretization()->GetStencil();
    97 
    98   this->CommToGhosts(sol());
    99 
    100   if (sol().Global().BoundaryType() == LocallyRefined)
    101     MG::GetDiscretization()->SetInnerBoundary(sol(), rhs(), sol(sol.Level()-1));
    102 
    103   for (i.X()=rhs().Local().Begin().X(); i.X()<rhs().Local().End().X(); ++i.X())
    104     for (i.Y()=rhs().Local().Begin().Y(); i.Y()<rhs().Local().End().Y(); ++i.Y())
    105       for (i.Z()=rhs().Local().Begin().Z(); i.Z()<rhs().Local().End().Z(); ++i.Z()) {
    106 
    107         val = rhs().GetVal(i) - prefactor * A.GetDiag() * sol().GetVal(i);
    108 
    109         for (iter = A.begin(); iter != A.end(); iter++)
    110           val -= prefactor * iter->val * sol().GetVal(i+iter);
    111 
    112         norm += val*val;
    113 
    114       }
    115 
    116   norm = sqrt( pow_3(sol().MeshWidth()) * norm );
    117 
    118 #ifdef DEBUG
    119   std::ofstream out;
    120   char path_str[129];
    121 
    122   CreateOutputDirectory();
    123 
    124   sprintf(path_str, "%sresidual.txt", defects_path_str.c_str());
    125 
    126   out.open(path_str, std::ios::app);
    127 
    128   assert(out.good());
    129 
    130   out << residual_count++ << " " << norm << std::endl;
    131 
    132   out.close();
    133 #endif
    134 
    135   return norm;
    136 }
    137 
    13877void CommSerial::CommFromGhosts(Grid& mesh)
    13978{
    140   if (this->BoundaryConditions() == Periodic) {
    141 
    142     // Copy values in x-direction
    143     for (int i=mesh.Local().HaloBegin1().X(); i<mesh.Local().HaloEnd1().X(); i++)
    144       for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
    145         for (int k=0; k<mesh.Local().SizeTotal().Z(); k++) {
    146           assert(i+mesh.Local().Size().X() >= mesh.Local().Begin().X() && i+mesh.Local().Size().X() < mesh.Local().End().X());
    147           mesh(i+mesh.Local().Size().X(),j,k) += mesh.GetVal(i,j,k);
    148         }
    149 
    150     for (int i=mesh.Local().HaloBegin2().X(); i<mesh.Local().HaloEnd2().X(); i++)
    151       for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
    152         for (int k=0; k<mesh.Local().SizeTotal().Z(); k++) {
    153           assert(i-mesh.Local().Size().X() >= mesh.Local().Begin().X() && i-mesh.Local().Size().X() < mesh.Local().End().X());
    154           mesh(i-mesh.Local().Size().X(),j,k) += mesh.GetVal(i,j,k);
    155         }
    156 
    157     // Copy values in y-direction
    158     for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    159       for (int j=mesh.Local().HaloBegin1().Y(); j<mesh.Local().HaloEnd1().Y(); j++)
    160         for (int k=0; k<mesh.Local().SizeTotal().Z(); k++) {
    161 
    162           mesh(i,j+mesh.Local().Size().Y(),k) += mesh.GetVal(i,j,k);
    163         }
    164     for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    165       for (int j=mesh.Local().HaloBegin2().Y(); j<mesh.Local().HaloEnd2().Y(); j++)
    166         for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    167           mesh(i,j-mesh.Local().Size().Y(),k) += mesh.GetVal(i,j,k);
    168 
    169     // Copy values in z-direction
    170     for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    171       for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
    172         for (int k=mesh.Local().HaloBegin1().Z(); k<mesh.Local().HaloEnd1().Z(); k++)
    173           mesh(i,j,k+mesh.Local().Size().Z()) += mesh.GetVal(i,j,k);
    174 
    175     for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    176       for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
    177         for (int k=mesh.Local().HaloBegin2().Z(); k<mesh.Local().HaloEnd2().Z(); k++)
    178           mesh(i,j,k-mesh.Local().Size().Z()) += mesh.GetVal(i,j,k);
    179   }
    180 
    181 #ifdef DEBUG
     79  Grid::iterator iter;
     80
     81  const GridIteratorSuite& iterators = mesh.Iterators();
     82  const Index& size = mesh.Local().Size();
     83
     84  if (BoundaryConditions().X() == Periodic) {
     85
     86    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);
     88
     89    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);
     91  }
     92
     93  if (BoundaryConditions().Y() == Periodic) {
     94
     95    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);
     97
     98    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);
     100
     101  }
     102
     103  if (BoundaryConditions().Z() == Periodic) {
     104
     105    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);
     107
     108    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);
     110
     111  }
     112
     113#ifdef DEBUG_MATRIX_CHECKS
    182114  mesh.IsConsistent();
    183115#endif
    … …  
    186118void CommSerial::CommToGhosts(Grid& mesh)
    187119{
    188 #ifdef DEBUG
     120  Grid::iterator iter;
     121
     122  const GridIteratorSuite& iterators = mesh.Iterators();
     123  const Index& size = mesh.Local().Size();
     124
     125  if (BoundaryConditions().Z() == Periodic) {
     126
     127    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());
     129
     130    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());
     132
     133  }
     134
     135  if (BoundaryConditions().Y() == Periodic) {
     136
     137    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());
     139
     140    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());
     142
     143  }
     144
     145  if (BoundaryConditions().X() == Periodic) {
     146
     147    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());
     149
     150    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());
     152
     153  }
     154
     155#ifdef DEBUG_MATRIX_CHECKS
    189156  mesh.IsConsistent();
    190157#endif
    191 
    192   int i,j,k;
    193 
    194   if (this->BoundaryConditions() == Periodic) {
    195 
    196     // Copy values in z-direction
    197     for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    198       for (j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
    199         for (k=mesh.Local().HaloBegin1().Z(); k<mesh.Local().HaloEnd1().Z(); k++)
    200           mesh(i,j,k) = mesh.GetVal(i,j,k+mesh.Local().Size().Z());
    201 
    202     for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    203       for (j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
    204         for (k=mesh.Local().HaloBegin2().Z(); k<mesh.Local().HaloEnd2().Z(); k++)
    205           mesh(i,j,k) = mesh.GetVal(i,j,k-mesh.Local().Size().Z());
    206 
    207     // Copy values in y-direction
    208     for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    209       for (j=mesh.Local().HaloBegin1().Y(); j<mesh.Local().HaloEnd1().Y(); j++)
    210         for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
    211           mesh(i,j,k) = mesh.GetVal(i,j+mesh.Local().Size().Y(),k);
    212 
    213     for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
    214       for (j=mesh.Local().HaloBegin2().Y(); j<mesh.Local().HaloEnd2().Y(); j++)
    215         for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
    216           mesh(i,j,k) = mesh.GetVal(i,j-mesh.Local().Size().Y(),k);
    217 
    218     // Copy values in x-direction
    219     for (i=mesh.Local().HaloBegin1().X(); i<mesh.Local().HaloEnd1().X(); i++)
    220       for (j=0; j<mesh.Local().SizeTotal().Y(); j++)
    221         for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
    222           mesh(i,j,k) = mesh.GetVal(i+mesh.Local().Size().X(),j,k);
    223 
    224     for (i=mesh.Local().HaloBegin2().X(); i<mesh.Local().HaloEnd2().X(); i++)
    225       for (j=0; j<mesh.Local().SizeTotal().Y(); j++)
    226         for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
    227           mesh(i,j,k) = mesh.GetVal(i-mesh.Local().Size().X(),j,k);
    228   }
    229 
    230 #ifdef DEBUG
    231   mesh.IsConsistent();
    232 #endif
     158}
     159
     160Grid& CommSerial::GetGlobalCoarseGrid(Multigrid& multigrid)
     161{
     162  return multigrid(multigrid.MinLevel());
     163}
     164
     165void CommSerial::CommParticles(const Grid& grid, std::list<Particle::Particle>& particles)
     166{
     167  Factory& factory = MG::GetFactory();
     168  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
     169  vmg_float* x = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POS_ARRAY");
     170  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
     171
     172  particles.clear();
     173
     174  for (vmg_int i=0; i<num_particles_local; ++i)
     175    particles.push_back(Particle::Particle(&x[3*i], q[i], 0.0, 0.0, 0, i));
     176}
     177
     178void CommSerial::CommLCListGhosts(const Grid& grid, Particle::LinkedCellList& lc)
     179{
     180  std::list<Particle::Particle>::iterator p_iter;
     181  Grid::iterator g_iter;
     182
     183  for (int i=2; i>=0; --i) {
     184
     185    if (BoundaryConditions()[i] == Periodic) {
     186
     187      for (g_iter=lc.Iterators().Boundary1()[i].Begin(); g_iter!=lc.Iterators().Boundary1()[i].End(); ++g_iter)
     188        lc(*g_iter).clear();
     189
     190      for (g_iter=lc.Iterators().NearBoundary1()[i].Begin(); g_iter!=lc.Iterators().NearBoundary1()[i].End(); ++g_iter)
     191        for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
     192          lc.AddParticleToComplete(p_iter->Pos() + grid.Extent().Size(), p_iter->Charge(),
     193                                   p_iter->Pot(), p_iter->Force(),
     194                                   p_iter->Rank(), p_iter->Index());
     195
     196      for (g_iter=lc.Iterators().Boundary2()[i].Begin(); g_iter!=lc.Iterators().Boundary2()[i].End(); ++g_iter)
     197        lc(*g_iter).clear();
     198
     199      for (g_iter=lc.Iterators().NearBoundary2()[i].Begin(); g_iter!=lc.Iterators().NearBoundary2()[i].End(); ++g_iter)
     200        for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
     201          lc.AddParticleToComplete(p_iter->Pos() - grid.Extent().Size(), p_iter->Charge(),
     202                                   p_iter->Pot(), p_iter->Force(),
     203                                   p_iter->Rank(), p_iter->Index());
     204
     205    }
     206  }
     207}
     208
     209void CommSerial::CommAddPotential(std::list<Particle::Particle>& particles)
     210{
     211  std::list<Particle::Particle>::iterator iter;
     212  vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
     213
     214  for (iter=particles.begin(); iter!=particles.end(); ++iter)
     215    p[iter->Index()] += iter->Pot();
     216}
     217
     218void CommSerial::CommAddPotential(Particle::LinkedCellList& lc)
     219{
     220  Grid::iterator lc_iter;
     221  std::list<Particle::Particle>::iterator p_iter;
     222  vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
     223
     224  for (lc_iter=lc.Iterators().Local().Begin(); lc_iter!=lc.Iterators().Local().End(); ++lc_iter)
     225    for (p_iter=lc(*lc_iter).begin(); p_iter!=lc(*lc_iter).end(); ++p_iter)
     226      p[p_iter->Index()] += p_iter->Pot();
     227}
     228
     229void CommSerial::PrintString(const char* format, ...)
     230{
     231  va_list args;
     232  va_start(args, format);
     233  vsprintf(print_buffer, format, args);
     234  printf("VMG: %s\n", print_buffer);
     235  va_end(args);
     236}
     237
     238void CommSerial::PrintStringOnce(const char* format, ...)
     239{
     240  va_list args;
     241  va_start(args, format);
     242  vsprintf(print_buffer, format, args);
     243  printf("VMG: %s\n", print_buffer);
     244  va_end(args);
     245}
     246
     247void CommSerial::PrintXML(const std::string& filename, const std::string& xml_data)
     248{
     249  pugi::xml_document doc;
     250  doc.load(xml_data.c_str());
     251
     252  pugi::xml_node node_global = doc.prepend_child("Global");
     253  pugi::xml_node node_num_processes = node_global.append_child("NumProcesses");
     254  pugi::xml_node node_num_processes_data = node_num_processes.append_child(pugi::node_pcdata);
     255  node_num_processes_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
     256
     257  doc.save_file(filename.c_str());
     258}
     259
     260void CommSerial::PrintXMLAll(const std::string& filename, const std::string& xml_data)
     261{
     262  PrintXML(filename, xml_data);
    233263}
    234264
    235265void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& mesh, const char* information, bool inner)
    236266{
    237   char file_str[129], path_str[129];
    238   struct tm *timeinfo;
    239   time_t rawtime;
    240 
    241   CreateOutputDirectory();
    242 
    243   output_count++;
    244 
    245   time(&rawtime);
    246   timeinfo = localtime(&rawtime);
    247 
    248   strftime(file_str, 128, "%H_%M_%S", timeinfo);
    249 
    250   sprintf(path_str, "%s%04d_%s.dat", defects_path_str.c_str(), output_count, file_str);
     267  char path_str[129];
     268  int count = OutputCount();
     269
     270  sprintf(path_str, "%s%04d.dat", OutputPath().c_str(), count);
    251271
    252272  out.open(path_str, std::ios::trunc);
    253273
    254274  if (!out.good()) {
     275#ifdef DEBUG_OUTPUT
    255276    printf("Multigrid: File %s not accessible.\n", path_str);
     277#endif /* DEBUG_OUTPUT */
    256278    return;
    257279  }
    258280
    259281  out << "# vtk DataFile Version 2.0" << std::endl
    260       << output_count << ": " << information << std::endl
     282      << count << ": " << information << std::endl
    261283      << "ASCII" << std::endl
    262284      << "DATASET STRUCTURED_POINTS" << std::endl;
    … …  
    288310}
    289311
    290 void CommSerial::PrintDefect(const Grid& sol, const Grid& rhs, const char* information)
    291 {
    292   std::ofstream out;
     312void CommSerial::PrintDefect(Grid& sol, Grid& rhs, const char* information)
     313{
     314  std::ofstream out_file;
     315  std::stringstream out;
    293316
    294317  TempGrid *temp = MG::GetTempGrid();
    … …  
    296319  temp->ImportFromResidual(sol, rhs);
    297320
    298 #ifdef DEBUG
     321#ifdef DEBUG_MATRIX_CHECKS
    299322  temp->IsConsistent();
    300323#endif
    301324
    302   OpenFileAndPrintHeader(out, *temp, information, true);
    303 
    304   assert(out.good());
    305 
    306   if (!out.good())
     325  OpenFileAndPrintHeader(out_file, *temp, information, true);
     326
     327  assert(out_file.good());
     328
     329  if (!out_file.good())
    307330    return;
    308331
    … …  
    312335        out << std::scientific << temp->GetVal(i,j,k) << std::endl;
    313336
    314   out.close();
    315 }
    316 
    317 void CommSerial::PrintGrid(const Grid& mesh, const char* information)
    318 {
    319   std::ofstream out;
    320 
    321   OpenFileAndPrintHeader(out, mesh, information, false);
     337  out_file << out.str();
     338
     339  out_file.close();
     340}
     341
     342void CommSerial::PrintGrid(Grid& mesh, const char* information)
     343{
     344  std::stringstream out;
     345  std::ofstream out_file;
     346
     347  CommToGhosts(mesh);
     348
     349  OpenFileAndPrintHeader(out_file, mesh, information, false);
    322350
    323351  for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    … …  
    326354        out << std::scientific << mesh.GetVal(i,j,k) << std::endl;
    327355
    328   out.close();
    329 }
    330 
    331 void CommSerial::PrintCorrectedGrid(const Grid& mesh, const char* information)
    332 {
    333   std::ofstream out;
    334 
    335   OpenFileAndPrintHeader(out, mesh, information, false);
     356  out_file << out.str();
     357  out_file.close();
     358}
     359
     360void CommSerial::PrintCorrectedGrid(Grid& mesh, const char* information)
     361{
     362  std::ofstream out_file;
     363  std::stringstream out;
     364
     365  CommToGhosts(mesh);
     366
     367  OpenFileAndPrintHeader(out_file, mesh, information, false);
    336368
    337369  for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    … …  
    340372        out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;
    341373
    342   out.close();
    343 }
    344 
    345 void CommSerial::PrintInnerGrid(const Grid& mesh, const char* information)
    346 {
    347   std::ofstream out;
    348 
    349   OpenFileAndPrintHeader(out, mesh, information, true);
     374  out_file << out.str();
     375
     376  out_file.close();
     377}
     378
     379void CommSerial::PrintInnerGrid(Grid& mesh, const char* information)
     380{
     381  std::ofstream out_file;
     382  std::stringstream out;
     383
     384  CommToGhosts(mesh);
     385
     386  OpenFileAndPrintHeader(out_file, mesh, information, true);
    350387
    351388  for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    … …  
    354391        out << std::scientific << mesh.GetVal(i,j,k) << std::endl;
    355392
    356   out.close();
    357 }
    358 
    359 void CommSerial::PrintInnerCorrectedGrid(const Grid& mesh, const char* information)
    360 {
    361   std::ofstream out;
    362 
    363   OpenFileAndPrintHeader(out, mesh, information, true);
     393  out_file << out.str();
     394
     395  out_file.close();
     396}
     397
     398void CommSerial::PrintInnerCorrectedGrid(Grid& mesh, const char* information)
     399{
     400  std::ofstream out_file;
     401  std::stringstream out;
     402
     403  CommToGhosts(mesh);
     404
     405  OpenFileAndPrintHeader(out_file, mesh, information, true);
    364406
    365407  for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    … …  
    368410        out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;
    369411
     412  out_file << out.str();
     413
     414  out_file.close();
     415}
     416
     417void CommSerial::PrintAllSettings()
     418{
     419  Multigrid* mg = MG::GetFactory().Get("SOL")->Cast<Multigrid>();
     420  std::stringstream buf;
     421  std::ofstream out;
     422
     423  buf << OutputPath() << "settings.txt";
     424
     425  out.open(buf.str().c_str(), std::ios::trunc);
     426
     427  for (int i=mg->MinLevel(); i<=mg->MaxLevel(); ++i) {
     428
     429    out << "GRID LEVEL " << i << std::endl
     430
     431        << "GLOBAL" << std::endl
     432        << "BEGINFINEST: " << (*mg)(i).Global().BeginFinest() << std::endl
     433        << "ENDFINEST: " << (*mg)(i).Global().EndFinest() << std::endl
     434        << "SIZEFINEST: " << (*mg)(i).Global().SizeFinest() << std::endl
     435        << "BEGINLOCAL: " << (*mg)(i).Global().BeginLocal() << std::endl
     436        << "ENDLOCAL: " << (*mg)(i).Global().EndLocal() << std::endl
     437        << "SIZELOCAL: " << (*mg)(i).Global().SizeLocal() << std::endl
     438        << "SIZEGLOBAL: " << (*mg)(i).Global().SizeGlobal() << std::endl
     439        << "BOUNDARYTYPE: " << (*mg)(i).Global().BoundaryType() << std::endl
     440
     441        << "LOCAL" << std::endl
     442        << "SIZE: " << (*mg)(i).Local().Size() << std::endl
     443        << "SIZETOTAL: " << (*mg)(i).Local().SizeTotal() << std::endl
     444        << "BEGIN: " << (*mg)(i).Local().Begin() << std::endl
     445        << "END: " << (*mg)(i).Local().End() << std::endl
     446        << "HALOBEGIN1: " << (*mg)(i).Local().HaloBegin1() << std::endl
     447        << "HALOEND1: " << (*mg)(i).Local().HaloEnd1() << std::endl
     448        << "HALOBEGIN2: " << (*mg)(i).Local().HaloBegin2() << std::endl
     449        << "HALOEND2: " << (*mg)(i).Local().HaloEnd2() << std::endl
     450        << "BOUNDARYBEGIN1: " << (*mg)(i).Local().BoundaryBegin1() << std::endl
     451        << "BOUNDARYEND1: " << (*mg)(i).Local().BoundaryEnd1() << std::endl
     452        << "BOUNDARYBEGIN2: " << (*mg)(i).Local().BoundaryBegin2() << std::endl
     453        << "BOUNDARYEND2: " << (*mg)(i).Local().BoundaryEnd2() << std::endl
     454        << "ALIGNMENTBEGIN: " << (*mg)(i).Local().AlignmentBegin() << std::endl
     455        << "ALIGNMENTEND: " << (*mg)(i).Local().AlignmentEnd() << std::endl
     456
     457        << "EXTENT" << std::endl
     458        << "SIZE: " << (*mg)(i).Extent().Size() << std::endl
     459        << "SIZEFACTOR: " << (*mg)(i).Extent().SizeFactor() << std::endl
     460        << "BEGIN: " << (*mg)(i).Extent().Begin() << std::endl
     461        << "END: " << (*mg)(i).Extent().End() << std::endl
     462        << "MESHWIDTH: " << (*mg)(i).Extent().MeshWidth() << std::endl;
     463
     464  }
     465
     466  assert(out.good());
    370467  out.close();
    371468}
    … …  
    374471{
    375472  Index i;
    376   std::ofstream out;
    377 
    378   OpenFileAndPrintHeader(out, sol, information, false);
     473  std::ofstream out_file;
     474  std::stringstream out;
     475
     476  OpenFileAndPrintHeader(out_file, sol, information, false);
    379477
    380478  for (i.Z()=0; i.Z()<sol.Local().SizeTotal().Z(); ++i.Z())
    381479    for (i.Y()=0; i.Y()<sol.Local().SizeTotal().Y(); ++i.Y())
    382       for (i.X()=0; i.X()<sol.Local().SizeTotal().X(); ++i.X()) {
     480      for (i.X()=0; i.X()<sol.Local().SizeTotal().X(); ++i.X())
    383481        out << std::scientific << sol.DebugKnownSolution(i) - sol.GetVal(i) << std::endl;
    384482
    385       }
     483  out_file << out.str();
     484
     485  out_file.close();
     486}
     487
     488void CommSerial::DebugPrintErrorNorm(Grid& sol)
     489{
     490#ifdef HAVE_LIBGSL
     491  char path_str[129];
     492  std::ofstream out;
     493  const vmg_float sp = sol.MeshWidth();
     494  const vmg_float sp_3_inv = 1.0/(sp*sp*sp);
     495  const int numsamples = 1e6;
     496  const vmg_float numsamples_inv = 1.0 / numsamples;
     497  vmg_float val1, val2, diff;
     498
     499  vmg_float norm_L1 = 0.0;
     500  vmg_float norm_L2 = 0.0;
     501  vmg_float norm_Linfty = 0.0;
     502  vmg_float norm_u_L1 = 0.0;
     503  vmg_float norm_u_L2 = 0.0;
     504  vmg_float norm_u_Linfty = 0.0;
     505
     506  Vector x, x_le, x_ri;
     507  Index ind;
     508
     509  gsl_qrng *q = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 3);
     510
     511  sprintf(path_str, "%serror.txt", defects_path_str.c_str());
     512
     513  out.open(path_str, std::ios::app);
     514
     515  assert(out.good());
     516
     517  CommToGhosts(sol);
     518
     519  for (int i=0; i<numsamples; i++) {
     520
     521    gsl_qrng_get(q, x.vec());
     522
     523    ind = x / sp;
     524
     525    x_le = x - sp * static_cast<Vector>(ind);
     526    x_ri = sp - x_le;
     527
     528    for (int j=0; j<3; ++j)
     529      if (this->BoundaryConditions()[j] == Periodic)
     530        ind[j] += sol.Local().Begin()[j];
     531
     532    val1 = ( sol.GetVal(ind) * x_ri[0] * x_ri[1] * x_ri[2] +
     533             sol.GetVal(ind[0]+1, ind[1]  , ind[2]  ) * x_le[0] * x_ri[1] * x_ri[2] +
     534             sol.GetVal(ind[0]  , ind[1]+1, ind[2]  ) * x_ri[0] * x_le[1] * x_ri[2] +
     535             sol.GetVal(ind[0]  , ind[1]  , ind[2]+1) * x_ri[0] * x_ri[1] * x_le[2] +
     536             sol.GetVal(ind[0]+1, ind[1]+1, ind[2]  ) * x_le[0] * x_le[1] * x_ri[2] +
     537             sol.GetVal(ind[0]+1, ind[1]  , ind[2]+1) * x_le[0] * x_ri[1] * x_le[2] +
     538             sol.GetVal(ind[0]  , ind[1]+1, ind[2]+1) * x_ri[0] * x_le[1] * x_le[2] +
     539             sol.GetVal(ind[0]+1, ind[1]+1, ind[2]+1) * x_le[0] * x_le[1] * x_le[2] ) * sp_3_inv;
     540
     541    val2 = sol.DebugKnownSolution(x);
     542
     543    diff = fabs(val1 - val2);
     544
     545    norm_L1 += diff * numsamples_inv;
     546    norm_L2 += diff*diff * numsamples_inv;
     547    norm_Linfty = std::max(norm_Linfty, diff);
     548
     549    norm_u_L1 += fabs(val2) * numsamples_inv;
     550    norm_u_L2 += val2*val2 * numsamples_inv;
     551    norm_u_Linfty = std::max(norm_u_Linfty, fabs(val2));
     552  }
     553
     554  norm_L2 = sqrt(norm_L2);
     555  norm_u_L2 = sqrt(norm_u_L2);
     556
     557  std::cout << std::scientific << std::endl
     558            << "Error L1-norm:              " << norm_L1 << std::endl
     559            << "Error L2-norm:              " << norm_L2 << std::endl
     560            << "Error Linfty-norm:          " << norm_Linfty << std::endl
     561            << "Relative error L1-norm:     " << norm_L1 / norm_u_L1 << std::endl
     562            << "Relative error L2-norm:     " << norm_L2 / norm_u_L2 << std::endl
     563            << "Relative error Linfty-norm: " << norm_Linfty / norm_u_Linfty << std::endl;
     564
     565  gsl_qrng_free(q);
     566
     567  out << std::scientific << error_norm_count++ << " "
     568      << norm_L1 << " "
     569      << norm_L2 << " "
     570      << norm_Linfty << " "
     571      << norm_L1/norm_u_L1 << " "
     572      << norm_L2/norm_u_L2 << " "
     573      << norm_Linfty/norm_u_Linfty << std::endl;
     574
     575  out.close();
     576
     577#endif
    386578}
    387579
    … …  
    527719  char path_str[129];
    528720
    529   CreateOutputDirectory();
    530 
    531   sprintf(path_str, "%sgrid.vtp", defects_path_str.c_str());
     721  sprintf(path_str, "%sgrid.vtp", OutputPath().c_str());
    532722
    533723  out.open(path_str, std::ios::trunc);
    534724
    535725  if (!out.good()) {
     726#ifdef DEBUG_OUTPUT
    536727    printf("Multigrid: File %s not accessible.\n", path_str);
     728#endif /* DEBUG_OUTPUT */
    537729    return;
    538730  }
    … …  
    550742  out.close();
    551743}
     744
     745std::string CommSerial::CreateOutputDirectory()
     746{
     747#ifdef HAVE_BOOST_FILESYSTEM
     748  std::string path, unique_path;
     749  std::stringstream unique_suffix;
     750  int suffix_counter = 0;
     751  char buffer[129];
     752  time_t rawtime;
     753  struct tm *timeinfo;
     754
     755  time(&rawtime);
     756  timeinfo = localtime(&rawtime);
     757  strftime(buffer, 128, "./output/%Y_%m_%d_%H_%M_%S/", timeinfo);
     758  path = buffer;
     759
     760  if (!fs::exists("output"))
     761    fs::create_directory("output");
     762
     763  unique_path = path;
     764
     765  while (fs::exists(unique_path.c_str())) {
     766
     767    unique_suffix.str("");
     768    unique_suffix << "_" << suffix_counter++ << "/";
     769
     770    unique_path = path;
     771    unique_path.replace(unique_path.size()-1, 1, unique_suffix.str());
     772
     773  }
     774
     775  fs::create_directory(unique_path.c_str());
     776
     777  return unique_path;
     778
     779#else
     780
     781  return "./";
     782
     783#endif
     784}
Note: See TracChangeset for help on using the changeset viewer.