Merge branch 'development' into fix_exact_riemann_pstar_test

AMReX-Astro · Jul 9, 2024 · 51f2c9f · 51f2c9f
2 parents bbcd4d8 + b373e92
commit 51f2c9f
Show file tree

Hide file tree

Showing 9 changed files with 34 additions and 26 deletions.
diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H
@@ -16,10 +16,6 @@ using namespace castro;
 
 #include <fundamental_constants.H>
 
-constexpr int HISTORY_SIZE=40;
-constexpr int PSTAR_BISECT_FACTOR = 5;
-
-
 #include <AMReX_BC_TYPES.H>
 #include <AMReX_AmrLevel.H>
 #include <AMReX_iMultiFab.H>

diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
@@ -58,6 +58,8 @@
 #include <ambient.H>
 #include <castro_limits.H>
 
+#include <riemann_constants.H>
+
 using namespace amrex;
 
 bool         Castro::signalStopJob = false;
@@ -392,7 +394,7 @@ Castro::read_params ()
 #endif
 
     if (riemann_solver == 1) {
-        if (riemann_shock_maxiter > HISTORY_SIZE) {
+        if (riemann_shock_maxiter > riemann_constants::HISTORY_SIZE) {
             amrex::Error("riemann_shock_maxiter > HISTORY_SIZE");
         }
 

diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package
@@ -30,6 +30,7 @@ CEXE_sources += riemann.cpp
 CEXE_headers += HLL_solvers.H
 CEXE_headers += riemann_solvers.H
 CEXE_headers += riemann_2shock_solvers.H
+CEXE_headers += riemann_constants.H
 CEXE_headers += riemann_type.H
 CEXE_headers += slope.H
 CEXE_headers += reconstruction.H

diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H
@@ -5,6 +5,7 @@
 
 using namespace amrex::literals;
 
+#include <riemann_constants.H>
 #include <riemann_type.H>
 #ifdef RADIATION
 #include <Radiation.H>
@@ -60,7 +61,7 @@ namespace TwoShock {
                     const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax,
                     const int lriemann_shock_maxiter, const amrex::Real lriemann_pstar_tol,
                     amrex::Real& pstar, amrex::Real& gamstar,
-                    bool& converged, amrex::GpuArray<amrex::Real, PSTAR_BISECT_FACTOR*HISTORY_SIZE>& pstar_hist_extra) {
+                    bool& converged, amrex::GpuArray<amrex::Real, riemann_constants::PSTAR_BISECT_FACTOR * riemann_constants::HISTORY_SIZE>& pstar_hist_extra) {
 
       // we want to zero
       // f(p*) = u*_l(p*) - u*_r(p*)
@@ -106,7 +107,7 @@ namespace TwoShock {
       converged = false;
       amrex::Real pstar_c = 0.0;
 
-      for (int iter = 0; iter < PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) {
+      for (int iter = 0; iter < riemann_constants::PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) {
 
         pstar_c = 0.5_rt * (pstar_lo + pstar_hi);
         pstar_hist_extra[iter] = pstar_c;
@@ -167,7 +168,7 @@ namespace TwoShock {
       constexpr amrex::Real weakwv = 1.e-3_rt;
 
     #ifndef AMREX_USE_GPU
-      amrex::GpuArray<amrex::Real, HISTORY_SIZE> pstar_hist;
+      amrex::GpuArray<amrex::Real, riemann_constants::HISTORY_SIZE> pstar_hist;
     #endif
 
 
@@ -336,7 +337,7 @@ namespace TwoShock {
               pstarl = amrex::max(pstarl, small_pres);
               pstaru = amrex::max(pstaru, small_pres);
 
-              amrex::GpuArray<amrex::Real, PSTAR_BISECT_FACTOR*HISTORY_SIZE> pstar_hist_extra;
+              amrex::GpuArray<amrex::Real, riemann_constants::PSTAR_BISECT_FACTOR * riemann_constants::HISTORY_SIZE> pstar_hist_extra;
 
               pstar_bisection(pstarl, pstaru,
                               ql.un, ql.p, taul, gamel, clsql,
@@ -352,7 +353,7 @@ namespace TwoShock {
                       std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl;
                   }
                   std::cout << "pstar extra history: " << std::endl;
-                  for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) {
+                  for (int iter_l = 0; iter_l < riemann_constants::PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) {
                       std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl;
                   }
 

diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H
@@ -0,0 +1,16 @@
+#ifndef RIEMANN_CONSTANTS_H
+#define RIEMANN_CONSTANTS_H
+
+#include <AMReX_REAL.H>
+
+using namespace amrex::literals;
+
+namespace riemann_constants {
+    constexpr amrex::Real smlp1 = 1.e-10_rt;
+    constexpr amrex::Real small = 1.e-8_rt;
+    constexpr amrex::Real smallu = 1.e-12_rt;
+    constexpr int HISTORY_SIZE=40;
+    constexpr int PSTAR_BISECT_FACTOR = 5;
+}
+
+#endif
diff --git a/Source/hydro/riemann_type.H b/Source/hydro/riemann_type.H
@@ -5,13 +5,6 @@
 
 using namespace amrex::literals;
 
-namespace riemann_constants {
-    const amrex::Real smlp1 = 1.e-10_rt;
-    const amrex::Real small = 1.e-8_rt;
-    const amrex::Real smallu = 1.e-12_rt;
-}
-
-
 struct RiemannState
 {
     amrex::Real rho;

diff --git a/Util/exact_riemann/GNUmakefile b/Util/exact_riemann/GNUmakefile
@@ -33,6 +33,6 @@ Blocs   := .
 
 # we explicitly want runtime_parameters.H so we have access to
 # Castro's runtime parameter system
-Blocs += $(CASTRO_HOME)/Source/driver
+Blocs += $(CASTRO_HOME)/Source/driver $(CASTRO_HOME)/Source/hydro
 
 include $(CASTRO_HOME)/Exec/Make.Castro
diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp
@@ -41,9 +41,7 @@ int main(int argc, char *argv[]) {
     network_init();
 #endif
 
-    amrex::Real small_temp = 1.e-200;
-    amrex::Real small_dens = 1.e-200;
-    eos_init(small_temp, small_dens);
+    eos_init(castro::small_temp, castro::small_dens);
 
     exact_riemann();
 

diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H
@@ -6,6 +6,7 @@
 #include <eos_type.H>
 #include <eos.H>
 
+#include <riemann_constants.H>
 
 AMREX_INLINE
 void
@@ -62,7 +63,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
     converged = false;
 
     int iter = 1;
-    while (! converged && iter < riemann_solver::max_iters) {
+    while (! converged && iter < castro::riemann_shock_maxiter) {
 
         amrex::Real rhostar_s, f, fprime;
 
@@ -71,7 +72,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
 
         amrex::Real dW = -f / fprime;
 
-        if (std::abs(dW) < riemann_solver::tol * W_s) {
+        if (std::abs(dW) < castro::riemann_pstar_tol * W_s) {
             converged = true;
         }
 
@@ -96,7 +97,7 @@ shock(const amrex::Real pstar,
 
     amrex::ignore_unused(u_s);
 
-    amrex::Real rhostar_hist[riemann_solver::max_iters], Ws_hist[riemann_solver::max_iters];
+    amrex::Real rhostar_hist[riemann_constants::HISTORY_SIZE], Ws_hist[riemann_constants::HISTORY_SIZE];
 
     // compute the Z_s function for a shock following C&G Eq. 20 and
     // 23.  Here the "_s" variables are the state (L or R) that we are
@@ -147,7 +148,7 @@ shock(const amrex::Real pstar,
     amrex::Real rhostar_s;
 
     if (taustar_s < 0.0_rt) {
-        rhostar_s = riemann_solver::smallrho;
+        rhostar_s = riemann_constants::small;
         W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s));
     }
 
@@ -162,7 +163,7 @@ shock(const amrex::Real pstar,
     // now did we converge?
 
     if (! converged) {
-        for (int i = 0; i < riemann_solver::max_iters; ++i) {
+        for (int i = 0; i < castro::riemann_shock_maxiter; ++i) {
             std::cout << i << " " << rhostar_hist[i] << " " << Ws_hist[i] << std::endl;
         }