INTERPOLATIONS: Implement multidimentional GP reconstruction.
This is truely multidimentional reconstruction based on Gaussian Processes. All neighbors of the interpolation cell are taken into account and properly weighted. The size of the neighbor region can be controlled. The monotonicity of the interpolated interfaces is detected and if the case of lack it, the linear TVD reconstruction is applied. Therefore, only the monotonic regions are interpolated using the high-order Gauss Process. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
0a067f1f26
commit
670668c293
@ -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 :: 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:
|
||||
! ----------------------
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user