From 6e9b1cfeac4800f1435faa37e3d75236e1c3e75f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 May 2020 12:01:53 -0400 Subject: [PATCH 1/2] remove old headers --- Source/mhd/Castro_mhd_F.H | 47 --------------------------------------- 1 file changed, 47 deletions(-) diff --git a/Source/mhd/Castro_mhd_F.H b/Source/mhd/Castro_mhd_F.H index e0484d706e..e781a29616 100644 --- a/Source/mhd/Castro_mhd_F.H +++ b/Source/mhd/Castro_mhd_F.H @@ -9,11 +9,6 @@ extern "C" { #endif - void mhd_flatten - (const int lo[], const int hi[], - const BL_FORT_FAB_ARG_3D(q), - BL_FORT_FAB_ARG_3D(flatn)); - void plm (const int lo[], const int hi[], const int idir, @@ -28,48 +23,6 @@ extern "C" const BL_FORT_FAB_ARG_3D(srcq), const amrex::Real dx[], const amrex::Real dt); - void electric_edge_x(const int lo[], const int hi[], - const BL_FORT_FAB_ARG_3D(q), - BL_FORT_FAB_ARG_3D(E), - BL_FORT_FAB_ARG_3D(flxy), - BL_FORT_FAB_ARG_3D(flxz)); - - void electric_edge_y(const int lo[], const int hi[], - const BL_FORT_FAB_ARG_3D(q), - BL_FORT_FAB_ARG_3D(E), - BL_FORT_FAB_ARG_3D(flxx), - BL_FORT_FAB_ARG_3D(flxz)); - - void electric_edge_z(const int lo[], const int hi[], - const BL_FORT_FAB_ARG_3D(q), - BL_FORT_FAB_ARG_3D(E), - BL_FORT_FAB_ARG_3D(flxx), - BL_FORT_FAB_ARG_3D(flxy)); - - void corner_couple(const int lo[], const int hi[], - BL_FORT_FAB_ARG_3D(qr_out), - BL_FORT_FAB_ARG_3D(ql_out), - const BL_FORT_FAB_ARG_3D(ur), - const BL_FORT_FAB_ARG_3D(ul), - const BL_FORT_FAB_ARG_3D(flxd2), - const BL_FORT_FAB_ARG_3D(Ed1), - const BL_FORT_FAB_ARG_3D(Ed3), - const int d1, const int d2, const int d3, - const Real dx, const Real dt); - - void half_step(const int lo[], const int hi[], - BL_FORT_FAB_ARG_3D(qr_out), - BL_FORT_FAB_ARG_3D(ql_out), - const BL_FORT_FAB_ARG_3D(ur), - const BL_FORT_FAB_ARG_3D(ul), - const BL_FORT_FAB_ARG_3D(flxd1), - const BL_FORT_FAB_ARG_3D(flxd2), - const BL_FORT_FAB_ARG_3D(Ed), - const BL_FORT_FAB_ARG_3D(Ed1), - const BL_FORT_FAB_ARG_3D(Ed2), - const int d, const int d1, const int d2, - const Real dx, const Real dt); - void hlld(const int lo[], const int hi[], const BL_FORT_FAB_ARG_3D(qleft), const BL_FORT_FAB_ARG_3D(qright), From ad01d7cc398c4db56ab7a75e495d24080dd7f12a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 May 2020 15:41:43 -0400 Subject: [PATCH 2/2] Port MHD flattening to C++ and address box sizes (#984) Now the flattening is done over the same size box as we do reconstruction on. Previously we were not computing the flattening coefficient in all the zones where we applied it. This also reuses the C++ flattening routine. Note: this adds QPTOT to the primitive variable state so we can flattening on that. --- Source/driver/_variables | 1 + Source/hydro/advection_util.cpp | 8 +- Source/mhd/Castro_mhd.cpp | 127 ++++++++---------- Source/mhd/Make.package | 1 - Source/mhd/flatten_mhd.F90 | 231 -------------------------------- 5 files changed, 67 insertions(+), 301 deletions(-) delete mode 100644 Source/mhd/flatten_mhd.F90 diff --git a/Source/driver/_variables b/Source/driver/_variables index 97133c394c..6244512278 100644 --- a/Source/driver/_variables +++ b/Source/driver/_variables @@ -39,6 +39,7 @@ species QFS (NQSRC, PRIM_SPECIES_HAVE_SOURCES) (nspec, NumSpec) None auxiliary QFX (NQSRC, PRIM_SPECIES_HAVE_SOURCES) (naux, NumAux) None total-pressure QPTOT None 1 RADIATION + total-pressure QPTOT None 1 MHD total-reint QREITOT None 1 RADIATION radiation QRAD None (ngroups, NGROUPS) RADIATION diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index b63d78c9e1..dbaf946eb4 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -61,7 +61,6 @@ Castro::ctoprim(const Box& bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept { - #ifndef AMREX_USE_CUDA if (uin(i,j,k,URHO) <= 0.0_rt) { std::cout << std::endl; @@ -162,6 +161,13 @@ Castro::ctoprim(const Box& bx, q_arr(i,j,k,QGC) = eos_state.gam1; #endif +#ifdef MHD + q_arr(i,j,k,QPTOT) = q_arr(i,j,k,QPRES) + + 0.5_rt * (q_arr(i,j,k,QMAGX) * q_arr(i,j,k,QMAGX) + + q_arr(i,j,k,QMAGY) * q_arr(i,j,k,QMAGY) + + q_arr(i,j,k,QMAGZ) * q_arr(i,j,k,QMAGZ)); +#endif + #ifdef RADIATION qaux_arr(i,j,k,QGAMCG) = eos_state.gam1; qaux_arr(i,j,k,QCG) = eos_state.cs; diff --git a/Source/mhd/Castro_mhd.cpp b/Source/mhd/Castro_mhd.cpp index 84e509d184..c457c70ff1 100644 --- a/Source/mhd/Castro_mhd.cpp +++ b/Source/mhd/Castro_mhd.cpp @@ -49,13 +49,10 @@ Castro::just_the_mhd(Real time, Real dt) FArrayBox srcQ; FArrayBox flatn; + FArrayBox flatg; - FArrayBox qx_left; - FArrayBox qx_right; - FArrayBox qy_left; - FArrayBox qy_right; - FArrayBox qz_left; - FArrayBox qz_right; + FArrayBox qleft[AMREX_SPACEDIM]; + FArrayBox qright[AMREX_SPACEDIM]; FArrayBox flxx1D; FArrayBox flxy1D; @@ -197,87 +194,81 @@ Castro::just_the_mhd(Real time, Real dt) check_for_mhd_cfl_violation(bx, dt, q_arr, qaux_arr); - flatn.resize(bx_gc, 1); + // we need to compute the flattening coefficient for every zone + // center where we do reconstruction + + const Box& bxi = amrex::grow(bx, IntVect(3, 3, 3)); + + flatn.resize(bxi, 1); auto flatn_arr = flatn.array(); auto elix_flatn = flatn.elixir(); - amrex::ParallelFor(bx_gc, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept - { - flatn_arr(i,j,k) = 0.0; - }); + flatg.resize(bxi, 1); + auto flatg_arr = flatg.array(); + auto elix_flatg = flatg.elixir(); - mhd_flatten(lo1, hi1, - BL_TO_FORTRAN_ANYD(q), - BL_TO_FORTRAN_ANYD(flatn)); + if (use_flattening == 0) { + amrex::ParallelFor(bxi, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + flatn_arr(i,j,k) = 0.0; + }); + } else { - // Interpolate Cell centered values to faces - qx_left.resize(bx_gc, NQ); - auto qx_left_arr = qx_left.array(); - auto elix_qx_left = qx_left.elixir(); + uflatten(bxi, q_arr, flatn_arr, QPRES); + uflatten(bxi, q_arr, flatg_arr, QPTOT); - qx_right.resize(bx_gc, NQ); - auto qx_right_arr = qx_right.array(); - auto elix_qx_right = qx_right.elixir(); + amrex::ParallelFor(obx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + flatn_arr(i,j,k) = flatn_arr(i,j,k) * flatg_arr(i,j,k); + }); - qy_left.resize(bx_gc, NQ); - auto qy_left_arr = qy_left.array(); - auto elix_qy_left = qy_left.elixir(); + } - qy_right.resize(bx_gc, NQ); - auto qy_right_arr = qy_right.array(); - auto elix_qy_right = qy_right.elixir(); + // Interpolate Cell centered values to faces + qleft[0].resize(bx_gc, NQ); + auto qx_left_arr = qleft[0].array(); + auto elix_qx_left = qleft[0].elixir(); - qz_left.resize(bx_gc, NQ); - auto qz_left_arr = qz_left.array(); - auto elix_qz_left = qz_left.elixir(); + qright[0].resize(bx_gc, NQ); + auto qx_right_arr = qright[0].array(); + auto elix_qx_right = qright[0].elixir(); - qz_right.resize(bx_gc, NQ); - auto qz_right_arr = qz_right.array(); - auto elix_qz_right = qz_right.elixir(); + qleft[1].resize(bx_gc, NQ); + auto qy_left_arr = qleft[1].array(); + auto elix_qy_left = qleft[1].elixir(); - const Box& nbxi = amrex::grow(bx, IntVect(3, 3, 3)); + qright[1].resize(bx_gc, NQ); + auto qy_right_arr = qright[1].array(); + auto elix_qy_right = qright[1].elixir(); - plm(nbxi.loVect(), nbxi.hiVect(), 1, - BL_TO_FORTRAN_ANYD(q), - BL_TO_FORTRAN_ANYD(qaux), - BL_TO_FORTRAN_ANYD(flatn), - BL_TO_FORTRAN_ANYD(Bx), - BL_TO_FORTRAN_ANYD(By), - BL_TO_FORTRAN_ANYD(Bz), - BL_TO_FORTRAN_ANYD(qx_left), - BL_TO_FORTRAN_ANYD(qx_right), - BL_TO_FORTRAN_ANYD(srcQ), - dx_f, dt); + qleft[2].resize(bx_gc, NQ); + auto qz_left_arr = qleft[2].array(); + auto elix_qz_left = qleft[2].elixir(); - const Box& nbyi = amrex::grow(bx, IntVect(3, 3, 3)); + qright[2].resize(bx_gc, NQ); + auto qz_right_arr = qright[2].array(); + auto elix_qz_right = qright[2].elixir(); - plm(nbyi.loVect(), nbyi.hiVect(), 2, - BL_TO_FORTRAN_ANYD(q), - BL_TO_FORTRAN_ANYD(qaux), - BL_TO_FORTRAN_ANYD(flatn), - BL_TO_FORTRAN_ANYD(Bx), - BL_TO_FORTRAN_ANYD(By), - BL_TO_FORTRAN_ANYD(Bz), - BL_TO_FORTRAN_ANYD(qy_left), - BL_TO_FORTRAN_ANYD(qy_right), - BL_TO_FORTRAN_ANYD(srcQ), - dx_f, dt); - const Box& nbzi = amrex::grow(bx, IntVect(3, 3, 3)); + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + + const int idir_f = idir + 1; - plm(nbzi.loVect(), nbzi.hiVect(), 3, + plm(bxi.loVect(), bxi.hiVect(), idir_f, BL_TO_FORTRAN_ANYD(q), BL_TO_FORTRAN_ANYD(qaux), BL_TO_FORTRAN_ANYD(flatn), BL_TO_FORTRAN_ANYD(Bx), BL_TO_FORTRAN_ANYD(By), BL_TO_FORTRAN_ANYD(Bz), - BL_TO_FORTRAN_ANYD(qz_left), - BL_TO_FORTRAN_ANYD(qz_right), + BL_TO_FORTRAN_ANYD(qleft[idir]), + BL_TO_FORTRAN_ANYD(qright[idir]), BL_TO_FORTRAN_ANYD(srcQ), dx_f, dt); + } // Corner Couple and find the correct fluxes + electric fields @@ -296,8 +287,8 @@ Castro::just_the_mhd(Real time, Real dt) auto elix_flxx1D = flxx1D.elixir(); hlld(bfx.loVect(), bfx.hiVect(), - BL_TO_FORTRAN_ANYD(qx_left), - BL_TO_FORTRAN_ANYD(qx_right), + BL_TO_FORTRAN_ANYD(qleft[0]), + BL_TO_FORTRAN_ANYD(qright[0]), BL_TO_FORTRAN_ANYD(flxx1D), 1); // y-dir @@ -309,8 +300,8 @@ Castro::just_the_mhd(Real time, Real dt) auto elix_flxy1D = flxy1D.elixir(); hlld(bfy.loVect(), bfy.hiVect(), - BL_TO_FORTRAN_ANYD(qy_left), - BL_TO_FORTRAN_ANYD(qy_right), + BL_TO_FORTRAN_ANYD(qleft[1]), + BL_TO_FORTRAN_ANYD(qright[1]), BL_TO_FORTRAN_ANYD(flxy1D), 2); // z-dir @@ -322,8 +313,8 @@ Castro::just_the_mhd(Real time, Real dt) auto elix_flxz1D = flxz1D.elixir(); hlld(bfz.loVect(), bfz.hiVect(), - BL_TO_FORTRAN_ANYD(qz_left), - BL_TO_FORTRAN_ANYD(qz_right), + BL_TO_FORTRAN_ANYD(qleft[2]), + BL_TO_FORTRAN_ANYD(qright[2]), BL_TO_FORTRAN_ANYD(flxz1D), 3); diff --git a/Source/mhd/Make.package b/Source/mhd/Make.package index 489ae98acd..ad5121736f 100644 --- a/Source/mhd/Make.package +++ b/Source/mhd/Make.package @@ -7,7 +7,6 @@ CEXE_sources += Castro_mhd.cpp CEXE_sources += mhd_util.cpp CEXE_headers += mhd_util.H -ca_F90EXE_sources += flatten_mhd.F90 CEXE_sources += ct_upwind.cpp CEXE_sources += electric.cpp ca_f90EXE_sources += hlld.f90 diff --git a/Source/mhd/flatten_mhd.F90 b/Source/mhd/flatten_mhd.F90 deleted file mode 100644 index 4129a67773..0000000000 --- a/Source/mhd/flatten_mhd.F90 +++ /dev/null @@ -1,231 +0,0 @@ -module flatten_module_mhd - - use amrex_mempool_module, only : bl_allocate, bl_deallocate - use amrex_constants_module, only : ZERO - - use amrex_fort_module, only : rt => amrex_real - implicit none - - private - - public :: uflatten -#ifdef MHD - public :: mhd_flatten -#endif -contains - -! ::: -! ::: ------------------------------------------------------------------ -! ::: - - subroutine uflatten(lo, hi, q, q_lo, q_hi, flatn, f_lo, f_hi, ipres) - - ! here, ipres is the pressure variable we want to consider jumps on - ! passing it in allows - use meth_params_module, only : small_pres, QU, QV, QW, NQ - use prob_params_module, only : dg - use amrex_constants_module - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: q_lo(3), q_hi(3) - integer, intent(in) :: f_lo(3), f_hi(3) - integer, intent(in) :: ipres - - real(rt) , intent(in) :: q(q_lo(1):q_hi(1),q_lo(2):q_hi(2),q_lo(3):q_hi(3),NQ+1) - real(rt) , intent(inout) :: flatn(f_lo(1):f_hi(1),f_lo(2):f_hi(2),f_lo(3):f_hi(3)) - - integer :: i, j, k, ishft - - real(rt) :: denom, zeta, tst, tmp, ftmp - - ! Local arrays - real(rt) , pointer :: dp(:,:,:), z(:,:,:), chi(:,:,:) - - ! Knobs for detection of strong shock - real(rt) , parameter :: shktst = 0.33e0_rt, zcut1 = 0.75e0_rt, zcut2 = 0.85e0_rt, dzcut = ONE/(zcut2-zcut1) - - call bl_allocate(dp ,lo(1)-1,hi(1)+1,lo(2)-1,hi(2)+1,lo(3)-1,hi(3)+1) - call bl_allocate(z ,lo(1)-1,hi(1)+1,lo(2)-1,hi(2)+1,lo(3)-1,hi(3)+1) - call bl_allocate(chi,lo(1)-1,hi(1)+1,lo(2)-1,hi(2)+1,lo(3)-1,hi(3)+1) - - ! x-direction flattening coef - do k = lo(3), hi(3) - do j = lo(2), hi(2) - !dir$ ivdep - do i = lo(1)-1*dg(1),hi(1)+1*dg(1) - dp(i,j,k) = q(i+1*dg(1),j,k,ipres) - q(i-1*dg(1),j,k,ipres) - denom = max(small_pres, abs(q(i+2*dg(1),j,k,ipres) - q(i-2*dg(1),j,k,ipres))) - zeta = abs(dp(i,j,k))/denom - z(i,j,k) = min( ONE, max( ZERO, dzcut*(zeta - zcut1) ) ) - if (q(i-1*dg(1),j,k,QU) - q(i+1*dg(1),j,k,QU) >= ZERO) then - tst = ONE - else - tst = ZERO - endif - tmp = min(q(i+1*dg(1),j,k,ipres), q(i-1*dg(1),j,k,ipres)) - if ((abs(dp(i,j,k))/tmp) > shktst) then - chi(i,j,k) = tst - else - chi(i,j,k) = ZERO - endif - enddo - do i = lo(1), hi(1) - if(dp(i,j,k) > ZERO)then - ishft = 1 - else - ishft = -1 - endif - flatn(i,j,k) = ONE - & - max(chi(i-ishft*dg(1),j,k)*z(i-ishft*dg(1),j,k), & - chi(i,j,k)*z(i,j,k)) - enddo - enddo - enddo - - ! y-direction flattening coef - do k = lo(3), hi(3) - do j = lo(2)-1*dg(2), hi(2)+1*dg(2) - !dir$ ivdep - do i = lo(1), hi(1) - dp(i,j,k) = q(i,j+1*dg(2),k,ipres) - q(i,j-1*dg(2),k,ipres) - denom = max(small_pres, abs(q(i,j+2*dg(2),k,ipres) - q(i,j-2*dg(2),k,ipres))) - zeta = abs(dp(i,j,k))/denom - z(i,j,k) = min( ONE, max( ZERO, dzcut*(zeta - zcut1) ) ) - if (q(i,j-1*dg(2),k,QV) - q(i,j+1*dg(2),k,QV) >= ZERO) then - tst = ONE - else - tst = ZERO - endif - tmp = min(q(i,j+1*dg(2),k,ipres), q(i,j-1*dg(2),k,ipres)) - if ((abs(dp(i,j,k))/tmp) > shktst) then - chi(i,j,k) = tst - else - chi(i,j,k) = ZERO - endif - enddo - end do - do j = lo(2), hi(2) - do i = lo(1), hi(1) - if(dp(i,j,k) > ZERO)then - ishft = 1 - else - ishft = -1 - endif - ftmp = ONE - & - max(chi(i,j-ishft*dg(2),k)*z(i,j-ishft*dg(2),k), & - chi(i,j,k)*z(i,j,k)) - flatn(i,j,k) = min( flatn(i,j,k), ftmp ) - enddo - enddo - enddo - - ! z-direction flattening coef - do k = lo(3)-1*dg(3), hi(3)+1*dg(3) - do j = lo(2), hi(2) - !dir$ ivdep - do i = lo(1), hi(1) - dp(i,j,k) = q(i,j,k+1*dg(3),ipres) - q(i,j,k-1*dg(3),ipres) - denom = max(small_pres, abs(q(i,j,k+2*dg(3),ipres) - q(i,j,k-2*dg(3),ipres))) - zeta = abs(dp(i,j,k))/denom - z(i,j,k) = min( ONE, max( ZERO, dzcut*(zeta - zcut1) ) ) - if (q(i,j,k-1*dg(3),QW) - q(i,j,k+1*dg(3),QW) >= ZERO) then - tst = ONE - else - tst = ZERO - endif - tmp = min(q(i,j,k+1*dg(3),ipres), q(i,j,k-1*dg(3),ipres)) - if ((abs(dp(i,j,k))/tmp) > shktst) then - chi(i,j,k) = tst - else - chi(i,j,k) = ZERO - endif - enddo - enddo - enddo - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - if(dp(i,j,k) > ZERO)then - ishft = 1 - else - ishft = -1 - endif - ftmp = ONE - & - max(chi(i,j,k-ishft*dg(3))*z(i,j,k-ishft*dg(3)), & - chi(i,j,k)*z(i,j,k)) - flatn(i,j,k) = min( flatn(i,j,k), ftmp ) - enddo - enddo - enddo - - call bl_deallocate(dp ) - call bl_deallocate(z ) - call bl_deallocate(chi) - - end subroutine uflatten - -#ifdef MHD - subroutine mhd_flatten(lo, hi, & - q, q_lo, q_hi, & - flatn, f_lo, f_hi) bind(C, name="mhd_flatten") - - use meth_params_module, only : QPRES, QU, QV, QW, NQ, QMAGX, QMAGY, QMAGZ - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: q_lo(3), q_hi(3) - integer, intent(in) :: f_lo(3), f_hi(3) - - real(rt) , intent(in) :: q(q_lo(1):q_hi(1),q_lo(2):q_hi(2),q_lo(3):q_hi(3),NQ) - real(rt) , intent(inout) :: flatn(f_lo(1):f_hi(1),f_lo(2):f_hi(2),f_lo(3):f_hi(3)) - - integer :: i, j, k, QPTOT - - real(rt) :: q2(q_lo(1):q_hi(1), q_lo(2):q_hi(2), q_lo(3):q_hi(3), NQ+1) - real(rt) , pointer :: flatg(:,:,:) - - - QPTOT = NQ+1 - - q2(:,:,:,1:NQ) = q(:,:,:,:) - q2(:,:,:,QPTOT) = q(:,:,:,QPRES) + 0.5d0 * ( q(:,:,:,QMAGX)**2 + q(:,:,:,QMAGY)**2 + q(:,:,:,QMAGZ)**2 ) - - call uflatten(lo, hi, q2, q_lo, q_hi, flatn, f_lo, f_hi, QPTOT) - - call bl_allocate(flatg, f_lo, f_hi) - call uflatten(lo, hi, q2, q_lo, q_hi, flatg, f_lo, f_hi, QPRES) - - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - - flatn(i,j,k) = flatn(i,j,k) * flatg(i,j,k) - -! if (flatten_pp_threshold > ZERO) then -! if ( q(i-1,j,k,QU) + q(i,j-1,k,QV) + q(i,j,k-1,QW) > & -! q(i+1,j,k,QU) + q(i,j+1,k,QV) + q(i,j,k+1,QW) ) then - -! if (q(i,j,k,QPRES) < flatten_pp_threshold * q(i,j,k,QPTOT)) then -! flatn(i,j,k) = ZERO -! end if - -! end if -! endif - - end do - end do - end do - - call bl_deallocate(flatg) - - end subroutine mhd_flatten -#endif - - - -end module flatten_module_mhd