source: src/comm/comm_serial.cpp@ 6f05224

Last change on this file since 6f05224 was f003a9, checked in by Julian Iseringhausen <isering@…>, 14 years ago

Refactored vmg in order to separate the core library and the particle simulation part properly.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1798 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 17.4 KB
Line 
1/**
2 * @file comm_serial.cpp
3 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
4 * @date Mon Apr 18 12:28:12 2011
5 *
6 * @brief VMG::CommSerial
7 *
8 */
9
10#ifdef HAVE_CONFIG_H
11#include <config.h>
12#endif
13
14#ifdef HAVE_BOOST_FILESYSTEM
15#include <boost/filesystem.hpp>
16namespace fs = boost::filesystem;
17#endif
18
19#include <cstdarg>
20#include <cstdio>
21#include <sstream>
22
23#include "base/discretization.hpp"
24#include "base/helper.hpp"
25#include "base/index.hpp"
26#include "base/stencil.hpp"
27#include "base/vector.hpp"
28#include "comm/comm_serial.hpp"
29#include "grid/multigrid.hpp"
30#include "grid/tempgrid.hpp"
31#include "mg.hpp"
32
33using namespace VMG;
34
35static char print_buffer[512];
36
37void CommSerial::InitCommSerial()
38{
39 error_norm_count = 0;
40 residual_count = 0;
41}
42
43void CommSerial::CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
44{
45 Grid::iterator iter;
46
47 if (&grid_old == &grid_new)
48 return;
49
50 for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
51 grid_new(*iter) = grid_old.GetVal(*iter);
52}
53
54void CommSerial::CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
55{
56 Grid::iterator iter;
57
58 if (&grid_old == &grid_new)
59 return;
60
61 for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
62 grid_new(*iter) += grid_old.GetVal(*iter);
63}
64
65Grid& CommSerial::GetFinerGrid(Multigrid& multigrid)
66{
67 return multigrid(multigrid.Level()+1);
68}
69
70Grid& CommSerial::GetCoarserGrid(Multigrid& multigrid)
71{
72 return multigrid(multigrid.Level()-1);
73}
74
75void CommSerial::CommFromGhosts(Grid& grid)
76{
77 Grid::iterator iter;
78
79 const GridIteratorSuite& iterators = grid.Iterators();
80 const Index& size = grid.Local().Size();
81
82 if (BoundaryConditions().X() == Periodic) {
83
84 for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
85 grid((*iter).X()+size.X(), (*iter).Y(), (*iter).Z()) += grid.GetVal(*iter);
86
87 for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
88 grid((*iter).X()-size.X(), (*iter).Y(), (*iter).Z()) += grid.GetVal(*iter);
89 }
90
91 if (BoundaryConditions().Y() == Periodic) {
92
93 for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
94 grid((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z()) += grid.GetVal(*iter);
95
96 for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
97 grid((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z()) += grid.GetVal(*iter);
98
99 }
100
101 if (BoundaryConditions().Z() == Periodic) {
102
103 for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
104 grid((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z()) += grid.GetVal(*iter);
105
106 for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
107 grid((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z()) += grid.GetVal(*iter);
108
109 }
110
111#ifdef DEBUG_MATRIX_CHECKS
112 grid.IsConsistent();
113#endif
114}
115
116void CommSerial::CommToGhosts(Grid& grid)
117{
118 Grid::iterator iter;
119
120 const GridIteratorSuite& iterators = grid.Iterators();
121 const Index& size = grid.Local().Size();
122
123 if (BoundaryConditions().Z() == Periodic) {
124
125 for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
126 grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z());
127
128 for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
129 grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z());
130
131 }
132
133 if (BoundaryConditions().Y() == Periodic) {
134
135 for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
136 grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z());
137
138 for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
139 grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z());
140
141 }
142
143 if (BoundaryConditions().X() == Periodic) {
144
145 for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
146 grid(*iter) = grid.GetVal((*iter).X()+size.X(), (*iter).Y(), (*iter).Z());
147
148 for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
149 grid(*iter) = grid.GetVal((*iter).X()-size.X(), (*iter).Y(), (*iter).Z());
150
151 }
152
153#ifdef DEBUG_MATRIX_CHECKS
154 grid.IsConsistent();
155#endif
156}
157
158void CommSerial::CommToGhostsAsyncStart(Grid& grid)
159{
160 CommToGhosts(grid);
161}
162
163void CommSerial::CommToGhostsAsyncFinish(Grid& grid)
164{
165}
166
167void CommSerial::CommFromGhostsAsyncStart(Grid& grid)
168{
169 CommFromGhosts(grid);
170}
171
172void CommSerial::CommFromGhostsAsyncFinish(Grid& grid)
173{
174}
175
176void CommSerial::PrintString(const char* format, ...)
177{
178 va_list args;
179 va_start(args, format);
180 vsprintf(print_buffer, format, args);
181 printf("VMG: %s\n", print_buffer);
182 va_end(args);
183}
184
185void CommSerial::PrintStringOnce(const char* format, ...)
186{
187 va_list args;
188 va_start(args, format);
189 vsprintf(print_buffer, format, args);
190 printf("VMG: %s\n", print_buffer);
191 va_end(args);
192}
193
194void CommSerial::PrintXML(const std::string& filename, const std::string& xml_data)
195{
196 pugi::xml_document doc;
197 doc.load(xml_data.c_str());
198
199 pugi::xml_node node_global = doc.prepend_child("Global");
200 pugi::xml_node node_num_processes = node_global.append_child("NumProcesses");
201 pugi::xml_node node_num_processes_data = node_num_processes.append_child(pugi::node_pcdata);
202 node_num_processes_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
203
204 doc.save_file(filename.c_str());
205}
206
207void CommSerial::PrintXMLAll(const std::string& filename, const std::string& xml_data)
208{
209 PrintXML(filename, xml_data);
210}
211
212void CommSerial::PrintGrid(Grid& grid, const char* information)
213{
214 Index i;
215 std::stringstream out;
216 std::ofstream out_file;
217
218 OpenFileAndPrintHeader(out_file, grid, information);
219
220 for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z())
221 for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
222 for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
223 out << std::scientific << grid.GetVal(i) << std::endl;
224
225 out_file << out.str();
226 out_file.close();
227}
228
229void CommSerial::PrintDefect(Grid& sol, Grid& rhs, const char* information)
230{
231 Grid::iterator iter;
232 std::ofstream out_file;
233 std::stringstream out;
234
235 TempGrid *temp = MG::GetTempGrid();
236 temp->SetProperties(sol);
237 temp->ImportFromResidual(sol, rhs);
238
239 OpenFileAndPrintHeader(out_file, *temp, information);
240
241 for (iter=temp->Iterators().Local().Begin(); iter!=temp->Iterators().Local().End(); ++iter)
242 out << std::scientific << temp->GetVal(*iter) << std::endl;
243
244 out_file << out.str();
245 out_file.close();
246}
247
248void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& grid, const char* information)
249{
250 char path_str[129];
251 int count = OutputCount();
252
253 sprintf(path_str, "%s%04d.dat", OutputPath().c_str(), count);
254
255 out.open(path_str, std::ios::trunc);
256
257 out << "# vtk DataFile Version 2.0" << std::endl
258 << count << ": " << information << std::endl
259 << "ASCII" << std::endl
260 << "DATASET STRUCTURED_POINTS" << std::endl
261 << grid.Local().Size().X() << " "
262 << grid.Local().Size().Y() << " "
263 << grid.Local().Size().Z() << std::endl
264 << "ORIGIN 0 0 0" << std::endl
265 << "SPACING " << grid.Extent().MeshWidth().X() << " "
266 << grid.Extent().MeshWidth() << " "
267 << grid.Extent().MeshWidth() << std::endl
268 << "POINT_DATA " << grid.Local().Size().Product() << std::endl
269 << "SCALARS residual double 1" << std::endl
270 << "LOOKUP_TABLE default" << std::endl;
271}
272
273void CommSerial::PrintAllSettings()
274{
275 Multigrid* mg = MG::GetFactory().Get("SOL")->Cast<Multigrid>();
276 std::stringstream buf;
277 std::ofstream out;
278
279 buf << OutputPath() << "settings.txt";
280
281 out.open(buf.str().c_str(), std::ios::trunc);
282
283 for (int i=mg->MinLevel(); i<=mg->MaxLevel(); ++i) {
284
285 out << "###########################################################" << std::endl
286 << "LEVEL: " << i << std::endl
287 << "GLOBAL:" << std::endl
288 << " LOCAL_BEGIN: " << (*mg)(i).Global().LocalBegin() << std::endl
289 << " LOCAL_END: " << (*mg)(i).Global().LocalEnd() << std::endl
290 << " LOCAL_SIZE: " << (*mg)(i).Global().LocalSize() << std::endl
291 << " GLOBAL_FINER_BEGIN: " << (*mg)(i).Global().GlobalFinerBegin() << std::endl
292 << " GLOBAL_FINER_END: " << (*mg)(i).Global().GlobalFinerEnd() << std::endl
293 << " GLOBAL_FINER_SIZE: " << (*mg)(i).Global().GlobalFinerSize() << std::endl
294 << " LOCAL_FINER_BEGIN: " << (*mg)(i).Global().LocalFinerBegin() << std::endl
295 << " LOCAL_FINER_END: " << (*mg)(i).Global().LocalFinerEnd() << std::endl
296 << " LOCAL_FINER_SIZE: " << (*mg)(i).Global().LocalFinerSize() << std::endl
297 << " FINEST_ABS_BEGIN: " << (*mg)(i).Global().FinestAbsBegin() << std::endl
298 << " FINEST_ABS_END: " << (*mg)(i).Global().FinestAbsEnd() << std::endl
299 << " FINEST_ABS_SIZE: " << (*mg)(i).Global().FinestAbsSize() << std::endl
300 << " GLOBAL_SIZE: " << (*mg)(i).Global().GlobalSize() << std::endl
301 << " BOUNDARY_TYPE: " << (*mg)(i).Global().BoundaryType() << std::endl
302 << "LOCAL:" << std::endl
303 << " BEGIN: " << (*mg)(i).Local().Begin() << std::endl
304 << " END: " << (*mg)(i).Local().End() << std::endl
305 << " SIZE: " << (*mg)(i).Local().Size() << std::endl
306 << " SIZE_TOTAL: " << (*mg)(i).Local().SizeTotal() << std::endl
307 << " HALO_BEGIN_1: " << (*mg)(i).Local().HaloBegin1() << std::endl
308 << " HALO_END_1: " << (*mg)(i).Local().HaloEnd1() << std::endl
309 << " HALO_SIZE_1: " << (*mg)(i).Local().HaloSize1() << std::endl
310 << " HALO_BEGIN_2: " << (*mg)(i).Local().HaloBegin2() << std::endl
311 << " HALO_END_2: " << (*mg)(i).Local().HaloEnd2() << std::endl
312 << " HALO_SIZE_2: " << (*mg)(i).Local().HaloSize2() << std::endl
313 << " BOUNDARY_BEGIN_1: " << (*mg)(i).Local().BoundaryBegin1() << std::endl
314 << " BOUNDARY_END_1: " << (*mg)(i).Local().BoundaryEnd1() << std::endl
315 << " BOUNDARY_SIZE_1: " << (*mg)(i).Local().BoundarySize1() << std::endl
316 << " BOUNDARY_BEGIN_2: " << (*mg)(i).Local().BoundaryBegin2() << std::endl
317 << " BOUNDARY_END_2: " << (*mg)(i).Local().BoundaryEnd2() << std::endl
318 << " BOUNDARY_SIZE_2: " << (*mg)(i).Local().BoundarySize2() << std::endl
319 << " FINER_BEGIN: " << (*mg)(i).Local().FinerBegin() << std::endl
320 << " FINER_END: " << (*mg)(i).Local().FinerEnd() << std::endl
321 << " FINER_SIZE: " << (*mg)(i).Local().FinerSize() << std::endl
322 << "EXTENT:" << std::endl
323 << " BEGIN: " << (*mg)(i).Extent().Begin() << std::endl
324 << " END: " << (*mg)(i).Extent().End() << std::endl
325 << " SIZE: " << (*mg)(i).Extent().Size() << std::endl
326 << " MESH_WIDTH: " << (*mg)(i).Extent().MeshWidth() << std::endl
327 << "###########################################################" << std::endl;
328
329 }
330
331 assert(out.good());
332 out.close();
333}
334
335inline int GetIndex(const Grid& grid, int i, int j, int k)
336{
337 if (grid.Global().BoundaryType() == LocallyRefined)
338 return k + grid.Local().Size().Z() * (j + grid.Local().Size().Y() * i);
339 else
340 return k + grid.Local().SizeTotal().Z() * (j + grid.Local().SizeTotal().Y() * i);
341}
342
343void CommSerial::PrintGridStructureLevel(Grid& grid, std::ofstream& out)
344{
345 const Vector& sp = grid.Extent().MeshWidth();
346 int numLines;
347
348 if (grid.Global().BoundaryType() == LocallyRefined)
349 numLines = grid.Local().Size().X() * grid.Local().Size().Y() +
350 grid.Local().Size().Y() * grid.Local().Size().Z() +
351 grid.Local().Size().X() * grid.Local().Size().Z();
352 else
353 numLines = grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Y() +
354 grid.Local().SizeTotal().Y() * grid.Local().SizeTotal().Z() +
355 grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Z();
356
357 out << " <Piece";
358
359 if (grid.Global().BoundaryType() == LocallyRefined) {
360 out << " NumberOfPoints=\"" << grid.Local().Size().Product() << "\""
361 << " NumberOfVerts=\"" << grid.Local().Size().Product() << "\"";
362 }else {
363 out << " NumberOfPoints=\"" << grid.Local().SizeTotal().Product() << "\""
364 << " NumberOfVerts=\"" << grid.Local().SizeTotal().Product() << "\"";
365 }
366
367 out << " NumberOfLines=\"" << numLines << "\""
368 << " NumberOfStrips=\"0\""
369 << " NumberOfPolys=\"0\"" << ">" << std::endl
370 << " <Points>" << std::endl
371 << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl
372 << " ";
373
374 if (grid.Global().BoundaryType() == LocallyRefined) {
375 for (int i=0; i<grid.Local().Size().X(); i++)
376 for (int j=0; j<grid.Local().Size().Y(); j++)
377 for (int k=0; k<grid.Local().Size().Z(); k++)
378 out << grid.Extent().Begin().X() + i * sp.X() << " "
379 << grid.Extent().Begin().Y() + j * sp.Y() << " "
380 << grid.Extent().Begin().Z() + k * sp.Z() << " ";
381 }else {
382 for (int i=0; i<grid.Local().SizeTotal().X(); i++)
383 for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
384 for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
385 out << grid.Extent().Begin().X() + i * sp.X() << " "
386 << grid.Extent().Begin().Y() + j * sp.Y() << " "
387 << grid.Extent().Begin().Z() + k * sp.Z() << " ";
388 }
389
390 out << std::endl
391 << " </DataArray>" << std::endl
392 << " </Points>" << std::endl
393 << " <Verts>" << std::endl
394 << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
395 << " ";
396
397 if (grid.Global().BoundaryType() == LocallyRefined) {
398 for (int i=0; i<grid.Local().Size().Product(); i++)
399 out << i << " ";
400 }else {
401 for (int i=0; i<grid.Local().SizeTotal().Product(); i++)
402 out << i << " ";
403 }
404
405 out << std::endl
406 << " </DataArray>" << std::endl
407 << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
408 << " ";
409
410 if (grid.Global().BoundaryType() == LocallyRefined) {
411 for (int i=1; i<=grid.Local().Size().Product(); i++)
412 out << i << " ";
413 }else {
414 for (int i=1; i<=grid.Local().SizeTotal().Product(); i++)
415 out << i << " ";
416 }
417
418 out << std::endl
419 << " </DataArray>" << std::endl
420 << " </Verts>" << std::endl
421 << " <Lines>" << std::endl
422 << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
423 << " ";
424
425 if (grid.Global().BoundaryType() == LocallyRefined) {
426 for (int i=0; i<grid.Local().Size().X(); i++)
427 for (int j=0; j<grid.Local().Size().Y(); j++)
428 out << GetIndex(grid,i,j,0) << " "
429 << GetIndex(grid,i,j,grid.Local().Size().Z()-1) << " ";
430
431 for (int j=0; j<grid.Local().Size().Y(); j++)
432 for (int k=0; k<grid.Local().Size().Z(); k++)
433 out << GetIndex(grid,0,j,k) << " "
434 << GetIndex(grid,grid.Local().Size().X()-1,j,k) << " ";
435
436 for (int i=0; i<grid.Local().Size().X(); i++)
437 for (int k=0; k<grid.Local().Size().Z(); k++)
438 out << GetIndex(grid,i,0,k) << " "
439 << GetIndex(grid,i,grid.Local().Size().Y()-1,k) << " ";
440 }else {
441 for (int i=0; i<grid.Local().SizeTotal().X(); i++)
442 for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
443 out << GetIndex(grid,i,j,0) << " "
444 << GetIndex(grid,i,j,grid.Local().SizeTotal().Z()-1) << " ";
445
446 for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
447 for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
448 out << GetIndex(grid,0,j,k) << " "
449 << GetIndex(grid,grid.Local().SizeTotal().X()-1,j,k) << " ";
450
451 for (int i=0; i<grid.Local().SizeTotal().X(); i++)
452 for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
453 out << GetIndex(grid,i,0,k) << " "
454 << GetIndex(grid,i,grid.Local().SizeTotal().Y()-1,k) << " ";
455 }
456
457 out << std::endl
458 << " </DataArray>" << std::endl
459 << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
460 << " ";
461
462 for (int i=1; i<=numLines; i++)
463 out << 2*i << " ";
464
465 out << std::endl
466 << " </DataArray>" << std::endl
467 << " </Lines>" << std::endl
468 << " </Piece>" << std::endl;
469}
470
471void CommSerial::DebugPrintGridStructure(Multigrid& grid)
472{
473 std::ofstream out;
474 char path_str[129];
475
476 sprintf(path_str, "%sgrid.vtp", OutputPath().c_str());
477
478 out.open(path_str, std::ios::trunc);
479
480 if (!out.good()) {
481#ifdef DEBUG_OUTPUT
482 printf("Multigrid: File %s not accessible.\n", path_str);
483#endif /* DEBUG_OUTPUT */
484 return;
485 }
486
487 out << "<?xml version=\"1.0\"?>" << std::endl
488 << "<VTKFile type=\"PolyData\">" << std::endl
489 << " <PolyData>" << std::endl;
490
491 for (int i=grid.MinLevel(); i<=grid.MaxLevel(); i++)
492 PrintGridStructureLevel(grid(i), out);
493
494 out << " </PolyData>" << std::endl
495 << "</VTKFile>" << std::endl;
496
497 out.close();
498}
499
500std::string CommSerial::CreateOutputDirectory()
501{
502#ifdef HAVE_BOOST_FILESYSTEM
503 std::string path, unique_path;
504 std::stringstream unique_suffix;
505 int suffix_counter = 0;
506 char buffer[129];
507 time_t rawtime;
508 struct tm *timeinfo;
509
510 time(&rawtime);
511 timeinfo = localtime(&rawtime);
512 strftime(buffer, 128, "./output/%Y_%m_%d_%H_%M_%S/", timeinfo);
513 path = buffer;
514
515 if (!fs::exists("output"))
516 fs::create_directory("output");
517
518 unique_path = path;
519
520 while (fs::exists(unique_path.c_str())) {
521
522 unique_suffix.str("");
523 unique_suffix << "_" << suffix_counter++ << "/";
524
525 unique_path = path;
526 unique_path.replace(unique_path.size()-1, 1, unique_suffix.str());
527
528 }
529
530 fs::create_directory(unique_path.c_str());
531
532 return unique_path;
533
534#else
535
536 return "./";
537
538#endif
539}
Note: See TracBrowser for help on using the repository browser.