From ec25e62821c80cfcdfeb5c6331899df8e655b34c Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 18 Nov 2012 17:43:00 -0500 Subject: CarpetIOF5: Various improvements --- CarpetDev/CarpetIOF5/interface.ccl | 1 + CarpetDev/CarpetIOF5/par/iof5-array.par | 8 +- .../par/iof5-multipatch-kerrschild-input.par | 2 + .../CarpetIOF5/par/iof5-multipatch-kerrschild.par | 2 + CarpetDev/CarpetIOF5/param.ccl | 11 +- CarpetDev/CarpetIOF5/schedule.ccl | 28 --- CarpetDev/CarpetIOF5/src/distribute.cc | 125 ++++++++++-- CarpetDev/CarpetIOF5/src/distribute.hh | 26 ++- CarpetDev/CarpetIOF5/src/input.cc | 54 ++++- CarpetDev/CarpetIOF5/src/iof5.cc | 90 ++++----- CarpetDev/CarpetIOF5/src/iof5.hh | 8 +- CarpetDev/CarpetIOF5/src/output.cc | 34 +++- CarpetDev/CarpetIOF5/src/util.cc | 221 ++++++++++++++++++++- 13 files changed, 489 insertions(+), 121 deletions(-) (limited to 'CarpetDev') diff --git a/CarpetDev/CarpetIOF5/interface.ccl b/CarpetDev/CarpetIOF5/interface.ccl index de552a6b2..59cea71d5 100644 --- a/CarpetDev/CarpetIOF5/interface.ccl +++ b/CarpetDev/CarpetIOF5/interface.ccl @@ -4,6 +4,7 @@ IMPLEMENTS: CarpetIOF5 INHERITS: grid +USES INCLUDE: cacheinfo.hh USES INCLUDE: carpet.hh USES INCLUDE: CarpetTimers.hh diff --git a/CarpetDev/CarpetIOF5/par/iof5-array.par b/CarpetDev/CarpetIOF5/par/iof5-array.par index 8c9fd2a03..26de7054d 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-array.par +++ b/CarpetDev/CarpetIOF5/par/iof5-array.par @@ -29,8 +29,6 @@ Cactus::cctk_itlast = 2 IO::out_dir = $parfile -Carpet::processor_topology = "balanced" - Carpet::domain_from_coordbase = yes CartGrid3D::type = "coordbase" CoordBase::domainsize = "minmax" @@ -40,9 +38,9 @@ CoordBase::zmin = -1.0 CoordBase::xmax = +1.0 CoordBase::ymax = +1.0 CoordBase::zmax = +1.0 -CoordBase::dx = 0.25 -CoordBase::dy = 0.25 -CoordBase::dz = 0.25 +CoordBase::dx = 0.025 +CoordBase::dy = 0.025 +CoordBase::dz = 0.025 InterpToArray::nscalars = 1 InterpToArray::scalar_vars[0] = "grid::r" diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par index c28ade8a0..d6c5cbe5c 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par @@ -145,6 +145,8 @@ CarpetIOASCII::out1D_vars = " +CarpetIOF5::verbose = yes + CarpetIOF5::out_every = 1024 CarpetIOF5::out_vars = " ADMBase::metric diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par index 171ef919e..8ddd56f94 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par @@ -150,6 +150,8 @@ CarpetIOASCII::out1D_vars = " +CarpetIOF5::verbose = yes + CarpetIOF5::out_every = 1024 CarpetIOF5::out_vars = " ADMBase::metric diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl index b83b74675..bd8f14332 100644 --- a/CarpetDev/CarpetIOF5/param.ccl +++ b/CarpetDev/CarpetIOF5/param.ccl @@ -154,10 +154,19 @@ STRING checkpoint_dir "Checkpoint directory (overrides IO::checkpoint_dir)" STEE +BOOLEAN use_chunks "Always use chunked layout" STEERABLE=always +{ +} "yes" + +INT max_chunksize "Maximum chunk size (in bytes)" STEERABLE=always +{ + 1:* :: "" +} 131072 # 128 kB + INT compression_level "Compression level" STEERABLE=always { -1 :: "no compression" - 0:9 :: "higher numbers compress better" + 0:9 :: "smaller numbers compress faster, larger numbers compress better" } 1 BOOLEAN use_checksums "Include a checksum for the data" STEERABLE=always diff --git a/CarpetDev/CarpetIOF5/schedule.ccl b/CarpetDev/CarpetIOF5/schedule.ccl index da8745b56..e633aef18 100644 --- a/CarpetDev/CarpetIOF5/schedule.ccl +++ b/CarpetDev/CarpetIOF5/schedule.ccl @@ -46,31 +46,3 @@ schedule CarpetIOF5_TerminationCheckpoint at TERMINATE LANG: C OPTIONS: global } "Termination checkpoint routine" - - - -#SCHEDULE F5_Output AT initial -#{ -# LANG: C -# OPTIONS: global-late -#} "Create an output file" -# -# -# -#SCHEDULE F5_Poison AT initial AFTER F5_Output BEFORE F5_Input -#{ -# LANG: C -# OPTIONS: global-late -#} "Poison all variables" -# -#SCHEDULE F5_Input AT initial AFTER F5_Output -#{ -# LANG: C -# OPTIONS: global-late -#} "Read from file" -# -#SCHEDULE F5_Check AT initial AFTER F5_Input -#{ -# LANG: C -# OPTIONS: global-late -#} "Check all variables for poison" diff --git a/CarpetDev/CarpetIOF5/src/distribute.cc b/CarpetDev/CarpetIOF5/src/distribute.cc index 5fd8e5a29..ab7183df3 100644 --- a/CarpetDev/CarpetIOF5/src/distribute.cc +++ b/CarpetDev/CarpetIOF5/src/distribute.cc @@ -2,6 +2,7 @@ #include +#include #include #include @@ -52,10 +53,48 @@ namespace CarpetIOF5 { + /*** scatter_t::busy_tags_t *************************************************/ + + scatter_t::busy_tags_t::busy_tags_t() + { + // do nothing + } + + scatter_t::busy_tags_t::~busy_tags_t() + { + // Ensure all tags are free + for (int tag=0; tag<(int)tags.size(); ++tag) { + assert(not tags.at(tag)); + } + } + + int scatter_t::busy_tags_t::allocate_tag() + { + if (tags.empty()) { + tags.resize(max_concurrent_sends, false); + } + for (int tag=0; tag<(int)tags.size(); ++tag) { + if (not tags.at(tag)) { + tags.at(tag) = true; + return tag; + } + } + return -1; // no free tag available + } + + void scatter_t::busy_tags_t::free_tag(int const tag) + { + assert(tags.at(tag)); + tags.at(tag) = false; + } + + + /*** scatter_t **************************************************************/ scatter_t::scatter_t(cGH const *const cctkGH_) : cctkGH(cctkGH_), + busy_tags(CCTK_nProcs(cctkGH)), num_received(0), num_sent(0), bytes_allocated(0), did_send_all(false), did_receive_sent_all(0), did_receive_all_sent_all(false) @@ -136,7 +175,7 @@ namespace CarpetIOF5 { } if (verbose) { - CCTK_INFO("Destroying down global scatter object"); + CCTK_INFO("Destroying global scatter object"); } assert(public_recvs.empty()); @@ -152,7 +191,7 @@ namespace CarpetIOF5 { // consecutive) int child() { return dist::rank()*nsiblings() + 1; } // Get process id if my parent - int parent() { return (dist::rank()-1) / nsiblings(); } + int parent() { return dist::rank()==0 ? -1 : (dist::rank()-1) / nsiblings(); } // Check whether we should tell our parent that all our (ours and // our childrens') messages have been sent @@ -164,7 +203,8 @@ namespace CarpetIOF5 { if (p < dist::size()) ++need_receive; } if (did_send_all and did_receive_sent_all == need_receive) { - if (dist::rank() == 0) { + int const p = parent(); + if (p < 0) { // We are root; now broadcast this to all set_did_receive_all_sent_all(); } else { @@ -172,7 +212,6 @@ namespace CarpetIOF5 { CCTK_INFO("[Telling our parent that we and our children sent all our data]"); } to_parent.state = fragdesc_t::state_sent_all; - int const p = parent(); MPI_Request request; MPI_Isend(&to_parent, to_parent.num_ints(), MPI_INT, p, tag_desc, dist::comm(), &request); @@ -265,10 +304,25 @@ namespace CarpetIOF5 { while (not tosend.empty()) { list::iterator const tmi = tosend.begin(); transmission_t *const tm = *tmi; + tosend.erase(tmi); + + // Choose tag + int const proc = tm->fragdesc.process; + int tag; + while (true) { + tag = busy_tags.at(proc).allocate_tag(); + if (tag >= 0) break; + // No tag was available; wait for some progress before trying + // again (until the receiver has received some of our + // messages) + do_some_work(true); + } + tm->fragdesc.tag = tag_data_min + tag; if (verbose) { CCTK_VInfo(CCTK_THORNSTRING, - " Sending to process %d...", tm->fragdesc.process); + " Sending data to process %d with tag %d...", + tm->fragdesc.process, tm->fragdesc.tag); } // Send descriptor and data @@ -276,11 +330,10 @@ namespace CarpetIOF5 { MPI_Isend(&tm->fragdesc, tm->fragdesc.num_ints(), MPI_INT, tm->fragdesc.process, tag_desc, dist::comm(), &req); + assert(tm->request == MPI_REQUEST_NULL); MPI_Isend(&tm->data[0], tm->fragdesc.npoints(), tm->fragdesc.datatype(), - tm->fragdesc.process, tag_data, + tm->fragdesc.process, tm->fragdesc.tag, dist::comm(), &tm->request); - - tosend.erase(tmi); sends.push_back(tm); } @@ -352,8 +405,9 @@ namespace CarpetIOF5 { requests.push_back(tm->request); iterators.push_back(tmi); } + assert((int)requests.size() == nrequests); - // Wait for (or test for) some open requests + // Wait (or test) for some open requests int outcount; vector indices(nrequests); vector statuses(nrequests); @@ -377,32 +431,38 @@ namespace CarpetIOF5 { // Process all completed requests for (int n=0; n=0 and source::iterator const tmi = iterators.at(idx); transmission_t *const tm = *tmi; + assert(tm->request != MPI_REQUEST_NULL); + tm->request = MPI_REQUEST_NULL; public_recvs.erase(tmi); if (tm->fragdesc.state != fragdesc_t::state_normal) { handle_state_transition(tm->fragdesc); } else { + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + "Receiving data from process %d with tag %d", + source, tm->fragdesc.tag); + } + // Prepare receiving the data assert(tm->fragdesc.process == dist::rank()); tm->data.resize(tm->fragdesc.npoints() * tm->fragdesc.vartypesize()); bytes_allocated += tm->data.size(); + assert(tm->request == MPI_REQUEST_NULL); MPI_Irecv(&tm->data[0], tm->fragdesc.npoints(), tm->fragdesc.datatype(), - source, tag_data, + source, tm->fragdesc.tag, dist::comm(), &tm->request); recvs.push_back(tm); if (verbose) { @@ -418,6 +478,8 @@ namespace CarpetIOF5 { // We completed receiving a dataset; process it list::iterator const tmi = iterators.at(idx); transmission_t *const tm = *tmi; + assert(tm->request != MPI_REQUEST_NULL); + tm->request = MPI_REQUEST_NULL; if (verbose) { char *const fullname = CCTK_FullName(tm->fragdesc.varindex); @@ -441,6 +503,8 @@ namespace CarpetIOF5 { // We completed sending a dataset; forget it list::iterator const tmi = iterators.at(idx); transmission_t *const tm = *tmi; + assert(tm->request != MPI_REQUEST_NULL); + tm->request = MPI_REQUEST_NULL; if (verbose) { char *const fullname = CCTK_FullName(tm->fragdesc.varindex); @@ -449,6 +513,10 @@ namespace CarpetIOF5 { free(fullname); } + int const proc = tm->fragdesc.process; + int const tag = tm->fragdesc.tag - tag_data_min; + busy_tags.at(proc).free_tag(tag); + bytes_allocated -= tm->data.size(); delete tm; if (verbose) { @@ -484,11 +552,28 @@ namespace CarpetIOF5 { baseext.lower() + fd.imax * baseext.stride(), baseext.stride()); + // Distribute load: don't start with sending to process 0, instead + // start with sending to (self+1) + int coffset; + { + int const nlc = hh.local_components(fd.reflevel); + if (nlc > 0) { + // find component of next processor + int c0 = hh.get_component(fd.reflevel, 0); + coffset = c0 + nlc; + } else { + // estimate component of next processor + int const nc = hh.components(fd.reflevel); + coffset = nc * dist::rank() / dist::size(); + } + } + ibset done; list tosend; dh::light_cboxes const& light_cbox = dd.light_boxes.at(fd.mglevel).at(fd.reflevel); - for (int c=0; c tags; + public: + busy_tags_t(); + ~busy_tags_t(); + int allocate_tag(); + void free_tag(int tag); + }; + vector busy_tags; struct transmission_t { fragdesc_t fragdesc; vector data; MPI_Request request; + transmission_t(): request(MPI_REQUEST_NULL) {} + ~transmission_t() { assert(request==MPI_REQUEST_NULL); } }; - cGH const *const cctkGH; - // Desired number of public receives to keep posted at all times static int const num_public_recvs = 10; list public_recvs, recvs, sends; @@ -68,7 +86,7 @@ namespace CarpetIOF5 { fragdesc_t to_parent, to_children; public: - scatter_t(cGH const *const cctkGH_); + scatter_t(cGH const *cctkGH_); ~scatter_t(); void send(fragdesc_t const& fd, void const *data); diff --git a/CarpetDev/CarpetIOF5/src/input.cc b/CarpetDev/CarpetIOF5/src/input.cc index b7f25b59f..eb2db82e5 100644 --- a/CarpetDev/CarpetIOF5/src/input.cc +++ b/CarpetDev/CarpetIOF5/src/input.cc @@ -119,15 +119,41 @@ namespace CarpetIOF5 { // TODO: synchronise all read grid functions } - + + + void read_grid(F5Path *const path) { indent_t indent; cout << indent << "grid=" << gridname << "\n"; + + + + if (input_metadata) { + if (reflevel == 0) { + indent_t indent2; + cout << indent2 << "reading metadata\n"; + // hid_t const metadata_group = path->Grid_hid; + ostringstream pathname; + pathname << FIBER_CONTENT_GRIDS << "/" << gridname; + hid_t group; + group = H5Gopen(path->ContentsGroup_hid, pathname.str().c_str(), + H5P_DEFAULT); + assert(group >= 0); + read_metadata(cctkGH, group); + herr_t const herr = H5Gclose(group); + assert(not herr); + } + } + + + // F5iterate_vertex_fields(path, NULL, field_iterator, this, NULL, NULL); F5iterate_topologies(path, NULL, topology_iterator, this); } + + void read_topology(F5Path *const path) { indent_t indent; @@ -151,9 +177,9 @@ namespace CarpetIOF5 { herr = H5Gget_linkval (path->Grid_hid, topologyname, sizeof linkval, linkval); assert(not herr); - indent_t indent; - cout << indent << "alias for topology \"" << linkval << "\"\n" - << indent << "ignoring this topology\n"; + indent_t indent2; + cout << indent2 << "alias for topology \"" << linkval << "\"\n" + << indent2 << "ignoring this topology\n"; return; } @@ -185,6 +211,8 @@ namespace CarpetIOF5 { F5iterate_topology_fields(path, NULL, field_iterator, this, NULL, NULL); } + + void read_field(F5Path *const path) { indent_t indent; @@ -220,8 +248,8 @@ namespace CarpetIOF5 { F5Fclose(path); } else { - indent_t indent; - cout << indent << "ignoring this field\n"; + indent_t indent2; + cout << indent2 << "ignoring this field\n"; } // TODO: keep track of which fields have been read, and complain // about unread ones @@ -229,6 +257,8 @@ namespace CarpetIOF5 { // complain about unread parts } + + void read_fragment(F5Path *const path) { indent_t indent; @@ -293,6 +323,8 @@ namespace CarpetIOF5 { } } + + void read_variable(F5Path *const path, char const *const name, int const var) { @@ -340,6 +372,10 @@ namespace CarpetIOF5 { // fragment, and (b) there exists a group with this name. (We // probably could examine F5's internal state as well, but // looking for an HDF5 group is the most robust approach.) + + // Note: F5 has a function F5Fis_group(path) that returns + // whether the field is a group or a dataset. + bool fragment_is_group = false; if (fragmentname) { H5O_info_t info; @@ -432,8 +468,8 @@ namespace CarpetIOF5 { ibbox const fbox(foff, foff+flen-1, 1); { - indent_t indent; - cout << indent + indent_t indent2; + cout << indent2 << "dataset bbox is " << foff << ":" << foff+flen-1 << "\n"; } @@ -603,7 +639,7 @@ namespace CarpetIOF5 { // Grid structure string gs; ReadLargeAttribute(group, grid_structure, gs); - // TODO: set grid structure + deserialise_grid_structure(cctkGH, gs); herr = H5Gclose(group); assert(not herr); diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc b/CarpetDev/CarpetIOF5/src/iof5.cc index 1141a7471..1642bfe0e 100644 --- a/CarpetDev/CarpetIOF5/src/iof5.cc +++ b/CarpetDev/CarpetIOF5/src/iof5.cc @@ -86,7 +86,7 @@ namespace CarpetIOF5 { // A mechanism to keep the HDF5 output file open across multiple // write operations: - hid_t file = H5I_INVALID_HID; + hid_t open_file = H5I_INVALID_HID; int keep_file_open = 0; void enter_keep_file_open() { @@ -96,10 +96,10 @@ namespace CarpetIOF5 { { assert(keep_file_open > 0); --keep_file_open; - if (keep_file_open==0 and file>=0) { - herr_t const herr = H5Fclose(file); + if (keep_file_open==0 and open_file>=0) { + herr_t const herr = H5Fclose(open_file); assert(not herr); - file = H5I_INVALID_HID; + open_file = H5I_INVALID_HID; H5garbage_collect(); } } @@ -294,22 +294,22 @@ namespace CarpetIOF5 { enter_keep_file_open(); bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH) and myproc == myioproc; - if (file < 0) { + if (open_file < 0) { // Reuse file hid if file is already open hid_t fapl = H5Pcreate(H5P_FILE_ACCESS); H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG); - file = + open_file = truncate_file ? H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl) : H5Fopen (name.c_str(), H5F_ACC_RDWR , fapl); - assert(file >= 0); + assert(open_file >= 0); H5Pclose(fapl); } first_time = false; vector output_var(CCTK_NumVars()); output_var.at(vindex) = true; - output(cctkGH, file, output_var, false, true); + output(cctkGH, open_file, output_var, false, true); // Close file leave_keep_file_open(); @@ -471,16 +471,22 @@ namespace CarpetIOF5 { + assert(called_from == CP_RECOVER_PARAMETERS or + called_from == CP_RECOVER_DATA or + called_from == FILEREADER_DATA); + bool const in_recovery = + called_from == CP_RECOVER_PARAMETERS or + called_from == CP_RECOVER_DATA; + if (not in_recovery) { + assert(is_level_mode()); + if (reflevel > 0) return 0; // no error + } + + + BEGIN_GLOBAL_MODE(cctkGH) { DECLARE_CCTK_ARGUMENTS; - assert(called_from == CP_RECOVER_PARAMETERS or - called_from == CP_RECOVER_DATA or - called_from == FILEREADER_DATA); - bool const in_recovery = - called_from == CP_RECOVER_PARAMETERS or - called_from == CP_RECOVER_DATA; - if (in_recovery) { CCTK_VInfo(CCTK_THORNSTRING, "F5::Input: recovering iteration %d", recovery_iteration); @@ -505,35 +511,28 @@ namespace CarpetIOF5 { - scatter_t scatter(cctkGH); + // Keep track of which files could be read, and which could not + int foundproc = -1, notfoundproc = -1; - // Open file // string const basename = // in_recovery // ? "checkpoint" // : generate_basename(cctkGH, CCTK_VarIndex("grid::r")); string const basename = basefilename; - // Keep track of which files could be read, and which could not - int foundproc = -1, notfoundproc = -1; - - // TODO: Store how many processes contributed to the output, and - // expect exactly that many files - int const myproc = CCTK_MyProc(cctkGH); - int const nprocs = CCTK_nProcs(cctkGH); - - int const ioproc_every = - max_nioprocs == 0 ? 1 : (nprocs + max_nioprocs - 1) / max_nioprocs; - assert(ioproc_every > 0); - int const nioprocs = nprocs / ioproc_every; - assert(nioprocs > 0 and nioprocs <= max_nioprocs); - int const myioproc = myproc / ioproc_every * ioproc_every; - - if (myproc == myioproc) { + { + scatter_t scatter(cctkGH); + + // Iterate over files + + // TODO: Store how many processes contributed to the output, + // and expect exactly that many files + int const myproc = CCTK_MyProc(cctkGH); + int const nprocs = CCTK_nProcs(cctkGH); // Loop over all (possible) files - for (int ioproc = myioproc / ioproc_every; ; ioproc += nioprocs) { + for (int proc=myproc; ; proc+=nprocs) { string const name = - create_filename(cctkGH, basename, cctkGH->cctk_iteration, ioproc, + create_filename(cctkGH, basename, cctkGH->cctk_iteration, proc, in_recovery ? io_dir_recover : io_dir_input, false); bool file_exists; @@ -541,24 +540,21 @@ namespace CarpetIOF5 { file_exists = H5Fis_hdf5(name.c_str()) > 0; } H5E_END_TRY; if (not file_exists) { - notfoundproc = ioproc; + notfoundproc = proc; break; } - foundproc = ioproc; + foundproc = proc; indent_t indent; - cout << indent << "I/O process=" << ioproc << "\n"; + cout << indent << "process=" << proc << "\n"; - hid_t fapl = H5Pcreate(H5P_FILE_ACCESS); - H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG); - hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, fapl); + hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); assert(file >= 0); - H5Pclose(fapl); // Iterate over all time slices bool const input_past_timelevels = in_recovery; - // TODO: read metadata when recoverying parameters - bool const input_metadata = false; + // Read metadata when recoverying parameters + bool const input_metadata = called_from == CP_RECOVER_PARAMETERS; input(cctkGH, file, input_var, input_past_timelevels, input_metadata, scatter); @@ -567,6 +563,8 @@ namespace CarpetIOF5 { assert(not herr); H5garbage_collect(); } + + // Destroy scatter object } { @@ -644,10 +642,10 @@ namespace CarpetIOF5 { int const iter = strtol(p, &p, 10); if (!*p) continue; - // Read the process number + // Read (and ignore) the process number if (infix.compare(0, infix.length(), p, infix.length()) != 0) continue; p += infix.length(); - int const proc = strtol(p, &p, 10); + /* int const proc = */ strtol(p, &p, 10); if (!*p) continue; // Finally check the file extension diff --git a/CarpetDev/CarpetIOF5/src/iof5.hh b/CarpetDev/CarpetIOF5/src/iof5.hh index c082206c3..a450ae812 100644 --- a/CarpetDev/CarpetIOF5/src/iof5.hh +++ b/CarpetDev/CarpetIOF5/src/iof5.hh @@ -264,13 +264,19 @@ namespace CarpetIOF5 { + // Write/read Cactus metadata to a particular location in an HDF5 + // file void - write_metadata(cGH const *const cctkGH, hid_t const file); + write_metadata(cGH const *const cctkGH, hid_t const group); + void + read_metadata(cGH const *const cctkGH, hid_t const group); // Handle Carpet's grid structure (this should move to Carpet and/or // CarpetLib) string serialise_grid_structure(cGH const *const cctkGH); + void + deserialise_grid_structure(cGH const *const cctkGH, string const buf); diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc index aad6174f2..d7c2f4649 100644 --- a/CarpetDev/CarpetIOF5/src/output.cc +++ b/CarpetDev/CarpetIOF5/src/output.cc @@ -104,8 +104,7 @@ namespace CarpetIOF5 { struct component_indices_t: map_indices_t { // elements >=dim remain undefined - ivect ash; - ivect lbnd, lsh, lghosts, ughosts; + ivect lbnd, lsh, ash, lghosts, ughosts; ivect imin, imax, ioff, ilen; rvect clower, cupper; @@ -758,14 +757,37 @@ namespace CarpetIOF5 { // Dataset properties hid_t const prop = H5Pcreate(H5P_DATASET_CREATE); assert(prop >= 0); + // Use chunked I/O if requested or needed + bool const set_chunks = + use_chunks or compression_level >= 0 or use_checksums; + if (set_chunks) { + assert(ci.dim > 0); + int const chunksize = max_chunksize / H5Tget_size(type); + assert(chunksize > 0); + int cs = 1, csd = 1; + while (csd < chunksize) { + cs <<= 1; csd <<= ci.dim; + assert(ci.dim > 0); + } + ivect chunkshape = cs; + chunkshape[ci.dim-1] = cs / (csd / chunksize); + chunkshape = max(chunkshape, 1); + chunkshape = min(chunkshape, ci.ilen); + FAILWARN(H5Pset_chunk(prop, ci.dim, + &v2h(chunkshape).reverse()[dim-ci.dim])); + } // Enable compression if requested if (compression_level >= 0) { - FAILWARN(H5Pset_chunk(prop, ci.dim, &v2h(ci.ilen).reverse()[0])); + assert(set_chunks); + indent_t indent2; + cout << indent2 + << "using chunk size " << ci.ilen << " " + << "(" << prod(ci.ilen) << ")\n"; FAILWARN(H5Pset_deflate(prop, compression_level)); } // Enable checksums if requested if (use_checksums) { - FAILWARN(H5Pset_chunk(prop, ci.dim, &v2h(ci.ilen).reverse()[0])); + assert(set_chunks); FAILWARN(H5Pset_fletcher32(prop)); } @@ -773,7 +795,8 @@ namespace CarpetIOF5 { // Write single-component tensors into non-separated fractions // for convenience (could also use a separated compound // instead) - // TODO: Extent F5 API to allow writing non-fragmented datasets + // TODO: Extend F5 API to allow writing non-fragmented datasets + // TODO: Extend F5 API to accept hyperslab descriptors FAILWARN (F5Fwrite_fraction(path, name.c_str(), ci.dim, @@ -798,6 +821,7 @@ namespace CarpetIOF5 { "This does not seem to be supported by F5ls " "(or is implemented wrong in the F5 library)."); } + // TODO: Extend F5 API to accept hyperslab descriptors FAILWARN (F5FSwrite_fraction(path, name.c_str(), ci.dim, diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc index 3d924e775..b2ef6fae2 100644 --- a/CarpetDev/CarpetIOF5/src/util.cc +++ b/CarpetDev/CarpetIOF5/src/util.cc @@ -381,7 +381,7 @@ namespace CarpetIOF5 { - char const *const grid_structure = "Grid Structure v5"; + char const *const grid_structure = "Grid Structure v6"; string serialise_grid_structure(cGH const *const cctkGH) @@ -398,14 +398,233 @@ namespace CarpetIOF5 { buf << "regions:" << vhh.AT(m)->regions.AT(0) << ","; buf << "ghost_widths:" << vdd.AT(m)->ghost_widths << ","; buf << "buffer_widths:" << vdd.AT(m)->buffer_widths << ","; + buf << "overlap_widths:" << vdd.AT(m)->overlap_widths << ","; buf << "prolongation_orders_space:" << vdd.AT(m)->prolongation_orders_space << ","; buf << "},"; } buf << "],"; + buf << "global_time:" << global_time << ","; + buf << "delta_time:" << delta_time << ","; buf << "times:" << *tt << ","; buf << "}."; return buf.str(); } + void + deserialise_grid_structure(cGH const *const cctkGH, string const str) + { + // Define variables for the metadata that will be read in + int const my_maps = maps; + vector my_superregions(my_maps); + vector my_regions(my_maps); + vector > my_ghost_widths(my_maps); + vector > my_buffer_widths(my_maps); + vector > my_overlap_widths(my_maps); + vector my_prolongation_orders_space(my_maps); + CCTK_REAL my_global_time, my_delta_time; + th my_tt(tt->h, tt->time_interpolation_during_regridding); + + // Read in the metadata + istringstream is(str); + skipws(is); + consume(is, grid_structure); + skipws(is); + consume(is, "{"); + skipws(is); + consume(is, "maps"); + skipws(is); + consume(is, '='); + skipws(is); + consume(is, '['); + for (int m=0; m> my_m; + assert(my_m == m); + skipws(is); + consume(is, ']'); + skipws(is); + consume(is, '='); + skipws(is); + consume(is, '{'); + skipws(is); + consume(is, "superregions:"); + skipws(is); + is >> my_superregions.AT(m); + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "regions:"); + skipws(is); + // ignore region specification + // is >> my_regions.AT(m); + gh::rregs tmp_regions; + is >> tmp_regions; + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "ghost_widths:"); + skipws(is); + is >> my_ghost_widths.AT(m); + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "buffer_widths:"); + skipws(is); + is >> my_buffer_widths.AT(m); + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "overlap_widths:"); + skipws(is); + is >> my_overlap_widths.AT(m); + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "prolongation_orders_space:"); + skipws(is); + is >> my_prolongation_orders_space.AT(m); + skipws(is); + consume(is, ','); + skipws(is); + consume(is, '}'); + skipws(is); + consume(is, ','); + } + skipws(is); + consume(is, ']'); + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "global_time:"); + is >> my_global_time; + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "delta_time:"); + is >> my_delta_time; + skipws(is); + consume(is, ','); + skipws(is); + consume(is, "times:"); + is >> my_tt; + skipws(is); + consume(is, ','); + skipws(is); + consume(is, '}'); + skipws(is); + consume(is, '.'); + + // Change Carpet's grid structure according to what has just been + // read + { + int type; + void const *const ptr = + CCTK_ParameterGet("regrid_in_level_mode", "Carpet", & type); + assert(ptr); + assert(type == PARAMETER_BOOLEAN); + CCTK_INT const regrid_in_level_mode = *static_cast(ptr); + assert(regrid_in_level_mode); + } + + // Count refinement levels + vector rls(maps); + for (int m=0; m superregss(maps); + for (int m=0; m regss(maps); + Carpet::SplitRegionsMaps(cctkGH, superregss, regss); + for (int m=0; m my_multigrid_regions(maps); + Carpet::MakeMultigridBoxesMaps(cctkGH, my_regions, my_multigrid_regions); + + // Regrid (ignoring the previous content of the grid hierarchy) + for (int m=0; mmglevels(); ++ ml) { + for (int rl=0; rlreflevels(); ++ rl) { + for (int tl=0; tltimelevels; ++ tl) { + tt->set_time(ml, rl, tl, my_tt.get_time(ml, rl, tl)); + } + tt->set_delta(ml, rl, my_tt.get_delta(ml, rl)); + } + } + + PostRegrid(cctkGH); + + // Ensure we are in global mode + assert(is_global_mode()); + global_time = my_global_time; + (const_cast(cctkGH))->cctk_time = global_time; + delta_time = my_delta_time; + timereflevelfact = timereffacts.AT(reflevels - 1); + spacereflevelfact = ivect(-get_deadbeef()); + ivect::ref(cctkGH->cctk_levfac) = spacereflevelfact; + (const_cast(cctkGH))->cctk_timefac = timereflevelfact; + + // Recompose the grid structure (ignoring the previous content of + // the grid hierarchy) + for (int rl=0; rlghost_widths.size() == my_ghost_widths.AT(m).size()); + for (int rl=0; rlghost_widths.AT(rl) == + my_ghost_widths.AT(m).AT(rl)))); + } + assert(vdd.AT(m)->buffer_widths.size() == my_buffer_widths.AT(m).size()); + for (int rl=0; rlbuffer_widths.AT(rl) == + my_buffer_widths.AT(m).AT(rl)))); + } + assert(vdd.AT(m)->overlap_widths.size() == + my_overlap_widths.AT(m).size()); + for (int rl=0; rloverlap_widths.AT(rl) == + my_overlap_widths.AT(m).AT(rl)))); + } + for (int rl=0; rlprolongation_orders_space.AT(rl) == + my_prolongation_orders_space.AT(rl)); + } + } + } + } // end namespace CarpetIOF5 -- cgit v1.2.3