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

improved source-slice implementation that handles symmetries #820

Closed
wants to merge 1 commit into from
Closed

improved source-slice implementation that handles symmetries #820

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
197 changes: 123 additions & 74 deletions src/array_slice.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,14 +46,17 @@ typedef struct {
direction ds[3];
size_t slice_size;

// if non-null, min_max_loc[0,1] are filled in by get_array_slice_dimensions_chunkloop
// with the (coordinate-wise) minimum and maximum grid points encountered
// in looping over the slice region.
vec *min_max_loc;

// 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<component> components;
component source_slice_component;
bool get_source_slice;

void *vslice;

Expand Down Expand Up @@ -113,48 +116,78 @@ static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, co
data->min_corner = min(data->min_corner, min(isS, ieS));
data->max_corner = max(data->max_corner, max(isS, ieS));
data->num_chunks++;

if (data->min_max_loc) {
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;
data->min_max_loc[0] = min(loc, data->min_max_loc[0]);
data->min_max_loc[1] = max(loc, data->min_max_loc[1]);
}
}
}

/*****************************************************************/
/* populate the array slice with information about sources with */
/* the component specified in source_slice_component. */
/*****************************************************************/
void fill_chunk_source_slice(fields_chunk *fc, array_slice_data *data) {
ndim dim = fc->gv.dim;
if (dim == Dcyl) {
fprintf(stderr, "warning: source slices not implemented for cylindrical coordinates; array "
"will be all zeros\n");
return;
}
/***************************************************************/
/* chunkloop for get_source_slice ******************************/
/***************************************************************/
typedef struct {
component source_component;
ivec slice_imin, slice_imax;
cdouble *slice;
} source_slice_data;

bool in_range(int imin, int i, int imax)
{ return (imin<=i && i<=imax); }

bool in_subgrid(ivec ivmin, ivec iv, ivec ivmax)
{ LOOP_OVER_DIRECTIONS(iv.dim,d)
if ( !in_range(ivmin.in_direction(d),iv.in_direction(d),ivmax.in_direction(d)) )
return false;
return true;
}

static void get_source_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<double> shift_phase, const symmetry &S,
int sn, void *data_) {

component slice_component = data->source_slice_component;
cdouble *slice = (cdouble *)data->vslice;
ivec slice_imin = data->min_corner, slice_imax = data->max_corner;
UNUSED(ichnk); UNUSED(cgrid); UNUSED(is); UNUSED(ie); UNUSED(s0); UNUSED(s1);
UNUSED(e0); UNUSED(e1); UNUSED(dV0); UNUSED(dV1); UNUSED(shift_phase);

source_slice_data *data = (source_slice_data *)data_;
ivec slice_imin = data->slice_imin, slice_imax=data->slice_imax;
ndim dim=fc->gv.dim;

// the following works in all cases except cylindrical coordinates
ptrdiff_t NY = 1, NZ = 1;
if (has_direction(fc->gv.dim, Z)) NZ = ((slice_imax - slice_imin).in_direction(Z) / 2) + 1;
if (has_direction(fc->gv.dim, Y)) NY = ((slice_imax - slice_imin).in_direction(Y) / 2) + 1;
if (has_direction(dim, Z)) NZ = ((slice_imax - slice_imin).in_direction(Z) / 2) + 1;
if (has_direction(dim, Y)) NY = ((slice_imax - slice_imin).in_direction(Y) / 2) + 1;

for (int ft = 0; ft < NUM_FIELD_TYPES; ft++)
for (src_vol *s = fc->sources[ft]; s; s = s->next) {
if (slice_component != s->c) continue;
component cS = S.transform(data->source_component, -sn);
if (s->c != cS) continue;

// loop over point sources in this src_vol. for each point source,
// the src_vol stores the amplitude and the global index of the grid point,
// from which we need to compute the local index of the grid point within the
// the src_vol stores the amplitude and the global index of the
// symmetry-parent grid point, from which we need to compute the
// 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];
ptrdiff_t chunk_index = s->index[npt];
ivec iloc = fc->gv.iloc(Dielectric, chunk_index);
if (iloc < slice_imin || iloc > slice_imax) continue; // source point outside slice
ivec slice_offset = iloc - slice_imin;
ivec iloc_parent = fc->gv.iloc(Dielectric, chunk_index);
ivec iloc_child = S.transform(iloc_parent,sn) + shift;
if (!in_subgrid(slice_imin, iloc_child, slice_imax)) continue; // source point outside slice
ivec slice_offset = iloc_child - slice_imin;
ptrdiff_t slice_index = 0;
// the following works to set the slice_index in all cases except cylindrical coordinates
if (has_direction(dim, Z)) slice_index += slice_offset.in_direction(Z) / 2;
if (has_direction(dim, Y)) slice_index += NZ * slice_offset.in_direction(Y) / 2;
if (has_direction(dim, X)) slice_index += NY * NZ * slice_offset.in_direction(X) / 2;

slice[slice_index] = amp;
data->slice[slice_index] = amp;
}
}
}
Expand All @@ -176,14 +209,6 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
UNUSED(dV1);
array_slice_data *data = (array_slice_data *)data_;

//-----------------------------------------------------------------------//
//- If we're fetching a 'source slice,' branch off here to handle that. //
//-----------------------------------------------------------------------//
if (data->get_source_slice) {
fill_chunk_source_slice(fc, data);
return;
}

//-----------------------------------------------------------------------//
// Find output chunk dimensions and strides, etc.

Expand Down Expand Up @@ -331,6 +356,30 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
} // LOOP_OVER_IVECS
}

/***************************************************************/
/* repeatedly call sum_to_all to consolidate all entries of */
/* an array on all cores. */
/***************************************************************/
#define BUFSIZE 1 << 16 // use 64k buffer
double *array_to_all(double *array, size_t array_size) {
double *buffer = new double[BUFSIZE];
ptrdiff_t offset = 0;
size_t remaining = array_size;
while (remaining != 0) {
size_t xfer_size = (remaining > BUFSIZE ? BUFSIZE : remaining);
sum_to_all(array + offset, buffer, xfer_size);
memcpy(array + offset, buffer, xfer_size*sizeof(double));
remaining -= xfer_size;
offset += xfer_size;
}
delete[] buffer;
return array;
}

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

/***************************************************************/
/* given a volume, fill in the dims[] and dirs[] arrays */
/* describing the array slice needed to store field data for */
Expand All @@ -344,7 +393,8 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
/* get_array_slice. */
/***************************************************************/
int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], direction dirs[3],
bool collapse_empty_dimensions, void *caller_data) {
bool collapse_empty_dimensions,
vec *min_max_loc, void *caller_data) {
am_now_working_on(FieldOutput);

// use a local data structure if the caller didn't provide one
Expand All @@ -356,10 +406,26 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], dire
data->max_corner = gv.round_vec(where.get_min_corner()) - one_ivec(gv.dim);
data->num_chunks = 0;

data->min_max_loc = min_max_loc;
vec *min_loc = 0, *max_loc = 0;
if (min_max_loc)
{ min_loc = min_max_loc + 0;
max_loc = min_max_loc + 1;
LOOP_OVER_DIRECTIONS(gv.dim,d) {
min_loc->set_direction(d,+infinity);
max_loc->set_direction(d,-infinity);
}
}

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
if (min_max_loc)
LOOP_OVER_DIRECTIONS(gv.dim,d) {
min_loc->set_direction(d, -1.0*max_to_all(-1.0*min_loc->in_direction(d)));
max_loc->set_direction(d, max_to_all( max_loc->in_direction(d)));
}
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;
Expand All @@ -386,12 +452,11 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], dire
}

/**********************************************************************/
/* precisely one of fun, rfun, source_slice should be non-NULL / true */
/* precisely one of fun, rfun should be non-NULL */
/**********************************************************************/
void *fields::do_get_array_slice(const volume &where, std::vector<component> components,
field_function fun, field_rfunction rfun, void *fun_data,
void *vslice, component source_slice_component,
bool get_source_slice) {
void *vslice) {
am_now_working_on(FieldOutput);

/***************************************************************/
Expand All @@ -404,7 +469,7 @@ void *fields::do_get_array_slice(const volume &where, std::vector<component> com
size_t dims[3];
direction dirs[3];
array_slice_data data;
int rank = get_array_slice_dimensions(where, dims, dirs, collapse, &data);
int rank = get_array_slice_dimensions(where, dims, dirs, collapse, 0, &data);
size_t slice_size = data.slice_size;
if (rank == 0 || slice_size == 0) return 0; // no data to write

Expand All @@ -427,8 +492,6 @@ void *fields::do_get_array_slice(const volume &where, std::vector<component> com
data.fun = fun;
data.rfun = rfun;
data.fun_data = fun_data;
data.source_slice_component = source_slice_component;
data.get_source_slice = get_source_slice;
data.components = components;
int num_components = components.size();
data.cS = new component[num_components];
Expand Down Expand Up @@ -470,37 +533,9 @@ void *fields::do_get_array_slice(const volume &where, std::vector<component> com
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 */
/*consolidate full array on all cores */
/***************************************************************/
#define BUFSIZE 1 << 16 // use 64k buffer
if (complex_data) {
cdouble *buffer = new cdouble[BUFSIZE];
cdouble *slice = (cdouble *)vslice;
ptrdiff_t offset = 0;
size_t remaining = slice_size;
while (remaining != 0) {
size_t 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;
ptrdiff_t offset = 0;
size_t remaining = slice_size;
while (remaining != 0) {
size_t 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;
}
vslice=(void *)array_to_all( (double *)vslice, (complex_data ? 2:1)*slice_size);

delete[] data.offsets;
delete[] data.fields;
Expand Down Expand Up @@ -546,9 +581,23 @@ cdouble *fields::get_complex_array_slice(const volume &where, component c, cdoub

cdouble *fields::get_source_slice(const volume &where, component source_slice_component,
cdouble *slice) {
vector<component> cs; // empty
return (cdouble *)do_get_array_slice(where, cs, 0, 0, 0, (void *)slice, source_slice_component,
true);

size_t dims[3];
direction dirs[3];
vec min_max_loc[2];
bool collapse=false;
int rank=get_array_slice_dimensions(where, dims, dirs, collapse, min_max_loc);
size_t slice_size = dims[0] * (rank>=2 ? dims[1] : 1) * (rank==3 ? dims[2] : 1);

source_slice_data data;
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 = slice ? slice : new cdouble[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, true);
return array_to_all(data.slice,slice_size);
}

} // namespace meep
6 changes: 3 additions & 3 deletions src/meep.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1517,7 +1517,8 @@ class fields {
// 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, size_t dims[3], direction dirs[3],
bool collapse_empty_dimensions = false, void *data = 0);
bool collapse_empty_dimensions = false,
vec *min_max_loc = NULL, void *data = 0);

int get_dft_array_dimensions(const volume &where, size_t dims[3], direction dirs[3]) {
return get_array_slice_dimensions(where, dims, dirs, true);
Expand Down Expand Up @@ -1551,8 +1552,7 @@ class fields {

// master routine for all above entry points
void *do_get_array_slice(const volume &where, std::vector<component> components,
field_function fun, field_rfunction rfun, void *fun_data, void *vslice,
component source_slice_component = Ex, bool get_source_slice = false);
field_function fun, field_rfunction rfun, void *fun_data, void *vslice);

// utility routine in loop_in_chunks.cpp to construct metadata for
// the arrays returned by get_array_slice and get_dft_array
Expand Down