update pressure via eos when transverse_use_eos=1 (#3004)

AMReX-Astro · Dec 13, 2024 · 4b1162f · 4b1162f
1 parent 50ed979
commit 4b1162f
Showing 1 changed file with 55 additions and 13 deletions.
diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp
@@ -437,18 +437,39 @@ Castro::actual_trans_single(const Box& bx,  // NOLINT(readability-convert-member
                 qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT);
             }
 
-            // Pretend QREINT has been fixed and transverse_use_eos != 1.
-            // If we are wrong, we will fix it later.
+            // Pretend QREINT has been fixed
+            // We can get pressure update via eos or using p-evolution equation
 
-            // Add the transverse term to the p evolution eq here.
+            if (transverse_use_eos) {
+
+                // With the EOS route:
+                eos_rep_t eos_state;
+                eos_state.rho = rrnewn;
+                eos_state.e = qo_arr(i,j,k,QREINT) / rrnewn;
+                for (int n = 0; n < NumSpec; n++) {
+                    eos_state.xn[n] = qo_arr(i,j,k,QFS+n);
+                }
+#if NAUX_NET > 0
+                for (int n = 0; n < NumAux; n++) {
+                    eos_state.aux[n] = qo_arr(i,j,k,QFX+n);
+                }
+#endif
+                eos(eos_input_re, eos_state);
+
+                Real pnewn = eos_state.p;
+                qo_arr(i,j,k,QPRES) = amrex::max(pnewn, small_p);
+
+            } else {
+
+                // Add the transverse term to the p evolution eq here.
 #if AMREX_SPACEDIM == 2
-            // the divergences here, dup and du, already have area factors
-            Real pnewn = q_arr(i,j,k,QPRES) - hdt * (dup + pav * du * (gamc - 1.0_rt)) * volinv;
+                // the divergences here, dup and du, already have area factors
+                Real pnewn = q_arr(i,j,k,QPRES) - hdt * (dup + pav * du * (gamc - 1.0_rt)) * volinv;
 #else
-            Real pnewn = q_arr(i,j,k,QPRES) - cdtdx * (dup + pav * du * (gamc - 1.0_rt));
+                Real pnewn = q_arr(i,j,k,QPRES) - cdtdx * (dup + pav * du * (gamc - 1.0_rt));
 #endif
-            qo_arr(i,j,k,QPRES) = amrex::max(pnewn, small_p);
-
+                qo_arr(i,j,k,QPRES) = amrex::max(pnewn, small_p);
+            }
         }
         else {
             qo_arr(i,j,k,QPRES) = q_arr(i,j,k,QPRES);
@@ -887,12 +908,33 @@ Castro::actual_trans_final(const Box& bx,  // NOLINT(readability-convert-member-
                 qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT);
             }
 
-            // Pretend QREINT has been fixed and transverse_use_eos != 1.
-            // If we are wrong, we will fix it later.
+            // Pretend QREINT has been fixed
+            // We can get pressure update via eos or using p-evolution equation
+
+            if (transverse_use_eos) {
+
+                // With the EOS route:
+                eos_rep_t eos_state;
+                eos_state.rho = rrnewn;
+                eos_state.e = qo_arr(i,j,k,QREINT) / rrnewn;
+                for (int n = 0; n < NumSpec; n++) {
+                    eos_state.xn[n] = qo_arr(i,j,k,QFS+n);
+                }
+#if NAUX_NET > 0
+                for (int n = 0; n < NumAux; n++) {
+                    eos_state.aux[n] = qo_arr(i,j,k,QFX+n);
+                }
+#endif
+                eos(eos_input_re, eos_state);
+                qo_arr(i,j,k,QPRES) = eos_state.p;
+
+            } else {
+
+                // add the transverse term to the p evolution eq here
+                Real pnewn = q_arr(i,j,k,QPRES) - pt1new - pt2new;
+                qo_arr(i,j,k,QPRES) = pnewn;
+            }
 
-            // add the transverse term to the p evolution eq here
-            Real pnewn = q_arr(i,j,k,QPRES) - pt1new - pt2new;
-            qo_arr(i,j,k,QPRES) = pnewn;
         }
         else {
             qo_arr(i,j,k,QPRES) = q_arr(i,j,k,QPRES);