diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package index 7e15b5b645..8a13baf408 100644 --- a/Source/radiation/Make.package +++ b/Source/radiation/Make.package @@ -41,8 +41,7 @@ ca_F90EXE_sources += rad_params.F90 ca_F90EXE_sources += Rad_nd.F90 CEXE_headers += fluxlimiter.H CEXE_headers += RadHydro.H -ca_F90EXE_sources += fluxlimiter.F90 -ca_F90EXE_sources += filter.F90 + CEXE_sources += RadDerive.cpp CEXE_headers += RadDerive.H CEXE_headers += rad_util.H diff --git a/Source/radiation/RAD_F.H b/Source/radiation/RAD_F.H index b1c28c0a07..6c3f1e4225 100644 --- a/Source/radiation/RAD_F.H +++ b/Source/radiation/RAD_F.H @@ -53,15 +53,6 @@ BL_FORT_PROC_DECL(CA_INITGROUPS3,ca_initgroups3) 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 -void ca_initfluxlimiter - (const int* limiter, const int* closure); -#ifdef __cplusplus -} -#endif extern "C" { diff --git a/Source/radiation/Radiation.cpp b/Source/radiation/Radiation.cpp index ace2a6969b..d2eb374b1b 100644 --- a/Source/radiation/Radiation.cpp +++ b/Source/radiation/Radiation.cpp @@ -333,8 +333,6 @@ Radiation::Radiation(Amr* Parent, Castro* castro, int restart) amrex::Abort("MGFLDSolver does not support limiter = 1"); } - ca_initfluxlimiter(&radiation::limiter, &radiation::closure); - inner_update_limiter = 0; pp.query("inner_update_limiter", inner_update_limiter); diff --git a/Source/radiation/filter.F90 b/Source/radiation/filter.F90 deleted file mode 100644 index f2f26ce47d..0000000000 --- a/Source/radiation/filter.F90 +++ /dev/null @@ -1,124 +0,0 @@ -! reference: R. J. Purser (J. of Clim. and Apld. Meteorology, 1987) - -module filter_module - - use amrex_fort_module, only: rt => amrex_real - - implicit none - - ! 3-point filter: T=1, R=S=0 - real(rt), dimension(0:1), parameter :: ff1 = [0.5e0_rt, 0.25e0_rt] - - ! For boundary cell - real(rt), dimension(0:1), parameter :: ff1b = [0.75e0_rt, 0.25e0_rt] - - ! 5-point filter: T=2, R+S+1=T - ! S = 0, 1 - real(rt), dimension(0:2,0:1), parameter :: ff2 = reshape( & - [ 0.625e0_rt, 0.25e0_rt, & - -0.0625e0_rt, 0.375e0_rt, & - 0.25e0_rt, 0.0625e0_rt ], & - [3,2] ) - - ! For first boundary cell - real(rt), dimension(0:2), parameter :: ff2b0 = [17.e0_rt/16.e0_rt, & - -2.e0_rt/16.e0_rt, & - 1.e0_rt/16.e0_rt] - - ! For second boundary cell - real(rt), dimension(-1:2), parameter :: ff2b1 = [-2.e0_rt/16.e0_rt, & - 21.e0_rt/16.e0_rt, & - -4.e0_rt/16.e0_rt, & - 1.e0_rt/16.e0_rt] - - ! 7-point filter: T=3, R+S+1=T - ! S = 0, 1, 2 - real(rt), dimension(0:3,0:2), parameter :: ff3 = reshape( & - [44.e0_rt/64.e0_rt, & - 15.e0_rt/64.e0_rt, & - -6.e0_rt/64.e0_rt, & - 1.e0_rt/64.e0_rt, & - 32.e0_rt/64.e0_rt, & - 18.e0_rt/64.e0_rt, & - 0.e0_rt, -2.e0_rt/64.e0_rt, & - & 20.e0_rt/64.e0_rt, & - 15.e0_rt/64.e0_rt, & - 6.e0_rt/64.e0_rt, & - 1.e0_rt/64.e0_rt ], & - [4,3] ) - ! For boundary cells - real(rt), dimension(0:3), parameter :: ff3b0 = [63.e0_rt/64.e0_rt, & - 3.e0_rt/64.e0_rt, & - -3.e0_rt/64.e0_rt, & - 1.e0_rt/64.e0_rt] - - real(rt), dimension(-1:3), parameter :: ff3b1 = [3.e0_rt/64.e0_rt, & - 54.e0_rt/64.e0_rt, & - 12.e0_rt/64.e0_rt, & - -6.e0_rt/64.e0_rt, & - 1.e0_rt/64.e0_rt] - - real(rt), dimension(-2:3), parameter :: ff3b2 = [-3.e0_rt/64.e0_rt, & - 12.e0_rt/64.e0_rt, & - 45.e0_rt/64.e0_rt, & - 15.e0_rt/64.e0_rt, & - -6.e0_rt/64.e0_rt, & - 1.e0_rt/64.e0_rt] - - ! 9-point filter: T=4, R+S+1=T - ! S = 0, 1, 2, 3 - real(rt), dimension(0:4,0:3), parameter :: ff4 = reshape( & - [ 186.e0_rt/256.e0_rt, & - 56.e0_rt/256.e0_rt, & - -28.e0_rt/256.e0_rt, & - 8.e0_rt/256.e0_rt, & - -1.e0_rt/256.e0_rt, & - 146.e0_rt/256.e0_rt, & - 72.e0_rt/256.e0_rt, & - -12.e0_rt/256.e0_rt, & - -8.e0_rt/256.e0_rt, & - 3.e0_rt/256.e0_rt, & - 110.e0_rt/256.e0_rt, & - 72.e0_rt/256.e0_rt, & - 12.e0_rt/256.e0_rt, & - -8.e0_rt/256.e0_rt, & - -3.e0_rt/256.e0_rt, & - 70.e0_rt/256.e0_rt, & - 56.e0_rt/256.e0_rt, & - 28.e0_rt/256.e0_rt, & - 8.e0_rt/256.e0_rt, & - 1.e0_rt/256.e0_rt ], & - [5,4] ) - - ! For boundary cells - real(rt), dimension(0:4), parameter :: ff4b0 = [257.e0_rt/256.e0_rt, & - -4.e0_rt/256.e0_rt, & - 6.e0_rt/256.e0_rt, & - -4.e0_rt/256.e0_rt, & - 1.e0_rt/256.e0_rt] - - real(rt), dimension(-1:4), parameter :: ff4b1 = [-4.e0_rt/256.e0_rt, & - 273.e0_rt/256.e0_rt, & - -28.e0_rt/256.e0_rt, & - 22.e0_rt/256.e0_rt, & - -8.e0_rt/256.e0_rt, & - 1.e0_rt/256.e0_rt] - - real(rt), dimension(-2:4), parameter :: ff4b2 = [6.e0_rt/256.e0_rt, & - -28.e0_rt/256.e0_rt, & - 309.e0_rt/256.e0_rt, & - -52.e0_rt/256.e0_rt, & - 28.e0_rt/256.e0_rt, & - -8.e0_rt/256.e0_rt, & - 1.e0_rt/256.e0_rt] - - real(rt), dimension(-3:4), parameter :: ff4b3 = [-4.e0_rt/256.e0_rt, & - 22.e0_rt/256.e0_rt, & - -52.e0_rt/256.e0_rt, & - 325.e0_rt/256.e0_rt, & - -56.e0_rt/256.e0_rt, & - 28.e0_rt/256.e0_rt, & - -8.e0_rt/256.e0_rt, & - 1.e0_rt/256.e0_rt] - -end module filter_module diff --git a/Source/radiation/fluxlimiter.F90 b/Source/radiation/fluxlimiter.F90 deleted file mode 100644 index 1c29118b36..0000000000 --- a/Source/radiation/fluxlimiter.F90 +++ /dev/null @@ -1,140 +0,0 @@ -module fluxlimiter_module - - use amrex_fort_module, only: rt => amrex_real - - implicit none - - integer, allocatable, save :: limiter, closure - - real(rt), parameter :: OneThird = 1.e0_rt/3.e0_rt, TwoThirds=2.e0_rt/3.e0_rt, & - OneSixth=1.e0_rt/6.e0_rt, TwoNinths=2.e0_rt/9.e0_rt, & - FiveNinths=5.e0_rt/9.e0_rt - -contains - - subroutine init_fluxlimiter_module(limiter_in, closure_in) - - use amrex_fort_module, only: rt => amrex_real - - implicit none - - integer, intent(in) :: limiter_in, closure_in - - allocate(limiter) - allocate(closure) - - limiter = limiter_in - closure = closure_in - - end subroutine init_fluxlimiter_module - - function FLDlambda(r, limiter_in) - - use amrex_fort_module, only: rt => amrex_real - - integer, intent(in), optional :: limiter_in - real(rt), intent(in) :: r - - real(rt) :: FLDlambda - integer :: limiter_local - - if (present(limiter_in))then - limiter_local = limiter_in - else - limiter_local = limiter - end if - - if (limiter_local .eq. 0) then ! no limiter - FLDlambda = OneThird - else if (limiter_local < 10) then ! approximate LP, [123] - FLDlambda = (2.e0_rt + r) / (6.e0_rt + r * (3.e0_rt + r)) - else if (limiter_local < 20) then ! Bruenn, 1[123] - FLDlambda = 1.e0_rt / (3.e0_rt + r) - else if (limiter_local < 30) then ! Larsen's square root, 2[123] - FLDlambda = 1.e0_rt / sqrt(9.e0_rt + r**2) - else if (limiter_local < 40) then ! Minerbo, 3[123] - if (r .lt. 1.5e0_rt) then - FLDlambda = 2.e0_rt/(3.e0_rt + sqrt(9.e0_rt+12.e0_rt*r**2)) - else - FLDlambda = 1.e0_rt/(1.e0_rt+r+sqrt(1.e0_rt+2.e0_rt*r)) - end if - else - print *, "Unknown limiter ", limiter_local - stop - endif - - end function FLDlambda - - function Edd_factor(lambda) result(f) - - use amrex_fort_module, only: rt => amrex_real - - implicit none - - real(rt), intent(in) :: lambda - real(rt) :: f - - if (closure .eq. 0) then - f = lambda - else if (closure .eq. 1) then - f = OneThird - else if (closure .eq. 2) then - f = 1.e0_rt - 2.e0_rt * lambda - else if (closure .eq. 3) then ! lambda + (lambda*R)**2 - if (limiter .eq. 0) then ! no limiter - f = OneThird - else if (limiter < 10) then ! approximate LP, [123] - f = lambda + 0.25e0_rt*(max(0.0_rt, (1.e0_rt-3.e0_rt*lambda)) + & - sqrt(max(0.0_rt, (1.e0_rt-3.e0_rt*lambda))*(1.e0_rt+5.e0_rt*lambda)))**2 - else if (limiter < 20) then ! Bruenn, 1[123] - f = 1.0e0_rt - 5.e0_rt*lambda + 9.e0_rt*lambda**2 - else if (limiter < 30) then ! Larsen's square root, 2[123] - f = 1.0e0_rt + lambda - 9.e0_rt*lambda**2 - else if (limiter < 40) then ! Minerbo - if (lambda .gt. TwoNinths) then - f = OneThird - else - f = 1.e0_rt + 3.e0_rt*lambda - 2.e0_rt*sqrt(2.e0_rt*lambda) - end if - else - print *, "Unknown limiter ", limiter - stop - endif - else if (closure .eq. 4) then ! 1/3 + 2/3*(lambda*R)**2 - if (limiter .eq. 0) then ! no limiter - f = OneThird - else if (limiter < 10) then ! approximate LP, [123] - f = OneThird + OneSixth*(max(0.0_rt, (1.e0_rt-3.e0_rt*lambda)) + & - sqrt(max(0.0_rt, (1.e0_rt-3.e0_rt*lambda))*(1.e0_rt+5.e0_rt*lambda)))**2 - else if (limiter < 20) then ! Bruenn, 1[123] - f = OneThird + TwoThirds*(1.0e0_rt - 6.e0_rt*lambda + 9.e0_rt*lambda**2) - else if (limiter < 30) then ! Larsen's square root, 2[123] - f = OneThird + TwoThirds*(1.0e0_rt - 9.e0_rt*lambda**2) - else if (limiter < 40) then ! Minerbo - if (lambda .gt. TwoNinths) then - f = FiveNinths - TwoThirds*lambda - else - f = OneThird + TwoThirds*(1.e0_rt + 2.e0_rt*lambda - 2.e0_rt*sqrt(2.e0_rt*lambda)) - end if - else - print *, "Unknown limiter ", limiter - stop - endif - end if - - end function Edd_factor - -end module fluxlimiter_module - -subroutine ca_initfluxlimiter(limiter, closure) bind(C, name="ca_initfluxlimiter") - - use fluxlimiter_module, only: init_fluxlimiter_module - use amrex_fort_module, only: rt => amrex_real - - implicit none - - integer, intent(in) :: limiter, closure - - call init_fluxlimiter_module(limiter, closure) - -end subroutine ca_initfluxlimiter diff --git a/Source/sdc/sdc_util.cpp b/Source/sdc/sdc_util.cpp index 6d471ab9a3..6d2cadaf39 100644 --- a/Source/sdc/sdc_util.cpp +++ b/Source/sdc/sdc_util.cpp @@ -20,7 +20,8 @@ Castro::ca_sdc_update_advection_o2_lobatto(const Box& bx, // Gauss-Lobatto / trapezoid - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + 0.5_rt * dt * (A_0_old(i,j,k,n) + A_1_old(i,j,k,n)); }); @@ -48,7 +49,8 @@ Castro::ca_sdc_update_advection_o2_radau(const Box& bx, if (m_start == 0) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_0_old(i,j,k,n)) + @@ -57,7 +59,8 @@ Castro::ca_sdc_update_advection_o2_radau(const Box& bx, } else if (m_start == 1) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_1_old(i,j,k,n)) + @@ -86,7 +89,8 @@ Castro::ca_sdc_update_advection_o4_lobatto(const Box& bx, if (m_start == 0) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_0_old(i,j,k,n)) + @@ -95,7 +99,8 @@ Castro::ca_sdc_update_advection_o4_lobatto(const Box& bx, } else if (m_start == 1) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_1_old(i,j,k,n)) + @@ -129,7 +134,8 @@ Castro::ca_sdc_update_advection_o4_radau(const Box& bx, if (m_start == 0) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_0_old(i,j,k,n)) + @@ -140,7 +146,8 @@ Castro::ca_sdc_update_advection_o4_radau(const Box& bx, } else if (m_start == 1) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_1_old(i,j,k,n)) + @@ -151,7 +158,8 @@ Castro::ca_sdc_update_advection_o4_radau(const Box& bx, } else if (m_start == 2) { - AMREX_PARALLEL_FOR_4D(bx, k_n.nComp(), i, j, k, n, + amrex::ParallelFor(bx, k_n.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { k_n(i,j,k,n) = k_m(i,j,k,n) + dt_m * (A_m(i,j,k,n) - A_2_old(i,j,k,n)) + @@ -182,7 +190,8 @@ Castro::ca_sdc_compute_C2_lobatto(const Box& bx, // Here, dt_m is the timestep between time-nodes m and m+1 - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // construct the source term to the update for 2nd order // Lobatto, there is no advective correction, and we have @@ -216,7 +225,8 @@ Castro::ca_sdc_compute_C2_radau(const Box& bx, if (m_start == 0) { - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { C(i,j,k,n) = -R_1_old(i,j,k,n) + (A_m(i,j,k,n) - A_0_old(i,j,k,n)) + @@ -227,7 +237,8 @@ Castro::ca_sdc_compute_C2_radau(const Box& bx, } else if (m_start == 1) { - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { C(i,j,k,n) = -R_2_old(i,j,k,n) + (A_m(i,j,k,n) - A_1_old(i,j,k,n)) + @@ -260,7 +271,8 @@ Castro::ca_sdc_compute_C4_lobatto(const Box& bx, if (m_start == 0) { // compute the integral from [t_m, t_{m+1}], normalized by dt_m - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real integral = 1.0_rt/12.0_rt * (5.0_rt*(A_0_old(i,j,k,n) + R_0_old(i,j,k,n)) + 8.0_rt*(A_1_old(i,j,k,n) + R_1_old(i,j,k,n)) - @@ -272,7 +284,8 @@ Castro::ca_sdc_compute_C4_lobatto(const Box& bx, else if (m_start == 1) { // compute the integral from [t_m, t_{m+1}], normalized by dt_m - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real integral = 1.0_rt/12.0_rt * (-(A_0_old(i,j,k,n) + R_0_old(i,j,k,n)) + 8.0_rt*(A_1_old(i,j,k,n) + R_1_old(i,j,k,n)) + @@ -312,7 +325,8 @@ Castro::ca_sdc_compute_C4_radau(const Box& bx, if (m_start == 0) { // compute the integral from [t_m, t_{m+1}], normalized by dt_m - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real integral = (dt/dt_m) * (1.0_rt/1800.0_rt) * ((-35.0_rt*std::sqrt(6.0_rt) + 440.0_rt)*(A_1_old(i,j,k,n) + R_1_old(i,j,k,n)) + @@ -325,7 +339,8 @@ Castro::ca_sdc_compute_C4_radau(const Box& bx, else if (m_start == 1) { // compute the integral from [t_m, t_{m+1}], normalized by dt_m - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real integral = (dt/dt_m) * (1.0_rt/150.0_rt) * ((-12.0_rt + 17.0_rt*std::sqrt(6.0_rt))*(A_1_old(i,j,k,n) + R_1_old(i,j,k,n)) + @@ -338,7 +353,8 @@ Castro::ca_sdc_compute_C4_radau(const Box& bx, else if (m_start == 2) { // compute the integral from [t_m, t_{m+1}], normalized by dt_m - AMREX_PARALLEL_FOR_4D(bx, C.nComp(), i, j, k, n, + amrex::ParallelFor(bx, C.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real integral = (dt/dt_m) * (1.0_rt/600.0_rt) * ((168.0_rt - 73.0_rt*std::sqrt(6.0_rt))*(A_1_old(i,j,k,n) + R_1_old(i,j,k,n)) + @@ -364,7 +380,8 @@ Castro::ca_sdc_conservative_update(const Box& bx, Real const dt_m, // given _old, _new, and , compute _new // now consider the reacting system - AMREX_PARALLEL_FOR_4D(bx, U_new.nComp(), i, j, k, n, + amrex::ParallelFor(bx, U_new.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { U_new(i,j,k,n) = U_old(i,j,k,n) + dt_m * R_new(i,j,k,n) + dt_m * C(i,j,k,n); }); @@ -385,14 +402,16 @@ void Castro::ca_sdc_compute_initial_guess(const Box& bx, if (sdc_iteration == 0) { - AMREX_PARALLEL_FOR_4D(bx, U_guess.nComp(), i, j, k, n, + amrex::ParallelFor(bx, U_guess.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { U_guess(i,j,k,n) = U_old(i,j,k,n) + dt_m * A_old(i,j,k,n) + dt_m * R_old(i,j,k,n); }); } else { - AMREX_PARALLEL_FOR_4D(bx, U_guess.nComp(), i, j, k, n, + amrex::ParallelFor(bx, U_guess.nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { U_guess(i,j,k,n) = U_new(i,j,k,n); }); @@ -412,20 +431,23 @@ void Castro::ca_store_reaction_state(const Box& bx, // Reactions_Type if (store_omegadot) { - AMREX_PARALLEL_FOR_4D(bx, NumSpec, i, j, k, n, + amrex::ParallelFor(bx, NumSpec, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { R_store(i,j,k,1+n) = R_old(i,j,k,UFS+n); }); #if NAUX_NET > 0 - AMREX_PARALLEL_FOR_4D(bx, NumAux, i, j, k, n, + amrex::ParallelFor(bx, NumAux, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { R_store(i,j,k,1+NumSpec+n) = R_old(i,j,k,UFX+n); }); #endif } - AMREX_PARALLEL_FOR_3D(bx, i, j, k, + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { R_store(i,j,k,0) = R_old(i,j,k,UEDEN); });