diff options
Diffstat (limited to 'CarpetDev/CarpetIOF5/src/output.cc')
| -rw-r--r-- | CarpetDev/CarpetIOF5/src/output.cc | 573 |
1 files changed, 573 insertions, 0 deletions
diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc new file mode 100644 index 000000000..0966e5dfb --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/output.cc @@ -0,0 +1,573 @@ +#include <cctk.h> +#include <cctk_Arguments.h> +#include <cctk_Parameters.h> +#include <cctk_Version.h> + +#include <cassert> +#include <cstdlib> +#include <cstring> +#include <iomanip> +#include <iostream> +#include <limits> +#include <sstream> +#include <string> +#include <vector> + +#include <hdf5.h> +#include <F5/F5F.h> +#include <F5/F5R.h> +#include <F5/F5X.h> +#include <F5/F5iterate.h> +#include <F5/F5uniform.h> + +#include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h" + +#include <bbox.hh> +#include <vect.hh> + +#include <carpet.hh> + +#include "iof5.hh" + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + class output_iterator_t { + cGH *const cctkGH; + + bool const is_multipatch; + + string gridname; + string chartname; + string fragmentname; + string topologyname; + + public: + output_iterator_t (cGH *const cctkGH_) + : cctkGH(cctkGH_), + is_multipatch + (CCTK_IsFunctionAliased("MultiPatch_GetSystemSpecification")) + { + } + + void iterate (hid_t const file) + { + int const rl = 0; + ENTER_LEVEL_MODE(cctkGH, rl) { + output_simulation (file); + } LEAVE_LEVEL_MODE; + } + + + + private: + + void output_simulation (hid_t const file) + { + assert (is_level_mode() and reflevel==0); + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + gridname = generate_gridname(cctkGH); + chartname = generate_chartname(cctkGH); + + cout << indent << "simulation=" << gridname << "\n"; + + ivect const reffact = spacereffacts.AT(reflevel); + F5Path *const globalpath = F5Rcreate_vertex_refinement3D + (file, cctk_time, gridname.c_str(), &v2h(reffact)[0], + chartname.c_str()); + + // Choose (arbitrarily) the root level as default topology, for + // readers which don't understand AMR + F5Rlink_default_vertex_topology (globalpath, &v2h(reffact)[0]); + + // Define iteration + F5Rset_timestep (globalpath, cctk_iteration); + + // Attach Cactus/Carpet metadata + write_metadata (cctkGH, globalpath->Grid_hid); + + // Close topology + F5close (globalpath); + + // Iterate over all maps (on the root level) + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + output_map (file); + } END_MAP_LOOP; + } + + + + void output_map (hid_t const file) + { + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + herr_t herr; + + assert (is_singlemap_mode() and reflevel==0); + + cout << indent << "map=" << Carpet::map << "\n"; + + fragmentname = generate_fragmentname(cctkGH, Carpet::map); + + // Iterate over all refinement levels of this map + int const m = Carpet::map; + BEGIN_GLOBAL_MODE(cctkGH) { + BEGIN_REFLEVEL_LOOP(cctkGH) { + if (vhh.AT(m)->local_components(reflevel) > 0) { + // Continue only if this process owns a component of this + // map (on this level) + ENTER_SINGLEMAP_MODE(cctkGH, m, CCTK_GF) { + output_reflevel (file); + } LEAVE_SINGLEMAP_MODE; + } + } END_REFLEVEL_LOOP; + } END_GLOBAL_MODE; + } + + + + void output_reflevel (hid_t const file) + { + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + cout << indent << "reflevel=" << reflevel << "\n"; + + ivect const reffact = spacereffacts.AT(reflevel); + topologyname = generate_topologyname(cctkGH, Carpet::map, reffact); + + // Define grid hierarchy + int const indexdepth = 0; // vertices + F5Path *const path = + F5Rcreate_coordinate_topology (file, &cctk_time, + gridname.c_str(), chartname.c_str(), + topologyname.c_str(), + indexdepth, + dim, dim, &v2h(reffact)[0]); + + // Determine level coordinates + ivect gsh; + rvect origin, delta; + rvect lower, upper; + for (int d=0; d<dim; ++d) { + gsh[d] = cctk_gsh[d]; + origin[d] = CCTK_ORIGIN_SPACE(d); + delta[d] = CCTK_DELTA_SPACE(d); + lower[d] = origin[d]; + upper[d] = lower[d] + (gsh[d]-1) * delta[d]; + } + + if (not is_multipatch) { + // Define level geometry + F5Fwrite_linear (path, FIBER_HDF5_POSITIONS_STRING, + dim, &v2h(gsh)[0], + F5T_COORD3_DOUBLE, + &v2d(lower), &v2d(delta)); + F5Fset_range (path, &v2d(lower), &v2d(upper)); + } + + BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { + output_component (path); + } END_LOCAL_COMPONENT_LOOP; + + // Close topology + F5close (path); + } + + + + void output_component (F5Path *const path) + { + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + cout << indent << "component=" << component << "\n"; + + // Determine component coordinates + ivect gsh; + ivect lbnd, lsh, lghosts, ughosts; + rvect origin, delta; + rvect clower, cupper; + for (int d=0; d<dim; ++d) { + gsh[d] = cctk_gsh[d]; + lbnd[d] = cctk_lbnd[d]; + lsh[d] = cctk_lsh[d]; + // F5 counts the total overlap, which is the sum of the ghost + // zones on this and the adjacent component + lghosts[d] = cctk_bbox[2*d ] ? 0 : 2*cctk_nghostzones[d]; + ughosts[d] = cctk_bbox[2*d+1] ? 0 : 2*cctk_nghostzones[d]; + origin[d] = CCTK_ORIGIN_SPACE(d); + delta[d] = CCTK_DELTA_SPACE(d); + clower[d] = origin[d] + cctk_lbnd[d] * delta[d]; + cupper[d] = clower[d] + (lsh[d]-1) * delta[d]; + } + + // Define coordinates + if (not is_multipatch) { + // (This is redundant, since the level's overall bounding box + // was already defined above, but it provides the individual + // components' bounding boxes.) + F5Fwrite_linear_fraction (path, FIBER_HDF5_POSITIONS_STRING, + cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0], + F5T_COORD3_DOUBLE, + &clower, &delta, &v2h(lbnd)[0], + fragmentname.c_str()); + } else { + output_variable (path, CCTK_VarIndex("grid::x"), true); + } + + // Output some variables + output_variable (path, CCTK_VarIndex("grid::r")); + + output_variable (path, CCTK_VarIndex("grid::x")); + output_variable (path, CCTK_VarIndex("grid::y")); + output_variable (path, CCTK_VarIndex("grid::z")); + } + + + + void output_variable (F5Path *const path, int const var, + bool const write_positions = false) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + indent_t indent; + + int ierr; + + assert (var >= 0); + { + char *const fullname = CCTK_FullName(var); + cout << indent << "variable=" << fullname << "\n"; + free (fullname); + } + + int const group = CCTK_GroupIndexFromVarI(var); + int const v0 = CCTK_FirstVarIndexI(group); + + cGroup groupdata; + ierr = CCTK_GroupData (group, &groupdata); + assert (not ierr); + + // TODO + assert (groupdata.grouptype == CCTK_GF); + assert (groupdata.vartype == CCTK_VARIABLE_REAL); + assert (groupdata.disttype == CCTK_DISTRIB_DEFAULT); + assert (groupdata.stagtype == 0); + assert (groupdata.dim == cctk_dim); + + int const timelevel = 0; + + // Determine component coordinates + int const dim = cctk_dim; + ivect gsh; + ivect lbnd, lsh, lghosts, ughosts; + ivect imin, imax, ioff, ilen; + ivect const izero(0); + for (int d=0; d<dim; ++d) { + gsh[d] = cctk_gsh[d]; + lbnd[d] = cctk_lbnd[d]; + lsh[d] = cctk_lsh[d]; + // F5 counts the total overlap, which is the sum of the ghost + // zones on this and the adjacent component + lghosts[d] = cctk_bbox[2*d ] ? 0 : 2*cctk_nghostzones[d]; + ughosts[d] = cctk_bbox[2*d+1] ? 0 : 2*cctk_nghostzones[d]; + imin[d] = 0; + imax[d] = lsh[d]; + if (not output_ghost_points) { + imin[d] += lghosts[d] / 2; + imax[d] -= ughosts[d] / 2; + } + ioff[d] = lbnd[d] + imin[d]; + ilen[d] = imax[d] - imin[d]; + } + +#warning "TODO: Do not output symmetry zones (unless requested by the user)" +#warning "TODO: Do not output buffer zones (is that easily possible?)" + + + + enum tensortype_t { + tt_scalar, tt_vector, tt_error + }; + + tensortype_t tensortype = tt_error; + + int const coordinates_group = CCTK_GroupIndex ("grid::coordinates"); + assert (coordinates_group >= 0); + + if (group == coordinates_group) { + + // Special case + if (var >= v0 and var < v0+dim) { + tensortype = tt_vector; + } else if (var == v0+dim) { + tensortype = tt_scalar; + } else { + assert (0); + } + + } else { + + // TODO: use the tensor type instead + switch (groupdata.numvars) { + case 1: + tensortype = tt_scalar; break; + case 3: + tensortype = tt_vector; break; + default: + // fallback: output group as scalars + tensortype = tt_scalar; break; + } + + } + + if (tensortype != tt_scalar and var != v0) { + // TODO: don't do this; instead, keep track of which groups + // have been output + char *const fullname = CCTK_FullName(var); + indent_t indent; + cout << indent << "skipping output since it is not the first variable in its group\n"; + free (fullname); + return; + } + + + + switch (tensortype) { + + case tt_scalar: { + // Scalar field + CCTK_REAL const *const rdata = + (CCTK_REAL const*)CCTK_VarDataPtrI (cctkGH, timelevel, var); + assert (rdata); + + int const vartype = CCTK_VarTypeI(var); + assert (vartype == CCTK_VARIABLE_REAL); + vector<CCTK_REAL> idata (prod(ilen)); + for (int k=0; k<ilen[2]; ++k) { + for (int j=0; j<ilen[1]; ++j) { + for (int i=0; i<ilen[0]; ++i) { + int const isrc = + CCTK_GFINDEX3D(cctkGH, imin[0]+i,imin[1]+j,imin[2]+k); + int const idst = i + ilen[0] * (j + ilen[1] * k); + idata[idst] = rdata[isrc]; + } + } + } + void const *const data = &idata.front(); + + char *const groupname = CCTK_GroupName(group); + char *const fullname = CCTK_FullName(var); + // Use the variable name instead of the group name if we may + // output several variables per group + assert (not write_positions); + char const *const name = groupdata.numvars == 1 ? groupname : fullname; +#if 0 + F5Fwrite_fraction (path, name, + dim, + is_multipatch ? NULL : &v2h(gsh)[0], &v2h(lsh)[0], + H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, + data, + &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0], + fragmentname.c_str(), H5P_DEFAULT); +#endif + F5Fwrite_fraction (path, name, + dim, + is_multipatch ? NULL : &v2h(gsh)[0], &v2h(ilen)[0], + H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, + data, + &v2h(ioff)[0], &v2h(izero)[0], &v2h(izero)[0], + fragmentname.c_str(), H5P_DEFAULT); + free (groupname); + free (fullname); + break; + } + + case tt_vector: { + // Vector field + CCTK_REAL const* rdata[cctk_dim]; + vector<CCTK_REAL> idata[cctk_dim]; + void const* data[cctk_dim]; + for (int d=0; d<cctk_dim; ++d) { + rdata[d] = + (CCTK_REAL const*)CCTK_VarDataPtrI (cctkGH, timelevel, v0+d); + assert (rdata[d]); + + int const vartype = CCTK_VarTypeI(var); + assert (vartype == CCTK_VARIABLE_REAL); + idata[d].resize (prod(ilen)); + for (int k=0; k<ilen[2]; ++k) { + for (int j=0; j<ilen[1]; ++j) { + for (int i=0; i<ilen[0]; ++i) { + int const isrc = + CCTK_GFINDEX3D(cctkGH, imin[0]+i,imin[1]+j,imin[2]+k); + int const idst = i + ilen[0] * (j + ilen[1] * k); + idata[d][idst] = rdata[d][isrc]; + } + } + } + data[d] = &idata[d].front(); + + } + char *const groupname = CCTK_GroupName(group); + char const *const name = + write_positions ? FIBER_HDF5_POSITIONS_STRING : groupname; + hid_t const type = + write_positions ? F5T_COORD3_DOUBLE : F5T_VEC3_DOUBLE; + int const will_cover_complete_domain = 0; +#if 0 + F5FSwrite_fraction (path, name, + dim, + is_multipatch ? NULL : &v2h(gsh)[0], &v2h(lsh)[0], + type, type, + data, + &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0], + fragmentname.c_str(), H5P_DEFAULT, + will_cover_complete_domain); +#endif + F5FSwrite_fraction (path, name, + dim, + is_multipatch ? NULL : &v2h(gsh)[0], &v2h(ilen)[0], + type, type, + data, + &v2h(ioff)[0], &v2h(izero)[0], &v2h(izero)[0], + fragmentname.c_str(), H5P_DEFAULT, + will_cover_complete_domain); + free (groupname); + break; + } + + default: + assert (0); + } + } + + }; // class output_iterator_t + + + + void write_metadata (cGH *const cctkGH, hid_t const file) + { + DECLARE_CCTK_PARAMETERS; + + assert (cctkGH); + assert (file >= 0); + + herr_t herr; + + + CCTK_INFO ("Writing simulation metadata..."); + + // NOTE: We could write the metadata into only one of the process + // files (to save space), or write it into a separate metadata + // file + + // Create a group to hold all metadata + hid_t const group = + H5Gcreate (file, metadata_group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert (group >= 0); + + // General information + WriteAttribute (group, "Cactus version", CCTK_FullVersion()); + + // Unique identifiers + if (CCTK_IsFunctionAliased ("UniqueConfigID")) { + WriteAttribute + (group, "config id", (char const*) UniqueConfigID(cctkGH)); + } + if (CCTK_IsFunctionAliased ("UniqueBuildID")) { + WriteAttribute + (group, "build id", (char const*) UniqueBuildID(cctkGH)); + } + if (CCTK_IsFunctionAliased ("UniqueSimulationID")) { + WriteAttribute + (group, "simulation id", (char const*) UniqueSimulationID(cctkGH)); + } + if (CCTK_IsFunctionAliased ("UniqueRunID")) { + WriteAttribute + (group, "run id", (char const*) UniqueRunID(cctkGH)); + } + + // Number of I/O processes (i.e. the number of output files) + int const nprocs = CCTK_nProcs(cctkGH); + int const nioprocs = nprocs; + WriteAttribute (group, "nioprocs", nioprocs); + + // Write parameters into a separate dataset (they may be too large + // for an attribute) + int const get_all = 1; + char *const parameters = IOUtil_GetAllParameters (cctkGH, get_all); + assert (parameters); + WriteLargeAttribute (group, all_parameters, parameters); + free (parameters); + + // Grid structure + string const gs = serialise_grid_structure (cctkGH); + WriteLargeAttribute (group, grid_structure, gs.c_str()); + + herr = H5Gclose (group); + assert (not herr); + } + + + + extern "C" + void F5_Output (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Output: iteration=%d", cctk_iteration); + + + + // We don't know how to open multiple files yet + assert (strcmp (file_content, "everything") == 0); + + // Open file + static bool first_time = true; + + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + int const myproc = CCTK_MyProc(cctkGH); + int const proc = myproc; + string const name = + create_filename (cctkGH, basename, proc, first_time); + + indent_t indent; + cout << indent << "process=" << proc << "\n"; + + bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH); + hid_t const file = + truncate_file ? + H5Fcreate (name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : + H5Fopen (name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + assert (file >= 0); + first_time = false; + + // write_metadata (cctkGH, file); + + output_iterator_t iterator(cctkGH); + iterator.iterate (file); + + // Close file + herr = H5Fclose (file); + assert (not herr); + } + +} // end namespace CarpetIOF5 |