diff options
author | schnetter <> | 2004-03-31 10:50:00 +0000 |
---|---|---|
committer | schnetter <> | 2004-03-31 10:50:00 +0000 |
commit | dc5ef8bd7548e3cca7906524f69804d25c9ff478 (patch) | |
tree | 6f62166db0663e0aec4204578636b17a1740cad6 | |
parent | b6ff83da7b4afff959a15886b23504fbbb476601 (diff) |
Handle looping over refinement levels more elegantly.
darcs-hash:20040331105058-07bb3-0c0495c75cd248152699703e2e6d1520c9c80052.gz
-rw-r--r-- | Carpet/CarpetInterp/src/interp.cc | 309 |
1 files changed, 154 insertions, 155 deletions
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc index 4c74c67c0..73ee23fc1 100644 --- a/Carpet/CarpetInterp/src/interp.cc +++ b/Carpet/CarpetInterp/src/interp.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.24 2004/03/08 09:12:44 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.25 2004/03/31 12:50:58 schnetter Exp $ #include <assert.h> #include <math.h> @@ -21,7 +21,7 @@ #include "interp.hh" extern "C" { - static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.24 2004/03/08 09:12:44 schnetter Exp $"; + static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.25 2004/03/31 12:50:58 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc); } @@ -300,186 +300,185 @@ namespace CarpetInterp { // int overall_ierr = 0; BEGIN_GLOBAL_MODE(cgh) { - - BEGIN_REFLEVEL_LOOP(cgh) { - if (reflevel>=minrl && reflevel<maxrl) { - - // Number of necessary time levels - CCTK_REAL const time1 = vtt.at(m)->time (0, reflevel, mglevel); - CCTK_REAL const time2 = cgh->cctk_time / cgh->cctk_delta_time; - bool const need_time_interp - = fabs((time1 - time2) / (fabs(time1) + fabs(time2) - + fabs(cgh->cctk_delta_time))) > 1e-12; - int const num_tl - = need_time_interp ? prolongation_order_time + 1 : 1; + for (int rl=minrl; rl<maxrl; ++rl) { + enter_level_mode(cgh, rl); + + // Number of necessary time levels + CCTK_REAL const time1 = vtt.at(m)->time (0, reflevel, mglevel); + CCTK_REAL const time2 = cgh->cctk_time / cgh->cctk_delta_time; + bool const need_time_interp + = fabs((time1 - time2) / (fabs(time1) + fabs(time2) + + fabs(cgh->cctk_delta_time))) > 1e-12; + int const num_tl + = need_time_interp ? prolongation_order_time + 1 : 1; + + BEGIN_MAP_LOOP(cgh, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { + // Find out about the local geometry + ivect lsh; + rvect coord_origin, coord_delta; + for (int d=0; d<dim; ++d) { + lsh[d] = cgh->cctk_lsh[d]; + coord_delta[d] = cgh->cctk_delta_space[d] / cgh->cctk_levfac[d]; + coord_origin[d] = cgh->cctk_origin_space[d] + (cgh->cctk_lbnd[d] + 1.0 * cgh->cctk_levoff[d] / cgh->cctk_levoffdenom[d]) * coord_delta[d]; - // Find out about the local geometry - ivect lsh; - rvect coord_origin, coord_delta; - for (int d=0; d<dim; ++d) { - lsh[d] = cgh->cctk_lsh[d]; - coord_delta[d] = cgh->cctk_delta_space[d] / cgh->cctk_levfac[d]; - coord_origin[d] = cgh->cctk_origin_space[d] + (cgh->cctk_lbnd[d] + 1.0 * cgh->cctk_levoff[d] / cgh->cctk_levoffdenom[d]) * coord_delta[d]; + } + + + + // Find out about the grid functions + vector<CCTK_INT> input_array_type_codes(N_input_arrays); + vector<const void *> input_arrays(N_input_arrays); + for (int n=0; n<N_input_arrays; ++n) { + if (input_array_variable_indices[n] == -1) { - } + // Ignore this entry + input_array_type_codes.at(n) = -1; + } else { - - // Find out about the grid functions - vector<CCTK_INT> input_array_type_codes(N_input_arrays); - vector<const void *> input_arrays(N_input_arrays); - for (int n=0; n<N_input_arrays; ++n) { - if (input_array_variable_indices[n] == -1) { - - // Ignore this entry - input_array_type_codes.at(n) = -1; - - } else { - - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - - int const gi = CCTK_GroupIndexFromVarI (vi); - assert (gi>=0 && gi<CCTK_NumGroups()); - - cGroup group; - ierr = CCTK_GroupData (gi, &group); - assert (!ierr); - assert (group.grouptype == CCTK_GF); - assert (group.dim == dim); - assert (group.disttype == CCTK_DISTRIB_DEFAULT); - assert (group.stagtype == 0); // not staggered - - assert (group.numtimelevels >= num_tl); - - input_array_type_codes.at(n) = group.vartype; - - } - } // for input arrays - - - - // Work on the data from all processors - for (int p=0; p<nprocs; ++p) { - assert (allcoords.at(ind_prc(p,m,reflevel,component)).owns_storage()); - assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == allcoords.at(ind_prc(p,m,reflevel,component)).shape()[0]); - assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == alloutputs.at(ind_prc(p,m,reflevel,component)).shape()[0]); + int const vi = input_array_variable_indices[n]; + assert (vi>=0 && vi<CCTK_NumVars()); - int const npoints = allhomecnts.at(ind_prc(p,m,reflevel,component)); + int const gi = CCTK_GroupIndexFromVarI (vi); + assert (gi>=0 && gi<CCTK_NumGroups()); - // Do the processor-local interpolation - vector<const void *> tmp_interp_coords (dim); - for (int d=0; d<dim; ++d) { - tmp_interp_coords.at(d) = &allcoords.at(ind_prc(p,m,reflevel,component))[ivect(0,d,0)]; - } + cGroup group; + ierr = CCTK_GroupData (gi, &group); + assert (!ierr); + assert (group.grouptype == CCTK_GF); + assert (group.dim == dim); + assert (group.disttype == CCTK_DISTRIB_DEFAULT); + assert (group.stagtype == 0); // not staggered + assert (group.numtimelevels >= num_tl); + input_array_type_codes.at(n) = group.vartype; + + } + } // for input arrays + + + + // Work on the data from all processors + for (int p=0; p<nprocs; ++p) { + assert (allcoords.at(ind_prc(p,m,reflevel,component)).owns_storage()); + assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == allcoords.at(ind_prc(p,m,reflevel,component)).shape()[0]); + assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == alloutputs.at(ind_prc(p,m,reflevel,component)).shape()[0]); + + int const npoints = allhomecnts.at(ind_prc(p,m,reflevel,component)); + + // Do the processor-local interpolation + vector<const void *> tmp_interp_coords (dim); + for (int d=0; d<dim; ++d) { + tmp_interp_coords.at(d) = &allcoords.at(ind_prc(p,m,reflevel,component))[ivect(0,d,0)]; + } + + + + vector<vector<void *> > tmp_output_arrays (num_tl); + + for (int tl=0; tl<num_tl; ++tl) { - vector<vector<void *> > tmp_output_arrays (num_tl); + for (int n=0; n<N_input_arrays; ++n) { + if (input_array_variable_indices[n] == -1) { + input_arrays.at(n) = 0; + } else { + int const vi = input_array_variable_indices[n]; + assert (vi>=0 && vi<CCTK_NumVars()); + input_arrays.at(n) = CCTK_VarDataPtrI (cgh, tl, vi); + } + } // for input arrays - for (int tl=0; tl<num_tl; ++tl) { - - for (int n=0; n<N_input_arrays; ++n) { - if (input_array_variable_indices[n] == -1) { - input_arrays.at(n) = 0; - } else { - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - input_arrays.at(n) = CCTK_VarDataPtrI (cgh, tl, vi); - } - } // for input arrays - - tmp_output_arrays.at(tl).resize (N_output_arrays); - for (int j=0; j<N_output_arrays; ++j) { - assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL); - if (need_time_interp) { - tmp_output_arrays.at(tl).at(j) = new CCTK_REAL [npoints]; - } else { - tmp_output_arrays.at(tl).at(j) = &alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(0,j,0)]; - } + tmp_output_arrays.at(tl).resize (N_output_arrays); + for (int j=0; j<N_output_arrays; ++j) { + assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL); + if (need_time_interp) { + tmp_output_arrays.at(tl).at(j) = new CCTK_REAL [npoints]; + } else { + tmp_output_arrays.at(tl).at(j) = &alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(0,j,0)]; } + } - ierr = CCTK_InterpLocalUniform - (N_dims, local_interp_handle, param_table_handle, - &coord_origin[0], &coord_delta[0], - npoints, - interp_coords_type_code, &tmp_interp_coords[0], - N_input_arrays, &lsh[0], - &input_array_type_codes.at(0), &input_arrays[0], - N_output_arrays, - output_array_type_codes, &tmp_output_arrays.at(tl)[0]); - if (ierr) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "The local interpolator returned the error code %d", ierr); - } - overall_ierr = min(overall_ierr, ierr); + ierr = CCTK_InterpLocalUniform + (N_dims, local_interp_handle, param_table_handle, + &coord_origin[0], &coord_delta[0], + npoints, + interp_coords_type_code, &tmp_interp_coords[0], + N_input_arrays, &lsh[0], + &input_array_type_codes.at(0), &input_arrays[0], + N_output_arrays, + output_array_type_codes, &tmp_output_arrays.at(tl)[0]); + if (ierr) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "The local interpolator returned the error code %d", ierr); + } + overall_ierr = min(overall_ierr, ierr); - } // for tl + } // for tl // Interpolate in time, if necessary - if (need_time_interp) { + if (need_time_interp) { - // Get interpolation times - CCTK_REAL const time = cgh->cctk_time / cgh->cctk_delta_time; - vector<CCTK_REAL> times(num_tl); - for (int tl=0; tl<num_tl; ++tl) { - times.at(tl) = vtt.at(m)->time (tl, reflevel, mglevel); - } + // Get interpolation times + CCTK_REAL const time = cgh->cctk_time / cgh->cctk_delta_time; + vector<CCTK_REAL> times(num_tl); + for (int tl=0; tl<num_tl; ++tl) { + times.at(tl) = vtt.at(m)->time (tl, reflevel, mglevel); + } - // Calculate interpolation weights - vector<CCTK_REAL> tfacs(num_tl); - switch (num_tl) { - case 1: - // no interpolation - assert (fabs((time - times.at(0)) / fabs(time + times.at(0) + cgh->cctk_delta_time)) < 1e-12); - tfacs.at(0) = 1.0; - break; - case 2: - // linear (2-point) interpolation - tfacs.at(0) = (time - times.at(1)) / (times.at(0) - times.at(1)); - tfacs.at(1) = (time - times.at(0)) / (times.at(1) - times.at(0)); - break; - case 3: - // quadratic (3-point) interpolation - tfacs.at(0) = (time - times.at(1)) * (time - times.at(2)) / ((times.at(0) - times.at(1)) * (times.at(0) - times.at(2))); - tfacs.at(1) = (time - times.at(0)) * (time - times.at(2)) / ((times.at(1) - times.at(0)) * (times.at(1) - times.at(2))); - tfacs.at(2) = (time - times.at(0)) * (time - times.at(1)) / ((times.at(2) - times.at(0)) * (times.at(2) - times.at(1))); - break; - default: - assert (0); - } + // Calculate interpolation weights + vector<CCTK_REAL> tfacs(num_tl); + switch (num_tl) { + case 1: + // no interpolation + assert (fabs((time - times.at(0)) / fabs(time + times.at(0) + cgh->cctk_delta_time)) < 1e-12); + tfacs.at(0) = 1.0; + break; + case 2: + // linear (2-point) interpolation + tfacs.at(0) = (time - times.at(1)) / (times.at(0) - times.at(1)); + tfacs.at(1) = (time - times.at(0)) / (times.at(1) - times.at(0)); + break; + case 3: + // quadratic (3-point) interpolation + tfacs.at(0) = (time - times.at(1)) * (time - times.at(2)) / ((times.at(0) - times.at(1)) * (times.at(0) - times.at(2))); + tfacs.at(1) = (time - times.at(0)) * (time - times.at(2)) / ((times.at(1) - times.at(0)) * (times.at(1) - times.at(2))); + tfacs.at(2) = (time - times.at(0)) * (time - times.at(1)) / ((times.at(2) - times.at(0)) * (times.at(2) - times.at(1))); + break; + default: + assert (0); + } - // Interpolate - for (int j=0; j<N_output_arrays; ++j) { - assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL); - for (int k=0; k<npoints; ++k) { - CCTK_REAL & dest = alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(k,j,0)]; - dest = 0; - for (int tl=0; tl<num_tl; ++tl) { - CCTK_REAL const src = ((CCTK_REAL const *)tmp_output_arrays.at(tl).at(j))[k]; - dest += tfacs[tl] * src; - } + // Interpolate + for (int j=0; j<N_output_arrays; ++j) { + assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL); + for (int k=0; k<npoints; ++k) { + CCTK_REAL & dest = alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(k,j,0)]; + dest = 0; + for (int tl=0; tl<num_tl; ++tl) { + CCTK_REAL const src = ((CCTK_REAL const *)tmp_output_arrays.at(tl).at(j))[k]; + dest += tfacs[tl] * src; } } + } - for (int tl=0; tl<num_tl; ++tl) { - for (int j=0; j<N_output_arrays; ++j) { - delete [] (CCTK_REAL *)tmp_output_arrays.at(tl).at(j); - } + for (int tl=0; tl<num_tl; ++tl) { + for (int j=0; j<N_output_arrays; ++j) { + delete [] (CCTK_REAL *)tmp_output_arrays.at(tl).at(j); } + } - } // if need_time_interp + } // if need_time_interp - } // for processors + } // for processors - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; - } // if reflevel active - } END_REFLEVEL_LOOP; + leave_level_mode (cgh); + } // for rl } END_GLOBAL_MODE; |