Grzegorz Kowal
b589ffede9
The coefficients were tuned to keep the same maximum disspersion error while shifting it toward higher frequencies with the scheme order. At the same time the dissipation decreses at high frequencies with the scheme order. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
5763 lines
157 KiB
Fortran
5763 lines
157 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-2020 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
|
||
|
||
#ifdef PROFILE
|
||
! import external subroutines
|
||
!
|
||
use timers, only : set_timer, start_timer, stop_timer
|
||
#endif /* PROFILE */
|
||
|
||
! module variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
#ifdef PROFILE
|
||
! timer indices
|
||
!
|
||
integer , save :: imi, imr, imf, imc
|
||
#endif /* PROFILE */
|
||
|
||
! 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 = .false.
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use iso_fortran_env, only : error_unit
|
||
use parameters , only : get_parameter
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(inout) :: nghosts
|
||
logical, intent(in) :: verbose
|
||
integer, intent(out) :: status
|
||
|
||
! local variables
|
||
!
|
||
character(len=255) :: sreconstruction = "tvd"
|
||
character(len=255) :: tlimiter = "mm"
|
||
character(len=255) :: plimiter = "mm"
|
||
character(len=255) :: climiter = "mm"
|
||
character(len=255) :: mlp_limiting = "off"
|
||
character(len=255) :: positivity_fix = "off"
|
||
character(len=255) :: clip_extrema = "off"
|
||
character(len= 16) :: stmp
|
||
real(kind=8) :: cfl = 0.5d+00
|
||
|
||
! local parameters
|
||
!
|
||
character(len=*), parameter :: loc = 'INTERPOLATIONS:initialize_interpolation()'
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! set timer descriptions
|
||
!
|
||
call set_timer('interpolations:: initialization', imi)
|
||
call set_timer('interpolations:: reconstruction', imr)
|
||
call set_timer('interpolations:: fix positivity', imf)
|
||
call set_timer('interpolations:: clip extrema' , imc)
|
||
|
||
! start accounting time for module initialization/finalization
|
||
!
|
||
call start_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
! reset the status flag
|
||
!
|
||
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) then
|
||
write(error_unit,"('[', a, ']: ', a)") trim(loc), &
|
||
"The parameter ngp has to be an odd integer >= 3. "// &
|
||
"Resetting to default value of ngp = 5."
|
||
end if
|
||
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, 4)
|
||
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) then
|
||
write(error_unit,"('[', a, ']: ', a)") trim(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'."
|
||
end if
|
||
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 default
|
||
name_tlim = "zero derivative"
|
||
limiter_tvd => limiter_zero
|
||
order = max(1, order)
|
||
end select
|
||
|
||
! select the prolongation limiter
|
||
!
|
||
select case(trim(plimiter))
|
||
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 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 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
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for module initialization/finalization
|
||
!
|
||
call stop_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(out) :: status
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! start accounting time for module initialization/finalization
|
||
!
|
||
call start_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
! reset the status flag
|
||
!
|
||
status = 0
|
||
|
||
! 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)
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for module initialization/finalization
|
||
!
|
||
call stop_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! import external procedures and variables
|
||
!
|
||
use helpers, only : print_section, print_parameter
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use algebra , only : max_prec, invert
|
||
use constants , only : pi
|
||
use iso_fortran_env, only : error_unit
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(inout) :: iret
|
||
|
||
! local variables
|
||
!
|
||
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 */
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=max_prec), dimension(:,:) , allocatable :: cov, inv
|
||
real(kind=max_prec), dimension(:,:,:), allocatable :: xgp
|
||
|
||
! local parameters
|
||
!
|
||
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
|
||
write(error_unit,"('[',a,']: ',a)") trim(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)
|
||
|
||
! include external procedures
|
||
!
|
||
use coordinates , only : nb, ne, nbl, neu
|
||
use iso_fortran_env, only : error_unit
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
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
|
||
|
||
! local variables
|
||
!
|
||
integer :: i , j , k = 1
|
||
integer :: im1, jm1, ip1, jp1
|
||
#if NDIMS == 3
|
||
integer :: km1, kp1
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! local vectors
|
||
!
|
||
real(kind=8), dimension(NDIMS) :: dql, dqr, dq
|
||
|
||
! local parameters
|
||
!
|
||
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(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||
, "Positive variable is not positive at ( ", i, j, k, " )"
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use coordinates, only : nb, ne, nbl, neu
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
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
|
||
|
||
! local variables
|
||
!
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use coordinates , only : nb, ne, nbl, neu
|
||
use iso_fortran_env, only : error_unit
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
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
|
||
|
||
! local variables
|
||
!
|
||
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 */
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=8), dimension(NDIMS) :: dql, dqr, dq
|
||
real(kind=8), dimension(dgp) :: u
|
||
|
||
! local parameters
|
||
!
|
||
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(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||
, "Positive variable is not positive at ( ", i, j, k, " )"
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use coordinates, only : nn => bcells
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8), dimension(:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:), intent(inout) :: qi
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local vectors
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: f
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! start accounting time for reconstruction
|
||
!
|
||
call start_timer(imr)
|
||
#endif /* PROFILE */
|
||
|
||
! 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(:))
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for reconstruction
|
||
!
|
||
call stop_timer(imr)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
integer :: n, i, im1, ip1
|
||
real(kind=8) :: df, dfl, dfr
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
integer :: n, i, im1, ip1
|
||
real(kind=8) :: bp, bm, ap, am, wp, wm, ww
|
||
real(kind=8) :: dfl, dfr, fp, fm, h2
|
||
|
||
! selection weights
|
||
!
|
||
real(kind=8), parameter :: dp = 2.0d+00 / 3.0d+00, dm = 1.0d+00 / 3.0d+00
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
integer :: n, i, im1, ip1
|
||
real(kind=8) :: dfl, dfr
|
||
real(kind=8) :: th, et, f1, f2, xl, xi, rdx, rdx2
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays
|
||
!
|
||
real(kind=8), dimension(size(fc)) :: fi, df2
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
|
||
! smoothness indicator coefficients
|
||
!
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
! weight coefficients
|
||
!
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
! interpolation coefficients
|
||
!
|
||
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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
|
||
! smoothness indicator coefficients
|
||
!
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
! weight coefficients
|
||
!
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
! interpolation coefficients
|
||
!
|
||
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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=8), dimension(size(fc)) :: dfm, dfp, df2
|
||
|
||
! weight coefficients
|
||
!
|
||
real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01
|
||
|
||
! interpolation coefficients
|
||
!
|
||
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
|
||
|
||
! the free parameter for smoothness indicators (see Eq. 3.6 in [1])
|
||
!
|
||
real(kind=8), parameter :: xi = 4.0d-01
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use algebra , only : tridiag
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
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
|
||
|
||
! smoothness indicator coefficients
|
||
!
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
! weight coefficients for implicit (c) and explicit (d) interpolations
|
||
!
|
||
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
|
||
|
||
! implicit method coefficients
|
||
!
|
||
real(kind=8), parameter :: dq = 5.0d-01
|
||
|
||
! interpolation coefficients
|
||
!
|
||
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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use algebra , only : tridiag
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
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
|
||
|
||
! smoothness indicator coefficients
|
||
!
|
||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01
|
||
|
||
! weight coefficients for implicit (c) and explicit (d) interpolations
|
||
!
|
||
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
|
||
|
||
! implicit method coefficients
|
||
!
|
||
real(kind=8), parameter :: dq = 5.0d-01
|
||
|
||
! interpolation coefficients
|
||
!
|
||
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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use algebra , only : tridiag
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
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
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
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
|
||
|
||
! the free parameter for smoothness indicators (see eq. 3.6 in [3])
|
||
!
|
||
real(kind=8), parameter :: xi = 4.0d-01
|
||
|
||
! weight coefficients for implicit (c) and explicit (d) interpolations
|
||
!
|
||
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
|
||
|
||
! implicit method coefficients
|
||
!
|
||
real(kind=8), parameter :: dq = 5.0d-01
|
||
|
||
! 3rd order interpolation coefficients for three stencils
|
||
!
|
||
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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! include external procedures
|
||
!
|
||
use algebra , only : max_prec, invert
|
||
use constants , only : pi
|
||
use iso_fortran_env, only : error_unit
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(inout) :: iret
|
||
|
||
! local variables
|
||
!
|
||
logical :: flag
|
||
integer :: i, j, ip, jp
|
||
real(kind=max_prec) :: sig, zl, zr, fc
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=max_prec), dimension(:,:), allocatable :: cov, agp
|
||
real(kind=max_prec), dimension(:) , allocatable :: xgp
|
||
|
||
! local parameters
|
||
!
|
||
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
|
||
write(error_unit,"('[',a,']: ',a)") trim(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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8) , intent(in) :: h
|
||
real(kind=8), dimension(:), intent(in) :: fc
|
||
real(kind=8), dimension(:), intent(out) :: fl, fr
|
||
|
||
! local variables
|
||
!
|
||
integer :: n, i
|
||
|
||
! local arrays for derivatives
|
||
!
|
||
real(kind=8), dimension(size(fc)) :: fi
|
||
real(kind=8), dimension(size(fc)) :: u
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! get the input vector length
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = 0.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_zero
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! 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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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)))
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function limiter_vanleer
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! 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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
real(kind=8), intent(in) :: x, a, b
|
||
real(kind=8) :: c
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
c = x * a * b * (a + b) / max(eps, a * a + b * b)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input arguments
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), dimension(:), intent(in) :: qi
|
||
real(kind=8), dimension(:), intent(inout) :: ql, qr
|
||
|
||
! local variables
|
||
!
|
||
integer :: n, i, im1, ip1
|
||
!
|
||
!------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! start accounting time for positivity fix
|
||
!
|
||
call start_timer(imf)
|
||
#endif /* PROFILE */
|
||
|
||
! check positivity only if desired
|
||
!
|
||
if (positivity) then
|
||
|
||
! get the input vector size
|
||
!
|
||
n = size(qi)
|
||
|
||
! look for negative values in the states along the vector
|
||
!
|
||
do i = 1, n
|
||
|
||
! check if the left state has a negative value
|
||
!
|
||
if (ql(i) <= 0.0d+00) then
|
||
|
||
! calculate the left neighbour index
|
||
!
|
||
im1 = max(1, i - 1)
|
||
|
||
! limit the states using the zeroth-order reconstruction
|
||
!
|
||
ql(i ) = qi(i)
|
||
qr(im1) = qi(i)
|
||
|
||
end if ! fl ≤ 0
|
||
|
||
! check if the right state has a negative value
|
||
!
|
||
if (qr(i) <= 0.0d+00) then
|
||
|
||
! calculate the right neighbour index
|
||
!
|
||
ip1 = min(n, i + 1)
|
||
|
||
! limit the states using the zeroth-order reconstruction
|
||
!
|
||
ql(ip1) = qi(ip1)
|
||
qr(i ) = qi(ip1)
|
||
|
||
end if ! fr ≤ 0
|
||
|
||
end do ! i = 1, n
|
||
|
||
end if ! positivity == .true.
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for positivity fix
|
||
!
|
||
call stop_timer(imf)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8), dimension(:), intent(in) :: qi
|
||
real(kind=8), dimension(:), intent(inout) :: ql, qr
|
||
|
||
! local variables
|
||
!
|
||
integer :: n, i, im1, ip1, ip2
|
||
real(kind=8) :: qmn, qmx
|
||
real(kind=8) :: dql, dqr, dq
|
||
!
|
||
!------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! start accounting time for extrema clipping
|
||
!
|
||
call start_timer(imc)
|
||
#endif /* PROFILE */
|
||
|
||
! get the input vector size
|
||
!
|
||
n = size(qi)
|
||
|
||
! iterate over all points
|
||
!
|
||
do i = 1, n
|
||
|
||
! calculate indices
|
||
!
|
||
im1 = max(1, i - 1)
|
||
ip1 = min(n, i + 1)
|
||
|
||
! estimate the bounds of the allowed interval for reconstructed states
|
||
!
|
||
qmn = min(qi(i), qi(ip1))
|
||
qmx = max(qi(i), qi(ip1))
|
||
|
||
! check if the left state lays in the allowed range
|
||
!
|
||
if (ql(i) < qmn .or. ql(i) > qmx) then
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
dql = qi(i ) - qi(im1)
|
||
dqr = qi(ip1) - qi(i )
|
||
|
||
! get the limited slope
|
||
!
|
||
dq = limiter_clip(0.5d+00, dql, dqr)
|
||
|
||
! calculate new states
|
||
!
|
||
ql(i ) = qi(i) + dq
|
||
qr(im1) = qi(i) - dq
|
||
|
||
end if
|
||
|
||
! check if the right state lays in the allowed range
|
||
!
|
||
if (qr(i) < qmn .or. qr(i) > qmx) then
|
||
|
||
! calculate the missing index
|
||
!
|
||
ip2 = min(n, i + 2)
|
||
|
||
! calculate the left and right derivatives
|
||
!
|
||
dql = qi(ip1) - qi(i )
|
||
dqr = qi(ip2) - qi(ip1)
|
||
|
||
! get the limited slope
|
||
!
|
||
dq = limiter_clip(0.5d+00, dql, dqr)
|
||
|
||
! calculate new states
|
||
!
|
||
ql(ip1) = qi(ip1) + dq
|
||
qr(i ) = qi(ip1) - dq
|
||
|
||
end if
|
||
|
||
end do ! i = 1, n
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for extrema clipping
|
||
!
|
||
call stop_timer(imc)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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
|
||
! [3] 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 mp_limiting(qc, qi)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8), dimension(:), intent(in) :: qc
|
||
real(kind=8), dimension(:), intent(inout) :: qi
|
||
|
||
! local variables
|
||
!
|
||
logical :: test
|
||
integer :: n, i, im2, im1, ip1, ip2
|
||
real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr, bt
|
||
real(kind=8) :: qlc, qmd, qmp, qmn, qmx, qul
|
||
|
||
! local vectors
|
||
!
|
||
real(kind=8), dimension(0:size(qc)+2) :: dm
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
n = size(qc)
|
||
|
||
! 1st order derivatives
|
||
!
|
||
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
|
||
|
||
! check the monotonicity condition and apply limiting if necessary
|
||
!
|
||
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
|
||
|
||
im2 = i - 2
|
||
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
|
||
|
||
test = qc(i) > max(qc(im1), qc(ip1)) .and. &
|
||
min(qc(im1), qc(ip1)) > max(qc(im2), qc(ip2))
|
||
test = test .or. qc(i) < min(qc(im1), qc(ip1)) .and. &
|
||
max(qc(im1), qc(ip1)) < min(qc(im2), qc(ip2))
|
||
if (test) then
|
||
if (qc(im2) <= qc(ip1) .and. qc(ip2) <= qc(im1)) then
|
||
bt = dc0 / (qc(ip2) + qc(im2) - 2.0d+00 * qc(i))
|
||
if (bt >= 3.0d-01) qlc = qc(im2)
|
||
end if
|
||
end if
|
||
|
||
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
|