gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Implement magnetic field prolongation with boundary.

Browse Source 
INTERPOLATION

 - add subroutine expand_1d_tvd to 1D TVD prolongation;

 - implement magnetic field boundary prolongation with a given boundary;
   this is used for boundary update of magnetic components in a block at
   higher level, where the boundary is updated from the block at lower
   level; the prolongation is divergence free; for now only 2D case is
   implemented;

BOUNDARY CONDITIONS

 - use the boundary prolongation subroutine expand_mag_bnd in the
   updating the magnetic field boundaries in a blocks which neighbors
   are at lower level; this should be ready for 3D too;
This commit is contained in:
Grzegorz Kowal 2010-07-03 23:37:02 -03:00
parent 4dea14df91
commit c057d098de
2 changed files with 588 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

298
src/boundaries.F90
View File

@ -859,7 +859,7 @@ module boundaries
use error , only : print_warning
use interpolation, only : expand
#if defined MHD && defined FLUXCT
use interpolation, only : expand_mag
use interpolation, only : expand_mag_bnd
#endif /* MHD & FLUXCT */
implicit none
@ -872,21 +872,22 @@ module boundaries
! local variables
!
integer :: q, i, j, k
integer :: q, i, j, k, nh
integer :: i1, i2, j1, j2, k1, k2
integer :: il, iu, jl, ju, kl, ku
integer :: is, it, js, jt, ks, kt
integer :: if, jf, kf
integer :: if, jf, kf, id, jd, kd
integer :: ih, jh, kh
! local arrays
!
integer, dimension(3) :: dm, cm, pm
real , dimension(:,:,:) , allocatable :: u
real , dimension(:,:,:,:), allocatable :: b
real , dimension(:,:,:) , allocatable :: u, bx, by, bz
real , dimension(:,:) , allocatable :: bn
! parameters
!
integer :: del = 1
integer :: dl = 2
!
!-------------------------------------------------------------------------------
!
@ -895,9 +896,9 @@ module boundaries
dm(:) = (/ im, jm, km /)
dm(idir) = ng
cm(:) = dm(:) / 2 + ng
cm(idir) = 3 * ng / 2 + 2 * del
cm(idir) = 3 * ng / 2 + 2 * dl
pm(:) = dm(:)
pm(idir) = ng + 4 * del
pm(idir) = ng + 4 * dl
#if NDIMS == 2
dm(3) = 1
cm(3) = 1
@ -917,11 +918,11 @@ module boundaries
! indices of the input array
!
if (if .eq. 0) then
i1 = ie - ng + 1 - del
i2 = ie + ng / 2 + del
i1 = ie - ng + 1 - dl
i2 = ie + ng / 2 + dl
else
i1 = ib - ng / 2 - del
i2 = ib + ng - 1 + del
i1 = ib - ng / 2 - dl
i2 = ib + ng - 1 + dl
end if
if (jf .eq. 0) then
j1 = 1
@ -945,8 +946,8 @@ module boundaries
! indices for the output array
!
il = 1 + 2 * del
iu = ng + 2 * del
il = 1 + 2 * dl
iu = ng + 2 * dl
jl = 1
ju = jm
kl = 1
@ -982,11 +983,11 @@ module boundaries
i2 = im
end if
if (jf .eq. 0) then
j1 = je - ng + 1 - del
j2 = je + ng / 2 + del
j1 = je - ng + 1 - dl
j2 = je + ng / 2 + dl
else
j1 = jb - ng / 2 - del
j2 = jb + ng - 1 + del
j1 = jb - ng / 2 - dl
j2 = jb + ng - 1 + dl
end if
#if NDIMS == 3
if (kf .eq. 0) then
@ -1005,8 +1006,8 @@ module boundaries
!
il = 1
iu = im
jl = 1 + 2 * del
ju = ng + 2 * del
jl = 1 + 2 * dl
ju = ng + 2 * dl
kl = 1
ku = km
@ -1048,11 +1049,11 @@ module boundaries
j2 = jm
end if
if (kf .eq. 0) then
k1 = ke - ng + 1 - del
k2 = ke + ng / 2 + del
k1 = ke - ng + 1 - dl
k2 = ke + ng / 2 + dl
else
k1 = kb - ng / 2 - del
k2 = kb + ng - 1 + del
k1 = kb - ng / 2 - dl
k2 = kb + ng - 1 + dl
end if
! indices for the output array
@ -1061,8 +1062,8 @@ module boundaries
iu = im
jl = 1
ju = jm
kl = 1 + 2 * del
ku = ng + 2 * del
kl = 1 + 2 * dl
ku = ng + 2 * dl
! indices of the destination array
!
@ -1115,45 +1116,240 @@ module boundaries
end do
#endif /* FIELDCD */
#ifdef FLUXCT
! prepare useful variables
!
nh = max(1, ng / 2)
ih = max(1, in / 2)
jh = max(1, jn / 2)
kh = max(1, kn / 2)
! prepare dimensions of the input arrays
!
dm(:) = (/ im, jm, km /)
dm(idir) = ng
pm(:) = (/ ih, jh, kh /) + ng
pm(idir) = nh
cm(:) = pm(:) + 2 * dl
! allocate space for the prolongated magnetic field components
!
allocate(b(3,pm(1),pm(2),pm(3)))
allocate(bx(dm(1)+1,dm(2) ,dm(3) ))
allocate(by(dm(1) ,dm(2)+1,dm(3) ))
allocate(bz(dm(1) ,dm(2) ,dm(3)+1))
! depending on the direction perform the right prolongation
!
select case(idir)
!! boundary update along the X direction
!!
case(1)
! prepare the index along the longitudinal direction from which the boundary
! should be copied
!
id = ibl + if * in
! prepare the indices in ourder to determine the input array
!
if (if .eq. 0) then
is = ie - nh
it = ie
else
is = ib - 1
it = ib + nh - 1
end if
il = is - dl + 1
iu = it + dl
jl = jb - nh - dl + jf * jh
js = jl - 1
ju = jl + cm(2) - 1
#if NDIMS == 2
kl = 1
ks = 1
ku = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
kl = kb - nh - dl + kf * kh
ks = kl - 1
ku = kl + cm(3) - 1
#endif /* NDIMS == 3 */
! allocate the temporary boundary array
!
allocate(bn(jm,km))
! copy the boundary from the higher refinement level neighbor
!
bn(1:jm,1:km) = pdata%u(ibx,id,1:jm,1:km)
! prolongate magnetic field preserving the divergence-free condition
!
call expand_mag(cm, pm, ng, ub(ibx,i1:i2,j1:j2,k1:k2), ub(iby,i1:i2,j1:j2,k1:k2), ub(ibz,i1:i2,j1:j2,k1:k2), b(1,:,:,:), b(2,:,:,:), b(3,:,:,:))
call expand_mag_bnd(idir, if, dl, pm, bn(:,:) &
, ub(ibx,is:it,jl:ju,kl:ku), ub(iby,il:iu,js:ju,kl:ku) &
, ub(ibz,il:iu,jl:ju,ks:ku), bx(:,:,:), by(:,:,:), bz(:,:,:))
! update the X component of magnetic field
! deallocate the boundary array
!
if (it .eq. ng) then
pdata%u(ibx,is:it-1,js:jt,ks:kt) = b(1,il:iu-1,jl:ju,kl:ku)
else
pdata%u(ibx,is:it ,js:jt,ks:kt) = b(1,il:iu ,jl:ju,kl:ku)
end if
deallocate(bn)
! update the Y component of magnetic field
!
if (jt .eq. ng) then
pdata%u(iby,is:it,js:jt-1,ks:kt) = b(2,il:iu,jl:ju-1,kl:ku)
else
pdata%u(iby,is:it,js:jt ,ks:kt) = b(2,il:iu,jl:ju ,kl:ku)
end if
! update the Z component of magnetic field
! update the magnetic field components
!
if (if .eq. 0) then
pdata%u(ibx, 1:ibl-1,1:jm,1:km) = bx(2:ng ,1:jm ,1:km )
pdata%u(iby, 1:ibl ,1:jm,1:km) = by(1:ng ,2:jm+1,1:km )
#if NDIMS == 3
if (kt .eq. ng) then
pdata%u(ibz,is:it,js:jt,ks:kt-1) = b(3,il:iu,jl:ju,kl:ku-1)
else
pdata%u(ibz,is:it,js:jt,ks:kt ) = b(3,il:iu,jl:ju,kl:ku )
end if
#else /* NDIMS == 3 */
pdata%u(ibz,is:it,js:jt,ks:kt) = b(3,il:iu,jl:ju,kl:ku)
pdata%u(ibz, 1:ibl ,1:jm,1:km) = bz(1:ng ,1:jm ,2:km+1)
#endif /* NDIMS == 3 */
else
pdata%u(ibx,ieu:im ,1:jm,1:km) = bx(2:ng+1,1:jm ,1:km )
pdata%u(iby,ieu:im ,1:jm,1:km) = by(1:ng ,2:jm+1,1:km )
#if NDIMS == 3
pdata%u(ibz,ieu:im ,1:jm,1:km) = bz(1:ng ,1:jm ,2:km+1)
#endif /* NDIMS == 3 */
end if
!! boundary update along the Y direction
!!
case(2)
! prepare the index along the longitudinal direction from which the boundary
! should be copied
!
jd = jbl + jf * jn
! prepare the indices in ourder to determine the input array
!
il = ib - nh - dl + if * ih
is = il - 1
iu = il + cm(1) - 1
if (jf .eq. 0) then
js = je - nh
jt = je
else
js = jb - 1
jt = jb + nh - 1
end if
jl = js - dl + 1
ju = jt + dl
#if NDIMS == 2
kl = 1
ks = 1
ku = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
kl = kb - nh - dl + kf * kh
ks = kl - 1
ku = kl + cm(3) - 1
#endif /* NDIMS == 3 */
! allocate the temporary boundary array
!
allocate(bn(im,km))
! copy the boundary from the higher refinement level neighbor
!
bn(1:im,1:km) = pdata%u(iby,1:im,jd,1:km)
! prolongate magnetic field preserving the divergence-free condition
!
call expand_mag_bnd(idir, jf, dl, pm, bn(:,:) &
, ub(ibx,is:iu,jl:ju,kl:ku), ub(iby,il:iu,js:jt,kl:ku) &
, ub(ibz,il:iu,jl:ju,ks:ku), bx(:,:,:), by(:,:,:), bz(:,:,:))
! deallocate the boundary array
!
deallocate(bn)
! update the magnetic field components
!
if (jf .eq. 0) then
pdata%u(ibx,1:im, 1:jbl ,1:km) = bx(2:im+1,1:ng ,1:km )
pdata%u(iby,1:im, 1:jbl-1,1:km) = by(1:im ,2:ng ,1:km )
#if NDIMS == 3
pdata%u(ibz,1:im, 1:jbl ,1:km) = bz(1:im ,1:ng ,2:km+1)
#endif /* NDIMS == 3 */
else
pdata%u(ibx,1:im,jeu:jm ,1:km) = bx(2:im+1,1:ng ,1:km )
pdata%u(iby,1:im,jeu:jm ,1:km) = by(1:im ,2:ng+1,1:km )
#if NDIMS == 3
pdata%u(ibz,1:im,jeu:jm ,1:km) = bz(1:im ,1:ng ,2:km+1)
#endif /* NDIMS == 3 */
end if
#if NDIMS == 3
!! boundary update along the Y direction
!!
case(3)
! prepare the index along the longitudinal direction from which the boundary
! should be copied
!
kd = kbl + kf * kn
! prepare the indices in ourder to determine the input array
!
il = ib - nh - dl + if * ih
is = il - 1
iu = il + cm(1) - 1
jl = jb - nh - dl + jf * jh
js = jl - 1
ju = jl + cm(2) - 1
if (kf .eq. 0) then
ks = ke - nh
kt = ke
else
ks = kb - 1
kt = kb + nh - 1
end if
kl = ks - dl + 1
ku = kt + dl
! allocate the temporary boundary array
!
allocate(bn(im,jm))
! copy the boundary from the higher refinement level neighbor
!
bn(1:im,1:jm) = pdata%u(ibz,1:im,1:jm,kd)
! prolongate magnetic field preserving the divergence-free condition
!
call expand_mag_bnd(idir, kf, dl, pm, bn(:,:) &
, ub(ibx,is:iu,jl:ju,kl:ku), ub(iby,il:iu,js:ju,kl:ku) &
, ub(ibz,il:iu,jl:ju,ks:kt), bx(:,:,:), by(:,:,:), bz(:,:,:))
! deallocate the boundary array
!
deallocate(bn)
! update the magnetic field components
!
if (kf .eq. 0) then
pdata%u(ibx,1:im,1:jm, 1:kbl ) = bx(2:im+1,1:jm ,1:ng )
pdata%u(iby,1:im,1:jm, 1:kbl ) = by(1:im ,2:jm+1,1:ng )
pdata%u(ibz,1:im,1:jm, 1:kbl-1) = bz(1:im ,1:jm ,2:ng )
else
pdata%u(ibx,1:im,1:jm,keu:km ) = bx(2:im+1,1:jm ,1:ng )
pdata%u(iby,1:im,1:jm,keu:km ) = by(1:im ,2:jm+1,1:ng )
pdata%u(ibz,1:im,1:jm,keu:km ) = bz(1:im ,1:jm ,2:ng+1)
end if
#endif /* NDIMS == 3 */
case default
end select
! deallocate space for the prolongated magnetic field components
!
deallocate(b)
deallocate(bx)
deallocate(by)
deallocate(bz)
#endif /* FLUXCT */
#endif /* MHD */

341
src/interpolation.F90
View File

@ -348,6 +348,56 @@ module interpolation
!
!===============================================================================
!
! expand_1d_tvd: one dimensional expansion using the TVD interpolation
!
!===============================================================================
!
subroutine expand_1d_tvd(nu, nv, ng, u, v)
implicit none
! input parameters
!
integer , intent(in) :: nu, nv, ng
real , dimension(nu), intent(in) :: u
real , dimension(nv), intent(inout) :: v
! local variables
!
integer :: i, ib, ie, il, ir
real :: dum, dup, ds, du
!
!-------------------------------------------------------------------------------
!
ib = 1 + ng
ie = nu - ng
do i = ib, ie
ir = 2 * (i - ng)
il = ir - 1
dup = u(i+1) - u(i)
dum = u(i-1) - u(i)
ds = - dup * dum
if (ds .gt. 0.0) then
du = 0.5 * ds / (dup - dum)
v(il) = u(i) - du
v(ir) = u(i) + du
else
v(il) = u(i)
v(ir) = u(i)
end if
end do
!
!-------------------------------------------------------------------------------
!
end subroutine expand_1d_tvd
!
!===============================================================================
!
! shrink_1d: one dimensional shrinking using different interpolations for
! different location of the new points
!
@ -933,6 +983,297 @@ module interpolation
#endif /* NDIMS == 3 */
!
end subroutine expand_mag
!
!===============================================================================
!
! expand_mag_bnd: subroutine for prolongation of the magnetic field with
! preserving the divergence free condition and applied
! boundaries (Li & Li, 2004, JCoPh, 199, 1)
!
! arguments:
!
! id - the direction of the expansion
! if - the side of the boundary values
! dl - number of safety zones
! dm(3) - the dimensions of the expanded region without boundaries
! bn(:,:) - the finer values of the boundary
! ux(:,:,:) - X component of the coarse field
! uy(:,:,:) - Y component of the coarse field
! uz(:,:,:) - Z component of the coarse field
! wx(:,:,:) - X component of the fine field
! wy(:,:,:) - Y component of the fine field
! wz(:,:,:) - Z component of the fine field
!
!===============================================================================
!
subroutine expand_mag_bnd(id, if, dl, dm, bn, ux, uy, uz, wx, wy, wz)
implicit none
! input and output variables
!
integer , intent(in) :: id, if, dl
integer, dimension(3) , intent(in) :: dm
real , dimension(:,:) , intent(in) :: bn
real , dimension(:,:,:), intent(in) :: ux, uy, uz
real , dimension(:,:,:), intent(inout) :: wx, wy, wz
! local variables
!
integer :: i, j, k, im1, jm1, km1, ip1, jp1, kp1
integer :: ib, ie, jb, je, kb, ke
integer :: il, jl, kl
! local arrays
!
integer, dimension(3) :: cm, fm
! temporary arrays
!
real, dimension(:,:,:), allocatable :: tx, ty
!
!------------------------------------------------------------------------------
!
! notes:
!
! - the dimensions of parallel and perpendicular components are different;
! the parallel components is [ dm(1) + 1, dm(2) + 2 * ng, dm(3) + 2 * ng];
! the perpendicular is [ dm(1) + 2 * ng, dm(2) + 2 * ng + 1, dm(3) + 2 * ng]
! or [ dm(1) + 2 * ng, dm(2) + 2 * ng, dm(3) + 2 * ng + 1]
!
! - the boundary is an array of the dimensions
! [ 4 * dm(2) + 2 * ng, 4 * dm(3) + 2 * ng ]
!
! - the parallel output vector is of the dimensions
! [ 2 * dm(1) + 1, 2 * dm(2) + 2 * ng, 2 * dm(2) + 2 * ng ];
!
! - the perpendicular output vector is of the dimensions
! [ 2 * dm(1), 2 * dm(2) + 2 * ng + 1, 2 * dm(2) + 2 * ng ] or
! [ 2 * dm(1), 2 * dm(2) + 2 * ng, 2 * dm(2) + 2 * ng + 1 ];
!
! - depending on the direction and side of the boundary vector choose the right
! interpolation;
!
! - perpendicular components to the direction of boundary should be larger by
! at least 2 ghost zones in order to handle properly TVD interpolation;
!
! - boundary array should be two dimensional and have perpendicular dimensions
! of the fine output array, e.g. if idir = 1, bn(fm(2),fm(3));
! prepare dimensions
!
cm(:) = dm(:) + 2 * dl
#if NDIMS == 2
cm(3) = 1
#endif /* NDIMS == 2 */
fm(:) = 2 * dm(:)
#if NDIMS == 2
fm(3) = 1
#endif /* NDIMS == 2 */
select case(id)
!! boundary perpendicular to the X direction
!!
case(1)
! allocate temporary arrays
!
allocate(tx(fm(1)+1,cm(2) ,cm(3)))
allocate(ty(fm(1) ,cm(2)+1,cm(3)))
! prepare indices
!
il = 1 + (1 - if) * fm(1)
jb = 1 + dl
je = cm(2) - dl
! steps:
!
! 1. check if the input arrays have the right dimensions
!
! 2. interpolate components perpendicular to the direction of boundary using the
! TVD interpolation; place the calculated values in the output array;
!
do k = 1, cm(3)
do j = 1, cm(2) + 1
call expand_1d_tvd(cm(1), fm(1), dl, uy(1:cm(1),j,k), ty(1:fm(1),j,k))
end do
end do
! 3. fill the already known values of the parallel component
!
tx(1:fm(1)+1:2,1:cm(2),1:cm(3)) = ux(1:dm(1)+1,1:cm(2),1:cm(3))
! 4. average the fine boundary in order to get the coarse values and substitute
! these values in the coarse parallel component;
!
tx(il,jb:je,1:cm(3)) = 0.50 * (bn(1:fm(2):2,1:cm(3)) &
+ bn(2:fm(2):2,1:cm(3)))
! 5. calculate the even values of the parallel component using the divergence
! free condition;
!
do k = 1, cm(3)
do j = 1, cm(2)
jp1 = j + 1
do i = 2, fm(1) + 1, 2
im1 = i - 1
ip1 = i + 1
tx(i,j,k) = 0.50 * (tx(im1,j ,k) + tx(ip1,j ,k)) &
+ 0.25 * (ty(i ,jp1,k) - ty(im1,jp1,k)) &
- 0.25 * (ty(i ,j ,k) - ty(im1,j ,k))
end do
end do
end do
! 6. interpolate components perpendicular to the direction of boundary using the
! TVD interpolation; place the calculated values in the output array;
!
do k = 1, cm(3)
do i = 1, fm(1) + 1
call expand_1d_tvd(cm(2), fm(2), dl, tx(i,1:cm(2),k), wx(i,1:fm(2),k))
end do
end do
! 7. substitute the original values of the fine boundary in the parallel
! component;
!
wx(il,1:fm(2),1:fm(3)) = bn(1:fm(2),1:fm(3))
! 8. substitute the odd values of the perpendicular component from the previous
! array;
!
wy(1:fm(1),1:fm(2)+1:2,1:fm(3)) = ty(1:fm(1),jb:je+1,1:fm(3))
! 9. calculate the even values of the perpendicular component using the
! divergence free condition;
!
do k = 1, fm(3)
do i = 1, fm(1)
ip1 = i + 1
do j = 2, fm(2), 2
jm1 = j - 1
jp1 = j + 1
wy(i,j,k) = 0.50 * (wy(i ,jm1,k) + wy(i ,jp1,k)) &
+ 0.50 * (wx(ip1,j ,k) - wx(ip1,jm1,k)) &
- 0.50 * (wx(i ,j ,k) - wx(i ,jm1,k))
end do
end do
end do
! deallocate temporary arrays
!
deallocate(tx)
deallocate(ty)
!! boundary perpendicular to the Y direction
!!
case(2)
! allocate temporary arrays
!
allocate(tx(cm(1)+1,fm(2) ,cm(3)))
allocate(ty(cm(1) ,fm(2)+1,cm(3)))
! prepare indices
!
jl = 1 + (1 - if) * fm(2)
ib = 1 + dl
ie = cm(1) - dl
! steps:
!
! 1. check if the input arrays have the right dimensions
!
! 2. interpolate components perpendicular to the direction of boundary using the
! TVD interpolation; place the calculated values in the output array;
!
do k = 1, cm(3)
do i = 1, cm(1) + 1
call expand_1d_tvd(cm(2), fm(2), dl, ux(i,1:cm(2),k), tx(i,1:fm(2),k))
end do
end do
! 3. fill the already known values of the parallel component
!
ty(1:cm(1),1:fm(2)+1:2,1:cm(3)) = uy(1:cm(1),1:dm(2)+1,1:cm(3))
! 4. average the fine boundary in order to get the coarse values and substitute
! these values in the coarse parallel component;
!
ty(ib:ie,jl,1:cm(3)) = 0.50 * (bn(1:fm(1):2,1:cm(3)) &
+ bn(2:fm(1):2,1:cm(3)))
! 5. calculate the even values of the parallel component using the divergence
! free condition;
!
do k = 1, cm(3)
do i = 1, cm(1)
ip1 = i + 1
do j = 2, fm(2) + 1, 2
jm1 = j - 1
jp1 = j + 1
ty(i,j,k) = 0.50 * (ty(i ,jm1,k) + ty(i ,jp1,k)) &
+ 0.25 * (tx(ip1,j ,k) - tx(ip1,jm1,k)) &
- 0.25 * (tx(i ,j ,k) - tx(i ,jm1,k))
end do
end do
end do
! 6. interpolate components perpendicular to the direction of boundary using the
! TVD interpolation; place the calculated values in the output array;
!
do k = 1, cm(3)
do j = 1, fm(2) + 1
call expand_1d_tvd(cm(1), fm(1), dl, ty(1:cm(1),j,k), wy(1:fm(1),j,k))
end do
end do
! 7. substitute the original values of the fine boundary in the parallel
! component;
!
wy(1:fm(1),jl,1:fm(3)) = bn(1:fm(1),1:fm(3))
! 8. substitute the odd values of the perpendicular component from the previous
! array;
!
wx(1:fm(1)+1:2,1:fm(2),1:fm(3)) = tx(ib:ie+1,1:fm(2),1:fm(3))
! 9. calculate the even values of the perpendicular component using the
! divergence free condition;
!
do k = 1, fm(3)
do j = 1, fm(2)
jp1 = j + 1
do i = 2, fm(1), 2
im1 = i - 1
ip1 = i + 1
wx(i,j,k) = 0.50 * (wx(im1,j ,k) + wx(ip1,j ,k)) &
+ 0.50 * (wy(i ,jp1,k) - wy(im1,jp1,k)) &
- 0.50 * (wy(i ,j ,k) - wy(im1,j ,k))
end do
end do
end do
! deallocate temporary arrays
!
deallocate(tx)
deallocate(ty)
case default
end select
!
end subroutine expand_mag_bnd
#endif /* FLUXCT */
!
!===============================================================================