From ca387a997fbd4624f0e5928c3192b6024c71dcd8 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 11 Aug 2023 16:22:33 -0700 Subject: [PATCH] fgradient: A new plotfile tool for computing gradient --- .../source/Post_Processing.rst | 19 + Src/Base/AMReX_PhysBCFunct.H | 9 + Tools/Plotfile/CMakeLists.txt | 1 + Tools/Plotfile/GNUmakefile | 3 + Tools/Plotfile/fgradient.cpp | 344 ++++++++++++++++++ 5 files changed, 376 insertions(+) create mode 100644 Tools/Plotfile/fgradient.cpp diff --git a/Docs/sphinx_documentation/source/Post_Processing.rst b/Docs/sphinx_documentation/source/Post_Processing.rst index c18d00359e7..855dfda0de2 100644 --- a/Docs/sphinx_documentation/source/Post_Processing.rst +++ b/Docs/sphinx_documentation/source/Post_Processing.rst @@ -316,3 +316,22 @@ Typing ``./faverage.gnu.ex`` without inputs will bring up usage and options. will compute the average density as a function of height, outputting a data file ``plt0000000.slice``. + +fgradient +--------- + +Compute the gradient of variables in a plotfile. + +**How to build and run** + +In ``amrex/Tools/Plotfile``, type ``make programs=fgradient`` and then ``./fgradient.gnu.ex`` to run. +Typing ``./fgradient.gnu.ex`` without inputs will bring up usage and options. + +**Example** + +.. code-block:: console + + user@:~/AMReX/amrex/Tools/Plotfile$ ./fgradient.gnu.ex -v density plt0000000 + +will compute the gradient of density as a function of height, outputting a +new plotfile ``grad.plt0000000``. diff --git a/Src/Base/AMReX_PhysBCFunct.H b/Src/Base/AMReX_PhysBCFunct.H index 0589d30a217..2d0906f27ad 100644 --- a/Src/Base/AMReX_PhysBCFunct.H +++ b/Src/Base/AMReX_PhysBCFunct.H @@ -90,6 +90,15 @@ protected: F m_user_f; }; +struct FabFillNoOp +{ + AMREX_GPU_DEVICE + void operator() (const amrex::IntVect&, amrex::Array4 const&, + int, int, amrex::GeometryData const&, amrex::Real, + const amrex::BCRec*, int, int) const + {} +}; + //! This cpu version calls function working on FArrayBox class CpuBndryFuncFab { diff --git a/Tools/Plotfile/CMakeLists.txt b/Tools/Plotfile/CMakeLists.txt index 9f8f066fbbb..0efc55c9a16 100644 --- a/Tools/Plotfile/CMakeLists.txt +++ b/Tools/Plotfile/CMakeLists.txt @@ -11,6 +11,7 @@ set(_exe_names fcompare fextract fextrema + fgradient fnan fsnapshot ftime diff --git a/Tools/Plotfile/GNUmakefile b/Tools/Plotfile/GNUmakefile index 336dc59f212..afd35e688f5 100644 --- a/Tools/Plotfile/GNUmakefile +++ b/Tools/Plotfile/GNUmakefile @@ -17,6 +17,7 @@ ifeq ($(strip $(programs)),) programs += fcompare programs += fextract programs += fextrema + programs += fgradient programs += fnan programs += fsnapshot programs += ftime @@ -31,6 +32,8 @@ default: $(multiple_executables) include ./Make.package include $(AMREX_HOME)/Src/Base/Make.package +include $(AMREX_HOME)/Src/Boundary/Make.package +include $(AMREX_HOME)/Src/AmrCore/Make.package include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Tools/Plotfile/fgradient.cpp b/Tools/Plotfile/fgradient.cpp new file mode 100644 index 00000000000..9347f7595e2 --- /dev/null +++ b/Tools/Plotfile/fgradient.cpp @@ -0,0 +1,344 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace amrex; + +void PrintUsage () +{ + amrex::Print() + << "\n" + << " Compute gradient of variables in a given plofile. The results are saved in a\n" + << " new plotfile.\n" + << "\n" + << " Usage:\n" + << " fgradient [-o outfile] [-v varnames] [-xlo xlo_bc] [-ylo ylo_bc]\n" + << " [-zlo zlo_bc] [-xhi xhi_bc] [-yhi yhi_bc] [-zhi zhi_bc]\n" + << " [-i interp] plotfile\n" + << "\n" + << " optional arguments:\n" + << " -o|--output outfile\n" + << " Gradient plotfile name. If it's not provided, it is grad.PLOTFILE,\n" + << " where PLOTFILE is the original plotfile name.\n" + << "\n" + << " -v|--variables varnames\n" + << " Variable names separated by spaces. If there are multiple names,\n" + << " they must be inside a pair of double quotes. If not provided, all\n" + << " variables in the plotfile are included.\n" + << "\n" + << " -xlo|--xlo_boundary xlo_bc\n" + << " Domain boundary conditions at x-lo. Choices are\n" + << " reflect_odd : reflecting odd\n" + << " reflect_even : reflecting even\n" + << " foextrap : first-order extrapolation\n" + << " hoextrap : second-order extrapolation.\n" + << " periodic : periodic boundary.\n" + << " One could provide either different BC types for different\n" + << " variables or just one type for all variables. If not provided, the\n" + << " default is hoextrap, equivalent to one-sided derivative.\n" + << "\n" + << " -ylo|--ylo_boundary ylo_bc\n" + << " Domain boundary conditions at y-lo.\n" + << "\n" + << " -zlo|--zlo_boundary zlo_bc\n" + << " Domain boundary conditions at z-lo.\n" + << "\n" + << " -xhi|--xhi_boundary xhi_bc\n" + << " Domain boundary conditions at x-hi.\n" + << "\n" + << " -yhi|--yhi_boundary yhi_bc\n" + << " Domain boundary conditions at y-hi.\n" + << "\n" + << " -zhi|--zhi_boundary zhi_bc\n" + << " Domain boundary conditions at z-hi.\n" + << "\n" + << " -i|--interpolater interp\n" + << " Interpolation from coarse to fine level. Choices are 0 (default)\n" + << " or 1, where 0 is amrex::cell_cons_interp (conservative linear)\n" + << " and 1 is amrex::pc_interp (piecewise constant).\n" + << "\n" + << " -h|--help\n" + << " Print usage.\n" + << "\n"; +} + +void main_main() +{ + const int narg = amrex::command_argument_count(); + + std::string pltfile; + std::string outfile; + std::string varnames_arg; + std::string xlobc_arg; + std::string xhibc_arg; + std::string ylobc_arg; + std::string yhibc_arg; + std::string zlobc_arg; + std::string zhibc_arg; + int interpolater = 0; + bool print_help = false; + + int farg = 1; + while (farg <= narg) { + const std::string& name = amrex::get_command_argument(farg); + if (name == "-o" || name == "--output") { + outfile = amrex::get_command_argument(++farg); + } else if (name == "-v" || name == "--variables") { + varnames_arg = amrex::get_command_argument(++farg); + } else if (name == "-xlo" || name == "--xlo_boundary") { + xlobc_arg = amrex::get_command_argument(++farg); + } else if (name == "-xhi" || name == "--xhi_boundary") { + xhibc_arg = amrex::get_command_argument(++farg); + } else if (name == "-ylo" || name == "--ylo_boundary") { + ylobc_arg = amrex::get_command_argument(++farg); + } else if (name == "-yhi" || name == "--yhi_boundary") { + yhibc_arg = amrex::get_command_argument(++farg); + } else if (name == "-zlo" || name == "--zlo_boundary") { + zlobc_arg = amrex::get_command_argument(++farg); + } else if (name == "-zhi" || name == "--zhi_boundary") { + zhibc_arg = amrex::get_command_argument(++farg); + } else if (name == "-i" || name == "--interpolater") { + interpolater = std::stoi(amrex::get_command_argument(++farg)); + } else if (name == "-h" || name == "--help") { + print_help = true; + } else { + break; + } + ++farg; + } + + if (farg <= narg) { + pltfile = amrex::get_command_argument(farg); + } + + if (pltfile.empty() || print_help) { + PrintUsage(); + return; + } + + if (pltfile.back() == '/') { + pltfile.pop_back(); + } + + if (outfile.empty()) { + outfile = "grad."+std::filesystem::path(pltfile).filename().string(); + } + + PlotFileData pf(pltfile); + + const int ndims = pf.spaceDim(); + AMREX_ALWAYS_ASSERT(ndims <= AMREX_SPACEDIM); + + const int nlevs = pf.finestLevel() + 1; + + Vector varnames; + if (varnames_arg.empty()) { + varnames = pf.varNames(); + } else { + auto const& all_names = pf.varNames(); + std::istringstream is(varnames_arg); + Vector args(std::istream_iterator{is}, + std::istream_iterator{ }); + for (auto const& vn : args) { + if (std::find(std::begin(all_names),std::end(all_names), vn) + != std::end(all_names)) { + varnames.push_back(vn); + } else { + amrex::Abort("fgradient: Unknown variable "+vn); + } + } + } + + Vector gvarnames; + for (auto const& name : varnames) { + gvarnames.push_back(name+"_x"); +#if (AMREX_SPACEDIM >= 2) + if (ndims > 1) { + gvarnames.push_back(name+"_y"); + } +#endif +#if (AMREX_SPACEDIM == 3) + if (ndims > 2) { + gvarnames.push_back(name+"_z"); + } +#endif + } + + const auto nvars = int(varnames.size()); + + BCRec bcr_default; + Array is_periodic{AMREX_D_DECL(0,0,0)}; + IntVect ng(1); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (idim < ndims) { + bcr_default.setLo(idim, BCType::hoextrapcc); + bcr_default.setHi(idim, BCType::hoextrapcc); + } else { + bcr_default.setLo(idim, BCType::int_dir); + bcr_default.setHi(idim, BCType::int_dir); + is_periodic[idim] = 1; + ng[idim] = 0; + } + } + + auto f_bcarg = [] (std::string const& bcarg) -> Vector + { + Vector r; + std::istringstream is(bcarg); + Vector args(std::istream_iterator{is}, + std::istream_iterator{ }); + for (auto const& bc : args) { + if (bc == "reflect_odd") { + r.push_back(BCType::reflect_odd); + } else if (bc == "reflect_even") { + r.push_back(BCType::reflect_even); + } else if (bc == "foextrap") { + r.push_back(BCType::foextrap); + } else if (bc == "hoextrap") { + r.push_back(BCType::hoextrapcc); + } else if (bc == "periodic") { + r.push_back(BCType::int_dir); + } else { + amrex::Abort("fgradient: Unknown BC type "+bc); + } + } + return r; + }; + + Array,3> lobc{f_bcarg(xlobc_arg), + f_bcarg(ylobc_arg), + f_bcarg(zlobc_arg)}; + Array,3> hibc{f_bcarg(xhibc_arg), + f_bcarg(yhibc_arg), + f_bcarg(zhibc_arg)}; + + Vector gmf(nlevs); + Vector geom; + for (int ilev = 0; ilev < nlevs; ++ilev) + { + gmf[ilev].define(pf.boxArray(ilev), pf.DistributionMap(ilev), nvars*ndims, 0); + + for (int ivar = 0; ivar < nvars; ++ivar) + { + Vector bcr{bcr_default}; + auto is_per = is_periodic; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + int bc_type = BCType::hoextrapcc; + auto nbcs = int(lobc[idim].size()); + if (nbcs == 1) { + bc_type = lobc[idim][0]; + } else if (ivar < nbcs) { + bc_type = lobc[idim][ivar]; + } + if (bc_type == BCType::int_dir) { + is_per[idim] = 1; + } + bcr[0].setLo(idim, bc_type); + // + bc_type = BCType::hoextrapcc; + nbcs = int(hibc[idim].size()); + if (nbcs == 1) { + bc_type = hibc[idim][0]; + } else if (ivar < nbcs) { + bc_type = hibc[idim][ivar]; + } + if (bc_type == BCType::int_dir) { + AMREX_ALWAYS_ASSERT(is_per[idim] == 1); + } + bcr[0].setHi(idim, bc_type); + } + + Geometry vargeom(pf.probDomain(ilev), RealBox(pf.probLo(),pf.probHi()), + pf.coordSys(), is_per); + if (ivar == 0) { + geom.push_back(vargeom); + } + + MultiFab mf(pf.boxArray(ilev), pf.DistributionMap(ilev), 1, ng); + if (ilev == 0) { + MultiFab smf = pf.get(ilev, varnames[ivar]); + PhysBCFunct> physbcf + (vargeom, bcr, GpuBndryFuncFab(FabFillNoOp{})); + + FillPatchSingleLevel(mf, ng, Real(0.0), {&smf}, {Real(0.0)}, + 0, 0, 1, vargeom, physbcf, 0); + } else { + MultiFab cmf = pf.get(ilev-1, varnames[ivar]); + MultiFab fmf = pf.get(ilev , varnames[ivar]); + Geometry cgeom(pf.probDomain(ilev-1), RealBox(pf.probLo(),pf.probHi()), + pf.coordSys(), is_per); + PhysBCFunct> cphysbcf + (cgeom, bcr, GpuBndryFuncFab(FabFillNoOp{})); + PhysBCFunct> fphysbcf + (vargeom, bcr, GpuBndryFuncFab(FabFillNoOp{})); + auto* mapper = (interpolater == 0) + ? (Interpolater*)(&cell_cons_interp) + : (Interpolater*)(&pc_interp); + + IntVect ratio(pf.refRatio(ilev-1)); + for (int idim = ndims; idim < AMREX_SPACEDIM; ++idim) { + ratio[idim] = 1; + } + FillPatchTwoLevels(mf, ng, Real(0.0), {&cmf}, {Real(0.0)}, + {&fmf}, {Real(0.0)}, 0, 0, 1, cgeom, vargeom, + cphysbcf, 0, fphysbcf, 0, ratio, mapper, bcr, 0); + } + + auto const& dx = pf.cellSize(ilev); + +#ifdef AMREX_USE_OMP +#pragma omp parallel +#endif + for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + auto const& ga = gmf[ilev].array(mfi,ivar*ndims); + auto const& a = mf.const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + ga(i,j,k,0) = (a(i+1,j,k) - a(i-1,j,k)) / (Real(2.0)*dx[0]); +#if (AMREX_SPACEDIM >= 2) + if (ndims > 1) { + ga(i,j,k,1) = (a(i,j+1,k) - a(i,j-1,k)) / (Real(2.0)*dx[1]); + } +#endif +#if (AMREX_SPACEDIM == 3) + if (ndims > 2) { + ga(i,j,k,2) = (a(i,j,k+1) - a(i,j,k-1)) / (Real(2.0)*dx[2]); + } +#endif + }); + } + } + } + + Vector level_steps; + Vector ref_ratio; + for (int ilev = 0; ilev < nlevs; ++ilev) { + level_steps.push_back(pf.levelStep(ilev)); + if (ilev < pf.finestLevel()) { + ref_ratio.push_back(IntVect(pf.refRatio(ilev))); + for (int idim = ndims; idim < AMREX_SPACEDIM; ++idim) { + ref_ratio[ilev][idim] = 1; + } + } + } + + WriteMultiLevelPlotfile(outfile, nlevs, GetVecOfConstPtrs(gmf), gvarnames, + geom, pf.time(), level_steps, ref_ratio); +} + +int main (int argc, char* argv[]) +{ + amrex::SetVerbose(0); + amrex::Initialize(argc, argv, false); + main_main(); + amrex::Finalize(); +}