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 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..8de577910 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,48 @@ 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) + print("Dumped structure to file: %s (%s)" % (fname, str(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) + print("Loaded structure from file: %s (%s)" % (fname, str(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) + print("Dumped fields to file: %s (%s)" % (fname, str(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) + print("Loaded fields from file: %s (%s)" % (fname, str(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,18 +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, 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) - def load(self, dirname, structure=True, single_parallel_file=True): + 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, load_structure=True, load_fields=True, single_parallel_file=True): """ Loads simulation state. @@ -1927,8 +1951,22 @@ 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 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: @@ -1977,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]) diff --git a/python/tests/test_dump_load.py b/python/tests/test_dump_load.py new file mode 100644 index 000000000..95c843f72 --- /dev/null +++ b/python/tests/test_dump_load.py @@ -0,0 +1,199 @@ +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() + print("Saving temp files to dir: {}".format(cls.temp_dir)) + + @classmethod + def tearDownClass(cls): + 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): + 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_structure') + sim1.dump(dump_dirname, dump_structure=True, dump_fields=False, 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, 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.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) + + # Tests dumping/loading of fields & structure. + def _load_dump_fields(self, single_parallel_file=True): + 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) + 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, + cell_size=cell, + boundary_layers=pml_layers, + geometry=geometry, + symmetries=symmetries, + sources=[sources]) + + sample_point = mp.Vector3(0.12, -0.29) + + 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=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 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, + 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 + + # 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) + + 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") + 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() 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/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.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 new file mode 100644 index 000000000..567a989e8 --- /dev/null +++ b/src/fields_dump.cpp @@ -0,0 +1,286 @@ +/* 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++) { + 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) { + realnum **f = field_ptr_getter(chunks[i], c, d); + if (*f) { + my_ntot += (num_f_[(chunk_i * NUM_FIELD_COMPONENTS + c) * 2 + d] = 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); + log(); + } + + 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 _t[1] = {(size_t)t}; + file.create_data("t", 1, dims); + if (am_master() || !single_parallel_file) file.write_chunk(1, start, dims, _t); + + 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]); }); + 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++) { + 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); + } + } +} + +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++) { + 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; + } + } + } + } + chunk_i += (chunks[i]->is_mine() || single_parallel_file); + } + + /* 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); + } + + /* 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); + } + 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; + } + } + } + } + } +} + +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); + + // Read in the current time 't' + int rank; + size_t dims[1] = {1}; + size_t start[1] = {0}; + 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"); + 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", + [](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]); }); + load_fields_chunk_field( + &file, single_parallel_file, "f_w_prev", + [](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++) { + 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 ec5b8db96..b74ead14a 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 { @@ -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); @@ -1736,6 +1737,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, @@ -2209,6 +2219,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..78685a4e9 --- /dev/null +++ b/tests/dump_load.cpp @@ -0,0 +1,232 @@ +/* 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 s(gv, eps, no_pml(), identity(), splitting); + + std::string filename_prefix = + std::string(tmpdir) + "/test_metal_" + std::to_string(splitting); + 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); + + while (f.time() < ttot) f.step(); + + 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); + + 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); + + 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; +} + +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 = + 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); + + while (f.time() < ttot) f.step(); + + 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); + + 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); + + 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; + + 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()); + + 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"); + + delete_directory(temp_dir.get()); + return 0; +}