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 #828

Merged
merged 1 commit into from
Apr 19, 2019
Merged

Improved source slice #828

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
173 changes: 90 additions & 83 deletions src/array_slice.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,17 +51,12 @@ typedef struct {
// in looping over the slice region.
vec *min_max_loc;

// if the component of the field array is NO_COMPONENT,
// the array is populated with the grid weights

// 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 @@ -133,54 +128,66 @@ static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, co
}
}

/*****************************************************************/
/* 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);
bool skip = false;
LOOP_OVER_DIRECTIONS(iloc.dim, d) {
if (iloc.in_direction(d) < slice_imin.in_direction(d) ||
iloc.in_direction(d) > slice_imax.in_direction(d)) {
skip = true;
break;
}
}
if (skip) 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 @@ -202,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 @@ -361,6 +360,31 @@ 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 Down Expand Up @@ -434,12 +458,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 Down Expand Up @@ -475,8 +498,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 @@ -518,37 +539,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 @@ -594,9 +587,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, snap=false;
int rank=get_array_slice_dimensions(where, dims, dirs, collapse, snap, 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);
}

/**********************************************************************/
Expand Down
3 changes: 1 addition & 2 deletions src/meep.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1555,8 +1555,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);

/* fetch and return coordinates and integration weights of grid points covered by an array slice,
*/
Expand Down