xrb layered updates (#396)

This syncs up to the version that was used in the xrb_layered paper

Exec/science/xrb_layered/.gitignore

@@ -0,0 +1,21 @@
+*.ipynb_checkpoints/
+**/.vscode
+*slurm*
+*submit*
+*out*
+*debug*
+log.md
+analysis
+
+# some useful analysis scripts
+scripts/*
+!scripts/average_vars.sh
+!scripts/avg.submit
+!scripts/chainslurm.sh
+!scripts/conv_grad.submit
+!scripts/derived.submit
+!scripts/get_extrema.sh
+!scripts/extrema.submit
+!scripts/rms.submit
+!scripts/rms_vars.sh
+!scripts/run.submit.template

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,121 +39,60 @@ 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 ? problem_rp::xrb_pert_height : center[1];
-#if AMREX_SPACEDIM == 3
-        const auto zcen = problem_rp::xrb_pert_height;
-#endif
-
-        const auto rad_pert =
-            -problem_rp::xrb_pert_size * problem_rp::xrb_pert_size / (4.0 * std::log(0.5));
-
-        const bool perturb_temp_true = problem_rp::xrb_pert_type == 1;
-
-        const auto xrb_pert_factor_loc = problem_rp::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 (problem_rp::apply_vel_field) {
-        const auto velpert_amplitude_loc = problem_rp::velpert_amplitude;
-        const auto velpert_scale_loc = problem_rp::velpert_scale;
-        const auto num_vortices_loc = problem_rp::num_vortices;
-        const auto velpert_height = problem_rp::velpert_height_loc;
-
-        const Real offset = (prob_hi[0] - prob_lo[0]) / (problem_rp::num_vortices);
-
-        // vortex x-coords
-        RealVector vortices_xloc(problem_rp::num_vortices);
-        for (auto i = 0; i < problem_rp::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, 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];
+                    }
+                }
             }
-            vel(i, j, k, 0) += upert;
-            vel(i, j, k, 1) += vpert;
-        });
+        }
     }
 }
 

Exec/science/xrb_layered/MaestroTagging.cpp

@@ -0,0 +1,71 @@
+
+#include <Maestro.H>
+
+using namespace amrex;
+
+void Maestro::RetagArray(const Box& bx, const int lev) {
+    // timer for profiling
+    BL_PROFILE_VAR("Maestro::RetagArray()", RetagArray);
+
+    // re-compute tag_array since the actual grid structure changed due to buffering
+    // this is required in order to compute numdisjointchunks, r_start_coord, r_end_coord
+
+    // Tag on regions including buffer cells
+    auto lo = bx.loVect3d()[AMREX_SPACEDIM - 1];
+    auto hi = bx.hiVect3d()[AMREX_SPACEDIM - 1];
+    const auto max_lev = base_geom.max_radial_level + 1;
+
+    for (auto r = lo; r <= hi; ++r) {
+        tag_array[lev - 1 + max_lev * (r / 2)] = TagBox::SET;
+    }
+}
+
+void Maestro::TagBoxes(TagBoxArray& tags, const MFIter& mfi, const int lev,
+                       const Real time) {
+    // timer for profiling
+    BL_PROFILE_VAR("Maestro::TagBoxes()", TagBoxes);
+
+    // tag all cells at a given height if any cells at that height were tagged
+
+    const Array4<TagBox::TagType> tag = tags.array(mfi);
+    const int* AMREX_RESTRICT tag_array_p = tag_array.dataPtr();
+    const int max_lev = base_geom.max_radial_level + 1;
+
+    const Box& tilebox = mfi.tilebox();
+
+    ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+        int r = AMREX_SPACEDIM == 2 ? j : k;
+
+        if (tag_array_p[lev + max_lev * r] > 0) {
+            tag(i, j, k) = TagBox::SET;
+        }
+    });
+}
+
+void Maestro::StateError(TagBoxArray& tags, const MultiFab& state_mf,
+                         const MFIter& mfi, const int lev, const Real time) {
+    // timer for profiling
+    BL_PROFILE_VAR("Maestro::StateError()", StateError);
+
+    const Box& tileBox = mfi.tilebox();
+
+    const auto tag = tags.array(mfi);
+    const Array4<const Real> state = state_mf.array(mfi);
+    int* AMREX_RESTRICT tag_array_p = tag_array.dataPtr();
+    const int max_lev = base_geom.max_radial_level + 1;
+
+    const Real dr_lev = base_geom.dr(lev);
+
+    // Tag based on height
+    if (problem_rp::do_height_based_refinement) {
+        ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+            int r = AMREX_SPACEDIM == 2 ? j : k;
+            Real height = (Real(r) + 0.5) * dr_lev;
+
+            if (height > problem_rp::height_refine_bottom && height < problem_rp::height_refine_top) {
+                tag(i, j, k) = TagBox::SET;
+                tag_array_p[lev + max_lev * r] = TagBox::SET;
+            }
+        });
+    } 
+}

Exec/science/xrb_layered/_parameters

@@ -1,63 +1,15 @@
 @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               integer            0
 
 # minimum value to use for the sponge
 sponge_min                          real               0.01d0
 
-# Mass fraction used in {\tt diag.f90} to define where the burning layer @@
-# begins.
-diag_define_layer                real               0.1
-
-# Velocity field initialization
-apply_vel_field                     integer            0
-
-# 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                      integer            0
-
-# minimum value to use for species in tag_boxes
-tag_minval                          real               1.0d-4
-
-# maximum value to use for species in tag_boxes
-tag_maxval                          real               0.99d0
-
-# factor to multiply the species by; for tagging off of spec
-# instead of EGR
-tag_xfac                           real                1.d-16
-
-# density tagging
-do_dens_tagging                    integer             0
-
-lo_dens_tag                        real                5.4d5
-hi_dens_tag                        real                1.75d6
-dens_tag_lev_fac                   real                1.0d0
+# Do refinement above and below given heights (given in cm)
+do_height_based_refinement          integer            0
+height_refine_bottom                real               1.0e3
+height_refine_top                   real               1.2e3