diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 8d8b0b4b7d..a56c7eb1b0 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -176,6 +176,26 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) shk.resize(obx, 1); Elixir elix_shk = shk.elixir(); fab_size += shk.nBytes(); + + // Multidimensional shock detection + // Used for the hybrid Riemann solver + +#ifdef SHOCK_VAR + bool compute_shock = true; +#else + bool compute_shock = false; +#endif + + if (hybrid_riemann == 1 || compute_shock) { +#pragma gpu box(obx) + ca_shock(AMREX_INT_ANYD(obx.loVect()), AMREX_INT_ANYD(obx.hiVect()), + BL_TO_FORTRAN_ANYD(q[mfi]), + BL_TO_FORTRAN_ANYD(shk), + AMREX_REAL_ANYD(dx)); + } + else { + shk.setVal(0.0); + } qxm.resize(obx, NQ); Elixir elix_qxm = qxm.elixir(); @@ -218,7 +238,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) BL_TO_FORTRAN_ANYD(flatn), BL_TO_FORTRAN_ANYD(qaux[mfi]), BL_TO_FORTRAN_ANYD(src_q[mfi]), - BL_TO_FORTRAN_ANYD(shk), BL_TO_FORTRAN_ANYD(dq), BL_TO_FORTRAN_ANYD(qxm), BL_TO_FORTRAN_ANYD(qxp), @@ -277,7 +296,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) BL_TO_FORTRAN_ANYD(flatn), BL_TO_FORTRAN_ANYD(qaux[mfi]), BL_TO_FORTRAN_ANYD(src_q[mfi]), - BL_TO_FORTRAN_ANYD(shk), BL_TO_FORTRAN_ANYD(Ip), BL_TO_FORTRAN_ANYD(Im), BL_TO_FORTRAN_ANYD(Ip_src), diff --git a/Source/hydro/Castro_ctu_nd.F90 b/Source/hydro/Castro_ctu_nd.F90 index ef1edfc3d3..44947c4aa1 100644 --- a/Source/hydro/Castro_ctu_nd.F90 +++ b/Source/hydro/Castro_ctu_nd.F90 @@ -16,7 +16,6 @@ subroutine ctu_ppm_states(lo, hi, & flatn, f_lo, f_hi, & qaux, qa_lo, qa_hi, & srcQ, src_lo, src_hi, & - shk, sk_lo, sk_hi, & Ip, Ip_lo, Ip_hi, & Im, Im_lo, Im_hi, & Ip_src, Ips_lo, Ips_hi, & @@ -58,7 +57,6 @@ subroutine ctu_ppm_states(lo, hi, & #else use trace_ppm_module, only : trace_ppm #endif - use advection_util_module, only : ca_shock use prob_params_module, only : dg implicit none @@ -69,7 +67,6 @@ subroutine ctu_ppm_states(lo, hi, & integer, intent(in) :: f_lo(3), f_hi(3) integer, intent(in) :: qa_lo(3), qa_hi(3) integer, intent(in) :: src_lo(3), src_hi(3) - integer, intent(in) :: sk_lo(3), sk_hi(3) integer, intent(in) :: Ip_lo(3), Ip_hi(3) integer, intent(in) :: Im_lo(3), Im_hi(3) integer, intent(in) :: Ips_lo(3), Ips_hi(3) @@ -100,7 +97,6 @@ subroutine ctu_ppm_states(lo, hi, & real(rt), intent(in) :: flatn(f_lo(1):f_hi(1),f_lo(2):f_hi(2),f_lo(3):f_hi(3)) ! flattening parameter real(rt), intent(in) :: srcQ(src_lo(1):src_hi(1),src_lo(2):src_hi(2),src_lo(3):src_hi(3),NQSRC) ! primitive variable source - real(rt), intent(inout) :: shk(sk_lo(1):sk_hi(1), sk_lo(2):sk_hi(2), sk_lo(3):sk_hi(3)) real(rt), intent(inout) :: Ip(Ip_lo(1):Ip_hi(1),Ip_lo(2):Ip_hi(2),Ip_lo(3):Ip_hi(3),1:3,NQ) real(rt), intent(inout) :: Im(Im_lo(1):Im_hi(1),Im_lo(2):Im_hi(2),Im_lo(3):Im_hi(3),1:3,NQ) real(rt), intent(inout) :: Ip_src(Ips_lo(1):Ips_hi(1),Ips_lo(2):Ips_hi(2),Ips_lo(3):Ips_hi(3),1:3,NQSRC) @@ -130,31 +126,10 @@ subroutine ctu_ppm_states(lo, hi, & logical :: source_nonzero(NQSRC) logical :: reconstruct_state(NQ) - logical :: compute_shock - !$gpu hdt = HALF*dt - ! multidimensional shock detection - -#ifdef SHOCK_VAR - compute_shock = .true. -#else - compute_shock = .false. -#endif - - ! multidimensional shock detection -- this will be used to do the - ! hybrid Riemann solver - if (hybrid_riemann == 1 .or. compute_shock) then - call ca_shock(lo, hi, & - q, qd_lo, qd_hi, & - shk, sk_lo, sk_hi, & - dx) - else - shk(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) = ZERO - endif - ! we don't need to reconstruct all of the NQ state variables, ! depending on how we are tracing reconstruct_state(:) = .true. @@ -385,7 +360,6 @@ subroutine ctu_plm_states(lo, hi, & flatn, f_lo, f_hi, & qaux, qa_lo, qa_hi, & srcQ, src_lo, src_hi, & - shk, sk_lo, sk_hi, & dq, dq_lo, dq_hi, & qxm, qxm_lo, qxm_hi, & qxp, qxp_lo, qxp_hi, & @@ -418,7 +392,6 @@ subroutine ctu_plm_states(lo, hi, & plm_iorder, use_pslope, hybrid_riemann use trace_plm_module, only : trace_plm use slope_module, only : uslope, pslope - use advection_util_module, only : ca_shock use prob_params_module, only : dg implicit none @@ -429,7 +402,6 @@ subroutine ctu_plm_states(lo, hi, & integer, intent(in) :: f_lo(3), f_hi(3) integer, intent(in) :: qa_lo(3), qa_hi(3) integer, intent(in) :: src_lo(3), src_hi(3) - integer, intent(in) :: sk_lo(3), sk_hi(3) integer, intent(in) :: dq_lo(3), dq_hi(3) integer, intent(in) :: qxm_lo(3), qxm_hi(3) integer, intent(in) :: qxp_lo(3), qxp_hi(3) @@ -453,7 +425,6 @@ subroutine ctu_plm_states(lo, hi, & real(rt), intent(in) :: flatn(f_lo(1):f_hi(1),f_lo(2):f_hi(2),f_lo(3):f_hi(3)) ! flattening parameter real(rt), intent(in) :: srcQ(src_lo(1):src_hi(1),src_lo(2):src_hi(2),src_lo(3):src_hi(3),NQSRC) ! primitive variable source - real(rt), intent(inout) :: shk(sk_lo(1):sk_hi(1), sk_lo(2):sk_hi(2), sk_lo(3):sk_hi(3)) real(rt), intent(inout) :: dq(dq_lo(1):dq_hi(1), dq_lo(2):dq_hi(2), dq_lo(3):dq_hi(3), NQ) real(rt), intent(inout) :: qxm(qxm_lo(1):qxm_hi(1), qxm_lo(2):qxm_hi(2), qxm_lo(3):qxm_hi(3), NQ) @@ -474,31 +445,10 @@ subroutine ctu_plm_states(lo, hi, & logical :: reconstruct_state(NQ) - logical :: compute_shock - !$gpu hdt = HALF*dt - ! multidimensional shock detection - -#ifdef SHOCK_VAR - compute_shock = .true. -#else - compute_shock = .false. -#endif - - ! multidimensional shock detection -- this will be used to do the - ! hybrid Riemann solver - if (hybrid_riemann == 1 .or. compute_shock) then - call ca_shock(lo, hi, & - q, qd_lo, qd_hi, & - shk, sk_lo, sk_hi, & - dx) - else - shk(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) = ZERO - endif - ! we don't need to reconstruct all of the NQ state variables, ! depending on how we are tracing reconstruct_state(:) = .true. diff --git a/Source/hydro/Castro_hydro_F.H b/Source/hydro/Castro_hydro_F.H index e639f63ad2..1a5bfcc5b5 100644 --- a/Source/hydro/Castro_hydro_F.H +++ b/Source/hydro/Castro_hydro_F.H @@ -162,6 +162,12 @@ extern "C" (const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(flux)); + void ca_shock + (const int* lo, const int* hi, + const BL_FORT_FAB_ARG_3D(q), + BL_FORT_FAB_ARG_3D(shk), + const amrex::Real* dx); + void ca_ctu_update (const int* lo, const int* hi, const int* is_finest_level, @@ -225,7 +231,6 @@ extern "C" const BL_FORT_FAB_ARG_3D(flatn), const BL_FORT_FAB_ARG_3D(qaux), const BL_FORT_FAB_ARG_3D(srcQ), - const BL_FORT_FAB_ARG_3D(shk), BL_FORT_FAB_ARG_3D(Ip), BL_FORT_FAB_ARG_3D(Im), BL_FORT_FAB_ARG_3D(Ip_src), @@ -257,7 +262,6 @@ extern "C" const BL_FORT_FAB_ARG_3D(flatn), const BL_FORT_FAB_ARG_3D(qaux), const BL_FORT_FAB_ARG_3D(srcQ), - const BL_FORT_FAB_ARG_3D(shk), BL_FORT_FAB_ARG_3D(dq), BL_FORT_FAB_ARG_3D(qxm), BL_FORT_FAB_ARG_3D(qxp), @@ -415,7 +419,6 @@ extern "C" (const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(q), const BL_FORT_FAB_ARG_3D(flatn), - BL_FORT_FAB_ARG_3D(shk), BL_FORT_FAB_ARG_3D(dq), BL_FORT_FAB_ARG_3D(qm), BL_FORT_FAB_ARG_3D(qp), @@ -425,7 +428,6 @@ extern "C" (const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(q), const BL_FORT_FAB_ARG_3D(flatn), - BL_FORT_FAB_ARG_3D(shk), BL_FORT_FAB_ARG_3D(qm), BL_FORT_FAB_ARG_3D(qp), const amrex::Real* dx); diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 8ecd96179c..d52e799e03 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -180,6 +180,26 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) shk.resize(obx, 1); Elixir elix_shk = shk.elixir(); + + // Multidimensional shock detection + // Used for the hybrid Riemann solver + +#ifdef SHOCK_VAR + bool compute_shock = true; +#else + bool compute_shock = false; +#endif + + if (hybrid_riemann == 1 || compute_shock) { +#pragma gpu box(obx) + ca_shock(AMREX_INT_ANYD(obx.loVect()), AMREX_INT_ANYD(obx.hiVect()), + BL_TO_FORTRAN_ANYD(q[mfi]), + BL_TO_FORTRAN_ANYD(shk), + AMREX_REAL_ANYD(dx)); + } + else { + shk.setVal(0.0); + } qm.resize(tbx, NQ*AMREX_SPACEDIM); Elixir elix_qm = qm.elixir(); @@ -199,7 +219,6 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) (AMREX_INT_ANYD(obx.loVect()), AMREX_INT_ANYD(obx.hiVect()), BL_TO_FORTRAN_ANYD(q[mfi]), BL_TO_FORTRAN_ANYD(flatn), - BL_TO_FORTRAN_ANYD(shk), BL_TO_FORTRAN_ANYD(dq), BL_TO_FORTRAN_ANYD(qm), BL_TO_FORTRAN_ANYD(qp), @@ -212,7 +231,6 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) (AMREX_INT_ANYD(obx.loVect()), AMREX_INT_ANYD(obx.hiVect()), BL_TO_FORTRAN_ANYD(q[mfi]), BL_TO_FORTRAN_ANYD(flatn), - BL_TO_FORTRAN_ANYD(shk), BL_TO_FORTRAN_ANYD(qm), BL_TO_FORTRAN_ANYD(qp), AMREX_REAL_ANYD(dx)); diff --git a/Source/hydro/Castro_mol_nd.F90 b/Source/hydro/Castro_mol_nd.F90 index fde63db007..baa19a9877 100644 --- a/Source/hydro/Castro_mol_nd.F90 +++ b/Source/hydro/Castro_mol_nd.F90 @@ -4,7 +4,6 @@ subroutine ca_mol_plm_reconstruct(lo, hi, & q, q_lo, q_hi, & flatn, fl_lo, fl_hi, & - shk, shk_lo, shk_hi, & dq, dq_lo, dq_hi, & qm, qm_lo, qm_hi, & qp, qp_lo, qp_hi, & @@ -12,7 +11,7 @@ subroutine ca_mol_plm_reconstruct(lo, hi, & use amrex_error_module use meth_params_module, only : NQ, NVAR, NGDNV, GDPRES, & - UTEMP, USHK, UMX, & + UTEMP, UMX, & use_flattening, QPRES, & QTEMP, QFS, QFX, QREINT, QRHO, & first_order_hydro, hybrid_riemann, & @@ -24,48 +23,28 @@ subroutine ca_mol_plm_reconstruct(lo, hi, & use eos_module, only : eos use network, only : nspec, naux use prob_params_module, only : dg, coord_type - use advection_util_module, only : ca_shock implicit none integer, intent(in) :: lo(3), hi(3) integer, intent(in) :: q_lo(3), q_hi(3) integer, intent(in) :: fl_lo(3), fl_hi(3) - integer, intent(in) :: shk_lo(3), shk_hi(3) integer, intent(in) :: dq_lo(3), dq_hi(3) integer, intent(in) :: qm_lo(3), qm_hi(3) integer, intent(in) :: qp_lo(3), qp_hi(3) real(rt), intent(inout) :: q(q_lo(1):q_hi(1), q_lo(2):q_hi(2), q_lo(3):q_hi(3), NQ) real(rt), intent(in) :: flatn(fl_lo(1):fl_hi(1), fl_lo(2):fl_hi(2), fl_lo(3):fl_hi(3)) - real(rt), intent(inout) :: shk(shk_lo(1):shk_hi(1), shk_lo(2):shk_hi(2), shk_lo(3):shk_hi(3)) real(rt), intent(inout) :: dq(dq_lo(1):dq_hi(1), dq_lo(2):dq_hi(2), dq_lo(3):dq_hi(3), NQ) real(rt), intent(inout) :: qm(qm_lo(1):qm_hi(1), qm_lo(2):qm_hi(2), qm_lo(3):qm_hi(3), NQ, AMREX_SPACEDIM) real(rt), intent(inout) :: qp(qp_lo(1):qp_hi(1), qp_lo(2):qp_hi(2), qp_lo(3):qp_hi(3), NQ, AMREX_SPACEDIM) real(rt), intent(in) :: dx(3) - integer :: idir, i, j, k, n - logical :: compute_shock type (eos_t) :: eos_state !$gpu -#ifdef SHOCK_VAR - compute_shock = .true. -#else - compute_shock = .false. -#endif - - if (hybrid_riemann == 1 .or. compute_shock) then - call ca_shock(lo, hi, & - q, q_lo, q_hi, & - shk, shk_lo, shk_hi, & - dx) - else - shk(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) = ZERO - endif - do idir = 1, AMREX_SPACEDIM do n = 1, NQ @@ -175,14 +154,13 @@ end subroutine ca_mol_plm_reconstruct subroutine ca_mol_ppm_reconstruct(lo, hi, & q, q_lo, q_hi, & flatn, fl_lo, fl_hi, & - shk, shk_lo, shk_hi, & qm, qm_lo, qm_hi, & qp, qp_lo, qp_hi, & dx) bind(C, name="ca_mol_ppm_reconstruct") use amrex_error_module use meth_params_module, only : NQ, NVAR, NGDNV, GDPRES, & - UTEMP, USHK, UMX, & + UTEMP, UMX, & use_flattening, QPRES, & QTEMP, QFS, QFX, QREINT, QRHO, & first_order_hydro, hybrid_riemann, & @@ -194,46 +172,26 @@ subroutine ca_mol_ppm_reconstruct(lo, hi, & use eos_module, only : eos use network, only : nspec, naux use prob_params_module, only : dg, coord_type - use advection_util_module, only : ca_shock implicit none integer, intent(in) :: lo(3), hi(3) integer, intent(in) :: q_lo(3), q_hi(3) integer, intent(in) :: fl_lo(3), fl_hi(3) - integer, intent(in) :: shk_lo(3), shk_hi(3) integer, intent(in) :: qm_lo(3), qm_hi(3) integer, intent(in) :: qp_lo(3), qp_hi(3) real(rt), intent(inout) :: q(q_lo(1):q_hi(1), q_lo(2):q_hi(2), q_lo(3):q_hi(3), NQ) real(rt), intent(in) :: flatn(fl_lo(1):fl_hi(1), fl_lo(2):fl_hi(2), fl_lo(3):fl_hi(3)) - real(rt), intent(inout) :: shk(shk_lo(1):shk_hi(1), shk_lo(2):shk_hi(2), shk_lo(3):shk_hi(3)) real(rt), intent(inout) :: qm(qm_lo(1):qm_hi(1), qm_lo(2):qm_hi(2), qm_lo(3):qm_hi(3), NQ, AMREX_SPACEDIM) real(rt), intent(inout) :: qp(qp_lo(1):qp_hi(1), qp_lo(2):qp_hi(2), qp_lo(3):qp_hi(3), NQ, AMREX_SPACEDIM) real(rt), intent(in) :: dx(3) integer :: idir, i, j, k, n - logical :: compute_shock type (eos_t) :: eos_state !$gpu -#ifdef SHOCK_VAR - compute_shock = .true. -#else - compute_shock = .false. -#endif - - if (hybrid_riemann == 1 .or. compute_shock) then - call ca_shock(lo, hi, & - q, q_lo, q_hi, & - shk, shk_lo, shk_hi, & - dx) - else - shk(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) = ZERO - endif - - do idir = 1, AMREX_SPACEDIM call ca_ppm_reconstruct(lo, hi, 1, idir, & q, q_lo, q_hi, NQ, 1, NQ, & @@ -310,7 +268,7 @@ subroutine ca_mol_consup(lo, hi, & use amrex_error_module use meth_params_module, only : NQ, NVAR, NGDNV, GDPRES, & - UTEMP, USHK, UMX, & + UTEMP, UMX, & use_flattening, QPRES, & QTEMP, QFS, QFX, QREINT, QRHO, & first_order_hydro, difmag, hybrid_riemann, &