source: src/comm/comm_serial.cpp@ f57182

Last change on this file since f57182 was f57182, checked in by Julian Iseringhausen <isering@…>, 13 years ago

Open boundary conditions.

Conflicts:

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