Skip to content

Commit

Permalink
enable single precision floating points for fields arrays (NanoComp#1544
Browse files Browse the repository at this point in the history
)

* remove cdouble type

* replace double,complex<double> with realnum,complex<realnum>

* always store DFT fields using double precision floating point

* fixes

* store all floating point parameters of MaterialGrid as doubles

* specify HDF5 read/write type format using single_precision parameter

* read_chunk/write_chunk in structure_dump.cpp use correct precision

* switch dft_ldos arrays to type complex<double> from complex<realnum>

* convert PML arrays sig,kap,siginv to realnum from double

* use single precision by default for h5file member functions write/read and related

* create two versions of read_chunk/write_chunk with input arrays of type float and double

* convert gyrotropic_susceptibility parameters to realnum

* minor formatting fixes

* update comment for single precision usage in meep.hpp
  • Loading branch information
oskooi authored Apr 14, 2021
1 parent 007b90d commit 2ea29d5
Show file tree
Hide file tree
Showing 32 changed files with 543 additions and 532 deletions.
33 changes: 15 additions & 18 deletions python/meep.i
Original file line number Diff line number Diff line change
Expand Up @@ -322,17 +322,14 @@ PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) {
}

// Wrapper around meep::dft_near2far::get_farfields_array
PyObject *_get_farfields_array(meep::dft_near2far *n2f, const meep::volume &where,
double resolution) {
PyObject *_get_farfields_array(meep::dft_near2far *n2f, const meep::volume &where,
double resolution) {
// Return value: New reference
size_t dims[4] = {1, 1, 1, 1};
int rank = 0;
size_t N = 1;

// TODO: Support single precision?
if (sizeof(realnum) == sizeof(float)) abort("Single precision not supported for get_farfields");

meep::realnum *EH = n2f->get_farfields_array(where, rank, dims, N, resolution);
double *EH = n2f->get_farfields_array(where, rank, dims, N, resolution);

if (!EH) return PyArray_SimpleNew(0, 0, NPY_CDOUBLE);

Expand All @@ -348,7 +345,7 @@ PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) {
}

PyObject *py_arr = PyArray_SimpleNew(rank, arr_dims, NPY_DOUBLE);
memcpy(PyArray_DATA((PyArrayObject*)py_arr), EH, sizeof(meep::realnum) * 2 * N * 6 * n2f->freq.size());
memcpy(PyArray_DATA((PyArrayObject*)py_arr), EH, sizeof(double) * 2 * N * 6 * n2f->freq.size());

delete[] arr_dims;
delete[] EH;
Expand Down Expand Up @@ -455,7 +452,7 @@ size_t _get_dft_data_size(meep::dft_chunk *dc) {
return meep::dft_chunks_Ntotal(dc, &istart) / 2;
}

void _get_dft_data(meep::dft_chunk *dc, std::complex<meep::realnum> *cdata, int size) {
void _get_dft_data(meep::dft_chunk *dc, std::complex<double> *cdata, int size) {
size_t istart;
size_t n = meep::dft_chunks_Ntotal(dc, &istart) / 2;
istart /= 2;
Expand All @@ -473,7 +470,7 @@ void _get_dft_data(meep::dft_chunk *dc, std::complex<meep::realnum> *cdata, int
}
}

void _load_dft_data(meep::dft_chunk *dc, std::complex<meep::realnum> *cdata, int size) {
void _load_dft_data(meep::dft_chunk *dc, std::complex<double> *cdata, int size) {
size_t istart;
size_t n = meep::dft_chunks_Ntotal(dc, &istart) / 2;
istart /= 2;
Expand Down Expand Up @@ -597,10 +594,10 @@ void _get_eigenmode(meep::fields *f, double frequency, meep::direction d, const
#endif
%}

%numpy_typemaps(std::complex<meep::realnum>, NPY_CDOUBLE, int);
%numpy_typemaps(std::complex<double>, NPY_CDOUBLE, int);
%numpy_typemaps(std::complex<double>, NPY_CDOUBLE, size_t);

%apply (std::complex<meep::realnum> *INPLACE_ARRAY1, int DIM1) {(std::complex<meep::realnum> *cdata, int size)};
%apply (std::complex<double> *INPLACE_ARRAY1, int DIM1) {(std::complex<double> *cdata, int size)};

// add_volume_source
%apply (std::complex<double> *INPLACE_ARRAY3, size_t DIM1, size_t DIM2, size_t DIM3) {
Expand Down Expand Up @@ -630,8 +627,8 @@ PyObject *_dft_ldos_J(meep::dft_ldos *f);
template<typename dft_type>
PyObject *_get_dft_array(meep::fields *f, dft_type dft, meep::component c, int num_freq);
size_t _get_dft_data_size(meep::dft_chunk *dc);
void _get_dft_data(meep::dft_chunk *dc, std::complex<meep::realnum> *cdata, int size);
void _load_dft_data(meep::dft_chunk *dc, std::complex<meep::realnum> *cdata, int size);
void _get_dft_data(meep::dft_chunk *dc, std::complex<double> *cdata, int size);
void _load_dft_data(meep::dft_chunk *dc, std::complex<double> *cdata, int size);
meep::volume_list *make_volume_list(const meep::volume &v, int c,
std::complex<double> weight,
meep::volume_list *next);
Expand Down Expand Up @@ -831,22 +828,22 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
if (!PyArray_Check(pao_grad)) meep::abort("grad parameter must be numpy array.");
if (!PyArray_ISCARRAY(pao_grad)) meep::abort("Numpy grad array must be C-style contiguous.");
if (PyArray_NDIM(pao_grad) !=2) {meep::abort("Numpy grad array must have 2 dimensions.");}
meep::realnum * grad_c = (meep::realnum *)PyArray_DATA(pao_grad);
double *grad_c = (double *)PyArray_DATA(pao_grad);
int ng = PyArray_DIMS(pao_grad)[1]; // number of design parameters

// clean the adjoint fields array
PyArrayObject *pao_fields_a = (PyArrayObject *)fields_a;
if (!PyArray_Check(pao_fields_a)) meep::abort("adjoint fields parameter must be numpy array.");
if (!PyArray_ISCARRAY(pao_fields_a)) meep::abort("Numpy adjoint fields array must be C-style contiguous.");
if (PyArray_NDIM(pao_fields_a) !=1) {meep::abort("Numpy adjoint fields array must have 1 dimension.");}
std::complex<double> * fields_a_c = (std::complex<double> *)PyArray_DATA(pao_fields_a);
std::complex<double> *fields_a_c = (std::complex<double> *)PyArray_DATA(pao_fields_a);

// clean the forward fields array
PyArrayObject *pao_fields_f = (PyArrayObject *)fields_f;
if (!PyArray_Check(pao_fields_f)) meep::abort("forward fields parameter must be numpy array.");
if (!PyArray_ISCARRAY(pao_fields_f)) meep::abort("Numpy forward fields array must be C-style contiguous.");
if (PyArray_NDIM(pao_fields_f) !=1) {meep::abort("Numpy forward fields array must have 1 dimension.");}
std::complex<double> * fields_f_c = (std::complex<double> *)PyArray_DATA(pao_fields_f);
std::complex<double> *fields_f_c = (std::complex<double> *)PyArray_DATA(pao_fields_f);

// scalegrad not currently used
double scalegrad = 1.0;
Expand All @@ -862,12 +859,12 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
PyArrayObject *pao_freqs = (PyArrayObject *)frequencies;
if (!PyArray_Check(pao_freqs)) meep::abort("frequencies parameter must be numpy array.");
if (!PyArray_ISCARRAY(pao_freqs)) meep::abort("Numpy fields array must be C-style contiguous.");
meep::realnum* frequencies_c = (meep::realnum *)PyArray_DATA(pao_freqs);
double *frequencies_c = (double *)PyArray_DATA(pao_freqs);
int nf = PyArray_DIMS(pao_freqs)[0];
if (PyArray_DIMS(pao_grad)[0] != nf) meep::abort("Numpy grad array is allocated for %d frequencies; it should be allocated for %d.",PyArray_DIMS(pao_grad)[0],nf);

// prepare a geometry_tree
//TODO eventually it would be nice to store the geom tree within the structure object so we don't have to recreate it here
// TODO eventually it would be nice to store the geom tree within the structure object so we don't have to recreate it here
geometric_object_list *l;
l = new geometric_object_list();
if (!py_list_to_gobj_list(py_geom_list,l)) meep::abort("Unable to convert geometry tree.");
Expand Down
8 changes: 4 additions & 4 deletions python/typemap_utils.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -488,8 +488,8 @@ static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) {
if (!PyArray_ISCARRAY(pao)) {
meep::abort("Numpy array weights must be C-style contiguous.");
}
md->weights = new realnum[PyArray_SIZE(pao)];
memcpy(md->weights, (realnum *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(realnum));
md->weights = new double[PyArray_SIZE(pao)];
memcpy(md->weights, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));

// if needed, combine sus structs to main object
PyObject *py_e_sus_m1 = PyObject_GetAttrString(po_medium1, "E_susceptibilities");
Expand Down Expand Up @@ -567,8 +567,8 @@ static int pymaterial_to_material(PyObject *po, material_type *mt) {
md = new material_data();
md->which_subclass = material_data::MATERIAL_FILE;
md->epsilon_dims[0] = md->epsilon_dims[1] = md->epsilon_dims[2] = 1;
md->epsilon_data = new realnum[PyArray_SIZE(pao)];
memcpy(md->epsilon_data, (realnum *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(realnum));
md->epsilon_data = new double[PyArray_SIZE(pao)];
memcpy(md->epsilon_data, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));

for (int i = 0; i < PyArray_NDIM(pao); ++i) {
md->epsilon_dims[i] = (size_t)PyArray_DIMS(pao)[i];
Expand Down
4 changes: 2 additions & 2 deletions scheme/structure.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@ static geom_box gv2box(const meep::volume &v) {

/***********************************************************************/

static meep::realnum *epsilon_data = NULL;
static double *epsilon_data = NULL;
static size_t epsilon_dims[3] = {0, 0, 0};

static void read_epsilon_file(const char *eps_input_file) {
Expand All @@ -175,7 +175,7 @@ static void read_epsilon_file(const char *eps_input_file) {
if (dataname) *(dataname++) = 0;
meep::h5file eps_file(fname, meep::h5file::READONLY, false);
int rank; // ignored since rank < 3 is equivalent to singleton dims
epsilon_data = eps_file.read(dataname, &rank, epsilon_dims, 3);
epsilon_data = (double *)eps_file.read(dataname, &rank, epsilon_dims, 3, false);
master_printf("read in %zdx%zdx%zd epsilon-input-file \"%s\"\n", epsilon_dims[0],
epsilon_dims[1], epsilon_dims[2], eps_input_file);
delete[] fname;
Expand Down
58 changes: 28 additions & 30 deletions src/array_slice.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,6 @@

using namespace std;

typedef complex<double> cdouble;

namespace meep {

/***************************************************************************/
Expand Down Expand Up @@ -62,8 +60,8 @@ typedef struct {

// temporary internal storage buffers
component *cS;
cdouble *ph;
cdouble *fields;
complex<double> *ph;
complex<double> *fields;
ptrdiff_t *offsets;

double frequency;
Expand All @@ -83,13 +81,13 @@ typedef struct {
#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_) {
static complex<double> default_field_func(const complex<double> *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_) {
static double default_field_rfunc(const complex<double> *fields, const vec &loc, void *data_) {
(void)loc; // unused
(void)data_; // unused
return real(fields[0]);
Expand Down Expand Up @@ -138,7 +136,7 @@ static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, co
typedef struct {
component source_component;
ivec slice_imin, slice_imax;
cdouble *slice;
complex<double> *slice;
} source_slice_data;

bool in_range(int imin, int i, int imax) { return (imin <= i && i <= imax); }
Expand Down Expand Up @@ -186,7 +184,7 @@ static void get_source_slice_chunkloop(fields_chunk *fc, int ichnk, component cg
// local index of the symmetry-child grid point within this
// slice (that is, if it even lies within the slice)
for (size_t npt = 0; npt < s->npts; npt++) {
cdouble amp = s->A[npt];
complex<double> amp = s->A[npt];
ptrdiff_t chunk_index = s->index[npt];
ivec iloc_parent = fc->gv.iloc(Dielectric, chunk_index);
ivec iloc_child = S.transform(iloc_parent, sn) + shift;
Expand Down Expand Up @@ -275,10 +273,10 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
// tabulated on the slice, exactly as in fields::integrate. //
//-----------------------------------------------------------------------//
double *slice = 0;
cdouble *zslice = 0;
complex<double> *zslice = 0;
bool complex_data = (data->rfun == 0);
if (complex_data)
zslice = (cdouble *)data->vslice;
zslice = (complex<double> *)data->vslice;
else
slice = (double *)data->vslice;

Expand Down Expand Up @@ -410,8 +408,8 @@ double *array_to_all(double *array, size_t array_size) {
return array;
}

cdouble *array_to_all(cdouble *array, size_t array_size) {
return (cdouble *)array_to_all((double *)array, 2 * array_size);
complex<double> *array_to_all(complex<double> *array, size_t array_size) {
return (complex<double> *)array_to_all((double *)array, 2 * array_size);
}

/***************************************************************/
Expand Down Expand Up @@ -569,9 +567,9 @@ double *collapse_array(double *array, int *rank, size_t dims[3], direction dirs[
return reduced_array;
}

cdouble *collapse_array(cdouble *array, int *rank, size_t dims[3], direction dirs[3],
volume where) {
return (cdouble *)collapse_array((double *)array, rank, dims, dirs, where, 2);
complex<double> *collapse_array(complex<double> *array, int *rank, size_t dims[3], direction dirs[3],
volume where) {
return (complex<double> *)collapse_array((double *)array, rank, dims, dirs, where, 2);
}

/**********************************************************************/
Expand Down Expand Up @@ -606,8 +604,8 @@ void *fields::do_get_array_slice(const volume &where, std::vector<component> com
data.frequency = frequency;
int num_components = components.size();
data.cS = new component[num_components];
data.ph = new cdouble[num_components];
data.fields = new cdouble[num_components];
data.ph = new complex<double>[num_components];
data.fields = new complex<double>[num_components];
data.offsets = new ptrdiff_t[2 * num_components];
memset(data.offsets, 0, 2 * num_components * sizeof(ptrdiff_t));
data.empty_dim[0] = data.empty_dim[1] = data.empty_dim[2] = data.empty_dim[3] =
Expand Down Expand Up @@ -680,11 +678,11 @@ double *fields::get_array_slice(const volume &where, std::vector<component> comp
frequency, snap);
}

cdouble *fields::get_complex_array_slice(const volume &where, std::vector<component> components,
field_function fun, void *fun_data, cdouble *slice,
double frequency, bool snap) {
return (cdouble *)do_get_array_slice(where, components, fun, 0, fun_data, (void *)slice,
frequency, snap);
complex<double> *fields::get_complex_array_slice(const volume &where, std::vector<component> components,
field_function fun, void *fun_data, complex<double> *slice,
double frequency, bool snap) {
return (complex<double> *)do_get_array_slice(where, components, fun, 0, fun_data, (void *)slice,
frequency, snap);
}

double *fields::get_array_slice(const volume &where, component c, double *slice, double frequency,
Expand All @@ -705,16 +703,16 @@ double *fields::get_array_slice(const volume &where, derived_component c, double
frequency, snap);
}

cdouble *fields::get_complex_array_slice(const volume &where, component c, cdouble *slice,
double frequency, bool snap) {
complex<double> *fields::get_complex_array_slice(const volume &where, component c, complex<double> *slice,
double frequency, bool snap) {
std::vector<component> components(1);
components[0] = c;
return (cdouble *)do_get_array_slice(where, components, default_field_func, 0, 0, (void *)slice,
return (complex<double> *)do_get_array_slice(where, components, default_field_func, 0, 0, (void *)slice,
frequency, snap);
}

cdouble *fields::get_source_slice(const volume &where, component source_slice_component,
cdouble *slice) {
complex<double> *fields::get_source_slice(const volume &where, component source_slice_component,
complex<double> *slice) {
size_t dims[3];
direction dirs[3];
vec min_max_loc[2];
Expand All @@ -725,18 +723,18 @@ cdouble *fields::get_source_slice(const volume &where, component source_slice_co
data.source_component = source_slice_component;
data.slice_imin = gv.round_vec(min_max_loc[0]);
data.slice_imax = gv.round_vec(min_max_loc[1]);
data.slice = new cdouble[slice_size];
data.slice = new complex<double>[slice_size];
if (!data.slice) abort("%s:%i: out of memory (%zu)", __FILE__, __LINE__, slice_size);

loop_in_chunks(get_source_slice_chunkloop, (void *)&data, where, Centered, true, false);

cdouble *slice_collapsed = collapse_array(data.slice, &rank, dims, dirs, where);
complex<double> *slice_collapsed = collapse_array(data.slice, &rank, dims, dirs, where);
rank = get_array_slice_dimensions(where, dims, dirs, true, false);
slice_size = dims[0] * (rank >= 2 ? dims[1] : 1) * (rank == 3 ? dims[2] : 1);

if (slice) {
memcpy(slice, slice_collapsed, 2 * slice_size * sizeof(double));
delete[] (cdouble*) slice_collapsed;
delete[] (complex<double>*) slice_collapsed;
}
else
slice = slice_collapsed;
Expand Down
Loading

0 comments on commit 2ea29d5

Please sign in to comment.