diff --git a/Source/radiation/CastroRad_1d.f90 b/Source/radiation/CastroRad_1d.f90 index ed6293e7e8..94732d9c25 100644 --- a/Source/radiation/CastroRad_1d.f90 +++ b/Source/radiation/CastroRad_1d.f90 @@ -1,223 +1,3 @@ -! This subroutine cannot be tiled -subroutine ca_compute_lamborder(Er, Er_l1, Er_h1, & - kap, kap_l1, kap_h1, & - lam, lam_l1, lam_h1, & - dx, ngrow, limiter, filter_T, S) - - use rad_params_module, only : ngroups - use fluxlimiter_module, only : FLDlambda - use filter_module - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: Er_l1, Er_h1, kap_l1, kap_h1, lam_l1, lam_h1 - integer, intent(in) :: ngrow, limiter, filter_T, S - real(rt) , intent(in) :: dx - real(rt) , intent(in) :: kap(kap_l1:kap_h1, 0:ngroups-1) - real(rt) , intent(in) :: Er(Er_l1:Er_h1, 0:ngroups-1) - real(rt) , intent(out) :: lam(lam_l1:lam_h1, 0:ngroups-1) - - integer :: i, reg_l1, reg_h1, g - real(rt) :: r - - real(rt) , allocatable :: lamfil(:) - - lam = -1.e50_rt - - reg_l1 = lam_l1 + ngrow - reg_h1 = lam_h1 - ngrow - - if (filter_T .gt. 0) then - allocate(lamfil(reg_l1:reg_h1)) - end if - - do g = 0, ngroups-1 - do i=lam_l1, lam_h1 - r = abs(Er(i+1,g) - Er(i-1,g)) / (2.*dx) - r = r / (kap(i,g) * max(Er(i,g), 1.e-50_rt)) - lam(i,g) = FLDlambda(r, limiter) - end do - - if (Er(reg_l1-1,g) .eq. -1.e0_rt) then - r = abs(Er(reg_l1+1,g) - Er(reg_l1,g)) / dx - r = r / (kap(reg_l1,g) * max(Er(reg_l1,g), 1.e-50_rt)) - lam(reg_l1,g) = FLDlambda(r, limiter) - end if - - if (Er(reg_h1+1,g) .eq. -1.e0_rt) then - r = abs(Er(reg_h1,g) - Er(reg_h1-1,g)) / dx - r = r / (kap(reg_h1,g) * max(Er(reg_h1,g), 1.e-50_rt)) - lam(reg_h1,g) = FLDlambda(r, limiter) - end if - - ! filter - if (filter_T .eq. 1) then - - do i=reg_l1, reg_h1 - lamfil(i) = ff1(0) * lam(i,g) & - & + ff1(1) * (lam(i-1,g)+lam(i+1,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end do - - if (Er(reg_l1-1,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i) = dot_product(ff1b, lam(i:i+1,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - if (Er(reg_h1+1,g) .eq. -1.e0_rt) then - i = reg_h1 - lamfil(i) = dot_product(ff1b(1:0:-1), lam(i-1:i,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - lam(reg_l1:reg_h1,g) = lamfil(reg_l1:reg_h1) - - else if (filter_T .eq. 2) then - - do i=reg_l1, reg_h1 - lamfil(i) = ff2(0,S) * lam(i,g) & - & + ff2(1,S) * (lam(i-1,g)+lam(i+1,g)) & - & + ff2(2,S) * (lam(i-2,g)+lam(i+2,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end do - - if (Er(reg_l1-1,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i) = dot_product(ff2b0, lam(i:i+2,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_l1 + 1 - lamfil(i) = dot_product(ff2b1, lam(i-1:i+2,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - if (Er(reg_h1+1,g) .eq. -1.e0_rt) then - i = reg_h1-1 - lamfil(i) = dot_product(ff2b1(2:-1:-1), lam(i-2:i+1,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_h1 - lamfil(i) = dot_product(ff2b0(2:0:-1), lam(i-2:i,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - lam(reg_l1:reg_h1,g) = lamfil(reg_l1:reg_h1) - - else if (filter_T .eq. 3) then - - do i=reg_l1, reg_h1 - lamfil(i) = ff3(0,S) * lam(i,g) & - & + ff3(1,S) * (lam(i-1,g)+lam(i+1,g)) & - & + ff3(2,S) * (lam(i-2,g)+lam(i+2,g)) & - & + ff3(3,S) * (lam(i-3,g)+lam(i+3,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end do - - if (Er(reg_l1-1,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i) = dot_product(ff3b0, lam(i:i+3,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_l1 + 1 - lamfil(i) = dot_product(ff3b1, lam(i-1:i+3,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_l1 + 2 - lamfil(i) = dot_product(ff3b2, lam(i-2:i+3,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - if (Er(reg_h1+1,g) .eq. -1.e0_rt) then - i = reg_h1 - 2 - lamfil(i) = dot_product(ff3b2(3:-2:-1), lam(i-3:i+2,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_h1 - 1 - lamfil(i) = dot_product(ff3b1(3:-1:-1), lam(i-3:i+1,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_h1 - lamfil(i) = dot_product(ff3b0(3:0:-1), lam(i-3:i,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - lam(reg_l1:reg_h1,g) = lamfil(reg_l1:reg_h1) - - else if (filter_T .eq. 4) then - - do i=reg_l1, reg_h1 - lamfil(i) = ff4(0,S) * lam(i,g) & - & + ff4(1,S) * (lam(i-1,g)+lam(i+1,g)) & - & + ff4(2,S) * (lam(i-2,g)+lam(i+2,g)) & - & + ff4(3,S) * (lam(i-3,g)+lam(i+3,g)) & - & + ff4(4,S) * (lam(i-4,g)+lam(i+4,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end do - - if (Er(reg_l1-1,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i) = dot_product(ff4b0, lam(i:i+4,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_l1 + 1 - lamfil(i) = dot_product(ff4b1, lam(i-1:i+4,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_l1 + 2 - lamfil(i) = dot_product(ff4b2, lam(i-2:i+4,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_l1 + 3 - lamfil(i) = dot_product(ff4b3, lam(i-3:i+4,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - if (Er(reg_h1+1,g) .eq. -1.e0_rt) then - i = reg_h1 - 3 - lamfil(i) = dot_product(ff4b3(4:-3:-1), lam(i-4:i+3,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_h1 - 2 - lamfil(i) = dot_product(ff4b2(4:-2:-1), lam(i-4:i+2,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_h1 - 1 - lamfil(i) = dot_product(ff4b1(4:-1:-1), lam(i-4:i+1,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - - i = reg_h1 - lamfil(i) = dot_product(ff4b0(4:0:-1), lam(i-4:i,g)) - lamfil(i) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i))) - end if - - lam(reg_l1:reg_h1,g) = lamfil(reg_l1:reg_h1) - - end if - - ! boundary - - if (Er(reg_l1-1,g) .eq. -1.e0_rt) then - do i=lam_l1, reg_l1-1 - lam(i,g) = lam(reg_l1,g) - end do - end if - - if (Er(reg_h1+1,g) .eq. -1.e0_rt) then - do i=reg_h1+1, lam_h1 - lam(i,g) = lam(reg_h1,g) - end do - end if - - end do - - if (filter_T .gt. 0) then - deallocate(lamfil) - end if - -end subroutine ca_compute_lamborder - ! no tiling subroutine ca_correct_dterm(dfx, dfx_l1, dfx_h1, & re, rc) bind(C, name="ca_correct_dterm") diff --git a/Source/radiation/CastroRad_2d.f90 b/Source/radiation/CastroRad_2d.f90 index 6fa40f49be..892aadac3b 100644 --- a/Source/radiation/CastroRad_2d.f90 +++ b/Source/radiation/CastroRad_2d.f90 @@ -1,467 +1,3 @@ -! This subroutine cannot be tiled -subroutine ca_compute_lamborder(Er, Er_l1, Er_l2, Er_h1, Er_h2, & - kap, kap_l1, kap_l2, kap_h1, kap_h2, & - lam, lam_l1, lam_l2, lam_h1, lam_h2, & - dx, ngrow, limiter, filter_T, S) - - use rad_params_module, only : ngroups - use fluxlimiter_module, only : FLDlambda - use filter_module - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: Er_l1, Er_l2, Er_h1, Er_h2, kap_l1, kap_l2, kap_h1, kap_h2, & - lam_l1, lam_l2, lam_h1, lam_h2 - integer, intent(in) :: ngrow, limiter, filter_T, S - real(rt) , intent(in) :: dx(2) - real(rt) , intent(in) :: kap(kap_l1:kap_h1, kap_l2:kap_h2) - real(rt) , intent(in) :: Er(Er_l1:Er_h1, Er_l2:Er_h2, 0:ngroups-1) - real(rt) , intent(out) :: lam(lam_l1:lam_h1, lam_l2:lam_h2, 0:ngroups-1) - - integer :: i, j, reg_l1, reg_l2, reg_h1, reg_h2, g - real(rt) :: r, r1, r2 - - real(rt) , allocatable :: lamfil(:,:) - - lam = -1.e50_rt - - reg_l1 = lam_l1 + ngrow - reg_l2 = lam_l2 + ngrow - reg_h1 = lam_h1 - ngrow - reg_h2 = lam_h2 - ngrow - - if (filter_T .gt. 0) then - allocate(lamfil(reg_l1:reg_h1,lam_l2:lam_h2)) - end if - - do g = 0, ngroups-1 - - do j=lam_l2, lam_h2 - do i=lam_l1, lam_h1 - if (Er(i,j,g) .eq. -1.e0_rt) then - cycle - end if - - if (Er(i-1,j,g) .eq. -1.e0_rt) then - r1 = (Er(i+1,j,g) - Er(i,j,g)) / (dx(1)) - else if (Er(i+1,j,g) .eq. -1.e0_rt) then - r1 = (Er(i,j,g) - Er(i-1,j,g)) / (dx(1)) - else - r1 = (Er(i+1,j,g) - Er(i-1,j,g)) / (2.e0_rt*dx(1)) - end if - - if (Er(i,j-1,g) .eq. -1.e0_rt) then - r2 = (Er(i,j+1,g) - Er(i,j,g)) / (dx(2)) - else if (Er(i,j+1,g) .eq. -1.e0_rt) then - r2 = (Er(i,j,g) - Er(i,j-1,g)) / (dx(2)) - else - r2 = (Er(i,j+1,g) - Er(i,j-1,g)) / (2.e0_rt*dx(2)) - end if - - r = sqrt(r1**2 + r2**2) - r = r / (kap(i,j) * max(Er(i,j,g), 1.e-50_rt)) - - lam(i,j,g) = FLDlambda(r, limiter) - end do - end do - - ! filter - if (filter_T .eq. 1) then - - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,g) .eq. -1.e0_rt) then - lamfil(:,j) = -1.e-50_rt - cycle - endif - - do i=reg_l1, reg_h1 - lamfil(i,j) = ff1(0) * lam(i,j,g) & - & + ff1(1) * (lam(i-1,j,g)+lam(i+1,j,g)) - end do - - if (Er(reg_l1-1,j,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j) = dot_product(ff1b, lam(i:i+1,j,g)) - end if - - if (Er(reg_h1+1,j,g) .eq. -1.e0_rt) then - i = reg_h1 - lamfil(i,j) = dot_product(ff1b(1:0:-1), lam(i-1:i,j,g)) - end if - end do - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,g) = ff1(0) * lamfil(i,j) & - & + ff1(1) * (lamfil(i,j-1)+lamfil(i,j+1)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end do - - if (Er(reg_l1,reg_l2-1,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff1b, lamfil(i,j:j+1)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - if (Er(reg_l1,reg_h2+1,g) .eq. -1.e0_rt) then - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff1b(1:0:-1), lamfil(i,j-1:j)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - else if (filter_T .eq. 2) then - - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,g) .eq. -1.e0_rt) then - lamfil(:,j) = -1.e-50_rt - cycle - endif - - do i=reg_l1, reg_h1 - lamfil(i,j) = ff2(0,S) * lam(i,j,g) & - & + ff2(1,S) * (lam(i-1,j,g)+lam(i+1,j,g)) & - & + ff2(2,S) * (lam(i-2,j,g)+lam(i+2,j,g)) - end do - - if (Er(reg_l1-1,j,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j) = dot_product(ff2b0, lam(i:i+2,j,g)) - - i = reg_l1 + 1 - lamfil(i,j) = dot_product(ff2b1, lam(i-1:i+2,j,g)) - end if - - if (Er(reg_h1+1,j,g) .eq. -1.e0_rt) then - i = reg_h1 - 1 - lamfil(i,j) = dot_product(ff2b1(2:-1:-1), lam(i-2:i+1,j,g)) - - i = reg_h1 - lamfil(i,j) = dot_product(ff2b0(2:0:-1), lam(i-2:i,j,g)) - end if - end do - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,g) = ff2(0,S) * lamfil(i,j) & - & + ff2(1,S) * (lamfil(i,j-1)+lamfil(i,j+1)) & - & + ff2(2,S) * (lamfil(i,j-2)+lamfil(i,j+2)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end do - - if (Er(reg_l1,reg_l2-1,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff2b0, lamfil(i,j:j+2)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_l2 + 1 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff2b1, lamfil(i,j-1:j+2)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - if (Er(reg_l1,reg_h2+1,g) .eq. -1.e0_rt) then - j = reg_h2 - 1 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff2b1(2:-1:-1), lamfil(i,j-2:j+1)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff2b0(2:0:-1), lamfil(i,j-2:j)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - else if (filter_T .eq. 3) then - - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,g) .eq. -1.e0_rt) then - lamfil(:,j) = -1.e-50_rt - cycle - endif - - do i=reg_l1, reg_h1 - lamfil(i,j) = ff3(0,S) * lam(i,j,g) & - & + ff3(1,S) * (lam(i-1,j,g)+lam(i+1,j,g)) & - & + ff3(2,S) * (lam(i-2,j,g)+lam(i+2,j,g)) & - & + ff3(3,S) * (lam(i-3,j,g)+lam(i+3,j,g)) - end do - - if (Er(reg_l1-1,j,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j) = dot_product(ff3b0, lam(i:i+3,j,g)) - - i = reg_l1 + 1 - lamfil(i,j) = dot_product(ff3b1, lam(i-1:i+3,j,g)) - - i = reg_l1 + 2 - lamfil(i,j) = dot_product(ff3b2, lam(i-2:i+3,j,g)) - end if - - if (Er(reg_h1+1,j,g) .eq. -1.e0_rt) then - i = reg_h1 - 2 - lamfil(i,j) = dot_product(ff3b2(3:-2:-1), lam(i-3:i+2,j,g)) - - i = reg_h1 - 1 - lamfil(i,j) = dot_product(ff3b1(3:-1:-1), lam(i-3:i+1,j,g)) - - i = reg_h1 - lamfil(i,j) = dot_product(ff3b0(3:0:-1), lam(i-3:i,j,g)) - end if - end do - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,g) = ff3(0,S) * lamfil(i,j) & - & + ff3(1,S) * (lamfil(i,j-1)+lamfil(i,j+1)) & - & + ff3(2,S) * (lamfil(i,j-2)+lamfil(i,j+2)) & - & + ff3(3,S) * (lamfil(i,j-3)+lamfil(i,j+3)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end do - - if (Er(reg_l1,reg_l2-1,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff3b0, lamfil(i,j:j+3)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_l2 + 1 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff3b1, lamfil(i,j-1:j+3)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_l2 + 2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff3b2, lamfil(i,j-2:j+3)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - if (Er(reg_l1,reg_h2+1,g) .eq. -1.e0_rt) then - j = reg_h2 - 2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff3b2(3:-2:-1), lamfil(i,j-3:j+2)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_h2 - 1 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff3b1(3:-1:-1), lamfil(i,j-3:j+1)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff3b0(3:0:-1), lamfil(i,j-3:j)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - else if (filter_T .eq. 4) then - - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,g) .eq. -1.e0_rt) then - lamfil(:,j) = -1.e-50_rt - cycle - endif - - do i=reg_l1, reg_h1 - lamfil(i,j) = ff4(0,S) * lam(i,j,g) & - & + ff4(1,S) * (lam(i-1,j,g)+lam(i+1,j,g)) & - & + ff4(2,S) * (lam(i-2,j,g)+lam(i+2,j,g)) & - & + ff4(3,S) * (lam(i-3,j,g)+lam(i+3,j,g)) & - & + ff4(4,S) * (lam(i-4,j,g)+lam(i+4,j,g)) - end do - - if (Er(reg_l1-1,j,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j) = dot_product(ff4b0, lam(i:i+4,j,g)) - - i = reg_l1 + 1 - lamfil(i,j) = dot_product(ff4b1, lam(i-1:i+4,j,g)) - - i = reg_l1 + 2 - lamfil(i,j) = dot_product(ff4b2, lam(i-2:i+4,j,g)) - - i = reg_l1 + 3 - lamfil(i,j) = dot_product(ff4b3, lam(i-3:i+4,j,g)) - end if - - if (Er(reg_h1+1,j,g) .eq. -1.e0_rt) then - i = reg_h1 - 3 - lamfil(i,j) = dot_product(ff4b3(4:-3:-1), lam(i-4:i+3,j,g)) - - i = reg_h1 - 2 - lamfil(i,j) = dot_product(ff4b2(4:-2:-1), lam(i-4:i+2,j,g)) - - i = reg_h1 - 1 - lamfil(i,j) = dot_product(ff4b1(4:-1:-1), lam(i-4:i+1,j,g)) - - i = reg_h1 - lamfil(i,j) = dot_product(ff4b0(4:0:-1), lam(i-4:i,j,g)) - end if - end do - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,g) = ff4(0,S) * lamfil(i,j) & - & + ff4(1,S) * (lamfil(i,j-1)+lamfil(i,j+1)) & - & + ff4(2,S) * (lamfil(i,j-2)+lamfil(i,j+2)) & - & + ff4(3,S) * (lamfil(i,j-3)+lamfil(i,j+3)) & - & + ff4(4,S) * (lamfil(i,j-4)+lamfil(i,j+4)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end do - - if (Er(reg_l1,reg_l2-1,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b0, lamfil(i,j:j+4)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_l2 + 1 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b1, lamfil(i,j-1:j+4)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_l2 + 2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b2, lamfil(i,j-2:j+4)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_l2 + 3 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b3, lamfil(i,j-3:j+4)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - if (Er(reg_l1,reg_h2+1,g) .eq. -1.e0_rt) then - j = reg_h2 - 3 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b3(4:-3:-1), lamfil(i,j-4:j+3)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_h2 - 2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b2(4:-2:-1), lamfil(i,j-4:j+2)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_h2 - 1 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b1(4:-1:-1), lamfil(i,j-4:j+1)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,g) = dot_product(ff4b0(4:0:-1), lamfil(i,j-4:j)) - lam(i,j,g) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lam(i,j,g))) - end do - end if - - end if - - ! lo-x lo-y - do j=lam_l2,reg_l2-1 - do i=lam_l1,reg_l1-1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(reg_l1,reg_l2,g) - end if - end do - end do - - ! reg-x lo-y - do j=lam_l2,reg_l2-1 - do i=reg_l1,reg_h1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(i,reg_l2,g) - end if - end do - end do - - ! hi-x lo-y - do j=lam_l2,reg_l2-1 - do i=reg_h1+1,lam_h1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(reg_h1,reg_l2,g) - end if - end do - end do - - ! lo-x reg-y - do j=reg_l2,reg_h2 - do i=lam_l1,reg_l1-1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(reg_l1,j,g) - end if - end do - end do - - ! hi-x reg-y - do j=reg_l2,reg_h2 - do i=reg_h1+1,lam_h1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(reg_h1,j,g) - end if - end do - end do - - ! lo-x hi-y - do j=reg_h2+1,lam_h2 - do i=lam_l1,reg_l1-1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(reg_l1,reg_h2,g) - end if - end do - end do - - ! reg-x hi-y - do j=reg_h2+1,lam_h2 - do i=reg_l1,reg_h1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(i,reg_h2,g) - end if - end do - end do - - ! hi-x hi-y - do j=reg_h2+1,lam_h2 - do i=reg_h1+1,lam_h1 - if (Er(i,j,g).eq.-1.e0_rt) then - lam(i,j,g) = lam(reg_h1,reg_h2,g) - end if - end do - end do - - end do - - if (filter_T .gt. 0) then - deallocate(lamfil) - end if - - return -end subroutine ca_compute_lamborder - ! no tiling subroutine ca_correct_dterm( & dfx, dfx_l1, dfx_l2, dfx_h1, dfx_h2, & diff --git a/Source/radiation/CastroRad_3d.f90 b/Source/radiation/CastroRad_3d.f90 index 4685964646..c21d2945ac 100644 --- a/Source/radiation/CastroRad_3d.f90 +++ b/Source/radiation/CastroRad_3d.f90 @@ -1,928 +1,3 @@ -! This subroutine cannot be tiled -subroutine ca_compute_lamborder(Er, Er_l1, Er_l2, Er_l3, Er_h1, Er_h2, Er_h3, & - kap, kap_l1, kap_l2, kap_l3, kap_h1, kap_h2, kap_h3, & - lam, lam_l1, lam_l2, lam_l3, lam_h1, lam_h2, lam_h3, & - dx, ngrow, limiter, filter_T, S) - - use rad_params_module, only : ngroups - use fluxlimiter_module, only : FLDlambda - use filter_module - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: Er_l1, Er_l2, Er_l3, Er_h1, Er_h2, Er_h3, & - kap_l1, kap_l2, kap_l3, kap_h1, kap_h2, kap_h3, & - lam_l1, lam_l2, lam_l3, lam_h1, lam_h2, lam_h3 - integer, intent(in) :: ngrow, limiter, filter_T, S - real(rt) , intent(in) :: dx(3) - real(rt) , intent(in) :: kap(kap_l1:kap_h1, kap_l2:kap_h2, kap_l3:kap_h3, 0:ngroups-1) - real(rt) , intent(in) :: Er(Er_l1:Er_h1, Er_l2:Er_h2, Er_l3:Er_h3, 0:ngroups-1) - real(rt) , intent(out) :: lam(lam_l1:lam_h1, lam_l2:lam_h2, lam_l3:lam_h3, 0:ngroups-1) - - integer :: i, j, k, reg_l1, reg_l2, reg_l3, reg_h1, reg_h2, reg_h3, g - real(rt) :: r, r1, r2, r3 - - real(rt) , allocatable :: lamfil(:,:,:) - - reg_l1 = lam_l1 + ngrow - reg_l2 = lam_l2 + ngrow - reg_l3 = lam_l3 + ngrow - reg_h1 = lam_h1 - ngrow - reg_h2 = lam_h2 - ngrow - reg_h3 = lam_h3 - ngrow - - if (filter_T .gt. 0) then - allocate(lamfil(reg_l1:reg_h1,lam_l2:lam_h2,lam_l3:lam_h3)) - end if - - do g = 0, ngroups-1 - - do k=lam_l3, lam_h3 - do j=lam_l2, lam_h2 - do i=lam_l1, lam_h1 - - lam(i,j,k,g) = -1.e50_rt - - if (Er(i,j,k,g) .eq. -1.e0_rt) then - cycle - end if - - if (Er(i-1,j,k,g) .eq. -1.e0_rt) then - r1 = (Er(i+1,j,k,g) - Er(i,j,k,g)) / (dx(1)) - else if (Er(i+1,j,k,g) .eq. -1.e0_rt) then - r1 = (Er(i,j,k,g) - Er(i-1,j,k,g)) / (dx(1)) - else - r1 = (Er(i+1,j,k,g) - Er(i-1,j,k,g)) / (2.e0_rt*dx(1)) - end if - - if (Er(i,j-1,k,g) .eq. -1.e0_rt) then - r2 = (Er(i,j+1,k,g) - Er(i,j,k,g)) / (dx(2)) - else if (Er(i,j+1,k,g) .eq. -1.e0_rt) then - r2 = (Er(i,j,k,g) - Er(i,j-1,k,g)) / (dx(2)) - else - r2 = (Er(i,j+1,k,g) - Er(i,j-1,k,g)) / (2.e0_rt*dx(2)) - end if - - if (Er(i,j,k-1,g) .eq. -1.e0_rt) then - r3 = (Er(i,j,k+1,g) - Er(i,j,k,g)) / dx(3) - else if (Er(i,j,k+1,g) .eq. -1.e0_rt) then - r3 = (Er(i,j,k,g) - Er(i,j,k-1,g)) / dx(3) - else - r3 = (Er(i,j,k+1,g) - Er(i,j,k-1,g)) / (2.e0_rt*dx(3)) - end if - - r = sqrt(r1**2 + r2**2 + r3**2) - r = r / (kap(i,j,k,g) * max(Er(i,j,k,g), 1.e-50_rt)) - - lam(i,j,k,g) = FLDlambda(r, limiter) - end do - end do - end do - - ! filter - if (filter_T .eq. 1) then - do k=lam_l3, lam_h3 - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,k,g) .eq. -1.e0_rt) then - lamfil(:,j,k) = -1.e-50_rt - cycle - end if - - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff1(0) * lam(i,j,k,g) & - & + ff1(1) * (lam(i-1,j,k,g)+lam(i+1,j,k,g)) - end do - - if (Er(reg_l1-1,j,k,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j,k) = dot_product(ff1b, lam(i:i+1,j,k,g)) - end if - - if (Er(reg_h1+1,j,k,g) .eq. -1.e0_rt) then - i = reg_h1 - lamfil(i,j,k) = dot_product(ff1b(1:0:-1), lam(i-1:i,j,k,g)) - end if - end do - end do - - do k=lam_l3, lam_h3 - if (Er(reg_l1,reg_l2,k,g) .eq. -1.e0_rt) then - cycle - end if - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,k,g) = ff1(0) * lamfil(i,j,k) & - & + ff1(1) * (lamfil(i,j-1,k)+lamfil(i,j+1,k)) - end do - end do - - if (Er(reg_l1,reg_l2-1,k,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff1b, lamfil(i,j:j+1,k)) - end do - end if - - if (Er(reg_l1,reg_h2+1,k,g) .eq. -1.e0_rt) then - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff1b(1:0:-1), lamfil(i,j-1:j,k)) - end do - end if - end do - - do k=reg_l3, reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff1(0) * lam(i,j,k,g) & - & + ff1(1) * (lam(i,j,k-1,g)+lam(i,j,k+1,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end do - - if (Er(reg_l1,reg_l2,reg_l3-1,g) .eq. -1.e0_rt) then - k = reg_l3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff1b, lam(i,j,k:k+1,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - if (Er(reg_l1,reg_l2,reg_h3+1,g) .eq. -1.e0_rt) then - k = reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff1b(1:0:-1), lam(i,j,k-1:k,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - lam(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3,g) = & - lamfil(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3) - - else if (filter_T .eq. 2) then - - do k=lam_l3, lam_h3 - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,k,g) .eq. -1.e0_rt) then - lamfil(:,j,k) = -1.e-50_rt - cycle - end if - - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff2(0,S) * lam(i,j,k,g) & - & + ff2(1,S) * (lam(i-1,j,k,g)+lam(i+1,j,k,g)) & - & + ff2(2,S) * (lam(i-2,j,k,g)+lam(i+2,j,k,g)) - end do - - if (Er(reg_l1-1,j,k,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j,k) = dot_product(ff2b0, lam(i:i+2,j,k,g)) - - i = reg_l1 + 1 - lamfil(i,j,k) = dot_product(ff2b1, lam(i-1:i+2,j,k,g)) - end if - - if (Er(reg_h1+1,j,k,g) .eq. -1.e0_rt) then - i = reg_h1 - 1 - lamfil(i,j,k) = dot_product(ff2b1(2:-1:-1), lam(i-2:i+1,j,k,g)) - - i = reg_h1 - lamfil(i,j,k) = dot_product(ff2b0(2:0:-1), lam(i-2:i,j,k,g)) - end if - end do - end do - - do k=lam_l3, lam_h3 - if (Er(reg_l1,reg_l2,k,g) .eq. -1.e0_rt) then - cycle - end if - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,k,g) = ff2(0,S) * lamfil(i,j,k) & - & + ff2(1,S) * (lamfil(i,j-1,k)+lamfil(i,j+1,k)) & - & + ff2(2,S) * (lamfil(i,j-2,k)+lamfil(i,j+2,k)) - end do - end do - - if (Er(reg_l1,reg_l2-1,k,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff2b0, lamfil(i,j:j+2,k)) - end do - - j = reg_l2 + 1 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff2b1, lamfil(i,j-1:j+2,k)) - end do - end if - - if (Er(reg_l1,reg_h2+1,k,g) .eq. -1.e0_rt) then - j = reg_h2 - 1 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff2b1(2:-1:-1), lamfil(i,j-2:j+1,k)) - end do - - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff2b0(2:0:-1), lamfil(i,j-2:j,k)) - end do - end if - end do - - do k=reg_l3, reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff2(0,S) * lam(i,j,k,g) & - & + ff2(1,S) * (lam(i,j,k-1,g)+lam(i,j,k+1,g)) & - & + ff2(2,S) * (lam(i,j,k-2,g)+lam(i,j,k+2,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end do - - if (Er(reg_l1,reg_l2,reg_l3-1,g) .eq. -1.e0_rt) then - k = reg_l3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff2b0, lam(i,j,k:k+2,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_l3 + 1 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff2b1, lam(i,j,k-1:k+2,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - if (Er(reg_l1,reg_l2,reg_h3+1,g) .eq. -1.e0_rt) then - k = reg_h3 - 1 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff2b1(2:-1:-1), lam(i,j,k-2:k+1,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff2b0(2:0:-1), lam(i,j,k-2:k,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - lam(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3,g) = & - lamfil(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3) - - else if (filter_T .eq. 3) then - - do k=lam_l3, lam_h3 - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,k,g) .eq. -1.e0_rt) then - lamfil(:,j,k) = -1.e-50_rt - cycle - end if - - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff3(0,S) * lam(i,j,k,g) & - & + ff3(1,S) * (lam(i-1,j,k,g)+lam(i+1,j,k,g)) & - & + ff3(2,S) * (lam(i-2,j,k,g)+lam(i+2,j,k,g)) & - & + ff3(3,S) * (lam(i-3,j,k,g)+lam(i+3,j,k,g)) - end do - - if (Er(reg_l1-1,j,k,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j,k) = dot_product(ff3b0, lam(i:i+3,j,k,g)) - - i = reg_l1 + 1 - lamfil(i,j,k) = dot_product(ff3b1, lam(i-1:i+3,j,k,g)) - - i = reg_l1 + 2 - lamfil(i,j,k) = dot_product(ff3b2, lam(i-2:i+3,j,k,g)) - end if - - if (Er(reg_h1+1,j,k,g) .eq. -1.e0_rt) then - i = reg_h1 - 2 - lamfil(i,j,k) = dot_product(ff3b2(3:-2:-1), lam(i-3:i+2,j,k,g)) - - i = reg_h1 - 1 - lamfil(i,j,k) = dot_product(ff3b1(3:-1:-1), lam(i-3:i+1,j,k,g)) - - i = reg_h1 - lamfil(i,j,k) = dot_product(ff3b0(3:0:-1), lam(i-3:i,j,k,g)) - end if - end do - end do - - do k=lam_l3, lam_h3 - if (Er(reg_l1,reg_l2,k,g) .eq. -1.e0_rt) then - cycle - end if - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,k,g) = ff3(0,S) * lamfil(i,j,k) & - & + ff3(1,S) * (lamfil(i,j-1,k)+lamfil(i,j+1,k)) & - & + ff3(2,S) * (lamfil(i,j-2,k)+lamfil(i,j+2,k)) & - & + ff3(3,S) * (lamfil(i,j-3,k)+lamfil(i,j+3,k)) - end do - end do - - if (Er(reg_l1,reg_l2-1,k,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff3b0, lamfil(i,j:j+3,k)) - end do - - j = reg_l2 + 1 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff3b1, lamfil(i,j-1:j+3,k)) - end do - - j = reg_l2 + 2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff3b2, lamfil(i,j-2:j+3,k)) - end do - end if - - if (Er(reg_l1,reg_h2+1,k,g) .eq. -1.e0_rt) then - j = reg_h2 - 2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff3b2(3:-2:-1), lamfil(i,j-3:j+2,k)) - end do - - j = reg_h2 - 1 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff3b1(3:-1:-1), lamfil(i,j-3:j+1,k)) - end do - - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff3b0(3:0:-1), lamfil(i,j-3:j,k)) - end do - end if - end do - - do k=reg_l3, reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff3(0,S) * lam(i,j,k,g) & - & + ff3(1,S) * (lam(i,j,k-1,g)+lam(i,j,k+1,g)) & - & + ff3(2,S) * (lam(i,j,k-2,g)+lam(i,j,k+2,g)) & - & + ff3(3,S) * (lam(i,j,k-3,g)+lam(i,j,k+3,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end do - - if (Er(reg_l1,reg_l2,reg_l3-1,g) .eq. -1.e0_rt) then - k = reg_l3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff3b0, lam(i,j,k:k+3,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_l3 + 1 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff3b1, lam(i,j,k-1:k+3,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_l3 + 2 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff3b2, lam(i,j,k-2:k+3,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - if (Er(reg_l1,reg_l2,reg_h3+1,g) .eq. -1.e0_rt) then - k = reg_h3 - 2 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff3b2(3:-2:-1), lam(i,j,k-3:k+2,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_h3 - 1 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff3b1(3:-1:-1), lam(i,j,k-3:k+1,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff3b0(3:0:-1), lam(i,j,k-3:k,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - lam(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3,g) = & - lamfil(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3) - - else if (filter_T .eq. 4) then - - do k=lam_l3, lam_h3 - do j=lam_l2, lam_h2 - if (Er(reg_l1,j,k,g) .eq. -1.e0_rt) then - lamfil(:,j,k) = -1.e-50_rt - cycle - end if - - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff4(0,S) * lam(i,j,k,g) & - & + ff4(1,S) * (lam(i-1,j,k,g)+lam(i+1,j,k,g)) & - & + ff4(2,S) * (lam(i-2,j,k,g)+lam(i+2,j,k,g)) & - & + ff4(3,S) * (lam(i-3,j,k,g)+lam(i+3,j,k,g)) & - & + ff4(4,S) * (lam(i-4,j,k,g)+lam(i+4,j,k,g)) - end do - - if (Er(reg_l1-1,j,k,g) .eq. -1.e0_rt) then - i = reg_l1 - lamfil(i,j,k) = dot_product(ff4b0, lam(i:i+4,j,k,g)) - - i = reg_l1 + 1 - lamfil(i,j,k) = dot_product(ff4b1, lam(i-1:i+4,j,k,g)) - - i = reg_l1 + 2 - lamfil(i,j,k) = dot_product(ff4b2, lam(i-2:i+4,j,k,g)) - - i = reg_l1 + 3 - lamfil(i,j,k) = dot_product(ff4b3, lam(i-3:i+4,j,k,g)) - end if - - if (Er(reg_h1+1,j,k,g) .eq. -1.e0_rt) then - i = reg_h1 - 3 - lamfil(i,j,k) = dot_product(ff4b3(4:-3:-1), lam(i-4:i+3,j,k,g)) - - i = reg_h1 - 2 - lamfil(i,j,k) = dot_product(ff4b2(4:-2:-1), lam(i-4:i+2,j,k,g)) - - i = reg_h1 - 1 - lamfil(i,j,k) = dot_product(ff4b1(4:-1:-1), lam(i-4:i+1,j,k,g)) - - i = reg_h1 - lamfil(i,j,k) = dot_product(ff4b0(4:0:-1), lam(i-4:i,j,k,g)) - end if - end do - end do - - do k=lam_l3, lam_h3 - if (Er(reg_l1,reg_l2,k,g) .eq. -1.e0_rt) then - cycle - end if - - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lam(i,j,k,g) = ff4(0,S) * lamfil(i,j,k) & - & + ff4(1,S) * (lamfil(i,j-1,k)+lamfil(i,j+1,k)) & - & + ff4(2,S) * (lamfil(i,j-2,k)+lamfil(i,j+2,k)) & - & + ff4(3,S) * (lamfil(i,j-3,k)+lamfil(i,j+3,k)) & - & + ff4(4,S) * (lamfil(i,j-4,k)+lamfil(i,j+4,k)) - end do - end do - - if (Er(reg_l1,reg_l2-1,k,g) .eq. -1.e0_rt) then - j = reg_l2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b0, lamfil(i,j:j+4,k)) - end do - - j = reg_l2 + 1 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b1, lamfil(i,j-1:j+4,k)) - end do - - j = reg_l2 + 2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b2, lamfil(i,j-2:j+4,k)) - end do - - j = reg_l2 + 3 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b3, lamfil(i,j-3:j+4,k)) - end do - end if - - if (Er(reg_l1,reg_h2+1,k,g) .eq. -1.e0_rt) then - j = reg_h2 - 3 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b3(4:-3:-1), lamfil(i,j-4:j+3,k)) - end do - - j = reg_h2 - 2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b2(4:-2:-1), lamfil(i,j-4:j+2,k)) - end do - - j = reg_h2 - 1 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b1(4:-1:-1), lamfil(i,j-4:j+1,k)) - end do - - j = reg_h2 - do i=reg_l1,reg_h1 - lam(i,j,k,g) = dot_product(ff4b0(4:0:-1), lamfil(i,j-4:j,k)) - end do - end if - end do - - do k=reg_l3, reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = ff4(0,S) * lam(i,j,k,g) & - & + ff4(1,S) * (lam(i,j,k-1,g)+lam(i,j,k+1,g)) & - & + ff4(2,S) * (lam(i,j,k-2,g)+lam(i,j,k+2,g)) & - & + ff4(3,S) * (lam(i,j,k-3,g)+lam(i,j,k+3,g)) & - & + ff4(4,S) * (lam(i,j,k-4,g)+lam(i,j,k+4,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end do - - if (Er(reg_l1,reg_l2,reg_l3-1,g) .eq. -1.e0_rt) then - k = reg_l3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b0, lam(i,j,k:k+4,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_l3 + 1 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b1, lam(i,j,k-1:k+4,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_l3 + 2 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b2, lam(i,j,k-2:k+4,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_l3 + 3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b3, lam(i,j,k-3:k+4,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - if (Er(reg_l1,reg_l2,reg_h3+1,g) .eq. -1.e0_rt) then - k = reg_h3 - 3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b3(4:-3:-1), lam(i,j,k-4:k+3,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_h3 - 2 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b2(4:-2:-1), lam(i,j,k-4:k+2,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_h3 - 1 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b1(4:-1:-1), lam(i,j,k-4:k+1,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - - k = reg_h3 - do j=reg_l2, reg_h2 - do i=reg_l1, reg_h1 - lamfil(i,j,k) = dot_product(ff4b0(4:0:-1), lam(i,j,k-4:k,g)) - lamfil(i,j,k) = min(1.e0_rt/3.e0_rt, max(1.e-25_rt, lamfil(i,j,k))) - end do - end do - end if - - lam(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3,g) = & - lamfil(reg_l1:reg_h1,reg_l2:reg_h2,reg_l3:reg_h3) - - end if - - ! lo-x lo-y lo-z - do k=lam_l3,reg_l3-1 - do j=lam_l2,reg_l2-1 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,reg_l2,reg_l3,g) - end if - end do - end do - end do - - ! reg-x lo-y lo-z - do k=lam_l3,reg_l3-1 - do j=lam_l2,reg_l2-1 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,reg_l2,reg_l3,g) - end if - end do - end do - end do - - ! hi-x lo-y lo-z - do k=lam_l3,reg_l3-1 - do j=lam_l2,reg_l2-1 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,reg_l2,reg_l3,g) - end if - end do - end do - end do - - ! lo-x reg-y lo-z - do k=lam_l3,reg_l3-1 - do j=reg_l2,reg_h2 - do i=lam_l1,reg_l1-1 - lam(i,j,k,g) = lam(reg_l1,j,reg_l3,g) - end do - end do - end do - - ! reg-x reg-y lo-z side - do k=lam_l3,reg_l3-1 - do j=reg_l2,reg_h2 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,j,reg_l3,g) - end if - end do - end do - end do - - ! hi-x reg-y lo-z - do k=lam_l3,reg_l3-1 - do j=reg_l2,reg_h2 - do i=reg_h1+1,lam_h1 - lam(i,j,k,g) = lam(reg_h1,j,reg_l3,g) - end do - end do - end do - - ! lo-x hi-y lo-z - do k=lam_l3,reg_l3-1 - do j=reg_h2+1,lam_h2 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,reg_h2,reg_l3,g) - end if - end do - end do - end do - - ! reg-x hi-y lo-z - do k=lam_l3,reg_l3-1 - do j=reg_h2+1,lam_h2 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,reg_h2,reg_l3,g) - end if - end do - end do - end do - - ! hi-x hi-y lo-z - do k=lam_l3,reg_l3-1 - do j=reg_h2+1,lam_h2 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,reg_h2,reg_l3,g) - end if - end do - end do - end do - - ! lo-x lo-y reg-z - do k=reg_l3,reg_h3 - do j=lam_l2,reg_l2-1 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,reg_l2,k,g) - end if - end do - end do - end do - - ! reg-x lo-y reg-z - do k=reg_l3,reg_h3 - do j=lam_l2,reg_l2-1 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,reg_l2,k,g) - end if - end do - end do - end do - - ! hi-x lo-y reg-z - do k=reg_l3,reg_h3 - do j=lam_l2,reg_l2-1 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,reg_l2,k,g) - end if - end do - end do - end do - - ! lo-x reg-y reg-z - do k=reg_l3,reg_h3 - do j=reg_l2,reg_h2 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,j,k,g) - end if - end do - end do - end do - - ! hi-x reg-y reg-z - do k=reg_l3,reg_h3 - do j=reg_l2,reg_h2 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,j,k,g) - end if - end do - end do - end do - - ! lo-x hi-y reg-z - do k=reg_l3,reg_h3 - do j=reg_h2+1,lam_h2 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,reg_h2,k,g) - end if - end do - end do - end do - - ! reg-x hi-y reg-z - do k=reg_l3,reg_h3 - do j=reg_h2+1,lam_h2 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,reg_h2,k,g) - end if - end do - end do - end do - - ! hi-x hi-y reg-z - do k=reg_l3,reg_h3 - do j=reg_h2+1,lam_h2 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,reg_h2,k,g) - end if - end do - end do - end do - - ! lo-x lo-y hi-z - do k=reg_h3+1,lam_h3 - do j=lam_l2,reg_l2-1 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,reg_l2,reg_h3,g) - end if - end do - end do - end do - - ! reg-x lo-y hi-z - do k=reg_h3+1,lam_h3 - do j=lam_l2,reg_l2-1 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,reg_l2,reg_h3,g) - end if - end do - end do - end do - - ! hi-x lo-y hi-z - do k=reg_h3+1,lam_h3 - do j=lam_l2,reg_l2-1 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,reg_l2,reg_h3,g) - end if - end do - end do - end do - - ! lo-x reg-y hi-z - do k=reg_h3+1,lam_h3 - do j=reg_l2,reg_h2 - do i=lam_l1,reg_l1-1 - lam(i,j,k,g) = lam(reg_l1,j,reg_h3,g) - end do - end do - end do - - ! reg-x reg-y hi-z - do k=reg_h3+1,lam_h3 - do j=reg_l2,reg_h2 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,j,reg_h3,g) - end if - end do - end do - end do - - ! hi-x reg-y hi-z - do k=reg_h3+1,lam_h3 - do j=reg_l2,reg_h2 - do i=reg_h1+1,lam_h1 - lam(i,j,k,g) = lam(reg_h1,j,reg_h3,g) - end do - end do - end do - - ! lo-x hi-y hi-z - do k=reg_h3+1,lam_h3 - do j=reg_h2+1,lam_h2 - do i=lam_l1,reg_l1-1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_l1,reg_h2,reg_h3,g) - end if - end do - end do - end do - - ! reg-x hi-y hi-z - do k=reg_h3+1,lam_h3 - do j=reg_h2+1,lam_h2 - do i=reg_l1,reg_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(i,reg_h2,reg_h3,g) - end if - end do - end do - end do - - ! hi-x hi-y hi-z - do k=reg_h3+1,lam_h3 - do j=reg_h2+1,lam_h2 - do i=reg_h1+1,lam_h1 - if (Er(i,j,k,g).eq.-1.e0_rt) then - lam(i,j,k,g) = lam(reg_h1,reg_h2,reg_h3,g) - end if - end do - end do - end do - - end do ! do g= - - if (filter_T .gt. 0) then - deallocate(lamfil) - end if - - return -end subroutine ca_compute_lamborder - subroutine ca_correct_dterm( & dfx, dfx_l1, dfx_l2, dfx_l3, dfx_h1, dfx_h2, dfx_h3, & dfy, dfy_l1, dfy_l2, dfy_l3, dfy_h1, dfy_h2, dfy_h3, & diff --git a/Source/radiation/MGFLD.cpp b/Source/radiation/MGFLD.cpp index cc42d401bc..595f14283a 100644 --- a/Source/radiation/MGFLD.cpp +++ b/Source/radiation/MGFLD.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -1001,16 +1002,256 @@ void Radiation::compute_limiter(int level, const BoxArray& grids, const Real* dx = parent->Geom(level).CellSize(); + using namespace filter; + #ifdef _OPENMP #pragma omp parallel #endif for (MFIter mfi(Er_wide,false); mfi.isValid(); ++mfi) { - BL_FORT_PROC_CALL(CA_COMPUTE_LAMBORDER, ca_compute_lamborder) - (BL_TO_FORTRAN(Er_wide[mfi]), - BL_TO_FORTRAN(kpr[mfi]), - BL_TO_FORTRAN(lamborder[mfi]), - dx, &ngrow, &radiation::limiter, &filter_lambda_T, &filter_lambda_S); - } + FArrayBox lamfilxfab; +#if AMREX_SPACEDIM >= 2 + FArrayBox lamfilyfab; +#endif +#if AMREX_SPACEDIM == 3 + FArrayBox lamfilzfab; +#endif + + Box lambx = lamborder[mfi].box(); + Box bx = grow(lambx, -ngrow); + + if (filter_lambda_T > 0) { + lamfilxfab.resize(bx); +#if AMREX_SPACEDIM >= 2 + lamfilyfab.resize(bx); +#endif +#if AMREX_SPACEDIM == 3 + lamfilzfab.resize(bx); +#endif + } + + Array4 const lam = lamborder[mfi].array(); + Array4 const Er = Er_wide[mfi].array(); + Array4 const kap = kpr[mfi].array(); + Array4 const lamfilx = lamfilxfab.array(); +#if AMREX_SPACEDIM >= 2 + Array4 const lamfily = lamfilyfab.array(); +#endif +#if AMREX_SPACEDIM == 3 + Array4 const lamfilz = lamfilzfab.array(); +#endif + + for (int g = 0; g < Radiation::nGroups; ++g) { + + amrex::Loop(lambx, [=] (int i, int j, int k) noexcept + { + lam(i,j,k,g) = -1.0e50_rt; + + // The radiation energy Er being equal to -1 is our arbitrary + // convention that we're on a boundary zone. This is meaningless + // from the perspective of computing the limiter, so we skip this + // step on all boundary zones. On all zones adjacent to the boundary + // we use a one-sided difference. + + if (Er(i,j,k,g) == -1.0_rt) { + return; + } + + Real r1 = 0.0; + Real r2 = 0.0; + Real r3 = 0.0; + + if (Er(i-1,j,k,g) == -1.0_rt) { + r1 = (Er(i+1,j,k,g) - Er(i,j,k,g)) / (dx[0]); + } + else if (Er(i+1,j,k,g) == -1.0_rt) { + r1 = (Er(i,j,k,g) - Er(i-1,j,k,g)) / (dx[0]); + } + else { + r1 = (Er(i+1,j,k,g) - Er(i-1,j,k,g)) / (2.e0_rt*dx[0]); + } + +#if AMREX_SPACEDIM >= 2 + if (Er(i,j-1,k,g) == -1.0_rt) { + r2 = (Er(i,j+1,k,g) - Er(i,j,k,g)) / (dx[1]); + } + else if (Er(i,j+1,k,g) == -1.0_rt) { + r2 = (Er(i,j,k,g) - Er(i,j-1,k,g)) / (dx[1]); + } + else { + r2 = (Er(i,j+1,k,g) - Er(i,j-1,k,g)) / (2.e0_rt*dx[1]); + } +#endif + +#if AMREX_SPACEDIM == 3 + if (Er(i,j,k-1,g) == -1.0_rt) { + r3 = (Er(i,j,k+1,g) - Er(i,j,k,g)) / dx[2]; + } + else if (Er(i,j,k+1,g) == -1.0_rt) { + r3 = (Er(i,j,k,g) - Er(i,j,k-1,g)) / dx[2]; + } + else { + r3 = (Er(i,j,k+1,g) - Er(i,j,k-1,g)) / (2.e0_rt*dx[2]); + } +#endif + + Real r = std::sqrt(r1 * r1 + r2 * r2 + r3 * r3); + r = r / (kap(i,j,k,g) * std::max(Er(i,j,k,g), 1.e-50_rt)); + + lam(i,j,k,g) = FLDlambda(r); + }); + + // filter + + int lam_ilo = lambx.loVect3d()[0]; + int lam_ihi = lambx.hiVect3d()[0]; + int lam_jlo = lambx.loVect3d()[1]; + int lam_jhi = lambx.hiVect3d()[1]; + int lam_klo = lambx.loVect3d()[2]; + int lam_khi = lambx.hiVect3d()[2]; + + int reg_ilo = bx.loVect3d()[0]; + int reg_ihi = bx.hiVect3d()[0]; + int reg_jlo = bx.loVect3d()[1]; + int reg_jhi = bx.hiVect3d()[1]; + int reg_klo = bx.loVect3d()[2]; + int reg_khi = bx.hiVect3d()[2]; + + if (filter_lambda_T > 0) { + // initialize + + amrex::Loop(bx, [=] (int i, int j, int k) noexcept + { + lamfilx(i,j,k) = -1.0e-50_rt; +#if AMREX_SPACEDIM >= 2 + lamfily(i,j,k) = -1.0e-50_rt; +#endif +#if AMREX_SPACEDIM == 3 + lamfilz(i,j,k) = -1.0e-50_rt; +#endif + }); + + // filter in the x direction + + int dir = 0; + + amrex::Loop(bx, [=] (int i, int j, int k) noexcept + { + if (filter_lambda_T == 1) { + lamfilx(i,j,k) = apply_filter<1>(Er, lam, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 2) { + lamfilx(i,j,k) = apply_filter<2>(Er, lam, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 3) { + lamfilx(i,j,k) = apply_filter<3>(Er, lam, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 4) { + lamfilx(i,j,k) = apply_filter<4>(Er, lam, dir, filter_lambda_S, i, j, k, g); + } + }); + + // filter in the y direction + +#if AMREX_SPACEDIM >= 2 + dir = 1; + + amrex::Loop(bx, [=] (int i, int j, int k) noexcept + { + if (filter_lambda_T == 1) { + lamfily(i,j,k) = apply_filter<1>(Er, lamfilx, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 2) { + lamfily(i,j,k) = apply_filter<2>(Er, lamfilx, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 3) { + lamfily(i,j,k) = apply_filter<3>(Er, lamfilx, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 4) { + lamfily(i,j,k) = apply_filter<4>(Er, lamfilx, dir, filter_lambda_S, i, j, k, g); + } + }); +#endif + + // filter in the z direction + +#if AMREX_SPACEDIM == 3 + dir = 2; + + amrex::Loop(bx, [=] (int i, int j, int k) noexcept + { + if (filter_lambda_T == 1) { + lamfilz(i,j,k) = apply_filter<1>(Er, lamfily, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 2) { + lamfilz(i,j,k) = apply_filter<2>(Er, lamfily, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 3) { + lamfilz(i,j,k) = apply_filter<3>(Er, lamfily, dir, filter_lambda_S, i, j, k, g); + } + else if (filter_lambda_T == 4) { + lamfilz(i,j,k) = apply_filter<4>(Er, lamfily, dir, filter_lambda_S, i, j, k, g); + } + }); +#endif + + // store the final results + + amrex::Loop(bx, [=] (int i, int j, int k) noexcept + { +#if AMREX_SPACEDIM == 1 + lam(i,j,k,g) = lamfilx(i,j,k); +#elif AMREX_SPACEDIM == 2 + lam(i,j,k,g) = lamfily(i,j,k); +#else + lam(i,j,k,g) = lamfilz(i,j,k); +#endif + }); + } // filter_lambda_T + + // For all zones outside the valid region, do a piecewise constant copy + // from the nearest valid edge or corner. + + amrex::Loop(lambx, [=] (int i, int j, int k) noexcept + { + // Only apply this step if we're outside the valid region. + + if (Er(i,j,k,g) != -1.0_rt) { + return; + } + + // Assume in each dimension that we're copying from the same index, + // then shift to the nearest edge in that dimension if we're outside + // the box. + + int ic = i; + int jc = j; + int kc = k; + + if (i < reg_ilo) { + ic = reg_ilo; + } + else if (i > reg_ihi) { + ic = reg_ihi; + } + + if (j < reg_jlo) { + jc = reg_jlo; + } + else if (j > reg_jhi) { + jc = reg_jhi; + } + + if (k < reg_klo) { + kc = reg_klo; + } + else if (k > reg_khi) { + kc = reg_khi; + } + + lam(i,j,k,g) = lam(ic,jc,kc,g); + }); + } // nGroups + } // mfiter if (filter_lambda_T) { lamborder.FillBoundary(parent->Geom(level).periodicity()); diff --git a/Source/radiation/RAD_F.H b/Source/radiation/RAD_F.H index bb8e1f5d86..f77512daa0 100644 --- a/Source/radiation/RAD_F.H +++ b/Source/radiation/RAD_F.H @@ -53,11 +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)); -BL_FORT_PROC_DECL(CA_COMPUTE_LAMBORDER, ca_compute_lamborder) - (const BL_FORT_FAB_ARG(Er), const BL_FORT_FAB_ARG(kap), - BL_FORT_FAB_ARG(lam), const amrex::Real* dx, const int* ngrow, const int* limiter, - const int* filter_lambda_T, const int* filter_lambda_S); - #ifdef __cplusplus extern "C" { #endif diff --git a/Source/radiation/filter.H b/Source/radiation/filter.H index 252d9b8185..3191599cad 100644 --- a/Source/radiation/filter.H +++ b/Source/radiation/filter.H @@ -5,7 +5,6 @@ namespace filter { - // 3-point filter: T=1, R=S=0 AMREX_GPU_HOST_DEVICE AMREX_INLINE Real ff1 (int i) @@ -400,6 +399,277 @@ namespace filter } } + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void load_filter_stencil(Array4 const input, + Array1D& stencil, + int dir, int i, int j, int k, int g) + { + stencil(0) = input(i,j,k,g); + + if constexpr (width >= 1) { + if (dir == 0) { + stencil(-1) = input(i-1,j,k,g); + stencil( 1) = input(i+1,j,k,g); + } + else if (dir == 1) { + stencil(-1) = input(i,j-1,k,g); + stencil( 1) = input(i,j+1,k,g); + } + else { + stencil(-1) = input(i,j,k-1,g); + stencil( 1) = input(i,j,k+1,g); + } + } + + if constexpr (width >= 2) { + if (dir == 0) { + stencil(-2) = input(i-2,j,k,g); + stencil( 2) = input(i+2,j,k,g); + } + else if (dir == 1) { + stencil(-2) = input(i,j-2,k,g); + stencil( 2) = input(i,j+2,k,g); + } + else { + stencil(-2) = input(i,j,k-2,g); + stencil( 2) = input(i,j,k+2,g); + } + } + + if constexpr (width >= 3) { + if (dir == 0) { + stencil(-3) = input(i-3,j,k,g); + stencil( 3) = input(i+3,j,k,g); + } + else if (dir == 1) { + stencil(-3) = input(i,j-3,k,g); + stencil( 3) = input(i,j+3,k,g); + } + else { + stencil(-3) = input(i,j,k-3,g); + stencil( 3) = input(i,j,k+3,g); + } + } + + if constexpr (width >= 4) { + if (dir == 0) { + stencil(-4) = input(i-4,j,k,g); + stencil( 4) = input(i+4,j,k,g); + } + else if (dir == 1) { + stencil(-4) = input(i,j-4,k,g); + stencil( 4) = input(i,j+4,k,g); + } + else { + stencil(-4) = input(i,j,k-4,g); + stencil( 4) = input(i,j,k+4,g); + } + } + } + + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + Real apply_filter (Array4 const Er, + Array4 const lam, + int dir, int S, int i, int j, int k, int g) + { + Real lamfil; + + // In whatever direction we're filtering, retrieve a stencil for both + // lambda and the radiation energy, Er. The energy does not directly + // factor into the calculation, but helps us understand when we have + // gone over the edge of the grid. In the code below, we'll look at + // the stencil only for those zones where Er != -1, which constitutes + // our valid data. The approach we will use is to look progressively + // one zone further out to the left or right; if we detect that any of + // the zones we would normally use for our centered stencil are an + // invalid zone, we use a one-sided or skewed stencil with only zones + // that are available. Note that this relies on the assumption that if + // Er(-1) is -1, then all zones to the left of that are also -1, so we + // don't bother actually examining the contents of those zones (and + // the analogous assumption holds on the right edge of the grid). + + Array1D Er_stencil, lam_stencil; + + load_filter_stencil(Er, Er_stencil, dir, i, j, k, g); + load_filter_stencil(lam, lam_stencil, dir, i, j, k, g); + + if constexpr (width == 1) { + if (Er_stencil(-1) == -1.0_rt) { + lamfil = ff1b(0) * lam_stencil(0) + + ff1b(1) * lam_stencil(1); + } + else if (Er_stencil(1) == -1.0_rt) { + lamfil = ff1b(1) * lam_stencil(-1) + + ff1b(0) * lam_stencil( 0); + } + else { + lamfil = ff1(0) * lam_stencil(0) + + ff1(1) * (lam_stencil(-1) + lam_stencil(1)); + } + } + else if constexpr (width == 2) { + if (Er_stencil(-1) == -1.0_rt) { + lamfil = ff2b0(0) * lam_stencil(0) + + ff2b0(1) * lam_stencil(1) + + ff2b0(2) * lam_stencil(2); + } + else if (Er_stencil(-2) == -1.0_rt) { + lamfil = ff2b1(-1) * lam_stencil(-1) + + ff2b1( 0) * lam_stencil( 0) + + ff2b1( 1) * lam_stencil( 1) + + ff2b1( 2) * lam_stencil( 2); + } + else if (Er_stencil(1) == -1.0_rt) { + lamfil = ff2b0(2) * lam_stencil(-2) + + ff2b0(1) * lam_stencil(-1) + + ff2b0(0) * lam_stencil( 0); + } + else if (Er_stencil(2) == -1.0_rt) { + lamfil = ff2b1( 2) * lam_stencil(-2) + + ff2b1( 1) * lam_stencil(-1) + + ff2b1( 0) * lam_stencil( 0) + + ff2b1(-1) * lam_stencil( 1); + } + else { + lamfil = ff2(0,S) * lam_stencil(0) + + ff2(1,S) * (lam_stencil(-1) + lam_stencil(1)) + + ff2(2,S) * (lam_stencil(-2) + lam_stencil(2)); + } + } + else if constexpr (width == 3) { + if (Er_stencil(-1) == -1.0_rt) { + lamfil = ff3b0(0) * lam_stencil(0) + + ff3b0(1) * lam_stencil(1) + + ff3b0(2) * lam_stencil(2) + + ff3b0(3) * lam_stencil(3); + } + else if (Er_stencil(-2) == -1.0_rt) { + lamfil = ff3b1(-1) * lam_stencil(-1) + + ff3b1( 0) * lam_stencil( 0) + + ff3b1( 1) * lam_stencil( 1) + + ff3b1( 2) * lam_stencil( 2) + + ff3b1( 3) * lam_stencil( 3); + } + else if (Er_stencil(-3) == -1.0_rt) { + lamfil = ff3b2(-2) * lam_stencil(-2) + + ff3b2(-1) * lam_stencil(-1) + + ff3b2( 0) * lam_stencil( 0) + + ff3b2( 1) * lam_stencil( 1) + + ff3b2( 2) * lam_stencil( 2) + + ff3b2( 3) * lam_stencil( 3); + } + else if (Er_stencil(1) == -1.0_rt) { + lamfil = ff3b0(3) * lam_stencil(-3) + + ff3b0(2) * lam_stencil(-2) + + ff3b0(1) * lam_stencil(-1) + + ff3b0(0) * lam_stencil( 0); + } + else if (Er_stencil(2) == -1.0_rt) { + lamfil = ff3b1( 3) * lam_stencil(-3) + + ff3b1( 2) * lam_stencil(-2) + + ff3b1( 1) * lam_stencil(-1) + + ff3b1( 0) * lam_stencil( 0) + + ff3b1(-1) * lam_stencil( 1); + } + else if (Er_stencil(3) == -1.0_rt) { + lamfil = ff3b2( 3) * lam_stencil(-3) + + ff3b2( 2) * lam_stencil(-2) + + ff3b2( 1) * lam_stencil(-1) + + ff3b2( 0) * lam_stencil( 0) + + ff3b2(-1) * lam_stencil( 1) + + ff3b2(-2) * lam_stencil( 2); + } + else { + lamfil = ff3(0,S) * lam_stencil(0) + + ff3(1,S) * (lam_stencil(-1) + lam_stencil(1)) + + ff3(2,S) * (lam_stencil(-2) + lam_stencil(2)) + + ff3(3,S) * (lam_stencil(-3) + lam_stencil(3)); + } + } + else if constexpr (width == 4) { + if (Er_stencil(-1) == -1.0_rt) { + lamfil = ff4b0(0) * lam_stencil(0) + + ff4b0(1) * lam_stencil(1) + + ff4b0(2) * lam_stencil(2) + + ff4b0(3) * lam_stencil(3) + + ff4b0(4) * lam_stencil(4); + } + else if (Er_stencil(-2) == -1.0_rt) { + lamfil = ff4b1(-1) * lam_stencil(-1) + + ff4b1( 0) * lam_stencil( 0) + + ff4b1( 1) * lam_stencil( 1) + + ff4b1( 2) * lam_stencil( 2) + + ff4b1( 3) * lam_stencil( 3) + + ff4b1( 4) * lam_stencil( 4); + } + else if (Er_stencil(-3) == -1.0_rt) { + lamfil = ff4b2(-2) * lam_stencil(-2) + + ff4b2(-1) * lam_stencil(-1) + + ff4b2( 0) * lam_stencil( 0) + + ff4b2( 1) * lam_stencil( 1) + + ff4b2( 2) * lam_stencil( 2) + + ff4b2( 3) * lam_stencil( 3) + + ff4b2( 4) * lam_stencil( 4); + } + else if (Er_stencil(-4) == -1.0_rt) { + lamfil = ff4b3(-3) * lam_stencil(-3) + + ff4b3(-2) * lam_stencil(-2) + + ff4b3(-1) * lam_stencil(-1) + + ff4b3( 0) * lam_stencil( 0) + + ff4b3( 1) * lam_stencil( 1) + + ff4b3( 2) * lam_stencil( 2) + + ff4b3( 3) * lam_stencil( 3) + + ff4b3( 4) * lam_stencil( 4); + } + else if (Er_stencil(1) == -1.0_rt) { + lamfil = ff4b0(4) * lam_stencil(-4) + + ff4b0(3) * lam_stencil(-3) + + ff4b0(2) * lam_stencil(-2) + + ff4b0(1) * lam_stencil(-1) + + ff4b0(0) * lam_stencil( 0); + } + else if (Er_stencil(2) == -1.0_rt) { + lamfil = ff4b1( 4) * lam_stencil(-4) + + ff4b1( 3) * lam_stencil(-3) + + ff4b1( 2) * lam_stencil(-2) + + ff4b1( 1) * lam_stencil(-1) + + ff4b1( 0) * lam_stencil( 0) + + ff4b1(-1) * lam_stencil( 1); + } + else if (Er_stencil(3) == -1.0_rt) { + lamfil = ff4b2( 4) * lam_stencil(-4) + + ff4b2( 3) * lam_stencil(-3) + + ff4b2( 2) * lam_stencil(-2) + + ff4b2( 1) * lam_stencil(-1) + + ff4b2( 0) * lam_stencil( 0) + + ff4b2(-1) * lam_stencil( 1) + + ff4b2(-2) * lam_stencil( 2); + } + else if (Er_stencil(4) == -1.0_rt) { + lamfil = ff4b3( 4) * lam_stencil(-4) + + ff4b3( 3) * lam_stencil(-3) + + ff4b3( 2) * lam_stencil(-2) + + ff4b3( 1) * lam_stencil(-1) + + ff4b3( 0) * lam_stencil( 0) + + ff4b3(-1) * lam_stencil( 1) + + ff4b3(-2) * lam_stencil( 2) + + ff4b3(-3) * lam_stencil( 3); + } + else { + lamfil = ff4(0,S) * lam_stencil(0) + + ff4(1,S) * (lam_stencil(-1) + lam_stencil(1)) + + ff4(2,S) * (lam_stencil(-2) + lam_stencil(2)) + + ff4(3,S) * (lam_stencil(-3) + lam_stencil(3)) + + ff4(4,S) * (lam_stencil(-4) + lam_stencil(4)); + } + } + + lamfil = std::min(1.0_rt / 3.0_rt, std::max(1.0e-25_rt, lamfil)); + + return lamfil; + } } #endif