Convective boundary

Exec/science/convective_boundary/GNUmakefile

@@ -0,0 +1,27 @@
+DEBUG      = FALSE
+DIM        = 2
+COMP	   = gnu
+USE_MPI    = FALSE
+USE_OMP    = FALSE
+USE_REACT  = TRUE
+
+TINY_PROFILE = FALSE
+PROFILE      = FALSE # TRUE overrides TINY_PROFILE
+
+# define the location of the MAESTROEX home directory
+MAESTROEX_HOME  := ../../..
+
+# Set the EOS, conductivity, and network directories
+# We first check if these exist in $(MAESTROEX_HOME)/Microphysics/(EOS/conductivity/networks)
+# If not we use the version in $(MICROPHYSICS_HOME)/Microphysics/(EOS/conductivity/networks)
+EOS_DIR          := gamma_law
+CONDUCTIVITY_DIR := constant
+NETWORK_DIR      := general_null
+NETWORK_INPUTS := simple.net
+INTEGRATOR_DIR   := VODE
+
+Bpack   := ./Make.package
+Blocs   := .
+PROBIN_PARAMETER_DIRS := .
+
+include $(MAESTROEX_HOME)/Exec/Make.Maestro

Exec/science/convective_boundary/MaestroBaseState.cpp

@@ -0,0 +1,104 @@
+// Initialize analytic hydrostatic background
+
+#include <Maestro.H>
+
+using namespace amrex;
+
+auto get_rho0(Real z);
+auto get_p0(Real rho_0, Real X0);
+
+void Maestro::InitBaseState(BaseState<Real>& rho0, BaseState<Real>& rhoh0,
+                            BaseState<Real>& p0, const int lev) {
+    // timer for profiling
+    BL_PROFILE_VAR("Maestro::InitBaseState()", InitBaseState);
+
+    auto rho0_arr = rho0.array();
+    auto rhoh0_arr = rhoh0.array();
+    auto p0_arr = p0.array();
+    auto p0_init_arr = p0_init.array();
+    auto tempbar_arr = tempbar.array();
+    auto tempbar_init_arr = tempbar_init.array();
+    auto s0_init_arr = s0_init.array();
+
+    const auto iH = network_spec_index("H2");
+    const auto iH2O = network_spec_index("H2O");
+
+    const auto g0    = problem_rp::g0;
+    const auto rho_t = problem_rp::rho_t;
+    const auto T_0   = problem_rp::T_0;
+    const auto k_B   = problem_rp::k_B;
+    const auto m_p   = problem_rp::m_p;
+    const auto X_b   = problem_rp::X_b;
+    const auto N_rho = problem_rp::N_rho;
+
+    const auto alpha = -N_rho / (std::log(9.0_rt / (9.0_rt - 8.0_rt * X_b)));
+    const auto rho_b = rho_t * std::exp(N_rho);
+    const auto R = k_B / m_p;
+    const auto D = (1.0_rt + alpha) * 4.0_rt * X_b * R * T_0 / (9.0_rt * g0);
+
+
+    // define some helper functions with lambdas
+    auto get_rho0 = [=](Real z) {
+        // return the background density
+        return rho_b * std::pow(1.0_rt + 8.0_rt * X_b * z / (9.0_rt - 8.0_rt * X_b) / D, alpha);
+    };
+
+    auto get_p0 = [=](Real rho_0, Real X0) {
+        // return the background pressure
+        return R * rho_0 * T_0 * (9.0_rt - 8.0_rt * X0) / 18.0_rt;
+    };
+
+    const int n = lev;
+    eos_rh_t eos_state;
+
+    for (auto r = 0; r < base_geom.nr(n); ++r) {
+        Real z = geom[lev].ProbLo(AMREX_SPACEDIM - 1) +
+                 (Real(r) + 0.5) * base_geom.dr[n];
+
+        eos_state.rho = get_rho0(z);
+        // Print() << "z/D = " << z/D << ",  rho/rho_t = " << eos_state.rho/rho_t << std::endl;
+
+        auto X0 = X_b * (1.0_rt - z / D);
+        eos_state.p = get_p0(eos_state.rho, X0);
+        for (auto comp = 0; comp < NumSpec; ++comp) {
+            eos_state.xn[comp] = 0.0_rt;
+        }
+        eos_state.xn[iH2O] = X0;
+        eos_state.xn[iH] = 1.0_rt - X0;
+        eos(eos_input_rp, eos_state);
+
+        // These prints validate that we are putting in the intended model for rho, p, composition
+        // Print() << "z = " << z << std::endl;
+        // Print() << "rho = " << eos_state.rho << std::endl;
+        // Print() << "p = " << eos_state.p << std::endl;
+        // Print() << "X(H20) = " << eos_state.xn[iH2O] << std::endl;
+        // Print() << "X(H2) = " << eos_state.xn[iH] << std::endl;
+        // Print() << "" << std::endl;
+
+        // Temperature (eos self-consistency check)
+        // Print() << "T_0 = " << T_0 << " " << "T_eos = " << eos_state.T << std::endl;
+
+        s0_init_arr(n, r, Rho) = eos_state.rho;
+        s0_init_arr(n, r, RhoH) = eos_state.rho * eos_state.h;
+        for (auto comp = 0; comp < NumSpec; ++comp) {
+            s0_init_arr(n, r, FirstSpec + comp) =
+                eos_state.rho * eos_state.xn[comp];
+        }
+        p0_init_arr(n, r) = eos_state.p;
+        s0_init_arr(n, r, Temp) = eos_state.T;
+    }
+
+    // copy s0_init and p0_init into rho0, rhoh0, p0, and tempbar
+    for (auto i = 0; i < base_geom.nr_fine; ++i) {
+        rho0_arr(lev, i) = s0_init_arr(lev, i, Rho);
+        rhoh0_arr(lev, i) = s0_init_arr(lev, i, RhoH);
+        tempbar_arr(lev, i) = s0_init_arr(lev, i, Temp);
+        tempbar_init_arr(lev, i) = s0_init_arr(lev, i, Temp);
+        p0_arr(lev, i) = p0_init_arr(lev, i);
+
+        // Print() << "p0 = " << p0_arr(lev,i) << std::endl;
+    }
+
+    // initialize any inlet BC parameters
+    SetInletBCs();
+}

Exec/science/convective_boundary/MaestroHeating.cpp

@@ -0,0 +1,60 @@
+// Apply cooling function at the top boundary
+
+#include <Maestro.H>
+
+using namespace amrex;
+
+// compute heating term, rho_Hext
+void Maestro::MakeHeating(Vector<MultiFab>& rho_Hext,
+                          const Vector<MultiFab>& scal) {
+    // timer for profiling
+    BL_PROFILE_VAR("Maestro::MakeHeating()", MakeHeating);
+
+    for (int lev = 0; lev <= finest_level; ++lev) {
+        // loop over boxes (make sure mfi takes a cell-centered multifab as an argument)
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+        for (MFIter mfi(rho_Hext[lev], TilingIfNotGPU()); mfi.isValid();
+             ++mfi) {
+            // Get the index space of the valid region
+            const Box& tileBox = mfi.tilebox();
+            const auto dx = geom[lev].CellSizeArray();
+            const auto prob_lo = geom[lev].ProbLoArray();
+            const auto prob_hi = geom[lev].ProbHiArray();
+
+            const Array4<const Real> scal_arr = scal[lev].array(mfi);
+            const Array4<Real> rho_Hext_arr = rho_Hext[lev].array(mfi);
+
+            ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                // Real r = (AMREX_SPACEDIM == 2) ? j : k;
+                // Real z = prob_lo[AMREX_SPACEDIM - 1] +
+                //          (Real(r) + 0.5) * dx[AMREX_SPACEDIM - 1];
+
+                rho_Hext_arr(i, j, k) = 0.0;
+
+                // if (z > 0.99 * prob_hi[AMREX_SPACEDIM - 1]) {
+                //     rho_Hext_arr(i, j, k) = scal_arr(i, j, k, Rho) * heat_flux_loc;
+                // }
+
+                if (problem_rp::do_cooling) {
+
+                    const Real L_x = prob_hi[0];
+
+                    Real x = (Real(i) + 0.5) * dx[0] + prob_lo[0];
+                    Real y = (Real(j) + 0.5) * dx[1] + prob_lo[1];
+                    Real top_y = prob_hi[1];
+                    Real delta_y = top_y - y;
+
+                    // Exponential decay from the top boundary
+                    Real er = std::exp(- delta_y * delta_y / 1.e11);
+                    rho_Hext_arr(i, j, k) = er * problem_rp::constant_heat_flux;
+                }
+
+            });
+        }
+    }
+
+    // average down
+    AverageDown(rho_Hext, 0, 1);
+}

Exec/science/convective_boundary/MaestroInitData.cpp

@@ -0,0 +1,108 @@
+// Initialize and perturb the temperature
+
+#include <Maestro.H>
+#include <random>
+
+using namespace amrex;
+
+// initializes data on a specific level
+void Maestro::InitLevelData(const int lev, const Real time, const MFIter& mfi,
+                            const Array4<Real> scal, const Array4<Real> vel) {
+    // timer for profiling
+    BL_PROFILE_VAR("Maestro::InitLevelData()", InitLevelData);
+
+    const auto tileBox = mfi.tilebox();
+
+    // set velocity to zero
+    ParallelFor(tileBox, AMREX_SPACEDIM,
+                [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) {
+                    vel(i, j, k, n) = 0.0;
+                });
+
+    const auto s0_arr = s0_init.const_array();
+    const auto p0_arr = p0_init.const_array();
+
+    ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+        int r = AMREX_SPACEDIM == 2 ? j : k;
+
+        // set the scalars using s0
+        scal(i, j, k, Rho) = s0_arr(lev, r, Rho);
+        scal(i, j, k, RhoH) = s0_arr(lev, r, RhoH);
+        scal(i, j, k, Temp) = s0_arr(lev, r, Temp);
+        for (auto comp = 0; comp < NumSpec; ++comp) {
+            scal(i, j, k, FirstSpec + comp) = s0_arr(lev, r, FirstSpec + comp);
+        }
+#if NAUX_NET > 0
+        for (auto comp = 0; comp < NumAux; ++comp) {
+            scal(i, j, k, FirstAux + comp) = s0_arr(lev, r, FirstAux + comp);
+        }
+#endif
+        // initialize pi to zero for now
+        scal(i, j, k, Pi) = 0.0;
+
+    });
+
+    if (perturb_model) {
+
+        // Generate random seed (or get from input), and RNG
+        int seed = 0;
+        if (problem_rp::pert_seed != -1) {
+            seed = problem_rp::pert_seed;
+        } else {
+            std::random_device r;
+            seed = r();
+        }
+        std::default_random_engine generator(seed);
+
+        // Gaussian distribution
+        std::normal_distribution<double> noise(1.0, 0.001); // mean of 1, std dev 10^-3
+
+        const auto lo = amrex::lbound(tileBox);
+        const auto hi = amrex::ubound(tileBox);
+
+        // Loop through data and perturb
+        // Because for and not parallelfor, need to run the first step not in serial (i.e. not with srun)
+        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++) {
+                    int r = AMREX_SPACEDIM == 2 ? j : k;  // j if 2D, k if 3D
+
+                    Real t0 = s0_arr(lev, r, Temp);
+                    Real temp = t0 * noise(generator);
+                    Real dens = s0_arr(lev, r, Rho); // no change
+
+                    // Create new eos object based on these modified values
+                    eos_t eos_state;
+                    eos_state.T = temp;
+                    eos_state.p = p0_arr(lev, r);
+                    eos_state.rho = dens;
+                    for (auto comp = 0; comp < NumSpec; comp++) {
+                        eos_state.xn[comp] =
+                            s0_arr(lev, r, FirstSpec + comp) / s0_arr(lev, r, Rho);
+                    }
+
+                    auto eos_input_flag = eos_input_tp; // temperature & pressure eos
+                    eos(eos_input_flag, eos_state);
+                    scal(i, j, k, Rho) = eos_state.rho;
+                    scal(i, j, k, RhoH) = eos_state.rho * eos_state.h; // re-compute enthalpy
+                    scal(i, j, k, Temp) = eos_state.T;
+                    for (auto comp=0; comp < NumSpec; comp++) {
+                        scal(i, j, k, FirstSpec + comp) =
+                            eos_state.rho * eos_state.xn[comp];
+                    }
+                }
+            }
+        }
+    }
+}
+
+void Maestro::InitLevelDataSphr(const int lev, const Real time, MultiFab& scal,
+                                MultiFab& vel) {
+
+    amrex::ignore_unused(lev);
+    amrex::ignore_unused(time);
+    amrex::ignore_unused(scal);
+    amrex::ignore_unused(vel);
+
+    Abort("InitLevelDataSphr not implemented");
+}

Exec/science/convective_boundary/_parameters

@@ -0,0 +1,17 @@
+@namespace: problem
+
+# flux boundary condition
+do_cooling               integer   1
+constant_heat_flux       real      -1.0e6
+
+# seed for the random noise (-1 for random seed)
+pert_seed                integer   2024
+
+# initial analytic model parameters
+g0                       real      -2479.0
+rho_t                    real      16.e-5
+T_0                      real      165.0
+k_B                      real      1.38e-16
+m_p                      real      1.67e-24
+X_b                      real      0.2
+N_rho                    real      3.0