Move dlognugroup to C++ (#2637)

AMReX-Astro · Oct 30, 2023 · ebd2255 · ebd2255
1 parent 13f63f2
commit ebd2255
Show file tree

Hide file tree

Showing 12 changed files with 30 additions and 294 deletions.
diff --git a/Exec/radiation_tests/RadSphere/Problem_Derive.cpp b/Exec/radiation_tests/RadSphere/Problem_Derive.cpp
@@ -2,11 +2,12 @@
 
 #include <Derive.H>
 #include <Castro.H>
+#include <Radiation.H>
 #include <prob_parameters.H>
 #include <extern_parameters.H>
+#include <global.H>
 #include <gravity_params.H>
 #include <problem_util.H>
-#include <RAD_F.H>
 
 using namespace amrex;
 
@@ -23,10 +24,12 @@ void deranalytic(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
   auto const der = derfab.array();
 
   GpuArray<Real, NGROUPS> nugroup = {0.0};
-  ca_get_nugroup(nugroup.begin());
-
   GpuArray<Real, NGROUPS> dnugroup = {0.0};
-  ca_get_dnugroup(dnugroup.begin());
+
+  for (int i = 0; i < NGROUPS; ++i) {
+      nugroup[i] = global::the_radiation_ptr->nugroup[i];
+      dnugroup[i] = global::the_radiation_ptr->dnugroup[i];
+  }
 
   amrex::ParallelFor(bx,
   [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)

diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
@@ -15,6 +15,7 @@
 #include <AMReX_Utility.H>
 #include <AMReX_CONSTANTS.H>
 #include <Castro.H>
+#include <global.H>
 #include <runtime_parameters.H>
 #include <AMReX_VisMF.H>
 #include <AMReX_TagBox.H>
@@ -729,6 +730,7 @@ Castro::Castro (Amr&            papa,
       if (radiation == nullptr) {
         // radiation is a static object, only alloc if not already there
         radiation = new Radiation(parent, this);
+        global::the_radiation_ptr = radiation;
       }
       radiation->regrid(level, grids, dmap);
 

diff --git a/Source/driver/global.H b/Source/driver/global.H
@@ -2,10 +2,16 @@
 #define GLOBAL_H
 
 #include <AMReX_Amr.H>
+#ifdef RADIATION
+#include <Radiation.H>
+#endif
 
 namespace global
 {
     AMREX_INLINE Amr* the_amr_ptr;
+#ifdef RADIATION
+    AMREX_INLINE Radiation* the_radiation_ptr;
+#endif
 }
 
 #endif
diff --git a/Source/hydro/Castro_ctu_rad.cpp b/Source/hydro/Castro_ctu_rad.cpp
@@ -36,24 +36,19 @@ Castro::ctu_rad_consup(const Box& bx,
   const int coord_type = geom.Coord();
 
   GpuArray<Real, NGROUPS> Erscale = {0.0};
-
   GpuArray<Real, NGROUPS> dlognu = {0.0};
 
-  GpuArray<Real, NGROUPS> nugroup = {0.0};
-
-
-
   if (NGROUPS > 1) {
-    ca_get_nugroup(nugroup.begin());
-    ca_get_dlognu(dlognu.begin());
-
+    for (int g = 0; g < NGROUPS; ++g) {
+      dlognu[g] = radiation->dlognugroup[g];
+    }
     if (radiation::fspace_advection_type == 1) {
       for (int g = 0; g < NGROUPS; g++) {
         Erscale[g] = dlognu[g];
       }
     } else {
       for (int g = 0; g < NGROUPS; g++) {
-        Erscale[g] = nugroup[g] * dlognu[g];
+        Erscale[g] = radiation->nugroup[g] * dlognu[g];
       }
     }
   }

diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package
@@ -33,8 +33,6 @@ FEXE_headers += RAD_F.H
 
 CEXE_sources += trace_ppm_rad.cpp
 
-ca_F90EXE_sources += rad_params.F90
-ca_F90EXE_sources += Rad_nd.F90
 CEXE_headers += fluxlimiter.H
 CEXE_headers += RadHydro.H
 

diff --git a/Source/radiation/RAD_F.H b/Source/radiation/RAD_F.H
@@ -14,46 +14,6 @@
 
 #include <AMReX_BLFort.H>
 
-#ifdef __cplusplus
-extern "C"
-{
-#endif
-
-void ca_init_fort_constants(const amrex::Real hplanck, const amrex::Real avogadro);
-
-void ca_initsinglegroup
-  (const int& ngroups);
-
-void ca_get_dlognu
-  (amrex::Real* dlognu_out);
-
-void ca_get_nugroup
-  (amrex::Real* nugroup);
-
-void ca_get_dnugroup
-  (amrex::Real* dnugroup);
-
-#ifdef __cplusplus
-}
-#endif
-
-
-BL_FORT_PROC_DECL(CA_INITGROUPS,ca_initgroups)
-  (const amrex::Real* nugroup, const amrex::Real* dnugroup,
-   const int& ngroups, const int& ng0, const int& ng1);
-
-BL_FORT_PROC_DECL(CA_INITGROUPS2,ca_initgroups2)
-  (const amrex::Real* nugroup, const amrex::Real* dnugroup,
-   const amrex::Real* xnu, const int& ngroups);
-
-BL_FORT_PROC_DECL(CA_INITGROUPS3,ca_initgroups3)
-  (const amrex::Real* nugroup, const amrex::Real* dnugroup, const amrex::Real* dlognugroup,
-   const amrex::Real* xnu, const int& ngroups, const int& ng0, const int& ng1);
-
-BL_FORT_PROC_DECL(CA_COMPUTE_KAPKAP, ca_compute_kapkap)
-     (BL_FORT_FAB_ARG(kapkap), const BL_FORT_FAB_ARG(kap_r)); 
-
-
 #ifdef __cplusplus
 extern "C" {
 #endif

diff --git a/Source/radiation/RadMultiGroup.cpp b/Source/radiation/RadMultiGroup.cpp
@@ -21,11 +21,10 @@ void Radiation::get_groups(int verbose)
   group_print_factor   = 1.0;
   group_units          = " (units are Hz)";
 
-  Vector<Real> dlognugroup;
-
   xnu.resize(nGroups+1, 0.0); // Bounds of the frequency group
   nugroup.resize(nGroups, 1.0); // Geometric center of the frequency group
   dnugroup.resize(nGroups, 0.0); // Width of the frequency group
+  lognugroup.resize(nGroups, 0.0); // Log of the center of the frequency group
   dlognugroup.resize(nGroups, 0.0); // Log of the width of the frequency group
 
   ParmParse pp("radiation");
@@ -79,22 +78,10 @@ void Radiation::get_groups(int verbose)
         dlognugroup[i] = log(xnu[i+1]) - log(xnu[i]);
       }
     }
-  }
-
-  int nG0 = 0, nG1 = 0;
 
-  if (SolverType == MGFLDSolver) { 
-    BL_FORT_PROC_CALL(CA_INITGROUPS3,ca_initgroups3)
-      (nugroup.dataPtr(), dnugroup.dataPtr(), dlognugroup.dataPtr(), xnu.dataPtr(), 
-       nGroups, nG0, nG1);
-  }
-  else if (xnu.size() > 0) {
-    BL_FORT_PROC_CALL(CA_INITGROUPS2,ca_initgroups2)
-      (nugroup.dataPtr(), dnugroup.dataPtr(), xnu.dataPtr(), nGroups);
-  }
-  else {
-    BL_FORT_PROC_CALL(CA_INITGROUPS,ca_initgroups)
-      (nugroup.dataPtr(), dnugroup.dataPtr(), nGroups, nG0, nG1);
+    for (int i = 0; i < nGroups; ++i) {
+        lognugroup[i] = std::log(nugroup[i]);
+    }
   }
 
   if (ParallelDescriptor::IOProcessor()) {

diff --git a/Source/radiation/RadPlotvar.cpp b/Source/radiation/RadPlotvar.cpp
@@ -72,7 +72,9 @@ void Radiation::save_lab_Er_in_plotvar(int level, const MultiFab& Snew,
     GpuArray<Real, NGROUPS> dlognu = {0.0};
 
     if (NGROUPS > 1) {
-        ca_get_dlognu(dlognu.begin());
+        for (int i = 0; i < NGROUPS; ++i) {
+            dlognu[i] = dlognugroup[i];
+        }
     }
 
 #ifdef _OPENMP
@@ -155,7 +157,9 @@ void Radiation::save_flux_in_plotvar(int level, const MultiFab& Snew,
     GpuArray<Real, NGROUPS> dlognu = {0.0};
 
     if (NGROUPS > 1) {
-        ca_get_dlognu(dlognu.begin());
+        for (int i = 0; i < NGROUPS; ++i) {
+            dlognu[i] = dlognugroup[i];
+        }
     }
 
 #ifdef _OPENMP

diff --git a/Source/radiation/Rad_nd.F90 b/Source/radiation/Rad_nd.F90