diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 5bc2245..67daed5 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -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 */ diff --git a/src/interpolation.F90 b/src/interpolation.F90 index a8ede8f..cd28db0 100644 --- a/src/interpolation.F90 +++ b/src/interpolation.F90 @@ -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 */ ! !===============================================================================