diff --git a/Exec/radiation_tests/RadSphere/Problem_Derive.cpp b/Exec/radiation_tests/RadSphere/Problem_Derive.cpp index 2203ba1a92..f781613a99 100644 --- a/Exec/radiation_tests/RadSphere/Problem_Derive.cpp +++ b/Exec/radiation_tests/RadSphere/Problem_Derive.cpp @@ -2,11 +2,12 @@ #include #include +#include #include #include +#include #include #include -#include using namespace amrex; @@ -23,10 +24,12 @@ void deranalytic(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, auto const der = derfab.array(); GpuArray nugroup = {0.0}; - ca_get_nugroup(nugroup.begin()); - GpuArray 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 index 45e03ba994..39e836b3f0 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -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 index 77c5e56302..3df820de66 100644 --- a/Source/driver/global.H +++ b/Source/driver/global.H @@ -2,10 +2,16 @@ #define GLOBAL_H #include +#ifdef RADIATION +#include +#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 index 0c22df74b8..00ee52923e 100644 --- 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 Erscale = {0.0}; - GpuArray dlognu = {0.0}; - GpuArray 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 index 3e0f97669b..93c5f935d9 100644 --- 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 index ef85a7ba73..a1768d1aa1 100644 --- a/Source/radiation/RAD_F.H +++ b/Source/radiation/RAD_F.H @@ -14,46 +14,6 @@ #include -#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 index fd2a1342ea..6d10bced88 100644 --- 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 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 index 8360c6083a..15ed08f304 100644 --- 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 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 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 deleted file mode 100644 index 2583066a99..0000000000 --- a/Source/radiation/Rad_nd.F90 +++ /dev/null @@ -1,190 +0,0 @@ -subroutine ca_init_fort_constants(hplanck_in, avogadro_in) bind(C, name="ca_init_fort_constants") - - use rad_params_module, only: hplanck, avogadro - use amrex_fort_module, only: rt => amrex_real - - implicit none - - real(rt), intent(in), value :: hplanck_in, avogadro_in - - hplanck = hplanck_in - avogadro = avogadro_in - -end subroutine ca_init_fort_constants - - - -! For single group, let set ngroups to 1. -subroutine ca_initsinglegroup(ngr) bind(C, name="ca_initsinglegroup") - - use rad_params_module, only: ngroups, nugroup, dnugroup, ng0, ng1 - use amrex_fort_module, only: rt => amrex_real - - implicit none - - integer :: ngr - - ! Local variables - integer :: i - - ng0 = 0 - ng1 = 0 - - allocate(nugroup( 0:ngroups-1)) - allocate(dnugroup(0:ngroups-1)) - - do i = 0, ngroups-1 - nugroup(i) = 1.e0_rt ! dummy - dnugroup(i) = 1.e0_rt - enddo - -end subroutine ca_initsinglegroup - -!! ----------------------------------------------------------- -!> @brief This routine is called at problem setup time and is used -!! to initialize the arrays nugroup and dnugroup in -!! probdata with the neutrino group energies and widths. -!! -!! The widths are used to derive neutrino spectrum for plot files -!! ----------------------------------------------------------- -subroutine ca_initgroups(nugr, dnugr, ngr, ngr0, ngr1) - - use rad_params_module, only: ngroups, ng0, ng1, nugroup, dnugroup, & - current_group - use amrex_fort_module, only: rt => amrex_real - - implicit none - - real(rt) :: nugr(0:ngr-1), dnugr(0:ngr-1) - integer :: ngr, ngr0, ngr1 - - ! Local variables - integer :: i - - allocate(current_group, ng0, ng1) - - ng0 = ngr0 - ng1 = ngr1 - - allocate(nugroup( 0:ngroups-1)) - allocate(dnugroup(0:ngroups-1)) - - do i = 0, ngroups-1 - nugroup(i) = nugr(i) - dnugroup(i) = dnugr(i) - enddo - -end subroutine ca_initgroups - -subroutine ca_initgroups2(nugr, dnugr, xnugr, ngr) - - use rad_params_module, only: ngroups, nugroup, dnugroup, xnu, dlognu, lognugroup, & - current_group, ng0, ng1 - use amrex_fort_module, only: rt => amrex_real - - implicit none - - real(rt), intent(in) :: nugr(0:ngr-1), dnugr(0:ngr-1), xnugr(0:ngr) - integer :: ngr - - ! Local variables - integer :: i - - allocate(current_group, ng0, ng1) - - allocate(nugroup( 0:ngroups-1)) - allocate(dnugroup(0:ngroups-1)) - allocate(xnu(0:ngroups)) - allocate(dlognu(0:ngroups-1)) - allocate(lognugroup(0:ngroups-1)) - - nugroup(:) = nugr(:) - dnugroup(:) = dnugr(:) - xnu(:) = xnugr(:) - lognugroup(:) = log(nugroup) - - dlognu(0:ngroups-1) = log(xnu(1:ngroups)) - log(xnu(0:ngroups-1)) - -end subroutine ca_initgroups2 - -subroutine ca_initgroups3(nugr, dnugr, dlognugr, xnugr, ngr, ngr0, ngr1) - ! used by MGFLDSolver - - use rad_params_module, only: ngroups, ng0, ng1, nugroup, dnugroup, & - xnu, dlognu, lognugroup, erg2rhoYe, avogadro, hplanck, & - current_group - use amrex_fort_module, only: rt => amrex_real - - implicit none - - real(rt), intent(in) :: nugr(0:ngr-1), dnugr(0:ngr-1), dlognugr(0:ngr-1), xnugr(0:ngr) - integer :: ngr, ngr0, ngr1 - - ! Local variables - integer :: i - - allocate(current_group, ng0, ng1) - - ng0 = ngr0 - ng1 = ngr1 - - allocate(nugroup( 0:ngroups-1)) - allocate(dnugroup(0:ngroups-1)) - allocate(xnu(0:ngroups)) - allocate(dlognu(0:ngroups-1)) - allocate(erg2rhoYe(0:ngroups-1)) - allocate(lognugroup( 0:ngroups-1)) - - nugroup(:) = nugr(:) - dnugroup(:) = dnugr(:) - xnu(:) = xnugr(:) - dlognu(:) = dlognugr(:) - lognugroup(:) = log(nugroup) - - erg2rhoYe = 0.e0_rt - if (ng0 > 0) then - erg2rhoYe(0:ng0-1) = 1.e0_rt / (avogadro*hplanck*nugroup(0:ng0-1)) - if (ng1 > 0) then - erg2rhoYe(ng0:ng0+ng1-1) = -1.e0_rt / (avogadro*hplanck*nugroup(ng0:ng0+ng1-1)) - end if - end if - -end subroutine ca_initgroups3 - -subroutine ca_get_dlognu(dlognu_out) bind(C, name="ca_get_dlognu") - - use amrex_fort_module, only: rt => amrex_real - use rad_params_module, only: ngroups, dlognu - implicit none - - real(rt), intent(out) :: dlognu_out(0:ngroups-1) - - dlognu_out(:) = dlognu(:) - -end subroutine ca_get_dlognu - -subroutine ca_get_nugroup(nugroup_out) bind(C, name="ca_get_nugroup") - - use amrex_fort_module, only: rt => amrex_real - use rad_params_module, only: ngroups, nugroup - implicit none - - real(rt), intent(out) :: nugroup_out(0:ngroups-1) - - nugroup_out(:) = nugroup(:) - -end subroutine ca_get_nugroup - -subroutine ca_get_dnugroup(dnugroup_out) bind(C, name="ca_get_dnugroup") - - use amrex_fort_module, only: rt => amrex_real - use rad_params_module, only: ngroups, dnugroup - implicit none - - real(rt), intent(out) :: dnugroup_out(0:ngroups-1) - - dnugroup_out(:) = dnugroup(:) - -end subroutine ca_get_dnugroup - - diff --git a/Source/radiation/Radiation.H b/Source/radiation/Radiation.H index 35b4f58f69..58be11061a 100644 --- a/Source/radiation/Radiation.H +++ b/Source/radiation/Radiation.H @@ -885,7 +885,7 @@ public: const amrex::Array& lambda, const amrex::MultiFab& Er, const amrex::MultiFab& F, int iflx); - amrex::Vector xnu, nugroup, dnugroup; + amrex::Vector xnu, nugroup, dnugroup, lognugroup, dlognugroup; protected: diff --git a/Source/radiation/Radiation.cpp b/Source/radiation/Radiation.cpp index 5534bc61b6..f1932f2a74 100644 --- a/Source/radiation/Radiation.cpp +++ b/Source/radiation/Radiation.cpp @@ -287,8 +287,6 @@ Radiation::Radiation(Amr* Parent, Castro* castro, int restart) aRad = 4.*C::sigma_SB / C::c_light; - ca_init_fort_constants(hPlanck, Avogadro); - c = clight; sigma = C::sigma_SB; @@ -460,13 +458,9 @@ Radiation::Radiation(Amr* Parent, Castro* castro, int restart) } if (do_multigroup) { - get_groups(verbose); - } else { - ca_initsinglegroup(nGroups); - // xnu is a dummy for single group xnu.resize(2, 1.0); nugroup.resize(1, 1.0); diff --git a/Source/radiation/rad_params.F90 b/Source/radiation/rad_params.F90 deleted file mode 100644 index 53fd9cca2b..0000000000 --- a/Source/radiation/rad_params.F90 +++ /dev/null @@ -1,23 +0,0 @@ - -! This module stores physical constants and radiation group information -! used for multigroup photon and neutrino radiation diffusion. -! These parameters are initialized in ca_initgroups? to match the -! values used in the C++ radiation code. - -module rad_params_module - - ! radiation energy group information - - use amrex_fort_module, only: rt => amrex_real - - implicit none - - integer, parameter :: ngroups = NGROUPS - - real :: hplanck, avogadro - - integer, allocatable, save :: current_group, ng0, ng1 - real(rt), save, allocatable :: nugroup(:), dnugroup(:), xnu(:), dlognu(:), & - erg2rhoYe(:), lognugroup(:) - -end module rad_params_module