In this tutorial you will be led through the process of modefying the code you created in the last 2 tutorials so that fluxes are 'fixed' at jumps in refinement so that conservation is maintained.
Since we use finite differences to represent derivatives, there is no guarantee that those derivatives so represented will be continuous at jumps in refinement in our mesh. In typical applications this manifests itself as a non-conservation of some conserved variable.
In this tutorial we describe how to make use of the PARAMESH subroutine 'amr_flux_conserve' which patches the fluxes at refinement jumps.
We start again with our simple 2D diffusion equation, but recast the
equation
in terms of fluxes.
If we define the fluxes at cell interfaces, then a finite difference representation of the above equations is,
U(i,j,t+dt) = U(i,j,t) + dt * { (Fx(i+1/2,j,t) - Fx(i-1/2,j,t)/dx + (Fy(i,j+1/2,t) - Fy(i,j-1/2,t))/dy }
where,
Fx(i+1/2,j,t) = (U(i+1,j,t) - U(i,j,t))/dx
and,
Fy(i,j+1/2,t) = (U(i,j+1,t) - U(i,j,t))/dy
Step 1. - Create a subroutine to compute the fluxes at block boundaries.
subroutine fluxes (ivar,idest,lb,fluxx_x,fluxx_y)
use physicaldata
use paramesh_dimensions
use tree
integer :: ivar, idest, i,j,k, ii,jj,kk
real :: fluxx_x(iu_bnd+1,ju_bnd,ku_bnd)
real :: fluxx_y(iu_bnd,ju_bnd+1,ku_bnd)
real :: dx,dy
dx = bsize(1,lb)/nxb
kk = nguard*k3d+1
do k = 1,ku_bnd
jj = nguard+1
do j = 1,ju_bnd
ii = nguard+1
do i = 1,iu_bnd+1
fluxx_x(i,j,k) = &
(unk1(ivar,ii,jj,kk,idest) - &
unk1(ivar,ii-1,jj,kk,idest))/dx
ii = ii + 1
end do
jj = jj +1
end do
kk = kk + 1
end do
dy = bsize(2,lb)/nyb
kk = nguard*k3d+1
do k = 1,ku_bnd
jj = nguard+1
do j = 1,ju_bnd+1
ii = nguard+1
do i = 1,iu_bnd
fluxx_y(i,j,k) = &
(unk1(ivar,ii,jj,kk,idest) - &
unk1(ivar,ii,jj-1,kk,idest))/dy
ii = ii + 1
end do
jj = jj + 1
end do
kk = kk + 1
end do
return
end subroutine fluxes
Explanation:
This subroutine simply computes the fluxes in the x
and y directions by finite differencing the variable 'U' (assumed to
be stored in the PARAMESH array unk1). The fluxes are computed at
cell faces such that Fx(i-1/2,j,t) -> fluxx_x(i,j) in the above
code.
Note that guardcell storage is not allocated for the fluxes.
Step 2. - Edit the file 'advance_soln_npgc.F90'.
real :: fluxx_x(iu_bnd+1,ju_bnd,ku_bnd)
real :: fluxx_y(iu_bnd,ju_bnd+1,ku_bnd)
call fluxes(1,idest,lb,fluxx_x,fluxx_y)
dx = bsize(1,lb)/nxb
dy = bsize(2,lb)/nyb
kk = 1
do k=1+nguard*k3d,nzb+nguard*k3d
jj = 1
do j=1+nguard*k2d,nyb+nguard*k2d
ii = 1
do i=1+nguard,nxb+nguard
unk1(1,i,j,k,idest) = unk1(1,i,j,k,idest) + dt*( &
(fluxx_x(ii+1,jj,kk) - fluxx_x(ii,jj,kk))/dx + &
(fluxx_y(ii,jj+1,kk) - fluxx_y(ii,jj,kk))/dy)
ii = ii + 1
enddo
jj = jj + 1
enddo
kk = kk + 1
enddo
Explanation:
This is the code to finite difference the fluxes and update the
solution. Note: We no longer need the array 'old_soln'.
Step 3. - Edit the makefile AMRDIR/your_tutorial/Makefile.gnu and add fluxes.F90 to the list of sources.
Step 4. - Re-compile and run it. The changes made should not change the output from the last time you ran it.
gmake -f make_tutor clean (NOTE: not needed if using LIBRARY)
gmake -f make_tutor your_tutorial
mpirun -np X your_tutorial/tutor
Step 5. - Now add in the flux-fix code.
! store fluxes in PARAMESH flux arrays for conservation step
do k = kl_bnd,ku_bnd
do j = jl_bnd,ju_bnd
ix = 1
do i = nguard0+1,nguard0+nxb+1,nxb
flux_x(1,ix,j,k,lb) = fluxx_x(i,j,k)
fluxx_x(i,j,k) = 0.
ix = ix + 1
end do
end do
end do
do k = kl_bnd,ku_bnd
iy = 1
do j = nguard0*k2d+1,nguard0+nyb+1,nyb
do i = il_bnd,iu_bnd
flux_y(1,i,iy,k,lb) = fluxx_y(i,j,k)
fluxx_y(i,j,k) = 0.
end do
iy = iy + 1
end do
end do
Explanation:
This code stores the user computed fluxes
(i.e. 'fluxx_x' and 'fluxx_y') into the PARAMESH arrays (flux_x and
flux_y) which will be altered by the flux-fixing routine,
'amr_flux_conserve'. Note that the arrays 'flux_x' and 'flux_y' are
only defined at the faces of the blocks since these will be the only
fluxes which need to be 'fixed'. For a more detailed explanation see
the Users Guide.
! FIX FLUXES AT BLOCK BOUNDARIES
call amr_flux_conserve(mype,0)
do lb = 1,lnblocks
if (nodetype(lb) == 1 .or. advance_all_levels) then
dx = bsize(1,lb)/real(nxb)
dy = bsize(2,lb)/real(nyb)
do k=1+nguard0*k3d,nzb+nguard0*k3d
do j=1+nguard0*k2d,nyb+nguard0*k2d
ix = 1
do i=1+nguard0,nxb+nguard0,nxb-1
if (ix == 1) then
unk(1,i,j,k,lb) = unk(1,i,j,k,lb) + dt*( &
(-flux_x(1,1,j,k,lb))/dx)
else
unk(1,i,j,k,lb) = unk(1,i,j,k,lb) + dt*( &
(flux_x(1,2,j,k,lb))/dx)
end if
ix = ix + 1
enddo
enddo
enddo
do k=1+nguard0*k3d,nzb+nguard0*k3d
iy = 1
do j=1+nguard0*k2d,nyb+nguard0*k2d,nyb-1
do i=1+nguard0,nxb+nguard0
if (iy == 1) then
unk(1,i,j,k,lb) = unk(1,i,j,k,lb) + dt*( &
(-flux_y(1,i,1,k,lb))/dy)
else
unk(1,i,j,k,lb) = unk(1,i,j,k,lb) + dt*( &
(flux_y(1,i,2,k,lb))/dy)
end if
enddo
iy = iy + 1
enddo
enddo
end if
end do
Explanation:
The fluxes at block faces are stored in 'flux_x' and
'flux_y' are corrected by the call to 'amr_flux_conserve' and then a
second loop over blocks is added. Inside this loop, the corrected
fluxes are applied at all block boundaries (This is perhaps somewhat
inefficient. We only really need to apply the fixed fluxes at jumps
in refinement). Here the fluxes are applied directly to the variables
'unk' rather then 'unk1' as in the previos loop over blocks.
Step 6. - Re-compile and run it. The changes made should not change the output from the last time you ran it.
gmake -f make_tutor clean (NOTE: not needed if using LIBRARY)
gmake -f make_tutor your_tutorial
mpirun -np X your_tutorial/tutor