diff options
Diffstat (limited to 'Source/Diagnostics/SliceDiagnostic.cpp')
| -rw-r--r-- | Source/Diagnostics/SliceDiagnostic.cpp | 252 |
1 files changed, 126 insertions, 126 deletions
diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp index 0ec528e32..79f03985b 100644 --- a/Source/Diagnostics/SliceDiagnostic.cpp +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -10,23 +10,23 @@ using namespace amrex; /* \brief * The functions creates the slice for diagnostics based on the user-input. - * The slice can be 1D, 2D, or 3D and it inherts the index type of the underlying data. - * The implementation assumes that the slice is aligned with the coordinate axes. - * The input parameters are modified if the user-input does not comply with requirements of coarsenability or if the slice extent is not contained within the simulation domain. - * First a slice multifab (smf) with cell size equal to that of the simulation grid is created such that it extends from slice.dim_lo to slice.dim_hi and shares the same index space as the source multifab (mf) + * The slice can be 1D, 2D, or 3D and it inherts the index type of the underlying data. + * The implementation assumes that the slice is aligned with the coordinate axes. + * The input parameters are modified if the user-input does not comply with requirements of coarsenability or if the slice extent is not contained within the simulation domain. + * First a slice multifab (smf) with cell size equal to that of the simulation grid is created such that it extends from slice.dim_lo to slice.dim_hi and shares the same index space as the source multifab (mf) * The values are copied from src mf to dst smf using amrex::ParallelCopy - * If interpolation is required, then on the smf, using data points stored in the ghost cells, the data in interpolated. - * If coarsening is required, then a coarse slice multifab is generated (cs_mf) and the - * values of the refined slice (smf) is averaged down to obtain the coarse slice. - * \param mf is the source multifab containing the field data - * \param dom_geom is the geometry of the domain and used in the function to obtain the - * CellSize of the underlying grid. - * \param slice_realbox defines the extent of the slice - * \param slice_cr_ratio provides the coarsening ratio for diagnostics + * If interpolation is required, then on the smf, using data points stored in the ghost cells, the data in interpolated. + * If coarsening is required, then a coarse slice multifab is generated (cs_mf) and the + * values of the refined slice (smf) is averaged down to obtain the coarse slice. + * \param mf is the source multifab containing the field data + * \param dom_geom is the geometry of the domain and used in the function to obtain the + * CellSize of the underlying grid. + * \param slice_realbox defines the extent of the slice + * \param slice_cr_ratio provides the coarsening ratio for diagnostics */ -std::unique_ptr<MultiFab> -CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, +std::unique_ptr<MultiFab> +CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, RealBox &slice_realbox, IntVect &slice_cr_ratio ) { std::unique_ptr<MultiFab> smf; @@ -36,10 +36,10 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, int nghost = 1; int nlevels = dom_geom.size(); int ncomp = (mf).nComp(); - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nlevels==1, + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nlevels==1, "Slice diagnostics does not work with mesh refinement yet (TO DO)."); - + const auto conversionType = (mf).ixType(); IntVect SliceType(AMREX_D_DECL(0,0,0)); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim ) @@ -50,7 +50,7 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, const RealBox& real_box = dom_geom[0].ProbDomain(); RealBox slice_cc_nd_box; int slice_grid_size = 32; - + bool interpolate = false; bool coarsen = false; @@ -60,45 +60,45 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, IntVect interp_lo(AMREX_D_DECL(0,0,0)); CheckSliceInput(real_box, slice_cc_nd_box, slice_realbox, slice_cr_ratio, - dom_geom, SliceType, slice_lo, + dom_geom, SliceType, slice_lo, slice_hi, interp_lo); int configuration_dim = 0; // Determine if interpolation is required and number of cells in slice // - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // Flag for interpolation if required // if ( interp_lo[idim] == 1) { - interpolate = 1; + interpolate = 1; } // For the case when a dimension is reduced // - if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { + if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { slice_ncells[idim] = 1; } - else { - slice_ncells[idim] = ( slice_hi[idim] - slice_lo[idim] + 1 ) + else { + slice_ncells[idim] = ( slice_hi[idim] - slice_lo[idim] + 1 ) / slice_cr_ratio[idim]; int refined_ncells = slice_hi[idim] - slice_lo[idim] + 1 ; - if ( slice_cr_ratio[idim] > 1) { + if ( slice_cr_ratio[idim] > 1) { coarsen = true; // modify slice_grid_size if >= refines_cells // if ( slice_grid_size >= refined_ncells ) { slice_grid_size = refined_ncells - 1; } - + } configuration_dim += 1; } } if (configuration_dim==1) { amrex::Warning("The slice configuration is 1D and cannot be visualized using yt."); - } + } // Slice generation with index type inheritance // Box slice(slice_lo, slice_hi); - + Vector<BoxArray> sba(1); sba[0].define(slice); sba[0].maxSize(slice_grid_size); @@ -106,7 +106,7 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, // Distribution mapping for slice can be different from that of domain // Vector<DistributionMapping> sdmap(1); sdmap[0] = DistributionMapping{sba[0]}; - + smf.reset(new MultiFab(amrex::convert(sba[0],SliceType), sdmap[0], ncomp, nghost)); @@ -115,12 +115,12 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, smf->ParallelCopy(mf, 0, 0, ncomp,nghost,nghost); // inteprolate if required on refined slice // - if (interpolate == 1 ) { - InterpolateSliceValues( *smf, interp_lo, slice_cc_nd_box, dom_geom, + if (interpolate == 1 ) { + InterpolateSliceValues( *smf, interp_lo, slice_cc_nd_box, dom_geom, ncomp, nghost, slice_lo, slice_hi, SliceType, real_box); } - + if (coarsen == false) { return smf; } @@ -131,14 +131,14 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, AMREX_ALWAYS_ASSERT(crse_ba[0].size() == sba[0].size()); - cs_mf.reset( new MultiFab(amrex::convert(crse_ba[0],SliceType), + cs_mf.reset( new MultiFab(amrex::convert(crse_ba[0],SliceType), sdmap[0], ncomp,nghost)); MultiFab& mfSrc = *smf; MultiFab& mfDst = *cs_mf; MFIter mfi_dst(mfDst); - for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { + for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { Array4<Real const> const& Src_fabox = mfSrc.const_array(mfi); @@ -150,7 +150,7 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, IntVect cctype(AMREX_D_DECL(0,0,0)); if( SliceType==cctype ) { - amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, + amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio); } IntVect ndtype(AMREX_D_DECL(1,1,1)); @@ -159,11 +159,11 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, scomp, ncomp, slice_cr_ratio); } if( SliceType == WarpX::Ex_nodal_flag ) { - amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 0); } if( SliceType == WarpX::Ey_nodal_flag) { - amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 1); } if( SliceType == WarpX::Ez_nodal_flag ) { @@ -171,19 +171,19 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, scomp, ncomp, slice_cr_ratio, 2); } if( SliceType == WarpX::Bx_nodal_flag) { - amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 0); } if( SliceType == WarpX::By_nodal_flag ) { - amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 1); } if( SliceType == WarpX::Bz_nodal_flag ) { - amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 2); } - if ( mfi_dst.isValid() ) { + if ( mfi_dst.isValid() ) { ++mfi_dst; } @@ -197,66 +197,66 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, /* \brief - * This function modifies the slice input parameters under certain conditions. - * The coarsening ratio, slice_cr_ratio is modified if the input is not an exponent of 2. - * for example, if the coarsening ratio is 3, 5 or 6, which is not an exponent of 2, - * then the value of coarsening ratio is modified to the nearest exponent of 2. - * The default value for coarsening ratio is 1. - * slice_realbox.lo and slice_realbox.hi are set equal to the simulation domain lo and hi - * if for the user-input for the slice lo and hi coordinates are outside the domain. - * If the slice_realbox.lo and slice_realbox.hi coordinates do not align with the data + * This function modifies the slice input parameters under certain conditions. + * The coarsening ratio, slice_cr_ratio is modified if the input is not an exponent of 2. + * for example, if the coarsening ratio is 3, 5 or 6, which is not an exponent of 2, + * then the value of coarsening ratio is modified to the nearest exponent of 2. + * The default value for coarsening ratio is 1. + * slice_realbox.lo and slice_realbox.hi are set equal to the simulation domain lo and hi + * if for the user-input for the slice lo and hi coordinates are outside the domain. + * If the slice_realbox.lo and slice_realbox.hi coordinates do not align with the data * points and the number of cells in that dimension is greater than 1, and if the extent of * the slice in that dimension is not coarsenable, then the value lo and hi coordinates are * shifted to the nearest coarsenable point to include some extra data points in the slice. * If slice_realbox.lo==slice_realbox.hi, then that dimension has only one cell and no - * modifications are made to the value. If the lo and hi do not align with a data point, - * then it is flagged for interpolation. - * \param real_box a Real box defined for the underlying domain. - * \param slice_realbox a Real box for defining the slice dimension. - * \param slice_cc_nd_box a Real box for defining the modified lo and hi of the slice + * modifications are made to the value. If the lo and hi do not align with a data point, + * then it is flagged for interpolation. + * \param real_box a Real box defined for the underlying domain. + * \param slice_realbox a Real box for defining the slice dimension. + * \param slice_cc_nd_box a Real box for defining the modified lo and hi of the slice * such that the coordinates align with the underlying data points. - * If the dimension is reduced to have only one cell, the slice_realbox is not modified and * instead the values are interpolated to the coordinate from the nearest data points. - * \param slice_cr_ratio contains values of the coarsening ratio which may be modified - * if the input values do not satisfy coarsenability conditions. + * If the dimension is reduced to have only one cell, the slice_realbox is not modified and * instead the values are interpolated to the coordinate from the nearest data points. + * \param slice_cr_ratio contains values of the coarsening ratio which may be modified + * if the input values do not satisfy coarsenability conditions. * \param slice_lo and slice_hi are the index values of the slice - * \param interp_lo are set to 0 or 1 if they are flagged for interpolation. - * The slice shares the same index space as that of the simulation domain. + * \param interp_lo are set to 0 or 1 if they are flagged for interpolation. + * The slice shares the same index space as that of the simulation domain. */ -void -CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, - RealBox &slice_realbox, IntVect &slice_cr_ratio, - Vector<Geometry> dom_geom, IntVect const SliceType, +void +CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, + RealBox &slice_realbox, IntVect &slice_cr_ratio, + Vector<Geometry> dom_geom, IntVect const SliceType, IntVect &slice_lo, IntVect &slice_hi, IntVect &interp_lo) { - + IntVect slice_lo2(AMREX_D_DECL(0,0,0)); - for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { // Modify coarsening ratio if the input value is not an exponent of 2 for AMR // - if ( slice_cr_ratio[idim] > 0 ) { + if ( slice_cr_ratio[idim] > 0 ) { int log_cr_ratio = floor ( log2( double(slice_cr_ratio[idim]))); slice_cr_ratio[idim] = exp2( double(log_cr_ratio) ); } - + //// Default coarsening ratio is 1 // // Modify lo if input is out of bounds // - if ( slice_realbox.lo(idim) < real_box.lo(idim) ) { + if ( slice_realbox.lo(idim) < real_box.lo(idim) ) { slice_realbox.setLo( idim, real_box.lo(idim)); - amrex::Print() << " slice lo is out of bounds. " << - " Modified it in dimension " << idim << + amrex::Print() << " slice lo is out of bounds. " << + " Modified it in dimension " << idim << " to be aligned with the domain box\n"; - } - + } + // Modify hi if input in out od bounds // - if ( slice_realbox.hi(idim) > real_box.hi(idim) ) { + if ( slice_realbox.hi(idim) > real_box.hi(idim) ) { slice_realbox.setHi( idim, real_box.hi(idim)); - amrex::Print() << " slice hi is out of bounds." << - " Modified it in dimension " << idim << + amrex::Print() << " slice hi is out of bounds." << + " Modified it in dimension " << idim << " to be aligned with the domain box\n"; } - + // Factor to ensure index values computation depending on index type // double fac = ( 1.0 - SliceType[idim] )*dom_geom[0].CellSize(idim)*0.5; // if dimension is reduced to one cell length // @@ -264,58 +264,58 @@ CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, { slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) ); slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) ); - + if ( slice_cr_ratio[idim] > 1) slice_cr_ratio[idim] = 1; - + // check for interpolation -- compute index lo with floor and ceil - if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) >= fac ) { - slice_lo[idim] = floor( ( (slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) + fac ) ) + if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) >= fac ) { + slice_lo[idim] = floor( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac ) ) / dom_geom[0].CellSize(idim)) + fac * 1E-10); - slice_lo2[idim] = ceil( ( (slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) + fac) ) - / dom_geom[0].CellSize(idim)) - fac * 1E-10 ); - } - else { - slice_lo[idim] = round( (slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) ) ) + slice_lo2[idim] = ceil( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac) ) + / dom_geom[0].CellSize(idim)) - fac * 1E-10 ); + } + else { + slice_lo[idim] = round( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) / dom_geom[0].CellSize(idim)); - slice_lo2[idim] = ceil((slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) ) ) + slice_lo2[idim] = ceil((slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) / dom_geom[0].CellSize(idim) ); } - + // flag for interpolation -- if reduced dimension location // // does not align with data point // - if ( slice_lo[idim] == slice_lo2[idim]) { + if ( slice_lo[idim] == slice_lo2[idim]) { if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) < fac ) { interp_lo[idim] = 1; } } - else { + else { interp_lo[idim] = 1; } - + // ncells = 1 if dimension is reduced // slice_hi[idim] = slice_lo[idim] + 1; } else { // moving realbox.lo and reabox.hi to nearest coarsenable grid point // - int index_lo = floor(((slice_realbox.lo(idim) + 1E-10 + int index_lo = floor(((slice_realbox.lo(idim) + 1E-10 - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); int index_hi = ceil(((slice_realbox.hi(idim) - 1E-10 - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); bool modify_cr = true; - - while ( modify_cr == true) { + + while ( modify_cr == true) { int lo_new = index_lo; - int hi_new = index_hi; + int hi_new = index_hi; int mod_lo = index_lo % slice_cr_ratio[idim]; int mod_hi = index_hi % slice_cr_ratio[idim]; modify_cr = false; - + // To ensure that the index.lo is coarsenable // if ( mod_lo > 0) { lo_new = index_lo - mod_lo; @@ -324,36 +324,36 @@ CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, if ( mod_hi > 0) { hi_new = index_hi + (slice_cr_ratio[idim] - mod_hi); } - - // If modified index.hi is > baselinebox.hi, move the point // + + // If modified index.hi is > baselinebox.hi, move the point // // to the previous coarsenable point // - if ( (hi_new * dom_geom[0].CellSize(idim)) + if ( (hi_new * dom_geom[0].CellSize(idim)) > real_box.hi(idim) - real_box.lo(idim) + dom_geom[0].CellSize(idim)*0.01 ) { hi_new = index_hi - mod_hi; } - - if ( (hi_new - lo_new) == 0 ){ + + if ( (hi_new - lo_new) == 0 ){ amrex::Print() << " Diagnostic Warning :: "; amrex::Print() << " Coarsening ratio "; amrex::Print() << slice_cr_ratio[idim] << " in dim "<< idim; amrex::Print() << "is leading to zero cells for slice."; amrex::Print() << " Thus reducing cr_ratio by half.\n"; - + slice_cr_ratio[idim] = slice_cr_ratio[idim]/2; modify_cr = true; } - - if ( modify_cr == false ) { + + if ( modify_cr == false ) { index_lo = lo_new; index_hi = hi_new; } slice_lo[idim] = index_lo; - slice_hi[idim] = index_hi - 1; // since default is cell-centered + slice_hi[idim] = index_hi - 1; // since default is cell-centered } - slice_realbox.setLo( idim, index_lo * dom_geom[0].CellSize(idim) + slice_realbox.setLo( idim, index_lo * dom_geom[0].CellSize(idim) + real_box.lo(idim) ); - slice_realbox.setHi( idim, index_hi * dom_geom[0].CellSize(idim) + slice_realbox.setHi( idim, index_hi * dom_geom[0].CellSize(idim) + real_box.lo(idim) ); slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) + fac ); slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) - fac ); @@ -363,14 +363,14 @@ CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, /* \brief - * This function is called if the coordinates of the slice do not align with data points - * \param interp_lo is an IntVect which is flagged as 1, if interpolation - is required in the dimension. + * This function is called if the coordinates of the slice do not align with data points + * \param interp_lo is an IntVect which is flagged as 1, if interpolation + is required in the dimension. */ -void +void InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, - Vector<Geometry> geom, int ncomp, int nghost, - IntVect slice_lo, IntVect slice_hi, IntVect SliceType, + Vector<Geometry> geom, int ncomp, int nghost, + IntVect slice_lo, IntVect slice_hi, IntVect SliceType, const RealBox real_box) { for (MFIter mfi(smf); mfi.isValid(); ++mfi) @@ -381,9 +381,9 @@ InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, const auto hi = amrex::ubound(bx); FArrayBox& fabox = smf[mfi]; - for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if ( interp_lo[idim] == 1 ) { - InterpolateLo( bx, fabox, slice_lo, geom, idim, SliceType, + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( interp_lo[idim] == 1 ) { + InterpolateLo( bx, fabox, slice_lo, geom, idim, SliceType, slice_realbox, 0, ncomp, nghost, real_box); } } @@ -391,10 +391,10 @@ InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, } -void -InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, - Vector<Geometry> geom, int idir, IntVect IndType, - RealBox slice_realbox, int srccomp, int ncomp, +void +InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, + Vector<Geometry> geom, int idir, IntVect IndType, + RealBox slice_realbox, int srccomp, int ncomp, int nghost, const RealBox real_box ) { auto fabarr = fabox.array(); @@ -428,7 +428,7 @@ InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, break; } case 1: - { + { if ( imin >= lo.y && imin <= lo.y) { for (int n = srccomp; n < srccomp+ncomp; ++n) { for (int k = lo.z; k <= hi.z; ++k) { |