Skip to content

Commit

Permalink
Dimension agnostic, GPU-friendly put_radial_phi
Browse files Browse the repository at this point in the history
  • Loading branch information
@maxpkatz
maxpkatz committed Dec 31, 2019
1 parent b386f6f commit 4c613a8
Show file tree
Hide file tree
Showing 6 changed files with 136 additions and 299 deletions.
14 changes: 9 additions & 5 deletions Source/gravity/Gravity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1452,11 +1452,15 @@ Gravity::make_radial_phi(int level, const MultiFab& Rhs, MultiFab& phi, int fill
for (MFIter mfi(phi, TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.growntilebox();
ca_put_radial_phi(bx.loVect(), bx.hiVect(),
domain.loVect(), domain.hiVect(),
dx,&dr, BL_TO_FORTRAN(phi[mfi]),
radial_phi.dataPtr(),geom.ProbLo(),
&n1d,&fill_interior);

#pragma gpu box(bx)
ca_put_radial_phi(AMREX_INT_ANYD(bx.loVect()), AMREX_INT_ANYD(bx.hiVect()),
AMREX_INT_ANYD(domain.loVect()), AMREX_INT_ANYD(domain.hiVect()),
AMREX_REAL_ANYD(dx), dr,
BL_TO_FORTRAN_ANYD(phi[mfi]),
radial_phi.dataPtr(),
AMREX_REAL_ANYD(geom.ProbLo()),
n1d, fill_interior);
}

if (verbose)
Expand Down
88 changes: 0 additions & 88 deletions Source/gravity/Gravity_1d.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,92 +158,4 @@ subroutine ca_put_radial_grav(lo,hi,dx,dr,&

end subroutine ca_put_radial_grav



subroutine ca_put_radial_phi (lo,hi,domlo,domhi,dx,dr,&
phi,p_l1,p_h1, &
radial_phi,problo, &
numpts_1d,fill_interior) bind(C, name="ca_put_radial_phi")

use prob_params_module, only: center
use amrex_constants_module, only: HALF, TWO

use amrex_fort_module, only : rt => amrex_real
implicit none
integer , intent(in ) :: lo(1), hi(1)
integer , intent(in ) :: domlo(1), domhi(1)
real(rt), intent(in ) :: dx(1), dr
real(rt), intent(in ) :: problo(1)

integer , intent(in ) :: numpts_1d
real(rt), intent(in ) :: radial_phi(0:numpts_1d-1)
integer , intent(in ) :: fill_interior

integer , intent(in ) :: p_l1, p_h1
real(rt), intent(inout) :: phi(p_l1:p_h1)

integer :: i, index
real(rt) :: r
real(rt) :: cen, xi, slope, p_lo, p_md, p_hi, minvar, maxvar

! Note that when we interpolate into the ghost cells we use the
! location of the edge, not the cell center.

do i = lo(1), hi(1)
if (i .gt. domhi(1)) then
r = problo(1) + (dble(i ) ) * dx(1) - center(1)
else if (i .lt. domlo(1)) then
r = problo(1) + (dble(i+1) ) * dx(1) - center(1)
else
r = problo(1) + (dble(i )+HALF) * dx(1) - center(1)
end if

index = int(r / dr)

if (index .gt. numpts_1d-1) then
print *,'PUT_RADIAL_PHI: INDEX TOO BIG ', index, ' > ', numpts_1d-1
print *,'AT (i) ',i
print *,'R / DR IS ', r, dr
call castro_error("Error:: Gravity_1d.f90 :: ca_put_radial_phi")
end if

if ( (fill_interior .eq. 1) .or. (i .lt. domlo(1) .or. i.gt.domhi(1)) ) then

cen = (dble(index) + HALF) * dr
xi = r - cen

if (index == 0) then

! Linear interpolation or extrapolation
slope = ( radial_phi(index+1) - radial_phi(index) ) / dr
phi(i) = radial_phi(index) + slope * xi

else if (index == numpts_1d-1) then

! Linear interpolation or extrapolation
slope = ( radial_phi(index) - radial_phi(index-1) ) / dr
phi(i) = radial_phi(index) + slope * xi

else

! Quadratic interpolation
p_hi = radial_phi(index+1)
p_md = radial_phi(index )
p_lo = radial_phi(index-1)
phi(i) = ( p_hi - TWO*p_md + p_lo)*xi**2/(TWO*dr**2) + &
( p_hi - p_lo )*xi /(TWO*dr ) + &
(-p_hi + 26.e0_rt*p_md - p_lo)/24.e0_rt
minvar = min(p_md, min(p_lo,p_hi))
maxvar = max(p_md, max(p_lo,p_hi))
phi(i) = max(phi(i),minvar)
phi(i) = min(phi(i),maxvar)

end if

end if

enddo

end subroutine ca_put_radial_phi

end module gravity_1D_module
91 changes: 0 additions & 91 deletions Source/gravity/Gravity_2d.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,95 +152,4 @@ subroutine ca_put_radial_grav (lo,hi,dx,dr,&

end subroutine ca_put_radial_grav



subroutine ca_put_radial_phi (lo,hi,domlo,domhi,dx,dr,&
phi,p_l1,p_l2,p_h1,p_h2, &
radial_phi,problo, &
numpts_1d,fill_interior) bind(C, name="ca_put_radial_phi")

use prob_params_module, only: center
use amrex_constants_module

use amrex_fort_module, only : rt => amrex_real
implicit none
integer , intent(in ) :: lo(2),hi(2)
integer , intent(in ) :: domlo(2),domhi(2)
real(rt), intent(in ) :: dx(2),dr
real(rt), intent(in ) :: problo(2)

integer , intent(in ) :: numpts_1d
real(rt), intent(in ) :: radial_phi(0:numpts_1d-1)
integer , intent(in ) :: fill_interior

integer , intent(in ) :: p_l1,p_l2,p_h1,p_h2
real(rt), intent(inout) :: phi(p_l1:p_h1,p_l2:p_h2)

integer :: i,j,index
real(rt) :: x,y,r
real(rt) :: cen,xi,slope,p_lo,p_md,p_hi,minvar,maxvar

! Note that when we interpolate into the ghost cells we use the
! location of the edge, not the cell center

do j = lo(2), hi(2)
if (j .gt. domhi(2)) then
y = problo(2) + (dble(j ) ) * dx(2) - center(2)
else if (j .lt. domlo(2)) then
y = problo(2) + (dble(j+1) ) * dx(2) - center(2)
else
y = problo(2) + (dble(j )+HALF) * dx(2) - center(2)
end if
do i = lo(1), hi(1)
if (i .gt. domhi(1)) then
x = problo(1) + (dble(i ) ) * dx(1) - center(1)
else if (i .lt. domlo(1)) then
x = problo(1) + (dble(i+1) ) * dx(1) - center(1)
else
x = problo(1) + (dble(i )+HALF) * dx(1) - center(1)
end if
r = sqrt( x**2 + y**2)
index = int(r/dr)
if (index .gt. numpts_1d-1) then
print *,'PUT_RADIAL_PHI: INDEX TOO BIG ',index,' > ',numpts_1d-1
print *,'AT (i,j) ',i,j
print *,'R / DR IS ',r,dr
call castro_error("Error:: Gravity_2d.f90 :: ca_put_radial_phi")
end if

if ( (fill_interior .eq. 1) .or. &
( j.lt.domlo(2).or.j.gt.domhi(2) .or. &
i.lt.domlo(1).or.i.gt.domhi(1) ) ) then

cen = (dble(index)+HALF)*dr
xi = r - cen
if (index == 0) then
! Linear interpolation or extrapolation
slope = ( radial_phi(index+1) - radial_phi(index) ) / dr
phi(i,j) = radial_phi(index) + slope * xi
else if (index == numpts_1d-1) then
! Linear interpolation or extrapolation
slope = ( radial_phi(index) - radial_phi(index-1) ) / dr
phi(i,j) = radial_phi(index) + slope * xi
else
! Quadratic interpolation
p_hi = radial_phi(index+1)
p_md = radial_phi(index )
p_lo = radial_phi(index-1)
phi(i,j) = &
( p_hi - TWO*p_md + p_lo)*xi**2/(TWO*dr**2) + &
( p_hi - p_lo )*xi /(TWO*dr ) + &
(-p_hi + 26.e0_rt*p_md - p_lo)/24.e0_rt
minvar = min(p_md, min(p_lo,p_hi))
maxvar = max(p_md, max(p_lo,p_hi))
phi(i,j) = max(phi(i,j),minvar)
phi(i,j) = min(phi(i,j),maxvar)
end if

end if
enddo
enddo

end subroutine ca_put_radial_phi

end module gravity_2D_module
110 changes: 0 additions & 110 deletions Source/gravity/Gravity_3d.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,114 +143,4 @@ subroutine ca_put_radial_grav (lo,hi,dx,dr,&

end subroutine ca_put_radial_grav



subroutine ca_put_radial_phi (lo,hi,domlo,domhi,dx,dr,&
phi,p_l1,p_l2,p_l3,p_h1,p_h2,p_h3, &
radial_phi,problo,&
numpts_1d,fill_interior) bind(C, name="ca_put_radial_phi")

use amrex_constants_module
use prob_params_module, only: center

use amrex_fort_module, only : rt => amrex_real
implicit none

integer , intent(in ) :: lo(3),hi(3)
integer , intent(in ) :: domlo(3),domhi(3)
real(rt), intent(in ) :: dx(3),dr
real(rt), intent(in ) :: problo(3)

integer , intent(in ) :: numpts_1d
real(rt), intent(in ) :: radial_phi(0:numpts_1d-1)
integer , intent(in ) :: fill_interior

integer , intent(in ) :: p_l1,p_l2,p_l3,p_h1,p_h2,p_h3
real(rt), intent(inout) :: phi(p_l1:p_h1,p_l2:p_h2,p_l3:p_h3)

integer :: i,j,k,index
real(rt) :: x,y,z,r
real(rt) :: cen,xi,slope,p_lo,p_md,p_hi,minvar,maxvar
!
! Note that when we interpolate into the ghost cells we use the
! location of the edge, not the cell center
!
do k = lo(3), hi(3)
if (k .gt. domhi(3)) then
z = problo(3) + (dble(k ) ) * dx(3) - center(3)
else if (k .lt. domlo(3)) then
z = problo(3) + (dble(k+1) ) * dx(3) - center(3)
else
z = problo(3) + (dble(k )+HALF) * dx(3) - center(3)
end if

do j = lo(2), hi(2)
if (j .gt. domhi(2)) then
y = problo(2) + (dble(j ) ) * dx(2) - center(2)
else if (j .lt. domlo(2)) then
y = problo(2) + (dble(j+1) ) * dx(2) - center(2)
else
y = problo(2) + (dble(j )+HALF) * dx(2) - center(2)
end if

do i = lo(1), hi(1)
if (i .gt. domhi(1)) then
x = problo(1) + (dble(i ) ) * dx(1) - center(1)
else if (i .lt. domlo(1)) then
x = problo(1) + (dble(i+1) ) * dx(1) - center(1)
else
x = problo(1) + (dble(i )+HALF) * dx(1) - center(1)
end if

r = sqrt( x**2 + y**2 + z**2 )
index = int(r/dr)

if (index .gt. numpts_1d-1) then
print *,'PUT_RADIAL_PHI: INDEX TOO BIG ',index,' > ',numpts_1d-1
print *,'AT (i,j) ',i,j,k
print *,'R / DR IS ',r,dr
call castro_error("Error:: Gravity_3d.f90 :: ca_put_radial_phi")
end if

if ( (fill_interior .eq. 1) .or. &
( i.lt.domlo(1).or.i.gt.domhi(1) .or. &
j.lt.domlo(2).or.j.gt.domhi(2) .or. &
k.lt.domlo(3).or.k.gt.domhi(3) ) ) then
cen = (dble(index)+HALF)*dr
xi = r - cen
if (index == 0) then
!
! Linear interpolation or extrapolation
!
slope = ( radial_phi(index+1) - radial_phi(index) ) / dr
phi(i,j,k) = radial_phi(index) + slope * xi
else if (index == numpts_1d-1) then
!
! Linear interpolation or extrapolation
!
slope = ( radial_phi(index) - radial_phi(index-1) ) / dr
phi(i,j,k) = radial_phi(index) + slope * xi
else
!
! Quadratic interpolation
!
p_hi = radial_phi(index+1)
p_md = radial_phi(index )
p_lo = radial_phi(index-1)
phi(i,j,k) = &
( p_hi - TWO*p_md + p_lo)*xi**2/(TWO*dr**2) + &
( p_hi - p_lo )*xi /(TWO*dr ) + &
(-p_hi + 26.e0_rt*p_md - p_lo)/24.e0_rt
minvar = min(p_md, min(p_lo,p_hi))
maxvar = max(p_md, max(p_lo,p_hi))
phi(i,j,k) = max(phi(i,j,k),minvar)
phi(i,j,k) = min(phi(i,j,k),maxvar)
end if
end if
enddo
enddo
enddo

end subroutine ca_put_radial_phi

end module gravity_3D_module
10 changes: 5 additions & 5 deletions Source/gravity/Gravity_F.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,13 @@ extern "C"
const int* n1d, const int* level);

void ca_put_radial_phi
(const int lo[], const int hi[],
const int domlo[], const int domhi[],
const amrex::Real* dx, const amrex::Real* dr,
BL_FORT_FAB_ARG(phi),
(const int* lo, const int* hi,
const int* domlo, const int* domhi,
const amrex::Real* dx, amrex::Real dr,
BL_FORT_FAB_ARG_3D(phi),
const amrex::Real* radial_phi,
const amrex::Real* problo,
const int* numpts_1d, const int* fill_interior);
int numpts_1d, int fill_interior);

void init_multipole_gravity
(const int* lnum, const int* lo_bc, const int* hi_bc);
Expand Down
Loading

0 comments on commit 4c613a8

Please sign in to comment.