From f007420ae6d59572ffec2ecaa20ab807f57a94ea Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Wed, 28 Jul 2021 18:59:54 +0000 Subject: [PATCH 01/17] Add support for fully local HDF5 files and shared dumping of meep::structure --- src/h5file.cpp | 31 +++--- src/meep.hpp | 21 +++- src/structure_dump.cpp | 219 +++++++++++++++++++++++------------------ 3 files changed, 160 insertions(+), 111 deletions(-) diff --git a/src/h5file.cpp b/src/h5file.cpp index 484640942..e2f511583 100644 --- a/src/h5file.cpp +++ b/src/h5file.cpp @@ -150,7 +150,7 @@ void h5file::close_id() { /* note: if parallel is true, then *all* processes must call this, and all processes will use I/O. */ -h5file::h5file(const char *filename_, access_mode m, bool parallel_) { +h5file::h5file(const char *filename_, access_mode m, bool parallel_, bool local_) { cur_dataname = NULL; id = (void *)malloc(sizeof(hid_t)); cur_id = (void *)malloc(sizeof(hid_t)); @@ -160,7 +160,12 @@ h5file::h5file(const char *filename_, access_mode m, bool parallel_) { filename = new char[strlen(filename_) + 1]; strcpy(filename, filename_); mode = m; + + if (parallel_ && local_) { + meep::abort("Can not open h5file (%s) in both parallel and local mode.", filename); + } parallel = parallel_; + local = local_; } h5file::~h5file() { @@ -191,7 +196,9 @@ void h5file::remove() { extending = 0; IF_EXCLUSIVE(if (parallel) all_wait(), (void)0); - if (am_master() && std::remove(filename)) meep::abort("error removing file %s", filename); + if (am_master() || local) { + if (std::remove(filename)) meep::abort("error removing file %s", filename); + } } h5file::extending_s *h5file::get_extending(const char *dataname) const { @@ -226,7 +233,7 @@ void h5file::set_cur(const char *dataname, void *data_id) { void h5file::read_size(const char *dataname, int *rank, size_t *dims, int maxrank) { #ifdef HAVE_HDF5 - if (parallel || am_master()) { + if (parallel || am_master() || local) { hid_t file_id = HID(get_id()), space_id, data_id; CHECK(file_id >= 0, "error opening HDF5 input file"); @@ -291,7 +298,7 @@ void *h5file::read(const char *dataname, int *rank, size_t *dims, int maxrank, bool single_precision) { #ifdef HAVE_HDF5 void *data = 0; - if (parallel || am_master()) { + if (parallel || am_master() || local) { int i, N; hid_t file_id = HID(get_id()), space_id, data_id; @@ -345,7 +352,7 @@ void *h5file::read(const char *dataname, int *rank, size_t *dims, int maxrank, if (close_data_id) H5Dclose(data_id); } - if (!parallel) { + if (!parallel && !local) { *rank = broadcast(0, *rank); broadcast(0, dims, *rank); size_t N = 1; @@ -375,7 +382,7 @@ char *h5file::read(const char *dataname) { #ifdef HAVE_HDF5 char *data = 0; int len = 0; - if (parallel || am_master()) { + if (parallel || am_master() || local) { hid_t file_id = HID(get_id()), space_id, data_id, type_id; CHECK(file_id >= 0, "error opening HDF5 input file"); @@ -404,7 +411,7 @@ char *h5file::read(const char *dataname) { H5Dclose(data_id); } - if (!parallel) { + if (!parallel && !local) { len = broadcast(0, len); if (!am_master()) data = new char[len]; broadcast(0, data, len); @@ -442,7 +449,7 @@ void h5file::remove_data(const char *dataname) { if (dataset_exists(dataname)) { /* this is hackish ...need to pester HDF5 developers to make H5Gunlink a collective operation for parallel mode */ - if (!parallel || am_master()) { + if (!parallel || am_master() || local) { H5Gunlink(file_id, dataname); /* delete it */ H5Fflush(file_id, H5F_SCOPE_GLOBAL); } @@ -471,7 +478,7 @@ void h5file::create_data(const char *dataname, int rank, const size_t *dims, boo unset_cur(); remove_data(dataname); // HDF5 gives error if we H5Dcreate existing dataset - if (IF_EXCLUSIVE(!parallel || am_master(), 1)) { + if (local || IF_EXCLUSIVE(!parallel || am_master(), 1)) { hsize_t *dims_copy = new hsize_t[rank1 + append_data]; hsize_t *maxdims = new hsize_t[rank1 + append_data]; hsize_t N = 1; @@ -708,12 +715,12 @@ void h5file::done_writing_chunks() { I'm assuming(?) that non-extensible datasets will use different files, etcetera, for different timesteps. All of this hackery goes away if we just use an MPI-compiled version of HDF5. */ - if (parallel && cur_dataname && get_extending(cur_dataname)) prevent_deadlock(); // closes id + if (parallel && !local && cur_dataname && get_extending(cur_dataname)) prevent_deadlock(); // closes id } void h5file::write(const char *dataname, int rank, const size_t *dims, void *data, bool single_precision) { - if (parallel || am_master()) { + if (parallel || am_master() || local) { size_t *start = new size_t[rank + 1]; for (int i = 0; i < rank; i++) start[i] = 0; @@ -732,7 +739,7 @@ void h5file::write(const char *dataname, int rank, const size_t *dims, void *dat void h5file::write(const char *dataname, const char *data) { #ifdef HAVE_HDF5 - if (IF_EXCLUSIVE(am_master(), parallel || am_master())) { + if (local || IF_EXCLUSIVE(am_master(), parallel || am_master())) { hid_t file_id = HID(get_id()), type_id, data_id, space_id; CHECK(file_id >= 0, "error opening HDF5 output file"); diff --git a/src/meep.hpp b/src/meep.hpp index 7bfcbe13b..8d0ecaa63 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -384,7 +384,13 @@ class h5file { public: typedef enum { READONLY, READWRITE, WRITE } access_mode; - h5file(const char *filename_, access_mode m = READWRITE, bool parallel_ = true); + // If 'parallel_' is true, then we assume that all processes will be doing + // I/O, else we assume that *only* the master is doing I/O and all other + // processes will send/receive data to/from the master. + // If 'local_' is true, then 'parallel_' *must* be false and assumes that + // each process is writing to a local non-shared file and the filename is + // unique to the process. + h5file(const char *filename_, access_mode m = READWRITE, bool parallel_ = true, bool local_ = false); ~h5file(); // closes the files (and any open dataset) bool ok(); @@ -423,6 +429,7 @@ class h5file { access_mode mode; char *filename; bool parallel; + bool local; bool is_cur(const char *dataname); void unset_cur(); @@ -849,9 +856,15 @@ class structure { std::vector get_chunk_owners() const; // structure_dump.cpp - void dump(const char *filename); + // Dump structure to specified file. If 'single_parallel_file' + // is 'true' (the default) - then all processes write to the same/single file + // file after computing their respective offsets into this file. When set to + // 'false', each process writes data for the chunks it owns to a separate + // (process unique) file. + void dump(const char *filename, bool single_parallel_file=true); + void load(const char *filename, bool single_parallel_file=true); + void dump_chunk_layout(const char *filename); - void load(const char *filename); void load_chunk_layout(const char *filename, boundary_region &br); void load_chunk_layout(const std::vector &gvs, const std::vector &ids, @@ -889,7 +902,7 @@ class structure { void changing_chunks(); // Helper methods for dumping and loading susceptibilities void set_chiP_from_file(h5file *file, const char *dataset, field_type ft); - void write_susceptibility_params(h5file *file, const char *dname, int EorH); + void write_susceptibility_params(h5file *file, bool single_parallel_file, const char *dname, int EorH); std::unique_ptr bp; }; diff --git a/src/structure_dump.cpp b/src/structure_dump.cpp index c5bcca7c1..5cdf1d118 100644 --- a/src/structure_dump.cpp +++ b/src/structure_dump.cpp @@ -33,7 +33,9 @@ namespace meep { // Write the parameters required to reconstruct the susceptibility (id, noise_amp (for noisy), // omega_0, gamma, no_omega_0_denominator) -void structure::write_susceptibility_params(h5file *file, const char *dname, int EorH) { +void structure::write_susceptibility_params(h5file *file, + bool single_parallel_file, + const char *dname, int EorH) { // Get number of susceptibility params from first chunk, since all chunks will have // the same susceptibility list. size_t params_ntotal = 0; @@ -47,7 +49,7 @@ void structure::write_susceptibility_params(h5file *file, const char *dname, int size_t params_start = 0; file->create_data(dname, 1, ¶ms_ntotal, false /* append_data */, sizeof(realnum) == sizeof(float) /* single_precision */); - if (am_master()) { + if (am_master() || !single_parallel_file) { susceptibility *sus = chunks[0]->chiP[EorH]; while (sus) { sus->dump_params(file, ¶ms_start); @@ -85,36 +87,60 @@ void structure::dump_chunk_layout(const char *filename) { delete[] nums; } -void structure::dump(const char *filename) { +void structure::dump(const char *filename, bool single_parallel_file) { if (verbosity > 0) master_printf("creating epsilon output file \"%s\"...\n", filename); - // make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting - // the number of entries in the chi1inv array for each chunk. - size_t *num_chi1inv_ = new size_t[num_chunks * NUM_FIELD_COMPONENTS * 5]; - memset(num_chi1inv_, 0, sizeof(size_t) * size_t(num_chunks * NUM_FIELD_COMPONENTS * 5)); + /* + * make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting + * the number of entries in the chi1inv array for each chunk. + * + * When 'single_parallel_file' is true, we are creating a single block of data + * for ALL chunks (that are merged using MPI). Otherwise, we are just + * making a copy of just the chunks that are ours. + */ + int my_num_chunks = 0; + for (int i = 0; i < num_chunks; i++) { + my_num_chunks += (single_parallel_file || chunks[i]->is_mine()); + } + size_t num_chi1inv_size = my_num_chunks * NUM_FIELD_COMPONENTS * 5; + std::vector num_chi1inv_(num_chi1inv_size); size_t my_ntot = 0; - for (int i = 0; i < num_chunks; i++) + for (int i = 0, chunk_i = 0; i < num_chunks; i++) { if (chunks[i]->is_mine()) { size_t ntot = chunks[i]->gv.ntot(); - for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) - for (int d = 0; d < 5; ++d) + for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { + for (int d = 0; d < 5; ++d) { if (chunks[i]->chi1inv[c][d]) - my_ntot += (num_chi1inv_[(i * NUM_FIELD_COMPONENTS + c) * 5 + d] = ntot); + my_ntot += (num_chi1inv_[(chunk_i * NUM_FIELD_COMPONENTS + c) * 5 + d] = ntot); + } + } } - size_t *num_chi1inv = new size_t[num_chunks * NUM_FIELD_COMPONENTS * 5]; - sum_to_master(num_chi1inv_, num_chi1inv, num_chunks * NUM_FIELD_COMPONENTS * 5); - delete[] num_chi1inv_; + chunk_i += (chunks[i]->is_mine() || single_parallel_file); + } + + std::vector num_chi1inv; + if (single_parallel_file) { + num_chi1inv.resize(num_chi1inv_size); + sum_to_master(num_chi1inv_.data(), num_chi1inv.data(), num_chi1inv_size); + } else { + num_chi1inv = std::move(num_chi1inv_); + } // determine total dataset size and offset of this process's data - size_t my_start = partial_sum_to_all(my_ntot) - my_ntot; - size_t ntotal = sum_to_all(my_ntot); + size_t my_start = 0; + size_t ntotal = my_ntot; + if (single_parallel_file) { + my_start = partial_sum_to_all(my_ntot) - my_ntot; + ntotal = sum_to_all(my_ntot); + } h5file file(filename, h5file::WRITE, true); - size_t dims[3] = {(size_t)num_chunks, NUM_FIELD_COMPONENTS, 5}; + size_t dims[3] = {(size_t)my_num_chunks, NUM_FIELD_COMPONENTS, 5}; size_t start[3] = {0, 0, 0}; file.create_data("num_chi1inv", 3, dims); - if (am_master()) file.write_chunk(3, start, dims, num_chi1inv); - delete[] num_chi1inv; + if (am_master() || !single_parallel_file) { + file.write_chunk(3, start, dims, num_chi1inv.data()); + } // write the data file.create_data("chi1inv", 1, &ntotal, false /* append_data */, false /* single_precision */); @@ -147,7 +173,7 @@ void structure::dump(const char *filename) { // Write the number of susceptibilites size_t len = 2; file.create_data("num_sus", 1, &len); - if (am_master()) { + if (am_master() || !single_parallel_file) { size_t start = 0; size_t ntot = 2; file.write_chunk(1, &start, &ntot, num_sus); @@ -157,13 +183,11 @@ void structure::dump(const char *filename) { // Get number of non-null sigma entries for each chiP in each chunk. // Assumes each susceptibility in the chiP[E_stuff] list has the // same number of non-null sigma elements. Likewise for chiP[H_stuff] - size_t *my_num_sigmas[2]; - size_t *num_sigmas[2]; + std::vector my_num_sigmas[2]; + std::vector num_sigmas[2]; for (int ft = 0; ft < 2; ++ft) { - my_num_sigmas[ft] = new size_t[num_chunks]; - num_sigmas[ft] = new size_t[num_chunks]; - memset(my_num_sigmas[ft], 0, sizeof(size_t) * num_chunks); - memset(num_sigmas[ft], 0, sizeof(size_t) * num_chunks); + my_num_sigmas[ft].resize(num_chunks); + num_sigmas[ft].resize(num_chunks); } for (int i = 0; i < num_chunks; ++i) { @@ -199,7 +223,7 @@ void structure::dump(const char *filename) { file.prevent_deadlock(); for (int ft = 0; ft < 2; ++ft) { - sum_to_all(my_num_sigmas[ft], num_sigmas[ft], num_chunks); + sum_to_all(my_num_sigmas[ft].data(), num_sigmas[ft].data(), num_chunks); } size_t num_E_sigmas = 0; @@ -210,17 +234,13 @@ void structure::dump(const char *filename) { } // Allocate space for component and direction of non-null sigmas - size_t *my_sigma_cd[2] = {NULL, NULL}; - my_sigma_cd[E_stuff] = new size_t[num_E_sigmas * 2]; - memset(my_sigma_cd[E_stuff], 0, sizeof(size_t) * num_E_sigmas * 2); - my_sigma_cd[H_stuff] = new size_t[num_H_sigmas * 2]; - memset(my_sigma_cd[H_stuff], 0, sizeof(size_t) * num_H_sigmas * 2); - - size_t *sigma_cd[2] = {NULL, NULL}; - sigma_cd[E_stuff] = new size_t[num_E_sigmas * 2]; - memset(sigma_cd[E_stuff], 0, sizeof(size_t) * num_E_sigmas * 2); - sigma_cd[H_stuff] = new size_t[num_H_sigmas * 2]; - memset(sigma_cd[H_stuff], 0, sizeof(size_t) * num_H_sigmas * 2); + std::vector my_sigma_cd[2]; + my_sigma_cd[E_stuff].resize(num_E_sigmas * 2); + my_sigma_cd[H_stuff].resize(num_H_sigmas * 2); + + std::vector sigma_cd[2]; + sigma_cd[E_stuff].resize(num_E_sigmas * 2); + sigma_cd[H_stuff].resize(num_H_sigmas * 2); // Find component and direction of non-null sigmas { @@ -246,8 +266,8 @@ void structure::dump(const char *filename) { } } } - bw_or_to_all(my_sigma_cd[E_stuff], sigma_cd[E_stuff], num_E_sigmas * 2); - bw_or_to_all(my_sigma_cd[H_stuff], sigma_cd[H_stuff], num_H_sigmas * 2); + bw_or_to_all(my_sigma_cd[E_stuff].data(), sigma_cd[E_stuff].data(), num_E_sigmas * 2); + bw_or_to_all(my_sigma_cd[H_stuff].data(), sigma_cd[H_stuff].data(), num_H_sigmas * 2); } // Write location (component and direction) data of non-null sigmas (sigma[c][d]) @@ -256,9 +276,9 @@ void structure::dump(const char *filename) { file.create_data("sigma_cd", 1, &len); size_t start = 0; for (int ft = 0; ft < 2; ++ft) { - if (am_master()) { - size_t count = (ft == 0 ? num_E_sigmas : num_H_sigmas) * 2; - file.write_chunk(1, &start, &count, sigma_cd[ft]); + if (am_master() || !single_parallel_file) { + size_t count = sigma_cd[ft].size(); + file.write_chunk(1, &start, &count, sigma_cd[ft].data()); start += count; } } @@ -292,15 +312,8 @@ void structure::dump(const char *filename) { } } - write_susceptibility_params(&file, "E_params", E_stuff); - write_susceptibility_params(&file, "H_params", H_stuff); - - for (int ft = 0; ft < 2; ++ft) { - delete[] sigma_cd[ft]; - delete[] my_sigma_cd[ft]; - delete[] num_sigmas[ft]; - delete[] my_num_sigmas[ft]; - } + write_susceptibility_params(&file, single_parallel_file, "E_params", E_stuff); + write_susceptibility_params(&file, single_parallel_file, "H_params", H_stuff); } // Reconstruct the chiP lists of susceptibilities from the params hdf5 data @@ -542,35 +555,49 @@ void structure::load_chunk_layout(const std::vector &gvs, check_chunks(); } -void structure::load(const char *filename) { +void structure::load(const char *filename, bool single_parallel_file) { h5file file(filename, h5file::READONLY, true); if (verbosity > 0) master_printf("reading epsilon from file \"%s\"...\n", filename); - // make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting - // the number of entries in the chi1inv array for each chunk. - size_t *num_chi1inv = new size_t[num_chunks * NUM_FIELD_COMPONENTS * 5]; + /* + * make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting + * the number of entries in the chi1inv array for each chunk. + * + * When 'single_parallel_file' is true, we are creating a single block of data + * for ALL chunks (that are merged using MPI). Otherwise, we are just + * making a copy of just the chunks that are ours. + */ + int my_num_chunks = 0; + for (int i = 0; i < num_chunks; i++) { + my_num_chunks += (single_parallel_file || chunks[i]->is_mine()); + } + size_t num_chi1inv_size = my_num_chunks * NUM_FIELD_COMPONENTS * 5; + std::vector num_chi1inv(num_chi1inv_size); + int rank; - size_t dims[3], _dims[3] = {(size_t)num_chunks, NUM_FIELD_COMPONENTS, 5}; + size_t dims[3], _dims[3] = {(size_t)my_num_chunks, NUM_FIELD_COMPONENTS, 5}; size_t start[3] = {0, 0, 0}; file.read_size("num_chi1inv", &rank, dims, 3); if (rank != 3 || _dims[0] != dims[0] || _dims[1] != dims[1] || _dims[2] != dims[2]) meep::abort("chunk mismatch in structure::load"); - if (am_master()) file.read_chunk(3, start, dims, num_chi1inv); + if (am_master() || !single_parallel_file) file.read_chunk(3, start, dims, num_chi1inv.data()); - file.prevent_deadlock(); - broadcast(0, num_chi1inv, dims[0] * dims[1] * dims[2]); + if (single_parallel_file) { + file.prevent_deadlock(); + broadcast(0, num_chi1inv.data(), dims[0] * dims[1] * dims[2]); + } changing_chunks(); // allocate data as needed and check sizes size_t my_ntot = 0; - for (int i = 0; i < num_chunks; i++) + for (int i = 0, chunk_i = 0; i < num_chunks; i++) { if (chunks[i]->is_mine()) { size_t ntot = chunks[i]->gv.ntot(); for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) for (int d = 0; d < 5; ++d) { - size_t n = num_chi1inv[(i * NUM_FIELD_COMPONENTS + c) * 5 + d]; + size_t n = num_chi1inv[(chunk_i * NUM_FIELD_COMPONENTS + c) * 5 + d]; if (n == 0) { delete[] chunks[i]->chi1inv[c][d]; chunks[i]->chi1inv[c][d] = NULL; @@ -582,10 +609,16 @@ void structure::load(const char *filename) { } } } + chunk_i += (chunks[i]->is_mine() || single_parallel_file); + } // determine total dataset size and offset of this process's data - size_t my_start = partial_sum_to_all(my_ntot) - my_ntot; - size_t ntotal = sum_to_all(my_ntot); + size_t my_start = 0; + size_t ntotal = my_ntot; + if (single_parallel_file) { + my_start = partial_sum_to_all(my_ntot) - my_ntot; + ntotal = sum_to_all(my_ntot); + } // read the data file.read_size("chi1inv", &rank, dims, 1); @@ -611,22 +644,22 @@ void structure::load(const char *filename) { size_t dims[] = {0, 0, 0}; file.read_size("num_sus", &rank, dims, 1); if (dims[0] > 0) { - if (am_master()) { + if (am_master() || !single_parallel_file) { size_t start = 0; size_t count = 2; file.read_chunk(rank, &start, &count, num_sus); } } - file.prevent_deadlock(); - broadcast(0, num_sus, 2); + if (single_parallel_file) { + file.prevent_deadlock(); + broadcast(0, num_sus, 2); + } } // Allocate and read non-null sigma entry data - size_t *num_sigmas[2] = {NULL, NULL}; - num_sigmas[E_stuff] = new size_t[num_chunks]; - memset(num_sigmas[E_stuff], 0, sizeof(size_t) * num_chunks); - num_sigmas[H_stuff] = new size_t[num_chunks]; - memset(num_sigmas[H_stuff], 0, sizeof(size_t) * num_chunks); + std::vector num_sigmas[2]; + num_sigmas[E_stuff].resize(num_chunks); + num_sigmas[H_stuff].resize(num_chunks); // Read num_sigmas data { @@ -634,21 +667,21 @@ void structure::load(const char *filename) { size_t dims[] = {0, 0, 0}; file.read_size("num_sigmas", &rank, dims, 1); if (dims[0] != (size_t)num_chunks * 2) { meep::abort("inconsistent data size in structure::load"); } - if (am_master()) { + if (am_master() || !single_parallel_file) { size_t start = 0; size_t count = num_chunks; - file.read_chunk(rank, &start, &count, num_sigmas[E_stuff]); + file.read_chunk(rank, &start, &count, num_sigmas[E_stuff].data()); start += count; - file.read_chunk(rank, &start, &count, num_sigmas[H_stuff]); + file.read_chunk(rank, &start, &count, num_sigmas[H_stuff].data()); + } + if (single_parallel_file) { + file.prevent_deadlock(); + broadcast(0, num_sigmas[E_stuff].data(), num_chunks); + broadcast(0, num_sigmas[H_stuff].data(), num_chunks); } - file.prevent_deadlock(); - broadcast(0, num_sigmas[E_stuff], num_chunks); - broadcast(0, num_sigmas[H_stuff], num_chunks); } // Allocate space for component and direction data of the non-null susceptibilities - size_t *sigma_cd[2] = {NULL, NULL}; - size_t nsig[2] = {0, 0}; for (int ft = 0; ft < 2; ++ft) { for (int i = 0; i < num_chunks; ++i) { @@ -656,10 +689,9 @@ void structure::load(const char *filename) { } } - sigma_cd[E_stuff] = new size_t[nsig[E_stuff] * 2]; - memset(sigma_cd[E_stuff], 0, sizeof(size_t) * nsig[E_stuff] * 2); - sigma_cd[H_stuff] = new size_t[nsig[H_stuff] * 2]; - memset(sigma_cd[H_stuff], 0, sizeof(size_t) * nsig[H_stuff] * 2); + std::vector sigma_cd[2]; + sigma_cd[E_stuff].resize(nsig[E_stuff] * 2); + sigma_cd[H_stuff].resize(nsig[H_stuff] * 2); // Read the component/direction data { @@ -670,17 +702,19 @@ void structure::load(const char *filename) { meep::abort("inconsistent data size in structure::load"); } - if (am_master()) { + if (am_master() || !single_parallel_file) { size_t start = 0; for (int ft = 0; ft < 2; ++ft) { - size_t count = nsig[ft] * 2; - file.read_chunk(rank, &start, &count, sigma_cd[ft]); + size_t count = sigma_cd[ft].size(); + file.read_chunk(rank, &start, &count, sigma_cd[ft].data()); start += count; } } - file.prevent_deadlock(); - broadcast(0, sigma_cd[E_stuff], nsig[E_stuff] * 2); - broadcast(0, sigma_cd[H_stuff], nsig[H_stuff] * 2); + if (single_parallel_file) { + file.prevent_deadlock(); + broadcast(0, sigma_cd[E_stuff].data(), nsig[E_stuff] * 2); + broadcast(0, sigma_cd[H_stuff].data(), nsig[H_stuff] * 2); + } } for (int ft = 0; ft < 2; ++ft) { @@ -730,10 +764,5 @@ void structure::load(const char *filename) { } } } - - for (int ft = 0; ft < 2; ++ft) { - delete[] num_sigmas[ft]; - delete[] sigma_cd[ft]; - } } } // namespace meep From 344fe551f7d7cd2571a8df545873f2f65356147b Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 29 Jul 2021 17:11:01 +0000 Subject: [PATCH 02/17] Add support for fully local HDF5 files and shared dumping of meep::structure --- python/simulation.py | 53 +++++++++++++++++++++++++-------- python/tests/test_simulation.py | 19 +++++++----- src/h5file.cpp | 2 +- src/structure_dump.cpp | 21 ++++++++----- 4 files changed, 67 insertions(+), 28 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index e8ef07d83..6bbb0de09 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -959,7 +959,6 @@ def __init__(self, filename_prefix=None, output_volume=None, output_single_precision=False, - load_structure='', geometry_center=mp.Vector3(), force_all_components=False, split_chunks_evenly=True, @@ -1157,11 +1156,6 @@ def __init__(self, that are specified by geometric objects. You should list any materials other than scalar dielectrics that are returned by `material_function` here. - + **`load_structure` [`string`]** — If not empty, Meep will load the structure - file specified by this string. The file must have been created by - `mp.dump_structure`. Defaults to an empty string. See [Load and Dump - Structure](#load-and-dump-structure) for more information. - + **`chunk_layout` [`string` or `Simulation` instance or `BinaryPartition` class]** — This will cause the `Simulation` to use the chunk layout described by either (1) an `.h5` file (created using `Simulation.dump_chunk_layout`), (2) another @@ -1233,7 +1227,6 @@ def __init__(self, self.is_cylindrical = False self.material_function = material_function self.epsilon_func = epsilon_func - self.load_structure_file = load_structure self.dft_objects = [] self._is_initialized = False self.force_all_components = force_all_components @@ -1244,6 +1237,9 @@ def __init__(self, self.fragment_stats = None self._output_stats = os.environ.get('MEEP_STATS', None) + self.load_single_parallel_file = True + self.load_structure_file = None + self.special_kz = False if self.cell_size.z == 0 and self.k_point and self.k_point.z != 0: if kz_2d == "complex": @@ -1715,7 +1711,8 @@ def _init_structure(self, k=False): self.num_chunks = self.chunk_layout.numchunks() if self.load_structure_file: - self.load_structure(self.load_structure_file) + self.load_structure( + self.load_structure_file, self.load_single_parallel_file) def _is_outer_boundary(self, vol, direction, side): @@ -1859,22 +1856,22 @@ def set_materials(self, geometry=None, default_material=None): None ) - def dump_structure(self, fname): + def dump_structure(self, fname, single_parallel_file=True): """ Dumps the structure to the file `fname`. """ if self.structure is None: raise ValueError("Fields must be initialized before calling dump_structure") - self.structure.dump(fname) + self.structure.dump(fname, single_parallel_file) - def load_structure(self, fname): + def load_structure(self, fname, single_parallel_file=True): """ Loads a structure from the file `fname`. A file name to load can also be passed to the `Simulation` constructor via the `load_structure` keyword argument. """ if self.structure is None: raise ValueError("Fields must be initialized before calling load_structure") - self.structure.load(fname) + self.structure.load(fname, single_parallel_file) def dump_chunk_layout(self, fname): """ @@ -1896,6 +1893,38 @@ def load_chunk_layout(self, br, source): ## source is either filename (string) self.structure.load_chunk_layout(source, br) + def _get_load_dump_dirname(self, dirname, single_parallel_file): + """ + """ + if single_parallel_file: + dump_dirname = dirname + else: + dump_dirname = os.path.join(dirname, 'rank%02d' % mp.my_rank()) + return dump_dirname + + def dump(self, dirname, structure=True, single_parallel_file=True): + """ + Dumps simulation state. + """ + dump_dirname = self._get_load_dump_dirname(dirname, single_parallel_file) + os.makedirs(dump_dirname, exist_ok=True) + + if structure: + structure_dump_filename = os.path.join(dump_dirname, 'structure.h5') + self.dump_structure(structure_dump_filename, single_parallel_file) + + def load(self, dirname, structure=True, single_parallel_file=True): + """ + Loads simulation state. + + This should called right after the Simulation object has been created + but before 'init_sim' is called. + """ + dump_dirname = self._get_load_dump_dirname(dirname, single_parallel_file) + self.load_single_parallel_file = single_parallel_file + if structure: + self.load_structure_file = os.path.join(dump_dirname, 'structure.h5') + def init_sim(self): if self._is_initialized: return diff --git a/python/tests/test_simulation.py b/python/tests/test_simulation.py index d5a9799da..b13cc2e68 100644 --- a/python/tests/test_simulation.py +++ b/python/tests/test_simulation.py @@ -299,7 +299,7 @@ def test_in_box_volumes(self): sim.field_energy_in_box(tv) sim.field_energy_in_box(v) - def _load_dump_structure(self, chunk_file=False, chunk_sim=False): + def _load_dump_structure(self, chunk_file=False, chunk_sim=False, single_parallel_file=True): from meep.materials import Al resolution = 50 cell = mp.Vector3(5, 5) @@ -326,12 +326,14 @@ def get_ref_field_point(sim): ref_field_points.append(p.real) sim1.run(mp.at_every(5, get_ref_field_point), until=50) - dump_fn = os.path.join(self.temp_dir, 'test_load_dump_structure.h5') + + dump_dirname = os.path.join(self.temp_dir, 'test_load_dump') + sim1.dump(dump_dirname, structure=True, single_parallel_file=single_parallel_file) + dump_chunk_fname = None chunk_layout = None - sim1.dump_structure(dump_fn) if chunk_file: - dump_chunk_fname = os.path.join(self.temp_dir, 'test_load_dump_structure_chunks.h5') + dump_chunk_fname = os.path.join(dump_dirname, 'chunk_layout.h5') sim1.dump_chunk_layout(dump_chunk_fname) chunk_layout = dump_chunk_fname if chunk_sim: @@ -342,9 +344,8 @@ def get_ref_field_point(sim): boundary_layers=pml_layers, sources=[sources], symmetries=symmetries, - chunk_layout=chunk_layout, - load_structure=dump_fn) - + chunk_layout=chunk_layout) + sim.load(dump_dirname, structure=True, single_parallel_file=single_parallel_file) field_points = [] def get_field_point(sim): @@ -359,6 +360,10 @@ def get_field_point(sim): def test_load_dump_structure(self): self._load_dump_structure() + @unittest.skipIf(not mp.with_mpi(), "MPI specific test") + def test_load_dump_structure_sharded(self): + self._load_dump_structure(single_parallel_file=False) + def test_load_dump_chunk_layout_file(self): self._load_dump_structure(chunk_file=True) diff --git a/src/h5file.cpp b/src/h5file.cpp index e2f511583..366e5d63c 100644 --- a/src/h5file.cpp +++ b/src/h5file.cpp @@ -260,7 +260,7 @@ void h5file::read_size(const char *dataname, int *rank, size_t *dims, int maxran H5Sclose(space_id); } - if (!parallel) { + if (!parallel && !local) { *rank = broadcast(0, *rank); broadcast(0, dims, *rank); diff --git a/src/structure_dump.cpp b/src/structure_dump.cpp index 5cdf1d118..8e19a5812 100644 --- a/src/structure_dump.cpp +++ b/src/structure_dump.cpp @@ -88,7 +88,7 @@ void structure::dump_chunk_layout(const char *filename) { } void structure::dump(const char *filename, bool single_parallel_file) { - if (verbosity > 0) master_printf("creating epsilon output file \"%s\"...\n", filename); + if (verbosity > 0) printf("creating epsilon from file \"%s\" (%d)...\n", filename, single_parallel_file); /* * make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting @@ -134,7 +134,7 @@ void structure::dump(const char *filename, bool single_parallel_file) { ntotal = sum_to_all(my_ntot); } - h5file file(filename, h5file::WRITE, true); + h5file file(filename, h5file::WRITE, single_parallel_file, !single_parallel_file); size_t dims[3] = {(size_t)my_num_chunks, NUM_FIELD_COMPONENTS, 5}; size_t start[3] = {0, 0, 0}; file.create_data("num_chi1inv", 3, dims); @@ -419,7 +419,7 @@ void structure::set_chiP_from_file(h5file *file, const char *dataset, field_type size_t dims[3] = {0, 0, 0}; file->read_size(dataset, &rank, dims, 1); - if (rank != 1) meep::abort("inconsistent data size in structure::load"); + if (rank != 1) meep::abort("inconsistent data size in structure::set_chiP_from_file"); if (dims[0] != 0) { for (int i = 0; i < num_chunks; ++i) { @@ -556,9 +556,9 @@ void structure::load_chunk_layout(const std::vector &gvs, } void structure::load(const char *filename, bool single_parallel_file) { - h5file file(filename, h5file::READONLY, true); + h5file file(filename, h5file::READONLY, single_parallel_file, !single_parallel_file); - if (verbosity > 0) master_printf("reading epsilon from file \"%s\"...\n", filename); + if (verbosity > 0) printf("reading epsilon from file \"%s\" (%d)...\n", filename, single_parallel_file); /* * make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting @@ -622,7 +622,12 @@ void structure::load(const char *filename, bool single_parallel_file) { // read the data file.read_size("chi1inv", &rank, dims, 1); - if (rank != 1 || dims[0] != ntotal) meep::abort("inconsistent data size in structure::load"); + if (rank != 1 || dims[0] != ntotal) { + meep::abort( + "inconsistent data size for chi1inv in structure::load (rank, dims[0]): " + "(%d, %zu) != (1, %zu)", + rank, dims[0], ntotal); + } for (int i = 0; i < num_chunks; i++) if (chunks[i]->is_mine()) { size_t ntot = chunks[i]->gv.ntot(); @@ -666,7 +671,7 @@ void structure::load(const char *filename, bool single_parallel_file) { int rank = 0; size_t dims[] = {0, 0, 0}; file.read_size("num_sigmas", &rank, dims, 1); - if (dims[0] != (size_t)num_chunks * 2) { meep::abort("inconsistent data size in structure::load"); } + if (dims[0] != (size_t)num_chunks * 2) { meep::abort("inconsistent data size for num_sigmas in structure::load"); } if (am_master() || !single_parallel_file) { size_t start = 0; size_t count = num_chunks; @@ -699,7 +704,7 @@ void structure::load(const char *filename, bool single_parallel_file) { size_t dims[] = {0, 0, 0}; file.read_size("sigma_cd", &rank, dims, 1); if (dims[0] != 2 * (nsig[E_stuff] + nsig[H_stuff])) { - meep::abort("inconsistent data size in structure::load"); + meep::abort("inconsistent data size for sigma_cd in structure::load"); } if (am_master() || !single_parallel_file) { From a33c0623c5daca16f017bf0d3db0a14f796cf205 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 29 Jul 2021 17:49:01 +0000 Subject: [PATCH 03/17] Update python func docs --- python/simulation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/python/simulation.py b/python/simulation.py index 6bbb0de09..dc02f9a05 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1895,10 +1895,13 @@ def load_chunk_layout(self, br, source): def _get_load_dump_dirname(self, dirname, single_parallel_file): """ + Get the dirname to dump simulation state to. """ if single_parallel_file: dump_dirname = dirname else: + # When doing a sharded dump (each process to its own file), use + # the process rank to get a unique name. dump_dirname = os.path.join(dirname, 'rank%02d' % mp.my_rank()) return dump_dirname From 3b04bb1cca01ba5f84632f5e6df1a6854e70632a Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 12 Aug 2021 22:39:42 +0000 Subject: [PATCH 04/17] Update python API documentation --- doc/docs/Python_User_Interface.md | 59 ++++++++++++++++++++++------ doc/docs/Python_User_Interface.md.in | 2 + 2 files changed, 48 insertions(+), 13 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index ad351dc69..c141ccd9c 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -114,7 +114,6 @@ def __init__(self, filename_prefix=None, output_volume=None, output_single_precision=False, - load_structure='', geometry_center=Vector3<0.0, 0.0, 0.0>, force_all_components=False, split_chunks_evenly=True, @@ -315,11 +314,6 @@ Python. `Vector3` is a `meep` class. that are specified by geometric objects. You should list any materials other than scalar dielectrics that are returned by `material_function` here. -+ **`load_structure` [`string`]** — If not empty, Meep will load the structure - file specified by this string. The file must have been created by - `mp.dump_structure`. Defaults to an empty string. See [Load and Dump - Structure](#load-and-dump-structure) for more information. - + **`chunk_layout` [`string` or `Simulation` instance or `BinaryPartition` class]** — This will cause the `Simulation` to use the chunk layout described by either (1) an `.h5` file (created using `Simulation.dump_chunk_layout`), (2) another @@ -534,10 +528,29 @@ def mean_time_spent_on(self, time_sink):
Return the mean time spent by all processes for a type of work `time_sink` which -can be one of ten integer values `0`-`9`: (`0`) connecting chunks, (`1`) time stepping, -(`2`) copying boundaries, (`3`) MPI all-to-all communication/synchronization, -(`4`) MPI one-to-one communication, (`5`) field output, (`6`) Fourier transforming, -(`7`) MPB mode solver, (`8`) near-to-far field transformation, and (`9`) other. +can be one of the following integer constants: +* meep.Stepping ("time stepping") +* meep.Connecting ("connecting chunks") +* meep.Boundaries ("copying boundaries") +* meep.MpiAllTime ("all-all communication") +* meep.MpiOneTime ("1-1 communication") +* meep.FieldOutput ("outputting fields") +* meep.FourierTransforming ("Fourier transforming") +* meep.MPBTime ("MPB mode solver") +* meep.GetFarfieldsTime ("far-field transform") +* meep.FieldUpdateB ("updating B field") +* meep.FieldUpdateH ("updating H field") +* meep.FieldUpdateD ("updating D field") +* meep.FieldUpdateE ("updating E field") +* meep.BoundarySteppingB ("boundary stepping B") +* meep.BoundarySteppingWH ("boundary stepping WH") +* meep.BoundarySteppingPH ("boundary stepping PH") +* meep.BoundarySteppingH ("boundary stepping H") +* meep.BoundarySteppingD ("boundary stepping D") +* meep.BoundarySteppingWE ("boundary stepping WE") +* meep.BoundarySteppingPE ("boundary stepping PE") +* meep.BoundarySteppingE ("boundary stepping E") +* meep.Other ("everything else")
@@ -1009,7 +1022,9 @@ def reset_meep(self):
Reset all of Meep's parameters, deleting the fields, structures, etcetera, from -memory as if you had not run any computations. +memory as if you had not run any computations. If the num_chunks or chunk_layout +attributes have been modified internally, they are reset to their original +values passed in at instantiation.
@@ -2279,13 +2294,15 @@ is a 1d array of `nfreq` dimensions. These functions dump the raw ε and μ data to disk and load it back for doing multiple simulations with the same materials but different sources etc. The only prerequisite is that the dump/load simulations have the same [chunks](Chunks_and_Symmetry.md) (i.e. the same grid, number of processors, symmetries, and PML). When using `split_chunks_evenly=False`, you must also dump the original chunk layout using `dump_chunk_layout` and load it into the new `Simulation` using the `chunk_layout` parameter. Currently only stores dispersive and non-dispersive $\varepsilon$ and $\mu$ but not nonlinearities. Note that loading data from a file in this way overwrites any `geometry` data passed to the `Simulation` constructor. +For dump_structure and load_structure, when `single_parallel_file=True` (the default) - all chunks write to the same/single file after computing their respective offsets into this file. When set to 'False', each chunk writes writes its data to a separate file that is appropriately named by its chunk index. +
```python -def dump_structure(self, fname): +def dump_structure(self, fname, single_parallel_file=True): ```
@@ -2301,7 +2318,7 @@ Dumps the structure to the file `fname`.
```python -def load_structure(self, fname): +def load_structure(self, fname, single_parallel_file=True): ```
@@ -7366,6 +7383,22 @@ of processes.
+ + +
+ +```python +def print(self): +``` + +
+ +Pretty-prints the tree structure of the BinaryPartition object. + +
+ +
+ Miscellaneous Functions Reference --------------------------------- diff --git a/doc/docs/Python_User_Interface.md.in b/doc/docs/Python_User_Interface.md.in index ac31a14d7..91805d8d5 100644 --- a/doc/docs/Python_User_Interface.md.in +++ b/doc/docs/Python_User_Interface.md.in @@ -351,6 +351,8 @@ And this [`DftNear2Far`](#DftNear2Far) method: These functions dump the raw ε and μ data to disk and load it back for doing multiple simulations with the same materials but different sources etc. The only prerequisite is that the dump/load simulations have the same [chunks](Chunks_and_Symmetry.md) (i.e. the same grid, number of processors, symmetries, and PML). When using `split_chunks_evenly=False`, you must also dump the original chunk layout using `dump_chunk_layout` and load it into the new `Simulation` using the `chunk_layout` parameter. Currently only stores dispersive and non-dispersive $\varepsilon$ and $\mu$ but not nonlinearities. Note that loading data from a file in this way overwrites any `geometry` data passed to the `Simulation` constructor. +For dump_structure and load_structure, when `single_parallel_file=True` (the default) - all chunks write to the same/single file after computing their respective offsets into this file. When set to 'False', each chunk writes writes its data to a separate file that is appropriately named by its chunk index. + @@ Simulation.dump_structure @@ @@ Simulation.load_structure @@ @@ Simulation.dump_chunk_layout @@ From 0f11d4118b44c235b94db4fc1c944e3e0ebd71bf Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Mon, 23 Aug 2021 17:37:25 +0000 Subject: [PATCH 05/17] Dump/Load of 'fields' The PR is incomplete (does not include the dft fields) and also has a lot of debugging statements. --- python/Makefile.am | 11 +- python/simulation.py | 55 ++++- python/tests/test_dump_load.py | 107 +++++++++ python/tests/test_simulation.py | 71 ------ src/Makefile.am | 2 +- src/fields_dump.cpp | 244 +++++++++++++++++++++ src/meep.hpp | 17 ++ tests/Makefile.am | 15 +- tests/dump_load.cpp | 369 ++++++++++++++++++++++++++++++++ 9 files changed, 800 insertions(+), 91 deletions(-) create mode 100644 python/tests/test_dump_load.py create mode 100644 src/fields_dump.cpp create mode 100644 tests/dump_load.cpp diff --git a/python/Makefile.am b/python/Makefile.am index 1870a70ef..49c600769 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -63,14 +63,15 @@ TESTS = \ $(TEST_DIR)/test_get_point.py \ $(TEST_DIR)/test_holey_wvg_bands.py \ $(TEST_DIR)/test_holey_wvg_cavity.py \ - $(KDOM_TEST) \ + $(KDOM_TEST) \ $(TEST_DIR)/test_ldos.py \ - $(MDPYTEST) \ + $(TEST_DIR)/test_dump_load.py \ + $(MDPYTEST) \ $(TEST_DIR)/test_material_grid.py \ $(TEST_DIR)/test_medium_evaluations.py \ - $(MPBPYTEST) \ - $(MODE_COEFFS_TEST) \ - $(MODE_DECOMPOSITION_TEST) \ + $(MPBPYTEST) \ + $(MODE_COEFFS_TEST) \ + $(MODE_DECOMPOSITION_TEST) \ $(TEST_DIR)/test_multilevel_atom.py \ $(TEST_DIR)/test_n2f_periodic.py \ $(TEST_DIR)/test_oblique_source.py \ diff --git a/python/simulation.py b/python/simulation.py index aa059198b..1e338aa08 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1241,6 +1241,7 @@ def __init__(self, self.load_single_parallel_file = True self.load_structure_file = None + self.load_fields_file = None self.special_kz = False if self.cell_size.z == 0 and self.k_point and self.k_point.z != 0: @@ -1863,29 +1864,44 @@ def dump_structure(self, fname, single_parallel_file=True): Dumps the structure to the file `fname`. """ if self.structure is None: - raise ValueError("Fields must be initialized before calling dump_structure") + raise ValueError("Structure must be initialized before calling dump_structure") self.structure.dump(fname, single_parallel_file) def load_structure(self, fname, single_parallel_file=True): """ - Loads a structure from the file `fname`. A file name to load can also be passed to - the `Simulation` constructor via the `load_structure` keyword argument. + Loads a structure from the file `fname`. """ if self.structure is None: - raise ValueError("Fields must be initialized before calling load_structure") + raise ValueError("Structure must be initialized before loading structure from file '%s'" % fname) self.structure.load(fname, single_parallel_file) + def dump_fields(self, fname, single_parallel_file=True): + """ + Dumps the fields to the file `fname`. + """ + if self.fields is None: + raise ValueError("Fields must be initialized before calling dump_fields") + self.fields.dump(fname, single_parallel_file) + + def load_fields(self, fname, single_parallel_file=True): + """ + Loads a fields from the file `fname`. + """ + if self.fields is None: + raise ValueError("Fields must be initialized before loading fields from file '%s'" % fname) + self.fields.load(fname, single_parallel_file) + def dump_chunk_layout(self, fname): """ Dumps the chunk layout to file `fname`. """ if self.structure is None: - raise ValueError("Fields must be initialized before calling load_structure") + raise ValueError("Structure must be initialized before calling dump_chunk_layout") self.structure.dump_chunk_layout(fname) def load_chunk_layout(self, br, source): if self.structure is None: - raise ValueError("Fields must be initialized before calling load_structure") + raise ValueError("Structure must be initialized before loading chunk layout from file '%s'" % fname) if isinstance(source, Simulation): vols = source.structure.get_chunk_volumes() @@ -1907,7 +1923,7 @@ def _get_load_dump_dirname(self, dirname, single_parallel_file): dump_dirname = os.path.join(dirname, 'rank%02d' % mp.my_rank()) return dump_dirname - def dump(self, dirname, structure=True, single_parallel_file=True): + def dump(self, dirname, structure=True, fields=True, single_parallel_file=True): """ Dumps simulation state. """ @@ -1918,7 +1934,11 @@ def dump(self, dirname, structure=True, single_parallel_file=True): structure_dump_filename = os.path.join(dump_dirname, 'structure.h5') self.dump_structure(structure_dump_filename, single_parallel_file) - def load(self, dirname, structure=True, single_parallel_file=True): + if fields: + fields_dump_filename = os.path.join(dump_dirname, 'fields.h5') + self.dump_fields(fields_dump_filename, single_parallel_file) + + def load(self, dirname, structure=True, fields=True, single_parallel_file=True): """ Loads simulation state. @@ -1927,9 +1947,13 @@ def load(self, dirname, structure=True, single_parallel_file=True): """ dump_dirname = self._get_load_dump_dirname(dirname, single_parallel_file) self.load_single_parallel_file = single_parallel_file + if structure: self.load_structure_file = os.path.join(dump_dirname, 'structure.h5') + if fields: + self.load_fields_file = os.path.join(dump_dirname, 'fields.h5') + def init_sim(self): if self._is_initialized: return @@ -1958,6 +1982,12 @@ def init_sim(self): self.loop_tile_base ) + if self.load_fields_file: + print("Loading fields from file: %s (%s)" % (self.load_fields_file, str(self.load_single_parallel_file))) + self.load_fields( + self.load_fields_file, self.load_single_parallel_file) + print("LOADED fields from file: %s (%s)" % (self.load_fields_file, str(self.load_single_parallel_file))) + if self.force_all_components and self.dimensions != 1: self.fields.require_component(mp.Ez) self.fields.require_component(mp.Hz) @@ -1977,6 +2007,7 @@ def init_sim(self): hook() self._is_initialized = True + print("Sim INITIALIZED.") def using_real_fields(self): cond1 = self.is_cylindrical and self.m != 0 @@ -2220,6 +2251,7 @@ def _run_until(self, cond, step_funcs): if isinstance(cond[i], numbers.Number): stop_time = cond[i] t0 = self.round_time() + print('_run_until: stop_time = %s, t0 = %s' % (str(stop_time), str(t0))) def stop_cond(sim): return sim.round_time() >= t0 + stop_time @@ -2237,6 +2269,7 @@ def stop_cond(sim): while not any([x(self) for x in cond]): for func in step_funcs: + print('_run_until: _eval_step_func: ' + str(func)) _eval_step_func(self, func, 'step') self.fields.step() @@ -3732,20 +3765,26 @@ def run(self, *step_funcs, **kwargs): until_after_sources = kwargs.pop('until_after_sources', None) if self.fields is None: + print("init_sim") self.init_sim() + print("_evaluate_dft_objects") self._evaluate_dft_objects() + print("_check_material_frequencies") self._check_material_frequencies() if kwargs: raise ValueError("Unrecognized keyword arguments: {}".format(kwargs.keys())) if until_after_sources is not None: + print("_run_sources_until") self._run_sources_until(until_after_sources, step_funcs) elif until is not None: + print("_run_until") self._run_until(until, step_funcs) else: raise ValueError("Invalid run configuration") + print("run_sim done.") def print_times(self): """ diff --git a/python/tests/test_dump_load.py b/python/tests/test_dump_load.py new file mode 100644 index 000000000..d018218ae --- /dev/null +++ b/python/tests/test_dump_load.py @@ -0,0 +1,107 @@ +import itertools +import os +import re +import sys +import unittest +import warnings +import h5py +import numpy as np +import meep as mp + +try: + unicode +except NameError: + unicode = str + + +class TestLoadDump(unittest.TestCase): + + fname_base = re.sub(r'\.py$', '', os.path.split(sys.argv[0])[1]) + fname = fname_base + '-ez-000200.00.h5' + + def setUp(self): + print("Running {}".format(self._testMethodName)) + + @classmethod + def setUpClass(cls): + cls.temp_dir = mp.make_output_directory() + + @classmethod + def tearDownClass(cls): + # mp.delete_directory(cls.temp_dir) + pass + + def _load_dump_structure(self, chunk_file=False, chunk_sim=False, single_parallel_file=True): + from meep.materials import Al + resolution = 50 + cell = mp.Vector3(5, 5) + sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.2), center=mp.Vector3(), component=mp.Ez) + one_by_one = mp.Vector3(1, 1, mp.inf) + geometry = [mp.Block(material=Al, center=mp.Vector3(), size=one_by_one), + mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)] + pml_layers = [mp.PML(0.5)] + + symmetries = [mp.Mirror(mp.Y)] + + sim1 = mp.Simulation(resolution=resolution, + cell_size=cell, + boundary_layers=pml_layers, + geometry=geometry, + symmetries=symmetries, + sources=[sources]) + + sample_point = mp.Vector3(0.12, -0.29) + ref_field_points = [] + + def get_ref_field_point(sim): + p = sim.get_field_point(mp.Ez, sample_point) + ref_field_points.append(p.real) + + sim1.run(mp.at_every(5, get_ref_field_point), until=50) + + dump_dirname = os.path.join(self.temp_dir, 'test_load_dump') + sim1.dump(dump_dirname, structure=True, fields=True, single_parallel_file=single_parallel_file) + + dump_chunk_fname = None + chunk_layout = None + if chunk_file: + dump_chunk_fname = os.path.join(dump_dirname, 'chunk_layout.h5') + sim1.dump_chunk_layout(dump_chunk_fname) + chunk_layout = dump_chunk_fname + if chunk_sim: + chunk_layout = sim1 + + sim = mp.Simulation(resolution=resolution, + cell_size=cell, + boundary_layers=pml_layers, + sources=[sources], + symmetries=symmetries, + chunk_layout=chunk_layout) + sim.load(dump_dirname, structure=True, fields=True, single_parallel_file=single_parallel_file) + field_points = [] + + def get_field_point(sim): + p = sim.get_field_point(mp.Ez, sample_point) + field_points.append(p.real) + + sim.run(mp.at_every(5, get_field_point), until=50) + + for ref_pt, pt in zip(ref_field_points, field_points): + self.assertAlmostEqual(ref_pt, pt) + + def test_load_dump_structure(self): + self._load_dump_structure() + + @unittest.skipIf(not mp.with_mpi(), "MPI specific test") + def test_load_dump_structure_sharded(self): + self._load_dump_structure(single_parallel_file=False) + + def test_load_dump_chunk_layout_file(self): + self._load_dump_structure(chunk_file=True) + + def test_load_dump_chunk_layout_sim(self): + self._load_dump_structure(chunk_sim=True) + + +if __name__ == '__main__': + unittest.main() diff --git a/python/tests/test_simulation.py b/python/tests/test_simulation.py index f85afeb7f..24d183467 100644 --- a/python/tests/test_simulation.py +++ b/python/tests/test_simulation.py @@ -302,77 +302,6 @@ def test_in_box_volumes(self): sim.field_energy_in_box(tv) sim.field_energy_in_box(v) - def _load_dump_structure(self, chunk_file=False, chunk_sim=False, single_parallel_file=True): - from meep.materials import Al - resolution = 50 - cell = mp.Vector3(5, 5) - sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.2), center=mp.Vector3(), component=mp.Ez) - one_by_one = mp.Vector3(1, 1, mp.inf) - geometry = [mp.Block(material=Al, center=mp.Vector3(), size=one_by_one), - mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)] - pml_layers = [mp.PML(0.5)] - - symmetries = [mp.Mirror(mp.Y)] - - sim1 = mp.Simulation(resolution=resolution, - cell_size=cell, - boundary_layers=pml_layers, - geometry=geometry, - symmetries=symmetries, - sources=[sources]) - - sample_point = mp.Vector3(0.12, -0.29) - ref_field_points = [] - - def get_ref_field_point(sim): - p = sim.get_field_point(mp.Ez, sample_point) - ref_field_points.append(p.real) - - sim1.run(mp.at_every(5, get_ref_field_point), until=50) - - dump_dirname = os.path.join(self.temp_dir, 'test_load_dump') - sim1.dump(dump_dirname, structure=True, single_parallel_file=single_parallel_file) - - dump_chunk_fname = None - chunk_layout = None - if chunk_file: - dump_chunk_fname = os.path.join(dump_dirname, 'chunk_layout.h5') - sim1.dump_chunk_layout(dump_chunk_fname) - chunk_layout = dump_chunk_fname - if chunk_sim: - chunk_layout = sim1 - - sim = mp.Simulation(resolution=resolution, - cell_size=cell, - boundary_layers=pml_layers, - sources=[sources], - symmetries=symmetries, - chunk_layout=chunk_layout) - sim.load(dump_dirname, structure=True, single_parallel_file=single_parallel_file) - field_points = [] - - def get_field_point(sim): - p = sim.get_field_point(mp.Ez, sample_point) - field_points.append(p.real) - - sim.run(mp.at_every(5, get_field_point), until=50) - - for ref_pt, pt in zip(ref_field_points, field_points): - self.assertAlmostEqual(ref_pt, pt) - - def test_load_dump_structure(self): - self._load_dump_structure() - - @unittest.skipIf(not mp.with_mpi(), "MPI specific test") - def test_load_dump_structure_sharded(self): - self._load_dump_structure(single_parallel_file=False) - - def test_load_dump_chunk_layout_file(self): - self._load_dump_structure(chunk_file=True) - - def test_load_dump_chunk_layout_sim(self): - self._load_dump_structure(chunk_sim=True) - def test_get_array_output(self): sim = self.init_simple_simulation() sim.use_output_directory(self.temp_dir) diff --git a/src/Makefile.am b/src/Makefile.am index 469f8aec2..bf2111a57 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -12,7 +12,7 @@ bicgstab.hpp meepgeom.hpp material_data.hpp adjust_verbosity.hpp libmeep_la_SOURCES = array_slice.cpp anisotropic_averaging.cpp \ bands.cpp boundaries.cpp bicgstab.cpp casimir.cpp \ cw_fields.cpp dft.cpp dft_ldos.cpp energy_and_flux.cpp \ -fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp \ +fields.cpp fields_dump.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp \ initialize.cpp integrate.cpp integrate2.cpp material_data.cpp monitor.cpp mympi.cpp \ multilevel-atom.cpp near2far.cpp output_directory.cpp random.cpp \ sources.cpp step.cpp step_db.cpp stress.cpp structure.cpp structure_dump.cpp \ diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp new file mode 100644 index 000000000..a19f41b52 --- /dev/null +++ b/src/fields_dump.cpp @@ -0,0 +1,244 @@ +/* Copyright (C) 2005-2021 Massachusetts Institute of Technology +% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2, or (at your option) +% any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software Foundation, +% Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +*/ + +// Dump/load raw fields data to/from an HDF5 file. Only +// works if the number of processors/chunks is the same. + +#include +#include +#include +#include + +#include + +#include "meep.hpp" +#include "meep_internals.hpp" + +namespace meep { + +void fields::dump_fields_chunk_field(h5file *h5f, bool single_parallel_file, + const std::string &field_name, + FieldPtrGetter field_ptr_getter) { + /* + * make/save a num_chunks x NUM_FIELD_COMPONENTS x 2 array counting + * the number of entries in the 'field_name' array for each chunk. + * + * When 'single_parallel_file' is true, we are creating a single block of data + * for ALL chunks (that are merged using MPI). Otherwise, we are just + * making a copy of just the chunks that are ours. + */ + int my_num_chunks = 0; + for (int i = 0; i < num_chunks; i++) { + my_num_chunks += (single_parallel_file || chunks[i]->is_mine()); + } + size_t num_f_size = my_num_chunks * NUM_FIELD_COMPONENTS * 2; + std::vector num_f_(num_f_size); + size_t my_ntot = 0; + for (int i = 0, chunk_i = 0; i < num_chunks; i++) { + printf( + "dump: i = %d, chunk_i = %d, is_mine = %d, num_chunks = %d, my_num_chunks = " + "%d\n", + i, chunk_i, chunks[i]->is_mine(), num_chunks, my_num_chunks); + if (chunks[i]->is_mine()) { + size_t ntot = chunks[i]->gv.ntot(); + for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { + for (int d = 0; d < 2; ++d) { + realnum **f = field_ptr_getter(chunks[i], c, d); + if (*f) { + my_ntot += (num_f_[(chunk_i * NUM_FIELD_COMPONENTS + c) * 2 + d] = ntot); + printf("dump: c=%d d=%d ntot=%zu my_ntot=%zu\n", c, d, ntot, my_ntot); + } + } + } + } + chunk_i += (chunks[i]->is_mine() || single_parallel_file); + } + + std::vector num_f; + if (single_parallel_file) { + num_f.resize(num_f_size); + sum_to_master(num_f_.data(), num_f.data(), num_f_size); + } else { + num_f = std::move(num_f_); + } + + /* determine total dataset size and offset of this process's data */ + size_t my_start = 0; + size_t ntotal = my_ntot; + if (single_parallel_file) { + my_start = partial_sum_to_all(my_ntot) - my_ntot; + ntotal = sum_to_all(my_ntot); + } + + size_t dims[3] = {(size_t)my_num_chunks, NUM_FIELD_COMPONENTS, 2}; + size_t start[3] = {0, 0, 0}; + std::string num_f_name = std::string("num_") + field_name; + h5f->create_data(num_f_name.c_str(), 3, dims); + if (am_master() || !single_parallel_file) { + h5f->write_chunk(3, start, dims, num_f.data()); + } + + /* write the data */ + h5f->create_data(field_name.c_str(), 1, &ntotal, false /* append_data */, false /* single_precision */); + for (int i = 0; i < num_chunks; i++) { + if (chunks[i]->is_mine()) { + size_t ntot = chunks[i]->gv.ntot(); + for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { + for (int d = 0; d < 2; ++d) { + realnum **f = field_ptr_getter(chunks[i], c, d); + if (*f) { + h5f->write_chunk(1, &my_start, &ntot, *f); + my_start += ntot; + } + } + } + } + } +} + +void fields::dump(const char *filename, bool single_parallel_file) { + if (verbosity > 0) { + printf("creating fields output file \"%s\" (%d)...\n", filename, single_parallel_file); + } + + h5file file(filename, h5file::WRITE, single_parallel_file, !single_parallel_file); + + dump_fields_chunk_field( + &file, single_parallel_file, "f", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f[c][d]); }); + dump_fields_chunk_field( + &file, single_parallel_file, "f_u", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_u[c][d]); }); + dump_fields_chunk_field( + &file, single_parallel_file, "f_w", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_w[c][d]); }); + dump_fields_chunk_field( + &file, single_parallel_file, "f_cond", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); +} + +void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, + const std::string &field_name, + FieldPtrGetter field_ptr_getter) { + int my_num_chunks = 0; + for (int i = 0; i < num_chunks; i++) { + my_num_chunks += (single_parallel_file || chunks[i]->is_mine()); + } + size_t num_f_size = my_num_chunks * NUM_FIELD_COMPONENTS * 2; + std::vector num_f(num_f_size); + + int rank; + size_t dims[3], _dims[3] = {(size_t)my_num_chunks, NUM_FIELD_COMPONENTS, 2}; + size_t start[3] = {0, 0, 0}; + + std::string num_f_name = std::string("num_") + field_name; + h5f->read_size(num_f_name.c_str(), &rank, dims, 3); + if (rank != 3 || _dims[0] != dims[0] || _dims[1] != dims[1] || _dims[2] != dims[2]) + meep::abort("chunk mismatch in fields::load"); + if (am_master() || !single_parallel_file) h5f->read_chunk(3, start, dims, num_f.data()); + + if (single_parallel_file) { + h5f->prevent_deadlock(); + broadcast(0, num_f.data(), dims[0] * dims[1] * dims[2]); + } + + /* allocate data as needed and check sizes */ + size_t my_ntot = 0; + for (int i = 0, chunk_i = 0; i < num_chunks; i++) { + printf( + "load: i = %d, chunk_i = %d, is_mine = %d, num_chunks = %d, my_num_chunks = " + "%d\n", + i, chunk_i, chunks[i]->is_mine(), num_chunks, my_num_chunks); + if (chunks[i]->is_mine()) { + size_t ntot = chunks[i]->gv.ntot(); + for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { + for (int d = 0; d < 2; ++d) { + size_t n = num_f[(chunk_i * NUM_FIELD_COMPONENTS + c) * 2 + d]; + realnum **f = field_ptr_getter(chunks[i], c, d); + if (n == 0) { + delete[] * f; + *f = NULL; + } else { + if (n != ntot) + meep::abort("grid size mismatch %zd vs %zd in fields::load", n, ntot); + *f = new realnum[ntot]; + my_ntot += ntot; + printf("load: c=%d d=%d ntot=%zu my_ntot=%zu\n", c, d, ntot, my_ntot); + } + } + } + } + chunk_i += (chunks[i]->is_mine() || single_parallel_file); + } + + printf("000\n"); + /* determine total dataset size and offset of this process's data */ + size_t my_start = 0; + size_t ntotal = my_ntot; + if (single_parallel_file) { + my_start = partial_sum_to_all(my_ntot) - my_ntot; + ntotal = sum_to_all(my_ntot); + } + + printf("001\n"); + /* read the data */ + h5f->read_size(field_name.c_str(), &rank, dims, 1); + if (rank != 1 || dims[0] != ntotal) { + meep::abort( + "inconsistent data size for '%s' in fields::load (rank, dims[0]): " + "(%d, %zu) != (1, %zu)", + field_name.c_str(), rank, dims[0], ntotal); + } + printf("002\n"); + for (int i = 0; i < num_chunks; i++) { + if (chunks[i]->is_mine()) { + size_t ntot = chunks[i]->gv.ntot(); + for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { + for (int d = 0; d < 2; ++d) { + realnum **f = field_ptr_getter(chunks[i], c, d); + if (*f) { + h5f->read_chunk(1, &my_start, &ntot, *f); + my_start += ntot; + } + } + } + } + } + printf("003\n"); +} + +void fields::load(const char *filename, bool single_parallel_file) { + if (verbosity > 0) + printf("reading fields from file \"%s\" (%d)...\n", filename, single_parallel_file); + + h5file file(filename, h5file::READONLY, single_parallel_file, !single_parallel_file); + + load_fields_chunk_field( + &file, single_parallel_file, "f", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f[c][d]); }); + load_fields_chunk_field( + &file, single_parallel_file, "f_u", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_u[c][d]); }); + load_fields_chunk_field( + &file, single_parallel_file, "f_w", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_w[c][d]); }); + load_fields_chunk_field( + &file, single_parallel_file, "f_cond", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); +} + +} // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index 7f336280f..b09a04e1b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1734,6 +1734,15 @@ class fields { volume total_volume(void) const; + // fields_dump.cpp + // Dump fields to specified file. If 'single_parallel_file' + // is 'true' (the default) - then all processes write to the same/single file + // file after computing their respective offsets into this file. When set to + // 'false', each process writes data for the chunks it owns to a separate + // (process unique) file. + void dump(const char *filename, bool single_parallel_file=true); + void load(const char *filename, bool single_parallel_file=true); + // h5fields.cpp: // low-level function: void output_hdf5(h5file *file, const char *dataname, int num_fields, const component *components, @@ -2206,6 +2215,14 @@ class fields { std::complex A(const vec &), std::complex amp, component c0, direction d, int has_tm, int has_te); + // fields_dump.cpp + // Helper methods for dumping field chunks. + using FieldPtrGetter = std::function; + void dump_fields_chunk_field(h5file *h5f, bool single_parallel_file, + const std::string& field_name, FieldPtrGetter field_ptr_getter); + void load_fields_chunk_field(h5file *h5f, bool single_parallel_file, + const std::string& field_name, FieldPtrGetter field_ptr_getter); + public: // monitor.cpp std::complex get_field(component c, const ivec &iloc, bool parallel = true) const; diff --git a/tests/Makefile.am b/tests/Makefile.am index 53eae02ae..1ce61b95b 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,8 +1,8 @@ SRC = aniso_disp.cpp bench.cpp bragg_transmission.cpp \ -convergence_cyl_waveguide.cpp cylindrical.cpp flux.cpp harmonics.cpp \ -integrate.cpp known_results.cpp near2far.cpp one_dimensional.cpp \ -physical.cpp stress_tensor.cpp symmetry.cpp three_d.cpp \ -two_dimensional.cpp 2D_convergence.cpp h5test.cpp pml.cpp +convergence_cyl_waveguide.cpp cylindrical.cpp dump_load.cpp flux.cpp \ +harmonics.cpp integrate.cpp known_results.cpp near2far.cpp \ +one_dimensional.cpp physical.cpp stress_tensor.cpp symmetry.cpp \ +three_d.cpp two_dimensional.cpp 2D_convergence.cpp h5test.cpp pml.cpp EXTRA_DIST = $(SRC) @@ -16,7 +16,7 @@ AM_CPPFLAGS = -I$(top_srcdir)/src -DDATADIR=\"$(srcdir)/\" .SUFFIXES = .dac .done -check_PROGRAMS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml pw-source-ll ring-ll cyl-ellipsoid-ll absorber-1d-ll array-slice-ll user-defined-material dft-fields gdsII-3d bend-flux-ll array-metadata +check_PROGRAMS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical dump_load flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml pw-source-ll ring-ll cyl-ellipsoid-ll absorber-1d-ll array-slice-ll user-defined-material dft-fields gdsII-3d bend-flux-ll array-metadata array_metadata_SOURCES = array-metadata.cpp array_metadata_LDADD = $(MEEPLIBS) @@ -36,6 +36,9 @@ convergence_cyl_waveguide_LDADD = $(MEEPLIBS) cylindrical_SOURCES = cylindrical.cpp cylindrical_LDADD = $(MEEPLIBS) +dump_load_SOURCES = dump_load.cpp +dump_load_LDADD = $(MEEPLIBS) + flux_SOURCES = flux.cpp flux_LDADD = $(MEEPLIBS) @@ -107,7 +110,7 @@ user_defined_material_LDADD = $(MEEPLIBS) dist_noinst_DATA = cyl-ellipsoid-eps-ref.h5 array-slice-ll-ref.h5 gdsII-3d.gds -TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml +TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical dump_load flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml if WITH_MPI LOG_COMPILER = $(RUNCODE) diff --git a/tests/dump_load.cpp b/tests/dump_load.cpp new file mode 100644 index 000000000..3adaa72c4 --- /dev/null +++ b/tests/dump_load.cpp @@ -0,0 +1,369 @@ +/* Copyright (C) 2005-2021 Massachusetts Institute of Technology +% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2, or (at your option) +% any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software Foundation, +% Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +*/ + +#include +#include +#include + +#include +using namespace meep; +using std::complex; + +double one(const vec &) { return 1.0; } +double targets(const vec &pt) { + const double r = sqrt(pt.x() * pt.x() + pt.y() * pt.y()); + double dr = r; + while (dr > 1) dr -= 1; + if (dr > 0.7001) return 12.0; + return 1.0; +} + +#if MEEP_SINGLE +static const double tol = 1e-3, thresh = 1e-3; +#else +static const double tol = 1e-9, thresh = 1e-15; +#endif + +int compare(double a, double b, const char *n) { + if (fabs(a - b) > fabs(b) * tol && fabs(b) > thresh) { + master_printf("%s differs by\t%g out of\t%g\n", n, a - b, b); + master_printf("This gives a fractional error of %g\n", + fabs(a - b) / fabs(b)); + return 0; + } else { + return 1; + } +} + +int compare_point(fields &f1, fields &f2, const vec &p) { + monitor_point m1, m_test; + f1.get_point(&m_test, p); + f2.get_point(&m1, p); + for (int i = 0; i < 10; i++) { + component c = (component)i; + if (f1.gv.has_field(c)) { + complex v1 = m_test.get_component(c), v2 = m1.get_component(c); + if (abs(v1 - v2) > tol * abs(v2) && abs(v2) > thresh) { + master_printf("%s differs: %g %g out of %g %g\n", component_name(c), + real(v2 - v1), imag(v2 - v1), real(v2), imag(v2)); + master_printf("This comes out to a fractional error of %g\n", + abs(v1 - v2) / abs(v2)); + master_printf("Right now I'm looking at %g %g %g, time %g\n", p.x(), + p.y(), p.z(), f1.time()); + return 0; + } + } + } + return 1; +} + +int approx_point(fields &f1, fields &f2, const vec &p) { + monitor_point m1, m_test; + f1.get_point(&m_test, p); + f2.get_point(&m1, p); + for (int i = 0; i < 10; i++) { + component c = (component)i; + if (f1.gv.has_field(c)) { + complex v1 = m_test.get_component(c), v2 = m1.get_component(c); + if (abs(v1 - v2) > tol * abs(v2) && abs(v2) > thresh) { + master_printf("%s differs: %g %g out of %g %g\n", component_name(c), + real(v2 - v1), imag(v2 - v1), real(v2), imag(v2)); + master_printf("This comes out to a fractional error of %g\n", + abs(v1 - v2) / abs(v2)); + master_printf("Right now I'm looking at %g %g %g, time %g\n", p.x(), + p.y(), p.z(), f1.time()); + return 0; + } + } + } + return 1; +} + +std::string structure_dump(structure *s, const std::string &filename_prefix, + const std::string &output_name) { + std::string filename = filename_prefix + "-structure-" + output_name; + s->dump(filename.c_str()); + master_printf("Dumping structure: %s\n", filename.c_str()); + return filename; +} + +void structure_load(structure *s, const std::string &filename) { + master_printf("Loading structure: %s\n", filename.c_str()); + s->load(filename.c_str()); +} + +std::string fields_dump(fields *f, const std::string &filename_prefix, + const std::string &output_name) { + std::string filename = filename_prefix + "-fields-" + output_name; + f->dump(filename.c_str()); + master_printf("Dumping fields: %s\n", filename.c_str()); + return filename; +} + +void fields_load(fields *f, const std::string &filename) { + master_printf("Loading fields: %s\n", filename.c_str()); + f->load(filename.c_str()); +} + +int test_metal(double eps(const vec &), int splitting, const char *tmpdir) { + double a = 10.0; + double ttot = 17.0; + + grid_volume gv = vol3d(1.5, 0.5, 1.0, a); + structure s1(gv, eps); + structure s(gv, eps, no_pml(), identity(), splitting); + + std::string filename_prefix = + std::string(tmpdir) + "/test_metal_" + std::to_string(splitting); + std::string structure_filename_original = + structure_dump(&s, filename_prefix, "original"); + + master_printf("Metal test using %d chunks...\n", splitting); + fields f(&s); + f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.299, 0.401), 1.0); + + std::string fields_filename_original = + fields_dump(&f, filename_prefix, "original"); + + fields f1(&s1); + f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.299, 0.401), 1.0); + + double field_energy_check_time = 8.0; + while (f.time() < ttot) { + f.step(); + f1.step(); + if (!compare_point(f, f1, vec(0.5, 0.5, 0.01))) return 0; + if (!compare_point(f, f1, vec(0.46, 0.33, 0.33))) return 0; + if (!compare_point(f, f1, vec(1.301, 0.301, 0.399))) return 0; + if (f.time() >= field_energy_check_time) { + if (!compare(f.field_energy(), f1.field_energy(), " total energy")) + return 0; + if (!compare(f.electric_energy_in_box(gv.surroundings()), + f1.electric_energy_in_box(gv.surroundings()), + "electric energy")) + return 0; + if (!compare(f.magnetic_energy_in_box(gv.surroundings()), + f1.magnetic_energy_in_box(gv.surroundings()), + "magnetic energy")) + return 0; + field_energy_check_time += 5.0; + } + } + + std::string structure_filename_after_sim = + structure_dump(&s, filename_prefix, "after-sim"); + + std::string fields_filename_after_sim = + fields_dump(&f, filename_prefix, "after_sim"); + + structure s_load(gv, eps, no_pml(), identity(), splitting); + structure_load(&s_load, structure_filename_after_sim); + + std::string structure_filename_dump_loaded = + structure_dump(&s, filename_prefix, "dump-loaded"); + + fields f_load(&s_load); + f_load.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.299, 0.401), + 1.0); + fields_load(&f_load, fields_filename_after_sim); + + std::string fields_filename_dump_loaded = + fields_dump(&f_load, filename_prefix, "dump-loaded"); + return 1; +} + +int test_periodic(double eps(const vec &), int splitting, const char *tmpdir) { + double a = 10.0; + double ttot = 17.0; + + grid_volume gv = vol3d(1.5, 0.5, 1.0, a); + structure s1(gv, eps); + structure s(gv, eps, no_pml(), identity(), splitting); + + std::string filename_prefix = + std::string(tmpdir) + "/test_periodic_" + std::to_string(splitting); + std::string structure_filename_original = + structure_dump(&s, filename_prefix, "original"); + + master_printf("Periodic test using %d chunks...\n", splitting); + fields f(&s); + f.use_bloch(vec(0.1, 0.7, 0.3)); + f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.25, 0.5), 1.0); + std::string fields_filename_original = + fields_dump(&f, filename_prefix, "original"); + + fields f1(&s1); + f1.use_bloch(vec(0.1, 0.7, 0.3)); + f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.25, 0.5), 1.0); + + double field_energy_check_time = 8.0; + while (f.time() < ttot) { + f.step(); + f1.step(); + if (!compare_point(f, f1, vec(0.5, 0.01, 0.5))) return 0; + if (!compare_point(f, f1, vec(0.46, 0.33, 0.2))) return 0; + if (!compare_point(f, f1, vec(1.0, 0.25, 0.301))) return 0; + if (f.time() >= field_energy_check_time) { + if (!compare(f.field_energy(), f1.field_energy(), " total energy")) + return 0; + if (!compare(f.electric_energy_in_box(gv.surroundings()), + f1.electric_energy_in_box(gv.surroundings()), + "electric energy")) + return 0; + if (!compare(f.magnetic_energy_in_box(gv.surroundings()), + f1.magnetic_energy_in_box(gv.surroundings()), + "magnetic energy")) + return 0; + field_energy_check_time += 5.0; + } + } + + std::string structure_filename_after_sim = + structure_dump(&s, filename_prefix, "after-sim"); + + std::string fields_filename_after_sim = + fields_dump(&f, filename_prefix, "after_sim"); + + structure s_load(gv, eps, no_pml(), identity(), splitting); + structure_load(&s_load, structure_filename_after_sim); + + std::string structure_filename_dump_loaded = + structure_dump(&s, filename_prefix, "dump-loaded"); + + fields f_load(&s_load); + f_load.use_bloch(vec(0.1, 0.7, 0.3)); + f_load.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.25, 0.5), 1.0); + fields_load(&f_load, fields_filename_after_sim); + + std::string fields_filename_dump_loaded = + fields_dump(&f_load, filename_prefix, "dump-loaded"); + return 1; +} + +int test_pml(double eps(const vec &), const char *tmpdir) { + master_printf("Testing pml quality...\n"); + + double a = 10.0; + grid_volume gv = vol3d(1.5, 1.0, 1.2, a); + structure s(gv, eps, pml(0.401)); + + std::string filename_prefix = std::string(tmpdir) + "/test_pml"; + std::string structure_filename_original = + structure_dump(&s, filename_prefix, "original"); + + fields f(&s); + f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.751, 0.5, 0.601), 1.0); + const double deltaT = 10.0; + const double ttot = 3.1 * deltaT; + double field_energy_check_time = deltaT; + + std::string fields_filename_original = + fields_dump(&f, filename_prefix, "original"); + + while (f.time() < f.last_source_time()) f.step(); + + double last_energy = f.field_energy(); + while (f.time() < ttot) { + f.step(); + if (f.time() >= field_energy_check_time) { + const double new_energy = f.field_energy(); + master_printf("Got newE/oldE of %g\n", new_energy / last_energy); + if (new_energy > last_energy * 4e-3) { + master_printf("Energy decaying too slowly: from %g to %g (%g)\n", + last_energy, new_energy, new_energy / last_energy); + return 0; + } + field_energy_check_time += deltaT; + } + } + + std::string structure_filename_after_sim = + structure_dump(&s, filename_prefix, "after-sim"); + + std::string fields_filename_after_sim = + fields_dump(&f, filename_prefix, "after_sim"); + + structure s_load(gv, eps, pml(0.401)); + structure_load(&s_load, structure_filename_after_sim); + + std::string structure_filename_dump_loaded = + structure_dump(&s, filename_prefix, "dump-loaded"); + + fields f_load(&s_load); + fields_load(&f_load, fields_filename_after_sim); + + std::string fields_filename_dump_loaded = + fields_dump(&f_load, filename_prefix, "dump-loaded"); + return 1; +} + +int test_pml_splitting(double eps(const vec &), int splitting, + const char *tmpdir) { + double a = 10.0; + + grid_volume gv = vol3d(1.5, 1.0, 1.2, a); + structure s1(gv, eps, pml(0.3)); + structure s(gv, eps, pml(0.3), identity(), splitting); + + master_printf("Testing pml while splitting into %d chunks...\n", splitting); + fields f(&s); + f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.099, 0.499, 0.501), 1.0); + fields f1(&s1); + f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.099, 0.499, 0.501), 1.0); + const double ttot = 31.0; + + double next_energy_time = 10.0; + while (f.time() < ttot) { + f.step(); + f1.step(); + if (!approx_point(f, f1, vec(0.5, 0.01, 1.0))) return 0; + if (!approx_point(f, f1, vec(0.46, 0.33, 0.33))) return 0; + if (!approx_point(f, f1, vec(1.0, 1.0, 0.33))) return 0; + if (!approx_point(f, f1, vec(1.3, 0.3, 0.15))) return 0; + if (f.time() > next_energy_time) { + if (!compare(f.field_energy(), f1.field_energy(), " total energy")) + return 0; + next_energy_time += 10.0; + } + } + return 1; +} + +int main(int argc, char **argv) { + initialize mpi(argc, argv); + verbosity = 0; + + std::unique_ptr temp_dir(make_output_directory()); + master_printf("Testing 3D dump/load: temp_dir = %s...\n", temp_dir.get()); + + if (!test_pml(one, temp_dir.get())) abort("error in test_pml vacuum\n"); + + for (int s = 2; s < 7; s++) + if (!test_periodic(targets, s, temp_dir.get())) + abort("error in test_periodic targets\n"); + + for (int s = 2; s < 8; s++) + if (!test_metal(one, s, temp_dir.get())) + abort("error in test_metal vacuum\n"); + + for (int s = 2; s < 4; s++) + if (!test_pml_splitting(one, s, temp_dir.get())) + abort("error in test_pml_splitting vacuum\n"); + + // delete_directory(temp_dir.get()); + return 0; +} From 16eb66112779996ab7d65108d2ea5ea048f643df Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Tue, 24 Aug 2021 01:13:15 +0000 Subject: [PATCH 06/17] Save dft chunks --- python/simulation.py | 2 +- python/tests/test_dump_load.py | 5 +++++ src/dft.cpp | 35 ++++++++++++++++++++++++---------- src/fields_dump.cpp | 18 +++++++++++++++++ src/meep.hpp | 8 ++++---- 5 files changed, 53 insertions(+), 15 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 1acd21253..37ec45b46 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -2269,7 +2269,7 @@ def stop_cond(sim): while not any([x(self) for x in cond]): for func in step_funcs: - print('_run_until: _eval_step_func: ' + str(func)) + print('_run_until: _eval_step_func: ' + str(func) + ', time=' + str(self.round_time())) _eval_step_func(self, func, 'step') self.fields.step() diff --git a/python/tests/test_dump_load.py b/python/tests/test_dump_load.py index d018218ae..74c5367bf 100644 --- a/python/tests/test_dump_load.py +++ b/python/tests/test_dump_load.py @@ -7,6 +7,7 @@ import h5py import numpy as np import meep as mp +import pdb try: unicode @@ -90,15 +91,19 @@ def get_field_point(sim): self.assertAlmostEqual(ref_pt, pt) def test_load_dump_structure(self): + # pdb.set_trace() self._load_dump_structure() @unittest.skipIf(not mp.with_mpi(), "MPI specific test") + @unittest.skip("foo") def test_load_dump_structure_sharded(self): self._load_dump_structure(single_parallel_file=False) + @unittest.skip("foo") def test_load_dump_chunk_layout_file(self): self._load_dump_structure(chunk_file=True) + @unittest.skip("foo") def test_load_dump_chunk_layout_sim(self): self._load_dump_structure(chunk_sim=True) diff --git a/src/dft.cpp b/src/dft.cpp index cf5177caa..4ab6a70c0 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -284,18 +284,33 @@ void dft_chunk::operator-=(const dft_chunk &chunk) { } } -size_t dft_chunks_Ntotal(dft_chunk *dft_chunks, size_t *my_start) { - size_t n = 0; +size_t my_dft_chunks_Ntotal(dft_chunk *dft_chunks, size_t *my_start) { + // When writing to a sharded file, we write out only the chunks we own. +size_t n = 0; for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_dft) n += cur->N * cur->omega.size() * 2; + + *my_start = 0; + return n; +} + +size_t dft_chunks_Ntotal(dft_chunk *dft_chunks, size_t *my_start) { + // If writing to a single parallel file, we are compute our chunks offset + // into the single-parallel-file that has all the data. + size_t n = my_dft_chunks_Ntotal(dft_chunks, my_start); *my_start = partial_sum_to_all(n) - n; // sum(n) for processes before this return sum_to_all(n); } +size_t dft_chunks_Ntotal(dft_chunk *dft_chunks, size_t *my_start, bool single_parallel_file) { + return single_parallel_file ? dft_chunks_Ntotal(dft_chunks, my_start) : + my_dft_chunks_Ntotal(dft_chunks, my_start); +} + // Note: the file must have been created in parallel mode, typically via fields::open_h5file. -void save_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix) { +void save_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix, bool single_parallel_file) { size_t istart; - size_t n = dft_chunks_Ntotal(dft_chunks, &istart); + size_t n = dft_chunks_Ntotal(dft_chunks, &istart, single_parallel_file); char dataname[1024]; snprintf(dataname, 1024, @@ -312,13 +327,13 @@ void save_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const file->done_writing_chunks(); } -void save_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix) { - save_dft_hdf5(dft_chunks, component_name(c), file, dprefix); +void save_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix, bool single_parallel_file) { + save_dft_hdf5(dft_chunks, component_name(c), file, dprefix, single_parallel_file); } -void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix) { +void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix, bool single_parallel_file) { size_t istart; - size_t n = dft_chunks_Ntotal(dft_chunks, &istart); + size_t n = dft_chunks_Ntotal(dft_chunks, &istart, single_parallel_file); char dataname[1024]; snprintf(dataname, 1024, @@ -339,8 +354,8 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const } } -void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix) { - load_dft_hdf5(dft_chunks, component_name(c), file, dprefix); +void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix, bool single_parallel_file) { + load_dft_hdf5(dft_chunks, component_name(c), file, dprefix, single_parallel_file); } dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index a19f41b52..811523248 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -129,6 +129,15 @@ void fields::dump(const char *filename, bool single_parallel_file) { dump_fields_chunk_field( &file, single_parallel_file, "f_cond", [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); + + // Dump DFT chunks. + for (int i = 0; i < num_chunks; i++) { + if (chunks[i]->is_mine()) { + char dataname[1024]; + snprintf(dataname, 1024, "chunk%02d", i); + save_dft_hdf5(chunks[i]->dft_chunks, dataname, &file, 0, single_parallel_file); + } + } } void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, @@ -239,6 +248,15 @@ void fields::load(const char *filename, bool single_parallel_file) { load_fields_chunk_field( &file, single_parallel_file, "f_cond", [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); + + // Load DFT chunks. + for (int i = 0; i < num_chunks; i++) { + if (chunks[i]->is_mine()) { + char dataname[1024]; + snprintf(dataname, 1024, "chunk%02d", i); + load_dft_hdf5(chunks[i]->dft_chunks, dataname, &file, 0, single_parallel_file); + } + } } } // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index 0ac0ad556..d522781f7 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1193,10 +1193,10 @@ class dft_chunk { int decimation_factor; }; -void save_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix = 0); -void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix = 0); -void save_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix = 0); -void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix = 0); +void save_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix = 0, bool single_parallel_file=true); +void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char *dprefix = 0, bool single_parallel_file=true); +void save_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix = 0, bool single_parallel_file=true); +void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const char *dprefix = 0, bool single_parallel_file=true); // dft.cpp (normally created with fields::add_dft_flux) class dft_flux { From 6cab1250fd5365a7385dcfe5ef4b0b125dc47851 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Tue, 24 Aug 2021 01:20:34 +0000 Subject: [PATCH 07/17] Remove debug log stmts --- python/simulation.py | 3 --- src/fields_dump.cpp | 14 -------------- 2 files changed, 17 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 37ec45b46..80ab7d34f 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1983,7 +1983,6 @@ def init_sim(self): ) if self.load_fields_file: - print("Loading fields from file: %s (%s)" % (self.load_fields_file, str(self.load_single_parallel_file))) self.load_fields( self.load_fields_file, self.load_single_parallel_file) print("LOADED fields from file: %s (%s)" % (self.load_fields_file, str(self.load_single_parallel_file))) @@ -2007,7 +2006,6 @@ def init_sim(self): hook() self._is_initialized = True - print("Sim INITIALIZED.") def using_real_fields(self): cond1 = self.is_cylindrical and self.m != 0 @@ -2251,7 +2249,6 @@ def _run_until(self, cond, step_funcs): if isinstance(cond[i], numbers.Number): stop_time = cond[i] t0 = self.round_time() - print('_run_until: stop_time = %s, t0 = %s' % (str(stop_time), str(t0))) def stop_cond(sim): return sim.round_time() >= t0 + stop_time diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index 811523248..79dd55380 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -49,10 +49,6 @@ void fields::dump_fields_chunk_field(h5file *h5f, bool single_parallel_file, std::vector num_f_(num_f_size); size_t my_ntot = 0; for (int i = 0, chunk_i = 0; i < num_chunks; i++) { - printf( - "dump: i = %d, chunk_i = %d, is_mine = %d, num_chunks = %d, my_num_chunks = " - "%d\n", - i, chunk_i, chunks[i]->is_mine(), num_chunks, my_num_chunks); if (chunks[i]->is_mine()) { size_t ntot = chunks[i]->gv.ntot(); for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { @@ -60,7 +56,6 @@ void fields::dump_fields_chunk_field(h5file *h5f, bool single_parallel_file, realnum **f = field_ptr_getter(chunks[i], c, d); if (*f) { my_ntot += (num_f_[(chunk_i * NUM_FIELD_COMPONENTS + c) * 2 + d] = ntot); - printf("dump: c=%d d=%d ntot=%zu my_ntot=%zu\n", c, d, ntot, my_ntot); } } } @@ -168,10 +163,6 @@ void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, /* allocate data as needed and check sizes */ size_t my_ntot = 0; for (int i = 0, chunk_i = 0; i < num_chunks; i++) { - printf( - "load: i = %d, chunk_i = %d, is_mine = %d, num_chunks = %d, my_num_chunks = " - "%d\n", - i, chunk_i, chunks[i]->is_mine(), num_chunks, my_num_chunks); if (chunks[i]->is_mine()) { size_t ntot = chunks[i]->gv.ntot(); for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { @@ -186,7 +177,6 @@ void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, meep::abort("grid size mismatch %zd vs %zd in fields::load", n, ntot); *f = new realnum[ntot]; my_ntot += ntot; - printf("load: c=%d d=%d ntot=%zu my_ntot=%zu\n", c, d, ntot, my_ntot); } } } @@ -194,7 +184,6 @@ void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, chunk_i += (chunks[i]->is_mine() || single_parallel_file); } - printf("000\n"); /* determine total dataset size and offset of this process's data */ size_t my_start = 0; size_t ntotal = my_ntot; @@ -203,7 +192,6 @@ void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, ntotal = sum_to_all(my_ntot); } - printf("001\n"); /* read the data */ h5f->read_size(field_name.c_str(), &rank, dims, 1); if (rank != 1 || dims[0] != ntotal) { @@ -212,7 +200,6 @@ void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, "(%d, %zu) != (1, %zu)", field_name.c_str(), rank, dims[0], ntotal); } - printf("002\n"); for (int i = 0; i < num_chunks; i++) { if (chunks[i]->is_mine()) { size_t ntot = chunks[i]->gv.ntot(); @@ -227,7 +214,6 @@ void fields::load_fields_chunk_field(h5file *h5f, bool single_parallel_file, } } } - printf("003\n"); } void fields::load(const char *filename, bool single_parallel_file) { From c28eb2c93984ac45d40c27a2be517696c5ed68a6 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Tue, 24 Aug 2021 21:41:18 +0000 Subject: [PATCH 08/17] Add saving of time value and also reorg tests. --- python/simulation.py | 34 ++++++++----- python/tests/test_dump_load.py | 88 ++++++++++++++++++++++++++++++---- src/fields_dump.cpp | 17 +++++++ 3 files changed, 120 insertions(+), 19 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 80ab7d34f..262987925 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1866,6 +1866,7 @@ def dump_structure(self, fname, single_parallel_file=True): if self.structure is None: raise ValueError("Structure must be initialized before calling dump_structure") self.structure.dump(fname, single_parallel_file) + print("Dumped structure to file: %s (%s)" % (fname, str(single_parallel_file))) def load_structure(self, fname, single_parallel_file=True): """ @@ -1874,6 +1875,7 @@ def load_structure(self, fname, single_parallel_file=True): if self.structure is None: raise ValueError("Structure must be initialized before loading structure from file '%s'" % fname) self.structure.load(fname, single_parallel_file) + print("Loaded structure from file: %s (%s)" % (fname, str(single_parallel_file))) def dump_fields(self, fname, single_parallel_file=True): """ @@ -1882,6 +1884,7 @@ def dump_fields(self, fname, single_parallel_file=True): if self.fields is None: raise ValueError("Fields must be initialized before calling dump_fields") self.fields.dump(fname, single_parallel_file) + print("Dumped fields to file: %s (%s)" % (fname, str(single_parallel_file))) def load_fields(self, fname, single_parallel_file=True): """ @@ -1890,6 +1893,7 @@ def load_fields(self, fname, single_parallel_file=True): if self.fields is None: raise ValueError("Fields must be initialized before loading fields from file '%s'" % fname) self.fields.load(fname, single_parallel_file) + print("Loaded fields from file: %s (%s)" % (fname, str(single_parallel_file))) def dump_chunk_layout(self, fname): """ @@ -1923,22 +1927,22 @@ def _get_load_dump_dirname(self, dirname, single_parallel_file): dump_dirname = os.path.join(dirname, 'rank%02d' % mp.my_rank()) return dump_dirname - def dump(self, dirname, structure=True, fields=True, single_parallel_file=True): + def dump(self, dirname, dump_structure=True, dump_fields=True, single_parallel_file=True): """ Dumps simulation state. """ dump_dirname = self._get_load_dump_dirname(dirname, single_parallel_file) os.makedirs(dump_dirname, exist_ok=True) - if structure: + if dump_structure: structure_dump_filename = os.path.join(dump_dirname, 'structure.h5') self.dump_structure(structure_dump_filename, single_parallel_file) - if fields: + if dump_fields: fields_dump_filename = os.path.join(dump_dirname, 'fields.h5') self.dump_fields(fields_dump_filename, single_parallel_file) - def load(self, dirname, structure=True, fields=True, single_parallel_file=True): + def load(self, dirname, load_structure=True, load_fields=True, single_parallel_file=True): """ Loads simulation state. @@ -1948,11 +1952,21 @@ def load(self, dirname, structure=True, fields=True, single_parallel_file=True): dump_dirname = self._get_load_dump_dirname(dirname, single_parallel_file) self.load_single_parallel_file = single_parallel_file - if structure: - self.load_structure_file = os.path.join(dump_dirname, 'structure.h5') - - if fields: - self.load_fields_file = os.path.join(dump_dirname, 'fields.h5') + if load_structure: + load_structure_file = os.path.join(dump_dirname, 'structure.h5') + # If structure is already initialized, load it straight away. + # Otherwise, do a delayed load. + if self.structure: + self.load_structure(load_structure_file, self.load_single_parallel_file) + else: + self.load_structure_file = load_structure_file + + if load_fields: + load_fields_file = os.path.join(dump_dirname, 'fields.h5') + if self.fields: + self.load_fields(load_fields_file, self.load_single_parallel_file) + else: + self.load_fields_file = load_fields_file def init_sim(self): if self._is_initialized: @@ -1985,7 +1999,6 @@ def init_sim(self): if self.load_fields_file: self.load_fields( self.load_fields_file, self.load_single_parallel_file) - print("LOADED fields from file: %s (%s)" % (self.load_fields_file, str(self.load_single_parallel_file))) if self.force_all_components and self.dimensions != 1: self.fields.require_component(mp.Ez) @@ -2266,7 +2279,6 @@ def stop_cond(sim): while not any([x(self) for x in cond]): for func in step_funcs: - print('_run_until: _eval_step_func: ' + str(func) + ', time=' + str(self.round_time())) _eval_step_func(self, func, 'step') self.fields.step() diff --git a/python/tests/test_dump_load.py b/python/tests/test_dump_load.py index 74c5367bf..ac6460ccb 100644 --- a/python/tests/test_dump_load.py +++ b/python/tests/test_dump_load.py @@ -7,7 +7,6 @@ import h5py import numpy as np import meep as mp -import pdb try: unicode @@ -26,12 +25,14 @@ def setUp(self): @classmethod def setUpClass(cls): cls.temp_dir = mp.make_output_directory() + print("Saving temp files to dir: {}".format(cls.temp_dir)) @classmethod def tearDownClass(cls): # mp.delete_directory(cls.temp_dir) pass + # Tests various combinations of dumping/loading structure & chunk layout. def _load_dump_structure(self, chunk_file=False, chunk_sim=False, single_parallel_file=True): from meep.materials import Al resolution = 50 @@ -60,8 +61,8 @@ def get_ref_field_point(sim): sim1.run(mp.at_every(5, get_ref_field_point), until=50) - dump_dirname = os.path.join(self.temp_dir, 'test_load_dump') - sim1.dump(dump_dirname, structure=True, fields=True, single_parallel_file=single_parallel_file) + dump_dirname = os.path.join(self.temp_dir, 'test_load_dump_structure') + sim1.dump(dump_dirname, dump_structure=True, dump_fields=False, single_parallel_file=single_parallel_file) dump_chunk_fname = None chunk_layout = None @@ -78,7 +79,7 @@ def get_ref_field_point(sim): sources=[sources], symmetries=symmetries, chunk_layout=chunk_layout) - sim.load(dump_dirname, structure=True, fields=True, single_parallel_file=single_parallel_file) + sim.load(dump_dirname, load_structure=True, load_fields=False, single_parallel_file=single_parallel_file) field_points = [] def get_field_point(sim): @@ -90,23 +91,94 @@ def get_field_point(sim): for ref_pt, pt in zip(ref_field_points, field_points): self.assertAlmostEqual(ref_pt, pt) + @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_structure(self): - # pdb.set_trace() self._load_dump_structure() @unittest.skipIf(not mp.with_mpi(), "MPI specific test") - @unittest.skip("foo") + @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_structure_sharded(self): self._load_dump_structure(single_parallel_file=False) - @unittest.skip("foo") + @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_chunk_layout_file(self): self._load_dump_structure(chunk_file=True) - @unittest.skip("foo") + @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_chunk_layout_sim(self): self._load_dump_structure(chunk_sim=True) + # Tests dumping/loading of fields & structure. + def _load_dump_fields(self, single_parallel_file=True): + from meep.materials import Al + resolution = 50 + cell = mp.Vector3(5, 5) + sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.2), center=mp.Vector3(), component=mp.Ez) + one_by_one = mp.Vector3(1, 1, mp.inf) + geometry = [mp.Block(material=Al, center=mp.Vector3(), size=one_by_one), + mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)] + pml_layers = [mp.PML(0.5)] + + symmetries = [mp.Mirror(mp.Y)] + + sim1 = mp.Simulation(resolution=resolution, + cell_size=cell, + boundary_layers=pml_layers, + geometry=geometry, + symmetries=symmetries, + sources=[sources]) + + sample_point = mp.Vector3(0.12, -0.29) + ref_field_points = {} + + def get_ref_field_point(sim): + p = sim.get_field_point(mp.Ez, sample_point) + ref_field_points[sim.meep_time()] = p.real + + # First run until t=25 and save structure/fields + sim1.run(mp.at_every(5, get_ref_field_point), until=25) + + dump_dirname = os.path.join(self.temp_dir, 'test_load_dump_fields') + sim1.dump(dump_dirname, dump_structure=True, dump_fields=True, single_parallel_file=single_parallel_file) + + # Then continue running another 25 until t=50 + sim1.run(mp.at_every(5, get_ref_field_point), until=25) + print('ref_field_points = ' + str(ref_field_points)) + + # Now create a new simulation and try restoring state. + sim = mp.Simulation(resolution=resolution, + cell_size=cell, + boundary_layers=pml_layers, + sources=[sources], + symmetries=symmetries, + chunk_layout=sim1) + # Just restore structure first. + sim.load(dump_dirname, load_structure=True, load_fields=False, single_parallel_file=single_parallel_file) + field_points = {} + + def get_field_point(sim): + p = sim.get_field_point(mp.Ez, sample_point) + field_points[sim.meep_time()] = p.real + + # Run until t=1 so that fields are initialized. + sim.run(mp.at_every(5, get_field_point), until=1) + + # Now load the fields. + sim.load(dump_dirname, load_structure=False, load_fields=True, single_parallel_file=single_parallel_file) + sim.run(mp.at_every(5, get_field_point), until=25) + print('field_points = ' + str(field_points)) + + for t, v in field_points.items(): + self.assertAlmostEqual(ref_field_points[t], v) + + def test_load_dump_fields(self): + self._load_dump_fields() + + @unittest.skipIf(not mp.with_mpi(), "MPI specific test") + @unittest.skip("disable_temporarily_for_debugging") + def test_load_dump_fields_sharded(self): + self._load_dump_fields(single_parallel_file=False) + if __name__ == '__main__': unittest.main() diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index 79dd55380..7fd5d9d5c 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -112,6 +112,13 @@ void fields::dump(const char *filename, bool single_parallel_file) { h5file file(filename, h5file::WRITE, single_parallel_file, !single_parallel_file); + // Write out the current time 't' + size_t dims[1] = {1}; + size_t start[1] = {0}; + size_t time[1] = {(size_t)t}; + file.create_data("t", 1, dims); + file.write_chunk(1, start, dims, time); + dump_fields_chunk_field( &file, single_parallel_file, "f", [](fields_chunk *chunk, int c, int d) { return &(chunk->f[c][d]); }); @@ -222,6 +229,16 @@ void fields::load(const char *filename, bool single_parallel_file) { h5file file(filename, h5file::READONLY, single_parallel_file, !single_parallel_file); + // Read in the current time 't' + int rank; + size_t dims[1] = {1}; + size_t start[1] = {0}; + size_t time[1]; + file.read_size("t", &rank, dims, 1); + if (rank != 1 || dims[0] != 1) meep::abort("time size mismatch in fields::load"); + file.read_chunk(1, start, dims, time); + t = static_cast(time[0]); + load_fields_chunk_field( &file, single_parallel_file, "f", [](fields_chunk *chunk, int c, int d) { return &(chunk->f[c][d]); }); From 42b90e0f1413ec187bab71a74e1d2d4cf60edbf6 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Wed, 25 Aug 2021 01:37:20 +0000 Subject: [PATCH 09/17] Fix dft-chunk saving for the single-parallel-file mode. --- src/fields.cpp | 9 +++++++++ src/fields_dump.cpp | 25 ++++++++++++++++++------- src/meep.hpp | 1 + 3 files changed, 28 insertions(+), 7 deletions(-) diff --git a/src/fields.cpp b/src/fields.cpp index f5d0d9767..01117e968 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -718,6 +718,15 @@ void fields::unset_solve_cw_omega() { chunks[i]->unset_solve_cw_omega(); } +void fields::log(const char* prefix) { + master_printf("%sFields State:\n", prefix); + master_printf("%s a = %g, dt = %g\n", prefix, a, dt); + master_printf("%s m = %g, beta = %g\n", prefix, m, beta); + master_printf("%s t = %d, phasein_time = %d, is_real = %d\n", prefix, t, phasein_time, is_real); + master_printf("\n"); + master_printf("%s num_chunks = %d (shared=%d)\n", prefix, num_chunks, shared_chunks); +} + /* implement mirror boundary conditions for i outside 0..n-1: */ int mirrorindex(int i, int n) { return i >= n ? 2 * n - 1 - i : (i < 0 ? -1 - i : i); } diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index 7fd5d9d5c..1f804786b 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -108,6 +108,7 @@ void fields::dump_fields_chunk_field(h5file *h5f, bool single_parallel_file, void fields::dump(const char *filename, bool single_parallel_file) { if (verbosity > 0) { printf("creating fields output file \"%s\" (%d)...\n", filename, single_parallel_file); + log(); } h5file file(filename, h5file::WRITE, single_parallel_file, !single_parallel_file); @@ -115,9 +116,9 @@ void fields::dump(const char *filename, bool single_parallel_file) { // Write out the current time 't' size_t dims[1] = {1}; size_t start[1] = {0}; - size_t time[1] = {(size_t)t}; + size_t _t[1] = {(size_t)t}; file.create_data("t", 1, dims); - file.write_chunk(1, start, dims, time); + if (am_master() || !single_parallel_file) file.write_chunk(1, start, dims, _t); dump_fields_chunk_field( &file, single_parallel_file, "f", @@ -134,7 +135,7 @@ void fields::dump(const char *filename, bool single_parallel_file) { // Dump DFT chunks. for (int i = 0; i < num_chunks; i++) { - if (chunks[i]->is_mine()) { + if (single_parallel_file || chunks[i]->is_mine()) { char dataname[1024]; snprintf(dataname, 1024, "chunk%02d", i); save_dft_hdf5(chunks[i]->dft_chunks, dataname, &file, 0, single_parallel_file); @@ -233,11 +234,18 @@ void fields::load(const char *filename, bool single_parallel_file) { int rank; size_t dims[1] = {1}; size_t start[1] = {0}; - size_t time[1]; + size_t _t[1]; file.read_size("t", &rank, dims, 1); if (rank != 1 || dims[0] != 1) meep::abort("time size mismatch in fields::load"); - file.read_chunk(1, start, dims, time); - t = static_cast(time[0]); + if (am_master() || !single_parallel_file) file.read_chunk(1, start, dims, _t); + + if (single_parallel_file) { + file.prevent_deadlock(); + broadcast(0, _t, dims[0]); + } + + t = static_cast(_t[0]); + calc_sources(time()); load_fields_chunk_field( &file, single_parallel_file, "f", @@ -254,12 +262,15 @@ void fields::load(const char *filename, bool single_parallel_file) { // Load DFT chunks. for (int i = 0; i < num_chunks; i++) { - if (chunks[i]->is_mine()) { + if (single_parallel_file || chunks[i]->is_mine()) { char dataname[1024]; snprintf(dataname, 1024, "chunk%02d", i); load_dft_hdf5(chunks[i]->dft_chunks, dataname, &file, 0, single_parallel_file); } } + + if (verbosity > 0) + log(); } } // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index d522781f7..b74ead14a 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1720,6 +1720,7 @@ class fields { void remove_susceptibilities(); void remove_fluxes(); void reset(); + void log(const char* prefix = ""); // time.cpp std::vector time_spent_on(time_sink sink); From bdef5e86276489697e496c875c992acca0521127 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 9 Sep 2021 14:36:54 +0000 Subject: [PATCH 10/17] Abort when trying to dump fields with non-null polarization state --- src/fields_dump.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index 1f804786b..3c400be9f 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -50,6 +50,10 @@ void fields::dump_fields_chunk_field(h5file *h5f, bool single_parallel_file, size_t my_ntot = 0; for (int i = 0, chunk_i = 0; i < num_chunks; i++) { if (chunks[i]->is_mine()) { + FOR_FIELD_TYPES(ft) { + if (chunks[i]->pol[ft] != NULL) + meep::abort("non-null polarization_state in fields::dump (unsupported)"); + } size_t ntot = chunks[i]->gv.ntot(); for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) { for (int d = 0; d < 2; ++d) { From e500897dd099d6f891cbeecaa0bab6ec781dc1a2 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 9 Sep 2021 17:13:50 +0000 Subject: [PATCH 11/17] Clean up test_dump_load.py and add test_dump_fails_for_non_null_polarization_state test --- python/tests/test_dump_load.py | 65 +++++++++++++++++++++------------- 1 file changed, 41 insertions(+), 24 deletions(-) diff --git a/python/tests/test_dump_load.py b/python/tests/test_dump_load.py index ac6460ccb..5602eff4d 100644 --- a/python/tests/test_dump_load.py +++ b/python/tests/test_dump_load.py @@ -29,8 +29,7 @@ def setUpClass(cls): @classmethod def tearDownClass(cls): - # mp.delete_directory(cls.temp_dir) - pass + mp.delete_directory(cls.temp_dir) # Tests various combinations of dumping/loading structure & chunk layout. def _load_dump_structure(self, chunk_file=False, chunk_sim=False, single_parallel_file=True): @@ -91,34 +90,28 @@ def get_field_point(sim): for ref_pt, pt in zip(ref_field_points, field_points): self.assertAlmostEqual(ref_pt, pt) - @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_structure(self): self._load_dump_structure() @unittest.skipIf(not mp.with_mpi(), "MPI specific test") - @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_structure_sharded(self): self._load_dump_structure(single_parallel_file=False) - @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_chunk_layout_file(self): self._load_dump_structure(chunk_file=True) - @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_chunk_layout_sim(self): self._load_dump_structure(chunk_sim=True) # Tests dumping/loading of fields & structure. def _load_dump_fields(self, single_parallel_file=True): - from meep.materials import Al resolution = 50 cell = mp.Vector3(5, 5) - sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.2), center=mp.Vector3(), component=mp.Ez) + sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.4), center=mp.Vector3(), component=mp.Ez) one_by_one = mp.Vector3(1, 1, mp.inf) - geometry = [mp.Block(material=Al, center=mp.Vector3(), size=one_by_one), + geometry = [mp.Block(material=mp.Medium(index=3.2), center=mp.Vector3(), size=one_by_one), mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)] pml_layers = [mp.PML(0.5)] - symmetries = [mp.Mirror(mp.Y)] sim1 = mp.Simulation(resolution=resolution, @@ -129,21 +122,21 @@ def _load_dump_fields(self, single_parallel_file=True): sources=[sources]) sample_point = mp.Vector3(0.12, -0.29) - ref_field_points = {} + dump_dirname = os.path.join(self.temp_dir, 'test_load_dump_fields') + os.makedirs(dump_dirname, exist_ok=True) + + ref_field_points = {} def get_ref_field_point(sim): p = sim.get_field_point(mp.Ez, sample_point) ref_field_points[sim.meep_time()] = p.real - # First run until t=25 and save structure/fields - sim1.run(mp.at_every(5, get_ref_field_point), until=25) - - dump_dirname = os.path.join(self.temp_dir, 'test_load_dump_fields') + # First run until t=15 and save structure/fields + sim1.run(mp.at_every(1, get_ref_field_point), until=15) sim1.dump(dump_dirname, dump_structure=True, dump_fields=True, single_parallel_file=single_parallel_file) - # Then continue running another 25 until t=50 - sim1.run(mp.at_every(5, get_ref_field_point), until=25) - print('ref_field_points = ' + str(ref_field_points)) + # Then continue running another 5 until t=20 + sim1.run(mp.at_every(1, get_ref_field_point), until=5) # Now create a new simulation and try restoring state. sim = mp.Simulation(resolution=resolution, @@ -160,13 +153,11 @@ def get_field_point(sim): p = sim.get_field_point(mp.Ez, sample_point) field_points[sim.meep_time()] = p.real - # Run until t=1 so that fields are initialized. - sim.run(mp.at_every(5, get_field_point), until=1) + sim.init_sim() - # Now load the fields. + # Now load the fields (at t=15) and then continue to t=20 sim.load(dump_dirname, load_structure=False, load_fields=True, single_parallel_file=single_parallel_file) - sim.run(mp.at_every(5, get_field_point), until=25) - print('field_points = ' + str(field_points)) + sim.run(mp.at_every(1, get_field_point), until=5) for t, v in field_points.items(): self.assertAlmostEqual(ref_field_points[t], v) @@ -175,10 +166,36 @@ def test_load_dump_fields(self): self._load_dump_fields() @unittest.skipIf(not mp.with_mpi(), "MPI specific test") - @unittest.skip("disable_temporarily_for_debugging") def test_load_dump_fields_sharded(self): self._load_dump_fields(single_parallel_file=False) + # This assertRaisesRegex check does not play well with MPI due to the + # underlying call to meep::abort + @unittest.skipIf(mp.with_mpi(), "MPI specific test") + def test_dump_fails_for_non_null_polarization_state(self): + resolution = 50 + cell = mp.Vector3(5, 5) + sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.4), center=mp.Vector3(), component=mp.Ez) + one_by_one = mp.Vector3(1, 1, mp.inf) + from meep.materials import Al + geometry = [mp.Block(material=Al, center=mp.Vector3(), size=one_by_one), + mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)] + + sim = mp.Simulation(resolution=resolution, + cell_size=cell, + boundary_layers=[], + geometry=geometry, + symmetries=[], + sources=[sources]) + + dump_dirname = os.path.join(self.temp_dir, 'test_load_dump_fields') + os.makedirs(dump_dirname, exist_ok=True) + + sim.run(until=1) + # NOTE: We do not yet support checkpoint/restore when there is a + # non-null polarization_state + with self.assertRaisesRegex(RuntimeError, 'meep: non-null polarization_state in fields::dump'): + sim.dump(dump_dirname, dump_structure=True, dump_fields=True) if __name__ == '__main__': unittest.main() From cd5bd6c341c4369c38c44aa8f97bbce03ca49f5d Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 9 Sep 2021 17:18:30 +0000 Subject: [PATCH 12/17] Add new tests/dump_load to set of files to ignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 1a667d31c..4ed1ef7e1 100644 --- a/.gitignore +++ b/.gitignore @@ -90,6 +90,7 @@ tests/convergence_cyl_waveguide tests/cyl-ellipsoid-ll tests/cylindrical tests/dft-fields +tests/dump_load tests/flux tests/gdsII-3d tests/h5test From 004e303aed96035e2e3ed5eed3f2bac3d9527e13 Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 9 Sep 2021 19:22:43 +0000 Subject: [PATCH 13/17] Clean up the test to remove unnecessary stuff from the copied over test --- tests/dump_load.cpp | 215 ++++++++------------------------------------ 1 file changed, 39 insertions(+), 176 deletions(-) diff --git a/tests/dump_load.cpp b/tests/dump_load.cpp index 3adaa72c4..78685a4e9 100644 --- a/tests/dump_load.cpp +++ b/tests/dump_load.cpp @@ -124,65 +124,43 @@ int test_metal(double eps(const vec &), int splitting, const char *tmpdir) { double ttot = 17.0; grid_volume gv = vol3d(1.5, 0.5, 1.0, a); - structure s1(gv, eps); structure s(gv, eps, no_pml(), identity(), splitting); std::string filename_prefix = std::string(tmpdir) + "/test_metal_" + std::to_string(splitting); - std::string structure_filename_original = + std::string structure_filename = structure_dump(&s, filename_prefix, "original"); master_printf("Metal test using %d chunks...\n", splitting); fields f(&s); f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.299, 0.401), 1.0); - std::string fields_filename_original = - fields_dump(&f, filename_prefix, "original"); - - fields f1(&s1); - f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.299, 0.401), 1.0); - - double field_energy_check_time = 8.0; - while (f.time() < ttot) { - f.step(); - f1.step(); - if (!compare_point(f, f1, vec(0.5, 0.5, 0.01))) return 0; - if (!compare_point(f, f1, vec(0.46, 0.33, 0.33))) return 0; - if (!compare_point(f, f1, vec(1.301, 0.301, 0.399))) return 0; - if (f.time() >= field_energy_check_time) { - if (!compare(f.field_energy(), f1.field_energy(), " total energy")) - return 0; - if (!compare(f.electric_energy_in_box(gv.surroundings()), - f1.electric_energy_in_box(gv.surroundings()), - "electric energy")) - return 0; - if (!compare(f.magnetic_energy_in_box(gv.surroundings()), - f1.magnetic_energy_in_box(gv.surroundings()), - "magnetic energy")) - return 0; - field_energy_check_time += 5.0; - } - } + while (f.time() < ttot) f.step(); - std::string structure_filename_after_sim = - structure_dump(&s, filename_prefix, "after-sim"); - - std::string fields_filename_after_sim = - fields_dump(&f, filename_prefix, "after_sim"); + std::string fields_filename = + fields_dump(&f, filename_prefix, "original"); structure s_load(gv, eps, no_pml(), identity(), splitting); - structure_load(&s_load, structure_filename_after_sim); - - std::string structure_filename_dump_loaded = - structure_dump(&s, filename_prefix, "dump-loaded"); + structure_load(&s_load, structure_filename); fields f_load(&s_load); f_load.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.299, 0.401), 1.0); - fields_load(&f_load, fields_filename_after_sim); + fields_load(&f_load, fields_filename); - std::string fields_filename_dump_loaded = - fields_dump(&f_load, filename_prefix, "dump-loaded"); + if (!compare_point(f, f_load, vec(0.5, 0.5, 0.01))) return 0; + if (!compare_point(f, f_load, vec(0.46, 0.33, 0.33))) return 0; + if (!compare_point(f, f_load, vec(1.301, 0.301, 0.399))) return 0; + if (!compare(f.field_energy(), f_load.field_energy(), " total energy")) + return 0; + if (!compare(f.electric_energy_in_box(gv.surroundings()), + f_load.electric_energy_in_box(gv.surroundings()), + "electric energy")) + return 0; + if (!compare(f.magnetic_energy_in_box(gv.surroundings()), + f_load.magnetic_energy_in_box(gv.surroundings()), + "magnetic energy")) + return 0; return 1; } @@ -196,150 +174,41 @@ int test_periodic(double eps(const vec &), int splitting, const char *tmpdir) { std::string filename_prefix = std::string(tmpdir) + "/test_periodic_" + std::to_string(splitting); - std::string structure_filename_original = + std::string structure_filename = structure_dump(&s, filename_prefix, "original"); master_printf("Periodic test using %d chunks...\n", splitting); fields f(&s); f.use_bloch(vec(0.1, 0.7, 0.3)); f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.25, 0.5), 1.0); - std::string fields_filename_original = - fields_dump(&f, filename_prefix, "original"); - fields f1(&s1); - f1.use_bloch(vec(0.1, 0.7, 0.3)); - f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.25, 0.5), 1.0); - - double field_energy_check_time = 8.0; - while (f.time() < ttot) { - f.step(); - f1.step(); - if (!compare_point(f, f1, vec(0.5, 0.01, 0.5))) return 0; - if (!compare_point(f, f1, vec(0.46, 0.33, 0.2))) return 0; - if (!compare_point(f, f1, vec(1.0, 0.25, 0.301))) return 0; - if (f.time() >= field_energy_check_time) { - if (!compare(f.field_energy(), f1.field_energy(), " total energy")) - return 0; - if (!compare(f.electric_energy_in_box(gv.surroundings()), - f1.electric_energy_in_box(gv.surroundings()), - "electric energy")) - return 0; - if (!compare(f.magnetic_energy_in_box(gv.surroundings()), - f1.magnetic_energy_in_box(gv.surroundings()), - "magnetic energy")) - return 0; - field_energy_check_time += 5.0; - } - } - - std::string structure_filename_after_sim = - structure_dump(&s, filename_prefix, "after-sim"); + while (f.time() < ttot) f.step(); - std::string fields_filename_after_sim = - fields_dump(&f, filename_prefix, "after_sim"); + std::string fields_filename = + fields_dump(&f, filename_prefix, "original"); structure s_load(gv, eps, no_pml(), identity(), splitting); - structure_load(&s_load, structure_filename_after_sim); - - std::string structure_filename_dump_loaded = - structure_dump(&s, filename_prefix, "dump-loaded"); + structure_load(&s_load, structure_filename); fields f_load(&s_load); f_load.use_bloch(vec(0.1, 0.7, 0.3)); f_load.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.25, 0.5), 1.0); - fields_load(&f_load, fields_filename_after_sim); - - std::string fields_filename_dump_loaded = - fields_dump(&f_load, filename_prefix, "dump-loaded"); - return 1; -} - -int test_pml(double eps(const vec &), const char *tmpdir) { - master_printf("Testing pml quality...\n"); - - double a = 10.0; - grid_volume gv = vol3d(1.5, 1.0, 1.2, a); - structure s(gv, eps, pml(0.401)); - - std::string filename_prefix = std::string(tmpdir) + "/test_pml"; - std::string structure_filename_original = - structure_dump(&s, filename_prefix, "original"); - - fields f(&s); - f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.751, 0.5, 0.601), 1.0); - const double deltaT = 10.0; - const double ttot = 3.1 * deltaT; - double field_energy_check_time = deltaT; - - std::string fields_filename_original = - fields_dump(&f, filename_prefix, "original"); - - while (f.time() < f.last_source_time()) f.step(); - - double last_energy = f.field_energy(); - while (f.time() < ttot) { - f.step(); - if (f.time() >= field_energy_check_time) { - const double new_energy = f.field_energy(); - master_printf("Got newE/oldE of %g\n", new_energy / last_energy); - if (new_energy > last_energy * 4e-3) { - master_printf("Energy decaying too slowly: from %g to %g (%g)\n", - last_energy, new_energy, new_energy / last_energy); - return 0; - } - field_energy_check_time += deltaT; - } - } - - std::string structure_filename_after_sim = - structure_dump(&s, filename_prefix, "after-sim"); - - std::string fields_filename_after_sim = - fields_dump(&f, filename_prefix, "after_sim"); - - structure s_load(gv, eps, pml(0.401)); - structure_load(&s_load, structure_filename_after_sim); - - std::string structure_filename_dump_loaded = - structure_dump(&s, filename_prefix, "dump-loaded"); - - fields f_load(&s_load); - fields_load(&f_load, fields_filename_after_sim); - - std::string fields_filename_dump_loaded = - fields_dump(&f_load, filename_prefix, "dump-loaded"); - return 1; -} + fields_load(&f_load, fields_filename); -int test_pml_splitting(double eps(const vec &), int splitting, - const char *tmpdir) { - double a = 10.0; - - grid_volume gv = vol3d(1.5, 1.0, 1.2, a); - structure s1(gv, eps, pml(0.3)); - structure s(gv, eps, pml(0.3), identity(), splitting); + if (!compare_point(f, f_load, vec(0.5, 0.01, 0.5))) return 0; + if (!compare_point(f, f_load, vec(0.46, 0.33, 0.2))) return 0; + if (!compare_point(f, f_load, vec(1.0, 0.25, 0.301))) return 0; + if (!compare(f.field_energy(), f_load.field_energy(), " total energy")) + return 0; + if (!compare(f.electric_energy_in_box(gv.surroundings()), + f_load.electric_energy_in_box(gv.surroundings()), + "electric energy")) + return 0; + if (!compare(f.magnetic_energy_in_box(gv.surroundings()), + f_load.magnetic_energy_in_box(gv.surroundings()), + "magnetic energy")) + return 0; - master_printf("Testing pml while splitting into %d chunks...\n", splitting); - fields f(&s); - f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.099, 0.499, 0.501), 1.0); - fields f1(&s1); - f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.099, 0.499, 0.501), 1.0); - const double ttot = 31.0; - - double next_energy_time = 10.0; - while (f.time() < ttot) { - f.step(); - f1.step(); - if (!approx_point(f, f1, vec(0.5, 0.01, 1.0))) return 0; - if (!approx_point(f, f1, vec(0.46, 0.33, 0.33))) return 0; - if (!approx_point(f, f1, vec(1.0, 1.0, 0.33))) return 0; - if (!approx_point(f, f1, vec(1.3, 0.3, 0.15))) return 0; - if (f.time() > next_energy_time) { - if (!compare(f.field_energy(), f1.field_energy(), " total energy")) - return 0; - next_energy_time += 10.0; - } - } return 1; } @@ -350,8 +219,6 @@ int main(int argc, char **argv) { std::unique_ptr temp_dir(make_output_directory()); master_printf("Testing 3D dump/load: temp_dir = %s...\n", temp_dir.get()); - if (!test_pml(one, temp_dir.get())) abort("error in test_pml vacuum\n"); - for (int s = 2; s < 7; s++) if (!test_periodic(targets, s, temp_dir.get())) abort("error in test_periodic targets\n"); @@ -360,10 +227,6 @@ int main(int argc, char **argv) { if (!test_metal(one, s, temp_dir.get())) abort("error in test_metal vacuum\n"); - for (int s = 2; s < 4; s++) - if (!test_pml_splitting(one, s, temp_dir.get())) - abort("error in test_pml_splitting vacuum\n"); - - // delete_directory(temp_dir.get()); + delete_directory(temp_dir.get()); return 0; } From 6064894929ef418c8ad21f692b7491490f5e1c5f Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 9 Sep 2021 19:26:26 +0000 Subject: [PATCH 14/17] Also dump/load 'f_w_prev' --- src/fields_dump.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index 3c400be9f..bb416f54f 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -136,6 +136,9 @@ void fields::dump(const char *filename, bool single_parallel_file) { dump_fields_chunk_field( &file, single_parallel_file, "f_cond", [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); + dump_fields_chunk_field( + &file, single_parallel_file, "f_w_prev", + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_w_prev[c][d]); }); // Dump DFT chunks. for (int i = 0; i < num_chunks; i++) { @@ -263,6 +266,9 @@ void fields::load(const char *filename, bool single_parallel_file) { load_fields_chunk_field( &file, single_parallel_file, "f_cond", [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); + load_fields_chunk_field( + &file, single_parallel_file, "f_w_prev", + [](fields_chunk *chunk, int c, int d) { return &(chunk->w_prev[c][d]); }); // Load DFT chunks. for (int i = 0; i < num_chunks; i++) { From 6e3cb0d7b19bad5e753ff2bf7cc97040b1ca3c4b Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 9 Sep 2021 19:37:30 +0000 Subject: [PATCH 15/17] Fix typo causing build breakage --- src/fields_dump.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fields_dump.cpp b/src/fields_dump.cpp index bb416f54f..567a989e8 100644 --- a/src/fields_dump.cpp +++ b/src/fields_dump.cpp @@ -268,7 +268,7 @@ void fields::load(const char *filename, bool single_parallel_file) { [](fields_chunk *chunk, int c, int d) { return &(chunk->f_cond[c][d]); }); load_fields_chunk_field( &file, single_parallel_file, "f_w_prev", - [](fields_chunk *chunk, int c, int d) { return &(chunk->w_prev[c][d]); }); + [](fields_chunk *chunk, int c, int d) { return &(chunk->f_w_prev[c][d]); }); // Load DFT chunks. for (int i = 0; i < num_chunks; i++) { From 9294d6416c8718c3d226a34d912964f8cb917cdd Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 23 Sep 2021 19:18:48 +0000 Subject: [PATCH 16/17] load fields at the very end of init_sim after add_sources --- python/simulation.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 262987925..8de577910 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1996,10 +1996,6 @@ def init_sim(self): self.loop_tile_base ) - if self.load_fields_file: - self.load_fields( - self.load_fields_file, self.load_single_parallel_file) - if self.force_all_components and self.dimensions != 1: self.fields.require_component(mp.Ez) self.fields.require_component(mp.Hz) @@ -2019,7 +2015,11 @@ def init_sim(self): hook() self._is_initialized = True - + + if self.load_fields_file: + self.load_fields( + self.load_fields_file, self.load_single_parallel_file) + def using_real_fields(self): cond1 = self.is_cylindrical and self.m != 0 cond2 = any([s.phase.imag for s in self.symmetries]) From b8bcf9d69bb2b0a64318f68864ae85c1a192b7be Mon Sep 17 00:00:00 2001 From: Krishna Gadepalli Date: Thu, 23 Sep 2021 19:29:24 +0000 Subject: [PATCH 17/17] remove init_sim before loading fields since that now works. --- python/tests/test_dump_load.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/python/tests/test_dump_load.py b/python/tests/test_dump_load.py index 5602eff4d..95c843f72 100644 --- a/python/tests/test_dump_load.py +++ b/python/tests/test_dump_load.py @@ -153,8 +153,6 @@ def get_field_point(sim): p = sim.get_field_point(mp.Ez, sample_point) field_points[sim.meep_time()] = p.real - sim.init_sim() - # Now load the fields (at t=15) and then continue to t=20 sim.load(dump_dirname, load_structure=False, load_fields=True, single_parallel_file=single_parallel_file) sim.run(mp.at_every(1, get_field_point), until=5)