5397 lines
151 KiB
Fortran
5397 lines
151 KiB
Fortran
!!******************************************************************************
|
||
!!
|
||
!! This file is part of the AMUN source code, a program to perform
|
||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||
!! adaptive mesh.
|
||
!!
|
||
!! Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
|
||
!!
|
||
!! This program is free software: you can redistribute it and/or modify
|
||
!! it under the terms of the GNU General Public License as published by
|
||
!! the Free Software Foundation, either version 3 of the License, or
|
||
!! (at your option) any later version.
|
||
!!
|
||
!! This program is distributed in the hope that it will be useful,
|
||
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
!! GNU General Public License for more details.
|
||
!!
|
||
!! You should have received a copy of the GNU General Public License
|
||
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||
!!
|
||
!!******************************************************************************
|
||
!!
|
||
!! module: INTERPOLATIONS
|
||
!!
|
||
!! This module provides subroutines to interpolate variables and reconstruct
|
||
!! the Riemann states.
|
||
!!
|
||
!!
|
||
!!******************************************************************************
|
||
!
|
||
module interpolations
|
||
|
||
implicit none
|
||
|
||
! interfaces for procedure pointers
|
||
!
|
||
abstract interface
|
||
subroutine interfaces_iface(positive, h, q, qi)
|
||
logical , intent(in) :: positive
|
||
real(kind=8), dimension(:) , intent(in) :: h
|
||
real(kind=8), dimension(:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi
|
||
end subroutine
|
||
subroutine reconstruct_iface(h, fc, fl, fr)
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
end subroutine
|
||
function limiter_iface(x, a, b) result(c)
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
end function
|
||
end interface
|
||
|
||
! pointers to the reconstruction and limiter procedures
|
||
!
|
||
procedure(interfaces_iface) , pointer, save :: interfaces => null()
|
||
procedure(reconstruct_iface) , pointer, save :: reconstruct_states => null()
|
||
procedure(limiter_iface) , pointer, save :: limiter_tvd => null()
|
||
procedure(limiter_iface) , pointer, save :: limiter_prol => null()
|
||
procedure(limiter_iface) , pointer, save :: limiter_clip => null()
|
||
|
||
! method names
|
||
!
|
||
character(len=255), save :: name_rec = ""
|
||
character(len=255), save :: name_tlim = ""
|
||
character(len=255), save :: name_plim = ""
|
||
character(len=255), save :: name_clim = ""
|
||
|
||
! module parameters
|
||
!
|
||
real(kind=8), save :: eps = epsilon(1.0d+00)
|
||
real(kind=8), save :: rad = 0.5d+00
|
||
|
||
! monotonicity preserving reconstruction coefficients
|
||
!
|
||
real(kind=8), save :: kappa = 1.0d+00
|
||
real(kind=8), save :: kbeta = 1.0d+00
|
||
|
||
! number of ghost zones (required for compact schemes)
|
||
!
|
||
integer , save :: ng = 2
|
||
|
||
! order of the reconstruction
|
||
!
|
||
integer , save :: order = 1
|
||
|
||
! number of cells used in the Gaussian process reconstruction
|
||
!
|
||
integer , save :: ngp = 5
|
||
integer , save :: mgp = 1
|
||
integer , save :: dgp = 9
|
||
|
||
! normal distribution width in the Gaussian process reconstruction
|
||
!
|
||
real(kind=8), save :: sgp = 9.0d+00
|
||
|
||
! PPM constant
|
||
!
|
||
real(kind=8), save :: ppm_const = 1.25d+00
|
||
|
||
! weight toward central scheme for low-dissipation methods
|
||
!
|
||
real(kind=8), save :: cweight = 0.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
|
||
!
|
||
logical , save :: mlp = .false.
|
||
logical , save :: positivity = .true.
|
||
logical , save :: clip = .false.
|
||
|
||
! interpolation coefficients
|
||
!
|
||
real(kind=8), dimension(10), save :: ce9a
|
||
real(kind=8), dimension( 8), save :: ce7a
|
||
real(kind=8), dimension( 6), save :: ce5a
|
||
real(kind=8), dimension( 3), save :: di5a, di7a
|
||
real(kind=8), dimension( 5), save :: di9a
|
||
real(kind=8), dimension( 4), save :: ci5a
|
||
real(kind=8), dimension( 6), save :: ci7a, ci9a
|
||
|
||
real(kind=8), dimension(9), parameter :: &
|
||
ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
|
||
1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 ] / 2.520d+03
|
||
real(kind=8), dimension(7), parameter :: &
|
||
ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, &
|
||
2.140d+02,-3.800d+01, 4.000d+00 ] / 4.200d+02
|
||
real(kind=8), dimension(5), parameter :: &
|
||
ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01
|
||
real(kind=8), dimension(3), parameter :: &
|
||
ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
|
||
real(kind=8), dimension(2), parameter :: &
|
||
ce2 = [ 5.0d-01, 5.0d-01 ]
|
||
|
||
! by default everything is private
|
||
!
|
||
private
|
||
|
||
! declare public subroutines
|
||
!
|
||
public :: initialize_interpolations, finalize_interpolations
|
||
public :: print_interpolations
|
||
public :: interfaces, reconstruct, limiter_prol
|
||
public :: fix_positivity
|
||
public :: order
|
||
|
||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||
!
|
||
contains
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INITIALIZE_INTERPOLATIONS:
|
||
! ------------------------------------
|
||
!
|
||
! Subroutine initializes the interpolation module by reading the module
|
||
! parameters.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! nghosts - the number of ghosts cells at each side of a block;
|
||
! verbose - a logical flag turning the information printing;
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine initialize_interpolations(nghosts, verbose, status)
|
||
|
||
use helpers , only : print_message
|
||
use operators , only : initialize_operators
|
||
use parameters, only : get_parameter
|
||
|
||
implicit none
|
||
|
||
integer, intent(inout) :: nghosts
|
||
logical, intent(in) :: verbose
|
||
integer, intent(out) :: status
|
||
|
||
character(len=255) :: sreconstruction = "tvd"
|
||
character(len=255) :: tlimiter = "mm"
|
||
character(len=255) :: climiter = "mm"
|
||
character(len=255) :: plimiter = "avg"
|
||
character(len=255) :: mlp_limiting = "off"
|
||
character(len=255) :: positivity_fix = "on"
|
||
character(len=255) :: clip_extrema = "off"
|
||
character(len= 16) :: stmp
|
||
real(kind=8) :: cfl = 0.5d+00
|
||
|
||
character(len=*), parameter :: &
|
||
loc = 'INTERPOLATIONS:initialize_interpolations()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
! obtain the user defined interpolation methods and coefficients
|
||
!
|
||
call get_parameter("reconstruction" , sreconstruction)
|
||
call get_parameter("limiter" , tlimiter )
|
||
call get_parameter("fix_positivity" , positivity_fix )
|
||
call get_parameter("clip_extrema" , clip_extrema )
|
||
call get_parameter("extrema_limiter" , climiter )
|
||
call get_parameter("prolongation_limiter", plimiter )
|
||
call get_parameter("mlp_limiting" , mlp_limiting )
|
||
call get_parameter("ngp" , ngp )
|
||
call get_parameter("sgp" , sgp )
|
||
call get_parameter("eps" , eps )
|
||
call get_parameter("limo3_rad" , rad )
|
||
call get_parameter("kappa" , kappa )
|
||
call get_parameter("kbeta" , kbeta )
|
||
call get_parameter("ppm_const" , ppm_const )
|
||
call get_parameter("central_weight" , cweight )
|
||
call get_parameter("cfl" , cfl )
|
||
|
||
! calculate κ = (1 - ν) / ν
|
||
!
|
||
kappa = min(kappa, (1.0d+00 - cfl) / cfl)
|
||
|
||
! correct central weight
|
||
!
|
||
cweight = max(0.0d+00, min(1.0d+00, cweight))
|
||
|
||
! check ngp
|
||
!
|
||
if (mod(ngp,2) == 0 .or. ngp < 3) then
|
||
if (verbose) &
|
||
call print_message(loc, "The parameter ngp has to be an odd " // &
|
||
"integer >= 3. Resetting to default " // &
|
||
"value of ngp = 5.")
|
||
ngp = 5
|
||
end if
|
||
|
||
! calculate mgp
|
||
!
|
||
mgp = (ngp - 1) / 2
|
||
dgp = ngp**NDIMS
|
||
|
||
! select the reconstruction method
|
||
!
|
||
select case(trim(sreconstruction))
|
||
case ("tvd", "TVD")
|
||
name_rec = "2nd order TVD"
|
||
interfaces => interfaces_tvd
|
||
reconstruct_states => reconstruct_tvd
|
||
case ("weno3", "WENO3")
|
||
name_rec = "3rd order WENO"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_weno3
|
||
order = 3
|
||
case ("limo3", "LIMO3", "LimO3")
|
||
name_rec = "3rd order logarithmic limited"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_limo3
|
||
order = 3
|
||
eps = max(1.0d-12, eps)
|
||
case ("ppm", "PPM")
|
||
name_rec = "3rd order PPM"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_ppm
|
||
order = 3
|
||
nghosts = max(nghosts, 4)
|
||
case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z")
|
||
name_rec = "5th order WENO-Z (Borges et al. 2008)"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_weno5z
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("weno5yc", "weno5-yc", "WENO5YC", "WENO5-YC")
|
||
name_rec = "5th order WENO-YC (Yamaleev & Carpenter 2009)"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_weno5yc
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("weno5ns", "weno5-ns", "WENO5NS", "WENO5-NS")
|
||
name_rec = "5th order WENO-NS (Ha et al. 2013)"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_weno5ns
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("crweno5z", "crweno5-z", "CRWENO5Z", "CRWENO5-Z")
|
||
name_rec = "5th order Compact WENO-Z"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_crweno5z
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("crweno5yc", "crweno5-yc", "CRWENO5YC", "CRWENO5-YC")
|
||
name_rec = "5th order Compact WENO-YC"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_crweno5yc
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("crweno5ns", "crweno5-ns", "CRWENO5NS", "CRWENO5-NS")
|
||
name_rec = "5th order Compact WENO-NS"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_crweno5ns
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("mp5", "MP5")
|
||
name_rec = "5th order Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_mp5
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("mp7", "MP7")
|
||
name_rec = "7th order Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_mp7
|
||
order = 7
|
||
nghosts = max(nghosts, 4)
|
||
case ("mp9", "MP9")
|
||
name_rec = "9th order Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_mp9
|
||
order = 9
|
||
nghosts = max(nghosts, 6)
|
||
case ("crmp5", "CRMP5")
|
||
name_rec = "5th order Compact Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_crmp5
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("crmp7", "CRMP7")
|
||
name_rec = "7th order Compact Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_crmp7
|
||
order = 7
|
||
nghosts = max(nghosts, 4)
|
||
case ("crmp9", "CRMP9")
|
||
name_rec = "9th order Compact Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_crmp9
|
||
order = 9
|
||
nghosts = max(nghosts, 6)
|
||
case ("ocmp5", "OCMP5")
|
||
name_rec = "5th order Optimized Compact Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_ocmp5
|
||
order = 5
|
||
nghosts = max(nghosts, 4)
|
||
case ("ocmp7", "OCMP7")
|
||
name_rec = "7th order Optimized Compact Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_ocmp7
|
||
order = 7
|
||
nghosts = max(nghosts, 6)
|
||
case ("ocmp9", "OCMP9")
|
||
name_rec = "9th order Optimized Compact Monotonicity Preserving"
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_ocmp9
|
||
order = 9
|
||
nghosts = max(nghosts, 6)
|
||
case ("gp", "GP")
|
||
write(stmp, '(f16.1)') sgp
|
||
write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp &
|
||
, trim(adjustl(stmp))
|
||
|
||
! allocate the Gaussian process reconstruction matrix and position vector
|
||
!
|
||
allocate(cgp(ngp), stat = status)
|
||
|
||
! prepare matrix coefficients
|
||
!
|
||
if (status == 0) call prepare_gp(status)
|
||
|
||
interfaces => interfaces_dir
|
||
reconstruct_states => reconstruct_gp
|
||
order = ngp
|
||
|
||
ng = 2
|
||
do while(ng < (ngp + 1) / 2)
|
||
ng = ng + 2
|
||
end do
|
||
nghosts = max(nghosts, ng)
|
||
|
||
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), stat = status)
|
||
|
||
! prepare matrix coefficients
|
||
!
|
||
if (status == 0) call prepare_mgp(status)
|
||
|
||
interfaces => interfaces_mgp
|
||
|
||
ng = 2
|
||
do while(ng < (ngp + 1) / 2)
|
||
ng = ng + 2
|
||
end do
|
||
nghosts = max(nghosts, ng)
|
||
|
||
case default
|
||
|
||
if (verbose) &
|
||
call print_message(loc, "The selected reconstruction method is not "// &
|
||
"implemented: " // trim(sreconstruction) // "." // &
|
||
"Available methods: 'tvd' 'limo3', 'ppm'," // &
|
||
" 'weno3', 'weno5z', 'weno5yc', 'weno5ns'," // &
|
||
" 'crweno5z', 'crweno5yc', 'crweno5ns','mp5', " // &
|
||
" 'mp7', 'mp9', 'crmp5', 'crmp5ld', 'crmp7', " // &
|
||
"'gp', 'mgp'.")
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! select the TVD limiter
|
||
!
|
||
select case(trim(tlimiter))
|
||
case ("mm", "minmod")
|
||
name_tlim = "minmod"
|
||
limiter_tvd => limiter_minmod
|
||
order = max(2, order)
|
||
case ("mc", "monotonized_central")
|
||
name_tlim = "monotonized central"
|
||
limiter_tvd => limiter_monotonized_central
|
||
order = max(2, order)
|
||
case ("sb", "superbee")
|
||
name_tlim = "superbee"
|
||
limiter_tvd => limiter_superbee
|
||
order = max(2, order)
|
||
case ("vl", "vanleer")
|
||
name_tlim = "van Leer"
|
||
limiter_tvd => limiter_vanleer
|
||
order = max(2, order)
|
||
case ("va", "vanalbada")
|
||
name_tlim = "van Albada"
|
||
limiter_tvd => limiter_vanalbada
|
||
order = max(2, order)
|
||
case ("ospre")
|
||
name_tlim = "ospre"
|
||
limiter_tvd => limiter_ospre
|
||
order = max(2, order)
|
||
case default
|
||
name_tlim = "zero derivative"
|
||
limiter_tvd => limiter_zero
|
||
order = max(1, order)
|
||
end select
|
||
|
||
! select the prolongation limiter
|
||
!
|
||
select case(trim(plimiter))
|
||
case ("avg", "average", "mean")
|
||
name_plim = "average"
|
||
limiter_prol => limiter_average
|
||
case ("mm", "minmod")
|
||
name_plim = "minmod"
|
||
limiter_prol => limiter_minmod
|
||
case ("mc", "monotonized_central")
|
||
name_plim = "monotonized central"
|
||
limiter_prol => limiter_monotonized_central
|
||
case ("sb", "superbee")
|
||
name_plim = "superbee"
|
||
limiter_prol => limiter_superbee
|
||
case ("vl", "vanleer")
|
||
name_plim = "van Leer"
|
||
limiter_prol => limiter_vanleer
|
||
case ("va", "vanalbada")
|
||
name_plim = "van Albada"
|
||
limiter_prol => limiter_vanalbada
|
||
case ("ospre")
|
||
name_plim = "ospre"
|
||
limiter_prol => limiter_ospre
|
||
case default
|
||
name_plim = "zero derivative"
|
||
limiter_prol => limiter_zero
|
||
end select
|
||
|
||
! select the clipping limiter
|
||
!
|
||
select case(trim(climiter))
|
||
case ("mm", "minmod")
|
||
name_clim = "minmod"
|
||
limiter_clip => limiter_minmod
|
||
case ("mc", "monotonized_central")
|
||
name_clim = "monotonized central"
|
||
limiter_clip => limiter_monotonized_central
|
||
case ("sb", "superbee")
|
||
name_clim = "superbee"
|
||
limiter_clip => limiter_superbee
|
||
case ("vl", "vanleer")
|
||
name_clim = "van Leer"
|
||
limiter_clip => limiter_vanleer
|
||
case ("ospre")
|
||
name_clim = "ospre"
|
||
limiter_clip => limiter_ospre
|
||
case default
|
||
name_clim = "zero derivative"
|
||
limiter_clip => limiter_zero
|
||
end select
|
||
|
||
! check additional reconstruction limiting
|
||
!
|
||
select case(trim(mlp_limiting))
|
||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||
mlp = .true.
|
||
case default
|
||
mlp = .false.
|
||
end select
|
||
select case(trim(positivity_fix))
|
||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||
positivity = .true.
|
||
case default
|
||
positivity = .false.
|
||
end select
|
||
select case(trim(clip_extrema))
|
||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||
clip = .true.
|
||
case default
|
||
clip = .false.
|
||
end select
|
||
|
||
! set the internal number of ghost cells
|
||
!
|
||
ng = nghosts
|
||
|
||
! calculate coefficients of the low-dissipation methods
|
||
!
|
||
ce5a = [ - cweight +2.000d+00, 5.00d+00 * cweight -1.300d+01, &
|
||
-1.00d+01 * cweight +4.700d+01, 1.00d+01 * cweight +2.700d+01, &
|
||
-5.00d+00 * cweight -3.000d+00, cweight ] / 6.00d+01
|
||
ce7a = [ 3.00d+00 * cweight -6.000d+00, -2.10d+01 * cweight +5.000d+01, &
|
||
6.30d+01 * cweight -2.020d+02, -1.05d+02 * cweight +6.380d+02, &
|
||
1.05d+02 * cweight +4.280d+02, -6.30d+01 * cweight -7.600d+01, &
|
||
2.10d+01 * cweight +8.000d+00, -3.00d+00 * cweight ] / 8.40d+02
|
||
ce9a = [ -2.00d+00 * cweight +4.000d+00, 1.80d+01 * cweight -4.100d+01, &
|
||
-7.20d+01 * cweight +1.990d+02, 1.68d+02 * cweight -6.410d+02, &
|
||
-2.52d+02 * cweight +1.879d+03, 2.52d+02 * cweight +1.375d+03, &
|
||
-1.68d+02 * cweight -3.050d+02, 7.20d+01 * cweight +5.500d+01, &
|
||
-1.80d+01 * cweight -5.000d+00, 2.00d+00 * cweight ] / 2.52d+03
|
||
|
||
di5a = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00, &
|
||
(cweight +1.0d+00) / 6.0d+00 ]
|
||
ci5a = [ - cweight +2.0d+00, -9.0d+00 * cweight +3.8d+01, &
|
||
9.0d+00 * cweight +2.0d+01, cweight ] / 3.6d+01
|
||
|
||
di7a = [ -(cweight -4.0d+00) / 8.0d+00, 1.0d+00, &
|
||
(cweight +2.0d+00) / 8.0d+00 ]
|
||
ci7a = [ cweight -2.0d+00, -1.50d+01 * cweight +3.8d+01, &
|
||
-8.0d+01 * cweight +4.78d+02, 8.0d+01 * cweight +3.18d+02, &
|
||
1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
|
||
|
||
di9a = [ -2.0d+00 * cweight +5.0d+00, -1.0d+01 * cweight +4.0d+01, &
|
||
6.0d+01, 1.0d+01 * cweight +2.0d+01, &
|
||
2.0d+00 * cweight +1.0d+00 ] / 6.0d+01
|
||
ci9a = [ -3.00d+00 * cweight +6.000d+00,-1.75d+02 * cweight +4.810d+02, &
|
||
-3.00d+02 * cweight +1.881d+03, 3.00d+02 * cweight +1.281d+03, &
|
||
1.75d+02 * cweight +1.310d+02, 3.00d+00 * cweight ] / 1.8d+03
|
||
|
||
! initialize operators
|
||
!
|
||
call initialize_operators(order, status)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine initialize_interpolations
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FINALIZE_INTERPOLATIONS:
|
||
! ----------------------------------
|
||
!
|
||
! Subroutine finalizes the interpolation module by releasing all memory used
|
||
! by its module variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine finalize_interpolations(status)
|
||
|
||
use operators, only : finalize_operators
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(out) :: status
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
! finalize operators
|
||
!
|
||
call finalize_operators(status)
|
||
|
||
! deallocate Gaussian process reconstruction coefficient vector if used
|
||
!
|
||
if (allocated(cgp)) deallocate(cgp, stat = status)
|
||
if (allocated(ugp)) deallocate(ugp, stat = status)
|
||
|
||
! release the procedure pointers
|
||
!
|
||
nullify(reconstruct_states)
|
||
nullify(limiter_tvd)
|
||
nullify(limiter_prol)
|
||
nullify(limiter_clip)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine finalize_interpolations
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRINT_INTERPOLATIONS:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine prints module parameters and settings.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! verbose - a logical flag turning the information printing;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine print_interpolations(verbose)
|
||
|
||
use helpers, only : print_section, print_parameter
|
||
|
||
implicit none
|
||
|
||
logical, intent(in) :: verbose
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (verbose) then
|
||
|
||
call print_section(verbose, "Interpolations")
|
||
call print_parameter(verbose, "reconstruction" , name_rec )
|
||
call print_parameter(verbose, "TVD limiter" , name_tlim)
|
||
call print_parameter(verbose, "prolongation limiter", name_plim)
|
||
if (mlp) then
|
||
call print_parameter(verbose, "MLP limiting" , "on" )
|
||
else
|
||
call print_parameter(verbose, "MLP limiting" , "off" )
|
||
end if
|
||
if (positivity) then
|
||
call print_parameter(verbose, "fix positivity" , "on" )
|
||
else
|
||
call print_parameter(verbose, "fix positivity" , "off" )
|
||
end if
|
||
if (clip) then
|
||
call print_parameter(verbose, "clip extrema" , "on" )
|
||
call print_parameter(verbose, "extrema limiter" , name_clim)
|
||
else
|
||
call print_parameter(verbose, "clip extrema" , "off" )
|
||
end if
|
||
call print_parameter(verbose, "central weight", cweight)
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine print_interpolations
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PREPARE_MGP:
|
||
! ---------------------
|
||
!
|
||
! Subroutine prepares matrixes for the multidimensional
|
||
! Gaussian Process (GP) method.
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prepare_mgp(iret)
|
||
|
||
use algebra , only : max_prec, invert
|
||
use constants, only : pi
|
||
use helpers , only : print_message
|
||
|
||
implicit none
|
||
|
||
integer, intent(inout) :: iret
|
||
|
||
logical :: flag
|
||
integer :: i, j, i1, j1, i2, j2
|
||
#if NDIMS == 3
|
||
integer :: k1, k2
|
||
#endif /* NDIMS == 3 */
|
||
real(kind=max_prec) :: sig, fc, fx, fy, xl, xr, yl, yr
|
||
#if NDIMS == 3
|
||
real(kind=max_prec) :: fz, zl, zr
|
||
#endif /* NDIMS == 3 */
|
||
|
||
real(kind=max_prec), dimension(:,:) , allocatable :: cov, inv
|
||
real(kind=max_prec), dimension(:,:,:), allocatable :: xgp
|
||
|
||
character(len=*), parameter :: loc = 'INTERPOLATIONS::prepare_mgp()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! 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) = real(matmul(xgp(1:dgp,i,j), inv(1:dgp,1:dgp)), kind=8)
|
||
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_message(loc, "Could not invert the covariance matrix!")
|
||
iret = 201
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prepare_mgp
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INTERFACES_TVD:
|
||
! -------------------------
|
||
!
|
||
! Subroutine reconstructs both side interfaces of variable using TVD methods.
|
||
!
|
||
! 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_tvd(positive, h, q, qi)
|
||
|
||
use coordinates, only : nb, ne, nbl, neu
|
||
use helpers , only : print_message
|
||
|
||
implicit none
|
||
|
||
logical , intent(in) :: positive
|
||
real(kind=8), dimension(:) , intent(in) :: h
|
||
real(kind=8), dimension(:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi
|
||
|
||
integer :: i , j , k = 1
|
||
integer :: im1, jm1, ip1, jp1
|
||
#if NDIMS == 3
|
||
integer :: km1, kp1
|
||
#endif /* NDIMS == 3 */
|
||
|
||
character(len=80) :: msg
|
||
|
||
real(kind=8), dimension(NDIMS) :: dql, dqr, dq
|
||
|
||
character(len=*), parameter :: loc = 'INTERPOLATIONS::interfaces_tvd()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! copy ghost zones
|
||
!
|
||
do j = 1, NDIMS
|
||
do i = 1, 2
|
||
qi( :nb, : , : ,i,j) = q( :nb, : , : )
|
||
qi(ne: , : , : ,i,j) = q(ne: , : , : )
|
||
qi(nb:ne, :nb, : ,i,j) = q(nb:ne, :nb, : )
|
||
qi(nb:ne,ne: , : ,i,j) = q(nb:ne,ne: , : )
|
||
#if NDIMS == 3
|
||
qi(nb:ne,nb:ne, :nb,i,j) = q(nb:ne,nb:ne, :nb)
|
||
qi(nb:ne,nb:ne,ne: ,i,j) = q(nb:ne,nb:ne,ne: )
|
||
#endif /* NDIMS == 3 */
|
||
end do
|
||
end do
|
||
|
||
! interpolate interfaces
|
||
!
|
||
#if NDIMS == 3
|
||
do k = nbl, neu
|
||
km1 = k - 1
|
||
kp1 = k + 1
|
||
#endif /* NDIMS == 3 */
|
||
do j = nbl, neu
|
||
jm1 = j - 1
|
||
jp1 = j + 1
|
||
do i = nbl, neu
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! 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
|
||
if (q(i,j,k) > 0.0d+00) then
|
||
do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||
dq(:) = 0.5d+00 * dq(:)
|
||
end do
|
||
else
|
||
write(msg,"(a,3i0,a)") &
|
||
"Positive variable is not positive at ( ", i, j, k, " )"
|
||
call print_message(loc, msg)
|
||
dq(:) = 0.0d+00
|
||
end if
|
||
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 do ! i = nbl, neu
|
||
end do ! j = nbl, neu
|
||
#if NDIMS == 3
|
||
end do ! k = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! apply Multi-dimensional Limiting Process
|
||
!
|
||
if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine interfaces_tvd
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INTERFACES_DIR:
|
||
! -------------------------
|
||
!
|
||
! Subroutine reconstructs both side interfaces of variable separately
|
||
! along each direction.
|
||
!
|
||
! 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_dir(positive, h, q, qi)
|
||
|
||
use coordinates, only : nb, ne, nbl, neu
|
||
|
||
implicit none
|
||
|
||
logical , intent(in) :: positive
|
||
real(kind=8), dimension(:) , intent(in) :: h
|
||
real(kind=8), dimension(:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi
|
||
|
||
integer :: i, j, k = 1
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! copy ghost zones
|
||
!
|
||
do j = 1, NDIMS
|
||
do i = 1, 2
|
||
qi( :nb, : , : ,i,j) = q( :nb, : , : )
|
||
qi(ne: , : , : ,i,j) = q(ne: , : , : )
|
||
qi(nb:ne, :nb, : ,i,j) = q(nb:ne, :nb, : )
|
||
qi(nb:ne,ne: , : ,i,j) = q(nb:ne,ne: , : )
|
||
#if NDIMS == 3
|
||
qi(nb:ne,nb:ne, :nb,i,j) = q(nb:ne,nb:ne, :nb)
|
||
qi(nb:ne,nb:ne,ne: ,i,j) = q(nb:ne,nb:ne,ne: )
|
||
#endif /* NDIMS == 3 */
|
||
end do
|
||
end do
|
||
|
||
! interpolate interfaces
|
||
!
|
||
#if NDIMS == 3
|
||
do k = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
do j = nbl, neu
|
||
call reconstruct(h(1), q(:,j,k), qi(:,j,k,1,1), qi(:,j,k,2,1))
|
||
end do ! j = nbl, neu
|
||
do i = nbl, neu
|
||
call reconstruct(h(2), q(i,:,k), qi(i,:,k,1,2), qi(i,:,k,2,2))
|
||
end do ! i = nbl, neu
|
||
#if NDIMS == 3
|
||
end do ! k = nbl, neu
|
||
do j = nbl, neu
|
||
do i = nbl, neu
|
||
call reconstruct(h(3), q(i,j,:), qi(i,j,:,1,3), qi(i,j,:,2,3))
|
||
end do ! i = nbl, neu
|
||
end do ! j = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! make sure the interface states are positive for positive variables
|
||
!
|
||
if (positive) then
|
||
|
||
#if NDIMS == 3
|
||
do k = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
do j = nbl, neu
|
||
call fix_positivity(q(:,j,k), qi(:,j,k,1,1), qi(:,j,k,2,1))
|
||
end do ! j = nbl, neu
|
||
do i = nbl, neu
|
||
call fix_positivity(q(i,:,k), qi(i,:,k,1,2), qi(i,:,k,2,2))
|
||
end do ! i = nbl, neu
|
||
#if NDIMS == 3
|
||
end do ! k = nbl, neu
|
||
do j = nbl, neu
|
||
do i = nbl, neu
|
||
call fix_positivity(q(i,j,:), qi(i,j,:,1,3), qi(i,j,:,2,3))
|
||
end do ! i = nbl, neu
|
||
end do ! j = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
end if
|
||
|
||
! apply Multi-dimensional Limiting Process
|
||
!
|
||
if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine interfaces_dir
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INTERFACES_MGP:
|
||
! -------------------------
|
||
!
|
||
! Subroutine reconstructs both side interfaces of variable using
|
||
! multidimensional 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)
|
||
|
||
use coordinates, only : nb, ne, nbl, neu
|
||
use helpers , only : print_message
|
||
|
||
implicit none
|
||
|
||
logical , intent(in) :: positive
|
||
real(kind=8), dimension(:) , intent(in) :: h
|
||
real(kind=8), dimension(:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi
|
||
|
||
logical :: flag
|
||
integer :: i, il, iu, im1, ip1
|
||
integer :: j, jl, ju, jm1, jp1
|
||
integer :: k = 1
|
||
#if NDIMS == 3
|
||
integer :: kl, ku, km1, kp1
|
||
#endif /* NDIMS == 3 */
|
||
|
||
character(len=80) :: msg
|
||
|
||
real(kind=8), dimension(NDIMS) :: dql, dqr, dq
|
||
real(kind=8), dimension(dgp) :: u
|
||
|
||
character(len=*), parameter :: loc = 'INTERPOLATIONS::interfaces_mgp()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! copy ghost zones
|
||
!
|
||
do j = 1, NDIMS
|
||
do i = 1, 2
|
||
qi( :nb, : , : ,i,j) = q( :nb, : , : )
|
||
qi(ne: , : , : ,i,j) = q(ne: , : , : )
|
||
qi(nb:ne, :nb, : ,i,j) = q(nb:ne, :nb, : )
|
||
qi(nb:ne,ne: , : ,i,j) = q(nb:ne,ne: , : )
|
||
#if NDIMS == 3
|
||
qi(nb:ne,nb:ne, :nb,i,j) = q(nb:ne,nb:ne, :nb)
|
||
qi(nb:ne,nb:ne,ne: ,i,j) = q(nb:ne,nb:ne,ne: )
|
||
#endif /* NDIMS == 3 */
|
||
end do
|
||
end do
|
||
|
||
! interpolate interfaces using precomputed interpolation vectors
|
||
!
|
||
#if NDIMS == 3
|
||
do k = nbl, neu
|
||
kl = k - mgp
|
||
ku = k + mgp
|
||
km1 = k - 1
|
||
kp1 = k + 1
|
||
#endif /* NDIMS == 3 */
|
||
do j = nbl, neu
|
||
jl = j - mgp
|
||
ju = j + mgp
|
||
jm1 = j - 1
|
||
jp1 = j + 1
|
||
do i = nbl, neu
|
||
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 /)) - q(i,j,k)
|
||
|
||
qi(i ,j,k,1,1) = sum(ugp(1:dgp,1,1) * u(1:dgp)) + q(i,j,k)
|
||
qi(im1,j,k,2,1) = sum(ugp(1:dgp,2,1) * u(1:dgp)) + q(i,j,k)
|
||
qi(i,j ,k,1,2) = sum(ugp(1:dgp,1,2) * u(1:dgp)) + q(i,j,k)
|
||
qi(i,jm1,k,2,2) = sum(ugp(1:dgp,2,2) * u(1:dgp)) + q(i,j,k)
|
||
qi(i,j,k ,1,3) = sum(ugp(1:dgp,1,3) * u(1:dgp)) + q(i,j,k)
|
||
qi(i,j,km1,2,3) = sum(ugp(1:dgp,2,3) * u(1:dgp)) + q(i,j,k)
|
||
#else /* NDIMS == 3 */
|
||
u(:) = reshape(q(il:iu,jl:ju,k ), (/ dgp /)) - q(i,j,k)
|
||
|
||
qi(i ,j,k,1,1) = sum(ugp(1:dgp,1,1) * u(1:dgp)) + q(i,j,k)
|
||
qi(im1,j,k,2,1) = sum(ugp(1:dgp,2,1) * u(1:dgp)) + q(i,j,k)
|
||
qi(i,j ,k,1,2) = sum(ugp(1:dgp,1,2) * u(1:dgp)) + q(i,j,k)
|
||
qi(i,jm1,k,2,2) = sum(ugp(1:dgp,2,2) * u(1:dgp)) + q(i,j,k)
|
||
#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)) > eps)
|
||
flag = flag .or. ((qi(im1,j,k,2,1) - q(im1,j,k)) &
|
||
* (qi(im1,j,k,2,1) - q(i ,j,k)) > eps)
|
||
flag = flag .or. ((qi(i,j ,k,1,2) - q(i,jp1,k)) &
|
||
* (qi(i,j ,k,1,2) - q(i,j ,k)) > eps)
|
||
flag = flag .or. ((qi(i,jm1,k,2,2) - q(i,jm1,k)) &
|
||
* (qi(i,jm1,k,2,2) - q(i,j ,k)) > eps)
|
||
#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 )) > eps)
|
||
flag = flag .or. ((qi(i,j,km1,2,3) - q(i,j,km1)) &
|
||
* (qi(i,j,km1,2,3) - q(i,j,k )) > eps)
|
||
#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
|
||
if (q(i,j,k) > 0.0d+00) then
|
||
do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||
dq(:) = 0.5d+00 * dq(:)
|
||
end do
|
||
else
|
||
write(msg,"(a,3i0,a)") &
|
||
"Positive variable is not positive at ( ", i, j, k, " )"
|
||
call print_message(loc, msg)
|
||
dq(:) = 0.0d+00
|
||
end if
|
||
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 = nbl, neu
|
||
end do ! j = nbl, neu
|
||
#if NDIMS == 3
|
||
end do ! k = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! apply Multi-dimensional Limiting Process
|
||
!
|
||
if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine interfaces_mgp
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine MLP_LIMITING:
|
||
! -----------------------
|
||
!
|
||
! Subroutine applies the multi-dimensional limiting process to
|
||
! the reconstructed states in order to control oscillations due to
|
||
! the Gibbs phenomena near discontinuities.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the variable array;
|
||
! qi - the array of reconstructed interfaces (2 in each direction);
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Gerlinger, P.,
|
||
! "Multi-dimensional limiting for high-order schemes including
|
||
! turbulence and combustion",
|
||
! Journal of Computational Physics,
|
||
! 2012, vol. 231, pp. 2199-2228,
|
||
! http://dx.doi.org/10.1016/j.jcp.2011.10.024
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine mlp_limiting(q, qi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:), intent(inout) :: qi
|
||
|
||
integer :: i, im1, ip1
|
||
integer :: j, jm1, jp1
|
||
integer :: k = 1
|
||
#if NDIMS == 3
|
||
integer :: km1, kp1
|
||
#endif /* NDIMS == 3 */
|
||
integer :: m
|
||
#if NDIMS == 3
|
||
integer :: n, np1, np2
|
||
#endif /* NDIMS == 3 */
|
||
real(kind=8) :: qmn, qmx, dqc
|
||
real(kind=8) :: fc, fl
|
||
|
||
real(kind=8), dimension(NDIMS) :: dql, dqr, dq
|
||
real(kind=8), dimension(NDIMS) :: dqm, ap
|
||
#if NDIMS == 3
|
||
real(kind=8), dimension(NDIMS) :: hh, uu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
km1 = max( 1, k - 1)
|
||
kp1 = min(nn, k + 1)
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
jm1 = max( 1, j - 1)
|
||
jp1 = min(nn, j + 1)
|
||
do i = 1, nn
|
||
im1 = max( 1, i - 1)
|
||
ip1 = min(nn, i + 1)
|
||
|
||
! calculate the minmod 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_minmod(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_minmod(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_minmod(0.5d+00, dql(3), dqr(3))
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! calculate dqc
|
||
!
|
||
#if NDIMS == 3
|
||
qmn = minval(q(im1:ip1,jm1:jp1,km1:kp1))
|
||
qmx = maxval(q(im1:ip1,jm1:jp1,km1:kp1))
|
||
#else /* NDIMS == 3 */
|
||
qmn = minval(q(im1:ip1,jm1:jp1,k))
|
||
qmx = maxval(q(im1:ip1,jm1:jp1,k))
|
||
#endif /* NDIMS == 3 */
|
||
dqc = min(qmx - q(i,j,k), q(i,j,k) - qmn)
|
||
|
||
! if needed, apply the multi-dimensional limiting process
|
||
!
|
||
if (sum(abs(dq(1:NDIMS))) > dqc) then
|
||
|
||
dqm(1) = abs(q(ip1,j,k) - q(im1,j,k))
|
||
dqm(2) = abs(q(i,jp1,k) - q(i,jm1,k))
|
||
#if NDIMS == 3
|
||
dqm(3) = abs(q(i,j,kp1) - q(i,j,km1))
|
||
#endif /* NDIMS == 3 */
|
||
|
||
fc = dqc / max(1.0d-16, sum(abs(dqm(1:NDIMS))))
|
||
do m = 1, NDIMS
|
||
ap(m) = fc * abs(dqm(m))
|
||
end do
|
||
|
||
#if NDIMS == 3
|
||
do n = 1, NDIMS
|
||
hh(n) = max(ap(n) - abs(dq(n)), 0.0d+00)
|
||
np1 = mod(n , NDIMS) + 1
|
||
np2 = mod(n + 1, NDIMS) + 1
|
||
uu(n) = ap(n) - hh(n)
|
||
uu(np1) = ap(np1) + 0.5d+00 * hh(n)
|
||
uu(np2) = ap(np2) + 0.5d+00 * hh(n)
|
||
fc = hh(n) / (hh(n) + 1.0d-16)
|
||
fl = fc * (max(ap(np1) - abs(dq(np1)), 0.0d+00) &
|
||
- max(ap(np2) - abs(dq(np2)), 0.0d+00))
|
||
ap(n ) = uu(n)
|
||
ap(np1) = uu(np1) - fl
|
||
ap(np2) = uu(np2) + fl
|
||
end do
|
||
#else /* NDIMS == 3 */
|
||
fl = max(ap(1) - abs(dq(1)), 0.0d+00) &
|
||
- max(ap(2) - abs(dq(2)), 0.0d+00)
|
||
ap(1) = ap(1) - fl
|
||
ap(2) = ap(2) + fl
|
||
#endif /* NDIMS == 3 */
|
||
|
||
do m = 1, NDIMS
|
||
dq(m) = sign(ap(m), dq(m))
|
||
end do
|
||
|
||
end if
|
||
|
||
! update the interpolated states
|
||
!
|
||
dql(1) = qi(i ,j,k,1,1) - q(i,j,k)
|
||
dqr(1) = qi(im1,j,k,2,1) - q(i,j,k)
|
||
if (max(abs(dql(1)), abs(dqr(1))) > abs(dq(1))) then
|
||
qi(i ,j,k,1,1) = q(i,j,k) + dq(1)
|
||
qi(im1,j,k,2,1) = q(i,j,k) - dq(1)
|
||
end if
|
||
|
||
dql(2) = qi(i,j ,k,1,2) - q(i,j,k)
|
||
dqr(2) = qi(i,jm1,k,2,2) - q(i,j,k)
|
||
if (max(abs(dql(2)), abs(dqr(2))) > abs(dq(2))) then
|
||
qi(i,j ,k,1,2) = q(i,j,k) + dq(2)
|
||
qi(i,jm1,k,2,2) = q(i,j,k) - dq(2)
|
||
end if
|
||
|
||
#if NDIMS == 3
|
||
dql(3) = qi(i,j,k ,1,3) - q(i,j,k)
|
||
dqr(3) = qi(i,j,km1,2,3) - q(i,j,k)
|
||
if (max(abs(dql(3)), abs(dqr(3))) > abs(dq(3))) then
|
||
qi(i,j,k ,1,3) = q(i,j,k) + dq(3)
|
||
qi(i,j,km1,2,3) = q(i,j,k) - dq(3)
|
||
end if
|
||
#endif /* NDIMS == 3 */
|
||
|
||
end do ! i = 1, nn
|
||
end do ! j = 1, nn
|
||
#if NDIMS == 3
|
||
end do ! k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine mlp_limiting
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT:
|
||
! ----------------------
|
||
!
|
||
! Subroutine calls a reconstruction procedure, depending on the compilation
|
||
! flag SPACE, in order to interpolate the left and right states from their
|
||
! cell integrals. These states are required by any approximate Riemann
|
||
! solver.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! h - the spatial step; this is required for some reconstruction methods;
|
||
! f - the input vector of cell averaged values;
|
||
! fl - the left side state reconstructed for location (i+1/2);
|
||
! fr - the right side state reconstructed for location (i+1/2);
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct(h, f, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: f
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! reconstruct the states using the selected subroutine
|
||
!
|
||
call reconstruct_states(h, f(:), fl(:), fr(:))
|
||
|
||
! correct the reconstruction near extrema by clipping them in order to improve
|
||
! the stability of scheme
|
||
!
|
||
if (clip) call clip_extrema(f(:), fl(:), fr(:))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_TVD:
|
||
! --------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the second order TVD
|
||
! method with a selected limiter.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_tvd(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1
|
||
real(kind=8) :: df, dfl, dfr
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left- and right-side interface interpolations
|
||
!
|
||
do i = 2, n - 1
|
||
|
||
! calculate left and right indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate left and right side derivatives
|
||
!
|
||
dfl = fc(i ) - fc(im1)
|
||
dfr = fc(ip1) - fc(i )
|
||
|
||
! obtain the TVD limited derivative
|
||
!
|
||
df = limiter_tvd(0.5d+00, dfl, dfr)
|
||
|
||
! update the left and right-side interpolation states
|
||
!
|
||
fl(i ) = fc(i) + df
|
||
fr(im1) = fc(i) - df
|
||
|
||
end do ! i = 2, n - 1
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_tvd
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_WENO3:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the third order
|
||
! Weighted Essentially Non-Oscillatory (WENO) method.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Yamaleev & Carpenter, 2009, J. Comput. Phys., 228, 3025
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_weno3(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1
|
||
real(kind=8) :: bp, bm, ap, am, wp, wm, ww
|
||
real(kind=8) :: dfl, dfr, fp, fm, h2
|
||
|
||
real(kind=8), parameter :: dp = 2.0d+00 / 3.0d+00, dm = 1.0d+00 / 3.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! prepare common parameters
|
||
!
|
||
h2 = h * h
|
||
|
||
! iterate along the vector
|
||
!
|
||
do i = 2, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
dfl = fc(i ) - fc(im1)
|
||
dfr = fc(ip1) - fc(i )
|
||
|
||
! calculate coefficient omega
|
||
!
|
||
ww = (dfr - dfl)**2
|
||
|
||
! calculate corresponding betas
|
||
!
|
||
bp = dfr * dfr
|
||
bm = dfl * dfl
|
||
|
||
! calculate improved alphas
|
||
!
|
||
ap = 1.0d+00 + ww / (bp + h2)
|
||
am = 1.0d+00 + ww / (bm + h2)
|
||
|
||
! calculate weights
|
||
!
|
||
wp = dp * ap
|
||
wm = dm * am
|
||
ww = 2.0d+00 * (wp + wm)
|
||
|
||
! calculate central interpolation
|
||
!
|
||
fp = fc(i ) + fc(ip1)
|
||
|
||
! calculate left side interpolation
|
||
!
|
||
fm = - fc(im1) + 3.0d+00 * fc(i )
|
||
|
||
! calculate the left state
|
||
!
|
||
fl( i ) = (wp * fp + wm * fm) / ww
|
||
|
||
! calculate weights
|
||
!
|
||
wp = dp * am
|
||
wm = dm * ap
|
||
ww = 2.0d+00 * (wp + wm)
|
||
|
||
! calculate central interpolation
|
||
!
|
||
fp = fc(i ) + fc(im1)
|
||
|
||
! calculate right side interpolation
|
||
!
|
||
fm = - fc(ip1) + 3.0d+00 * fc(i )
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(im1) = (wp * fp + wm * fm) / ww
|
||
|
||
end do ! i = 2, n - 1
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_weno3
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_LIMO3:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the third order method
|
||
! with a limiter function LimO3.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Cada, M. & Torrilhon, M.,
|
||
! "Compact third-order limiter functions for finite volume methods",
|
||
! Journal of Computational Physics, 2009, 228, 4118-4145
|
||
! [2] Mignone, A., Tzeferacos, P., & Bodo, G.,
|
||
! "High-order conservative finite divergence GLM-MHD schemes for
|
||
! cell-centered MHD",
|
||
! Journal of Computational Physics, 2010, 229, 5896-5920
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_limo3(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1
|
||
real(kind=8) :: dfl, dfr
|
||
real(kind=8) :: th, et, f1, f2, xl, xi, rdx, rdx2
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! prepare parameters
|
||
!
|
||
rdx = rad * h
|
||
rdx2 = rdx * rdx
|
||
|
||
! iterate over positions and interpolate states
|
||
!
|
||
do i = 2, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! prepare left and right differences
|
||
!
|
||
dfl = fc(i ) - fc(im1)
|
||
dfr = fc(ip1) - fc(i )
|
||
|
||
! calculate the indicator function (eq. 3.17 in [1])
|
||
!
|
||
et = (dfl * dfl + dfr * dfr) / rdx2
|
||
|
||
! the switching function (embedded in eq. 3.22 in [1], eq. 32 in [2])
|
||
!
|
||
xi = max(0.0d+00, 0.5d+00 * min(2.0d+00, 1.0d+00 + (et - 1.0d+00) / eps))
|
||
xl = 1.0d+00 - xi
|
||
|
||
! calculate values at i + ½
|
||
!
|
||
if (abs(dfr) > eps) then
|
||
|
||
! calculate the slope ratio (eq. 2.8 in [1])
|
||
!
|
||
th = dfl / dfr
|
||
|
||
! calculate the quadratic reconstruction (eq. 3.8 in [1], divided by 2)
|
||
!
|
||
f1 = (2.0d+00 + th) / 6.0d+00
|
||
|
||
! calculate the third order limiter (eq. 3.13 in [1], cofficients divided by 2)
|
||
!
|
||
if (th >= 0.0d+00) then
|
||
f2 = max(0.0d+00, min(f1, th, 0.8d+00))
|
||
else
|
||
f2 = max(0.0d+00, min(f1, - 0.25d+00 * th))
|
||
end if
|
||
|
||
! interpolate the left state (eq. 3.5 in [1], eq. 30 in [2])
|
||
!
|
||
fl(i) = fc(i) + dfr * (xl * f1 + xi * f2)
|
||
|
||
else
|
||
|
||
fl(i) = fc(i)
|
||
|
||
end if
|
||
|
||
! calculate values at i - ½
|
||
!
|
||
if (abs(dfl) > eps) then
|
||
|
||
! calculate the slope ratio (eq. 2.8 in [1])
|
||
!
|
||
th = dfr / dfl
|
||
|
||
! calculate the quadratic reconstruction (eq. 3.8 in [1], divided by 2)
|
||
!
|
||
f1 = (2.0d+00 + th) / 6.0d+00
|
||
|
||
! calculate the third order limiter (eq. 3.13 in [1], cofficients divided by 2)
|
||
!
|
||
if (th >= 0.0d+00) then
|
||
f2 = max(0.0d+00, min(f1, th, 0.8d+00))
|
||
else
|
||
f2 = max(0.0d+00, min(f1, - 0.25d+00 * th))
|
||
end if
|
||
|
||
! interpolate the right state (eq. 3.5 in [1], eq. 30 in [2])
|
||
!
|
||
fr(im1) = fc(i) - dfl * (xl * f1 + xi * f2)
|
||
|
||
else
|
||
|
||
fr(im1) = fc(i)
|
||
|
||
end if
|
||
|
||
end do ! i = 2, n - 1
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_limo3
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_PPM:
|
||
! --------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the third order
|
||
! Piecewise Parabolic Method (PPM) with new limiters. This version lacks
|
||
! the support for flattening important for keeping the spurious oscillations
|
||
! near strong shocks/discontinuties under control.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Colella, P., Woodward, P. R.,
|
||
! "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
|
||
! Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
|
||
! https://dx.doi.org/10.1016/0021-9991(84)90143-8
|
||
! [2] Colella, P., Sekora, M. D.,
|
||
! "A limiter for PPM that preserves accuracy at smooth extrema",
|
||
! Journal of Computational Physics, 2008, vol. 227, pp. 7069-7076,
|
||
! https://doi.org/10.1016/j.jcp.2008.03.034
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_ppm(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
logical :: cm, cp, ext
|
||
integer :: n, i, im1, ip1, im2
|
||
real(kind=8) :: df2c, df2m, df2p, df2s, df2l
|
||
real(kind=8) :: dfm, dfp, dcm, dcp
|
||
real(kind=8) :: alm, alp, amx, sgn, del
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi, df2
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the high-order interface interpolation; eq. (16) in [2]
|
||
!
|
||
fi(1 ) = 5.0d-01 * (fc(1) + fc(2))
|
||
do i = 2, n - 2
|
||
fi(i) = (7.0d+00 * (fc(i) + fc(i+1)) - (fc(i-1) + fc(i+2))) / 1.2d+01
|
||
end do
|
||
fi(n-1) = 5.0d-01 * (fc(n-1) + fc(n))
|
||
fi(n ) = fc(n)
|
||
|
||
! calculate second-order central derivative
|
||
!
|
||
df2(1) = 0.0d+00
|
||
do i = 2, n - 1
|
||
df2(i) = (fc(i+1) + fc(i-1)) - 2.0d+00 * fc(i)
|
||
end do
|
||
df2(n) = 0.0d+00
|
||
|
||
! limit the interpolation; eqs. (18) and (19) in [2]
|
||
!
|
||
do i = 2, n - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
if ((fc(ip1) - fi(i)) * (fi(i) - fc(i)) < 0.0d+00) then
|
||
df2c = 3.0d+00 * ((fc(ip1) + fc(i )) - 2.0d+00 * fi(i))
|
||
df2m = df2(i )
|
||
df2p = df2(ip1)
|
||
if (min(df2c, df2m, df2p) * max(df2c, df2m, df2p) > 0.0d+00) then
|
||
df2l = sign(min(ppm_const * min(abs(df2m), abs(df2p)), abs(df2c)), df2c)
|
||
else
|
||
df2l = 0.0d+00
|
||
end if
|
||
fi(i) = 5.0d-01 * (fc(i) + fc(ip1)) - df2l / 6.0d+00
|
||
end if
|
||
end do
|
||
|
||
! iterate along the vector
|
||
!
|
||
do i = 2, n - 1
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! limit states if local extremum is detected or the interpolation is not
|
||
! monotone
|
||
!
|
||
alm = fi(im1) - fc(i)
|
||
alp = fi(i ) - fc(i)
|
||
cm = abs(alm) >= 2.0d+00 * abs(alp)
|
||
cp = abs(alp) >= 2.0d+00 * abs(alm)
|
||
ext = .false.
|
||
|
||
! check if we have local extremum here
|
||
!
|
||
if (alm * alp > 0.0d+00) then
|
||
ext = .true.
|
||
else if (cm .or. cp) then
|
||
im2 = max(1, i - 1)
|
||
|
||
dfm = fi(im1) - fi(im2)
|
||
dfp = fi(ip1) - fi(i )
|
||
dcm = fc(i ) - fc(im1)
|
||
dcp = fc(ip1) - fc(i )
|
||
if (min(abs(dfm),abs(dfp)) >= min(abs(dcm),abs(dcp))) then
|
||
ext = dfm * dfp < 0.0d+00
|
||
else
|
||
ext = dcm * dcp < 0.0d+00
|
||
end if
|
||
end if
|
||
|
||
! limit states if the local extremum is detected
|
||
!
|
||
if (ext) then
|
||
df2s = 6.0d+00 * (alm + alp)
|
||
df2m = df2(im1)
|
||
df2c = df2(i )
|
||
df2p = df2(ip1)
|
||
if (min(df2s, df2c, df2m, df2p) &
|
||
* max(df2s, df2c, df2m, df2p) > 0.0d+00) then
|
||
df2l = sign(min(ppm_const * min(abs(df2m), abs(df2c), abs(df2p)) &
|
||
, abs(df2s)), df2s)
|
||
if (abs(df2l) > 0.0d+00) then
|
||
alm = alm * df2l / df2s
|
||
alp = alp * df2l / df2s
|
||
else
|
||
alm = 0.0d+00
|
||
alp = 0.0d+00
|
||
end if
|
||
else
|
||
alm = 0.0d+00
|
||
alp = 0.0d+00
|
||
end if
|
||
else
|
||
|
||
! the interpolation is not monotonic so apply additional limiter
|
||
!
|
||
if (cp) then
|
||
sgn = sign(1.0d+00, alp)
|
||
amx = - 2.5d-01 * alp**2 / (alp + alm)
|
||
dfp = fc(ip1) - fc(i)
|
||
if (sgn * amx >= sgn * dfp) then
|
||
del = dfp * (dfp - alm)
|
||
if (del >= 0.0d+00) then
|
||
alp = - 2.0d+00 * (dfp + sgn * sqrt(del))
|
||
else
|
||
alp = - 2.0d+00 * alm
|
||
end if
|
||
end if
|
||
end if
|
||
|
||
if (cm) then
|
||
sgn = sign(1.0d+00, alm)
|
||
amx = - 2.5d-01 * alm**2 / (alp + alm)
|
||
dfm = fc(im1) - fc(i)
|
||
if (sgn * amx >= sgn * dfm) then
|
||
del = dfm * (dfm - alp)
|
||
if (del >= 0.0d+00) then
|
||
alm = - 2.0d+00 * (dfm + sgn * sqrt(del))
|
||
else
|
||
alm = - 2.0d+00 * alp
|
||
end if
|
||
end if
|
||
end if
|
||
end if
|
||
|
||
! update the states
|
||
!
|
||
fr(im1) = fc(i) + alm
|
||
fl(i ) = fc(i) + alp
|
||
|
||
end do ! i = 2, n
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
fl(1 ) = fi(1)
|
||
fr(n-1) = fi(n)
|
||
fl(n ) = fi(n)
|
||
fr(n ) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_ppm
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_WENO5Z:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the fifth order
|
||
! Explicit Weighted Essentially Non-Oscillatory (WENO5) method with
|
||
! stencil weights by Borges et al. (2008).
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Borges, R., Carmona, M., Costa, B., & Don, W.-S.,
|
||
! "An improved weighted essentially non-oscillatory scheme for
|
||
! hyperbolic conservation laws"
|
||
! Journal of Computational Physics,
|
||
! 2008, vol. 227, pp. 3191-3211,
|
||
! http://dx.doi.org/10.1016/j.jcp.2007.11.038
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_weno5z(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1, im2, ip2
|
||
real(kind=8) :: bl, bc, br, tt, df
|
||
real(kind=8) :: al, ac, ar
|
||
real(kind=8) :: wl, wc, wr, ww
|
||
real(kind=8) :: ql, qc, qr
|
||
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 &
|
||
, a12 = - 7.0d+00 / 6.0d+00 &
|
||
, a13 = 1.1d+01 / 6.0d+00
|
||
real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 &
|
||
, a22 = 5.0d+00 / 6.0d+00 &
|
||
, a23 = 2.0d+00 / 6.0d+00
|
||
real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 &
|
||
, a32 = 5.0d+00 / 6.0d+00 &
|
||
, a33 = - 1.0d+00 / 6.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
do i = 1, n - 1
|
||
ip1 = i + 1
|
||
dfp(i ) = fc(ip1) - fc(i)
|
||
dfm(ip1) = dfp(i)
|
||
end do
|
||
dfm(1) = dfp(1)
|
||
dfp(n) = dfm(n)
|
||
|
||
! calculate the absolute value of the second derivative
|
||
!
|
||
df2(:) = c1 * (dfp(:) - dfm(:))**2
|
||
|
||
! iterate along the vector
|
||
!
|
||
do i = 3, n - 2
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! calculate βₖ (eqs. 9-11 in [1])
|
||
!
|
||
bl = df2(im1) + c2 * (3.0d+00 * dfm(i ) - dfm(im1))**2
|
||
bc = df2(i ) + c2 * ( dfp(i ) + dfm(i ))**2
|
||
br = df2(ip1) + c2 * (3.0d+00 * dfp(i ) - dfp(ip1))**2
|
||
|
||
! calculate τ (below eq. 25 in [1])
|
||
!
|
||
tt = abs(br - bl)
|
||
|
||
! calculate αₖ (eq. 28 in [1])
|
||
!
|
||
al = 1.0d+00 + tt / (bl + eps)
|
||
ac = 1.0d+00 + tt / (bc + eps)
|
||
ar = 1.0d+00 + tt / (br + eps)
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al
|
||
wc = dc * ac
|
||
wr = dr * ar
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i ) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar
|
||
wc = dc * ac
|
||
wr = dr * al
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(im1) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
end do ! i = 3, n - 2
|
||
|
||
! update the interpolation of the first and last two points
|
||
!
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
df = limiter_tvd(0.5d+00, dfm(2), dfp(2))
|
||
fr(1) = fc(2) - df
|
||
fl(2) = fc(2) + df
|
||
i = n - 1
|
||
df = limiter_tvd(0.5d+00, dfm(i), dfp(i))
|
||
fr(i-1) = fc(i) - df
|
||
fl(i) = fc(i) + df
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_weno5z
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_WENO5YC:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the fifth order
|
||
! Explicit Weighted Essentially Non-Oscillatory (WENO5) method with
|
||
! stencil weights by Yamaleev & Carpenter (2009).
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Yamaleev, N. K. & Carpenter, H. C.,
|
||
! "A Systematic Methodology for Constructing High-Order Energy Stable
|
||
! WENO Schemes"
|
||
! Journal of Computational Physics,
|
||
! 2009, vol. 228, pp. 4248-4272,
|
||
! http://dx.doi.org/10.1016/j.jcp.2009.03.002
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_weno5yc(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1, im2, ip2
|
||
real(kind=8) :: bl, bc, br, tt, df
|
||
real(kind=8) :: al, ac, ar
|
||
real(kind=8) :: wl, wc, wr, ww
|
||
real(kind=8) :: ql, qc, qr
|
||
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 &
|
||
, a12 = - 7.0d+00 / 6.0d+00 &
|
||
, a13 = 1.1d+01 / 6.0d+00
|
||
real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 &
|
||
, a22 = 5.0d+00 / 6.0d+00 &
|
||
, a23 = 2.0d+00 / 6.0d+00
|
||
real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 &
|
||
, a32 = 5.0d+00 / 6.0d+00 &
|
||
, a33 = - 1.0d+00 / 6.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
do i = 1, n - 1
|
||
ip1 = i + 1
|
||
dfp(i ) = fc(ip1) - fc(i)
|
||
dfm(ip1) = dfp(i)
|
||
end do
|
||
dfm(1) = dfp(1)
|
||
dfp(n) = dfm(n)
|
||
|
||
! calculate the absolute value of the second derivative
|
||
!
|
||
df2(:) = c1 * (dfp(:) - dfm(:))**2
|
||
|
||
! iterate along the vector
|
||
!
|
||
do i = 3, n - 2
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! calculate βₖ (eq. 19 in [1])
|
||
!
|
||
bl = df2(im1) + c2 * (3.0d+00 * dfm(i ) - dfm(im1))**2
|
||
bc = df2(i ) + c2 * ( dfp(i ) + dfm(i ))**2
|
||
br = df2(ip1) + c2 * (3.0d+00 * dfp(i ) - dfp(ip1))**2
|
||
|
||
! calculate τ (below eq. 20 in [1])
|
||
!
|
||
tt = (6.0d+00 * fc(i) - 4.0d+00 * (fc(im1) + fc(ip1)) &
|
||
+ (fc(im2) + fc(ip2)))**2
|
||
|
||
! calculate αₖ (eqs. 18 or 58 in [1])
|
||
!
|
||
al = 1.0d+00 + tt / (bl + eps)
|
||
ac = 1.0d+00 + tt / (bc + eps)
|
||
ar = 1.0d+00 + tt / (br + eps)
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al
|
||
wc = dc * ac
|
||
wr = dr * ar
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i ) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar
|
||
wc = dc * ac
|
||
wr = dr * al
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(im1) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
end do ! i = 3, n - 2
|
||
|
||
! update the interpolation of the first and last two points
|
||
!
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
df = limiter_tvd(0.5d+00, dfm(2), dfp(2))
|
||
fr(1) = fc(2) - df
|
||
fl(2) = fc(2) + df
|
||
i = n - 1
|
||
df = limiter_tvd(0.5d+00, dfm(i), dfp(i))
|
||
fr(i-1) = fc(i) - df
|
||
fl(i) = fc(i) + df
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_weno5yc
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_WENO5NS:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the fifth order
|
||
! Explicit Weighted Essentially Non-Oscillatory (WENO5) method with new
|
||
! smoothness indicators and stencil weights by Ha et al. (2013).
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Ha, Y., Kim, C. H., Lee, Y. J., & Yoon, J.,
|
||
! "An improved weighted essentially non-oscillatory scheme with a new
|
||
! smoothness indicator",
|
||
! Journal of Computational Physics,
|
||
! 2013, vol. 232, pp. 68-86
|
||
! http://dx.doi.org/10.1016/j.jcp.2012.06.016
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_weno5ns(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1, im2, ip2
|
||
real(kind=8) :: bl, bc, br
|
||
real(kind=8) :: al, ac, ar, aa
|
||
real(kind=8) :: wl, wc, wr
|
||
real(kind=8) :: df, lq, l3, zt
|
||
real(kind=8) :: ql, qc, qr
|
||
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 &
|
||
, a12 = - 7.0d+00 / 6.0d+00 &
|
||
, a13 = 1.1d+01 / 6.0d+00
|
||
real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 &
|
||
, a22 = 5.0d+00 / 6.0d+00 &
|
||
, a23 = 2.0d+00 / 6.0d+00
|
||
real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 &
|
||
, a32 = 5.0d+00 / 6.0d+00 &
|
||
, a33 = - 1.0d+00 / 6.0d+00
|
||
|
||
real(kind=8), parameter :: xi = 4.0d-01
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
do i = 1, n - 1
|
||
ip1 = i + 1
|
||
dfp(i ) = fc(ip1) - fc(i)
|
||
dfm(ip1) = dfp(i)
|
||
end do
|
||
dfm(1) = dfp(1)
|
||
dfp(n) = dfm(n)
|
||
|
||
! calculate the absolute value of the second derivative
|
||
!
|
||
df2(:) = 0.5d+00 * abs(dfp(:) - dfm(:))
|
||
|
||
! iterate along the vector
|
||
!
|
||
do i = 3, n - 2
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! calculate βₖ (eq. 3.6 in [1])
|
||
!
|
||
df = abs(dfp(i))
|
||
lq = xi * df
|
||
bl = df2(im1) + xi * abs(2.0d+00 * dfm(i) - dfm(im1))
|
||
bc = df2(i ) + lq
|
||
br = df2(ip1) + lq
|
||
|
||
! calculate ζ (below eq. 3.6 in [1])
|
||
!
|
||
l3 = df**3
|
||
zt = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)
|
||
|
||
! calculate αₖ (eq. 3.9 in [4])
|
||
!
|
||
al = dl * (1.0d+00 + zt / (bl + eps)**2)
|
||
ac = dc * (1.0d+00 + zt / (bc + eps)**2)
|
||
ar = dr * (1.0d+00 + zt / (br + eps)**2)
|
||
|
||
! calculate weights
|
||
!
|
||
aa = (al + ar) + ac
|
||
wl = al / aa
|
||
wr = ar / aa
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i ) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! calculate βₖ (eq. 3.6 in [1])
|
||
!
|
||
df = abs(dfm(i))
|
||
lq = xi * df
|
||
bl = df2(ip1) + xi * abs(2.0d+00 * dfp(i) - dfp(ip1))
|
||
bc = df2(i ) + lq
|
||
br = df2(im1) + lq
|
||
|
||
! calculate ζ (below eq. 3.6 in [1])
|
||
|
||
l3 = df**3
|
||
zt = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)
|
||
|
||
! calculate αₖ (eq. 3.9 in [4])
|
||
!
|
||
al = dl * (1.0d+00 + zt / (bl + eps)**2)
|
||
ac = dc * (1.0d+00 + zt / (bc + eps)**2)
|
||
ar = dr * (1.0d+00 + zt / (br + eps)**2)
|
||
|
||
! normalize weights
|
||
!
|
||
aa = (al + ar) + ac
|
||
wl = al / aa
|
||
wr = ar / aa
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(im1) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
end do ! i = 3, n - 2
|
||
|
||
! update the interpolation of the first and last two points
|
||
!
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
df = limiter_tvd(0.5d+00, dfm(2), dfp(2))
|
||
fr(1) = fc(2) - df
|
||
fl(2) = fc(2) + df
|
||
i = n - 1
|
||
df = limiter_tvd(0.5d+00, dfm(i), dfp(i))
|
||
fr(i-1) = fc(i) - df
|
||
fl(i) = fc(i) + df
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_weno5ns
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_CRWENO5Z:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the fifth order
|
||
! Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
|
||
! method and smoothness indicators by Borges et al. (2008).
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Ghosh, D. & Baeder, J. D.,
|
||
! "Compact Reconstruction Schemes with Weighted ENO Limiting for
|
||
! Hyperbolic Conservation Laws"
|
||
! SIAM Journal on Scientific Computing,
|
||
! 2012, vol. 34, no. 3, pp. A1678-A1706,
|
||
! http://dx.doi.org/10.1137/110857659
|
||
! [2] Ghosh, D. & Baeder, J. D.,
|
||
! "Weighted Non-linear Compact Schemes for the Direct Numerical
|
||
! Simulation of Compressible, Turbulent Flows"
|
||
! Journal on Scientific Computing,
|
||
! 2014,
|
||
! http://dx.doi.org/10.1007/s10915-014-9818-0
|
||
! [3] Borges, R., Carmona, M., Costa, B., & Don, W.-S.,
|
||
! "An improved weighted essentially non-oscillatory scheme for
|
||
! hyperbolic conservation laws"
|
||
! Journal of Computational Physics,
|
||
! 2008, vol. 227, pp. 3191-3211,
|
||
! http://dx.doi.org/10.1016/j.jcp.2007.11.038
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_crweno5z(h, fc, fl, fr)
|
||
|
||
use algebra, only : tridiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1, im2, ip2
|
||
real(kind=8) :: bl, bc, br, tt
|
||
real(kind=8) :: wl, wc, wr, ww
|
||
real(kind=8) :: ql, qc, qr
|
||
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
real(kind=8), dimension(size(fc)) :: al, ac, ar
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
real(kind=8), dimension(size(fc),2) :: a, b, c, r
|
||
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
real(kind=8), parameter :: cl = 1.0d+00 / 9.0d+00
|
||
real(kind=8), parameter :: cc = 5.0d+00 / 9.0d+00
|
||
real(kind=8), parameter :: cr = 1.0d+00 / 3.0d+00
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
real(kind=8), parameter :: dq = 5.0d-01
|
||
|
||
real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 &
|
||
, a12 = - 7.0d+00 / 6.0d+00 &
|
||
, a13 = 1.1d+01 / 6.0d+00
|
||
real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 &
|
||
, a22 = 5.0d+00 / 6.0d+00 &
|
||
, a23 = 2.0d+00 / 6.0d+00
|
||
real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 &
|
||
, a32 = 5.0d+00 / 6.0d+00 &
|
||
, a33 = - 1.0d+00 / 6.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
do i = 1, n - 1
|
||
ip1 = i + 1
|
||
dfp(i ) = fc(ip1) - fc(i)
|
||
dfm(ip1) = dfp(i)
|
||
end do
|
||
dfm(1) = dfp(1)
|
||
dfp(n) = dfm(n)
|
||
|
||
! calculate the absolute value of the second derivative
|
||
!
|
||
df2(:) = c1 * (dfp(:) - dfm(:))**2
|
||
|
||
! prepare smoothness indicators
|
||
!
|
||
do i = 2, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate βₖ (eqs. 9-11 in [1])
|
||
!
|
||
bl = df2(im1) + c2 * (3.0d+00 * dfm(i ) - dfm(im1))**2
|
||
bc = df2(i ) + c2 * ( dfp(i ) + dfm(i ))**2
|
||
br = df2(ip1) + c2 * (3.0d+00 * dfp(i ) - dfp(ip1))**2
|
||
|
||
! calculate τ (below eq. 25 in [1])
|
||
!
|
||
tt = abs(br - bl)
|
||
|
||
! calculate αₖ (eq. 28 in [1])
|
||
!
|
||
al(i) = 1.0d+00 + tt / (bl + eps)
|
||
ac(i) = 1.0d+00 + tt / (bc + eps)
|
||
ar(i) = 1.0d+00 + tt / (br + eps)
|
||
|
||
end do ! i = 2, n - 1
|
||
|
||
! prepare tridiagonal system coefficients
|
||
!
|
||
do i = ng, n - ng + 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate weights
|
||
!
|
||
wl = cl * al(i)
|
||
wc = cc * ac(i)
|
||
wr = cr * ar(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate tridiagonal matrix coefficients
|
||
!
|
||
a(i,1) = 2.0d+00 * wl + wc
|
||
b(i,1) = wl + 2.0d+00 * (wc + wr)
|
||
c(i,1) = wr
|
||
|
||
! prepare right hand side of tridiagonal equation
|
||
!
|
||
r(i,1) = (wl * fc(im1) + (5.0d+00 * (wl + wc) + wr) * fc(i ) &
|
||
+ (wc + 5.0d+00 * wr) * fc(ip1)) * dq
|
||
|
||
! calculate weights
|
||
!
|
||
wl = cl * ar(i)
|
||
wc = cc * ac(i)
|
||
wr = cr * al(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate tridiagonal matrix coefficients
|
||
!
|
||
a(i,2) = wr
|
||
b(i,2) = wl + 2.0d+00 * (wc + wr)
|
||
c(i,2) = 2.0d+00 * wl + wc
|
||
|
||
! prepare right hand side of tridiagonal equation
|
||
!
|
||
r(i,2) = (wl * fc(ip1) + (5.0d+00 * (wl + wc) + wr) * fc(i ) &
|
||
+ (wc + 5.0d+00 * wr) * fc(im1)) * dq
|
||
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
! interpolate ghost zones using explicit solver (left-side reconstruction)
|
||
!
|
||
do i = 2, ng
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * ar(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,1) = 0.0d+00
|
||
b(i,1) = 1.0d+00
|
||
c(i,1) = 0.0d+00
|
||
r(i,1) = fl(i)
|
||
|
||
end do ! i = 2, ng
|
||
a(1,1) = 0.0d+00
|
||
b(1,1) = 1.0d+00
|
||
c(1,1) = 0.0d+00
|
||
r(1,1) = 0.5d+00 * (fc(1) + fc(2))
|
||
|
||
! interpolate ghost zones using explicit solver (left-side reconstruction)
|
||
!
|
||
do i = n - ng, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = min(n, i + 2)
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * ar(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,1) = 0.0d+00
|
||
b(i,1) = 1.0d+00
|
||
c(i,1) = 0.0d+00
|
||
r(i,1) = fl(i)
|
||
|
||
end do ! i = n - ng, n - 1
|
||
a(n,1) = 0.0d+00
|
||
b(n,1) = 1.0d+00
|
||
c(n,1) = 0.0d+00
|
||
r(n,1) = fc(n)
|
||
|
||
! interpolate ghost zones using explicit solver (right-side reconstruction)
|
||
!
|
||
do i = 2, ng + 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * al(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,2) = 0.0d+00
|
||
b(i,2) = 1.0d+00
|
||
c(i,2) = 0.0d+00
|
||
r(i,2) = fr(i)
|
||
|
||
end do ! i = 2, ng + 1
|
||
a(1,2) = 0.0d+00
|
||
b(1,2) = 1.0d+00
|
||
c(1,2) = 0.0d+00
|
||
r(1,2) = fc(1)
|
||
|
||
! interpolate ghost zones using explicit solver (right-side reconstruction)
|
||
!
|
||
do i = n - ng + 1, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = max(1, i - 1)
|
||
ip1 = min(n, i + 1)
|
||
ip2 = min(n, i + 2)
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * al(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,2) = 0.0d+00
|
||
b(i,2) = 1.0d+00
|
||
c(i,2) = 0.0d+00
|
||
r(i,2) = fr(i)
|
||
|
||
end do ! i = n - ng + 1, n - 1
|
||
a(n,2) = 0.0d+00
|
||
b(n,2) = 1.0d+00
|
||
c(n,2) = 0.0d+00
|
||
r(n,2) = 0.5d+00 * (fc(n-1) + fc(n))
|
||
|
||
! solve the tridiagonal system of equations for the left-side interpolation
|
||
!
|
||
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
|
||
|
||
! substitute the left-side values
|
||
!
|
||
fl(1:n ) = u(1:n)
|
||
|
||
! solve the tridiagonal system of equations for the left-side interpolation
|
||
!
|
||
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
|
||
|
||
! substitute the right-side values
|
||
!
|
||
fr(1:n-1) = u(2:n)
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_crweno5z
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_CRWENO5YC:
|
||
! --------------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the fifth order
|
||
! Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
|
||
! method and smoothness indicators by Yamaleev & Carpenter (2009).
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Ghosh, D. & Baeder, J. D.,
|
||
! "Compact Reconstruction Schemes with Weighted ENO Limiting for
|
||
! Hyperbolic Conservation Laws"
|
||
! SIAM Journal on Scientific Computing,
|
||
! 2012, vol. 34, no. 3, pp. A1678-A1706,
|
||
! http://dx.doi.org/10.1137/110857659
|
||
! [2] Ghosh, D. & Baeder, J. D.,
|
||
! "Weighted Non-linear Compact Schemes for the Direct Numerical
|
||
! Simulation of Compressible, Turbulent Flows"
|
||
! Journal on Scientific Computing,
|
||
! 2014,
|
||
! http://dx.doi.org/10.1007/s10915-014-9818-0
|
||
! [3] Yamaleev, N. K. & Carpenter, H. C.,
|
||
! "A Systematic Methodology for Constructing High-Order Energy Stable
|
||
! WENO Schemes"
|
||
! Journal of Computational Physics,
|
||
! 2009, vol. 228, pp. 4248-4272,
|
||
! http://dx.doi.org/10.1016/j.jcp.2009.03.002
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_crweno5yc(h, fc, fl, fr)
|
||
|
||
use algebra, only : tridiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1, im2, ip2
|
||
real(kind=8) :: bl, bc, br, tt
|
||
real(kind=8) :: wl, wc, wr, ww
|
||
real(kind=8) :: ql, qc, qr
|
||
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
real(kind=8), dimension(size(fc)) :: al, ac, ar
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
real(kind=8), dimension(size(fc),2) :: a, b, c, r
|
||
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
real(kind=8), parameter :: cl = 1.0d+00 / 9.0d+00
|
||
real(kind=8), parameter :: cc = 5.0d+00 / 9.0d+00
|
||
real(kind=8), parameter :: cr = 1.0d+00 / 3.0d+00
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
real(kind=8), parameter :: dq = 5.0d-01
|
||
|
||
real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 &
|
||
, a12 = - 7.0d+00 / 6.0d+00 &
|
||
, a13 = 1.1d+01 / 6.0d+00
|
||
real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 &
|
||
, a22 = 5.0d+00 / 6.0d+00 &
|
||
, a23 = 2.0d+00 / 6.0d+00
|
||
real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 &
|
||
, a32 = 5.0d+00 / 6.0d+00 &
|
||
, a33 = - 1.0d+00 / 6.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
do i = 1, n - 1
|
||
ip1 = i + 1
|
||
dfp(i ) = fc(ip1) - fc(i)
|
||
dfm(ip1) = dfp(i)
|
||
end do
|
||
dfm(1) = dfp(1)
|
||
dfp(n) = dfm(n)
|
||
|
||
! calculate the absolute value of the second derivative
|
||
!
|
||
df2(:) = c1 * (dfp(:) - dfm(:))**2
|
||
|
||
! prepare smoothness indicators
|
||
!
|
||
do i = 2, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = min(n, i + 2)
|
||
|
||
! calculate βₖ (eqs. 9-11 in [1])
|
||
!
|
||
bl = df2(im1) + c2 * (3.0d+00 * dfm(i ) - dfm(im1))**2
|
||
bc = df2(i ) + c2 * ( dfp(i ) + dfm(i ))**2
|
||
br = df2(ip1) + c2 * (3.0d+00 * dfp(i ) - dfp(ip1))**2
|
||
|
||
! calculate τ (below eq. 64 in [3])
|
||
!
|
||
tt = (6.0d+00 * fc(i) + (fc(im2) + fc(ip2)) &
|
||
- 4.0d+00 * (fc(im1) + fc(ip1)))**2
|
||
|
||
! calculate αₖ (eq. 28 in [1])
|
||
!
|
||
al(i) = 1.0d+00 + tt / (bl + eps)
|
||
ac(i) = 1.0d+00 + tt / (bc + eps)
|
||
ar(i) = 1.0d+00 + tt / (br + eps)
|
||
|
||
end do ! i = 2, n - 1
|
||
|
||
! prepare tridiagonal system coefficients
|
||
!
|
||
do i = ng, n - ng + 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate weights
|
||
!
|
||
wl = cl * al(i)
|
||
wc = cc * ac(i)
|
||
wr = cr * ar(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate tridiagonal matrix coefficients
|
||
!
|
||
a(i,1) = 2.0d+00 * wl + wc
|
||
b(i,1) = wl + 2.0d+00 * (wc + wr)
|
||
c(i,1) = wr
|
||
|
||
! prepare right hand side of tridiagonal equation
|
||
!
|
||
r(i,1) = (wl * fc(im1) + (5.0d+00 * (wl + wc) + wr) * fc(i ) &
|
||
+ (wc + 5.0d+00 * wr) * fc(ip1)) * dq
|
||
|
||
! calculate weights
|
||
!
|
||
wl = cl * ar(i)
|
||
wc = cc * ac(i)
|
||
wr = cr * al(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate tridiagonal matrix coefficients
|
||
!
|
||
a(i,2) = wr
|
||
b(i,2) = wl + 2.0d+00 * (wc + wr)
|
||
c(i,2) = 2.0d+00 * wl + wc
|
||
|
||
! prepare right hand side of tridiagonal equation
|
||
!
|
||
r(i,2) = (wl * fc(ip1) + (5.0d+00 * (wl + wc) + wr) * fc(i ) &
|
||
+ (wc + 5.0d+00 * wr) * fc(im1)) * dq
|
||
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
! interpolate ghost zones using explicit solver (left-side reconstruction)
|
||
!
|
||
do i = 2, ng
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * ar(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,1) = 0.0d+00
|
||
b(i,1) = 1.0d+00
|
||
c(i,1) = 0.0d+00
|
||
r(i,1) = fl(i)
|
||
|
||
end do ! i = 2, ng
|
||
a(1,1) = 0.0d+00
|
||
b(1,1) = 1.0d+00
|
||
c(1,1) = 0.0d+00
|
||
r(1,1) = 0.5d+00 * (fc(1) + fc(2))
|
||
|
||
! interpolate ghost zones using explicit solver (left-side reconstruction)
|
||
!
|
||
do i = n - ng, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = min(n, i + 2)
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * ar(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,1) = 0.0d+00
|
||
b(i,1) = 1.0d+00
|
||
c(i,1) = 0.0d+00
|
||
r(i,1) = fl(i)
|
||
|
||
end do ! i = n - ng, n - 1
|
||
a(n,1) = 0.0d+00
|
||
b(n,1) = 1.0d+00
|
||
c(n,1) = 0.0d+00
|
||
r(n,1) = fc(n)
|
||
|
||
! interpolate ghost zones using explicit solver (right-side reconstruction)
|
||
!
|
||
do i = 2, ng + 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * al(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,2) = 0.0d+00
|
||
b(i,2) = 1.0d+00
|
||
c(i,2) = 0.0d+00
|
||
r(i,2) = fr(i)
|
||
|
||
end do ! i = 2, ng + 1
|
||
a(1,2) = 0.0d+00
|
||
b(1,2) = 1.0d+00
|
||
c(1,2) = 0.0d+00
|
||
r(1,2) = fc(1)
|
||
|
||
! interpolate ghost zones using explicit solver (right-side reconstruction)
|
||
!
|
||
do i = n - ng + 1, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = min(n, i + 2)
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar(i)
|
||
wc = dc * ac(i)
|
||
wr = dr * al(i)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,2) = 0.0d+00
|
||
b(i,2) = 1.0d+00
|
||
c(i,2) = 0.0d+00
|
||
r(i,2) = fr(i)
|
||
|
||
end do ! i = n - ng + 1, n - 1
|
||
a(n,2) = 0.0d+00
|
||
b(n,2) = 1.0d+00
|
||
c(n,2) = 0.0d+00
|
||
r(n,2) = 0.5d+00 * (fc(n-1) + fc(n))
|
||
|
||
! solve the tridiagonal system of equations for the left-side interpolation
|
||
!
|
||
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
|
||
|
||
! substitute the left-side values
|
||
!
|
||
fl(1:n ) = u(1:n)
|
||
|
||
! solve the tridiagonal system of equations for the left-side interpolation
|
||
!
|
||
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
|
||
|
||
! substitute the right-side values
|
||
!
|
||
fr(1:n-1) = u(2:n)
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_crweno5yc
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_CRWENO5NS:
|
||
! --------------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the fifth order
|
||
! Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
|
||
! method combined with the smoothness indicators by Ha et al. (2013).
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Ghosh, D. & Baeder, J. D.,
|
||
! "Compact Reconstruction Schemes with Weighted ENO Limiting for
|
||
! Hyperbolic Conservation Laws"
|
||
! SIAM Journal on Scientific Computing,
|
||
! 2012, vol. 34, no. 3, pp. A1678-A1706,
|
||
! http://dx.doi.org/10.1137/110857659
|
||
! [2] Ghosh, D. & Baeder, J. D.,
|
||
! "Weighted Non-linear Compact Schemes for the Direct Numerical
|
||
! Simulation of Compressible, Turbulent Flows"
|
||
! Journal on Scientific Computing,
|
||
! 2014,
|
||
! http://dx.doi.org/10.1007/s10915-014-9818-0
|
||
! [3] Ha, Y., Kim, C. H., Lee, Y. J., & Yoon, J.,
|
||
! "An improved weighted essentially non-oscillatory scheme with a new
|
||
! smoothness indicator",
|
||
! Journal of Computational Physics,
|
||
! 2013, vol. 232, pp. 68-86
|
||
! http://dx.doi.org/10.1016/j.jcp.2012.06.016
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_crweno5ns(h, fc, fl, fr)
|
||
|
||
use algebra, only : tridiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i, im1, ip1, im2, ip2
|
||
real(kind=8) :: bl, bc, br
|
||
real(kind=8) :: wl, wc, wr, ww
|
||
real(kind=8) :: df, lq, l3, zt
|
||
real(kind=8) :: ql, qc, qr
|
||
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
real(kind=8), dimension(size(fc),2) :: al, ac, ar
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
real(kind=8), dimension(size(fc),2) :: a, b, c, r
|
||
|
||
real(kind=8), parameter :: xi = 4.0d-01
|
||
|
||
real(kind=8), parameter :: cl = 1.0d+00 / 9.0d+00
|
||
real(kind=8), parameter :: cc = 5.0d+00 / 9.0d+00
|
||
real(kind=8), parameter :: cr = 1.0d+00 / 3.0d+00
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
real(kind=8), parameter :: dq = 5.0d-01
|
||
|
||
real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 &
|
||
, a12 = - 7.0d+00 / 6.0d+00 &
|
||
, a13 = 1.1d+01 / 6.0d+00
|
||
real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 &
|
||
, a22 = 5.0d+00 / 6.0d+00 &
|
||
, a23 = 2.0d+00 / 6.0d+00
|
||
real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 &
|
||
, a32 = 5.0d+00 / 6.0d+00 &
|
||
, a33 = - 1.0d+00 / 6.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
do i = 1, n - 1
|
||
ip1 = i + 1
|
||
dfp(i ) = fc(ip1) - fc(i)
|
||
dfm(ip1) = dfp(i)
|
||
end do
|
||
dfm(1) = dfp(1)
|
||
dfp(n) = dfm(n)
|
||
|
||
! calculate the absolute value of the second derivative
|
||
!
|
||
df2(:) = 0.5d+00 * abs(dfp(:) - dfm(:))
|
||
|
||
! prepare smoothness indicators
|
||
!
|
||
do i = 2, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate βₖ
|
||
!
|
||
df = abs(dfp(i))
|
||
lq = xi * df
|
||
bl = df2(im1) + xi * abs(2.0d+00 * dfm(i) - dfm(im1))
|
||
bc = df2(i ) + lq
|
||
br = df2(ip1) + lq
|
||
|
||
! calculate ζ
|
||
!
|
||
l3 = df**3
|
||
zt = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)
|
||
|
||
! calculate αₖ
|
||
!
|
||
al(i,1) = 1.0d+00 + zt / (bl + eps)**2
|
||
ac(i,1) = 1.0d+00 + zt / (bc + eps)**2
|
||
ar(i,1) = 1.0d+00 + zt / (br + eps)**2
|
||
|
||
! calculate βₖ
|
||
!
|
||
df = abs(dfm(i))
|
||
lq = xi * df
|
||
bl = df2(im1) + lq
|
||
bc = df2(i ) + lq
|
||
br = df2(ip1) + xi * abs(2.0d+00 * dfp(i) - dfp(ip1))
|
||
|
||
! calculate ζ
|
||
|
||
l3 = df**3
|
||
zt = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)
|
||
|
||
! calculate αₖ
|
||
!
|
||
al(i,2) = 1.0d+00 + zt / (bl + eps)**2
|
||
ac(i,2) = 1.0d+00 + zt / (bc + eps)**2
|
||
ar(i,2) = 1.0d+00 + zt / (br + eps)**2
|
||
|
||
end do ! i = 2, n - 1
|
||
|
||
! prepare tridiagonal system coefficients
|
||
!
|
||
do i = ng, n - ng + 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
|
||
! calculate weights
|
||
!
|
||
wl = cl * al(i,1)
|
||
wc = cc * ac(i,1)
|
||
wr = cr * ar(i,1)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate tridiagonal matrix coefficients
|
||
!
|
||
a(i,1) = 2.0d+00 * wl + wc
|
||
b(i,1) = wl + 2.0d+00 * (wc + wr)
|
||
c(i,1) = wr
|
||
|
||
! prepare right hand side of tridiagonal equation
|
||
!
|
||
r(i,1) = (wl * fc(im1) + (5.0d+00 * (wl + wc) + wr) * fc(i ) &
|
||
+ (wc + 5.0d+00 * wr) * fc(ip1)) * dq
|
||
|
||
! calculate weights
|
||
!
|
||
wl = cl * ar(i,2)
|
||
wc = cc * ac(i,2)
|
||
wr = cr * al(i,2)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate tridiagonal matrix coefficients
|
||
!
|
||
a(i,2) = wr
|
||
b(i,2) = wl + 2.0d+00 * (wc + wr)
|
||
c(i,2) = 2.0d+00 * wl + wc
|
||
|
||
! prepare right hand side of tridiagonal equation
|
||
!
|
||
r(i,2) = (wl * fc(ip1) + (5.0d+00 * (wl + wc) + wr) * fc(i ) &
|
||
+ (wc + 5.0d+00 * wr) * fc(im1)) * dq
|
||
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
! interpolate ghost zones using explicit solver (left-side reconstruction)
|
||
!
|
||
do i = 2, ng
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al(i,1)
|
||
wc = dc * ac(i,1)
|
||
wr = dr * ar(i,1)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,1) = 0.0d+00
|
||
b(i,1) = 1.0d+00
|
||
c(i,1) = 0.0d+00
|
||
r(i,1) = fl(i)
|
||
|
||
end do ! i = 2, ng
|
||
a(1,1) = 0.0d+00
|
||
b(1,1) = 1.0d+00
|
||
c(1,1) = 0.0d+00
|
||
r(1,1) = 0.5d+00 * (fc(1) + fc(2))
|
||
|
||
! interpolate ghost zones using explicit solver (left-side reconstruction)
|
||
!
|
||
do i = n - ng, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = min(n, i + 2)
|
||
|
||
! calculate weights
|
||
!
|
||
wl = dl * al(i,1)
|
||
wc = dc * ac(i,1)
|
||
wr = dr * ar(i,1)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the left state
|
||
!
|
||
ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i )
|
||
qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1)
|
||
qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2)
|
||
|
||
! calculate the left state
|
||
!
|
||
fl(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,1) = 0.0d+00
|
||
b(i,1) = 1.0d+00
|
||
c(i,1) = 0.0d+00
|
||
r(i,1) = fl(i)
|
||
|
||
end do ! i = n - ng, n - 1
|
||
a(n,1) = 0.0d+00
|
||
b(n,1) = 1.0d+00
|
||
c(n,1) = 0.0d+00
|
||
r(n,1) = fc(n)
|
||
|
||
! interpolate ghost zones using explicit solver (right-side reconstruction)
|
||
!
|
||
do i = 2, ng + 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = max(1, i - 2)
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = i + 2
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar(i,2)
|
||
wc = dc * ac(i,2)
|
||
wr = dr * al(i,2)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,2) = 0.0d+00
|
||
b(i,2) = 1.0d+00
|
||
c(i,2) = 0.0d+00
|
||
r(i,2) = fr(i)
|
||
|
||
end do ! i = 2, ng + 1
|
||
a(1,2) = 0.0d+00
|
||
b(1,2) = 1.0d+00
|
||
c(1,2) = 0.0d+00
|
||
r(1,2) = fc(1)
|
||
|
||
! interpolate ghost zones using explicit solver (right-side reconstruction)
|
||
!
|
||
do i = n - ng + 1, n - 1
|
||
|
||
! prepare neighbour indices
|
||
!
|
||
im2 = i - 2
|
||
im1 = i - 1
|
||
ip1 = i + 1
|
||
ip2 = min(n, i + 2)
|
||
|
||
! normalize weights
|
||
!
|
||
wl = dl * ar(i,2)
|
||
wc = dc * ac(i,2)
|
||
wr = dr * al(i,2)
|
||
ww = (wl + wr) + wc
|
||
wl = wl / ww
|
||
wr = wr / ww
|
||
wc = 1.0d+00 - (wl + wr)
|
||
|
||
! calculate the interpolations of the right state
|
||
!
|
||
ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i )
|
||
qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1)
|
||
qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2)
|
||
|
||
! calculate the right state
|
||
!
|
||
fr(i) = (wl * ql + wr * qr) + wc * qc
|
||
|
||
! prepare coefficients of the tridiagonal system
|
||
!
|
||
a(i,2) = 0.0d+00
|
||
b(i,2) = 1.0d+00
|
||
c(i,2) = 0.0d+00
|
||
r(i,2) = fr(i)
|
||
|
||
end do ! i = n - ng + 1, n - 1
|
||
a(n,2) = 0.0d+00
|
||
b(n,2) = 1.0d+00
|
||
c(n,2) = 0.0d+00
|
||
r(n,2) = 0.5d+00 * (fc(n-1) + fc(n))
|
||
|
||
! solve the tridiagonal system of equations for the left-side interpolation
|
||
!
|
||
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
|
||
|
||
! substitute the left-side values
|
||
!
|
||
fl(1:n ) = u(1:n)
|
||
|
||
! solve the tridiagonal system of equations for the left-side interpolation
|
||
!
|
||
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
|
||
|
||
! substitute the right-side values
|
||
!
|
||
fr(1:n-1) = u(2:n)
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_crweno5ns
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_MP5:
|
||
! --------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 5th order
|
||
! Monotonicity Preserving (MP) method with adjustable dissipation.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Suresh, A. & Huynh, H. T.,
|
||
! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
|
||
! Time Stepping"
|
||
! Journal on Computational Physics,
|
||
! 1997, vol. 136, pp. 83-99,
|
||
! http://dx.doi.org/10.1006/jcph.1997.5745
|
||
! [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
|
||
! "Investigation of low-dissipation monotonicity-preserving scheme
|
||
! for direct numerical simulation of compressible turbulent flows",
|
||
! Computers & Fluids,
|
||
! 2014, vol. 104, pp. 55-72,
|
||
! http://dx.doi.org/10.1016/j.compfluid.2014.07.024
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_mp5(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = 3, n - 3
|
||
u(i) = sum(ce5a(:) * fc(i-2:i+3))
|
||
end do
|
||
|
||
u( 1) = sum(ce2(:) * fc( 1:2))
|
||
u( 2) = sum(ce3(:) * fc( 1:3))
|
||
u(n-2) = sum(ce5(:) * fc(n-4:n))
|
||
u(n-1) = sum(ce3(:) * fc(n-2:n))
|
||
u(n ) = fc(n )
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = 3, n - 3
|
||
u(i) = sum(ce5a(:) * fi(i-2:i+3))
|
||
end do
|
||
|
||
u( 1) = sum(ce2(:) * fi( 1:2))
|
||
u( 2) = sum(ce3(:) * fi( 1:3))
|
||
u(n-2) = sum(ce5(:) * fi(n-4:n))
|
||
u(n-1) = sum(ce3(:) * fi(n-2:n))
|
||
u(n ) = fi(n )
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_mp5
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_MP7:
|
||
! --------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 7th order
|
||
! Monotonicity Preserving (MP) method with adjustable dissipation.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! - See RECONSTRUCT_MP5
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_mp7(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = 4, n - 4
|
||
u(i) = sum(ce7a(:) * fc(i-3:i+4))
|
||
end do
|
||
|
||
u( 1) = sum(ce2(:) * fc( 1:2))
|
||
u( 2) = sum(ce3(:) * fc( 1:3))
|
||
u( 3) = sum(ce5(:) * fc( 1:5))
|
||
u(n-3) = sum(ce7(:) * fc(n-6:n))
|
||
u(n-2) = sum(ce5(:) * fc(n-4:n))
|
||
u(n-1) = sum(ce3(:) * fc(n-2:n))
|
||
u(n ) = fc(n )
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = 4, n - 4
|
||
u(i) = sum(ce7a(:) * fi(i-3:i+4))
|
||
end do
|
||
|
||
u( 1) = sum(ce2(:) * fi( 1:2))
|
||
u( 2) = sum(ce3(:) * fi( 1:3))
|
||
u( 3) = sum(ce5(:) * fi( 1:5))
|
||
u(n-3) = sum(ce7(:) * fi(n-6:n))
|
||
u(n-2) = sum(ce5(:) * fi(n-4:n))
|
||
u(n-1) = sum(ce3(:) * fi(n-2:n))
|
||
u(n ) = fi(n )
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_mp7
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_MP9:
|
||
! --------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 9th order
|
||
! Monotonicity Preserving (MP) method with adjustable dissipation.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! - See RECONSTRUCT_MP5LD
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_mp9(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = 5, n - 5
|
||
u(i) = sum(ce9a(:) * fc(i-4:i+5))
|
||
end do
|
||
|
||
u( 1) = sum(ce2(:) * fc( 1:2))
|
||
u( 2) = sum(ce3(:) * fc( 1:3))
|
||
u( 3) = sum(ce5(:) * fc( 1:5))
|
||
u( 4) = sum(ce7(:) * fc( 1:7))
|
||
u(n-4) = sum(ce9(:) * fc(n-8:n))
|
||
u(n-3) = sum(ce7(:) * fc(n-6:n))
|
||
u(n-2) = sum(ce5(:) * fc(n-4:n))
|
||
u(n-1) = sum(ce3(:) * fc(n-2:n))
|
||
u(n ) = fc(n )
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = 5, n - 5
|
||
u(i) = sum(ce9a(:) * fi(i-4:i+5))
|
||
end do
|
||
|
||
u( 1) = sum(ce2(:) * fi( 1:2))
|
||
u( 2) = sum(ce3(:) * fi( 1:3))
|
||
u( 3) = sum(ce5(:) * fi( 1:5))
|
||
u( 4) = sum(ce7(:) * fi( 1:7))
|
||
u(n-4) = sum(ce9(:) * fi(n-8:n))
|
||
u(n-3) = sum(ce7(:) * fi(n-6:n))
|
||
u(n-2) = sum(ce5(:) * fi(n-4:n))
|
||
u(n-1) = sum(ce3(:) * fi(n-2:n))
|
||
u(n ) = fi(n )
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_mp9
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_CRMP5:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 5th order
|
||
! Compact Reconstruction Monotonicity Preserving (CRMP) method with
|
||
! adjustable dissipation.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Suresh, A. & Huynh, H. T.,
|
||
! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
|
||
! Time Stepping"
|
||
! Journal on Computational Physics,
|
||
! 1997, vol. 136, pp. 83-99,
|
||
! http://dx.doi.org/10.1006/jcph.1997.5745
|
||
! [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
|
||
! "Investigation of low-dissipation monotonicity-preserving scheme
|
||
! for direct numerical simulation of compressible turbulent flows",
|
||
! Computers & Fluids,
|
||
! 2014, vol. 104, pp. 55-72,
|
||
! http://dx.doi.org/10.1016/j.compfluid.2014.07.024
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_crmp5(h, fc, fl, fr)
|
||
|
||
use algebra, only : tridiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: a, b, c
|
||
real(kind=8), dimension(size(fc)) :: r
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
do i = 1, ng
|
||
a(i) = 0.0d+00
|
||
b(i) = 1.0d+00
|
||
c(i) = 0.0d+00
|
||
end do
|
||
do i = ng + 1, n - ng - 1
|
||
a(i) = di5a(1)
|
||
b(i) = di5a(2)
|
||
c(i) = di5a(3)
|
||
end do
|
||
do i = n - ng, n
|
||
a(i) = 0.0d+00
|
||
b(i) = 1.0d+00
|
||
c(i) = 0.0d+00
|
||
end do
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci5a(:) * fc(i-1:i+2))
|
||
end do
|
||
|
||
r( 1) = sum(ce2(:) * fc( 1: 2))
|
||
r( 2) = sum(ce3(:) * fc( 1: 3))
|
||
do i = 3, ng
|
||
r(i) = sum(ce5(:) * fc(i-2:i+2))
|
||
end do
|
||
do i = n - ng, n - 2
|
||
r(i) = sum(ce5(:) * fc(i-2:i+2))
|
||
end do
|
||
r(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
r(n ) = fc(n )
|
||
|
||
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci5a(:) * fi(i-1:i+2))
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||
do i = 3, ng
|
||
r(i) = sum(ce5(:) * fi(i-2:i+2))
|
||
end do
|
||
do i = n - ng, n - 2
|
||
r(i) = sum(ce5(:) * fi(i-2:i+2))
|
||
end do
|
||
r(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
r(n ) = fi(n )
|
||
|
||
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_crmp5
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_CRMP7:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 7th order
|
||
! Compact Reconstruction Monotonicity Preserving (CRMP) method with
|
||
! adjustable dissipation.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! - See RECONSTRUCT_CRMP5
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_crmp7(h, fc, fl, fr)
|
||
|
||
use algebra, only : tridiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: a, b, c
|
||
real(kind=8), dimension(size(fc)) :: r
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
do i = 1, ng
|
||
a(i) = 0.0d+00
|
||
b(i) = 1.0d+00
|
||
c(i) = 0.0d+00
|
||
end do
|
||
do i = ng + 1, n - ng - 1
|
||
a(i) = di7a(1)
|
||
b(i) = di7a(2)
|
||
c(i) = di7a(3)
|
||
end do
|
||
do i = n - ng, n
|
||
a(i) = 0.0d+00
|
||
b(i) = 1.0d+00
|
||
c(i) = 0.0d+00
|
||
end do
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci7a(:) * fc(i-2:i+3))
|
||
end do
|
||
|
||
r( 1) = sum(ce2(:) * fc( 1: 2))
|
||
r( 2) = sum(ce3(:) * fc( 1: 3))
|
||
r( 3) = sum(ce5(:) * fc( 1: 5))
|
||
do i = 4, ng
|
||
r(i) = sum(ce7(:) * fc(i-3:i+3))
|
||
end do
|
||
do i = n - ng, n - 3
|
||
r(i) = sum(ce7(:) * fc(i-3:i+3))
|
||
end do
|
||
r(n-2) = sum(ce5(:) * fc(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
r(n ) = fc(n )
|
||
|
||
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci7a(:) * fi(i-2:i+3))
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||
r( 3) = sum(ce5(:) * fi( 1: 5))
|
||
do i = 4, ng
|
||
r(i) = sum(ce7(:) * fi(i-3:i+3))
|
||
end do
|
||
do i = n - ng, n - 3
|
||
r(i) = sum(ce7(:) * fi(i-3:i+3))
|
||
end do
|
||
r(n-2) = sum(ce5(:) * fi(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
r(n ) = fi(n )
|
||
|
||
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_crmp7
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_CRMP9:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 9th order
|
||
! Compact Reconstruction Monotonicity Preserving (CRMP) method with
|
||
! adjustable dissipation.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! - See RECONSTRUCT_CRMP5
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_crmp9(h, fc, fl, fr)
|
||
|
||
use algebra, only : pentadiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: e, c, d, a, b
|
||
real(kind=8), dimension(size(fc)) :: r
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
do i = 1, ng
|
||
e(i) = 0.0d+00
|
||
c(i) = 0.0d+00
|
||
d(i) = 1.0d+00
|
||
a(i) = 0.0d+00
|
||
b(i) = 0.0d+00
|
||
end do
|
||
do i = ng + 1, n - ng - 1
|
||
e(i) = di9a(1)
|
||
c(i) = di9a(2)
|
||
d(i) = di9a(3)
|
||
a(i) = di9a(4)
|
||
b(i) = di9a(5)
|
||
end do
|
||
do i = n - ng, n
|
||
e(i) = 0.0d+00
|
||
c(i) = 0.0d+00
|
||
d(i) = 1.0d+00
|
||
a(i) = 0.0d+00
|
||
b(i) = 0.0d+00
|
||
end do
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci9a(:) * fc(i-2:i+3))
|
||
end do
|
||
|
||
r( 1) = sum(ce2(:) * fc( 1: 2))
|
||
r( 2) = sum(ce3(:) * fc( 1: 3))
|
||
r( 3) = sum(ce5(:) * fc( 1: 5))
|
||
r( 4) = sum(ce7(:) * fc( 1: 7))
|
||
do i = 5, ng
|
||
r(i) = sum(ce9(:) * fc(i-4:i+4))
|
||
end do
|
||
do i = n - ng, n - 4
|
||
r(i) = sum(ce9(:) * fc(i-4:i+4))
|
||
end do
|
||
r(n-3) = sum(ce7(:) * fc(n-6: n))
|
||
r(n-2) = sum(ce5(:) * fc(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
r(n ) = fc(n )
|
||
|
||
call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci9a(:) * fi(i-2:i+3))
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||
r( 3) = sum(ce5(:) * fi( 1: 5))
|
||
r( 4) = sum(ce7(:) * fi( 1: 7))
|
||
do i = 5, ng
|
||
r(i) = sum(ce9(:) * fi(i-4:i+4))
|
||
end do
|
||
do i = n - ng, n - 4
|
||
r(i) = sum(ce9(:) * fi(i-4:i+4))
|
||
end do
|
||
r(n-3) = sum(ce7(:) * fi(n-6: n))
|
||
r(n-2) = sum(ce5(:) * fi(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
r(n ) = fi(n )
|
||
|
||
call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_crmp9
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_OCMP5:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 5th order Optimized
|
||
! Compact Reconstruction Monotonicity Preserving (CRMP) method.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Myeong-Hwan Ahn, Duck-Joo Lee,
|
||
! "Modified Monotonicity Preserving Constraints for High-Resolution
|
||
! Optimized Compact Scheme",
|
||
! Journal of Scientific Computing,
|
||
! 2020, vol. 83, p. 34
|
||
! https://doi.org/10.1007/s10915-020-01221-0
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_ocmp5(h, fc, fl, fr)
|
||
|
||
use algebra, only : tridiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: a, b, c
|
||
real(kind=8), dimension(size(fc)) :: r
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
real(kind=8), parameter :: a1 = 5.0160152430504783780454277648572d-01
|
||
real(kind=8), parameter :: a2 = 2.5395800542237252167985944322018d-01
|
||
real(kind=8), dimension(3), parameter :: &
|
||
di5 = [ a1, 1.0d+00, a2 ]
|
||
real(kind=8), dimension(5), parameter :: &
|
||
ci5 = [- 3.0d+00 * a1 - 3.0d+00 * a2 + 2.0d+00, &
|
||
2.7d+01 * a1 + 1.7d+01 * a2 - 1.3d+01, &
|
||
4.7d+01 * a1 - 4.3d+01 * a2 + 4.7d+01, &
|
||
- 1.3d+01 * a1 + 7.7d+01 * a2 + 2.7d+01, &
|
||
2.0d+00 * a1 + 1.2d+01 * a2 - 3.0d+00 ] / 6.0d+01
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
do i = 1, ng
|
||
a(i) = 0.0d+00
|
||
b(i) = 1.0d+00
|
||
c(i) = 0.0d+00
|
||
end do
|
||
do i = ng + 1, n - ng - 1
|
||
a(i) = di5(1)
|
||
b(i) = di5(2)
|
||
c(i) = di5(3)
|
||
end do
|
||
do i = n - ng, n
|
||
a(i) = 0.0d+00
|
||
b(i) = 1.0d+00
|
||
c(i) = 0.0d+00
|
||
end do
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci5(:) * fc(i-2:i+2))
|
||
end do
|
||
|
||
r( 1) = sum(ce2(:) * fc( 1: 2))
|
||
r( 2) = sum(ce3(:) * fc( 1: 3))
|
||
do i = 3, ng
|
||
r(i) = sum(ce5(:) * fc(i-2:i+2))
|
||
end do
|
||
do i = n - ng, n - 2
|
||
r(i) = sum(ce5(:) * fc(i-2:i+2))
|
||
end do
|
||
r(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
r(n ) = fc(n )
|
||
|
||
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci5(:) * fi(i-2:i+2))
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||
do i = 3, ng
|
||
r(i) = sum(ce5(:) * fi(i-2:i+2))
|
||
end do
|
||
do i = n - ng, n - 2
|
||
r(i) = sum(ce5(:) * fi(i-2:i+2))
|
||
end do
|
||
r(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
r(n ) = fi(n )
|
||
|
||
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_ocmp5
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_OCMP7:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 7th order Optimized
|
||
! Compact Reconstruction Monotonicity Preserving (CRMP) method.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Myeong-Hwan Ahn, Duck-Joo Lee,
|
||
! "Modified Monotonicity Preserving Constraints for High-Resolution
|
||
! Optimized Compact Scheme",
|
||
! Journal of Scientific Computing,
|
||
! 2020, vol. 83, p. 34
|
||
! https://doi.org/10.1007/s10915-020-01221-0
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_ocmp7(h, fc, fl, fr)
|
||
|
||
use algebra, only : pentadiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: e, c, d, a, b
|
||
real(kind=8), dimension(size(fc)) :: r
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
real(kind=8), parameter :: a1 = 6.7080731007136045597206649140123d-01
|
||
real(kind=8), parameter :: a2 = 3.4031334443871704940876966367699d-01
|
||
real(kind=8), dimension(5), parameter :: &
|
||
di7 = [ (3.0d+00 * a1 + 2.0d+00 * a2 - 2.0d+00) / 8.0d+00, &
|
||
a1, 1.0d+00, a2, &
|
||
( a1 + 6.0d+00 * a2 - 2.0d+00) / 4.0d+01 ]
|
||
real(kind=8), dimension(5), parameter :: &
|
||
ci7 = [ 1.80d+01 * a1 + 1.80d+01 * a2 - 1.60d+01, &
|
||
5.43d+02 * a1 + 2.68d+02 * a2 - 2.91d+02, &
|
||
3.43d+02 * a1 - 3.32d+02 * a2 + 5.09d+02, &
|
||
-1.07d+02 * a1 + 5.68d+02 * a2 + 3.09d+02, &
|
||
4.30d+01 * a1 + 3.18d+02 * a2 - 9.10d+01 ] / 6.0d+02
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! prepare the diagonals of the tridiagonal matrix
|
||
!
|
||
do i = 1, ng
|
||
e(i) = 0.0d+00
|
||
c(i) = 0.0d+00
|
||
d(i) = 1.0d+00
|
||
a(i) = 0.0d+00
|
||
b(i) = 0.0d+00
|
||
end do
|
||
do i = ng + 1, n - ng - 1
|
||
e(i) = di7(1)
|
||
c(i) = di7(2)
|
||
d(i) = di7(3)
|
||
a(i) = di7(4)
|
||
b(i) = di7(5)
|
||
end do
|
||
do i = n - ng, n
|
||
e(i) = 0.0d+00
|
||
c(i) = 0.0d+00
|
||
d(i) = 1.0d+00
|
||
a(i) = 0.0d+00
|
||
b(i) = 0.0d+00
|
||
end do
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci7(:) * fc(i-2:i+2))
|
||
end do
|
||
|
||
r( 1) = sum(ce2(:) * fc( 1: 2))
|
||
r( 2) = sum(ce3(:) * fc( 1: 3))
|
||
r( 3) = sum(ce5(:) * fc( 1: 5))
|
||
do i = 4, ng
|
||
r(i) = sum(ce7(:) * fc(i-3:i+3))
|
||
end do
|
||
do i = n - ng, n - 3
|
||
r(i) = sum(ce7(:) * fc(i-3:i+3))
|
||
end do
|
||
r(n-2) = sum(ce5(:) * fc(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
r(n ) = fc(n )
|
||
|
||
call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci7(:) * fi(i-2:i+2))
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||
r( 3) = sum(ce5(:) * fi( 1: 5))
|
||
do i = 4, ng
|
||
r(i) = sum(ce7(:) * fi(i-3:i+3))
|
||
end do
|
||
do i = n - ng, n - 3
|
||
r(i) = sum(ce7(:) * fi(i-3:i+3))
|
||
end do
|
||
r(n-2) = sum(ce5(:) * fi(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
r(n ) = fi(n )
|
||
|
||
call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_ocmp7
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_OCMP9:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using the 9th order Optimized
|
||
! Compact Reconstruction Monotonicity Preserving (CRMP) method.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Myeong-Hwan Ahn, Duck-Joo Lee,
|
||
! "Modified Monotonicity Preserving Constraints for High-Resolution
|
||
! Optimized Compact Scheme",
|
||
! Journal of Scientific Computing,
|
||
! 2020, vol. 83, p. 34
|
||
! https://doi.org/10.1007/s10915-020-01221-0
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_ocmp9(h, fc, fl, fr)
|
||
|
||
use algebra, only : pentadiag
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: e, c, d, a, b
|
||
real(kind=8), dimension(size(fc)) :: r
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
real(kind=8), parameter :: a1 = 6.7049430606778809680492458975703d-01
|
||
real(kind=8), parameter :: a2 = 4.0910274527161221589360057231066d-01
|
||
real(kind=8), dimension(5), parameter :: &
|
||
di9 = [ (9.0d+00 * a1 + 5.0d+00 * a2 - 6.0d+00) / 2.0d+01, &
|
||
a1, 1.0d+00, a2, &
|
||
( a1 + 5.0d+00 * a2 - 2.0d+00) / 2.0d+01 ]
|
||
real(kind=8), dimension(7), parameter :: &
|
||
ci9 = [-1.000d+01 * a1 -1.000d+01 * a2 +1.000d+01, &
|
||
2.420d+02 * a1 +2.000d+02 * a2 -2.140d+02, &
|
||
4.085d+03 * a1 +1.635d+03 * a2 -2.146d+03, &
|
||
2.335d+03 * a1 -1.865d+03 * a2 +3.454d+03, &
|
||
-8.150d+02 * a1 +3.385d+03 * a2 +2.404d+03, &
|
||
4.450d+02 * a1 +2.895d+03 * a2 -9.560d+02, &
|
||
1.800d+01 * a1 +6.000d+01 * a2 -3.200d+01 ] / 4.2d+03
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
! prepare the diagonals of the tridiagonal matrix
|
||
!
|
||
do i = 1, ng
|
||
e(i) = 0.0d+00
|
||
c(i) = 0.0d+00
|
||
d(i) = 1.0d+00
|
||
a(i) = 0.0d+00
|
||
b(i) = 0.0d+00
|
||
end do
|
||
do i = ng + 1, n - ng - 1
|
||
e(i) = di9(1)
|
||
c(i) = di9(2)
|
||
d(i) = di9(3)
|
||
a(i) = di9(4)
|
||
b(i) = di9(5)
|
||
end do
|
||
do i = n - ng, n
|
||
e(i) = 0.0d+00
|
||
c(i) = 0.0d+00
|
||
d(i) = 1.0d+00
|
||
a(i) = 0.0d+00
|
||
b(i) = 0.0d+00
|
||
end do
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci9(:) * fc(i-3:i+3))
|
||
end do
|
||
|
||
r( 1) = sum(ce2(:) * fc( 1: 2))
|
||
r( 2) = sum(ce3(:) * fc( 1: 3))
|
||
r( 3) = sum(ce5(:) * fc( 1: 5))
|
||
r( 4) = sum(ce7(:) * fc( 1: 7))
|
||
do i = 5, ng
|
||
r(i) = sum(ce9(:) * fc(i-4:i+4))
|
||
end do
|
||
do i = n - ng, n - 4
|
||
r(i) = sum(ce9(:) * fc(i-4:i+4))
|
||
end do
|
||
r(n-3) = sum(ce7(:) * fc(n-6: n))
|
||
r(n-2) = sum(ce5(:) * fc(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
r(n ) = fc(n )
|
||
|
||
call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
|
||
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
do i = ng, n - ng + 1
|
||
r(i) = sum(ci9(:) * fi(i-3:i+3))
|
||
end do ! i = ng, n - ng + 1
|
||
|
||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||
r( 3) = sum(ce5(:) * fi( 1: 5))
|
||
r( 4) = sum(ce7(:) * fi( 1: 7))
|
||
do i = 5, ng
|
||
r(i) = sum(ce9(:) * fi(i-4:i+4))
|
||
end do
|
||
do i = n - ng, n - 4
|
||
r(i) = sum(ce9(:) * fi(i-4:i+4))
|
||
end do
|
||
r(n-3) = sum(ce7(:) * fi(n-6: n))
|
||
r(n-2) = sum(ce5(:) * fi(n-4: n))
|
||
r(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
r(n ) = fi(n )
|
||
|
||
call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
|
||
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_ocmp9
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PREPARE_GP:
|
||
! ---------------------
|
||
!
|
||
! Subroutine prepares matrixes for the Gaussian Process (GP) method.
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prepare_gp(iret)
|
||
|
||
use algebra , only : max_prec, invert
|
||
use constants, only : pi
|
||
use helpers , only : print_message
|
||
|
||
implicit none
|
||
|
||
integer, intent(inout) :: iret
|
||
|
||
logical :: flag
|
||
integer :: i, j, ip, jp
|
||
real(kind=max_prec) :: sig, zl, zr, fc
|
||
|
||
real(kind=max_prec), dimension(:,:), allocatable :: cov, agp
|
||
real(kind=max_prec), dimension(:) , allocatable :: xgp
|
||
|
||
character(len=*), parameter :: loc = 'INTERPOLATIONS::prepare_gp()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate normal distribution sigma
|
||
!
|
||
sig = sqrt(2.0d+00) * sgp
|
||
|
||
! allocate the convariance matrix and interpolation position vector
|
||
!
|
||
allocate(cov(ngp,ngp))
|
||
allocate(agp(ngp,ngp))
|
||
allocate(xgp(ngp))
|
||
|
||
! prepare the covariance matrix
|
||
!
|
||
fc = 0.5d+00 * sqrt(pi) * sig
|
||
i = 0
|
||
do ip = - mgp, mgp
|
||
i = i + 1
|
||
j = 0
|
||
do jp = - mgp, mgp
|
||
j = j + 1
|
||
zl = (1.0d+00 * (ip - jp) - 0.5d+00) / sig
|
||
zr = (1.0d+00 * (ip - jp) + 0.5d+00) / sig
|
||
cov(i,j) = fc * (erf(zr) - erf(zl))
|
||
end do
|
||
end do
|
||
|
||
! invert the matrix
|
||
!
|
||
call invert(ngp, cov(:,:), agp(:,:), flag)
|
||
|
||
! continue if the matrix inversion was successful
|
||
!
|
||
if (flag) then
|
||
|
||
! make the inversed matrix symmetric
|
||
!
|
||
do j = 1, ngp
|
||
do i = 1, ngp
|
||
if (i /= j) then
|
||
fc = 0.5d+00 * (agp(i,j) + agp(j,i))
|
||
agp(i,j) = fc
|
||
agp(j,i) = fc
|
||
end if
|
||
end do
|
||
end do
|
||
|
||
! prepare the interpolation position vector
|
||
!
|
||
j = 0
|
||
do jp = - mgp, mgp
|
||
j = j + 1
|
||
zr = (0.5d+00 - 1.0d+00 * jp) / sig
|
||
xgp(j) = exp(- zr**2)
|
||
end do
|
||
|
||
! prepare the interpolation coefficients vector
|
||
!
|
||
cgp(1:ngp) = real(matmul(xgp(1:ngp), agp(1:ngp,1:ngp)), kind=8)
|
||
|
||
end if
|
||
|
||
! deallocate the convariance matrix and interpolation position vector
|
||
!
|
||
deallocate(cov)
|
||
deallocate(agp)
|
||
deallocate(xgp)
|
||
|
||
! check if the matrix was inverted successfully
|
||
!
|
||
if (.not. flag) then
|
||
call print_message(loc, "Could not invert the covariance matrix!")
|
||
iret = 202
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prepare_gp
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_GP:
|
||
! -------------------------
|
||
!
|
||
! Subroutine reconstructs the interface states using one dimensional
|
||
! high order Gaussian Process (GP) method.
|
||
!
|
||
! Arguments are described in subroutine reconstruct().
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_gp(h, fc, fl, fr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
integer :: n, i
|
||
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(fc)
|
||
|
||
!! === left-side interpolation ===
|
||
!!
|
||
! interpolate the interface state of the ghost zones using the interpolations
|
||
! of lower orders
|
||
!
|
||
u( 1) = sum(ce2(:) * fc( 1: 2))
|
||
u( 2) = sum(ce3(:) * fc( 1: 3))
|
||
do i = 3, mgp
|
||
u(i) = sum(ce5(:) * fc(i-2:i+2))
|
||
end do
|
||
do i = n + 1 - mgp, n - 2
|
||
u(i) = sum(ce5(:) * fc(i-2:i+2))
|
||
end do
|
||
u(n-1) = sum(ce3(:) * fc(n-2: n))
|
||
u(n ) = fc(n )
|
||
|
||
! reconstruct the interface state using the Gauss process interpolation
|
||
!
|
||
do i = 1 + mgp, n - mgp
|
||
u(i) = sum(cgp(1:ngp) * (fc(i-mgp:i+mgp) - fc(i))) + fc(i)
|
||
end do
|
||
|
||
! apply the monotonicity preserving limiting
|
||
!
|
||
call mp_limiting(fc(:), u(:))
|
||
|
||
! copy the interpolation to the respective vector
|
||
!
|
||
fl(1:n) = u(1:n)
|
||
|
||
!! === right-side interpolation ===
|
||
!!
|
||
! invert the cell-centered value vector
|
||
!
|
||
fi(1:n) = fc(n:1:-1)
|
||
|
||
! interpolate the interface state of the ghost zones using the interpolations
|
||
! of lower orders
|
||
!
|
||
u( 1) = sum(ce2(:) * fi( 1: 2))
|
||
u( 2) = sum(ce3(:) * fi( 1: 3))
|
||
do i = 3, mgp
|
||
u(i) = sum(ce5(:) * fi(i-2:i+2))
|
||
end do
|
||
do i = n + 1 - mgp, n - 2
|
||
u(i) = sum(ce5(:) * fi(i-2:i+2))
|
||
end do
|
||
u(n-1) = sum(ce3(:) * fi(n-2: n))
|
||
u(n ) = fi(n )
|
||
|
||
! reconstruct the interface state using the Gauss process interpolation
|
||
!
|
||
do i = 1 + mgp, n - mgp
|
||
u(i) = sum(cgp(1:ngp) * (fi(i-mgp:i+mgp) - fi(i))) + fi(i)
|
||
end do
|
||
|
||
! apply the monotonicity preserving limiting
|
||
!
|
||
call mp_limiting(fi(:), u(:))
|
||
|
||
! copy the interpolation to the respective vector
|
||
!
|
||
fr(1:n-1) = u(n-1:1:-1)
|
||
|
||
! update the interpolation of the first and last points
|
||
!
|
||
i = n - 1
|
||
fl(1) = 0.5d+00 * (fc(1) + fc(2))
|
||
fr(i) = 0.5d+00 * (fc(i) + fc(n))
|
||
fl(n) = fc(n)
|
||
fr(n) = fc(n)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_gp
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_ZERO:
|
||
! ---------------------
|
||
!
|
||
! Function returns zero.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_zero(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = 0.0d+00
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_zero
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_AVERAGE:
|
||
! ------------------------
|
||
!
|
||
! Function returns the average among two arguments.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_average(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = 5.0d-01 * (a + b) * x
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_average
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_MINMOD:
|
||
! -----------------------
|
||
!
|
||
! Function returns the minimum module value among two arguments using
|
||
! minmod limiter.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_minmod(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = 0.5d+00 * (sign(x, a) + sign(x, b)) * min(abs(a), abs(b))
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_minmod
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_MONOTONIZED_CENTRAL:
|
||
! ------------------------------------
|
||
!
|
||
! Function returns the minimum module value among two arguments using
|
||
! the monotonized central TVD limiter.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_monotonized_central(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = (sign(x, a) + sign(x, b)) * min(abs(a), abs(b), 2.5d-01 * abs(a + b))
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_monotonized_central
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_SUPERBEE:
|
||
! -------------------------
|
||
!
|
||
! Function returns the minimum module value among two arguments using
|
||
! superbee limiter.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_superbee(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = 0.5d+00 * (sign(x, a) + sign(x, b)) &
|
||
* max(min(2.0d+00 * abs(a), abs(b)), min(abs(a), 2.0d+00 * abs(b)))
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_superbee
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_VANLEER:
|
||
! ------------------------
|
||
!
|
||
! Function returns the minimum module value among two arguments using
|
||
! van Leer's limiter.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_vanleer(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = a * b
|
||
if (c > 0.0d+00) then
|
||
c = 2.0d+00 * x * c / (a + b)
|
||
else
|
||
c = 0.0d+00
|
||
end if
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_vanleer
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_OSPRE:
|
||
! ----------------------
|
||
!
|
||
! Function returns the minimum module value among two arguments using
|
||
! OSPRE limiter.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_ospre(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = a * b
|
||
c = 1.5d+00 * x * c * (a + b) / max(eps, a * a + c + b * b)
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_ospre
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LIMITER_VANALBADA:
|
||
! --------------------------
|
||
!
|
||
! Function returns the minimum module value among two arguments using
|
||
! van Albada's limiter.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! x - scaling factor;
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
function limiter_vanalbada(x, a, b) result(c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = x * a * b * (a + b) / max(eps, a * a + b * b)
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_vanalbada
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function MINMOD:
|
||
! ===============
|
||
!
|
||
! Function returns the minimum module value among two arguments.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! a, b - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function minmod(a, b)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: a, b
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
minmod = (sign(0.5d+00, a) + sign(0.5d+00, b)) * min(abs(a), abs(b))
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function minmod
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function MINMOD4:
|
||
! ================
|
||
!
|
||
! Function returns the minimum module value among four arguments.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! a, b, c, d - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function minmod4(a, b, c, d)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: a, b, c, d
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
minmod4 = minmod(minmod(a, b), minmod(c, d))
|
||
return
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function minmod4
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function MEDIAN:
|
||
! ===============
|
||
!
|
||
! Function returns the median of three argument values.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! a, b, c - the input values;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function median(a, b, c)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: a, b, c
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
median = a + minmod(b - a, c - a)
|
||
return
|
||
|
||
end function median
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FIX_POSITIVITY:
|
||
! -------------------------
|
||
!
|
||
! Subroutine scans the input arrays of the left and right states for
|
||
! negative values. If it finds one, it repeates the state reconstruction
|
||
! the integrated cell-centered values using the zeroth order interpolation.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qi - the cell centered integrals of variable;
|
||
! ql, qr - the left and right states of variable;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fix_positivity(qi, ql, qr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:), intent(in) :: qi
|
||
real(kind=8), dimension(:), intent(inout) :: ql, qr
|
||
|
||
integer :: n, i, im1, ip1
|
||
|
||
!------------------------------------------------------------------------------
|
||
!
|
||
if (positivity) then
|
||
n = size(qi)
|
||
|
||
do i = 1, n
|
||
if (ql(i) <= 0.0d+00) then
|
||
im1 = max(1, i - 1)
|
||
|
||
ql(i ) = qi(i)
|
||
qr(im1) = qi(i)
|
||
end if ! fl ≤ 0
|
||
|
||
if (qr(i) <= 0.0d+00) then
|
||
ip1 = min(n, i + 1)
|
||
|
||
ql(ip1) = qi(ip1)
|
||
qr(i ) = qi(ip1)
|
||
|
||
end if ! fr ≤ 0
|
||
end do ! i = 1, n
|
||
end if ! positivity == .true.
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fix_positivity
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CLIP_EXTREMA:
|
||
! -----------------------
|
||
!
|
||
! Subroutine scans the reconstructed states and check if they didn't leave
|
||
! the allowed limits. In the case where the limits where exceeded,
|
||
! the states are limited using constant reconstruction.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qi - the cell centered integrals of variable;
|
||
! ql, qr - the left and right states of variable;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine clip_extrema(qi, ql, qr)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:), intent(in) :: qi
|
||
real(kind=8), dimension(:), intent(inout) :: ql, qr
|
||
|
||
integer :: n, i, im1, ip1, ip2
|
||
real(kind=8) :: qmn, qmx
|
||
real(kind=8) :: dql, dqr, dq
|
||
|
||
!------------------------------------------------------------------------------
|
||
!
|
||
n = size(qi)
|
||
|
||
do i = 1, n
|
||
|
||
im1 = max(1, i - 1)
|
||
ip1 = min(n, i + 1)
|
||
|
||
qmn = min(qi(i), qi(ip1))
|
||
qmx = max(qi(i), qi(ip1))
|
||
|
||
if (ql(i) < qmn .or. ql(i) > qmx) then
|
||
dql = qi(i ) - qi(im1)
|
||
dqr = qi(ip1) - qi(i )
|
||
|
||
dq = limiter_clip(0.5d+00, dql, dqr)
|
||
|
||
ql(i ) = qi(i) + dq
|
||
qr(im1) = qi(i) - dq
|
||
end if
|
||
|
||
if (qr(i) < qmn .or. qr(i) > qmx) then
|
||
ip2 = min(n, i + 2)
|
||
|
||
dql = qi(ip1) - qi(i )
|
||
dqr = qi(ip2) - qi(ip1)
|
||
|
||
dq = limiter_clip(0.5d+00, dql, dqr)
|
||
|
||
ql(ip1) = qi(ip1) + dq
|
||
qr(i ) = qi(ip1) - dq
|
||
end if
|
||
|
||
end do ! i = 1, n
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine clip_extrema
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine MP_LIMITING:
|
||
! ----------------------
|
||
!
|
||
! Subroutine applies the monotonicity preserving (MP) limiter to a vector of
|
||
! high order reconstructed interface values.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qc - the vector of cell-centered values;
|
||
! qi - the vector of interface values obtained from the high order
|
||
! interpolation as input and its monotonicity limited values as output;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Suresh, A. & Huynh, H. T.,
|
||
! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
|
||
! Time Stepping"
|
||
! Journal on Computational Physics,
|
||
! 1997, vol. 136, pp. 83-99,
|
||
! http://dx.doi.org/10.1006/jcph.1997.5745
|
||
! [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
|
||
! "A 5th order monotonicity-preserving upwind compact difference
|
||
! scheme",
|
||
! Science China Physics, Mechanics and Astronomy,
|
||
! Volume 54, Issue 3, pp. 511-522,
|
||
! http://dx.doi.org/10.1007/s11433-010-4220-x
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine mp_limiting(qc, qi)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:), intent(in) :: qc
|
||
real(kind=8), dimension(:), intent(inout) :: qi
|
||
|
||
integer :: n, i, im1, ip1, ip2
|
||
real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr
|
||
real(kind=8) :: qlc, qmd, qmp, qmn, qmx, qul
|
||
|
||
real(kind=8), dimension(0:size(qc)+2) :: dm
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(qc)
|
||
|
||
dm(0 ) = 0.0d+00
|
||
dm(1 ) = 0.0d+00
|
||
dm(2:n) = qc(2:n) - qc(1:n-1)
|
||
dm(n+1) = 0.0d+00
|
||
dm(n+2) = 0.0d+00
|
||
|
||
do i = 1, n - 1
|
||
|
||
ip1 = i + 1
|
||
|
||
if (dm(i) * dm(ip1) >= 0.0d+00) then
|
||
dq = kappa * dm(i)
|
||
else
|
||
dq = kbeta * dm(i)
|
||
end if
|
||
|
||
qmp = qc(i) + minmod(dm(ip1), dq)
|
||
ds = (qi(i) - qc(i)) * (qi(i) - qmp)
|
||
|
||
if (ds > eps) then
|
||
|
||
im1 = i - 1
|
||
ip2 = i + 2
|
||
|
||
dm1 = dm(i ) - dm(im1)
|
||
dc0 = dm(ip1) - dm(i )
|
||
dp1 = dm(ip2) - dm(ip1)
|
||
dc4 = 4.0d+00 * dc0
|
||
|
||
dml = 0.5d+00 * minmod4(dc4 - dm1, 4.0d+00 * dm1 - dc0, dc0, dm1)
|
||
dmr = 0.5d+00 * minmod4(dc4 - dp1, 4.0d+00 * dp1 - dc0, dc0, dp1)
|
||
|
||
qmd = qc(i) + 0.5d+00 * dm(ip1) - dmr
|
||
qul = qc(i) + dq
|
||
qlc = qc(i) + 0.5d+00 * dq + dml
|
||
|
||
qmx = max(min(qc(i), qc(ip1), qmd), min(qc(i), qul, qlc))
|
||
qmn = min(max(qc(i), qc(ip1), qmd), max(qc(i), qul, qlc))
|
||
|
||
qi(i) = median(qi(i), qmn, qmx)
|
||
|
||
end if
|
||
|
||
end do ! i = 1, n - 1
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine mp_limiting
|
||
|
||
!===============================================================================
|
||
!
|
||
end module interpolations
|