Merge branch 'exact_riemann_castro_params' into move_history_params_t…

…o_riemann
AMReX-Astro · Jul 1, 2024 · a3f5802 · a3f5802
2 parents 03b9653 + 0d51f94
commit a3f5802
Show file tree

Hide file tree

Showing 18 changed files with 215 additions and 185 deletions.
diff --git a/Exec/Make.auto_source b/Exec/Make.auto_source
@@ -72,9 +72,15 @@ CPP_PARAMETERS := $(TOP)/Source/driver/_cpp_parameters
 
 $(CASTRO_AUTO_SOURCE_DIR)/runtime_params.cpp: $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H
 
+STRUCT_USE_CASTRO_CLASS ?= TRUE
+PARSE_ARGS :=
+ifneq ($(STRUCT_USE_CASTRO_CLASS), TRUE)
+   PARSE_ARGS += --without-castro-class
+endif
+
 $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H: $(CPP_PARAMETERS)
 	@if [ ! -d $(CASTRO_AUTO_SOURCE_DIR) ]; then mkdir -p $(CASTRO_AUTO_SOURCE_DIR); fi
-	PYTHONPATH=$(MICROPHYSICS_HOME)/util/build_scripts $(TOP)/Source/driver/parse_castro_params.py -o $(CASTRO_AUTO_SOURCE_DIR) $(CPP_PARAMETERS)
+	PYTHONPATH=$(MICROPHYSICS_HOME)/util/build_scripts $(TOP)/Source/driver/parse_castro_params.py $(PARSE_ARGS) -o $(CASTRO_AUTO_SOURCE_DIR) $(CPP_PARAMETERS)
 
 # for debugging
 test_cxx_params: $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H

diff --git a/Source/driver/parse_castro_params.py b/Source/driver/parse_castro_params.py
@@ -129,7 +129,7 @@ def read_param_file(infile):
 
     return params
 
-def write_headers_and_source(params, out_directory, struct_name):
+def write_headers_and_source(params, out_directory, struct_name, without_castro_class):
 
     # output
 
@@ -176,19 +176,23 @@ def write_headers_and_source(params, out_directory, struct_name):
 
         cq.write(CWARNING)
 
+        class_name = "Castro"
+        if without_castro_class:
+            class_name = None
+
         for ifdef in ifdefs:
             if ifdef is None:
                 for p in [q for q in params_nm if q.ifdef is None]:
                     cq.write(p.get_default_string())
                     cq.write(p.get_query_string())
-                    cq.write(p.get_query_struct_string(struct_name=struct_name, class_name="Castro"))
+                    cq.write(p.get_query_struct_string(struct_name=struct_name, class_name=class_name))
                     cq.write("\n")
             else:
                 cq.write(f"#ifdef {ifdef}\n")
                 for p in [q for q in params_nm if q.ifdef == ifdef]:
                     cq.write(p.get_default_string())
                     cq.write(p.get_query_string())
-                    cq.write(p.get_query_struct_string(struct_name=struct_name, class_name="Castro"))
+                    cq.write(p.get_query_struct_string(struct_name=struct_name, class_name=class_name))
                     cq.write("\n")
                 cq.write("#endif\n")
             cq.write("\n")
@@ -305,13 +309,14 @@ def main():
                         help="output directory for the generated files")
     parser.add_argument("-s", type=str, default="params",
                         help="name for the name struct that will hold the parameters")
+    parser.add_argument("--without-castro-class", action="store_true", help="don't include Castro:: in the struct namespace")
     parser.add_argument("input_file", type=str, nargs=1,
                         help="input file containing the list of parameters we will define")
 
     args = parser.parse_args()
 
     p = read_param_file(args.input_file[0])
-    write_headers_and_source(p, args.o, args.s)
+    write_headers_and_source(p, args.o, args.s, args.without_castro_class)
 
 if __name__ == "__main__":
     main()
diff --git a/Source/driver/runtime_parameters.H b/Source/driver/runtime_parameters.H
@@ -1,6 +1,7 @@
 #ifndef RUNTIME_PARAMETERS_H
 #define RUNTIME_PARAMETERS_H
 
+#include <AMReX_Print.H>
 #include <AMReX_ParmParse.H>
 #include <castro_limits.H>
 
@@ -25,34 +26,34 @@ initialize_cpp_runparams()
 {
 
     {
-        ParmParse pp("castro");
+        amrex::ParmParse pp("castro");
 #include <castro_queries.H>
     }
 
 #ifdef AMREX_PARTICLES
     {
-        ParmParse pp("particles");
+        amrex::ParmParse pp("particles");
 #include <particles_queries.H>
     }
 #endif
 
 #ifdef DIFFUSION
     {
-        ParmParse pp("diffusion");
+        amrex::ParmParse pp("diffusion");
 #include <diffusion_queries.H>
     }
 #endif
 
 #ifdef GRAVITY
     {
-        ParmParse pp("gravity");
+        amrex::ParmParse pp("gravity");
 #include <gravity_queries.H>
     }
 #endif
 
 #ifdef RADIATION
     {
-        ParmParse pp("radiation");
+        amrex::ParmParse pp("radiation");
 #include <radiation_queries.H>
     }
 #endif
@@ -102,14 +103,14 @@ validate_runparams()
     for (const auto& nm: check_namespaces)
     {
         // "castro"
-        if (ParmParse::hasUnusedInputs(nm)) {
+        if (amrex::ParmParse::hasUnusedInputs(nm)) {
             amrex::Print() << "Warning: the following " + nm + ".* parameters are ignored\n";
-            auto unused = ParmParse::getUnusedInputs(nm);
+            auto unused = amrex::ParmParse::getUnusedInputs(nm);
             for (const auto& p: unused) {
                 amrex::Print() << p << "\n";
             }
             amrex::Print() << std::endl;
-            if (abort_on_invalid_params) {
+            if (castro::abort_on_invalid_params) {
                 amrex::Error("Error: invalid parameters");
             }
         }

diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp
@@ -8,6 +8,10 @@
 #include <Gravity.H>
 #endif
 
+#ifdef HYBRID_MOMENTUM
+#include <hybrid.H>
+#endif
+
 #include <ppm.H>
 #include <slope.H>
 
@@ -321,3 +325,143 @@ Castro::mol_diffusive_flux(const Box& bx,
 
   });
 }
+
+
+void
+Castro::compute_flux_from_q(const Box& bx,
+                            Array4<Real const> const& qint,
+                            Array4<Real> const& F,
+                            const int idir) {
+
+  // given a primitive state, compute the flux in direction idir
+  //
+
+  int iu, iv1, iv2;
+  int im1, im2, im3;
+
+  auto coord = geom.Coord();
+  auto mom_check = mom_flux_has_p(idir, idir, coord);
+
+  if (idir == 0) {
+    iu = QU;
+    iv1 = QV;
+    iv2 = QW;
+    im1 = UMX;
+    im2 = UMY;
+    im3 = UMZ;
+
+  } else if (idir == 1) {
+    iu = QV;
+    iv1 = QU;
+    iv2 = QW;
+    im1 = UMY;
+    im2 = UMX;
+    im3 = UMZ;
+
+  } else {
+    iu = QW;
+    iv1 = QU;
+    iv2 = QV;
+    im1 = UMZ;
+    im2 = UMX;
+    im3 = UMY;
+  }
+
+#ifdef HYBRID_MOMENTUM
+  GeometryData geomdata = geom.data();
+#endif
+
+  amrex::ParallelFor(bx,
+  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+  {
+
+    Real u_adv = qint(i,j,k,iu);
+    Real rhoeint = qint(i,j,k,QREINT);
+
+    // Compute fluxes, order as conserved state (not q)
+    F(i,j,k,URHO) = qint(i,j,k,QRHO)*u_adv;
+
+    F(i,j,k,im1) = F(i,j,k,URHO)*qint(i,j,k,iu);
+    if (mom_check) {
+      F(i,j,k,im1) += qint(i,j,k,QPRES);
+    }
+    F(i,j,k,im2) = F(i,j,k,URHO)*qint(i,j,k,iv1);
+    F(i,j,k,im3) = F(i,j,k,URHO)*qint(i,j,k,iv2);
+
+    Real rhoetot = rhoeint + 0.5_rt * qint(i,j,k,QRHO)*
+      (qint(i,j,k,iu)*qint(i,j,k,iu) +
+       qint(i,j,k,iv1)*qint(i,j,k,iv1) +
+       qint(i,j,k,iv2)*qint(i,j,k,iv2));
+
+    F(i,j,k,UEDEN) = u_adv*(rhoetot + qint(i,j,k,QPRES));
+    F(i,j,k,UEINT) = u_adv*rhoeint;
+
+    F(i,j,k,UTEMP) = 0.0;
+#ifdef SHOCK_VAR
+    F(i,j,k,USHK) = 0.0;
+#endif
+
+#ifdef NSE_NET
+    F(i,j,k,UMUP) = 0.0;
+    F(i,j,k,UMUN) = 0.0;
+#endif
+    // passively advected quantities
+    for (int ipassive = 0; ipassive < npassive; ipassive++) {
+      int n  = upassmap(ipassive);
+      int nqp = qpassmap(ipassive);
+
+      F(i,j,k,n) = F(i,j,k,URHO)*qint(i,j,k,nqp);
+    }
+
+#ifdef HYBRID_MOMENTUM
+    // the hybrid routine uses the Godunov indices, not the full NQ state
+    GpuArray<Real, NGDNV> qgdnv_zone;
+    qgdnv_zone[GDRHO] = qint(i,j,k,QRHO);
+    qgdnv_zone[GDU] = qint(i,j,k,QU);
+    qgdnv_zone[GDV] = qint(i,j,k,QV);
+    qgdnv_zone[GDW] = qint(i,j,k,QW);
+    qgdnv_zone[GDPRES] = qint(i,j,k,QPRES);
+    GpuArray<Real, NUM_STATE> F_zone;
+    for (int n = 0; n < NUM_STATE; n++) {
+        F_zone[n] = F(i,j,k,n);
+    }
+    compute_hybrid_flux(qgdnv_zone, geomdata, idir, i, j, k, F_zone);
+    for (int n = 0; n < NUM_STATE; n++) {
+        F(i,j,k,n) = F_zone[n];
+    }
+#endif
+  });
+}
+
+void
+Castro::store_godunov_state(const Box& bx,
+                            Array4<Real const> const& qint,
+#ifdef RADIATION
+                            Array4<Real const> const& lambda,
+#endif
+                            Array4<Real> const& qgdnv) {
+
+  // this copies the full interface state (NQ -- one for each primitive
+  // variable) over to a smaller subset of size NGDNV for use later in the
+  // hydro advancement.
+
+  amrex::ParallelFor(bx,
+  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+  {
+
+
+#ifdef HYBRID_MOMENTUM
+    qgdnv(i,j,k,GDRHO) = qint(i,j,k,QRHO);
+#endif
+    qgdnv(i,j,k,GDU) = qint(i,j,k,QU);
+    qgdnv(i,j,k,GDV) = qint(i,j,k,QV);
+    qgdnv(i,j,k,GDW) = qint(i,j,k,QW);
+    qgdnv(i,j,k,GDPRES) = qint(i,j,k,QPRES);
+#ifdef RADIATION
+    for (int g = 0; g < NGROUPS; g++) {
+      qgdnv(i,j,k,GDLAMS+g) = lambda(i,j,k,g);
+      qgdnv(i,j,k,GDERADS+g) = qint(i,j,k,QRAD+g);
+    }
+#endif
+  });
+}
diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package
@@ -30,7 +30,6 @@ CEXE_sources += riemann.cpp
 CEXE_headers += HLL_solvers.H
 CEXE_headers += riemann_solvers.H
 CEXE_headers += riemann_2shock_solvers.H
-CEXE_sources += riemann_util.cpp
 CEXE_headers += riemann_type.H
 CEXE_headers += slope.H
 CEXE_headers += reconstruction.H