Merge branch 'development' into shock_paper_detonation

AMReX-Astro · Jun 24, 2024 · 1e03d10 · 1e03d10
2 parents 0a09d64 + 82b3d4e
commit 1e03d10
Show file tree

Hide file tree

Showing 12 changed files with 65 additions and 117 deletions.
diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H
@@ -522,7 +522,7 @@ public:
 /// @param state        MultiFab to check
 /// @param check_ghost  do we check the ghost cells?
 ///
-    static void check_for_nan(amrex::MultiFab& state, int check_ghost=0);
+    static void check_for_nan(const amrex::MultiFab& state, int check_ghost=0);
 
 
 #include <Castro_sources.H>
@@ -1138,7 +1138,7 @@ public:
 /// @param S        Current state
 /// @param time     current time
 ///
-    void define_new_center (amrex::MultiFab& S, amrex::Real time);
+    void define_new_center (const amrex::MultiFab& S, amrex::Real time);
     void write_center ();
 
 ///

diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
@@ -2782,8 +2782,8 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep)
             // Update the flux register now that we may have modified some of the flux corrections.
 
             for (OrientationIter fi; fi.isValid(); ++fi) {
-                FabSet& fs = (*reg)[fi()];
                 if (fi().coordDir() == idir) {
+                    FabSet& fs = (*reg)[fi()];
                     fs.copyFrom(temp_fluxes[idir], 0, 0, 0, temp_fluxes[idir].nComp());
                 }
             }
@@ -3551,7 +3551,7 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time
             {
                 bool outer_boundary_test[3] = {false};
 
-                int idx[3] = {i, j, k};
+                const int idx[3] = {i, j, k};
 
                 for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) {
 
@@ -4239,7 +4239,7 @@ Castro::get_numpts ()
 }
 
 void
-Castro::define_new_center(MultiFab& S, Real time)
+Castro::define_new_center(const MultiFab& S, Real time)
 {
     BL_PROFILE("Castro::define_new_center()");
 
@@ -4428,7 +4428,7 @@ Castro::expand_state(MultiFab& S, Real time, int ng)
 
 
 void
-Castro::check_for_nan(MultiFab& state_in, int check_ghost)
+Castro::check_for_nan(const MultiFab& state_in, int check_ghost)
 {
   BL_PROFILE("Castro::check_for_nan()");
 

diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H
@@ -97,9 +97,9 @@ void position(int i, int j, int k,
     // Given 3D indices (i,j,k), return the cell-centered spatial position.
     // Optionally we can also be edge-centered in any of the directions.
 
-    int idx[3] = {i, j, k};
+    const int idx[3] = {i, j, k};
 
-    bool cc[3] = {ccx, ccy, ccz};
+    const bool cc[3] = {ccx, ccy, ccz};
 
     Real offset[AMREX_SPACEDIM];
 

diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp
@@ -259,16 +259,13 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)  // NOLINT(readability-co
       shk.resize(obx, 1);
       fab_size += shk.nBytes();
 
-      Array4<Real> const shk_arr = shk.array();
-
       // Multidimensional shock detection
-      // Used for the hybrid Riemann solver
 
-#ifdef SHOCK_VAR
-      bool compute_shock = true;
-#else
-      bool compute_shock = false;
-#endif
+      // this is a local shock variable used only for the Riemann
+      // solver -- this will never be used to update the value in the
+      // conserved state.
+
+      Array4<Real> const shk_arr = shk.array();
 
       // get the primitive variable hydro sources
 
@@ -285,7 +282,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)  // NOLINT(readability-co
           hydro::src_to_prim(i, j, k, dt, U_old_arr, q_arr, old_src_arr, src_corr_arr, src_q_arr);
       });
 
-      if (hybrid_riemann == 1 || compute_shock) {
+      if (hybrid_riemann == 1) {
         shock(obx, q_arr, old_src_arr, shk_arr);
       }
       else {

diff --git a/Util/exact_riemann/_prob_params b/Util/exact_riemann/_prob_params
@@ -21,7 +21,3 @@ npts                integer  128          y
 use_Tinit           integer  0            y
 
 initial_temp_guess  real     1.0d5        y
-
-riemann_max_iter    integer  10           y
-
-co_moving_frame     integer  0            y
diff --git a/Util/exact_riemann/ci-benchmarks/test1-riemann.out b/Util/exact_riemann/ci-benchmarks/test1-riemann.out
diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H
@@ -47,12 +47,6 @@ exact_riemann() {
 
     }
 
-    if (problem::co_moving_frame) {
-        amrex::Real W_avg = 0.5_rt * (problem::u_l + problem::u_r);
-        problem::u_l -= W_avg;
-        problem::u_r -= W_avg;
-    }
-
     amrex::Real ustar, pstar, W_l, W_r;
 
     riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn,
@@ -86,12 +80,6 @@ exact_riemann() {
                        x, problem::xjump, problem::t,
                        rho, u, p, xn_s);
 
-        if (problem::co_moving_frame) {
-            amrex::Real W_avg = 0.5_rt * (problem::u_l + problem::u_r);
-            u += W_avg;
-            x += problem::t * W_avg;
-        }
-
         // get the thermodynamics for this state for output
 
         eos_t eos_state;

diff --git a/Util/exact_riemann/inputs.moving b/Util/exact_riemann/inputs.moving
diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp
@@ -10,16 +10,14 @@
 #include <prob_parameters.H>
 #include <eos.H>
 #include <network.H>
-#include <castro_params.H>
+
 #include <exact_riemann.H>
 
 int main(int argc, char *argv[]) {
 
     amrex::Initialize(argc, argv);
 
     std::cout << "starting the exact Riemann solver..." << std::endl;
-    std::cout << argv[1] << std::endl;
-    std::cout << strlen(argv[1]) << std::endl;
 
     // initialize the external runtime parameters in C++
 

diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/riemann_rarefaction.H
@@ -68,7 +68,6 @@ riemann_invariant_rhs2(const amrex::Real u, const amrex::Real tau, const amrex::
 
     amrex::ignore_unused(u);
 
-    std::cout << "in riemann_invariant_rhs2: " << u << " " << tau << " " << p << std::endl;
     eos_rep_t eos_state;
     eos_state.rho = 1.0_rt / tau;
     eos_state.p = p;
@@ -326,8 +325,6 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re
 
     bool finished = false;
 
-    std::cout << "integrating from u: " << u << " " << ustop << " " << xi << " " << c << std::endl;
-
     // this will be used as an initial guess to accelerate the EOS inversions
 
     amrex::Real T = eos_state.T;

diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H
@@ -129,10 +129,6 @@ shock(const amrex::Real pstar,
         2.0_rt * (1.0_rt - gammaE_bar / gammaC_bar) * (gammaE_bar - 1.0_rt) *
         (pstar - p_s) / (pstar + p_s);
 
-    std::cout << "pstar, ps = " << pstar << " " << p_s << " " << gammaE_s << " " << gammaE_star << std::endl;
-    std::cout << (pstar/rho_s - (gammaE_star - 1.0_rt)/(gammaE_s - 1.0_rt) * p_s/rho_s);
-    std::cout << (pstar + 0.5_rt * (gammaE_star - 1.0_rt) * (pstar + p_s));
-
     // there is a pathological case that if p_s - pstar ~ 0, the root finding
     // just doesn't work.  In this case, we use the ideas from CG, Eq. 35, and
     // take W = sqrt(gamma p rho)

diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H
@@ -80,9 +80,6 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
     // find the exact pstar and ustar
 
     std::cout << "solving for star state: " << std::endl;
-    std::cout << rho_l << " " << u_l << " " << p_l << std::endl;
-    std::cout << rho_r << " " << u_r << " " << p_r << std::endl;
-
 
     // this procedure follows directly from Colella & Glaz 1985, section 1
     // the basic idea is that we want to find the pstar that satisfies:
@@ -111,14 +108,13 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
     bool converged = false;
 
     int iter = 1;
-    while (! converged && iter < problem::riemann_max_iter) {
+    while (! converged && iter < riemann_solver::max_iters) {
 
         // compute Z_l, W_l and Z_r, W_r -- the form of these depend
         // on whether the wave is a shock or a rarefaction
 
         amrex::Real Z_l, Z_r;
 
-        std::cout << "iteration: " << iter << " " << pstar << " " << p_l << " " << p_r << std::endl;
         // left wave
 
         if (pstar - p_l > SMALL * p_l) {
@@ -146,9 +142,6 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
 
         amrex::Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r);
 
-        std::cout << "done with iteration " << iter << std::endl;
-        std::cout << "ustar_l/r, pstar: " << ustar_l << " " << ustar_r << " " << pstar_new << std::endl;
-
         // estimate the error in the current star solution
         amrex::Real err1 = std::abs(ustar_r - ustar_l);
         amrex::Real err2 = pstar_new - pstar;
@@ -164,6 +157,10 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
         iter++;
     }
 
+    if (! converged) {
+        amrex::Error("star state did not converge");
+    }
+
     ustar = 0.5_rt * (ustar_l + ustar_r);
 
     std::cout << "found pstar, ustar: " <<  pstar << " " << ustar << std::endl;