Skip to content

Commit

Permalink
updates (#828)
Browse files Browse the repository at this point in the history
  • Loading branch information
@HomerReid @stevengj
HomerReid authored and stevengj committed Apr 19, 2019
1 parent a49c030 commit 5afef0c
Show file tree
Hide file tree
Showing 2 changed files with 91 additions and 85 deletions.
173 changes: 90 additions & 83 deletions src/array_slice.cpp
View file
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
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

0 comments on commit 5afef0c

Please sign in to comment.