Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

move the fourth order center <-> average stuff to C++ #884

Merged
merged 11 commits into from
May 6, 2020
Next Next commit
move the fourth order center <-> average stuff to C++
zingale committed Apr 27, 2020
commit 85423e5315f2c8ec50c4ca4196f7b34a3834e5f4
240 changes: 240 additions & 0 deletions Source/hydro/fourth_center_average.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,240 @@


// Note: pretty much all of these routines below assume that dx(1) = dx(2) = dx(3)

Real compute_laplacian(int i, int j, int k, int n,
Array4<Real const> const a,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

Real lapx = 0.0;
Real lapy = 0.0;
Real lapz = 0.0;

// we use 2nd-order accurate one-sided stencils at the physical
// boundaries note: this differs from the suggestion in MC2011 --
// they just suggest using the Laplacian from +1 off the interior.
// I like the one-sided better.

if (i == domlo[0] && physbc_lo[0] != Interior) {
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i+1,j,k,n) + 4.0_rt*a(i+2,j,k,n) - a(i+3,j,k,n);

} else if (i == domhi[0] && physbc_hi[0] != Interior) {
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i-1,j,k,n) + 4.0_rt*a(i-2,j,k,n) - a(i-3,j,k,n);

} else {
lapx = a(i+1,j,k,n) - TWO*a(i,j,k,n) + a(i-1,j,k,n);
}

#if AMREX_SPACEDIM >= 2
if (j == domlo[1] && physbc_lo[1] != Interior) {
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j+1,k,n) + 4.0_rt*a(i,j+2,k,n) - a(i,j+3,k,n);

} else if (j == domhi[1] && physbc_hi[1] != Interior) {
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j-1,k,n) + 4.0_rt*a(i,j-2,k,n) - a(i,j-3,k,n);

} else {
lapy = a(i,j+1,k,n) - TWO*a(i,j,k,n) + a(i,j-1,k,n);
}
#endif

#if AMREX_SPACEDIM == 3
if (k == domlo[2] && physbc_lo[2] != Interior) {
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k+1,n) + 4.0_rt*a(i,j,k+2,n) - a(i,j,k+3,n);

} else if (k == domhi[2] && physbc_hi[2] != Interior) {
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k-1,n) + 4.0_rt*a(i,j,k-2,n) - a(i,j,k-3,n);

} else {
lapz = a(i,j,k+1,n) - TWO*a(i,j,k,n) + a(i,j,k-1,n);
}
#endif

return lapx + lapy + lapz;
}


Real trans_laplacian(int i, int j, int k, int n,
int idir,
Array4<Real const> const a,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

Real lapx = 0.0;
Real lapy = 0.0;
Real lapz = 0.0;

// we use 2nd-order accurate one-sided stencils at the physical boundaries
if (idir != 0) {

if (i == domlo[0] && physbc_lo[0] != Interior) {
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i+1,j,k,n) + 4.0_rt*a(i+2,j,k,n) - a(i+3,j,k,n);

} else if (i == domhi[0] && physbc_hi[0] != Interior) {
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i-1,j,k,n) + 4.0_rt*a(i-2,j,k,n) - a(i-3,j,k,n);

} else {
lapx = a(i+1,j,k,n) - TWO*a(i,j,k,n) + a(i-1,j,k,n);
}
}

if (idir != 1) {

if (j == domlo[1] && physbc_lo[1] != Interior) {
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j+1,k,n) + 4.0_rt*a(i,j+2,k,n) - a(i,j+3,k,n);

} else if (j == domhi[1] && physbc_hi[1] != Interior) {
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j-1,k,n) + 4.0_rt*a(i,j-2,k,n) - a(i,j-3,k,n);

} else {
lapy = a(i,j+1,k,n) - TWO*a(i,j,k,n) + a(i,j-1,k,n);
}
}

#if AMREX_SPACEDIM == 3
if (idir != 2) {

if (k == domlo[2] && physbc_lo[2] != Interior) {
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k+1,n) + 4.0_rt*a(i,j,k+2,n) - a(i,j,k+3,n);

} else if (k == domhi[2] && physbc_hi[2] != Interior) {
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k-1,n) + 4.0_rt*a(i,j,k-2,n) - a(i,j,k-3,n);

} else {
lapz = a(i,j,k+1,n) - TWO*a(i,j,k,n) + a(i,j,k-1,n);
}
}
#endif

return lapx + lapy + lapz;
}


void make_cell_center(const Box& bx,
Array4<Real const> const& U,
Array4<Real> const& U_cc,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

// Take a cell-average state U and a convert it to a cell-center
// state U_cc via U_cc = U - 1/24 L U

U_lo = U.lbound();
U_hi = U.ubound();

AMREX_ASSERT(U_lo[0] <= lo[0]-1 && U_hi[0] >= hi[0]+1 &&
(AMREX_SPACEDIM >= 2 && (U_lo[1] <= lo[1]-1 && U_hi[1] >= hi[1]+1)) &&
(AMREX_SPACEDIM == 3 && (U_lo[2] <= lo[2]-1 && U_hi[2] >= hi[2]+1)));

amrex::ParallelFor(bx, U.nComp(),
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
Real lap = compute_laplacian(i, j, k, n, U, domlo, domhi);

U_cc(i,j,k,n) = U(i,j,k,n) - (1.0_rt/24.0_rt) * lap;

});
}

void make_cell_center_in_place(const Box& bx,
Array4<Real> const& U,
Array4<Real> const& tmp,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

// Take a cell-average state U and make it cell-centered in place
// via U <- U - 1/24 L U. Note that this operation is not tile
// safe.

// here, tmp is a temporary memory space we use to store one-component's Laplacian

for (int n = 0; n < U.nComp(); n++) {

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
tmp(i,j,k) = compute_laplacian(i, j, k, n, U, domlo, domhi);
});

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
U(i,j,k,n) = U(i,j,k,n) - (1.0_rt/24.0_rt) * tmp(i,j,k);
});
}
}


void compute_lap_term(const Box& bx,
Array4<Real const> const& Y,
Array4<Real> const& lap, const int ncomp,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

// Computes the h**2/24 L U term that is used in correcting
// cell-center to averages (and back)

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
lap(i,j,k) = (1.0_rt/24.0_rt) *
compute_laplacian(i, j, k, ncomp, U, domlo, domhi);
});

}


void make_fourth_average(const Box& bx,
Array4<Real> const& q,
Array4<Real const> const& q_bar,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

// Take the cell-center state q and another state q_bar (e.g.,
// constructed from the cell-average U) and replace the cell-center
// q with a 4th-order accurate cell-average, q <- q + 1/24 L q_bar

amrex::ParallelFor(bx, q.nComp(),
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
Real lap = compute_laplacian(i, j, k, n, q_bar, domlo, domhi);

q(i,j,k,n) += (1.0_rt/24.0_rt) * lap;
});
}


void make_fourth_in_place(const Box& bx,
Array4<Real> const& q,
Array4<Real> const& tmp,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

// Take the cell-center q and makes it a cell-average q, in place
// (e.g. q is overwritten by its average), q <- q + 1/24 L q.
// Note: this routine is not tile safe.

for (int n = 0; n < q.nComp(); n++) {
make_fourth_in_place_n(bx, q, n, tmp, domlo, domhi);
}
}


void make_fourth_in_place_n(const Box& bx,
Array4<Real> const& q, const int ncomp,
Array4<Real> const& tmp,
GpuArray<Real, 3>& domlo, GpuArray<Real, 3>& domhi) {

// Take the cell-center q and makes it a cell-average q, in place
// (e.g. q is overwritten by its average), q <- q + 1/24 L q.
// Note: this routine is not tile safe.
//
// This version operates on a single component. Here ncomp is the
// component to update

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
tmp(i,j,k) = compute_laplacian(i, j, k, ncomp, q, domlo, domhi);
});

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
q(i,j,k,ncomp) += (1.0_rt/24.0_rt) * lap(i,j,k);
});

}
410 changes: 0 additions & 410 deletions Source/hydro/fourth_order_nd.F90
Original file line number Diff line number Diff line change
@@ -790,414 +790,4 @@ subroutine ca_states(lo, hi, &

end subroutine ca_states


! Note: pretty much all of these routines below assume that dx(1) = dx(2) = dx(3)
pure function compute_laplacian(i, j, k, n, &
a, a_lo, a_hi, nc, &
domlo, domhi) result (lap)

use prob_params_module, only : physbc_lo, physbc_hi, Interior
implicit none

integer, intent(in) :: i, j, k, n
integer, intent(in) :: a_lo(3), a_hi(3)
integer, intent(in) :: nc
real(rt), intent(in) :: a(a_lo(1):a_hi(1), a_lo(2):a_hi(2), a_lo(3):a_hi(3), nc)
integer, intent(in) :: domlo(3), domhi(3)

real(rt) :: lapx, lapy, lapz
real(rt) :: lap

lapx = ZERO
lapy = ZERO
lapz = ZERO

! we use 2nd-order accurate one-sided stencils at the physical
! boundaries note: this differs from the suggestion in MC2011 --
! they just suggest using the Laplacian from +1 off the interior.
! I like the one-sided better.

if (i == domlo(1) .and. physbc_lo(1) /= Interior) then
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i+1,j,k,n) + 4.0_rt*a(i+2,j,k,n) - a(i+3,j,k,n)

else if (i == domhi(1) .and. physbc_hi(1) /= Interior) then
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i-1,j,k,n) + 4.0_rt*a(i-2,j,k,n) - a(i-3,j,k,n)

else
lapx = a(i+1,j,k,n) - TWO*a(i,j,k,n) + a(i-1,j,k,n)
end if

#if AMREX_SPACEDIM >= 2
if (j == domlo(2) .and. physbc_lo(2) /= Interior) then
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j+1,k,n) + 4.0_rt*a(i,j+2,k,n) - a(i,j+3,k,n)

else if (j == domhi(2) .and. physbc_hi(2) /= Interior) then
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j-1,k,n) + 4.0_rt*a(i,j-2,k,n) - a(i,j-3,k,n)

else
lapy = a(i,j+1,k,n) - TWO*a(i,j,k,n) + a(i,j-1,k,n)
end if

#endif

#if AMREX_SPACEDIM == 3
if (k == domlo(3) .and. physbc_lo(3) /= Interior) then
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k+1,n) + 4.0_rt*a(i,j,k+2,n) - a(i,j,k+3,n)

else if (k == domhi(3) .and. physbc_hi(3) /= Interior) then
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k-1,n) + 4.0_rt*a(i,j,k-2,n) - a(i,j,k-3,n)

else
lapz = a(i,j,k+1,n) - TWO*a(i,j,k,n) + a(i,j,k-1,n)
end if
#endif

lap = lapx + lapy + lapz

end function compute_laplacian

subroutine trans_laplacian(i, j, k, n, &
idir, &
a, a_lo, a_hi, nc, &
lap, &
domlo, domhi) bind(C, name="trans_laplacian")

use prob_params_module, only : physbc_lo, physbc_hi, Interior
implicit none

integer, intent(in), value :: i, j, k, n
integer, intent(in), value :: idir
integer, intent(in) :: a_lo(3), a_hi(3)
integer, intent(in), value :: nc
real(rt), intent(in) :: a(a_lo(1):a_hi(1), a_lo(2):a_hi(2), a_lo(3):a_hi(3), nc)
integer, intent(in) :: domlo(3), domhi(3)
real(rt), intent(out) :: lap

real(rt) :: lapx, lapy, lapz

lapx = ZERO
lapy = ZERO
lapz = ZERO

! we use 2nd-order accurate one-sided stencils at the physical boundaries
if (idir /= 1) then

if (i == domlo(1) .and. physbc_lo(1) /= Interior) then
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i+1,j,k,n) + 4.0_rt*a(i+2,j,k,n) - a(i+3,j,k,n)

else if (i == domhi(1) .and. physbc_hi(1) /= Interior) then
lapx = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i-1,j,k,n) + 4.0_rt*a(i-2,j,k,n) - a(i-3,j,k,n)

else
lapx = a(i+1,j,k,n) - TWO*a(i,j,k,n) + a(i-1,j,k,n)
end if
end if

if (idir /= 2) then

if (j == domlo(2) .and. physbc_lo(2) /= Interior) then
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j+1,k,n) + 4.0_rt*a(i,j+2,k,n) - a(i,j+3,k,n)

else if (j == domhi(2) .and. physbc_hi(2) /= Interior) then
lapy = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j-1,k,n) + 4.0_rt*a(i,j-2,k,n) - a(i,j-3,k,n)

else
lapy = a(i,j+1,k,n) - TWO*a(i,j,k,n) + a(i,j-1,k,n)
end if
end if

#if AMREX_SPACEDIM == 3
if (idir /= 3) then

if (k == domlo(3) .and. physbc_lo(3) /= Interior) then
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k+1,n) + 4.0_rt*a(i,j,k+2,n) - a(i,j,k+3,n)

else if (k == domhi(3) .and. physbc_hi(3) /= Interior) then
lapz = 2.0_rt*a(i,j,k,n) - 5.0_rt*a(i,j,k-1,n) + 4.0_rt*a(i,j,k-2,n) - a(i,j,k-3,n)

else
lapz = a(i,j,k+1,n) - TWO*a(i,j,k,n) + a(i,j,k-1,n)
end if
end if
#endif

lap = lapx + lapy + lapz

end subroutine trans_laplacian


subroutine ca_make_cell_center(lo, hi, &
U, U_lo, U_hi, nc, &
U_cc, U_cc_lo, U_cc_hi, nc_cc, &
domlo, domhi) &
bind(C, name="ca_make_cell_center")
! Take a cell-average state U and a convert it to a cell-center
! state U_cc via U_cc = U - 1/24 L U

implicit none

integer, intent(in) :: lo(3), hi(3)
integer, intent(in) :: U_lo(3), U_hi(3)
integer, intent(in) :: U_cc_lo(3), U_cc_hi(3)
integer, intent(in) :: nc, nc_cc
real(rt), intent(in) :: U(U_lo(1):U_hi(1), U_lo(2):U_hi(2), U_lo(3):U_hi(3), nc)
real(rt), intent(inout) :: U_cc(U_cc_lo(1):U_cc_hi(1), U_cc_lo(2):U_cc_hi(2), U_cc_lo(3):U_cc_hi(3), nc_cc)
integer, intent(in) :: domlo(3), domhi(3)

integer :: i, j, k, n
real(rt) :: lap

if (U_lo(1) > lo(1)-1 .or. U_hi(1) < hi(1)+1 .or. &
(AMREX_SPACEDIM >= 2 .and. (U_lo(2) > lo(2)-1 .or. U_hi(2) < hi(2)+1)) .or. &
(AMREX_SPACEDIM == 3 .and. (U_lo(3) > lo(3)-1 .or. U_hi(3) < hi(3)+1))) then
call bl_error("insufficient ghostcells in ca_make_cell_center")
end if

do n = 1, nc

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)

lap = compute_laplacian(i, j, k, n, &
U, U_lo, U_hi, nc, &
domlo, domhi)

U_cc(i,j,k,n) = U(i,j,k,n) - TWENTYFOURTH * lap

end do
end do
end do
end do

end subroutine ca_make_cell_center

subroutine ca_make_cell_center_in_place(lo, hi, &
U, U_lo, U_hi, nc, &
domlo, domhi) &
bind(C, name="ca_make_cell_center_in_place")
! Take a cell-average state U and make it cell-centered in place
! via U <- U - 1/24 L U. Note that this operation is not tile
! safe.

use amrex_mempool_module, only : bl_allocate, bl_deallocate

implicit none

integer, intent(in) :: lo(3), hi(3)
integer, intent(in) :: U_lo(3), U_hi(3)
integer, intent(in) :: nc
real(rt), intent(inout) :: U(U_lo(1):U_hi(1), U_lo(2):U_hi(2), U_lo(3):U_hi(3), nc)
integer, intent(in) :: domlo(3), domhi(3)

integer :: i, j, k, n

real(rt), pointer :: lap(:,:,:)

call bl_allocate(lap, lo, hi)

do n = 1, nc

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)

lap(i,j,k) = compute_laplacian(i, j, k, n, &
U, U_lo, U_hi, nc, &
domlo, domhi)

end do
end do
end do

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)
U(i,j,k,n) = U(i,j,k,n) - TWENTYFOURTH * lap(i,j,k)
end do
end do
end do
end do

call bl_deallocate(lap)

end subroutine ca_make_cell_center_in_place

subroutine ca_compute_lap_term(lo, hi, &
U, U_lo, U_hi, nc, &
lap, lap_lo, lap_hi, ncomp, &
domlo, domhi) &
bind(C, name="ca_compute_lap_term")
! Computes the h**2/24 L U term that is used in correcting
! cell-center to averages (and back)

implicit none

integer, intent(in) :: lo(3), hi(3)
integer, intent(in) :: U_lo(3), U_hi(3)
integer, intent(in) :: lap_lo(3), lap_hi(3)
integer, intent(in) :: nc
real(rt), intent(in) :: U(U_lo(1):U_hi(1), U_lo(2):U_hi(2), U_lo(3):U_hi(3), nc)
real(rt), intent(inout) :: lap(lap_lo(1):lap_hi(1), lap_lo(2):lap_hi(2), lap_lo(3):lap_hi(3))
integer, intent(in), value :: ncomp
integer, intent(in) :: domlo(3), domhi(3)

! note: ncomp is C++ index (0-based)

integer :: i, j, k

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)

lap(i,j,k) = compute_laplacian(i, j, k, ncomp+1, &
U, U_lo, U_hi, nc, &
domlo, domhi)

end do
end do
end do

lap(lo(1):hi(1), lo(2):hi(2), lo(3):hi(3)) = &
TWENTYFOURTH * lap(lo(1):hi(1), lo(2):hi(2), lo(3):hi(3))

end subroutine ca_compute_lap_term

subroutine ca_make_fourth_average(lo, hi, &
q, q_lo, q_hi, nc, &
q_bar, q_bar_lo, q_bar_hi, nc_bar, &
domlo, domhi) &
bind(C, name="ca_make_fourth_average")
! Take the cell-center state q and another state q_bar (e.g.,
! constructed from the cell-average U) and replace the cell-center
! q with a 4th-order accurate cell-average, q <- q + 1/24 L q_bar

integer, intent(in) :: lo(3), hi(3)
integer, intent(in) :: q_lo(3), q_hi(3)
integer, intent(in) :: q_bar_lo(3), q_bar_hi(3)
integer, intent(in) :: nc, nc_bar
real(rt), intent(inout) :: q(q_lo(1):q_hi(1), q_lo(2):q_hi(2), q_lo(3):q_hi(3), nc)
real(rt), intent(in) :: q_bar(q_bar_lo(1):q_bar_hi(1), q_bar_lo(2):q_bar_hi(2), q_bar_lo(3):q_bar_hi(3), nc_bar)
integer, intent(in) :: domlo(3), domhi(3)

integer :: i, j, k, n
real(rt) :: lap

do n = 1, nc
do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)

lap = compute_laplacian(i, j, k, n, &
q_bar, q_bar_lo, q_bar_hi, nc, &
domlo, domhi)

q(i,j,k,n) = q(i,j,k,n) + TWENTYFOURTH * lap

end do
end do
end do
end do

end subroutine ca_make_fourth_average

subroutine ca_make_fourth_in_place(lo, hi, &
q, q_lo, q_hi, nc, &
domlo, domhi) &
bind(C, name="ca_make_fourth_in_place")
! Take the cell-center q and makes it a cell-average q, in place
! (e.g. q is overwritten by its average), q <- q + 1/24 L q.
! Note: this routine is not tile safe.

use amrex_mempool_module, only : bl_allocate, bl_deallocate

integer, intent(in) :: lo(3), hi(3)
integer, intent(in) :: q_lo(3), q_hi(3)
integer, intent(in) :: nc
real(rt), intent(inout) :: q(q_lo(1):q_hi(1), q_lo(2):q_hi(2), q_lo(3):q_hi(3), nc)
integer, intent(in) :: domlo(3), domhi(3)

integer :: i, j, k, n
real(rt), pointer :: lap(:,:,:)

call bl_allocate(lap, lo, hi)

do n = 1, nc

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)

lap(i,j,k) = compute_laplacian(i, j, k, n, &
q, q_lo, q_hi, nc, &
domlo, domhi)

end do
end do
end do

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)
q(i,j,k,n) = q(i,j,k,n) + TWENTYFOURTH * lap(i,j,k)
end do
end do
end do
end do

call bl_deallocate(lap)

end subroutine ca_make_fourth_in_place

subroutine ca_make_fourth_in_place_n(lo, hi, &
q, q_lo, q_hi, nc, ncomp, &
domlo, domhi) &
bind(C, name="ca_make_fourth_in_place_n")
! Take the cell-center q and makes it a cell-average q, in place
! (e.g. q is overwritten by its average), q <- q + 1/24 L q.
! Note: this routine is not tile safe.
!
! This version operates on a single component. Here ncomp is the
! component to update -- we expect this to come in 0-based from
! C++, so we add 1

use amrex_mempool_module, only : bl_allocate, bl_deallocate

integer, intent(in) :: lo(3), hi(3)
integer, intent(in) :: q_lo(3), q_hi(3)
integer, intent(in) :: nc
integer, intent(in), value :: ncomp
real(rt), intent(inout) :: q(q_lo(1):q_hi(1), q_lo(2):q_hi(2), q_lo(3):q_hi(3), nc)
integer, intent(in) :: domlo(3), domhi(3)

integer :: i, j, k
real(rt), pointer :: lap(:,:,:)

call bl_allocate(lap, lo, hi)

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)

lap(i,j,k) = compute_laplacian(i, j, k, ncomp+1, &
q, q_lo, q_hi, nc, &
domlo, domhi)

end do
end do
end do

do k = lo(3), hi(3)
do j = lo(2), hi(2)
do i = lo(1), hi(1)
q(i,j,k,ncomp+1) = q(i,j,k,ncomp+1) + TWENTYFOURTH * lap(i,j,k)
end do
end do
end do

call bl_deallocate(lap)

end subroutine ca_make_fourth_in_place_n


end module fourth_order