Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

enable single precision floating points for fields arrays #1544

Merged
merged 14 commits into from
Apr 14, 2021
Merged

enable single precision floating points for fields arrays #1544

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
584f584
remove cdouble type
oskooi Mar 31, 2021
fc603a3
replace double,complex<double> with realnum,complex<realnum>
oskooi Apr 1, 2021
abaad41
always store DFT fields using double precision floating point
oskooi Apr 2, 2021
3254c69
fixes
oskooi Apr 2, 2021
a1ba1bd
store all floating point parameters of MaterialGrid as doubles
oskooi Apr 4, 2021
99ec5a3
specify HDF5 read/write type format using single_precision parameter
oskooi Apr 4, 2021
a6673a4
read_chunk/write_chunk in structure_dump.cpp use correct precision
oskooi Apr 5, 2021
c800556
switch dft_ldos arrays to type complex<double> from complex<realnum>
oskooi Apr 5, 2021
01b6a2a
convert PML arrays sig,kap,siginv to realnum from double
oskooi Apr 8, 2021
394c359
use single precision by default for h5file member functions write/rea…
oskooi Apr 8, 2021
00c9e5e
create two versions of read_chunk/write_chunk with input arrays of ty…
oskooi Apr 8, 2021
ce17449
convert gyrotropic_susceptibility parameters to realnum
oskooi Apr 8, 2021
9beea0f
minor formatting fixes
oskooi Apr 13, 2021
3d67879
update comment for single precision usage in meep.hpp
oskooi Apr 14, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 15 additions & 18 deletions python/meep.i
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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