source: ThirdParty/vmg/src/comm/comm_serial.cpp@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 7faa5c, checked in by Frederik Heber <heber@…>, 9 years ago

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