[WIP] xrb layered

Exec/science/xrb_layered/.gitignore

@@ -0,0 +1,6 @@
+*.ipynb_checkpoints/
+**/.vscode
+*slurm*
+log.md
+analysis
+scripts

Exec/science/xrb_layered/MaestroInitData.cpp

@@ -1,5 +1,6 @@
-
 #include <Maestro.H>
+#include <random>
+
 using namespace amrex;
 
 // initializes data on a specific level
@@ -38,131 +39,59 @@ void Maestro::InitLevelData(const int lev, [[maybe_unused]] const Real time, con
     const auto dx = geom[lev].CellSizeArray();
 
     if (perturb_model) {
-        const auto xcen = center[0];
-        const auto ycen = AMREX_SPACEDIM == 2 ? xrb_pert_height : center[1];
-#if AMREX_SPACEDIM == 3
-        const auto zcen = xrb_pert_height;
-#endif
-
-        const auto rad_pert =
-            -xrb_pert_size * xrb_pert_size / (4.0 * std::log(0.5));
-
-        const bool perturb_temp_true = xrb_pert_type == 1;
-
-        const auto xrb_pert_factor_loc = xrb_pert_factor;
-        const auto rad_pert_loc = rad_pert;
-
-        ParallelFor(tileBox, [=] (int i, int j, int k) {
-            int r = AMREX_SPACEDIM == 2 ? j : k;
-
-            const auto x = prob_lo[0] + (Real(i) + 0.5) * dx[0] - xcen;
-            const auto y = prob_lo[1] + (Real(j) + 0.5) * dx[1] - ycen;
-#if AMREX_SPACEDIM == 3
-            const auto z = prob_lo[2] + (Real(k) + 0.5) * dx[2] - zcen;
-#else
-            const Real z = 0.0;
-#endif
-
-            const auto dist = std::sqrt(x * x + y * y + z * z);
-
-            Real temp = 0.0;
-            Real dens = 0.0;
-            auto eos_input_flag = eos_input_tp;
-
-            if (perturb_temp_true) {
-                Real t0 = s0_arr(lev, r, Temp);
-                temp = t0 * (1.0 + xrb_pert_factor_loc *
-                                       std::exp(-dist * dist / rad_pert_loc));
-                dens = s0_arr(lev, r, Rho);
-                eos_input_flag = eos_input_tp;
-            } else {
-                Real d0 = s0_arr(lev, r, Rho);
-                dens = d0 * (1.0 + xrb_pert_factor_loc *
-                                       std::exp(-dist * dist / rad_pert_loc));
-                temp = s0_arr(lev, r, Temp);
-                eos_input_flag = eos_input_rp;
-            }
-
-            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);
-            }
-
-            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;
-            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];
-            }
-        });
-    }
-
-    if (apply_vel_field) {
-        const auto velpert_amplitude_loc = velpert_amplitude;
-        const auto velpert_scale_loc = velpert_scale;
-        const auto num_vortices_loc = num_vortices;
-        const auto velpert_height = velpert_height_loc;
-
-        const Real offset = (prob_hi[0] - prob_lo[0]) / (num_vortices);
-
-        // vortex x-coords
-        RealVector vortices_xloc(num_vortices);
-        for (auto i = 0; i < num_vortices; ++i) {
-            vortices_xloc[i] = (static_cast<Real>(i) + 0.5_rt) * offset;
+
+        // 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();
         }
-
-        const Real* vortices_xloc_p = vortices_xloc.dataPtr();
-
-        ParallelFor(tileBox, [=] (int i, int j, int k) {
-            const auto x = prob_lo[0] + (Real(i) + 0.5) * dx[0];
-            const auto y = prob_lo[1] + (Real(j) + 0.5) * dx[1];
-
-            Real ydist = y - velpert_height;
-
-            Real upert = 0.0;
-            Real vpert = 0.0;
-
-            for (auto vortex = 0; vortex < num_vortices_loc; ++vortex) {
-                Real xdist = x - vortices_xloc_p[vortex];
-
-                Real rad = std::sqrt(xdist * xdist + ydist * ydist);
-
-                // e.g. Calder et al. ApJSS 143, 201-229 (2002)
-                // we set things up so that every other vortex has the same
-                // orientation
-                upert -=
-                    ydist * velpert_amplitude_loc *
-                    std::exp(-rad * rad /
-                             (2.0 * velpert_scale_loc * velpert_scale_loc)) *
-                    pow(-1.0, vortex + 1);
-
-                vpert +=
-                    xdist * velpert_amplitude_loc *
-                    std::exp(-rad * rad /
-                             (2.0 * velpert_scale_loc * velpert_scale_loc)) *
-                    pow(-1.0, vortex + 1);
+        std::default_random_engine generator(seed);
+
+        // Gaussian distribution
+        std::normal_distribution<double> noise(1.0, 0.001); // mean of 1, std dev 10^-3
+
+        // Loop through data and perturb
+        const auto lo = amrex::lbound(tileBox);
+        const auto hi = amrex::ubound(tileBox);
+
+        // Why for and not parallelfor?
+        // Because RNG's are not thread safe. so if we wanted to use ParallelFor, we would have to create a new RNG inside
+        // every iteration, but that loses the ability to keep a fixed seed and be fully reproducible.. maybe that wouldn't be so bad
+        // but (probably?) because of this standard for, this initialization can't be run in parallel (with e.g. srun)
+        // Have to run as a single process until at least get an initial checkfile, from which we can run in parallel.
+        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, j, 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];
+                    }
+                }
             }
-            vel(i, j, k, 0) += upert;
-            vel(i, j, k, 1) += vpert;
-        });
+        }
     }
 }
-
-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/xrb_layered/MaestroSponge.cpp

@@ -60,7 +60,7 @@ void Maestro::SpongeInit(const BaseState<Real>& rho0_s) {
     }
 
     if (maestro_verbose >= 1) {
-        if (xrb_use_bottom_sponge) {
+        if (problem_rp::xrb_use_bottom_sponge) {
             Print() << "inner sponge: r_sp      , r_tp      : " << r_sp << ", "
                     << r_tp << std::endl;
         }
@@ -78,7 +78,7 @@ void Maestro::MakeSponge(Vector<MultiFab>& sponge) {
     const Real botsponge_hi_r = r_tp;
     const Real topsponge_lo_r = r_sp_outer;
     const Real topsponge_hi_r = r_tp_outer;
-    const Real sponge_min_loc = sponge_min;
+    const Real sponge_min_loc = problem_rp::sponge_min;
 
     if (AMREX_SPACEDIM != 2) {
         Abort("ERROR: sponge only supported for 2d in xrb_mixed");
@@ -105,7 +105,7 @@ void Maestro::MakeSponge(Vector<MultiFab>& sponge) {
             ParallelFor(tileBox, [=] (int i, int j, int k)
                         { sponge_arr(i, j, k) = 1.0; });
 
-            if (xrb_use_bottom_sponge) {
+            if (problem_rp::xrb_use_bottom_sponge) {
 
                 // do both the top and bottom sponges
 

Exec/science/xrb_layered/_parameters

@@ -1,18 +1,7 @@
 @namespace: problem
 
-# Magnitude of the perturbation.
-xrb_pert_factor                     real               1.d-2               
-
-# Radius of the perturbation.
-xrb_pert_size                       real               5.0d1               
-
-# Type of perturbation @@
-# 1 = Temperature perturbation @@
-# 2 = Density perturbation
-xrb_pert_type                       integer            1                   
-
-# the height location for the perturbation
-xrb_pert_height                     real               -1.0d0
+# seed for the random noise (-1 for random seed)
+pert_seed                           integer            -1
 
 # Turn on (.true.) or off (.false.) sponging at the bottom of the domain.
 xrb_use_bottom_sponge               logical            .false.
@@ -24,24 +13,6 @@ sponge_min                          real               0.01d0
 # begins.
 diag_define_layer                real               0.1
 
-# Velocity field initialization
-apply_vel_field                     logical            .false.
-
-# Initial velocity field vortex scale (2-d); thickness of velocity layer (3-d)
-velpert_scale                       real               1.0d2
-
-# Initial velocity field vortex amplitude
-velpert_amplitude                   real               1.0d2
-
-# Initial velocity field vortex height location (2-d); base of perturbation layer (3-d)
-velpert_height_loc                  real               6.5d3
-
-# thickness of velocity pert transition region (3-d)
-velpert_steep                       real               12.d0
-
-# number of vortices (2-d)
-num_vortices                        integer            1
-
 # look at the deltap diagnostic
 do_deltap_diag                      logical            .false.