move the riemann_util functions into Castro_mol.cpp

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
+  });
+}

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