From a7dffefc31fff25f4e37fb0dc951d71cb5704e8f Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 23 Feb 2026 18:58:07 -0800 Subject: [PATCH 01/23] Reading plot file and interpolate directory --- .../ReadPltFileAndInterpolate/Makefile | 33 +++++ .../ReadPltFileAndInterpolate/ReadPlotFile.H | 23 ++++ .../ReadPlotFile.cpp | 122 ++++++++++++++++++ .../ReadPltFileAndInterpolate/main.cpp | 39 ++++++ 4 files changed, 217 insertions(+) create mode 100644 Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile create mode 100644 Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H create mode 100644 Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp create mode 100644 Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile new file mode 100644 index 0000000000..9f2c8cc9ac --- /dev/null +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile @@ -0,0 +1,33 @@ +#--------------------------- +# 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 + +SRCS = main.cpp ReadPlotFile.cpp +HDRS = ReadPlotFile.H + +OBJS = $(SRCS:.cpp=.o) +EXE = out + +# Add math, dl, pthread as required by AMReX +LIBS = -lamrex -lpthread -lgfortran -ldl -lm + +all: $(EXE) + +$(EXE): $(OBJS) + $(CXX) $(CXXFLAGS) -o $@ $^ -L$(AMREX_LIB) $(LIBS) + +# Compile rule with header dependency +%.o: %.cpp $(HDRS) + $(CXX) $(CXXFLAGS) -I$(AMREX_INC) -c $< -o $@ + +clean: + rm -f $(OBJS) $(EXE) + +.PHONY: all clean diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H new file mode 100644 index 0000000000..1711f48ac9 --- /dev/null +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -0,0 +1,23 @@ +#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); + +#endif diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp new file mode 100644 index 0000000000..8b093cbecc --- /dev/null +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -0,0 +1,122 @@ +#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; +} + +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 = 0; + + 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); + } +} diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp new file mode 100644 index 0000000000..e608029f6e --- /dev/null +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -0,0 +1,39 @@ +#include +#include +#include + +#include "ReadPlotFile.H" + +using namespace amrex; + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + + { + amrex::MultiFab mf; + + PlotFileData pf("plt000000"); + + ReadPlotFile("vars.txt", pf, mf); + + amrex::Print() << "MultiFab has " + << mf.nComp() << " components\n"; + + Array is_periodic{AMREX_D_DECL(0,0,0)}; + + Geometry geom(pf.probDomain(0), + RealBox(pf.probLo(), pf.probHi()), + pf.coordSys(), + is_periodic); + + // Read variable names + Vector varnames = ReadVarNames("vars.txt"); + + // Write plotfile + WriteSingleLevelPlotfile("plt_new", mf, varnames, geom, 0.0, 0); + + } + amrex::Finalize(); +} + From 3fbabec228ecca10df634da9175ccf8d7a207c41 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 23 Feb 2026 19:06:36 -0800 Subject: [PATCH 02/23] Updating routines --- .../ReadPltFileAndInterpolate/main.cpp | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp index e608029f6e..58e600dd71 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -8,32 +8,33 @@ using namespace amrex; int main (int argc, char* argv[]) { - amrex::Initialize(argc, argv); + Initialize(argc, argv); { - amrex::MultiFab mf; + MultiFab mf_coarse; + MultiFab mf_fine; - PlotFileData pf("plt000000"); + PlotFileData pf_coarse("plotfile_coarse"); - ReadPlotFile("vars.txt", pf, mf); + ReadPlotFile("vars.txt", pf_coarse, mf_coarse); - amrex::Print() << "MultiFab has " - << mf.nComp() << " components\n"; + PlotFileData pf_fine("plotfile_fine"); + ReadPlotFile("vars.txt", pf_fine, mf_fine); Array is_periodic{AMREX_D_DECL(0,0,0)}; - Geometry geom(pf.probDomain(0), - RealBox(pf.probLo(), pf.probHi()), - pf.coordSys(), + Geometry geom_coarse(pf_coarse.probDomain(0), + RealBox(pf_coarse.probLo(), pf_coarse.probHi()), + pf_coarse.coordSys(), is_periodic); // Read variable names Vector varnames = ReadVarNames("vars.txt"); // Write plotfile - WriteSingleLevelPlotfile("plt_new", mf, varnames, geom, 0.0, 0); + WriteSingleLevelPlotfile("plt_new", mf_coarse, varnames, geom_coarse, 0.0, 0); } - amrex::Finalize(); + Finalize(); } From e33c79402bbf1b8860c436718c2178424a42b7c0 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Tue, 24 Feb 2026 15:34:09 -0800 Subject: [PATCH 03/23] Verifying workflow --- .../ReadPltFileAndInterpolate/ReadPlotFile.H | 16 ++ .../ReadPlotFile.cpp | 248 +++++++++++++++++- .../ReadPltFileAndInterpolate/main.cpp | 38 ++- 3 files changed, 292 insertions(+), 10 deletions(-) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H index 1711f48ac9..9bcc738393 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -20,4 +20,20 @@ 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 + #endif diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp index 8b093cbecc..0c5fca4180 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -46,6 +46,8 @@ ReadVarNames(const std::string& 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, @@ -107,7 +109,7 @@ ReadPlotFile(const std::string& var_filename, DistributionMapping dm(ba); int ncomp = varnames.size(); - int ngrow = 0; + int ngrow = 1; mf.define(ba, dm, ncomp, ngrow); @@ -120,3 +122,247 @@ ReadPlotFile(const std::string& var_filename, MultiFab::Copy(mf, src, 0, comp, 1, 0); } } + +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 +} + +/*IntVect +find_bound_idx(Real x, Real y, Real z, + const Geometry& geom, + BoundType bound_type) +{ + const auto prob_lo = geom.ProbLoArray(); + const auto dx = geom.CellSizeArray(); + constexpr Real eps = 1.e-12; + + auto map = [&](Real x, int d) { + Real r = (x - prob_lo[d]) / dx[d]; + return (bound_type == BoundType::Lo) + ? static_cast(std::floor(r + eps)) + : static_cast(std::floor(r - eps)); + }; + + IntVect idx(map(x,0), map(y,1), map(z,2)); + + const Box& domain = geom.Domain(); + idx = amrex::max(idx, domain.smallEnd()); + idx = amrex::min(idx, domain.bigEnd()); + + return idx; +}*/ + + +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); + } + + 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); +} diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp index 58e600dd71..8e0d51d7aa 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -11,15 +11,19 @@ int main (int argc, char* argv[]) Initialize(argc, argv); { - MultiFab mf_coarse; - MultiFab mf_fine; + MultiFab mf_cc_coarse; + MultiFab mf_cc_fine; - PlotFileData pf_coarse("plotfile_coarse"); + PlotFileData pf_coarse("plt_coarse"); - ReadPlotFile("vars.txt", pf_coarse, mf_coarse); + ReadPlotFile("vars.txt", pf_coarse, mf_cc_coarse); - PlotFileData pf_fine("plotfile_fine"); - ReadPlotFile("vars.txt", pf_fine, mf_fine); + PlotFileData pf_fine("plt_fine"); + + // Read variable names + Vector varnames = ReadVarNames("vars.txt"); + + ReadPlotFile("vars.txt", pf_fine, mf_cc_fine); Array is_periodic{AMREX_D_DECL(0,0,0)}; @@ -28,11 +32,27 @@ int main (int argc, char* argv[]) pf_coarse.coordSys(), is_periodic); - // Read variable names - Vector varnames = ReadVarNames("vars.txt"); + 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); + 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); + + CreateCellCenteredMultiFabFromNodalMultiFab(mf_cc_tmp, coarse_multifab_on_fine_dmap); // Write plotfile - WriteSingleLevelPlotfile("plt_new", mf_coarse, varnames, geom_coarse, 0.0, 0); + WriteSingleLevelPlotfile("plt_2", mf_cc_tmp, varnames, geom_coarse, 0.0, 0); } Finalize(); From 55baf13f7ae1a70dcbcca86b13e03f94e23a2690 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Tue, 3 Mar 2026 12:51:48 -0800 Subject: [PATCH 04/23] Updating with some minor changes --- .../ReadPlotFile.cpp | 31 ++++--------------- .../ReadPltFileAndInterpolate/main.cpp | 8 ++--- 2 files changed, 9 insertions(+), 30 deletions(-) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp index 0c5fca4180..ad0a913a01 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -293,31 +293,6 @@ for (const auto& b : bl_coarse) { return IntVect::TheZeroVector(); // unreachable if Abort } -/*IntVect -find_bound_idx(Real x, Real y, Real z, - const Geometry& geom, - BoundType bound_type) -{ - const auto prob_lo = geom.ProbLoArray(); - const auto dx = geom.CellSizeArray(); - constexpr Real eps = 1.e-12; - - auto map = [&](Real x, int d) { - Real r = (x - prob_lo[d]) / dx[d]; - return (bound_type == BoundType::Lo) - ? static_cast(std::floor(r + eps)) - : static_cast(std::floor(r - eps)); - }; - - IntVect idx(map(x,0), map(y,1), map(z,2)); - - const Box& domain = geom.Domain(); - idx = amrex::max(idx, domain.smallEnd()); - idx = amrex::min(idx, domain.bigEnd()); - - return idx; -}*/ - void GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, @@ -359,6 +334,12 @@ GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, 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)); diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp index 8e0d51d7aa..40ad4e0e7c 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -15,14 +15,9 @@ int main (int argc, char* argv[]) MultiFab mf_cc_fine; PlotFileData pf_coarse("plt_coarse"); - ReadPlotFile("vars.txt", pf_coarse, mf_cc_coarse); PlotFileData pf_fine("plt_fine"); - - // Read variable names - Vector varnames = ReadVarNames("vars.txt"); - ReadPlotFile("vars.txt", pf_fine, mf_cc_fine); Array is_periodic{AMREX_D_DECL(0,0,0)}; @@ -42,6 +37,9 @@ int main (int argc, char* argv[]) 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; From 9f148d7ff413a698e54255ebd0bcd75aaafe35e7 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Wed, 4 Mar 2026 15:42:39 -0800 Subject: [PATCH 05/23] Correcting some errors in ERF_Prob.cpp --- Exec/DevTests/DataAssimilation/ERF_Prob.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Exec/DevTests/DataAssimilation/ERF_Prob.cpp b/Exec/DevTests/DataAssimilation/ERF_Prob.cpp index 8140346813..32bf66dbc8 100644 --- a/Exec/DevTests/DataAssimilation/ERF_Prob.cpp +++ b/Exec/DevTests/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 From 509c5f38992fdf434e1821d077f58fa7f2428e1d Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Wed, 4 Mar 2026 15:48:55 -0800 Subject: [PATCH 06/23] Interpolation works now --- .../ReadPltFileAndInterpolate/ReadPlotFile.H | 13 +++ .../ReadPlotFile.cpp | 93 +++++++++++++++++++ .../ReadPltFileAndInterpolate/main.cpp | 11 +++ 3 files changed, 117 insertions(+) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H index 9bcc738393..c64a4f1d43 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -36,4 +36,17 @@ 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); + #endif diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp index ad0a913a01..3fccc0f640 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -347,3 +347,96 @@ GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, 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; + }); + } +} diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp index 40ad4e0e7c..833edf7dbb 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -47,11 +47,22 @@ int main (int argc, char* argv[]) 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_from_coarse; + PopulateFineCellCenteredFromCoarseNodal(geom_coarse, geom_fine, coarse_multifab_on_fine_dmap, + mf_cc_fine, mf_cc_from_coarse); + + WriteSingleLevelPlotfile("plt_final", mf_cc_from_coarse, varnames, geom_fine, 0.0, 0); + } Finalize(); } From 87979bf9186c41e09d902a335da970aa2231e7bf Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Wed, 4 Mar 2026 15:53:59 -0800 Subject: [PATCH 07/23] Correcting inputs for coarse --- CMake/BuildERFExe.cmake | 1 + .../DataAssimilation/inputs_advecting_coarse | 4 +- .../ERF_InitCustomPertState.cpp | 100 ---------------- Source/Initialization/ERF_InitForEnsemble.cpp | 108 ++++++++++++++++++ Source/Initialization/Make.package | 1 + 5 files changed, 112 insertions(+), 102 deletions(-) create mode 100644 Source/Initialization/ERF_InitForEnsemble.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 54f7fcd3cf..341236c882 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -307,6 +307,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/DevTests/DataAssimilation/inputs_advecting_coarse b/Exec/DevTests/DataAssimilation/inputs_advecting_coarse index cbf8e27170..af3ff9bf31 100644 --- a/Exec/DevTests/DataAssimilation/inputs_advecting_coarse +++ b/Exec/DevTests/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/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 7b6af11648..8f4e8fa369 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -110,103 +110,3 @@ 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(3.0 * sigma / dmesh); // stencil radius - - // ---- Precompute Gaussian weights on host ---- - 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 = 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(); - - // 1. Define ngrow_big using the actual dimension macro - IntVect ngrow_big(AMREX_D_DECL(r, r, 0)); - - // 2. 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); - - // 3. 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 = 0.0; - 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; - }); - } -} diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp new file mode 100644 index 0000000000..784a683fe5 --- /dev/null +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -0,0 +1,108 @@ +/** + * \file ERF_InitForEnsemble.cpp + */ + +#include +#include + +using namespace amrex; + +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(3.0 * sigma / dmesh); // stencil radius + + // ---- Precompute Gaussian weights on host ---- + 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 = 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(); + + // 1. Define ngrow_big using the actual dimension macro + IntVect ngrow_big(AMREX_D_DECL(r, r, 0)); + + // 2. 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); + + // 3. 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 = 0.0; + 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; + }); + } +} diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 4aef3e9a0c..10c2ff9c48 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 From d7c3e674c3f6fda810a8bcde8c7e134058ce4fb3 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Wed, 4 Mar 2026 15:57:08 -0800 Subject: [PATCH 08/23] Correcting trailing whitespaces and tabs --- .../ReadPltFileAndInterpolate/ReadPlotFile.H | 2 +- .../ReadPltFileAndInterpolate/ReadPlotFile.cpp | 8 ++++---- .../DataAssimilation/ReadPltFileAndInterpolate/main.cpp | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H index c64a4f1d43..1d569d6046 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -46,7 +46,7 @@ PopulateFineCellCenteredFromCoarseNodal(const amrex::Geometry& geom_coarse, const amrex::MultiFab& mf_cc_fine, amrex::MultiFab& mf_cc_from_coarse); -void +void check_large_values(const amrex::MultiFab& mf); #endif diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp index 3fccc0f640..e7604bd10d 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -17,7 +17,7 @@ using namespace amrex; // ------------------------------------------------------------ // Read variable names from a file // ------------------------------------------------------------ -Vector +Vector ReadVarNames(const std::string& filename) { Vector varnames; @@ -294,7 +294,7 @@ for (const auto& b : bl_coarse) { } -void +void GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, const Geometry& geom_fine, const MultiFab& mf_nc_coarse, @@ -324,7 +324,7 @@ GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, 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]; @@ -339,7 +339,7 @@ GetCoarseMultiFabOnFineDMap(const Geometry& geom_coarse, 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)); diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp index 833edf7dbb..448b8acba8 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -43,7 +43,7 @@ int main (int argc, char* argv[]) WriteSingleLevelPlotfile("plt_1", mf_cc_tmp, varnames, geom_coarse, 0.0, 0); MultiFab coarse_multifab_on_fine_dmap; - GetCoarseMultiFabOnFineDMap(geom_coarse, geom_fine, + GetCoarseMultiFabOnFineDMap(geom_coarse, geom_fine, mf_nc_coarse, mf_cc_fine, coarse_multifab_on_fine_dmap); @@ -60,7 +60,7 @@ int main (int argc, char* argv[]) MultiFab mf_cc_from_coarse; PopulateFineCellCenteredFromCoarseNodal(geom_coarse, geom_fine, coarse_multifab_on_fine_dmap, mf_cc_fine, mf_cc_from_coarse); - + WriteSingleLevelPlotfile("plt_final", mf_cc_from_coarse, varnames, geom_fine, 0.0, 0); } From b59e46e95680d6243335c3e378bbb9413d7a6bcb Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Sun, 15 Mar 2026 15:37:32 -0700 Subject: [PATCH 09/23] Getting ready for ensemble simulations --- Source/DataStructs/ERF_DataStruct.H | 28 +- Source/Initialization/ERF_InitForEnsemble.cpp | 364 +++++++++++++++++- 2 files changed, 381 insertions(+), 11 deletions(-) diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 0f2f8be487..aac03a2d3f 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -750,14 +750,21 @@ 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 <= 0.0) { + pp.query("is_init_for_ensemble", is_init_for_ensemble); + if(is_init_for_ensemble) { + amrex::ParmParse pp_ens("ensemble"); + pp_ens.query("ens_pert_amplitude", ens_pert_amplitude); + pp_ens.query("ens_pert_correlated_radius", ens_pert_correlated_radius); + 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, " + "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"); } } @@ -1244,7 +1251,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 = -1.0; + bool is_init_for_ensemble = true; + amrex::Real ens_pert_correlated_radius = -1.0; + amrex::Real ens_pert_amplitude = -1.0; + + }; #endif diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp index 784a683fe5..b8e148db34 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -22,9 +22,9 @@ ERF::create_random_perturbations(const int 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 + ParallelForRNG(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const RandomEngine& engine) noexcept { - xvel_pert_arr(i, j, k) = amrex::Random(engine); + xvel_pert_arr(i, j, k) = Random(engine); }); } } @@ -106,3 +106,363 @@ ERF::apply_gaussian_smoothing_to_perturbations(const int lev, }); } } + +// Reads the plotfile data into cell cenetred multifab +// Does not fill ghost cells +void +read_plot_file(PlotFileData& pf, + const std::string varnames, + MultiFab& mf) +{ + // ------------------------------------------------------------ + // 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("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(); + 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 +create_nodal_mf_from_cc_mf (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) ); + }); + } +} + +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 +get_coarse_mf_on_fine_dmap(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 +populate_fine_cc_mf_from_coarse_nodal_mf(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 +ERF::create_background_state_for_ensemble () +{ + bckgnd_state.resize(max_level+1); + for (int lev = 0; lev < max_level+1; ++lev) { + bckgnd_state[lev].resize(vars_new[lev].size()+1); + for (int comp = 0; comp < vars_new[lev].size(); ++comp) { + const MultiFab& src = vars_new[lev][comp]; + bckgnd_state[lev][comp].define(src.boxArray(), src.DistributionMap(), + src.nComp(), src.nGrow()); + } + } + + MultiFab mf_cc_coarse; + PlotFileData pf_coarse(pltfile_bckgnd_coarse); + Vector varnames = {"density"}; + read_plot_file(pf_coarse, varnames, mf_cc_coarse); + + Geometry geom_coarse(pf_coarse.probDomain(0), + RealBox(pf_coarse.probLo(), pf_coarse.probHi()), + pf_coarse.coordSys(), + is_periodic); + + MultiFab mf_nc_coarse; + create_nodal_mf_from_cc_mf(mf_nc_coarse, + mf_cc_coarse, + geom_coarse); + + MultiFab coarse_multifab_on_fine_dmap; + get_coarse_mf_on_fine_dmap(geom_coarse, geom_fine, + mf_nc_coarse, mf_cc_fine, + coarse_multifab_on_fine_dmap); + + MultiFab mf_cc_fine_from_coarse; + populate_fine_cc_mf_from_coarse_nodal_mf(geom_coarse, + geom_fine, + coarse_multifab_on_fine_dmap, + mf_cc_fine, + mf_cc_fine_from_coarse); +} From e399f63e7384bc49ad1bcb49f19354325d81bb47 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Sun, 15 Mar 2026 15:39:33 -0700 Subject: [PATCH 10/23] Getting ready for ensemble simulations --- Source/ERF.H | 4 ++++ Source/Initialization/ERF_InitCustomPertState.cpp | 12 ++++++++---- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/Source/ERF.H b/Source/ERF.H index 7ea4dcf993..33eb4eeef6 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,6 +498,8 @@ public: amrex::Vector>& z_phys_nd, bool regrid_forces_file_read); + void create_background_state_for_ensemble (); + void create_random_perturbations (const int lev, amrex::MultiFab& cons_pert, amrex::MultiFab& xvel_pert, diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 8f4e8fa369..8480ca4573 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -42,10 +42,14 @@ 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) { + // If initializing for enembls simluations, then + // 1. read in the coarse background state from a plotfile + // 2. interpolate the state data onto the current mesh + // 3. Create spatially correlated perturbations + // 4. Add the perturbations to the background state and then populate the "pert variables + + if(solverChoice.is_init_for_ensemble) { + create_background_state (); 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); } From b9b9738f9aeaa34e5b59fbff04066810a4695c63 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Sun, 15 Mar 2026 15:40:07 -0700 Subject: [PATCH 11/23] Getting ready for ensemble simulations --- .../DataAssimilation/ReadPltFileAndInterpolate/main.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp index 448b8acba8..735e4d10d3 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp @@ -57,11 +57,11 @@ int main (int argc, char* argv[]) // Write plotfile WriteSingleLevelPlotfile("plt_2", mf_cc_tmp, varnames, geom_coarse, 0.0, 0); - MultiFab mf_cc_from_coarse; + MultiFab mf_cc_fine_from_coarse; PopulateFineCellCenteredFromCoarseNodal(geom_coarse, geom_fine, coarse_multifab_on_fine_dmap, - mf_cc_fine, mf_cc_from_coarse); + mf_cc_fine, mf_cc_fine_from_coarse); - WriteSingleLevelPlotfile("plt_final", mf_cc_from_coarse, varnames, geom_fine, 0.0, 0); + WriteSingleLevelPlotfile("plt_final", mf_cc_fine_from_coarse, varnames, geom_fine, 0.0, 0); } Finalize(); From 81d223ddc4183566ad7bf85e3b5c7187d14d7c5e Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 16 Mar 2026 15:14:52 -0700 Subject: [PATCH 12/23] Coarse basckground state to fine interpolation --- Source/DataStructs/ERF_DataStruct.H | 2 +- Source/Initialization/ERF_InitCustomPertState.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index aac03a2d3f..2f5ca320e0 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -1251,7 +1251,7 @@ struct SolverChoice { bool io_hurricane_eye_tracker = false; amrex::Real hurricane_eye_latitude = -1e10, hurricane_eye_longitude = -1e10; - bool is_init_for_ensemble = true; + bool is_init_for_ensemble = false; amrex::Real ens_pert_correlated_radius = -1.0; amrex::Real ens_pert_amplitude = -1.0; diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 8480ca4573..bfe26d6700 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -49,7 +49,7 @@ ERF::init_custom (int lev) // 4. Add the perturbations to the background state and then populate the "pert variables if(solverChoice.is_init_for_ensemble) { - create_background_state (); + create_background_state_for_ensemble(); 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); } From 7f1bbec187f9daade3d438140b022788010931c3 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Wed, 25 Mar 2026 17:43:44 -0700 Subject: [PATCH 13/23] Updating with coarse mesh to fine mesh interpolation routines --- .../ReadPltFileAndInterpolate/Makefile | 47 +- .../ReadPltFileAndInterpolate/ReadPlotFile.H | 8 + .../ReadPlotFile.cpp | 112 ++++ .../{main.cpp => main_parallel.cpp} | 0 .../ReadPltFileAndInterpolate/main_serial.cpp | 48 ++ Source/ERF.H | 6 +- .../ERF_InitCustomPertState.cpp | 6 +- Source/Initialization/ERF_InitForEnsemble.cpp | 576 ++++++++++-------- 8 files changed, 519 insertions(+), 284 deletions(-) rename Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/{main.cpp => main_parallel.cpp} (100%) create mode 100644 Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile index 9f2c8cc9ac..20b24cd94c 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile @@ -1,7 +1,8 @@ + + #--------------------------- # Makefile for AMReX project #--------------------------- - CXX = mpicxx CXXFLAGS = -O2 -std=c++17 -fopenmp -Wall @@ -9,25 +10,49 @@ CXXFLAGS = -O2 -std=c++17 -fopenmp -Wall AMREX_INC = /pscratch/sd/n/nataraj2/AIModCon/GitHub/amrex_repos/amrex/include AMREX_LIB = /pscratch/sd/n/nataraj2/AIModCon/GitHub/amrex_repos/amrex/lib -SRCS = main.cpp ReadPlotFile.cpp HDRS = ReadPlotFile.H +LIBS = -lamrex -lpthread -lgfortran -ldl -lm -lstdc++ -OBJS = $(SRCS:.cpp=.o) -EXE = out +# Default: do nothing if 'make' is called without target +all: + @echo "Use 'make serial' or 'make parallel'" -# Add math, dl, pthread as required by AMReX -LIBS = -lamrex -lpthread -lgfortran -ldl -lm +# ------------------------- +# Serial build +# ------------------------- +SRCS_SERIAL = main_serial.cpp ReadPlotFile.cpp +OBJS_SERIAL = $(SRCS_SERIAL:.cpp=.o) +EXE_SERIAL = out_serial -all: $(EXE) +serial: $(EXE_SERIAL) -$(EXE): $(OBJS) +$(EXE_SERIAL): $(OBJS_SERIAL) $(CXX) $(CXXFLAGS) -o $@ $^ -L$(AMREX_LIB) $(LIBS) -# Compile rule with header dependency +# ------------------------- +# 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 $(OBJS) $(EXE) + rm -f *.o out_serial out_parallel + +.PHONY: all serial parallel clean + -.PHONY: all clean diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H index 1d569d6046..94ada1c1c6 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -49,4 +49,12 @@ PopulateFineCellCenteredFromCoarseNodal(const amrex::Geometry& geom_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/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp index e7604bd10d..d53054d33e 100644 --- a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp +++ b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -123,6 +123,52 @@ ReadPlotFile(const std::string& var_filename, } } +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) @@ -440,3 +486,69 @@ PopulateFineCellCenteredFromCoarseNodal(const Geometry& geom_coarse, }); } } + +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/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp similarity index 100% rename from Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main.cpp rename to Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp b/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp new file mode 100644 index 0000000000..563dcb568a --- /dev/null +++ b/Exec/DevTests/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/Source/ERF.H b/Source/ERF.H index 5d1253b980..449ca463eb 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -498,7 +498,11 @@ public: amrex::Vector>& z_phys_nd, bool regrid_forces_file_read); - void create_background_state_for_ensemble (); + void create_background_state_for_ensemble (int lev, + 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, diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 1be91aef61..68e6ebfea0 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -49,9 +49,9 @@ ERF::init_custom (int lev) // 4. Add the perturbations to the background state and then populate the "pert variables if(solverChoice.is_init_for_ensemble) { - create_background_state_for_ensemble(); - 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); + create_background_state_for_ensemble(lev, cons_pert, xvel_pert, yvel_pert, zvel_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 diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp index b8e148db34..646dc4235f 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -4,6 +4,7 @@ #include #include +#include using namespace amrex; @@ -46,7 +47,7 @@ ERF::apply_gaussian_smoothing_to_perturbations(const int lev, const Real dmesh = std::min(dx, dy); // ---- User choices ---- - const Real sigma = solverChoice.pert_correlated_radius; // e.g. 2 km correlation length + const Real sigma = solverChoice.ens_pert_correlated_radius; // e.g. 2 km correlation length const int r = static_cast(3.0 * sigma / dmesh); // stencil radius // ---- Precompute Gaussian weights on host ---- @@ -107,77 +108,11 @@ ERF::apply_gaussian_smoothing_to_perturbations(const int lev, } } -// Reads the plotfile data into cell cenetred multifab -// Does not fill ghost cells -void -read_plot_file(PlotFileData& pf, - const std::string varnames, - MultiFab& mf) +void ApplyNeumannBCs(const Geometry& geom, + MultiFab& mf_cc) { - // ------------------------------------------------------------ - // 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("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(); - 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 -create_nodal_mf_from_cc_mf (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()); @@ -216,218 +151,313 @@ create_nodal_mf_from_cc_mf (MultiFab& mf_nc, // output nodal MF 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) ); - }); - } } -IntVect -find_bound_idx(const Real& x, const Real& y, const Real& z, - const BoxList& bl_coarse, const Geometry& geom_coarse, - BoundType bound_type) + +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) { - 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])); + std::ifstream ifs(filename_custom, std::ios::binary); + if (!ifs.is_open()) { + Abort("Failed to open file for reading"); } - IntVect idx(i, j, k); - - for (const auto& b : bl_coarse) { - if (b.contains(idx)) { - return idx; + // ---------------------------- + // 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; + } } } - amrex::Print() << x << " " << y << " " << z << " " << idx << std::endl; - - amrex::Print() << "Printing BoxList (coarse):\n"; -for (const auto& b : bl_coarse) { - amrex::Print() << b << "\n"; + ifs.close(); } - amrex::Abort("Bound index not found in any box in BoxList!"); - return IntVect::TheZeroVector(); // unreachable if Abort +void +populate_mf_cc_fine_from_mf_cc_coarse (const Geometry& geom_coarse, + const Geometry& geom_fine, + const MultiFab& coarse_mf_cc_on_fine_dmap, + const MultiFab& mf_cc_fine, + MultiFab& mf_cc_from_coarse) +{ } void -get_coarse_mf_on_fine_dmap(const Geometry& geom_coarse, - const Geometry& geom_fine, - const MultiFab& mf_nc_coarse, - const MultiFab& mf_cc_fine, - MultiFab& coarse_multifab_on_fine_dmap) +populate_mf_face_fine_from_mf_cc_coarse(const Geometry& geom_coarse, + const Geometry& geom_fine, + const MultiFab& mf_cc_coarse, + const MultiFab& mf_face_fine, + MultiFab& mf_face_from_coarse, + int dir) // 0=x,1=y,2=z { - 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); +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); +} - b.setSmall(idx_lo); - b.setBig(idx_hi); +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; +} - /*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;*/ +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); + }); } - - 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 -populate_fine_cc_mf_from_coarse_nodal_mf(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) +MakeFaceCenteredVelocities (const MultiFab& mf_cc_fine, + MultiFab& mf_xvel, + MultiFab& mf_yvel, + MultiFab& mf_zvel) { - AMREX_ALWAYS_ASSERT(coarse_multifab_on_fine_dmap.ixType().nodeCentered()); + BL_PROFILE("MakeFaceCenteredVelocities"); - 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); - } + const BoxArray& ba = mf_cc_fine.boxArray(); + const DistributionMapping& dm = mf_cc_fine.DistributionMap(); + int ng = mf_xvel.nGrow(); + + // --- Define face-centered MultiFabs --- + BoxArray ba_x = amrex::convert(ba, IntVect(1,0,0)); + BoxArray ba_y = amrex::convert(ba, IntVect(0,1,0)); + BoxArray ba_z = amrex::convert(ba, IntVect(0,0,1)); - 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()); + mf_xvel.define(ba_x, dm, 1, ng); + mf_yvel.define(ba_y, dm, 1, ng); + mf_zvel.define(ba_z, dm, 1, ng); - const auto prob_lo_coarse = geom_coarse.ProbLoArray(); - const auto dx_coarse = geom_coarse.CellSizeArray(); + // --- X-faces (component 2) --- + for (MFIter mfi(mf_xvel, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& uface = mf_xvel.array(mfi); + auto const& cc = mf_cc_fine.const_array(mfi); - const auto prob_lo_fine = geom_fine.ProbLoArray(); - const auto dx_fine = geom_fine.CellSizeArray(); + 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)); + }); + } - Box nodal_domain = amrex::convert(geom_coarse.Domain(), IntVect::TheNodeVector()); - const auto nd_lo = nodal_domain.smallEnd(); - const auto nd_hi = nodal_domain.bigEnd(); + // --- Y-faces (component 3) --- + for (MFIter mfi(mf_yvel, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& vface = mf_yvel.array(mfi); + auto const& cc = mf_cc_fine.const_array(mfi); - int ncomp = mf_cc_from_coarse.nComp(); + 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)); + }); + } - for (MFIter mfi(mf_cc_from_coarse, TilingIfNotGPU()); mfi.isValid(); ++mfi) + // --- Z-faces (component 4) --- + for (MFIter mfi(mf_zvel, 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); + const Box& bx = mfi.tilebox(); + auto const& wface = mf_zvel.array(mfi); + auto const& cc = mf_cc_fine.const_array(mfi); - ParallelFor(cbox, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - 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; + wface(i,j,k) = 0.5 * (cc(i,j,k-1,4) + cc(i,j,k,4)); }); } } void -ERF::create_background_state_for_ensemble () +ERF::create_background_state_for_ensemble (int lev, + MultiFab& cons_pert, + MultiFab& xvel_pert, + MultiFab& yvel_pert, + MultiFab& zvel_pert) { bckgnd_state.resize(max_level+1); for (int lev = 0; lev < max_level+1; ++lev) { @@ -439,30 +469,38 @@ ERF::create_background_state_for_ensemble () } } - MultiFab mf_cc_coarse; - PlotFileData pf_coarse(pltfile_bckgnd_coarse); - Vector varnames = {"density"}; - read_plot_file(pf_coarse, varnames, mf_cc_coarse); - - Geometry geom_coarse(pf_coarse.probDomain(0), - RealBox(pf_coarse.probLo(), pf_coarse.probHi()), - pf_coarse.coordSys(), - is_periodic); - - MultiFab mf_nc_coarse; - create_nodal_mf_from_cc_mf(mf_nc_coarse, - mf_cc_coarse, - geom_coarse); - - MultiFab coarse_multifab_on_fine_dmap; - get_coarse_mf_on_fine_dmap(geom_coarse, geom_fine, - mf_nc_coarse, mf_cc_fine, - coarse_multifab_on_fine_dmap); - - MultiFab mf_cc_fine_from_coarse; - populate_fine_cc_mf_from_coarse_nodal_mf(geom_coarse, - geom_fine, - coarse_multifab_on_fine_dmap, - mf_cc_fine, - mf_cc_fine_from_coarse); + std::string filename_custom = "coarse_data.bin"; + 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(filename_custom, + 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 ngrow = src.nGrow(); + 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"}; + WriteSingleLevelPlotfile("plt_final", mf_cc_fine, varnames, geom_fine, 0.0, 0); + + MakeFaceCenteredVelocities(mf_cc_fine, xvel_pert, yvel_pert, zvel_pert); } From 4603ab574e7ce34bc018355c08098190f07e8856 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 14:39:09 -0700 Subject: [PATCH 14/23] Update routines for ensemble simulations --- .Exec_dev/DataAssimilation/main.cpp | 8 +- Source/DataStructs/ERF_DataStruct.H | 9 +- Source/ERF.H | 11 +- .../ERF_InitCustomPertState.cpp | 26 +- Source/Initialization/ERF_InitForEnsemble.cpp | 445 +++++++++++++----- 5 files changed, 368 insertions(+), 131 deletions(-) diff --git a/.Exec_dev/DataAssimilation/main.cpp b/.Exec_dev/DataAssimilation/main.cpp index 6e2d761d19..dbda8d4095 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(".")) { @@ -147,10 +148,12 @@ int main (int argc, char* argv[]) } } // Optional: barrier after move to ensure rank 0 is done - ParallelDescriptor::Barrier(); - + 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/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 84cd2adada..128c64acc9 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -754,8 +754,14 @@ struct SolverChoice { pp.query("is_init_for_ensemble", is_init_for_ensemble); if(is_init_for_ensemble) { amrex::ParmParse pp_ens("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 (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 " @@ -1255,7 +1261,6 @@ struct SolverChoice { bool is_init_for_ensemble = false; 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 449ca463eb..039443ac65 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -499,22 +499,17 @@ public: 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); #ifdef ERF_USE_MULTIBLOCK // constructor used when ERF is created by a multiblock driver diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 68e6ebfea0..e8303fa753 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -42,17 +42,6 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - // If initializing for enembls simluations, then - // 1. read in the coarse background state from a plotfile - // 2. interpolate the state data onto the current mesh - // 3. Create spatially correlated perturbations - // 4. Add the perturbations to the background state and then populate the "pert variables - - if(solverChoice.is_init_for_ensemble) { - create_background_state_for_ensemble(lev, cons_pert, xvel_pert, yvel_pert, zvel_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()) @@ -90,6 +79,7 @@ ERF::init_custom (int lev) solverChoice, 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()); @@ -120,4 +110,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()); } + + // 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 index 646dc4235f..374a5693f0 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -7,105 +7,155 @@ #include using namespace amrex; +namespace fs = std::filesystem; void ERF::create_random_perturbations(const int lev, - MultiFab& cons_pert, - MultiFab& xvel_pert, - MultiFab& yvel_pert, - MultiFab& zvel_pert) + MultiFab& mf_cc_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 RandomEngine& engine) noexcept + 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 { - xvel_pert_arr(i, j, k) = Random(engine); + 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& cons_pert, - MultiFab& xvel_pert, - MultiFab& yvel_pert, - MultiFab& zvel_pert) + MultiFab& mf_cc_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.ens_pert_correlated_radius; // e.g. 2 km correlation length - const int r = static_cast(3.0 * sigma / dmesh); // stencil radius - // ---- Precompute Gaussian weights on host ---- + // ---- 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)); + 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 = v/Z; + v /= Z; } - Gpu::DeviceVector w_dev; - w_dev.resize(w_host.size()); + 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(); - // 1. Define ngrow_big using the actual dimension macro + // ---- Create a grown copy (for stencil access) ---- IntVect ngrow_big(AMREX_D_DECL(r, r, 0)); - // 2. 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); + MultiFab mf_copy(mf_cc_pert.boxArray(), + mf_cc_pert.DistributionMap(), + ncomp, ngrow_big); - // 3. 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()); + mf_copy.ParallelCopy(mf_cc_pert, + 0, 0, ncomp, + IntVect(0), ngrow_big, + gm.periodicity()); - for (MFIter mfi(xvel_pert, TileNoZ()); mfi.isValid(); ++mfi) + // ---- Apply smoothing ---- + for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& tbx = mfi.tilebox(); + const Box& bx = mfi.tilebox(); - auto const& in = xvel_pert_copy.array(mfi); - auto const& out = xvel_pert.array(mfi); + auto const& in = mf_copy.const_array(mfi); + auto const& out = mf_cc_pert.array(mfi); - ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + 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 n = -r; n <= r; ++n) { - Real wij = w[(m+r)*wsize + (n+r)]; - sum += wij * in(i+m, j+n, k); + 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) = sum; + + out(i,j,k,n) = sum; }); } + NormalizeMultiFabRMS_PerComponent(mf_cc_pert); } void ApplyNeumannBCs(const Geometry& geom, @@ -167,7 +217,7 @@ void ReadCustomDataFile(const std::string& filename_custom, { std::ifstream ifs(filename_custom, std::ios::binary); if (!ifs.is_open()) { - Abort("Failed to open file for reading"); + Abort("Failed to open file " + filename_custom + " for reading"); } // ---------------------------- @@ -231,26 +281,6 @@ void ReadCustomDataFile(const std::string& filename_custom, ifs.close(); } -void -populate_mf_cc_fine_from_mf_cc_coarse (const Geometry& geom_coarse, - const Geometry& geom_fine, - const MultiFab& coarse_mf_cc_on_fine_dmap, - const MultiFab& mf_cc_fine, - MultiFab& mf_cc_from_coarse) -{ -} - -void -populate_mf_face_fine_from_mf_cc_coarse(const Geometry& geom_coarse, - const Geometry& geom_fine, - const MultiFab& mf_cc_coarse, - const MultiFab& mf_face_fine, - MultiFab& mf_face_from_coarse, - int dir) // 0=x,1=y,2=z -{ -} - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int idx(int i, int j, int k, int nx, int ny) @@ -392,31 +422,38 @@ InterpolateToFineMF( } void -MakeFaceCenteredVelocities (const MultiFab& mf_cc_fine, - MultiFab& mf_xvel, - MultiFab& mf_yvel, - MultiFab& mf_zvel) +MakeFinalMultiFabs (const MultiFab& mf_cc_fine, + MultiFab& cons_pert, + MultiFab& xvel_pert, + MultiFab& yvel_pert, + MultiFab& zvel_pert) { - BL_PROFILE("MakeFaceCenteredVelocities"); - const BoxArray& ba = mf_cc_fine.boxArray(); - const DistributionMapping& dm = mf_cc_fine.DistributionMap(); - int ng = mf_xvel.nGrow(); + for (MFIter mfi(cons_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); - // --- Define face-centered MultiFabs --- - BoxArray ba_x = amrex::convert(ba, IntVect(1,0,0)); - BoxArray ba_y = amrex::convert(ba, IntVect(0,1,0)); - BoxArray ba_z = amrex::convert(ba, IntVect(0,0,1)); + auto const& mf_cc_fine_arr = mf_cc_fine.const_array(mfi); + auto const& cons_pert_arr = cons_pert.array(mfi); - mf_xvel.define(ba_x, dm, 1, ng); - mf_yvel.define(ba_y, dm, 1, ng); - mf_zvel.define(ba_z, dm, 1, ng); + 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; + }); + } + + const BoxArray& ba = mf_cc_fine.boxArray(); + const DistributionMapping& dm = mf_cc_fine.DistributionMap(); + int ng = xvel_pert.nGrow(); // --- X-faces (component 2) --- - for (MFIter mfi(mf_xvel, TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(xvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - auto const& uface = mf_xvel.array(mfi); + 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) @@ -426,10 +463,10 @@ MakeFaceCenteredVelocities (const MultiFab& mf_cc_fine, } // --- Y-faces (component 3) --- - for (MFIter mfi(mf_yvel, TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(yvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - auto const& vface = mf_yvel.array(mfi); + 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) @@ -439,44 +476,233 @@ MakeFaceCenteredVelocities (const MultiFab& mf_cc_fine, } // --- Z-faces (component 4) --- - for (MFIter mfi(mf_zvel, TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(zvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - auto const& wface = mf_zvel.array(mfi); + 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.5 * (cc(i,j,k-1,4) + cc(i,j,k,4)); + 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) { - bckgnd_state.resize(max_level+1); - for (int lev = 0; lev < max_level+1; ++lev) { - bckgnd_state[lev].resize(vars_new[lev].size()+1); - for (int comp = 0; comp < vars_new[lev].size(); ++comp) { - const MultiFab& src = vars_new[lev][comp]; - bckgnd_state[lev][comp].define(src.boxArray(), src.DistributionMap(), - src.nComp(), src.nGrow()); - } - } - - std::string filename_custom = "coarse_data.bin"; 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(filename_custom, + 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); @@ -500,7 +726,12 @@ ERF::create_background_state_for_ensemble (int lev, ApplyNeumannBCs(geom_fine, mf_cc_fine); Vector varnames = {"density","theta", "x_velocity","y_velocity","z_velocity"}; - WriteSingleLevelPlotfile("plt_final", mf_cc_fine, varnames, geom_fine, 0.0, 0); - MakeFaceCenteredVelocities(mf_cc_fine, xvel_pert, yvel_pert, zvel_pert); + // Add pertubrations stored in the "pert" variables in the function arguments + // (mutliplied 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); } From b8bd13c7e36ab2ad2995dc3d705d358073665679 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 14:40:11 -0700 Subject: [PATCH 15/23] Some minor changes --- Exec/ERF_Prob.cpp | 6 ++++++ Source/DataStructs/ERF_DataStruct.H | 7 +++++++ 2 files changed, 13 insertions(+) diff --git a/Exec/ERF_Prob.cpp b/Exec/ERF_Prob.cpp index fb070df8c5..fe980f162b 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" } amrex::Gpu::streamSynchronize(); @@ -224,6 +227,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" } amrex::Gpu::streamSynchronize(); diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 128c64acc9..38cb576424 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -754,9 +754,14 @@ struct SolverChoice { 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 " @@ -1259,8 +1264,10 @@ struct SolverChoice { amrex::Real hurricane_eye_latitude = -1e10, hurricane_eye_longitude = -1e10; 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 From 44d29150f5f1432c14699c64a847fe4069704813 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 14:40:36 -0700 Subject: [PATCH 16/23] updating ERF.H --- Source/ERF.H | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/ERF.H b/Source/ERF.H index 039443ac65..68c6fe6c0e 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -511,6 +511,8 @@ public: void apply_gaussian_smoothing_to_perturbations (const int lev, amrex::MultiFab& mf_cc_pert); + void ComputeAndWriteEnsemblePerturbations(); + #ifdef ERF_USE_MULTIBLOCK // constructor used when ERF is created by a multiblock driver // calls AmrCore constructor -> AmrMesh constructor From fc43788207977abe973dbaf4848eee6c4800571b Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 19:18:14 -0700 Subject: [PATCH 17/23] Routine for reading plotfile and wrting coarse data --- .../DataAssimilation/ReadPltFileAndInterpolate/Makefile | 0 .../DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H | 0 .../DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp | 0 .../DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp | 0 .../DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp | 0 5 files changed, 0 insertions(+), 0 deletions(-) rename {Exec/DevTests => .Exec_dev}/DataAssimilation/ReadPltFileAndInterpolate/Makefile (100%) rename {Exec/DevTests => .Exec_dev}/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H (100%) rename {Exec/DevTests => .Exec_dev}/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp (100%) rename {Exec/DevTests => .Exec_dev}/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp (100%) rename {Exec/DevTests => .Exec_dev}/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp (100%) diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/Makefile similarity index 100% rename from Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/Makefile rename to .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/Makefile diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H similarity index 100% rename from Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H rename to .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp similarity index 100% rename from Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp rename to .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp similarity index 100% rename from Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp rename to .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_parallel.cpp diff --git a/Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp similarity index 100% rename from Exec/DevTests/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp rename to .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp From 507a60fba466517abb1695c74879a6572b6fa76a Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 19:20:43 -0700 Subject: [PATCH 18/23] Adding header files for data assimilation simulations --- ..._InitCustomPertVels_DataAssimilation_ISV.H | 49 ++++++++++++ .../ERF_InitCustomPert_DataAssimilation_ISV.H | 77 +++++++++++++++++++ 2 files changed, 126 insertions(+) create mode 100644 Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H create mode 100644 Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H diff --git a/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H b/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H new file mode 100644 index 0000000000..7732f4b4e8 --- /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 0000000000..185cf07622 --- /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)); + }); From 4a27a0385dee198f206b3453dba5032b88c3ea99 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 19:24:14 -0700 Subject: [PATCH 19/23] Correcting trailing whitespaces and tabs --- .../ReadPltFileAndInterpolate/ReadPlotFile.H | 4 ++-- .../ReadPltFileAndInterpolate/ReadPlotFile.cpp | 2 +- .../ReadPltFileAndInterpolate/main_serial.cpp | 2 +- .../run_interactive_GPU.sh | 2 ++ .../ReadPltFileAndInterpolate/vars.txt | 5 +++++ .Exec_dev/DataAssimilation/main.cpp | 4 ++-- .../Initialization/ERF_InitCustomPertState.cpp | 6 +++--- Source/Initialization/ERF_InitForEnsemble.cpp | 18 +++++++++--------- 8 files changed, 25 insertions(+), 18 deletions(-) create mode 100644 .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/run_interactive_GPU.sh create mode 100644 .Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/vars.txt diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H index 94ada1c1c6..861203e4b4 100644 --- a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.H @@ -49,9 +49,9 @@ PopulateFineCellCenteredFromCoarseNodal(const amrex::Geometry& geom_coarse, void check_large_values(const amrex::MultiFab& mf); -void +void ApplyNeumannBCs(const amrex::Geometry& geom, - amrex::MultiFab& mf_cc); + amrex::MultiFab& mf_cc); void WriteCustomDataFile(const amrex::Geometry& geom, const amrex::MultiFab& mf_cc, diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp index d53054d33e..e8d6e21012 100644 --- a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/ReadPlotFile.cpp @@ -123,7 +123,7 @@ ReadPlotFile(const std::string& var_filename, } } -void ApplyNeumannBCs(const Geometry& geom, +void ApplyNeumannBCs(const Geometry& geom, MultiFab& mf_cc) { diff --git a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp index 563dcb568a..ad4fc47329 100644 --- a/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp +++ b/.Exec_dev/DataAssimilation/ReadPltFileAndInterpolate/main_serial.cpp @@ -33,7 +33,7 @@ int main (int argc, char* argv[]) is_periodic); - std::string filename_custom = "coarse_data.bin"; + std::string filename_custom = "coarse_data.bin"; ApplyNeumannBCs(geom_coarse, mf_cc_coarse); WriteCustomDataFile(geom_coarse, mf_cc_coarse, filename_custom); 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 0000000000..3374a93d45 --- /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 0000000000..5c98397ff7 --- /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/main.cpp b/.Exec_dev/DataAssimilation/main.cpp index dbda8d4095..370d4136aa 100644 --- a/.Exec_dev/DataAssimilation/main.cpp +++ b/.Exec_dev/DataAssimilation/main.cpp @@ -148,10 +148,10 @@ int main (int argc, char* argv[]) } } // Optional: barrier after move to ensure rank 0 is done - ParallelDescriptor::Barrier(); + ParallelDescriptor::Barrier(); } - ERF tmp_erf; + ERF tmp_erf; tmp_erf.ComputeAndWriteEnsemblePerturbations(); BL_PROFILE_VAR_STOP(pmain); diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp index 413786e5ba..74a7bf6928 100644 --- a/Source/Initialization/ERF_InitCustomPertState.cpp +++ b/Source/Initialization/ERF_InitCustomPertState.cpp @@ -147,15 +147,15 @@ ERF::init_custom (int lev) // If initializing for ensemble simluations, then // 1. Create cell-centered random perturbations - // 2. Create cell-centered spatially correlated 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 + // 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]); + 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 index 374a5693f0..266d999a59 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -409,7 +409,7 @@ InterpolateToFineMF( 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); + //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); @@ -672,7 +672,7 @@ for (const auto& pf_name : pltfiles) // 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" @@ -716,18 +716,18 @@ ERF::create_background_state_for_ensemble (int lev, 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, + + 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 + // Add pertubrations stored in the "pert" variables in the function arguments // (mutliplied by the corresponding amplitude) AddPertToBckgnd(mf_cc_fine, mf_cc_pert); ApplyNeumannBCs(geom_fine, mf_cc_fine); From 3111a62d0f7d687268cbf27bf9e0e39eafe1970a Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 19:25:55 -0700 Subject: [PATCH 20/23] Adding header file --- .Exec_dev/DataAssimilation/ERF_Prob.H | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.Exec_dev/DataAssimilation/ERF_Prob.H b/.Exec_dev/DataAssimilation/ERF_Prob.H index 1203051b7f..8ff4b61266 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: From 9d6313b8766652fce4aa28750c682bd1516bb37c Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Mon, 30 Mar 2026 19:40:20 -0700 Subject: [PATCH 21/23] Correcting spelling error --- Source/Initialization/ERF_InitForEnsemble.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp index 266d999a59..92be33eaf0 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -728,7 +728,7 @@ ERF::create_background_state_for_ensemble (int lev, Vector varnames = {"density","theta", "x_velocity","y_velocity","z_velocity"}; // Add pertubrations stored in the "pert" variables in the function arguments - // (mutliplied by the corresponding amplitude) + // (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); From f305641c18f34b2cc82b137085c4b1e0e05ddfd0 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Tue, 31 Mar 2026 08:28:11 -0700 Subject: [PATCH 22/23] Correcting unused variables errors --- Exec/ERF_Prob.cpp | 2 ++ Source/Initialization/ERF_InitForEnsemble.cpp | 7 ++----- Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H | 8 ++++---- Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H | 2 +- 4 files changed, 9 insertions(+), 10 deletions(-) diff --git a/Exec/ERF_Prob.cpp b/Exec/ERF_Prob.cpp index 030dce9efb..2fe030087b 100644 --- a/Exec/ERF_Prob.cpp +++ b/Exec/ERF_Prob.cpp @@ -142,6 +142,7 @@ Problem::init_custom_pert ( } 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" } @@ -232,6 +233,7 @@ Problem::init_custom_pert_vels ( } 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/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp index 92be33eaf0..316037f6f5 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -445,10 +445,6 @@ MakeFinalMultiFabs (const MultiFab& mf_cc_fine, }); } - const BoxArray& ba = mf_cc_fine.boxArray(); - const DistributionMapping& dm = mf_cc_fine.DistributionMap(); - int ng = xvel_pert.nGrow(); - // --- X-faces (component 2) --- for (MFIter mfi(xvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -696,6 +692,8 @@ ERF::create_background_state_for_ensemble (int lev, 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; @@ -712,7 +710,6 @@ ERF::create_background_state_for_ensemble (int lev, // distributed mapping, nGrow, but with 5 components MultiFab mf_cc_fine; const MultiFab& src = vars_new[0][0]; - int ngrow = src.nGrow(); int ncomp = 5; mf_cc_fine.define(src.boxArray(), src.DistributionMap(), ncomp, src.nGrow()); diff --git a/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H b/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H index 7732f4b4e8..90a2b948f8 100644 --- a/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H +++ b/Source/Prob/ERF_InitCustomPertVels_DataAssimilation_ISV.H @@ -30,8 +30,8 @@ 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; + x_vel_pert(i, j, k) = (M_inf * std::cos(alpha) + - (y - yc)/R * Omg) * a_inf; }); // Set the y-velocity @@ -44,6 +44,6 @@ 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; + 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 index 185cf07622..7300ef74b6 100644 --- a/Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H +++ b/Source/Prob/ERF_InitCustomPert_DataAssimilation_ISV.H @@ -70,7 +70,7 @@ 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 + 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)); From 9b132d45df88c397d55198bcc611f6018ba594c5 Mon Sep 17 00:00:00 2001 From: nataraj2 Date: Tue, 31 Mar 2026 11:05:57 -0700 Subject: [PATCH 23/23] Correcting unused variable error --- Source/Initialization/ERF_InitForEnsemble.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp index 316037f6f5..2d6b358322 100644 --- a/Source/Initialization/ERF_InitForEnsemble.cpp +++ b/Source/Initialization/ERF_InitForEnsemble.cpp @@ -476,7 +476,7 @@ MakeFinalMultiFabs (const MultiFab& mf_cc_fine, { const Box& bx = mfi.tilebox(); auto const& wface = zvel_pert.array(mfi); - auto const& cc = mf_cc_fine.const_array(mfi); + //auto const& cc = mf_cc_fine.const_array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {