diff --git a/libmeepgeom/Makefile.am b/libmeepgeom/Makefile.am index 819a774ab..7a67b7485 100644 --- a/libmeepgeom/Makefile.am +++ b/libmeepgeom/Makefile.am @@ -11,9 +11,9 @@ AM_CPPFLAGS = -I$(top_srcdir)/src -DDATADIR=\"$(srcdir)/\" libmeepgeom_la_SOURCES = \ meepgeom.cpp meepgeom.hpp material_data.hpp -libmeepgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ +libmeepgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ -check_PROGRAMS = pw-source-ll ring-ll cyl-ellipsoid-ll absorber-1d-ll +check_PROGRAMS = pw-source-ll ring-ll cyl-ellipsoid-ll absorber-1d-ll array-slice-ll absorber_1d_ll_SOURCES = absorber-1d-ll.cpp absorber_1d_ll_LDADD = libmeepgeom.la $(MEEPLIBS) @@ -27,10 +27,14 @@ pw_source_ll_LDADD = libmeepgeom.la $(MEEPLIBS) ring_ll_SOURCES = ring-ll.cpp ring_ll_LDADD = libmeepgeom.la $(MEEPLIBS) -TESTS = cyl-ellipsoid-ll +array_slice_ll_SOURCES = array-slice-ll.cpp +array_slice_ll_LDADD = libmeepgeom.la $(MEEPLIBS) + +TESTS = cyl-ellipsoid-ll array-slice-ll noinst_PROGRAMS = bend-flux-ll + bend_flux_ll_SOURCES = bend-flux-ll.cpp bend_flux_ll_LDADD = libmeepgeom.la $(MEEPLIBS) -noinst_DATA = cyl-ellipsoid-eps-ref.h5 +noinst_DATA = cyl-ellipsoid-eps-ref.h5 array-slice-ll.h5 diff --git a/libmeepgeom/array-slice-ll.cpp b/libmeepgeom/array-slice-ll.cpp new file mode 100644 index 000000000..0040669ea --- /dev/null +++ b/libmeepgeom/array-slice-ll.cpp @@ -0,0 +1,249 @@ +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +#include +#include +#include +#include + +#include "meep.hpp" + +#include "ctl-math.h" +#include "ctlgeom.h" + +#include "meepgeom.hpp" + +#ifndef DATADIR + #define DATADIR "./" +#endif + +using namespace meep; + +typedef std::complex cdouble; + +vector3 v3(double x, double y=0.0, double z=0.0) +{ + vector3 v; + v.x=x; v.y=y; v.z=z; + return v; +} + +// passthrough field function +cdouble default_field_function(const cdouble *fields, + const vec &loc, void *data_) +{ + (void) loc; // unused + (void) data_; // unused + return fields[0]; +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +#define RELTOL 1.0e-6 +double Compare(double *d1, double *d2, int N, const char *Name) +{ + double Norm1=0.0, Norm2=0.0, NormDelta=0.0; + for(int n=0; n RELTOL) + abort("fail: rel error in %s data = %e\n",Name,RelErr); + return RelErr; +} + +double Compare(cdouble *d1, cdouble *d2, int N, const char *Name) +{ + double Norm1=0.0, Norm2=0.0, NormDelta=0.0; + for(int n=0; n RELTOL) + abort("fail: rel error in %s data = %e\n",Name,RelErr); + return RelErr; +} + +/***************************************************************/ +/* dummy material function needed to pass to structure( ) */ +/* constructor as a placeholder before we can call */ +/* set_materials_from_geometry */ +/***************************************************************/ +double dummy_eps(const vec &) { return 1.0; } + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +void usage(char *progname) +{ master_printf("usage: %s [options]\n",progname); + master_printf("options: \n"); + master_printf(" --use-symmetry use geometric symmetries\n"); + master_printf(" --write-files write reference data files\n"); + abort(); +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +int main(int argc, char *argv[], char **envp) +{ + initialize mpi(argc, argv); + + /***************************************************************/ + /* parse command-line options **********************************/ + /***************************************************************/ + bool use_symmetry=false; + bool write_files=false; + for(int narg=1; nargread("hz.r", &rank, dims1D, 1); + if (rank!=1 || dims1D[0]!=NX) + abort("failed to read 1D data(hz.r) from file %s.h5",H5FILENAME); + double *idata = file->read("hz.i", &rank, dims1D, 1); + if (rank!=1 || dims1D[0]!=NX) + abort("failed to read 1D data(hz.i) from file %s.h5",H5FILENAME); + file_slice1d = new cdouble[dims1D[0]]; + for(int n=0; nread("sy", &rank, dims2D, 2); + if (rank!=2 || dims2D[0]!=NX || dims2D[1]!=NY) + abort("failed to read 2D reference data from file %s.h5",H5FILENAME); + delete file; + + // + // generate 1D and 2D array slices and compare to + // data read from file + // + rank=f.get_array_slice_dimensions(v1d, dims1D); + if (rank!=1 || dims1D[0]!=NX) + abort("incorrect dimensions for 1D slice"); + cdouble *slice1d=f.get_complex_array_slice(v1d, Hz); + double RelErr1D=Compare(slice1d, file_slice1d, NX, "Hz_1d"); + master_printf("1D: rel error %e\n",RelErr1D); + + rank=f.get_array_slice_dimensions(v2d, dims2D); + if (rank!=2 || dims2D[0]!=NX || dims2D[1]!=NY) + abort("incorrect dimensions for 2D slice"); + double *slice2d=f.get_array_slice(v2d, Sy); + double RelErr2D=Compare(slice2d, file_slice2d, NX*NY, "Sy_2d"); + master_printf("2D: rel error %e\n",RelErr2D); + + }; // if (write_files) ... else ... + + return 0; +} diff --git a/libmeepgeom/array-slice-ll.h5 b/libmeepgeom/array-slice-ll.h5 new file mode 100644 index 000000000..9d359db83 Binary files /dev/null and b/libmeepgeom/array-slice-ll.h5 differ diff --git a/python/examples/cavity_arrayslice.py b/python/examples/cavity_arrayslice.py new file mode 100644 index 000000000..bc97e9686 --- /dev/null +++ b/python/examples/cavity_arrayslice.py @@ -0,0 +1,92 @@ +import unittest +import meep as mp +import numpy as np +import matplotlib.pyplot as plt + +################################################## +# set up the geometry +################################################## +eps = 13 +w = 1.2 +r = 0.36 +d = 1.4 +N = 3 +sy = 6 +pad = 2 +dpml = 1 +sx = (2 * (pad + dpml + N)) + d - 1 +fcen = 0.25 +df = 0.2 +nfreq = 500 + +cell = mp.Vector3(sx, sy, 0) + +blk = mp.Block(size=mp.Vector3(1e20, w, 1e20), + material=mp.Medium(epsilon=eps)) + +geometry = [blk] + +for i in range(3): + geometry.append(mp.Cylinder(r, center=mp.Vector3(d / 2 + i))) + +for i in range(3): + geometry.append(mp.Cylinder(r, center=mp.Vector3(d / -2 - i))) + +sim = mp.Simulation(cell_size=cell, + geometry=geometry, + sources=[], + boundary_layers=[mp.pml(dpml)], + resolution=20) + +################################################## +# add sources +################################################## +sim.sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), + mp.Hz, mp.Vector3())] + +#sim.symmetries = [mp.Mirror(mp.Y, phase=-1), mp.Mirror(mp.X, phase=-1)] + +################################################## +# run until sources are finished (and no later) +################################################## +sim._run_sources_until(0,[]) + +################################################## +# get 1D and 2D array slices +################################################## +xMin = -0.25*sx; +xMax = +0.25*sx; +yMin = -0.15*sy; +yMax = +0.15*sy; + +# 1D slice of Hz data +dims1d=np.zeros(3,dtype=np.int32) +v1d=mp.volume( mp.vec(xMin, 0.0), mp.vec(xMax, 0.0) ) +rank1d = sim.fields.get_array_slice_dimensions(v1d, dims1d); +NX1 = dims1d[0]; +slice1d = np.zeros(NX1, dtype=np.double); +sim.fields.get_array_slice(v1d, mp.Hz, slice1d); + +# 2D slice of Hz data +dims2d=np.zeros(3,dtype=np.int32) +v2d=mp.volume( mp.vec(xMin, yMin), mp.vec(xMax, yMax) ) +rank2d = sim.fields.get_array_slice_dimensions(v2d, dims2d); +NX2 = dims2d[0]; +NY2 = dims2d[1]; +N2 = NX2*NY2; +slice2d = np.zeros(N2, dtype=np.double); +sim.fields.get_array_slice(v2d, mp.Hz, slice2d); + +# plot 1D slice +plt.subplot(1,2,1); +x1d=np.linspace(xMin, xMax, dims1d[0]); +plt.plot(x1d, slice1d); + +# plot 2D slice +plt.subplot(1,2,2); +dy = (yMax-yMin) / dims2d[1]; +dx = (xMax-xMin) / dims2d[0]; +(x2d,y2d)=np.mgrid[ slice(xMin, xMax, dx), slice(yMin, yMax, dy)]; +plt.contourf( x2d, y2d, np.reshape(slice2d, (dims2d[0], dims2d[1])) ) +plt.colorbar(); +plt.show() diff --git a/python/meep.i b/python/meep.i index 55627b0c8..a2dd0230d 100644 --- a/python/meep.i +++ b/python/meep.i @@ -261,6 +261,7 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i SWIG_exception_fail(SWIG_ArgError(tmp_res), "Couldn't convert Python object to meep::src_time"); } $1 = reinterpret_cast(tmp_ptr); + } // Typemap suite for boundary_region @@ -295,6 +296,11 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i } } +// Typemap suite for array_slice +// TODO: add (cdouble *, int) version +%apply (double* INPLACE_ARRAY1, int DIM1) {(double *slice, int slice_length)}; +%apply int INPLACE_ARRAY1[ANY] { int [3] }; + // Rename python builtins %rename(br_apply) meep::boundary_region::apply; %rename(_is) meep::dft_chunk::is; diff --git a/src/Makefile.am b/src/Makefile.am index 5be09c990..3f5deb098 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -9,15 +9,15 @@ BUILT_SOURCES = sphere-quad.h step_generic_stride1.cpp HDRS = meep.hpp meep_internals.hpp meep/mympi.hpp meep/vec.hpp \ bicgstab.hpp -libmeep@MEEP_SUFFIX@_la_SOURCES = anisotropic_averaging.cpp bands.cpp \ -boundaries.cpp bicgstab.cpp casimir.cpp control_c.cpp cw_fields.cpp \ -dft.cpp dft_ldos.cpp energy_and_flux.cpp fields.cpp loop_in_chunks.cpp \ -h5fields.cpp h5file.cpp initialize.cpp integrate.cpp \ -integrate2.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 susceptibility.cpp time.cpp update_eh.cpp \ -mpb.cpp update_pols.cpp vec.cpp step_generic.cpp $(HDRS) \ -$(BUILT_SOURCES) +libmeep@MEEP_SUFFIX@_la_SOURCES = array_slice.cpp anisotropic_averaging.cpp \ +bands.cpp boundaries.cpp bicgstab.cpp casimir.cpp \ +control_c.cpp cw_fields.cpp dft.cpp dft_ldos.cpp energy_and_flux.cpp \ +fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp \ +initialize.cpp integrate.cpp integrate2.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 \ +susceptibility.cpp time.cpp update_eh.cpp mpb.cpp update_pols.cpp \ +vec.cpp step_generic.cpp $(HDRS) $(BUILT_SOURCES) libmeep@MEEP_SUFFIX@_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ diff --git a/src/array_slice.cpp b/src/array_slice.cpp new file mode 100644 index 000000000..b2656ec2e --- /dev/null +++ b/src/array_slice.cpp @@ -0,0 +1,482 @@ +/* Copyright (C) 2005-2015 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. +*/ + +/* HDF5 output of fields and arbitrary functions thereof. Works + very similarly to integrate.cpp (using fields::loop_in_chunks). */ + +#include +#include +#include + +#include "meep_internals.hpp" + +using namespace std; + +typedef complex cdouble; + +namespace meep { + +/***************************************************************************/ + +typedef struct { + + // information related to the volume covered by the + // array slice (its size, etcetera) + // these fields are filled in by get_array_slice_dimensions + // if the data parameter is non-null + ivec min_corner, max_corner; + int num_chunks; + int rank; + direction ds[3]; + int slice_size; + + // the function to output and related info (offsets for averaging, etc.) + // note: either fun *or* rfun should be non-NULL (not both) + field_function fun; + field_rfunction rfun; + void *fun_data; + std::vector components; + + void *vslice; + + // temporary internal storage buffers + component *cS; + cdouble *ph; + cdouble *fields; + int *offsets; + + int ninveps; + component inveps_cs[3]; + direction inveps_ds[3]; + + int ninvmu; + component invmu_cs[3]; + direction invmu_ds[3]; + +} array_slice_data; + +#define UNUSED(x) (void) x // silence compiler warnings + +/* passthrough field function equivalent to component_fun in h5fields.cpp */ +static cdouble default_field_func(const cdouble *fields, + const vec &loc, void *data_) +{ + (void) loc; // unused + (void) data_; // unused + return fields[0]; +} + +static double default_field_rfunc(const cdouble *fields, + const vec &loc, void *data_) +{ + (void) loc; // unused + (void) data_; // unused + return real(fields[0]); +} + +/***************************************************************/ +/* callback function passed to loop_in_chunks to compute */ +/* dimensions of array slice */ +/***************************************************************/ +static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, component cgrid, + ivec is, ivec ie, + vec s0, vec s1, vec e0, vec e1, + double dV0, double dV1, + ivec shift, complex shift_phase, + const symmetry &S, int sn, + void *data_) +{ + UNUSED(ichnk);UNUSED(cgrid);UNUSED(s0);UNUSED(s1);UNUSED(e0);UNUSED(e1); + UNUSED(dV0);UNUSED(dV1);UNUSED(shift_phase); UNUSED(fc); + array_slice_data *data = (array_slice_data *) data_; + ivec isS = S.transform(is, sn) + shift; + ivec ieS = S.transform(ie, sn) + shift; + data->min_corner = min(data->min_corner, min(isS, ieS)); + data->max_corner = max(data->max_corner, max(isS, ieS)); + data->num_chunks++; + +} + +/***************************************************************/ +/* callback function passed to loop_in_chunks to fill array slice */ +/***************************************************************/ +static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgrid, + ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, + double dV0, double dV1, + ivec shift, complex shift_phase, + const symmetry &S, int sn, void *data_) +{ + UNUSED(ichnk);UNUSED(cgrid);UNUSED(s0);UNUSED(s1);UNUSED(e0);UNUSED(e1); + UNUSED(dV0);UNUSED(dV1); + array_slice_data *data = (array_slice_data *) data_; + + //-----------------------------------------------------------------------// + // Find output chunk dimensions and strides, etc. + + int count[3]={1,1,1}, offset[3]={0,0,0}, stride[3]={1,1,1}; + + ivec isS = S.transform(is, sn) + shift; + ivec ieS = S.transform(ie, sn) + shift; + + // figure out what yucky_directions (in LOOP_OVER_IVECS) + // correspond to what directions in the transformed vectors (in output). + ivec permute(zero_ivec(fc->gv.dim)); + for (int i = 0; i < 3; ++i) + permute.set_direction(fc->gv.yucky_direction(i), i); + permute = S.transform_unshifted(permute, sn); + LOOP_OVER_DIRECTIONS(permute.dim, d) + permute.set_direction(d, abs(permute.in_direction(d))); + + // compute the size of the chunk to output, and its strides etc. + for (int i = 0; i < data->rank; ++i) { + direction d = data->ds[i]; + int isd = isS.in_direction(d), ied = ieS.in_direction(d); + count[i] = abs(ied - isd) / 2 + 1; + if (ied < isd) offset[permute.in_direction(d)] = count[i] - 1; + } + for (int i = 0; i < data->rank; ++i) { + direction d = data->ds[i]; + int j = permute.in_direction(d); + for (int k = i + 1; k < data->rank; ++k) stride[j] *= count[k]; + offset[j] *= stride[j]; + if (offset[j]) stride[j] *= -1; + } + + //-----------------------------------------------------------------------// + // Compute the function to output, exactly as in fields::integrate. + int *off = data->offsets; + component *cS = data->cS; + complex *fields = data->fields, *ph = data->ph; + const component *iecs = data->inveps_cs; + const direction *ieds = data->inveps_ds; + int ieos[6]; + const component *imcs = data->invmu_cs; + const direction *imds = data->invmu_ds; + int imos[6]; + int num_components=data->components.size(); + + double *slice=0; + cdouble *zslice=0; + bool complex_data = (data->rfun==0); + if (complex_data) + zslice = (cdouble *)data->vslice; + else + slice = (double *)data->vslice; + + for (int i = 0; i < num_components; ++i) { + cS[i] = S.transform(data->components[i], -sn); + if (cS[i] == Dielectric || cS[i] == Permeability) + ph[i] = 1.0; + else { + fc->gv.yee2cent_offsets(cS[i], off[2*i], off[2*i+1]); + ph[i] = shift_phase * S.phase_shift(cS[i], sn); + } + } + for (int k = 0; k < data->ninveps; ++k) + fc->gv.yee2cent_offsets(iecs[k], ieos[2*k], ieos[2*k+1]); + for (int k = 0; k < data->ninvmu; ++k) + fc->gv.yee2cent_offsets(imcs[k], imos[2*k], imos[2*k+1]); + + vec rshift(shift * (0.5*fc->gv.inva)); + LOOP_OVER_IVECS(fc->gv, is, ie, idx) { + IVEC_LOOP_LOC(fc->gv, loc); + loc = S.transform(loc, sn) + rshift; + + for (int i = 0; i < num_components; ++i) { + if (cS[i] == Dielectric) { + double tr = 0.0; + for (int k = 0; k < data->ninveps; ++k) { + const realnum *ie = fc->s->chi1inv[iecs[k]][ieds[k]]; + if (ie) tr += (ie[idx] + ie[idx+ieos[2*k]] + ie[idx+ieos[1+2*k]] + + ie[idx+ieos[2*k]+ieos[1+2*k]]); + else tr += 4; // default inveps == 1 + } + fields[i] = (4 * data->ninveps) / tr; + } + else if (cS[i] == Permeability) { + double tr = 0.0; + for (int k = 0; k < data->ninvmu; ++k) { + const realnum *im = fc->s->chi1inv[imcs[k]][imds[k]]; + if (im) tr += (im[idx] + im[idx+imos[2*k]] + im[idx+imos[1+2*k]] + + im[idx+imos[2*k]+imos[1+2*k]]); + else tr += 4; // default invmu == 1 + } + fields[i] = (4 * data->ninvmu) / tr; + } + else { + double f[2]; + for (int k = 0; k < 2; ++k) + if (fc->f[cS[i]][k]) + f[k] = 0.25 * (fc->f[cS[i]][k][idx] + + fc->f[cS[i]][k][idx+off[2*i]] + + fc->f[cS[i]][k][idx+off[2*i+1]] + + fc->f[cS[i]][k][idx+off[2*i]+off[2*i+1]]); + else + f[k] = 0; + fields[i] = complex(f[0], f[1]) * ph[i]; + } + } + int idx2 = ((((offset[0] + offset[1] + offset[2]) + + loop_i1 * stride[0]) + + loop_i2 * stride[1]) + + loop_i3 * stride[2]); + + if (complex_data) + zslice[idx2] = data->fun(fields, loc, data->fun_data); + else + slice[idx2] = data->rfun(fields, loc, data->fun_data); + + }; + +} + +/***************************************************************/ +/* given a volume, fill in the dims[] and directions[] arrays */ +/* describing the array slice needed to store field data for */ +/* all grid points in the volume. */ +/* */ +/* return value is rank of array slice. */ +/* */ +/* if caller_data is non-NULL, it should point to a */ +/* caller-allocated array_slice_data structure which will be */ +/* initialized appopriately for subsequent use in */ +/* get_array_slice. */ +/***************************************************************/ +int fields::get_array_slice_dimensions(const volume &where, int dims[3], void *caller_data) +{ + am_now_working_on(ArraySlice); + + // use a local data structure if the caller didn't provide one + array_slice_data local_data; + array_slice_data *data=(array_slice_data *)caller_data; + if (data==0) + data=&local_data; + + data->min_corner = gv.round_vec(where.get_max_corner()) + one_ivec(gv.dim); + data->max_corner = gv.round_vec(where.get_min_corner()) - one_ivec(gv.dim); + data->num_chunks = 0; + + loop_in_chunks(get_array_slice_dimensions_chunkloop, + (void *) data, where, Centered, true, true); + + data->max_corner = max_to_all(data->max_corner); + data->min_corner = -max_to_all(-data->min_corner); // i.e., min_to_all + data->num_chunks = sum_to_all(data->num_chunks); + if (data->num_chunks == 0 || !(data->min_corner <= data->max_corner)) + return 0; // no data to write; + + int rank=0, slice_size=1; + LOOP_OVER_DIRECTIONS(gv.dim, d) { + if (rank >= 3) abort("too many dimensions in array_slice"); + int n = (data->max_corner.in_direction(d) + - data->min_corner.in_direction(d)) / 2 + 1; + if (n > 1) { + data->ds[rank] = d; + dims[rank++] = n; + slice_size *= n; + } + } + data->rank=rank; + data->slice_size=slice_size; + finished_working(); + + return rank; +} + +/***************************************************************/ +/* precisely one of fun, rfun should be non-NULL */ +/***************************************************************/ +void *fields::do_get_array_slice(const volume &where, + std::vector components, + field_function fun, + field_rfunction rfun, + void *fun_data, + void *vslice) { + + am_now_working_on(ArraySlice); + + /***************************************************************/ + /* call get_array_slice_dimensions to get slice dimensions and */ + /* partially initialze an array_slice_data struct */ + /***************************************************************/ + int dims[3]; + array_slice_data data; + int rank=get_array_slice_dimensions(where, dims, &data); + int slice_size=data.slice_size; + if (rank==0 || slice_size==0) return 0; // no data to write + + bool complex_data = (rfun==0); + cdouble *zslice; + double *slice; + if (vslice==0) + { if (complex_data) + { zslice = new cdouble[slice_size]; + vslice = (void *)zslice; + } + else + { slice = new double[slice_size]; + vslice = (void *)slice; + }; + }; + + data.vslice = vslice; + data.fun = fun; + data.rfun = rfun; + data.fun_data = fun_data; + data.components = components; + + int num_components = components.size(); + + data.cS = new component[num_components]; + data.ph = new cdouble[num_components]; + data.fields = new cdouble[num_components]; + + data.offsets = new int[2 * num_components]; + for (int i = 0; i < 2 * num_components; ++i) + data.offsets[i] = 0; + + /* compute inverse-epsilon directions for computing Dielectric fields */ + data.ninveps = 0; + bool needs_dielectric = false; + for (int i = 0; i < num_components; ++i) + if (components[i] == Dielectric) { needs_dielectric = true; break; } + if (needs_dielectric) + FOR_ELECTRIC_COMPONENTS(c) if (gv.has_field(c)) { + if (data.ninveps == 3) abort("more than 3 field components??"); + data.inveps_cs[data.ninveps] = c; + data.inveps_ds[data.ninveps] = component_direction(c); + ++data.ninveps; + } + + /* compute inverse-mu directions for computing Permeability fields */ + data.ninvmu = 0; + bool needs_permeability = false; + for (int i = 0; i < num_components; ++i) + if (components[i] == Permeability) { needs_permeability = true; break; } + if (needs_permeability) + FOR_MAGNETIC_COMPONENTS(c) if (gv.has_field(c)) { + if (data.ninvmu == 3) abort("more than 3 field components??"); + data.invmu_cs[data.ninvmu] = c; + data.invmu_ds[data.ninvmu] = component_direction(c); + ++data.ninvmu; + } + + loop_in_chunks(get_array_slice_chunkloop, (void *) &data, + where, Centered, true, true); + + /***************************************************************/ + /* repeatedly call sum_to_all to consolidate full array slice */ + /* on all cores */ + /***************************************************************/ +#define BUFSIZE 1<<16 // use 64k buffer + if (complex_data) + { cdouble *buffer = new cdouble[BUFSIZE]; + cdouble *slice = (cdouble *)vslice; + int offset=0, remaining=slice_size; + while(remaining!=0) + { + int size = (remaining > BUFSIZE ? BUFSIZE : remaining); + sum_to_all(slice + offset, buffer, size); + memcpy(slice+offset, buffer, size*sizeof(cdouble)); + remaining-=size; + offset+=size; + }; + delete[] buffer; + } + else + { double *buffer = new double[BUFSIZE]; + double *slice = (double *)vslice; + int offset=0, remaining=slice_size; + while(remaining!=0) + { int size = (remaining > BUFSIZE ? BUFSIZE : remaining); + sum_to_all(slice + offset, buffer, size); + memcpy(slice+offset, buffer, size*sizeof(double)); + remaining-=size; + offset+=size; + }; + delete[] buffer; + }; + + delete[] data.offsets; + delete[] data.fields; + delete[] data.ph; + delete[] data.cS; + finished_working(); + + return vslice; +} + +/***************************************************************/ +/* entry points to get_array_slice */ +/***************************************************************/ +double *fields::get_array_slice(const volume &where, + std::vector components, + field_rfunction rfun, void *fun_data, + double *slice) +{ + return (double *)do_get_array_slice(where, components, + 0, rfun, fun_data, + (void *)slice); +} + +cdouble *fields::get_complex_array_slice(const volume &where, + std::vector components, + field_function fun, void *fun_data, + cdouble *slice) +{ + return (cdouble *)do_get_array_slice(where, components, + fun, 0, fun_data, + (void *)slice); +} + +double *fields::get_array_slice(const volume &where, component c, + double *slice, int slice_length) +{ + (void) slice_length; + std::vector components(1); + components[0]=c; + return (double *)do_get_array_slice(where, components, + 0, default_field_rfunc, 0, + (void *)slice); +} + +double *fields::get_array_slice(const volume &where, + derived_component c, + double *slice, int slice_length) +{ + (void) slice_length; + int nfields; + component carray[12]; + field_rfunction rfun = derived_component_func(c, gv, nfields, carray); + std::vector cs(carray, carray+nfields); + return (double *)do_get_array_slice(where, cs, + 0, rfun, &nfields, + (void *)slice); +} + +cdouble *fields::get_complex_array_slice(const volume &where, component c, + cdouble *slice, int slice_length) +{ + (void) slice_length; + std::vector components(1); + components[0]=c; + return (cdouble *)do_get_array_slice(where, components, + default_field_func, 0, 0, + (void *)slice); +} + +} // namespace meep diff --git a/src/h5file.cpp b/src/h5file.cpp index 4c3ac1dcd..5f7b89b9d 100644 --- a/src/h5file.cpp +++ b/src/h5file.cpp @@ -309,7 +309,7 @@ void h5file::read_size(const char *dataname, int *rank, int *dims, int maxrank) #define REALNUM_H5T (sizeof(realnum) == sizeof(double) ? H5T_NATIVE_DOUBLE : H5T_NATIVE_FLOAT) realnum *h5file::read(const char *dataname, - int *rank, int *dims, int maxrank) + int *rank, int *dims, int maxrank) { #ifdef HAVE_HDF5 realnum *data = 0; diff --git a/src/meep.hpp b/src/meep.hpp index 9c9d4da16..6774568b1 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -23,6 +23,8 @@ #include "meep/vec.hpp" #include "meep/mympi.hpp" +#include + namespace meep { /* We use the type realnum for large arrays, e.g. the fields. @@ -1165,7 +1167,8 @@ class fields_chunk { enum boundary_condition { Periodic=0, Metallic, Magnetic, None }; enum time_sink { Connecting, Stepping, Boundaries, MpiTime, - FieldOutput, FourierTransforming, Other }; + FieldOutput, FourierTransforming, + ArraySlice, Other }; typedef void (*field_chunkloop)(fields_chunk *fc, int ichunk, component cgrid, ivec is, ivec ie, @@ -1285,6 +1288,50 @@ class fields { const char *h5file_name(const char *name, const char *prefix = NULL, bool timestamp = false); + // array_slice.cpp methods + + // given a subvolume, compute the dimensions of the array slice + // needed to store field data for that subvolume. + // the data parameter is used internally in get_array_slice + // and should be ignored by external callers. + int get_array_slice_dimensions(const volume &where, int dims[3], void *data=0); + + // given a subvolume, return a column-major array containing + // the given function of the field components in that subvolume + // if slice is non-null, it must be a user-allocated buffer + // of the correct size. + // otherwise, a new buffer is allocated and returned; it + // must eventually be caller-deallocated via delete[]. + double *get_array_slice(const volume &where, + std::vector components, + field_rfunction rfun, void *fun_data, + double *slice=0); + + std::complex *get_complex_array_slice(const volume &where, + std::vector components, + field_function fun, + void *fun_data, + std::complex *slice=0); + + // alternative entry points for when you have no field + // function, i.e. you want just a single component or + // derived component.) (The slice_length parameter is only + // present to facilitate SWIG-generated python wrappers; + // it is ignored by the C code). + double *get_array_slice(const volume &where, component c, double *slice=0, int slice_length=0); + double *get_array_slice(const volume &where, derived_component c, double *slice=0, int slice_length=0); + std::complex *get_complex_array_slice(const volume &where, + component c, + std::complex *slice=0, int slice_length=0); + + // master routine for all above entry points + void *do_get_array_slice(const volume &where, + std::vector components, + field_function fun, + field_rfunction rfun, + void *fun_data, + void *vslice); + // step.cpp methods: double last_step_output_wall_time; int last_step_output_t; diff --git a/tests/Makefile.am b/tests/Makefile.am index e03ce8fb2..cbd1bcd74 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -74,7 +74,7 @@ h5test_LDADD = $(LIBMEEP) pml_SOURCES = pml.cpp pml_LDADD = $(LIBMEEP) -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 flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml LOG_COMPILER = $(RUNCODE) diff --git a/tests/h5test.cpp b/tests/h5test.cpp index 247f38d5b..00bb79c4a 100644 --- a/tests/h5test.cpp +++ b/tests/h5test.cpp @@ -185,7 +185,7 @@ bool check_2d(double eps(const vec &), double a, int splitting, symfunc Sf, delete[] h5data; } - file->remove(); + //file->remove(); delete file; master_printf("Passed %s (%g..%g), err=%g\n", name, @@ -300,7 +300,7 @@ bool check_3d(double eps(const vec &), double a, int splitting, symfunc Sf, delete[] h5data; } - file->remove(); + //file->remove(); delete file; master_printf("Passed %s (%g..%g), err=%g\n", name, @@ -384,7 +384,7 @@ bool check_2d_monitor(double eps(const vec &), delete[] mon; - file->remove(); + //file->remove(); delete file; master_printf("Passed %s (%g..%g), err=%g\n", name,