diff options
author | Erik Schnetter <schnetter@cct.lsu.edu> | 2008-02-29 22:45:27 -0600 |
---|---|---|
committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2008-02-29 22:45:27 -0600 |
commit | 765661c0a8a480dbb2a05c0a4ea97b609d8ec7a1 (patch) | |
tree | f288c5c29bf698f1a8e5a7d369d93804be3c651a /Carpet | |
parent | 18f653b07145c5fccdf22a5cd7dc876411916144 (diff) |
Add new parameter Carpet::grid_coordinates_filename
Add a new parameter Carpet::grid_coordinates_filename which specifies
a file name to which the grid structure is output after regridding.
This is similar to the existing parameter Carpet::grid_structure_filename,
except that the output is in terms of coordinates, which is easier to
interpret.
Diffstat (limited to 'Carpet')
-rw-r--r-- | Carpet/Carpet/param.ccl | 5 | ||||
-rw-r--r-- | Carpet/Carpet/src/Initialise.cc | 3 | ||||
-rw-r--r-- | Carpet/Carpet/src/Recompose.cc | 197 | ||||
-rw-r--r-- | Carpet/Carpet/src/functions.hh | 5 |
4 files changed, 199 insertions, 11 deletions
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl index db2ad25e0..ab643e276 100644 --- a/Carpet/Carpet/param.ccl +++ b/Carpet/Carpet/param.ccl @@ -449,6 +449,11 @@ STRING grid_structure_filename "File name to output grid structure to (empty = n ".*" :: "must be a legal file name" } "" +STRING grid_coordinates_filename "File name to output grid coordinates to (empty = no output)" STEERABLE=recover +{ + ".*" :: "must be a legal file name" +} "" + private: diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc index 0d56f6a90..0afecaa5b 100644 --- a/Carpet/Carpet/src/Initialise.cc +++ b/Carpet/Carpet/src/Initialise.cc @@ -80,10 +80,13 @@ namespace Carpet { CCTKi_InitGHExtensions (cctkGH); +#if 0 // Write grid structure to file for (int m=0; m<maps; ++m) { OutputGridStructure (cctkGH, m, vhh.at(m)->regions); + OutputGridCoordinates (cctkGH, m, vhh.at(m)->regions); } // for m +#endif CallSetup (cctkGH); diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index d25d12e47..3027290de 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -219,6 +219,7 @@ namespace Carpet { // Write grid structure to file OutputGridStructure (cctkGH, m, regsss); + OutputGridCoordinates (cctkGH, m, regsss); if (verbose) OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m)); @@ -477,7 +478,7 @@ namespace Carpet { if (CCTK_MyProc(cctkGH) != 0) return; // Output only if output is desired - if (strcmp(grid_structure_filename, "") == 0) return; + if (CCTK_EQUALS (grid_structure_filename, "")) return; // Create the output directory CCTK_CreateDirectory (0755, out_dir); @@ -499,8 +500,8 @@ namespace Carpet { or stat(filename, &fileinfo)!=0) { file.open (filename, ios::out | ios::trunc); assert (file.good()); - file << "# grid structure" << endl - << "# format: map reflevel component mglevel processor bounding-box is-outer-boundary" << endl; + file << "# grid structure" << eol + << "# format: map reflevel component mglevel processor bounding-box is-outer-boundary" << eol; assert (file.good()); } } @@ -509,22 +510,196 @@ namespace Carpet { assert (file.good()); } - file << "iteration " << cctkGH->cctk_iteration << endl; - file << "maps " << maps << endl; - file << m << " mglevels " << regsss.size() << endl; + file << "iteration " << cctkGH->cctk_iteration << eol; + file << "maps " << maps << eol; + file << m << " mglevels " << regsss.size() << eol; for (int ml=0; ml<(int)regsss.size(); ++ml) { - file << m << " " << ml << " reflevels " << regsss.at(ml).size() << endl; + file << m << " " << ml << " reflevels " << regsss.at(ml).size() << eol; for (int rl=0; rl<(int)regsss.at(ml).size(); ++rl) { - file << m << " " << ml << " " << rl << " components " << regsss.at(ml).at(rl).size() << endl; + file << m << " " << ml << " " << rl << " components " << regsss.at(ml).at(rl).size() << eol; for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) { file << m << " " << ml << " " << rl << " " << c << " " << regsss.at(ml).at(rl).at(c).processor << " " - << regsss.at(ml).at(rl).at(c).extent - << regsss.at(ml).at(rl).at(c).outer_boundaries << endl; + << regsss.at(ml).at(rl).at(c).extent << " " + << regsss.at(ml).at(rl).at(c).outer_boundaries << eol; } } } - file << endl; + file << eol; + + file.close(); + assert (file.good()); + } + + + + void + OutputGridCoordinates (cGH const * const cctkGH, + int const m, + gh::mregs const & regsss) + { + DECLARE_CCTK_PARAMETERS; + + // Output only on the root processor + if (CCTK_MyProc(cctkGH) != 0) return; + + // Output only if output is desired + if (CCTK_EQUALS (grid_coordinates_filename, "")) return; + + // Create the output directory + CCTK_CreateDirectory (0755, out_dir); + + ostringstream filenamebuf; + filenamebuf << out_dir << "/" << grid_coordinates_filename; + // we need a persistent temporary here + string filenamestr = filenamebuf.str(); + const char * filename = filenamestr.c_str(); + + ofstream file; + + static bool do_truncate = true; + + if (do_truncate) { + do_truncate = false; + struct stat fileinfo; + if (IO_TruncateOutputFiles (cctkGH) + or stat(filename, &fileinfo)!=0) { + file.open (filename, ios::out | ios::trunc); + assert (file.good()); + file << "# grid coordinates" << eol + << "# format: map reflevel region mglevel bounding-box" << eol; + assert (file.good()); + } + } + if (not file.is_open()) { + file.open (filename, ios::app); + assert (file.good()); + } + + // + // Find extent of domain + // + + // This requires that CoordBase is used +#warning "TODO: check this (for Carpet, and maybe also for CartGrid3D)" +#warning "TODO: (the check that these two are consistent should be in Carpet)" + + typedef vect<vect<CCTK_INT,2>,3> jjvect; + jjvect nboundaryzones; + jjvect is_internal; + jjvect is_staggered; + jjvect shiftout; + { + CCTK_INT const ierr = GetBoundarySpecification + (2*dim, + & nboundaryzones[0][0], + & is_internal[0][0], + & is_staggered[0][0], + & shiftout[0][0]); + assert (not ierr); + } + + rvect physical_lower, physical_upper; + rvect interior_lower, interior_upper; + rvect exterior_lower, exterior_upper; + rvect spacing; + { + CCTK_INT const ierr = GetDomainSpecification + (dim, + & physical_lower[0], & physical_upper[0], + & interior_lower[0], & interior_upper[0], + & exterior_lower[0], & exterior_upper[0], + & spacing[0]); + assert (not ierr); + } + + // Adapt spacing for convergence level +#warning "TODO: take ml into account" + spacing *= ipow ((CCTK_REAL) mgfact, basemglevel); + + { + CCTK_INT const ierr = + ConvertFromPhysicalBoundary + (dim, + & physical_lower[0], & physical_upper[0], + & interior_lower[0], & interior_upper[0], + & exterior_lower[0], & exterior_upper[0], + & spacing[0]); + assert (not ierr); + } + + // Affine transformation between index space and coordinate space + rvect const origin = exterior_lower; + rvect const scale = + rvect (vhh.at(m)->baseextents.at(0).at(0).stride()) / spacing; + + + + file << "iteration " << cctkGH->cctk_iteration << eol; + file << "maps " << maps << eol; + file << m << " mglevels " << regsss.size() << eol; + for (int ml=0; ml<(int)regsss.size(); ++ml) { + file << m << " " << ml << " reflevels " << regsss.at(ml).size() << eol; + for (int rl=0; rl<(int)regsss.at(ml).size(); ++rl) { + ibset extents; + for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) { + extents += regsss.at(ml).at(rl).at(c).extent; + } + extents.normalize(); + file << m << " " << ml << " " << rl << " regions " << extents.setsize() << eol; + int c=0; + for (ibset::const_iterator + bi = extents.begin(); bi != extents.end(); ++ bi) + { +#if 0 + ibbox const & ext = * bi; + ibbox const & baseext = vhh.at(m)->baseextents.at(ml).at(rl); + ibbox const & coarseext = vhh.at(m)->baseextents.at(ml).at(0 ); + + // This is nice, but is wrong since CartGrid3D has not yet + // initialised the coordinates + ivect const cctk_levfac = spacereffacts.at(rl); + ivect const cctk_levoff = baseext.lower() - coarseext.lower(); + ivect const cctk_levoffdenom = baseext.stride(); + + ivect const cctk_lbnd = + (ext.lower() - baseext.lower()) / ext.stride(); + ivect const cctk_ubnd = + (ext.upper() - baseext.lower()) / ext.stride(); + + rvect const cctk_origin_space = + origin_space.at(m).at(ml); + rvect const cctk_delta_space = + delta_space.at(m) * rvect (ipow (mgfact, ml)); + + rvect const CCTK_ORIGIN_SPACE = + cctk_origin_space + + cctk_delta_space / rvect (cctk_levfac) * + rvect (cctk_levoff) / rvect (cctk_levoffdenom); + rvect const CCTK_DELTA_SPACE = + cctk_delta_space / rvect (cctk_levfac); + + rvect const rlower = + CCTK_ORIGIN_SPACE + rvect (cctk_lbnd) * CCTK_DELTA_SPACE; + rvect const rupper = + CCTK_ORIGIN_SPACE + rvect (cctk_ubnd) * CCTK_DELTA_SPACE; + rvect const rdelta = + CCTK_DELTA_SPACE; +#endif + + ibbox const & ext = * bi; + + rvect const rlower = origin + rvect (ext.lower()) / scale; + rvect const rupper = origin + rvect (ext.upper()) / scale; + rvect const rdelta = rvect (ext.stride()) / scale; + + rbbox const rb (rlower, rupper, rdelta); + file << m << " " << ml << " " << rl << " " << c << " " << rb << eol; + ++c; + } + } + } + file << eol; file.close(); assert (file.good()); diff --git a/Carpet/Carpet/src/functions.hh b/Carpet/Carpet/src/functions.hh index 44aee74c3..bbe2ff724 100644 --- a/Carpet/Carpet/src/functions.hh +++ b/Carpet/Carpet/src/functions.hh @@ -125,6 +125,11 @@ namespace Carpet { gh::mregs const & regsss); void + OutputGridCoordinates (cGH const * cctkGH, + int const m, + gh::mregs const & regsss); + + void OutputGridStatistics (cGH const * cctkGH); |