diff --git a/.Exec_dev/DataAssimilation/ERF_Prob.H b/.Exec_dev/DataAssimilation/ERF_Prob.H index 1203051b7..8ff4b6126 100644 --- a/.Exec_dev/DataAssimilation/ERF_Prob.H +++ b/.Exec_dev/DataAssimilation/ERF_Prob.H @@ -4,7 +4,6 @@ #include #include "AMReX_REAL.H" - #include "ERF_ProbCommon.H" struct ProbParm : ProbParmDefaults { @@ -29,6 +28,8 @@ class Problem : public ProblemBase public: Problem(const amrex_real* problo, const amrex_real* probhi); +#include "Prob/ERF_InitDensityHSE.H" +#include "Prob/ERF_InitRayleighDamping.H" #include "Prob/ERF_InitCustomPert.H" protected: diff --git a/.Exec_dev/DataAssimilation/ERF_Prob.cpp b/.Exec_dev/DataAssimilation/ERF_Prob.cpp index 814034681..32bf66dbc 100644 --- a/.Exec_dev/DataAssimilation/ERF_Prob.cpp +++ b/.Exec_dev/DataAssimilation/ERF_Prob.cpp @@ -129,13 +129,14 @@ Problem::init_custom_pert_vels ( Array4 const& x_vel_pert, Array4 const& y_vel_pert, Array4 const& z_vel_pert, - Array4 const& r_hse, + amrex::Array4 const& z_nd, GeometryData const& geomdata, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, const SolverChoice& /*sc*/, const int /*lev*/) { + ignore_unused(z_nd); // -------------------------------------------------------- // Per-ensemble perturbation controls // -------------------------------------------------------- @@ -143,6 +144,9 @@ Problem::init_custom_pert_vels ( ParmParse pp_ens("ensemble_pert"); pp_ens.query("amplitude", ens_pert_amplitude); + Real xc = parms.xc; Real yc = parms.yc; + Real R = parms.R ; Real beta = parms.beta; + Real sigma = parms.sigma; // Set the x-velocity ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/.Exec_dev/DataAssimilation/GNUmakefile b/.Exec_dev/DataAssimilation/GNUmakefile index 56fa6cabf..640374a14 100644 --- a/.Exec_dev/DataAssimilation/GNUmakefile +++ b/.Exec_dev/DataAssimilation/GNUmakefile @@ -28,6 +28,6 @@ USE_ASSERTION = TRUE Bpack := ./Make.package Blocs := . -ERF_HOME := ../../.. -ERF_PROBLEM_DIR = $(ERF_HOME)/.Exec_dev/Assimilation +ERF_HOME := ../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/.Exec_dev/DataAssimilation include $(ERF_HOME)/Exec/Make.ERF diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/Makefile b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/Makefile new file mode 100644 index 000000000..20b24cd94 --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/Makefile @@ -0,0 +1,58 @@ + + +#--------------------------- +# Makefile for AMReX project +#--------------------------- +CXX = mpicxx +CXXFLAGS = -O2 -std=c++17 -fopenmp -Wall + +# Replace with your AMReX paths +AMREX_INC = /pscratch/sd/n/nataraj2/AIModCon/GitHub/amrex_repos/amrex/include +AMREX_LIB = /pscratch/sd/n/nataraj2/AIModCon/GitHub/amrex_repos/amrex/lib + +HDRS = ReadPlotFile.H +LIBS = -lamrex -lpthread -lgfortran -ldl -lm -lstdc++ + +# Default: do nothing if 'make' is called without target +all: + @echo "Use 'make serial' or 'make parallel'" + +# ------------------------- +# Serial build +# ------------------------- +SRCS_SERIAL = main_serial.cpp ReadPlotFile.cpp +OBJS_SERIAL = $(SRCS_SERIAL:.cpp=.o) +EXE_SERIAL = out_serial + +serial: $(EXE_SERIAL) + +$(EXE_SERIAL): $(OBJS_SERIAL) + $(CXX) $(CXXFLAGS) -o $@ $^ -L$(AMREX_LIB) $(LIBS) + +# ------------------------- +# Parallel build +# ------------------------- +SRCS_PARALLEL = main_parallel.cpp ReadPlotFile.cpp +OBJS_PARALLEL = $(SRCS_PARALLEL:.cpp=.o) +EXE_PARALLEL = out_parallel + +parallel: $(EXE_PARALLEL) + +$(EXE_PARALLEL): $(OBJS_PARALLEL) + $(CXX) $(CXXFLAGS) -o $@ $^ -L$(AMREX_LIB) $(LIBS) + +# ------------------------- +# Compile rules for object files +# ------------------------- +%.o: %.cpp $(HDRS) + $(CXX) $(CXXFLAGS) -I$(AMREX_INC) -c $< -o $@ + +# ------------------------- +# Clean +# ------------------------- +clean: + rm -f *.o out_serial out_parallel + +.PHONY: all serial parallel clean + + diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H new file mode 100644 index 000000000..861203e4b --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -0,0 +1,60 @@ +#ifndef READ_PLOTFILE_H_ +#define READ_PLOTFILE_H_ + +#include +#include +#include +#include +#include + +amrex::Vector ReadVarNames(const std::string& filename); + +// ------------------------------------------------------------ +// ReadPlotFile +// +// var_filename : text file containing variable names (one per line) +// plot_filename : AMReX plotfile directory +// mf : MultiFab that will be defined and filled +// ------------------------------------------------------------ +void ReadPlotFile(const std::string& var_filename, + amrex::PlotFileData& pf, + amrex::MultiFab& mf); + +void +GetCoarseMultiFabOnFineDMap(const amrex::Geometry& geom_coarse, + const amrex::Geometry& geom_fine, + const amrex::MultiFab& multifab_coarse, + const amrex::MultiFab& multifab_fine, + amrex::MultiFab& coarse_multifab_on_fine_dmap); + +void +CreateNodalMultiFabFromCellCenteredMultiFab (amrex::MultiFab& mf_nc, // output nodal MF + amrex::MultiFab& mf_cc, // input cell-centered MF (coarse) + const amrex::Geometry& geom); + +void +CreateCellCenteredMultiFabFromNodalMultiFab(amrex::MultiFab& tmp_cc, // cell-centered output + amrex::MultiFab& mf_nc); // node-centered input + +// Populate a fine, cell-centered MultiFab from coarse nodal data using +// trilinear interpolation at fine cell centers. +// mf_cc_fine provides the fine BoxArray + DistributionMapping template. +void +PopulateFineCellCenteredFromCoarseNodal(const amrex::Geometry& geom_coarse, + const amrex::Geometry& geom_fine, + const amrex::MultiFab& coarse_multifab_on_fine_dmap, + const amrex::MultiFab& mf_cc_fine, + amrex::MultiFab& mf_cc_from_coarse); + +void +check_large_values(const amrex::MultiFab& mf); + +void +ApplyNeumannBCs(const amrex::Geometry& geom, + amrex::MultiFab& mf_cc); + +void WriteCustomDataFile(const amrex::Geometry& geom, + const amrex::MultiFab& mf_cc, + const std::string& filename_custom); + +#endif diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp new file mode 100644 index 000000000..e8d6e2101 --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -0,0 +1,554 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "ReadPlotFile.H" + +#include +#include // for std::remove_if, std::isspace + +using namespace amrex; +// ------------------------------------------------------------ +// Read variable names from a file +// ------------------------------------------------------------ +Vector +ReadVarNames(const std::string& filename) +{ + Vector varnames; + std::ifstream infile(filename); + if (!infile.is_open()) { + amrex::Abort("Cannot open variable file: " + filename); + } + + std::string line; + while (std::getline(infile, line)) { + // trim whitespace + line.erase(line.begin(), + std::find_if(line.begin(), line.end(), + [](unsigned char ch){ return !std::isspace(ch); })); + line.erase(std::find_if(line.rbegin(), line.rend(), + [](unsigned char ch){ return !std::isspace(ch); }).base(), + line.end()); + + if (!line.empty()) varnames.push_back(line); + } + + if (varnames.empty()) { + amrex::Abort("Variable file is empty: " + filename); + } + + return varnames; +} + +// Reads the plotfile data into cell cenetred multifab +// Does not fill ghost cells +void +ReadPlotFile(const std::string& var_filename, + PlotFileData& pf, + amrex::MultiFab& mf) +{ + // ------------------------------------------------------------ + // Read variable list from file + // ------------------------------------------------------------ + Vector varnames; + { + std::ifstream infile(var_filename); + if (!infile.is_open()) { + Abort("ReadPlotFile: Cannot open variable file: " + var_filename); + } + + std::string line; + while (std::getline(infile, line)) { + if (!line.empty()) { + // trim whitespace + line.erase(0, line.find_first_not_of(" \t\r\n")); + line.erase(line.find_last_not_of(" \t\r\n") + 1); + + if (!line.empty()) varnames.push_back(line); + } + } + } + + if (varnames.empty()) { + Abort("ReadPlotFile: No variables found in " + var_filename); + } + + // ------------------------------------------------------------ + // Open plotfile + // ------------------------------------------------------------ + const Vector& var_names_pf = pf.varNames(); + + // ------------------------------------------------------------ + // Validate requested variables + // ------------------------------------------------------------ + for (auto const& v : varnames) { + bool found = false; + for (auto const& vpf : var_names_pf) { + if (v == vpf) { + found = true; + break; + } + } + if (!found) { + Abort("ReadPlotFile: invalid variable name: " + v); + } + } + + // ------------------------------------------------------------ + // Define destination MultiFab (single level only) + // ------------------------------------------------------------ + const int level = 0; + + BoxArray ba = pf.boxArray(level); + DistributionMapping dm(ba); + + int ncomp = varnames.size(); + int ngrow = 1; + + mf.define(ba, dm, ncomp, ngrow); + + // ------------------------------------------------------------ + // Copy plotfile data → mf + // ------------------------------------------------------------ + for (int comp = 0; comp < ncomp; ++comp) + { + const MultiFab& src = pf.get(level, varnames[comp]); + MultiFab::Copy(mf, src, 0, comp, 1, 0); + } +} + +void ApplyNeumannBCs(const Geometry& geom, + MultiFab& mf_cc) +{ + + // ------------------------------------------------- + // 2. Fill interior + periodic ghost cells + // ------------------------------------------------- + mf_cc.FillBoundary(geom.periodicity()); + // ------------------------------------------------- + // 3. Apply FOExtrap (Neumann) at domain boundaries + // ------------------------------------------------- + const Box& domain = geom.Domain(); + + for (MFIter mfi(mf_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& gbx = mfi.growntilebox(); // includes ghost cells + const Box& vbx = mfi.validbox(); + + auto const& arr = mf_cc.array(mfi); + int ncomp = mf_cc.nComp(); + + ParallelFor(gbx, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + if (vbx.contains(i,j,k)) return; + + int ii = i; + int jj = j; + int kk = k; + + // Clamp to domain interior (FOExtrap) + ii = amrex::max(domain.smallEnd(0), + amrex::min(i, domain.bigEnd(0))); + + jj = amrex::max(domain.smallEnd(1), + amrex::min(j, domain.bigEnd(1))); + + kk = amrex::max(domain.smallEnd(2), + amrex::min(k, domain.bigEnd(2))); + + arr(i,j,k,n) = arr(ii,jj,kk,n); + }); + } +} + + +void +CreateNodalMultiFabFromCellCenteredMultiFab (MultiFab& mf_nc, // output nodal MF + MultiFab& mf_cc, // input cell-centered MF (coarse) + const Geometry& geom) +{ + + // ------------------------------------------------- + // 1. Build nodal MultiFab if not already defined + // ------------------------------------------------- + if (!mf_nc.isDefined()) + { + BoxArray ba_nd = amrex::convert(mf_cc.boxArray(), + IntVect::TheNodeVector()); + + mf_nc.define(ba_nd, + mf_cc.DistributionMap(), + mf_cc.nComp(), + 0); // nodal MF typically needs no ghosts + } + // ------------------------------------------------- + // 2. Fill interior + periodic ghost cells + // ------------------------------------------------- + mf_cc.FillBoundary(geom.periodicity()); + // ------------------------------------------------- + // 3. Apply FOExtrap (Neumann) at domain boundaries + // ------------------------------------------------- + const Box& domain = geom.Domain(); + + for (MFIter mfi(mf_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& gbx = mfi.growntilebox(); // includes ghost cells + const Box& vbx = mfi.validbox(); + + auto const& arr = mf_cc.array(mfi); + int ncomp = mf_cc.nComp(); + + ParallelFor(gbx, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + if (vbx.contains(i,j,k)) return; + + int ii = i; + int jj = j; + int kk = k; + + // Clamp to domain interior (FOExtrap) + ii = amrex::max(domain.smallEnd(0), + amrex::min(i, domain.bigEnd(0))); + + jj = amrex::max(domain.smallEnd(1), + amrex::min(j, domain.bigEnd(1))); + + kk = amrex::max(domain.smallEnd(2), + amrex::min(k, domain.bigEnd(2))); + + arr(i,j,k,n) = arr(ii,jj,kk,n); + }); + } + + AMREX_ALWAYS_ASSERT(mf_nc.nComp() == mf_cc.nComp()); + + int ncomp = mf_cc.nComp(); + + for (MFIter mfi(mf_nc, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); // nodal valid box + + auto const& arr_nc = mf_nc.array(mfi); + auto const& arr_cc = mf_cc.array(mfi); + + ParallelFor(bx, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) + { + // 3D: eight surrounding cells + arr_nc(i,j,k,n) = + 0.125 * ( arr_cc(i-1,j-1,k-1,n) + arr_cc(i ,j-1,k-1,n) + + arr_cc(i-1,j ,k-1,n) + arr_cc(i ,j ,k-1,n) + + arr_cc(i-1,j-1,k ,n) + arr_cc(i ,j-1,k ,n) + + arr_cc(i-1,j ,k ,n) + arr_cc(i ,j ,k ,n) ); + }); + } +} + + + + +void +CreateCellCenteredMultiFabFromNodalMultiFab(MultiFab& mf_cc_tmp, // cell-centered output + MultiFab& mf_nc) // node-centered input +{ + + // ------------------------------------------------- + // 1. Define cell-centered MF with same structure + // ------------------------------------------------- + BoxArray ba_cc = amrex::convert(mf_nc.boxArray(), + IntVect::TheCellVector()); + + mf_cc_tmp.define(ba_cc, + mf_nc.DistributionMap(), + mf_nc.nComp(), + 0); // no ghosts needed for simple averaging + + int ncomp = mf_cc_tmp.nComp(); + + for (MFIter mfi(mf_cc_tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& cbox = mfi.tilebox(); // cell-centered valid box + auto const& arr_cc = mf_cc_tmp.array(mfi); + auto const& arr_nc = mf_nc.array(mfi); + + ParallelFor(cbox, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) + { + // Average 8 surrounding nodes + arr_cc(i,j,k,n) = 0.125 * ( arr_nc(i,j,k,n) // corner 0 + + arr_nc(i+1,j,k,n) // corner 1 + + arr_nc(i,j+1,k,n) // corner 2 + + arr_nc(i+1,j+1,k,n) // corner 3 + + arr_nc(i,j,k+1,n) // corner 4 + + arr_nc(i+1,j,k+1,n) // corner 5 + + arr_nc(i,j+1,k+1,n) // corner 6 + + arr_nc(i+1,j+1,k+1,n) // corner 7 + ); + }); + } +} + + +enum class BoundType { Lo, Hi }; +enum class MultiFabType { CC, NC }; + +IntVect +find_bound_idx(const Real& x, const Real& y, const Real& z, + const BoxList& bl_coarse, const Geometry& geom_coarse, + BoundType bound_type) +{ + const auto prob_lo_coarse = geom_coarse.ProbLoArray(); + const auto dx_coarse = geom_coarse.CellSizeArray(); + + int i, j, k; + + if (bound_type == BoundType::Lo) { + i = static_cast(std::floor((x - prob_lo_coarse[0]) / dx_coarse[0])); + j = static_cast(std::floor((y - prob_lo_coarse[1]) / dx_coarse[1])); + k = static_cast(std::floor((z - prob_lo_coarse[2]) / dx_coarse[2])); + } else { // BoundType::Hi + i = static_cast(std::ceil((x - prob_lo_coarse[0]) / dx_coarse[0])); + j = static_cast(std::ceil((y - prob_lo_coarse[1]) / dx_coarse[1])); + k = static_cast(std::ceil((z - prob_lo_coarse[2]) / dx_coarse[2])); + } + + IntVect idx(i, j, k); + + for (const auto& b : bl_coarse) { + if (b.contains(idx)) { + return idx; + } + } + + amrex::Print() << x << " " << y << " " << z << " " << idx << std::endl; + + amrex::Print() << "Printing BoxList (coarse):\n"; +for (const auto& b : bl_coarse) { + amrex::Print() << b << "\n"; +} + + amrex::Abort("Bound index not found in any box in BoxList!"); + return IntVect::TheZeroVector(); // unreachable if Abort +} + + +void +GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, + const Geometry& geom_fine, + const MultiFab& mf_nc_coarse, + const MultiFab& mf_cc_fine, + MultiFab& coarse_multifab_on_fine_dmap) +{ + BoxList bl_coarse = mf_nc_coarse.boxArray().boxList(); + BoxList bl_fine = mf_cc_fine.boxArray().boxList(); + + const auto prob_lo_fine = geom_fine.ProbLoArray(); + const auto dx_fine = geom_fine.CellSizeArray(); + + for (auto& b : bl_fine) { + // You look at the lo corner of b, and find out the lowest cell in + // coarse mutlifab data you need for the interpolation. That gives + // you the lo corner of the new b. Similarly, you can find out the + // hi corner of the new b. For cells outside the coarse multifab data's + // bounding data, it's up to you. You probably want to use a biased + // interpolation stencil. + + // Get the cell indices of the bottom corner and top corner + const IntVect& lo_fine = b.smallEnd(); // Lower corner (inclusive) + const IntVect& hi_fine = b.bigEnd(); // Upper corner (inclusive) + + Real x = prob_lo_fine[0] + lo_fine[0] * dx_fine[0]; + Real y = prob_lo_fine[1] + lo_fine[1] * dx_fine[1]; + Real z = prob_lo_fine[2] + lo_fine[2] * dx_fine[2]; + + auto idx_lo = find_bound_idx(x, y, z, bl_coarse, geom_coarse, BoundType::Lo); + + + x = prob_lo_fine[0] + hi_fine[0] * dx_fine[0]; + y = prob_lo_fine[1] + hi_fine[1] * dx_fine[1]; + z = prob_lo_fine[2] + hi_fine[2] * dx_fine[2]; + + auto idx_hi = find_bound_idx(x, y, z, bl_coarse, geom_coarse, BoundType::Hi); + + b.setSmall(idx_lo); + b.setBig(idx_hi); + + /*Print() << "lo fine = " << lo_fine << std::endl; + Print() << "hi fine = " << hi_fine << std::endl; + Print() << " idx lo = " << idx_lo << std::endl; + Print() << "idx_hi = " << idx_hi << std::endl;*/ + + } + + BoxArray cba(std::move(bl_fine)); + cba.convert(IndexType::TheNodeType()); // <-- Make it nodal in all directions + coarse_multifab_on_fine_dmap.define(cba, mf_cc_fine.DistributionMap(), mf_nc_coarse.nComp(), 0); + coarse_multifab_on_fine_dmap.ParallelCopy(mf_nc_coarse); +} + + +void +PopulateFineCellCenteredFromCoarseNodal(const Geometry& geom_coarse, + const Geometry& geom_fine, + const MultiFab& coarse_multifab_on_fine_dmap, + const MultiFab& mf_cc_fine, + MultiFab& mf_cc_from_coarse) +{ + AMREX_ALWAYS_ASSERT(coarse_multifab_on_fine_dmap.ixType().nodeCentered()); + + if (!mf_cc_from_coarse.isDefined()) + { + mf_cc_from_coarse.define(mf_cc_fine.boxArray(), + mf_cc_fine.DistributionMap(), + coarse_multifab_on_fine_dmap.nComp(), + 0); + } + + AMREX_ALWAYS_ASSERT(mf_cc_from_coarse.boxArray() == mf_cc_fine.boxArray()); + AMREX_ALWAYS_ASSERT(mf_cc_from_coarse.DistributionMap() == mf_cc_fine.DistributionMap()); + AMREX_ALWAYS_ASSERT(mf_cc_from_coarse.nComp() == coarse_multifab_on_fine_dmap.nComp()); + + const auto prob_lo_coarse = geom_coarse.ProbLoArray(); + const auto dx_coarse = geom_coarse.CellSizeArray(); + + const auto prob_lo_fine = geom_fine.ProbLoArray(); + const auto dx_fine = geom_fine.CellSizeArray(); + + Box nodal_domain = amrex::convert(geom_coarse.Domain(), IntVect::TheNodeVector()); + const auto nd_lo = nodal_domain.smallEnd(); + const auto nd_hi = nodal_domain.bigEnd(); + + int ncomp = mf_cc_from_coarse.nComp(); + + for (MFIter mfi(mf_cc_from_coarse, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& cbox = mfi.tilebox(); + auto const& arr_cc = mf_cc_from_coarse.array(mfi); + auto const& arr_nc = coarse_multifab_on_fine_dmap.array(mfi); + + ParallelFor(cbox, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) + { + Real x = prob_lo_fine[0] + (static_cast(i) + 0.5_rt) * dx_fine[0]; + Real y = prob_lo_fine[1] + (static_cast(j) + 0.5_rt) * dx_fine[1]; + Real z = prob_lo_fine[2] + (static_cast(k) + 0.5_rt) * dx_fine[2]; + + Real fx = (x - prob_lo_coarse[0]) / dx_coarse[0]; + Real fy = (y - prob_lo_coarse[1]) / dx_coarse[1]; + Real fz = (z - prob_lo_coarse[2]) / dx_coarse[2]; + + int i0 = static_cast(amrex::Math::floor(fx)); + int j0 = static_cast(amrex::Math::floor(fy)); + int k0 = static_cast(amrex::Math::floor(fz)); + + Real wx = fx - static_cast(i0); + Real wy = fy - static_cast(j0); + Real wz = fz - static_cast(k0); + + int i1 = i0 + 1; + int j1 = j0 + 1; + int k1 = k0 + 1; + + i0 = amrex::max(nd_lo[0], amrex::min(i0, nd_hi[0]-1)); + j0 = amrex::max(nd_lo[1], amrex::min(j0, nd_hi[1]-1)); + k0 = amrex::max(nd_lo[2], amrex::min(k0, nd_hi[2]-1)); + + i1 = amrex::max(nd_lo[0], amrex::min(i1, nd_hi[0])); + j1 = amrex::max(nd_lo[1], amrex::min(j1, nd_hi[1])); + k1 = amrex::max(nd_lo[2], amrex::min(k1, nd_hi[2])); + + Real c000 = arr_nc(i0,j0,k0,n); + Real c100 = arr_nc(i1,j0,k0,n); + Real c010 = arr_nc(i0,j1,k0,n); + Real c110 = arr_nc(i1,j1,k0,n); + Real c001 = arr_nc(i0,j0,k1,n); + Real c101 = arr_nc(i1,j0,k1,n); + Real c011 = arr_nc(i0,j1,k1,n); + Real c111 = arr_nc(i1,j1,k1,n); + + Real c00 = c000 * (1.0_rt - wx) + c100 * wx; + Real c10 = c010 * (1.0_rt - wx) + c110 * wx; + Real c01 = c001 * (1.0_rt - wx) + c101 * wx; + Real c11 = c011 * (1.0_rt - wx) + c111 * wx; + + Real c0 = c00 * (1.0_rt - wy) + c10 * wy; + Real c1 = c01 * (1.0_rt - wy) + c11 * wy; + + arr_cc(i,j,k,n) = c0 * (1.0_rt - wz) + c1 * wz; + }); + } +} + +void WriteCustomDataFile(const Geometry& geom, + const MultiFab& mf_cc, + const std::string& filename_custom) +{ + AMREX_ALWAYS_ASSERT(mf_cc.nComp() >= 1); + AMREX_ALWAYS_ASSERT(mf_cc.local_size() == 1); // single box assumed + + const int ncomp = mf_cc.nComp(); + const int ng = mf_cc.nGrow(); + + const Box grown_domain = amrex::grow(geom.Domain(), ng); + const auto lo = lbound(grown_domain); + const auto hi = ubound(grown_domain); + + const auto problo = geom.ProbLoArray(); + const auto probhi = geom.ProbHiArray(); + const auto dx = geom.CellSizeArray(); + + // extend physical domain to include ghosts + amrex::GpuArray problo_ext, probhi_ext; + for (int d=0; d<3; ++d) { + problo_ext[d] = problo[d] - ng*dx[d]; + probhi_ext[d] = probhi[d] + ng*dx[d]; + } + + std::ofstream ofs(filename_custom, std::ios::binary); + if (!ofs.is_open()) Abort("Failed to open file for writing"); + + // header + const int nx = hi.x - lo.x + 1; + const int ny = hi.y - lo.y + 1; + const int nz = hi.z - lo.z + 1; + + ofs.write(reinterpret_cast(&nx), sizeof(int)); + ofs.write(reinterpret_cast(&ny), sizeof(int)); + ofs.write(reinterpret_cast(&nz), sizeof(int)); + ofs.write(reinterpret_cast(&ng), sizeof(int)); + ofs.write(reinterpret_cast(&ncomp), sizeof(int)); + + for (int d=0; d<3; ++d) ofs.write(reinterpret_cast(&problo_ext[d]), sizeof(double)); + for (int d=0; d<3; ++d) ofs.write(reinterpret_cast(&probhi_ext[d]), sizeof(double)); + + // data + const auto& arr = mf_cc.const_array(0); + for (int k=lo.z; k<=hi.z; ++k) + for (int j=lo.y; j<=hi.y; ++j) + for (int i=lo.x; i<=hi.x; ++i) + { + double x = problo[0] + (i + 0.5) * dx[0]; + double y = problo[1] + (j + 0.5) * dx[1]; + double z = problo[2] + (k + 0.5) * dx[2]; + + ofs.write(reinterpret_cast(&x), sizeof(double)); + ofs.write(reinterpret_cast(&y), sizeof(double)); + ofs.write(reinterpret_cast(&z), sizeof(double)); + + for (int n=0; n(&val), sizeof(double)); + } + } + + ofs.close(); +} diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp new file mode 100644 index 000000000..735e4d10d --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp @@ -0,0 +1,69 @@ +#include +#include +#include + +#include "ReadPlotFile.H" + +using namespace amrex; + +int main (int argc, char* argv[]) +{ + Initialize(argc, argv); + + { + MultiFab mf_cc_coarse; + MultiFab mf_cc_fine; + + PlotFileData pf_coarse("plt_coarse"); + ReadPlotFile("vars.txt", pf_coarse, mf_cc_coarse); + + PlotFileData pf_fine("plt_fine"); + ReadPlotFile("vars.txt", pf_fine, mf_cc_fine); + + Array is_periodic{AMREX_D_DECL(0,0,0)}; + + Geometry geom_coarse(pf_coarse.probDomain(0), + RealBox(pf_coarse.probLo(), pf_coarse.probHi()), + pf_coarse.coordSys(), + is_periodic); + + Geometry geom_fine(pf_fine.probDomain(0), + RealBox(pf_fine.probLo(), pf_fine.probHi()), + pf_fine.coordSys(), + is_periodic); + + MultiFab mf_nc_coarse; + CreateNodalMultiFabFromCellCenteredMultiFab(mf_nc_coarse, mf_cc_coarse, geom_coarse); + + MultiFab mf_cc_tmp; + CreateCellCenteredMultiFabFromNodalMultiFab(mf_cc_tmp, mf_nc_coarse); + // Read variable names + Vector varnames = ReadVarNames("vars.txt"); + + WriteSingleLevelPlotfile("plt_1", mf_cc_tmp, varnames, geom_coarse, 0.0, 0); + + MultiFab coarse_multifab_on_fine_dmap; + GetCoarseMultiFabOnFineDMap(geom_coarse, geom_fine, + mf_nc_coarse, mf_cc_fine, + coarse_multifab_on_fine_dmap); + + Print() << "Checking for large values on coarse_multifab_on_fine_dmap" << std::endl; + //check_large_values(coarse_multifab_on_fine_dmap); + + CreateCellCenteredMultiFabFromNodalMultiFab(mf_cc_tmp, coarse_multifab_on_fine_dmap); + Print() << "Checking for large values on mf_cc_tmp" << std::endl; + //check_large_values(mf_cc_tmp); + + // Write plotfile + WriteSingleLevelPlotfile("plt_2", mf_cc_tmp, varnames, geom_coarse, 0.0, 0); + + MultiFab mf_cc_fine_from_coarse; + PopulateFineCellCenteredFromCoarseNodal(geom_coarse, geom_fine, coarse_multifab_on_fine_dmap, + mf_cc_fine, mf_cc_fine_from_coarse); + + WriteSingleLevelPlotfile("plt_final", mf_cc_fine_from_coarse, varnames, geom_fine, 0.0, 0); + + } + Finalize(); +} + diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp new file mode 100644 index 000000000..ad4fc4732 --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +#include "ReadPlotFile.H" + +using namespace amrex; + +int main (int argc, char* argv[]) +{ + Initialize(argc, argv); + + { + MultiFab mf_cc_coarse; + MultiFab mf_cc_fine; + + PlotFileData pf_coarse("plt_coarse"); + ReadPlotFile("vars.txt", pf_coarse, mf_cc_coarse); + + PlotFileData pf_fine("plt_fine"); + ReadPlotFile("vars.txt", pf_fine, mf_cc_fine); + + Array is_periodic{AMREX_D_DECL(0,0,0)}; + + Geometry geom_coarse(pf_coarse.probDomain(0), + RealBox(pf_coarse.probLo(), pf_coarse.probHi()), + pf_coarse.coordSys(), + is_periodic); + + Geometry geom_fine(pf_fine.probDomain(0), + RealBox(pf_fine.probLo(), pf_fine.probHi()), + pf_fine.coordSys(), + is_periodic); + + + std::string filename_custom = "coarse_data.bin"; + ApplyNeumannBCs(geom_coarse, mf_cc_coarse); + WriteCustomDataFile(geom_coarse, mf_cc_coarse, filename_custom); + + // Read variable names + Vector varnames = ReadVarNames("vars.txt"); + + //WriteSingleLevelPlotfile("plt_final", mf_cc_fine_from_coarse, varnames, geom_fine, 0.0, 0); + + } + Finalize(); +} + diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/run_interactive_GPU.sh b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/run_interactive_GPU.sh new file mode 100644 index 000000000..3374a93d4 --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/run_interactive_GPU.sh @@ -0,0 +1,2 @@ +salloc --nodes 4 --qos interactive --time 04:00:00 --constraint gpu --gpus 16 --account=m5175 + diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/vars.txt b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/vars.txt new file mode 100644 index 000000000..5c98397ff --- /dev/null +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/vars.txt @@ -0,0 +1,5 @@ +density +theta +x_velocity +y_velocity +z_velocity diff --git a/.Exec_dev/DataAssimilation/inputs_advecting_coarse b/.Exec_dev/DataAssimilation/inputs_advecting_coarse index cbf8e2717..af3ff9bf3 100644 --- a/.Exec_dev/DataAssimilation/inputs_advecting_coarse +++ b/.Exec_dev/DataAssimilation/inputs_advecting_coarse @@ -6,8 +6,8 @@ amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_lo = 0 0 -1 -geometry.prob_hi = 50 25 1 +geometry.prob_lo = 0 0 -1.0 +geometry.prob_hi = 50.0 25.0 1.0 amr.n_cell = 64 32 4 geometry.is_periodic = 1 1 0 diff --git a/.Exec_dev/DataAssimilation/main.cpp b/.Exec_dev/DataAssimilation/main.cpp index 6e2d761d1..370d4136a 100644 --- a/.Exec_dev/DataAssimilation/main.cpp +++ b/.Exec_dev/DataAssimilation/main.cpp @@ -129,6 +129,7 @@ int main (int argc, char* argv[]) fs::create_directory(member_dir); fs::create_directory(member_dir + "/plotfiles"); fs::create_directory(member_dir + "/chkfiles"); + fs::create_directory(member_dir + "/pertfiles"); // Move plotfiles (plt*) from current directory into member_dir/plotfiles for (auto &p : fs::directory_iterator(".")) { @@ -148,9 +149,11 @@ int main (int argc, char* argv[]) } // Optional: barrier after move to ensure rank 0 is done ParallelDescriptor::Barrier(); - } + ERF tmp_erf; + tmp_erf.ComputeAndWriteEnsemblePerturbations(); + BL_PROFILE_VAR_STOP(pmain); amrex::Finalize(); @@ -161,4 +164,3 @@ int main (int argc, char* argv[]) return 0; } - diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 2373f2ed5..702a4daa8 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -319,6 +319,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Initialization/ERF_Init1D.cpp ${SRC_DIR}/Initialization/ERF_InitTurbPert.cpp ${SRC_DIR}/Initialization/ERF_InitImmersedForcing.cpp + ${SRC_DIR}/Initialization/ERF_InitForEnsemble.cpp ${SRC_DIR}/IO/ERF_Checkpoint.cpp ${SRC_DIR}/IO/ERF_ReadBndryPlanes.cpp ${SRC_DIR}/IO/ERF_WriteBndryPlanes.cpp diff --git a/Exec/ERF_Prob.cpp b/Exec/ERF_Prob.cpp index 151291eb5..2fe030087 100644 --- a/Exec/ERF_Prob.cpp +++ b/Exec/ERF_Prob.cpp @@ -139,6 +139,9 @@ Problem::init_custom_pert ( } else if (my_prob_name_ci == "supercell") { #include "Prob/ERF_InitCustomPert_SuperCell.H" + } + else if (my_prob_name_ci == "data_assimilation_isv") { +#include "Prob/ERF_InitCustomPert_DataAssimilation_ISV.H" } else if (my_prob_name_ci == "sinusoidalmassflux") { #include "Prob/ERF_InitCustomPert_Bomex.H" @@ -227,6 +230,9 @@ Problem::init_custom_pert_vels ( } else if (my_prob_name_ci == "userdefined") { #include "Prob/ERF_InitCustomPertVels_UserDefined.H" + } + else if (my_prob_name_ci == "data_assimilation_isv") { +#include "Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H" } else if (my_prob_name_ci == "sinusoidalmassflux") { #include "Prob/ERF_InitCustomPertVels_Bomex.H" diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index fd80f8c4d..7c975f723 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -754,15 +754,33 @@ struct SolverChoice { } } - amrex::ParmParse pp_ens("ensemble"); - pp_ens.query("is_init_with_correlated_pert", is_init_with_correlated_pert); - if(is_init_with_correlated_pert) { - pp_ens.query("pert_correlated_radius", pert_correlated_radius); - if(pert_correlated_radius <= zero) { + pp.query("is_init_for_ensemble", is_init_for_ensemble); + if(is_init_for_ensemble) { + amrex::ParmParse pp_ens("ensemble"); + pp_ens.query("n_members", n_ensemble); + pp_ens.query("coarse_bckgnd_data_file", coarse_bckgnd_data_file); + pp_ens.query("ens_pert_amplitude", ens_pert_amplitude); + pp_ens.query("ens_pert_correlated_radius", ens_pert_correlated_radius); + if (n_ensemble < 0) { + amrex::Abort("You are running an ensemble simulation, There needs to be an entry in the inputs " + "ensemble.n_ensemble which is the number of ensemble members, and should be greater than 0."); + } + if (coarse_bckgnd_data_file.empty()) { + amrex::Abort("coarse_bckgnd_data_file is empty! For ensmeble simulations, there needs to be a coarse background file which " + "contains the data which will be interpolated onto the fine mesh. There has to a entry ensemble.coarse_bckgnd_data_file " + "which contains the filename in the inputs."); + } + if(ens_pert_amplitude <= 0.0) { + amrex::Error("You are using initialization for ensemble simulations using the inputs option " + "ensemble.is_init_for_ensemble=true. In this case, there has to be an option " + "ensemble.ens_pert_amplitude which is the value of the amplitude of the perturbation " + "to be added to the background state and has to be greater than 0.0"); + } + if(ens_pert_correlated_radius <= 0.0) { amrex::Error("You are using initialization with spatially correlated perturbations using the inputs option " - "ensemble.is_init_with_correlated_pert=true. In this case, there has to be an option " - "ensemble.pert_correlated_radius which is the value of the the spatial correlation radius, " - "and has to be greater than zero"); + "ensemble.is_init_for_ensemble=true. In this case, there has to be an option " + "ensemble.ens_pert_correlated_radius which is the value of the the spatial correlation radius, " + "and has to be greater than 0.0"); } } } @@ -1282,7 +1300,10 @@ struct SolverChoice { bool io_hurricane_eye_tracker = false; amrex::Real hurricane_eye_latitude = -1e10, hurricane_eye_longitude = -1e10; - bool is_init_with_correlated_pert = false; - amrex::Real pert_correlated_radius = -one; + bool is_init_for_ensemble = false; + int n_ensemble = -1; + amrex::Real ens_pert_correlated_radius = -1.0; + amrex::Real ens_pert_amplitude = -1.0; + std::string coarse_bckgnd_data_file; }; #endif diff --git a/Source/ERF.H b/Source/ERF.H index 23e94e4b4..08bc5eb4f 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -167,6 +167,8 @@ public: surface_state_2, surface_state_interp; + amrex::Vector> bckgnd_state; + std::string MakeVTKFilename(int nstep); std::string MakeVTKFilename_TrackerCircle(int nstep); std::string MakeVTKFilename_EyeTracker_xy(int nstep); @@ -496,17 +498,20 @@ public: amrex::Vector>& z_phys_nd, bool regrid_forces_file_read); + void create_background_state_for_ensemble (int lev, + amrex::MultiFab& mf_cc_pert, + amrex::MultiFab& cons_pert, + amrex::MultiFab& xvel_pert, + amrex::MultiFab& yvel_pert, + amrex::MultiFab& zvel_pert); + void create_random_perturbations (const int lev, - amrex::MultiFab& cons_pert, - amrex::MultiFab& xvel_pert, - amrex::MultiFab& yvel_pert, - amrex::MultiFab& zvel_pert); + amrex::MultiFab& mf_cc_pert); void apply_gaussian_smoothing_to_perturbations (const int lev, - amrex::MultiFab& cons_pert, - amrex::MultiFab& xvel_pert, - amrex::MultiFab& yvel_pert, - amrex::MultiFab& zvel_pert); + amrex::MultiFab& mf_cc_pert); + + void ComputeAndWriteEnsemblePerturbations(); void init_custom (int lev); diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 81866facc..74a7bf692 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -42,13 +42,6 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - // If the initial condition needs to have spatially correlated perturbations superimposed onto - // the base state, then populate the "pert" multifabs with random perturbations, - // and then apply a Gaussian smoothing to make the perturbations spatially correlated - if(solverChoice.is_init_with_correlated_pert) { - create_random_perturbations(lev, cons_pert, xvel_pert, yvel_pert, zvel_pert); - apply_gaussian_smoothing_to_perturbations(lev, cons_pert, xvel_pert, yvel_pert, zvel_pert); - } #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -120,6 +113,7 @@ ERF::init_custom (int lev) } //mfi + // Add problem-specific perturbation to background flow if not doing anelastic with fixed-in-time density if (!solverChoice.fixed_density[lev]) { MultiFab::Add(lev_new[Vars::cons], cons_pert, Rho_comp, Rho_comp, 1, cons_pert.nGrow()); @@ -150,104 +144,18 @@ ERF::init_custom (int lev) MultiFab::Add(lev_new[Vars::yvel], yvel_pert, 0, 0, 1, yvel_pert.nGrowVect()); MultiFab::Add(lev_new[Vars::zvel], zvel_pert, 0, 0, 1, zvel_pert.nGrowVect()); } -} -void -ERF::create_random_perturbations(const int lev, - MultiFab& cons_pert, - MultiFab& xvel_pert, - MultiFab& yvel_pert, - MultiFab& zvel_pert) -{ - ignore_unused(cons_pert); - ignore_unused(yvel_pert); - ignore_unused(zvel_pert); - - auto& lev_new = vars_new[lev]; - for (MFIter mfi(lev_new[Vars::cons], TileNoZ()); mfi.isValid(); ++mfi) { - const auto &xvel_pert_arr = xvel_pert.array(mfi); - const Box &xbx = mfi.tilebox(IntVect(1,0,0)); - ParallelForRNG(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept - { - xvel_pert_arr(i, j, k) = amrex::Random(engine); - }); - } -} - -void -ERF::apply_gaussian_smoothing_to_perturbations(const int lev, - MultiFab& cons_pert, - MultiFab& xvel_pert, - MultiFab& yvel_pert, - MultiFab& zvel_pert) -{ - ignore_unused(cons_pert); - ignore_unused(yvel_pert); - ignore_unused(zvel_pert); - - const Geometry& gm = geom[lev]; - const Real dx = gm.CellSize(0); - const Real dy = gm.CellSize(1); - - const Real dmesh = std::min(dx, dy); - // ---- User choices ---- - const Real sigma = solverChoice.pert_correlated_radius; // e.g. 2 km correlation length - const int r = static_cast(three * sigma / dmesh); // stencil radius - - // ---- Precompute Gaussian weights on host ---- - const int wsize = 2*r + 1; - Vector w_host(wsize * wsize); - - Real Z = zero; - for (int m = -r; m <= r; ++m) { - for (int n = -r; n <= r; ++n) { - Real val = std::exp(-(m*m*dx*dx + n*n*dy*dy)/(two*sigma*sigma)); - w_host[(m+r)*wsize + (n+r)] = val; - Z += val; - } - } - for (auto& v : w_host) { - v = v/Z; - } - - Gpu::DeviceVector w_dev; - w_dev.resize(w_host.size()); - Gpu::copy(Gpu::hostToDevice, w_host.begin(), w_host.end(), w_dev.begin()); - - Real const* w = w_dev.data(); - - // one Define ngrow_big using the actual dimension macro - IntVect ngrow_big(AMREX_D_DECL(r, r, 0)); - - // two Create the copy - MultiFab xvel_pert_copy(xvel_pert.boxArray(), - xvel_pert.DistributionMap(), - 1, ngrow_big); - //MultiFab::Copy(xvel_pert_copy, xvel_pert, 0, 0, 1, 0); - - // three Use the built-in copy that includes ghost cell logic - // Copy(dst, src, src_comp, dst_comp, num_comp, ngrow) - // Setting ngrow to 0 ensures we only take valid data from the original - xvel_pert_copy.ParallelCopy(xvel_pert, 0, 0, 1, IntVect(0), ngrow_big, gm.periodicity()); - - for (MFIter mfi(xvel_pert, TileNoZ()); mfi.isValid(); ++mfi) - { - const Box& tbx = mfi.tilebox(); - - auto const& in = xvel_pert_copy.array(mfi); - auto const& out = xvel_pert.array(mfi); - - ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real sum = zero; - for (int m = -r; m <= r; ++m) { - for (int n = -r; n <= r; ++n) { - Real wij = w[(m+r)*wsize + (n+r)]; - sum += wij * in(i+m, j+n, k); - } - } - out(i,j,k) = sum; - }); + // If initializing for ensemble simluations, then + // 1. Create cell-centered random perturbations + // 2. Create cell-centered spatially correlated perturbations + // 3. Read in the coarse background state from the coarse data file, + // interpolate the state data onto the current mesh, and + // add the perturbations to the background state and then populate the "pert variables + + if(solverChoice.is_init_for_ensemble) { + MultiFab mf_cc_pert; + create_random_perturbations(lev, mf_cc_pert); + apply_gaussian_smoothing_to_perturbations(lev, mf_cc_pert); + create_background_state_for_ensemble(lev, mf_cc_pert, lev_new[Vars::cons], lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::zvel]); } } diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp new file mode 100644 index 000000000..2d6b35832 --- /dev/null +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -0,0 +1,734 @@ +/** + * \file ERF_InitForEnsemble.cpp + */ + +#include +#include +#include + +using namespace amrex; +namespace fs = std::filesystem; + +void +ERF::create_random_perturbations(const int lev, + MultiFab& mf_cc_pert) +{ + const MultiFab& src = vars_new[lev][Vars::cons]; + + int ncomp = 5; + mf_cc_pert.define(src.boxArray(), src.DistributionMap(), + ncomp, src.nGrow()); + + // Loop over cell-centered boxes + for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + auto const& pert_arr = mf_cc_pert.array(mfi); + + // Loop over all 5 components + amrex::ParallelForRNG(bx, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n, + const amrex::RandomEngine& engine) noexcept + { + pert_arr(i,j,k,n) = amrex::Random(engine); + }); + } +} + +void NormalizeMultiFabRMS_PerComponent(MultiFab& mf_cc_pert) +{ + const int ncomp = mf_cc_pert.nComp(); + + for (int n = 0; n < ncomp; ++n) + { + // 1. Set up AMReX reduction (sum of squares + count) + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + // 2. Loop over tiles and accumulate + for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& arr = mf_cc_pert.const_array(mfi); + + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> ReduceTuple + { + Real v = arr(i, j, k, n); + return { v * v, 1L }; + }); + } + + // 3. Retrieve results (includes implicit GPU sync + host copy) + auto rv = reduce_data.value(reduce_op); + Real h_sumsq = amrex::get<0>(rv); + Long h_count = amrex::get<1>(rv); + + // 4. Sum across MPI ranks + ParallelDescriptor::ReduceRealSum(h_sumsq); + ParallelDescriptor::ReduceLongSum(h_count); + + // 5. Compute RMS and normalize + if (h_count > 0) + { + Real rms = std::sqrt(h_sumsq / static_cast(h_count)); + if (rms > 0.0) { + mf_cc_pert.mult(1.0 / rms, n, 1); + } + } + } +} + +void +ERF::apply_gaussian_smoothing_to_perturbations(const int lev, + MultiFab& mf_cc_pert) +{ + const Geometry& gm = geom[lev]; + const Real dx = gm.CellSize(0); + const Real dy = gm.CellSize(1); + + const Real dmesh = std::min(dx, dy); + + // ---- User choice ---- + const Real sigma = solverChoice.ens_pert_correlated_radius; + const int r = static_cast(3.0 * sigma / dmesh); + + const int ncomp = mf_cc_pert.nComp(); + + // ---- Precompute Gaussian weights ---- + const int wsize = 2*r + 1; + Vector w_host(wsize * wsize); + + Real Z = 0.0; + for (int m = -r; m <= r; ++m) { + for (int n = -r; n <= r; ++n) { + Real val = std::exp(-(m*m*dx*dx + n*n*dy*dy) + /(2.0*sigma*sigma)); + w_host[(m+r)*wsize + (n+r)] = val; + Z += val; + } + } + + for (auto& v : w_host) { + v /= Z; + } + + Gpu::DeviceVector w_dev(w_host.size()); + Gpu::copy(Gpu::hostToDevice, w_host.begin(), w_host.end(), w_dev.begin()); + + Real const* w = w_dev.data(); + + // ---- Create a grown copy (for stencil access) ---- + IntVect ngrow_big(AMREX_D_DECL(r, r, 0)); + + MultiFab mf_copy(mf_cc_pert.boxArray(), + mf_cc_pert.DistributionMap(), + ncomp, ngrow_big); + + mf_copy.ParallelCopy(mf_cc_pert, + 0, 0, ncomp, + IntVect(0), ngrow_big, + gm.periodicity()); + + // ---- Apply smoothing ---- + for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + auto const& in = mf_copy.const_array(mfi); + auto const& out = mf_cc_pert.array(mfi); + + ParallelFor(bx, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + Real sum = 0.0; + + for (int m = -r; m <= r; ++m) { + for (int nn = -r; nn <= r; ++nn) { + Real wij = w[(m+r)*wsize + (nn+r)]; + sum += wij * in(i+m, j+nn, k, n); + } + } + + out(i,j,k,n) = sum; + }); + } + NormalizeMultiFabRMS_PerComponent(mf_cc_pert); +} + +void ApplyNeumannBCs(const Geometry& geom, + MultiFab& mf_cc) +{ + + // ------------------------------------------------- + // 2. Fill interior + periodic ghost cells + // ------------------------------------------------- + mf_cc.FillBoundary(geom.periodicity()); + // ------------------------------------------------- + // 3. Apply FOExtrap (Neumann) at domain boundaries + // ------------------------------------------------- + const Box& domain = geom.Domain(); + + for (MFIter mfi(mf_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& gbx = mfi.growntilebox(); // includes ghost cells + const Box& vbx = mfi.validbox(); + + auto const& arr = mf_cc.array(mfi); + int ncomp = mf_cc.nComp(); + + ParallelFor(gbx, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + if (vbx.contains(i,j,k)) return; + + int ii = i; + int jj = j; + int kk = k; + + // Clamp to domain interior (FOExtrap) + ii = amrex::max(domain.smallEnd(0), + amrex::min(i, domain.bigEnd(0))); + + jj = amrex::max(domain.smallEnd(1), + amrex::min(j, domain.bigEnd(1))); + + kk = amrex::max(domain.smallEnd(2), + amrex::min(k, domain.bigEnd(2))); + + arr(i,j,k,n) = arr(ii,jj,kk,n); + }); + } +} + + +void ReadCustomDataFile(const std::string& filename_custom, + int& nx, int& ny, int& nz, + int& ng, int& ncomp, + std::array& problo_ext, + std::array& probhi_ext, + Vector& data_rho, + Vector& data_theta, + Vector& data_xvel, + Vector& data_yvel, + Vector& data_zvel) +{ + std::ifstream ifs(filename_custom, std::ios::binary); + if (!ifs.is_open()) { + Abort("Failed to open file " + filename_custom + " for reading"); + } + + // ---------------------------- + // Read header + // ---------------------------- + ifs.read(reinterpret_cast(&nx), sizeof(int)); + ifs.read(reinterpret_cast(&ny), sizeof(int)); + ifs.read(reinterpret_cast(&nz), sizeof(int)); + + ifs.read(reinterpret_cast(&ng), sizeof(int)); + ifs.read(reinterpret_cast(&ncomp), sizeof(int)); + + ifs.read(reinterpret_cast(&problo_ext[0]), sizeof(Real)); + ifs.read(reinterpret_cast(&problo_ext[1]), sizeof(Real)); + ifs.read(reinterpret_cast(&problo_ext[2]), sizeof(Real)); + + ifs.read(reinterpret_cast(&probhi_ext[0]), sizeof(Real)); + ifs.read(reinterpret_cast(&probhi_ext[1]), sizeof(Real)); + ifs.read(reinterpret_cast(&probhi_ext[2]), sizeof(Real)); + + const std::size_t ncell = static_cast(nx) * ny * nz; + + // ---------------------------- + // Allocate storage + // ---------------------------- + data_rho.resize(ncell); + data_theta.resize(ncell); + data_xvel.resize(ncell); + data_yvel.resize(ncell); + data_zvel.resize(ncell); + + // ---------------------------- + // Read data + // ---------------------------- + std::size_t idx = 0; + + for (int k = 0; k < nz; ++k) + { + for (int j = 0; j < ny; ++j) + { + for (int i = 0; i < nx; ++i) + { + // Skip coordinates + Real x, y, z; + ifs.read(reinterpret_cast(&x), sizeof(Real)); + ifs.read(reinterpret_cast(&y), sizeof(Real)); + ifs.read(reinterpret_cast(&z), sizeof(Real)); + + // Read components (fixed order) + ifs.read(reinterpret_cast(&data_rho[idx]), sizeof(Real)); + ifs.read(reinterpret_cast(&data_theta[idx]), sizeof(Real)); + ifs.read(reinterpret_cast(&data_xvel[idx]), sizeof(Real)); + ifs.read(reinterpret_cast(&data_yvel[idx]), sizeof(Real)); + ifs.read(reinterpret_cast(&data_zvel[idx]), sizeof(Real)); + + ++idx; + } + } + } + + ifs.close(); +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +int idx(int i, int j, int k, int nx, int ny) +{ + return i + nx * (j + ny * k); +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +Real interp_trilinear( + const Real* f, // <-- raw pointer + int i, int j, int k, + Real tx, Real ty, Real tz, + int nx, int ny, int nz) +{ + int i1 = amrex::min(i+1, nx-1); + int j1 = amrex::min(j+1, ny-1); + int k1 = amrex::min(k+1, nz-1); + + Real c000 = f[idx(i ,j ,k ,nx,ny)]; + Real c100 = f[idx(i1,j ,k ,nx,ny)]; + Real c010 = f[idx(i ,j1,k ,nx,ny)]; + Real c110 = f[idx(i1,j1,k ,nx,ny)]; + Real c001 = f[idx(i ,j ,k1,nx,ny)]; + Real c101 = f[idx(i1,j ,k1,nx,ny)]; + Real c011 = f[idx(i ,j1,k1,nx,ny)]; + Real c111 = f[idx(i1,j1,k1,nx,ny)]; + + Real c00 = c000*(1-tx) + c100*tx; + Real c10 = c010*(1-tx) + c110*tx; + Real c01 = c001*(1-tx) + c101*tx; + Real c11 = c011*(1-tx) + c111*tx; + + Real c0 = c00*(1-ty) + c10*ty; + Real c1 = c01*(1-ty) + c11*ty; + + return c0*(1-tz) + c1*tz; +} + +void +InterpolateToFineMF( + const Vector& data_rho, + const Vector& data_theta, + const Vector& data_xvel, + const Vector& data_yvel, + const Vector& data_zvel, + int nx, int ny, int nz, + const std::array& problo, + const std::array& probhi, + MultiFab& mf_fine, + const Geometry& geom_fine) +{ + // coarse spacing + Real dx_c[3]; + dx_c[0] = (probhi[0] - problo[0]) / nx; + dx_c[1] = (probhi[1] - problo[1]) / ny; + dx_c[2] = (probhi[2] - problo[2]) / nz; + + const auto problo_f = geom_fine.ProbLoArray(); + const auto dx_f = geom_fine.CellSizeArray(); + + // Step 1: declare device vectors with correct size + amrex::Gpu::DeviceVector d_rho(data_rho.size()); + amrex::Gpu::DeviceVector d_theta(data_theta.size()); + amrex::Gpu::DeviceVector d_xvel(data_xvel.size()); + amrex::Gpu::DeviceVector d_yvel(data_yvel.size()); + amrex::Gpu::DeviceVector d_zvel(data_zvel.size()); + + // Step 2: copy data from host to device + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + data_rho.begin(), data_rho.end(), + d_rho.begin()); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + data_theta.begin(), data_theta.end(), + d_theta.begin()); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + data_xvel.begin(), data_xvel.end(), + d_xvel.begin()); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + data_yvel.begin(), data_yvel.end(), + d_yvel.begin()); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + data_zvel.begin(), data_zvel.end(), + d_zvel.begin()); + + const Real* rho_ptr = d_rho.data(); + const Real* theta_ptr = d_theta.data(); + const Real* xvel_ptr = d_xvel.data(); + const Real* yvel_ptr = d_yvel.data(); + const Real* zvel_ptr = d_zvel.data(); + // ------------------------------- + // GPU kernel over MultiFab + // ------------------------------- + for (MFIter mfi(mf_fine); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + auto arr = mf_fine.array(mfi); + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + // physical location (fine cell center) + Real x = problo_f[0] + (i + 0.5) * dx_f[0]; + Real y = problo_f[1] + (j + 0.5) * dx_f[1]; + Real z = problo_f[2] + (k + 0.5) * dx_f[2]; + + // map to coarse index space + Real rx = (x - problo[0]) / dx_c[0] - 0.5; + Real ry = (y - problo[1]) / dx_c[1] - 0.5; + Real rz = (z - problo[2]) / dx_c[2] - 0.5; + + int ic = static_cast(floor(rx)); + int jc = static_cast(floor(ry)); + int kc = static_cast(floor(rz)); + + Real tx = rx - ic; + Real ty = ry - jc; + Real tz = rz - kc; + + // clamp + ic = amrex::max(0, amrex::min(ic, nx-1)); + jc = amrex::max(0, amrex::min(jc, ny-1)); + kc = amrex::max(0, amrex::min(kc, nz-1)); + + //printf("The values are x, y, z, rx, ry, rz = %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", x, y, z, rx, ry, rz); + + // interpolate each component using device trilinear + arr(i,j,k,0) = interp_trilinear(rho_ptr, ic,jc,kc, tx,ty,tz, nx,ny,nz); + arr(i,j,k,1) = interp_trilinear(theta_ptr, ic,jc,kc, tx,ty,tz, nx,ny,nz); + arr(i,j,k,2) = interp_trilinear(xvel_ptr, ic,jc,kc, tx,ty,tz, nx,ny,nz); + arr(i,j,k,3) = interp_trilinear(yvel_ptr, ic,jc,kc, tx,ty,tz, nx,ny,nz); + arr(i,j,k,4) = interp_trilinear(zvel_ptr, ic,jc,kc, tx,ty,tz, nx,ny,nz); + }); + } +} + +void +MakeFinalMultiFabs (const MultiFab& mf_cc_fine, + MultiFab& cons_pert, + MultiFab& xvel_pert, + MultiFab& yvel_pert, + MultiFab& zvel_pert) +{ + + for (MFIter mfi(cons_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + auto const& mf_cc_fine_arr = mf_cc_fine.const_array(mfi); + auto const& cons_pert_arr = cons_pert.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real tmp_rho = mf_cc_fine_arr(i,j,k,0); + Real tmp_theta = mf_cc_fine_arr(i,j,k,1); + cons_pert_arr(i,j,k,Rho_comp) = tmp_rho; + cons_pert_arr(i,j,k,RhoTheta_comp) = tmp_rho*tmp_theta; + }); + } + + // --- X-faces (component 2) --- + for (MFIter mfi(xvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& uface = xvel_pert.array(mfi); + auto const& cc = mf_cc_fine.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + uface(i,j,k) = 0.5 * (cc(i-1,j,k,2) + cc(i,j,k,2)); + }); + } + + // --- Y-faces (component 3) --- + for (MFIter mfi(yvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& vface = yvel_pert.array(mfi); + auto const& cc = mf_cc_fine.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + vface(i,j,k) = 0.5 * (cc(i,j-1,k,3) + cc(i,j,k,3)); + }); + } + + // --- Z-faces (component 4) --- + for (MFIter mfi(zvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& wface = zvel_pert.array(mfi); + //auto const& cc = mf_cc_fine.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + wface(i,j,k) = 0.0;//0.5 * (cc(i,j,k-1,4) + cc(i,j,k,4)); + }); + } +} + +void +AddPertToBckgnd(MultiFab& mf_cc_fine, + const MultiFab& mf_cc_pert) +{ + const int ncomp = mf_cc_fine.nComp(); + + // Optional safety check (recommended) + AMREX_ALWAYS_ASSERT(mf_cc_pert.nComp() == ncomp); + + for (MFIter mfi(mf_cc_fine, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + auto const& bg = mf_cc_fine.array(mfi); + auto const& pert = mf_cc_pert.const_array(mfi); + + amrex::ParallelFor(bx, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + Real ens_amp = 0.02*std::abs(bg(i,j,k,n)); + bg(i,j,k,n) += ens_amp*pert(i,j,k,n); + }); + } +} + + +#include +namespace fs = std::filesystem; + +std::string +get_last_plotfile(const std::string& plotfile_dir) +{ + std::vector pltfiles; + + for (const auto& entry : fs::directory_iterator(plotfile_dir)) { + if (entry.is_directory()) { + std::string name = entry.path().filename().string(); + if (name.find("plt") == 0) { + pltfiles.push_back(entry.path().string()); + } + } + } + + std::sort(pltfiles.begin(), pltfiles.end()); + return pltfiles.back(); // last one +} + + +// Reads the plotfile data into cell cenetred multifab +// Does not fill ghost cells +void +read_plot_file(PlotFileData& pf, + const std::vector varnames, + MultiFab& mf) +{ + // ------------------------------------------------------------ + // Open plotfile + // ------------------------------------------------------------ + const std::vector& var_names_pf = pf.varNames(); + + // ------------------------------------------------------------ + // Validate requested variables + // ------------------------------------------------------------ + for (auto const& v : varnames) { + bool found = false; + for (auto const& vpf : var_names_pf) { + if (v == vpf) { + found = true; + break; + } + } + if (!found) { + Abort("read_plot_file: invalid variable name: " + v); + } + } + + // ------------------------------------------------------------ + // Define destination MultiFab (single level only) + // ------------------------------------------------------------ + const int level = 0; + + BoxArray ba = pf.boxArray(level); + DistributionMapping dm(ba); + + int ncomp = varnames.size(); + + mf.define(ba, dm, ncomp, 0); + + // ------------------------------------------------------------ + // Copy plotfile data → mf + // ------------------------------------------------------------ + for (int comp = 0; comp < ncomp; ++comp) + { + const MultiFab& src = pf.get(level, varnames[comp]); + MultiFab::Copy(mf, src, 0, comp, 1, 0); + } +} + +void +ERF::ComputeAndWriteEnsemblePerturbations() +{ + int Nens = solverChoice.n_ensemble; + Vector varnames = {"density","theta", "x_velocity","y_velocity","z_velocity"}; + const std::string member_prefix = "member_"; + // ------------------------------- + // Step 1: Read last plotfile for each member + // ------------------------------- + Vector> mf_ens(Nens); + std::vector pltfiles; + + // Step 1: find all plotfiles (assuming all members have same set) + { + std::string pf_dir = member_prefix + "00/plotfiles"; // member 0 as reference + for (const auto& entry : fs::directory_iterator(pf_dir)) { + if (entry.is_directory()) { + std::string name = entry.path().filename().string(); + if (name.find("plt") == 0) { + pltfiles.push_back(name); // store just the directory name + } + } + } + std::sort(pltfiles.begin(), pltfiles.end()); + } + + +// Step 2: loop over all plotfiles (timesteps) +for (const auto& pf_name : pltfiles) +{ + for (int n = 0; n < Nens; ++n) + { + std::string member_dir = member_prefix + amrex::Concatenate("", n, 2); + std::string pf_path = member_dir + "/plotfiles/" + pf_name; + + PlotFileData pf(pf_path); + + const BoxArray& ba = pf.boxArray(0); + const DistributionMapping& dm = pf.DistributionMap(0); + int ncomp = varnames.size(); + + mf_ens[n] = std::make_unique(ba, dm, ncomp, 0); + + read_plot_file(pf, varnames, *mf_ens[n]); + } + + // ------------------------------- + // Step 2: Compute ensemble mean + // ------------------------------- + MultiFab mf_mean(mf_ens[0]->boxArray(), + mf_ens[0]->DistributionMap(), + mf_ens[0]->nComp(), + mf_ens[0]->nGrow()); + + mf_mean.setVal(0.0); + + for (int n = 0; n < Nens; ++n) { + MultiFab::Add(mf_mean, *mf_ens[n], 0, 0, mf_mean.nComp(), mf_mean.nGrow()); + } + + mf_mean.mult(1.0 / Real(Nens), 0, mf_mean.nComp(), mf_mean.nGrow()); + + // ------------------------------- + // Step 3 & 4: Compute perturbations and write plotfiles + // ------------------------------- + for (int n = 0; n < Nens; ++n) + { + MultiFab mf_pert(mf_mean.boxArray(), + mf_mean.DistributionMap(), + mf_mean.nComp(), + mf_mean.nGrow()); + + // perturbation = member - mean + MultiFab::Copy(mf_pert, *mf_ens[n], 0, 0, mf_mean.nComp(), mf_mean.nGrow()); + MultiFab::Subtract(mf_pert, mf_mean, 0, 0, mf_mean.nComp(), mf_mean.nGrow()); + + // Create output directory + std::string member_dir = member_prefix + amrex::Concatenate("", n, 2); + std::string out_dir = member_dir + "/pertfiles"; + fs::create_directories(out_dir); + + std::string pf_path = member_dir + "/plotfiles/" + pf_name; + // Write perturbation plotfile + std::string last_pf_name = fs::path(pf_path).filename().string(); // "pltfile00020" + + + // Extract numeric suffix (everything after "pltfile") + std::string suffix = last_pf_name.substr(3); // "00020" + + // Construct perturbation plotfile name + std::string pltname = out_dir + "/plt_pert_" + suffix; + WriteSingleLevelPlotfile(pltname, + mf_pert, + varnames, + geom[0], + 0.0, // time + 0); // level + } +} +} + +void +ERF::create_background_state_for_ensemble (int lev, + MultiFab& mf_cc_pert, + MultiFab& cons_pert, + MultiFab& xvel_pert, + MultiFab& yvel_pert, + MultiFab& zvel_pert) +{ + + ignore_unused(lev); + int nx_crse, ny_crse, nz_crse, ng_crse, ncomp_crse; + Vector> data_crse; + std::array problo_ext, probhi_ext; + + Vector data_rho, data_theta, data_xvel, data_yvel, data_zvel; + + ReadCustomDataFile(solverChoice.coarse_bckgnd_data_file, + nx_crse, ny_crse, nz_crse, ng_crse, ncomp_crse, + problo_ext, probhi_ext, + data_rho, data_theta, data_xvel, data_yvel, data_zvel); + + Geometry& geom_fine = geom[0]; + // Create a cell-centered multifab on the fine mesh - ie. something with the same boxarray, + // distributed mapping, nGrow, but with 5 components + MultiFab mf_cc_fine; + const MultiFab& src = vars_new[0][0]; + int ncomp = 5; + mf_cc_fine.define(src.boxArray(), src.DistributionMap(), + ncomp, src.nGrow()); + + InterpolateToFineMF(data_rho, data_theta, data_xvel, data_yvel, data_zvel, + nx_crse, ny_crse, nz_crse, + problo_ext, probhi_ext, + mf_cc_fine, + geom_fine); + + ApplyNeumannBCs(geom_fine, mf_cc_fine); + + Vector varnames = {"density","theta", "x_velocity","y_velocity","z_velocity"}; + + // Add pertubrations stored in the "pert" variables in the function arguments + // (multiplied by the corresponding amplitude) + AddPertToBckgnd(mf_cc_fine, mf_cc_pert); + ApplyNeumannBCs(geom_fine, mf_cc_fine); + //WriteSingleLevelPlotfile("1_plt_final", mf_cc_fine, varnames, geom_fine, 0.0, 0); + + MakeFinalMultiFabs(mf_cc_fine, cons_pert, xvel_pert, yvel_pert, zvel_pert); +} diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 4aef3e9a0..10c2ff9c4 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -9,6 +9,7 @@ CEXE_sources += ERF_InitSponge.cpp CEXE_sources += ERF_Init1D.cpp CEXE_sources += ERF_InitTurbPert.cpp CEXE_sources += ERF_InitImmersedForcing.cpp +CEXE_sources += ERF_InitForEnsemble.cpp ifeq ($(USE_WINDFARM),TRUE) CEXE_sources += ERF_InitWindFarm.cpp diff --git a/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H b/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H new file mode 100644 index 000000000..90a2b948f --- /dev/null +++ b/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H @@ -0,0 +1,49 @@ + + // Parse params + ParmParse pp_for_pert_vels("prob"); + + amrex::Real M_inf = 0.2; pp_for_pert_vels.query("M_inf", M_inf);//freestream Mach number [-] + amrex::Real T_inf = 300.0; pp_for_pert_vels.query("T_inf", T_inf);//freestream temperature [K] + amrex::Real alpha = 0.0; pp_for_pert_vels.query("alpha", alpha);//inflow angle, 0 --> x-aligned [rad] + amrex::Real gamma = Gamma; pp_for_pert_vels.query("gamma", gamma);//specific heat ratio [-] + amrex::Real beta = 0.01; pp_for_pert_vels.query("beta", beta);//non-dimensional max perturbation strength [-] + amrex::Real sigma = 1.0; pp_for_pert_vels.query("sigma", sigma);//Gaussian standard deviation, i.e., spreading parameter [-] + amrex::Real R = 2.0; pp_for_pert_vels.query("R", R);//characteristic length scale for grid [m] + amrex::Real xc_frac = 0.5; pp_for_pert_vels.query("xc", xc_frac);//normalized x-location of vortex center [-] + amrex::Real yc_frac = 0.5; pp_for_pert_vels.query("yc", yc_frac);//normalized y-location of vortex center [-] + + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* prob_hi = geomdata.ProbHi(); + + amrex::Real xc = xc_frac * (prob_hi[0] - prob_lo[0]); + amrex::Real yc = yc_frac * (prob_hi[1] - prob_lo[1]); + + amrex::Real a_inf = std::sqrt(gamma * R_d * T_inf); + + // Set the x-velocity + ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real* dx = geomdata.CellSize(); + + // Note that since we only use (x-xc) and (y-yc), we neglect the prob_lo offset in these + const Real x = i * dx[0]; // face center + const Real y = (j + 0.5) * dx[1]; // cell center + const Real Omg = erf_vortex_Gaussian(x,y,xc,yc,R,beta,sigma); + + x_vel_pert(i, j, k) = (M_inf * std::cos(alpha) + - (y - yc)/R * Omg) * a_inf; + }); + + // Set the y-velocity + ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real* dx = geomdata.CellSize(); + + // Note that since we only use (x-xc) and (y-yc), we neglect the prob_lo offset in these + const Real x = (i + 0.5) * dx[0]; // cell center + const Real y = j * dx[1]; // face center + const Real Omg = erf_vortex_Gaussian(x,y,xc,yc,R,beta,sigma); + + y_vel_pert(i, j, k) = (M_inf * std::sin(alpha) + + (x - xc)/R * Omg) * a_inf; + }); diff --git a/Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H b/Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H new file mode 100644 index 000000000..7300ef74b --- /dev/null +++ b/Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H @@ -0,0 +1,77 @@ + + // Parse params + ParmParse pp("prob"); + + amrex::Real p_inf = p_0; pp.query("p_inf", p_inf); //freestream pressure [Pa] + amrex::Real T_inf = 300.0; pp.query("T_inf", T_inf);//freestream temperature [K] + amrex::Real M_inf = 0.2; pp.query("M_inf", M_inf);//freestream Mach number [-] + amrex::Real alpha = 0.0; pp.query("alpha", alpha);//inflow angle, 0 --> x-aligned [rad] + amrex::Real gamma = Gamma; pp.query("gamma", gamma);//specific heat ratio [-] + amrex::Real beta = 0.01; pp.query("beta", beta);//non-dimensional max perturbation strength [-] + amrex::Real sigma = 1.0; pp.query("sigma", sigma);//Gaussian standard deviation, i.e., spreading parameter [-] + amrex::Real R = 2.0; pp.query("R", R);//characteristic length scale for grid [m] + amrex::Real xc_frac = 0.5; pp.query("xc", xc_frac);//normalized x-location of vortex center [-] + amrex::Real yc_frac = 0.5; pp.query("yc", yc_frac);//normalized y-location of vortex center [-] + + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* prob_hi = geomdata.ProbHi(); + + amrex::Real xc = xc_frac * (prob_hi[0] - prob_lo[0]); + amrex::Real yc = yc_frac * (prob_hi[1] - prob_lo[1]); + + // amrex::Print() << " vortex initialized at (" + // << parms.xc << ", " + // << parms.yc << ")" + // << std::endl; + + amrex::Real inv_gm1 = 1.0 / (gamma - 1.0); + + // amrex::Print() << " reference pressure = " << parms.p_inf << " Pa" << std::endl; + // amrex::Print() << " reference temperature = " << parms.T_inf << " K" << std::endl; + // amrex::Print() << " reference potential temperature (not used) = " << parms.T_0 << " K" << std::endl; + + amrex::Real rho_0 = p_inf / (R_d * T_inf); + // amrex::Print() << " calculated freestream air density = " + // << parms.rho_0 << " kg/m^3" + // << std::endl; + + amrex::Real a_inf = std::sqrt(gamma * R_d * T_inf); + // amrex::Print() << " calculated speed of sound, a = " + // << parms.a_inf << " m/s" + // << std::endl; + + // amrex::Print() << " freestream u/a = " + // << parms.M_inf * std::cos(parms.alpha) + // << std::endl; + // amrex::Print() << " freestream v/a = " + // << parms.M_inf * std::sin(parms.alpha) + // << std::endl; + + const amrex::Real rdOcp = sc.rdOcp; + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real* dx = geomdata.CellSize(); + + // Note that since we only use (x-xc) and (y-yc), we neglect the prob_lo offset in these + const Real x = (i + 0.5) * dx[0]; // cell center + const Real y = (j + 0.5) * dx[1]; // cell center + + // Calculate perturbation temperature + const Real Omg = erf_vortex_Gaussian(x,y,xc,yc,R,beta,sigma); + const Real deltaT = -(gamma - 1.0)/(2.0*sigma*sigma) * Omg*Omg; + + // Set the perturbation density + const Real rho_norm = std::pow(1.0 + deltaT, inv_gm1); + //state_pert(i, j, k, Rho_comp) = rho_norm * rho_0 - r_hse(i,j,k); + + // Initial _potential_ temperature + const Real T = (1.0 + deltaT) * T_inf; + const Real p = std::pow(rho_norm, Gamma) / Gamma // isentropic relation + * rho_0*a_inf*a_inf; + const Real rho_theta = rho_0 * rho_norm * (T * std::pow(p_0 / p, rdOcp)); // T --> theta + state_pert(i, j, k, RhoTheta_comp) = rho_theta - getRhoThetagivenP(p_hse(i,j,k)); // Set the perturbation rho*theta + + const Real r2d_xy = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc)); + state_pert(i, j, k, RhoScalar_comp) = 0.25 * (1.0 + std::cos(PI * std::min(r2d_xy, R) / R)); + });