diff --git a/src/interpolations.F90 b/src/interpolations.F90 index ca62d08..2150319 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -71,15 +71,18 @@ module interpolations ! number of cells used in the Gaussian process reconstruction ! - integer , save :: ngp = 5 + integer , save :: ngp = 3 + integer , save :: mgp = 1 + integer , save :: dgp = 9 ! normal distribution width in the Gaussian process reconstruction ! - real(kind=8), save :: sgp = 1.0d+01 + real(kind=8), save :: sgp = 3.0d+00 ! Gaussian process reconstruction coefficients vector ! - real(kind=8), dimension(:) , allocatable, save :: cgp + real(kind=8), dimension(:) , allocatable, save :: cgp + real(kind=8), dimension(:,:,:), allocatable, save :: ugp ! flags for reconstruction corrections ! @@ -179,6 +182,11 @@ module interpolations ! kappa = min(kappa, (1.0d+00 - cfl) / cfl) +! calculate mgp +! + mgp = (ngp - 1) / 2 + dgp = ngp**NDIMS + ! select the reconstruction method ! select case(trim(sreconstruction)) @@ -288,6 +296,27 @@ module interpolations if (verbose .and. mod(ngp,2) == 0) & call print_warning("interpolations:initialize_interpolation" & , "The parameter ngp has to be integer with odd value.") + case ("mgp", "MGP") + write(stmp, '(f16.1)') sgp + write(name_rec, & + '("Multidimensional Gaussian Process (",i1,"-point, δ=",a,")")') & + ngp, trim(adjustl(stmp)) + +! allocate the Gaussian process reconstruction matrix and position vector +! + allocate(ugp(dgp,2,NDIMS)) + +! prepare matrix coefficients +! + call prepare_mgp() + + interfaces => interfaces_mgp + if (verbose .and. 2 * ng <= ngp - 1) & + call print_warning("interpolations:initialize_interpolation" & + , "Increase the number of ghost cells (at least (ngp+1)/2).") + if (verbose .and. mod(ngp,2) == 0) & + call print_warning("interpolations:initialize_interpolation" & + , "The parameter ngp has to be integer with odd value.") case default if (verbose) then write (*,"(1x,a)") "The selected reconstruction method is not " // & @@ -431,6 +460,7 @@ module interpolations ! deallocate Gaussian process reconstruction coefficient vector if used ! if (allocated(cgp)) deallocate(cgp) + if (allocated(ugp)) deallocate(ugp) ! release the procedure pointers ! @@ -451,6 +481,163 @@ module interpolations ! !=============================================================================== ! +! subroutine PREPARE_MGP: +! --------------------- +! +! Subroutine prepares matrixes for the multidimensional +! Gaussian Process (GP) method. +! +!=============================================================================== +! + subroutine prepare_mgp() + +! include external procedures +! + use algebra , only : invert + use constants , only : pi + use error , only : print_error + +! local variables are not implicit by default +! + implicit none + +! local variables +! + logical :: flag + integer :: i, j, i1, j1, k1, i2, j2, k2 + real(kind=16) :: sig, fc, fx, fy, fz, xl, xr, yl, yr, zl, zr + +! local arrays for derivatives +! + real(kind=16), dimension(:,:) , allocatable :: cov, inv + real(kind=16), dimension(:,:,:), allocatable :: xgp +! +!------------------------------------------------------------------------------- +! +! calculate normal distribution sigma +! + sig = sqrt(2.0d+00) * sgp + +! allocate the convariance matrix and interpolation position vector +! + allocate(cov(dgp,dgp)) + allocate(inv(dgp,dgp)) + allocate(xgp(dgp,2,NDIMS)) + +! prepare the covariance matrix +! + fc = 0.5d+00 * sqrt(pi) * sig + i = 0 +#if NDIMS == 3 + do k1 = - mgp, mgp +#endif /* NDIMS == 3 */ + do j1 = - mgp, mgp + do i1 = - mgp, mgp + i = i + 1 + j = 0 +#if NDIMS == 3 + do k2 = - mgp, mgp +#endif /* NDIMS == 3 */ + do j2 = - mgp, mgp + do i2 = - mgp, mgp + j = j + 1 + + xl = (1.0d+00 * (i1 - i2) - 0.5d+00) / sig + xr = (1.0d+00 * (i1 - i2) + 0.5d+00) / sig + yl = (1.0d+00 * (j1 - j2) - 0.5d+00) / sig + yr = (1.0d+00 * (j1 - j2) + 0.5d+00) / sig +#if NDIMS == 3 + zl = (1.0d+00 * (k1 - k2) - 0.5d+00) / sig + zr = (1.0d+00 * (k1 - k2) + 0.5d+00) / sig +#endif /* NDIMS == 3 */ + + cov(i,j) = fc * (erf(xr) - erf(xl)) * (erf(yr) - erf(yl)) +#if NDIMS == 3 + cov(i,j) = cov(i,j) * (erf(zr) - erf(zl)) +#endif /* NDIMS == 3 */ + end do + end do + end do + end do +#if NDIMS == 3 + end do + end do +#endif /* NDIMS == 3 */ + +! prepare the interpolation position vector +! + i = 0 +#if NDIMS == 3 + do k1 = - mgp, mgp +#endif /* NDIMS == 3 */ + do j1 = - mgp, mgp + do i1 = - mgp, mgp + i = i + 1 + + xl = (1.0d+00 * i1 - 0.5d+00) / sig + xr = (1.0d+00 * i1 + 0.5d+00) / sig + yl = (1.0d+00 * j1 - 0.5d+00) / sig + yr = (1.0d+00 * j1 + 0.5d+00) / sig +#if NDIMS == 3 + zl = (1.0d+00 * k1 - 0.5d+00) / sig + zr = (1.0d+00 * k1 + 0.5d+00) / sig +#endif /* NDIMS == 3 */ + + fx = erf(xr) - erf(xl) + fy = erf(yr) - erf(yl) +#if NDIMS == 3 + fz = erf(zr) - erf(zl) + + xgp(i,1,1) = exp(- xl**2) * fy * fz + xgp(i,2,1) = exp(- xr**2) * fy * fz + xgp(i,1,2) = exp(- yl**2) * fx * fz + xgp(i,2,2) = exp(- yr**2) * fx * fz + xgp(i,1,3) = exp(- zl**2) * fx * fy + xgp(i,2,3) = exp(- zr**2) * fx * fy +#else /* NDIMS == 3 */ + xgp(i,1,1) = exp(- xl**2) * fy + xgp(i,2,1) = exp(- xr**2) * fy + xgp(i,1,2) = exp(- yl**2) * fx + xgp(i,2,2) = exp(- yr**2) * fx +#endif /* NDIMS == 3 */ + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + +! invert the matrix +! + call invert(dgp, cov(1:dgp,1:dgp), inv(1:dgp,1:dgp), flag) + +! prepare the interpolation coefficients vector +! + do j = 1, NDIMS + do i = 1, 2 + ugp(1:dgp,i,j) = matmul(xgp(1:dgp,i,j), inv(1:dgp,1:dgp)) + end do + end do + +! deallocate the convariance matrix and interpolation position vector +! + deallocate(cov) + deallocate(inv) + deallocate(xgp) + +! check if the matrix was inverted successfully +! + if (.not. flag) then + call print_error("interpolations::prepare_mgp" & + , "Could not invert covariance matrix!") + stop + end if + +!------------------------------------------------------------------------------- +! + end subroutine prepare_mgp +! +!=============================================================================== +! ! subroutine INTERFACES_TVD: ! ------------------------- ! @@ -677,6 +864,176 @@ module interpolations ! !=============================================================================== ! +! subroutine INTERFACES_MGP: +! ------------------------- +! +! Subroutine reconstructs both side interfaces of variable using +! multidimentional Gaussian Process method. +! +! Arguments: +! +! positive - the variable positivity flag; +! h - the spatial step; +! q - the variable array; +! qi - the array of reconstructed interfaces (2 in each direction); +! +!=============================================================================== +! + subroutine interfaces_mgp(positive, h, q, qi) + +! include external procedures +! + use coordinates , only : im , jm , km + use coordinates , only : ib , jb , kb , ie , je , ke + use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + logical , intent(in) :: positive + real(kind=8), dimension(NDIMS) , intent(in) :: h + real(kind=8), dimension(im,jm,km) , intent(in) :: q + real(kind=8), dimension(im,jm,km,2,NDIMS), intent(out) :: qi + +! local variables +! + logical :: flag + integer :: i, il, iu, im1, ip1 + integer :: j, jl, ju, jm1, jp1 + integer :: k, kl, ku, km1, kp1 + +! local arrays for derivatives +! + real(kind=8), dimension(NDIMS) :: dql, dqr, dq + real(kind=8), dimension(dgp) :: u +! +!------------------------------------------------------------------------------- +! +! copy ghost zones +! + do k = 1, NDIMS + do j = 1, 2 + qi( 1:ib, 1:jm, 1:km,j,k) = q( 1:ib, 1:jm, 1:km) + qi(ie:im, 1:jm, 1:km,j,k) = q(ie:im, 1:jm, 1:km) + qi(ib:ie, 1:jb, 1:km,j,k) = q(ib:ie, 1:jb, 1:km) + qi(ib:ie,je:jm, 1:km,j,k) = q(ib:ie,je:jm, 1:km) +#if NDIMS == 3 + qi(ib:ie,jb:je, 1:kb,j,k) = q(ib:ie,jb:je, 1:kb) + qi(ib:ie,jb:je,ke:km,j,k) = q(ib:ie,jb:je,ke:km) +#endif /* NDIMS == 3 */ + end do + end do + +! interpolate interfaces using precomputed interpolation vectors +! + do k = kbl, keu +#if NDIMS == 3 + kl = k - mgp + ku = k + mgp + km1 = k - 1 + kp1 = k + 1 +#endif /* NDIMS == 3 */ + do j = jbl, jeu + jl = j - mgp + ju = j + mgp + jm1 = j - 1 + jp1 = j + 1 + do i = ibl, ieu + il = i - mgp + iu = i + mgp + im1 = i - 1 + ip1 = i + 1 + +#if NDIMS == 3 + u(:) = reshape(q(il:iu,jl:ju,kl:ku), (/ dgp /)) + + qi(i ,j,k,1,1) = sum(ugp(1:dgp,1,1) * u(1:dgp)) + qi(im1,j,k,2,1) = sum(ugp(1:dgp,2,1) * u(1:dgp)) + qi(i,j ,k,1,2) = sum(ugp(1:dgp,1,2) * u(1:dgp)) + qi(i,jm1,k,2,2) = sum(ugp(1:dgp,2,2) * u(1:dgp)) + qi(i,j,k ,1,3) = sum(ugp(1:dgp,1,3) * u(1:dgp)) + qi(i,j,km1,2,3) = sum(ugp(1:dgp,2,3) * u(1:dgp)) +#else /* NDIMS == 3 */ + u(:) = reshape(q(il:iu,jl:ju,k ), (/ dgp /)) + + qi(i ,j,k,1,1) = sum(ugp(1:dgp,1,1) * u(1:dgp)) + qi(im1,j,k,2,1) = sum(ugp(1:dgp,2,1) * u(1:dgp)) + qi(i,j ,k,1,2) = sum(ugp(1:dgp,1,2) * u(1:dgp)) + qi(i,jm1,k,2,2) = sum(ugp(1:dgp,2,2) * u(1:dgp)) +#endif /* NDIMS == 3 */ + +! if the interpolation is not monotonic, apply a TVD slope +! + flag = ((qi(i ,j,k,1,1) - q(ip1,j,k)) & + * (qi(i ,j,k,1,1) - q(i ,j,k)) > 0.0d+00) + flag = flag .or. ((qi(im1,j,k,2,1) - q(im1,j,k)) & + * (qi(im1,j,k,2,1) - q(i ,j,k)) > 0.0d+00) + flag = flag .or. ((qi(i,j ,k,1,2) - q(i,jp1,k)) & + * (qi(i,j ,k,1,2) - q(i,j ,k)) > 0.0d+00) + flag = flag .or. ((qi(i,jm1,k,2,2) - q(i,jm1,k)) & + * (qi(i,jm1,k,2,2) - q(i,j ,k)) > 0.0d+00) +#if NDIMS == 3 + flag = flag .or. ((qi(i,j,k ,1,3) - q(i,j,kp1)) & + * (qi(i,j,k ,1,3) - q(i,j,k )) > 0.0d+00) + flag = flag .or. ((qi(i,j,km1,2,3) - q(i,j,km1)) & + * (qi(i,j,km1,2,3) - q(i,j,k )) > 0.0d+00) +#endif /* NDIMS == 3 */ + + if (flag) then + +! calculate the TVD derivatives +! + dql(1) = q(i ,j,k) - q(im1,j,k) + dqr(1) = q(ip1,j,k) - q(i ,j,k) + dq (1) = limiter_tvd(0.5d+00, dql(1), dqr(1)) + + dql(2) = q(i,j ,k) - q(i,jm1,k) + dqr(2) = q(i,jp1,k) - q(i,j ,k) + dq (2) = limiter_tvd(0.5d+00, dql(2), dqr(2)) + +#if NDIMS == 3 + dql(3) = q(i,j,k ) - q(i,j,km1) + dqr(3) = q(i,j,kp1) - q(i,j,k ) + dq (3) = limiter_tvd(0.5d+00, dql(3), dqr(3)) +#endif /* NDIMS == 3 */ + +! limit the derivatives if they produce negative interpolation for positive +! variables +! + if (positive) then + do while (q(i,j,k) <= sum(abs(dq(1:NDIMS)))) + dq(:) = 0.5d+00 * dq(:) + end do + end if + +! interpolate states +! + qi(i ,j,k,1,1) = q(i,j,k) + dq(1) + qi(im1,j,k,2,1) = q(i,j,k) - dq(1) + + qi(i,j ,k,1,2) = q(i,j,k) + dq(2) + qi(i,jm1,k,2,2) = q(i,j,k) - dq(2) + +#if NDIMS == 3 + qi(i,j,k ,1,3) = q(i,j,k) + dq(3) + qi(i,j,km1,2,3) = q(i,j,k) - dq(3) +#endif /* NDIMS == 3 */ + + end if + + end do ! i = ibl, ieu + end do ! j = jbl, jeu + end do ! k = kbl, keu + +!------------------------------------------------------------------------------- +! + end subroutine interfaces_mgp +! +!=============================================================================== +! ! subroutine RECONSTRUCT: ! ---------------------- !