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 6d8a63795..b34526be4 100644
--- a/python/Makefile.am
+++ b/python/Makefile.am
@@ -64,14 +64,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 568995705..5485f6a94 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -1250,6 +1250,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:
@@ -1874,29 +1875,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()
@@ -1918,18 +1938,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.
@@ -1938,8 +1962,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:
@@ -1988,7 +2026,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 ae8e4f4e4..1f15ab3d3 100644
--- a/src/dft.cpp
+++ b/src/dft.cpp
@@ -377,18 +377,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,
@@ -405,13 +420,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,
@@ -432,8 +447,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 6d1ba792a..d9dcda7df 100644
--- a/src/fields.cpp
+++ b/src/fields.cpp
@@ -722,6 +722,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 <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <cassert>
+
+#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<size_t> 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<size_t> 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<size_t> 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<int>(_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 3caeb9366..9c1a95fd8 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -1199,10 +1199,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 {
@@ -1730,6 +1730,7 @@ class fields {
void remove_susceptibilities();
void remove_fluxes();
void reset();
+ void log(const char* prefix = "");
// time.cpp
std::vector<double> time_spent_on(time_sink sink);
@@ -1746,6 +1747,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,
@@ -2223,6 +2233,14 @@ class fields {
std::complex<double> A(const vec &), std::complex<double> amp,
component c0, direction d, int has_tm, int has_te);
+ // fields_dump.cpp
+ // Helper methods for dumping field chunks.
+ using FieldPtrGetter = std::function<realnum **(fields_chunk *, int, int)>;
+ 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<double> 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 <signal.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <meep.hpp>
+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<double> 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<double> 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<const char[]> 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;
+}