Merge branch 'master' into reconnection
This commit is contained in:
commit
5b1654d8f4
@ -50,7 +50,8 @@ module interpolations
|
||||
! pointers to the reconstruction and limiter procedures
|
||||
!
|
||||
procedure(reconstruct) , pointer, save :: reconstruct_states => null()
|
||||
procedure(limiter_zero), pointer, save :: limiter => null()
|
||||
procedure(stencil_weights_js), pointer, save :: stencil_weights => null()
|
||||
procedure(limiter_zero) , pointer, save :: limiter => null()
|
||||
|
||||
! module parameters
|
||||
!
|
||||
@ -105,10 +106,12 @@ module interpolations
|
||||
! local variables
|
||||
!
|
||||
character(len=255) :: sreconstruction = "tvd"
|
||||
character(len=255) :: sweights = "yc"
|
||||
character(len=255) :: slimiter = "mm"
|
||||
character(len=255) :: positivity_fix = "off"
|
||||
character(len=255) :: clip_extrema = "off"
|
||||
character(len=255) :: name_rec = ""
|
||||
character(len=255) :: name_wei = ""
|
||||
character(len=255) :: name_lim = ""
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -127,9 +130,10 @@ module interpolations
|
||||
|
||||
! obtain the user defined interpolation methods and coefficients
|
||||
!
|
||||
call get_parameter_string("reconstruction", sreconstruction)
|
||||
call get_parameter_string("reconstruction" , sreconstruction)
|
||||
call get_parameter_string("stencil_weights", sweights )
|
||||
call get_parameter_string("limiter" , slimiter )
|
||||
call get_parameter_string("fix_positivity", positivity_fix )
|
||||
call get_parameter_string("fix_positivity" , positivity_fix )
|
||||
call get_parameter_string("clip_extrema" , clip_extrema )
|
||||
call get_parameter_real ("eps" , eps )
|
||||
call get_parameter_real ("limo3_rad" , rad )
|
||||
@ -146,6 +150,9 @@ module interpolations
|
||||
case ("limo3", "LIMO3", "LimO3")
|
||||
name_rec = "3rd order logarithmic limited"
|
||||
reconstruct_states => reconstruct_limo3
|
||||
case ("crweno5", "CRWENO5")
|
||||
name_rec = "5th order Compact WENO"
|
||||
reconstruct_states => reconstruct_crweno5
|
||||
case default
|
||||
if (verbose) then
|
||||
write (*,"(1x,a)") "The selected reconstruction method is not " // &
|
||||
@ -154,6 +161,26 @@ module interpolations
|
||||
end if
|
||||
end select
|
||||
|
||||
! select the stencil weights
|
||||
!
|
||||
select case(trim(sweights))
|
||||
case ("js", "JS")
|
||||
name_wei = "Jiang-Shu"
|
||||
stencil_weights => stencil_weights_js
|
||||
case ("z", "Z")
|
||||
name_wei = "Borges et al."
|
||||
stencil_weights => stencil_weights_z
|
||||
case ("yc", "YC")
|
||||
name_wei = "Yamaleev-Carpenter"
|
||||
stencil_weights => stencil_weights_yc
|
||||
case default
|
||||
if (verbose) then
|
||||
write (*,"(1x,a)") "The selected stencil weight method is not " // &
|
||||
"implemented: " // trim(sweights)
|
||||
stop
|
||||
end if
|
||||
end select
|
||||
|
||||
! select the limiter
|
||||
!
|
||||
select case(trim(slimiter))
|
||||
@ -196,10 +223,15 @@ module interpolations
|
||||
!
|
||||
if (verbose) then
|
||||
|
||||
write (*,"(4x,a14,9x,'=',1x,a)") "reconstruction", trim(name_rec)
|
||||
write (*,"(4x,a14,9x,'=',1x,a)") "limiter ", trim(name_lim)
|
||||
write (*,"(4x,a14,9x,'=',1x,a)") "fix positivity", trim(positivity_fix)
|
||||
write (*,"(4x,a14,9x,'=',1x,a)") "clip extrema ", trim(clip_extrema)
|
||||
write (*,"(4x,a15,8x,'=',1x,a)") "reconstruction ", trim(name_rec)
|
||||
select case(trim(sreconstruction))
|
||||
case ("crweno5", "CRWENO5")
|
||||
write (*,"(4x,a15,8x,'=',1x,a)") "stencil weights", trim(name_wei)
|
||||
case default
|
||||
end select
|
||||
write (*,"(4x,a15,8x,'=',1x,a)") "limiter ", trim(name_lim)
|
||||
write (*,"(4x,a15,8x,'=',1x,a)") "fix positivity ", trim(positivity_fix)
|
||||
write (*,"(4x,a15,8x,'=',1x,a)") "clip extrema ", trim(clip_extrema)
|
||||
|
||||
end if
|
||||
|
||||
@ -648,6 +680,560 @@ module interpolations
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RECONSTRUCT_CRWENO5:
|
||||
! ------------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the fifth order
|
||||
! Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
|
||||
! method.
|
||||
!
|
||||
! 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
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine reconstruct_crweno5(n, h, f, fl, fr)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8) , intent(in) :: h
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(out) :: fl, fr
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(1:n) :: bl, bc, br
|
||||
real(kind=8), dimension(1:n) :: wl, wc, wr
|
||||
real(kind=8), dimension(1:n) :: u
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate smoothness indicators
|
||||
!
|
||||
call smoothness_indicators(n, f(1:n), bl(1:n), bc(1:n), br(1:n))
|
||||
|
||||
! calculate stencil weights
|
||||
!
|
||||
call stencil_weights(n, f(1:n), bl(1:n), bc(1:n), br(1:n) &
|
||||
, wl(1:n), wc(1:n), wr(1:n))
|
||||
|
||||
! find the left state interpolation implicitelly
|
||||
!
|
||||
call crweno5_implicit(n, f(1:n), wl(1:n), wc(1:n), wr(1:n), u(1:n))
|
||||
|
||||
! substitute the left state
|
||||
!
|
||||
fl(1:n) = u(1:n)
|
||||
|
||||
! find the right state interpolation implicitelly
|
||||
!
|
||||
call crweno5_implicit(n, f(n:1:-1), wr(n:1:-1), wc(n:1:-1), wl(n:1:-1) &
|
||||
, u(n:1:-1))
|
||||
|
||||
! substitute the right state
|
||||
!
|
||||
fr(1:n-1) = u (2:n)
|
||||
fr( n ) = fl( n)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine reconstruct_crweno5
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine CRWENO5_IMPLICIT:
|
||||
! ---------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the fifth order
|
||||
! Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
|
||||
! method.
|
||||
!
|
||||
! 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
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine crweno5_implicit(n, f, wl, wc, wr, u)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(in) :: wl, wc, wr
|
||||
real(kind=8), dimension(n), intent(out) :: u
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: im2, im1, i , ip1, ip2
|
||||
integer :: nm1, nm2, nm3, nm4
|
||||
real(kind=8) :: xl, xc, xr, xx
|
||||
real(kind=8) :: bl, bc, br, bt
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(1:n) :: fm, fp
|
||||
real(kind=8), dimension(1:n) :: a, b, c
|
||||
real(kind=8), dimension(1:n) :: r, g
|
||||
|
||||
! local constants
|
||||
!
|
||||
real(kind=8), parameter :: al = 1.0d-01, ac = 6.0d-01, ar = 3.0d-01
|
||||
real(kind=8), parameter :: cl = 2.0d-01, cc = 5.0d-01, cr = 3.0d-01
|
||||
real(kind=8), parameter :: dh = 5.0d-01, ds = 1.0d+00 / 6.0d+00
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate indices
|
||||
!
|
||||
nm1 = n - 1
|
||||
nm2 = n - 2
|
||||
nm3 = n - 3
|
||||
nm4 = n - 4
|
||||
|
||||
! prepare coefficients
|
||||
!
|
||||
do i = 4, nm3
|
||||
|
||||
! prepare indices
|
||||
!
|
||||
im1 = i - 1
|
||||
ip1 = i + 1
|
||||
|
||||
! normalize weigths
|
||||
!
|
||||
xc = cc * wc(i)
|
||||
xl = cl * wl(i)
|
||||
xr = cr * wr(i)
|
||||
xx = xc + (xl + xr)
|
||||
bl = xl / xx
|
||||
br = xr / xx
|
||||
bc = 1.0d+00 - (bl + br)
|
||||
|
||||
! calculate tridiagonal matrix coefficients
|
||||
!
|
||||
a(i) = 2.0d+00 * bl + bc
|
||||
b(i) = bl + 2.0d+00 * (bc + br)
|
||||
c(i) = br
|
||||
|
||||
! prepare right hand side of tridiagonal equation
|
||||
!
|
||||
r(i) = ( bl * f(im1) &
|
||||
+ (5.0d+00 * (bl + bc) + br) * f(i) &
|
||||
+ (bc + 5.0d+00 * br) * f(ip1)) * dh
|
||||
|
||||
end do
|
||||
|
||||
! at the left boundaries, we apply 5th order explicit WENO interpolation with
|
||||
! different weights
|
||||
!
|
||||
! normalize weigths
|
||||
!
|
||||
xc = ac * wc(3)
|
||||
xl = al * wl(3)
|
||||
xr = ar * wr(3)
|
||||
xx = xc + (xl + xr)
|
||||
bl = xl / xx
|
||||
br = xr / xx
|
||||
bc = 1.0d+00 - (bl + br)
|
||||
|
||||
! prepare right hand side of tridiagonal equation
|
||||
!
|
||||
r(1) = dh * (f(1) + f(2 ))
|
||||
r(2) = f(2) + limiter(dh, f(2) - f(1), f(3) - f(2))
|
||||
r(3) = (bl * (2.0d+00 * f(1) - 7.0d+00 * f(2) + 1.1d+01 * f(3)) &
|
||||
+ bc * ( - f(2) + 5.0d+00 * f(3) + 2.0d+00 * f(4)) &
|
||||
+ br * (2.0d+00 * f(3) + 5.0d+00 * f(4) - f(5))) * ds
|
||||
|
||||
! at the right boundaries, we do the similar thing
|
||||
!
|
||||
! normalize weigths
|
||||
!
|
||||
xc = ac * wc(nm2)
|
||||
xl = al * wl(nm2)
|
||||
xr = ar * wr(nm2)
|
||||
xx = xc + (xl + xr)
|
||||
bl = xl / xx
|
||||
br = xr / xx
|
||||
bc = 1.0d+00 - (bl + br)
|
||||
|
||||
! prepare right hand side of tridiagonal equation
|
||||
!
|
||||
r(nm2) = (bl * (2.0d+00 * f(nm4) - 7.0d+00 * f(nm3) + 1.1d+01 * f(nm2)) &
|
||||
+ bc * ( - f(nm3) + 5.0d+00 * f(nm2) + 2.0d+00 * f(nm1)) &
|
||||
+ br * (2.0d+00 * f(nm2) + 5.0d+00 * f(nm1) - f(n ))) &
|
||||
* ds
|
||||
|
||||
r(nm1) = f(nm1) + limiter(dh, f(nm1) - f(nm2), f(n) - f(nm1))
|
||||
r(n ) = f(n ) + dh * (f(n) - f(nm1))
|
||||
|
||||
! correct matrix coefficients for boundaries with explicit interpolation
|
||||
!
|
||||
do i = 1, 3
|
||||
a(i) = 0.0d+00
|
||||
b(i) = 1.0d+00
|
||||
c(i) = 0.0d+00
|
||||
end do
|
||||
|
||||
do i = nm2, n
|
||||
a(i) = 0.0d+00
|
||||
b(i) = 1.0d+00
|
||||
c(i) = 0.0d+00
|
||||
end do
|
||||
|
||||
! solve the tridiagonal system of equations
|
||||
!
|
||||
bt = b(1)
|
||||
u(1) = r(1) / bt
|
||||
do i = 2, n
|
||||
im1 = i - 1
|
||||
g(i) = c(im1) / bt
|
||||
bt = b(i) - a(i) * g(i)
|
||||
u(i) = (r(i) - a(i) * u(im1)) / bt
|
||||
end do
|
||||
do i = nm1, 1, -1
|
||||
ip1 = i + 1
|
||||
u(i) = u(i) - g(ip1) * u(ip1)
|
||||
end do
|
||||
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine crweno5_implicit
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine SMOOTHNESS_INDICATORS:
|
||||
! --------------------------------
|
||||
!
|
||||
! Subroutine calculates smoothness indicators for a given vector of variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! n - the length of the input vector;
|
||||
! f - the input vector of cell averaged values;
|
||||
! bl - the smoothness indicators for the left stencil;
|
||||
! bc - the smoothness indicators for the central stencil;
|
||||
! br - the smoothness indicators for the right stencil;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Jiang, G.-S., Shu, C.-W.,
|
||||
! "Efficient Implementation of Weighted ENO Schemes"
|
||||
! Journal of Computational Physics,
|
||||
! 1996, vol. 126, pp. 202-228,
|
||||
! http://dx.doi.org/10.1006/jcph.1996.0130
|
||||
!
|
||||
! [2] 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
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine smoothness_indicators(n, f, bl, bc, br)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(out) :: bl, bc, br
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: np2, np1, nm1, nm2, i, j
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(0:n+2) :: df2, df1
|
||||
real(kind=8), dimension(1:n) :: dfl, dfc, dfr
|
||||
|
||||
! local constants
|
||||
!
|
||||
real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 1.0d+00 / 4.0d+00
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate indices
|
||||
!
|
||||
np2 = n + 2
|
||||
np1 = n + 1
|
||||
nm1 = n - 1
|
||||
nm2 = n - 2
|
||||
|
||||
! calculate the second order derivative
|
||||
!
|
||||
df2(2:nm1) = ((f(3:n) + f(1:nm2)) - 2.0d+00 * f(2:nm1))**2
|
||||
df2(1 ) = (f(2 ) - f(1))**2
|
||||
df2( n ) = (f(nm1) - f(n))**2
|
||||
df2(0 ) = 0.0d+00
|
||||
df2( np1) = 0.0d+00
|
||||
df2( np2) = 0.0d+00
|
||||
|
||||
! calculate the first derivative
|
||||
!
|
||||
df1(2:n) = f(2:n) - f(1:nm1)
|
||||
df1(0 ) = 0.0d+00
|
||||
df1(1 ) = 0.0d+00
|
||||
df1(np1) = 0.0d+00
|
||||
df1(np2) = 0.0d+00
|
||||
|
||||
! calculate the first order derivatives for left, central and right stencils
|
||||
!
|
||||
dfl( 1:n ) = 3.0d+00 * df1(1:n) - df1(0:nm1)
|
||||
|
||||
dfc( 2:nm1) = f(1:nm2) - f(3:n)
|
||||
dfc( 1 ) = f(1 ) - f(2 )
|
||||
dfc( n ) = f( nm1) - f( n)
|
||||
|
||||
dfr( 1:n) = - 3.0d+00 * df1(2:np1) + df1(3:np2)
|
||||
|
||||
! calculate the left, central and right smoothness indicators
|
||||
!
|
||||
bl(1:n) = c1 * df2(0:nm1) + c2 * dfl(1:n)**2
|
||||
bc(1:n) = c1 * df2(1:n ) + c2 * dfc(1:n)**2
|
||||
br(1:n) = c1 * df2(2:np1) + c2 * dfr(1:n)**2
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine smoothness_indicators
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine STENCIL_WEIGHTS_JS:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine calculate the stencil weights using Jiang-Shu method.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! n - the length of the input vector;
|
||||
! f - the input vector of cell averaged values;
|
||||
! bl - the smoothness indicators for the left stencil;
|
||||
! bc - the smoothness indicators for the central stencil;
|
||||
! br - the smoothness indicators for the right stencil;
|
||||
! wl - the weights the left stencil;
|
||||
! wc - the weights for the central stencil;
|
||||
! wr - the weights for the right stencil;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Jiang, G.-S., Shu, C.-W.,
|
||||
! "Efficient Implementation of Weighted ENO Schemes"
|
||||
! Journal of Computational Physics,
|
||||
! 1996, vol. 126, pp. 202-228,
|
||||
! http://dx.doi.org/10.1006/jcph.1996.0130
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine stencil_weights_js(n, f, bl, bc, br, wl, wc, wr)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(in) :: bl, bc, br
|
||||
real(kind=8), dimension(n), intent(out) :: wl, wc, wr
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate the weights
|
||||
!
|
||||
wl(1:n) = 1.0d+00 / (bl(1:n) + eps)**2
|
||||
wc(1:n) = 1.0d+00 / (bc(1:n) + eps)**2
|
||||
wr(1:n) = 1.0d+00 / (br(1:n) + eps)**2
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine stencil_weights_js
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine STENCIL_WEIGHTS_Z:
|
||||
! ----------------------------
|
||||
!
|
||||
! Subroutine calculate the stencil weights using Borges et al. method.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! n - the length of the input vector;
|
||||
! f - the input vector of cell averaged values;
|
||||
! bl - the smoothness indicators for the left stencil;
|
||||
! bc - the smoothness indicators for the central stencil;
|
||||
! br - the smoothness indicators for the right stencil;
|
||||
! wl - the weights the left stencil;
|
||||
! wc - the weights for the central stencil;
|
||||
! wr - the weights for the right stencil;
|
||||
!
|
||||
! 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 stencil_weights_z(n, f, bl, bc, br, wl, wc, wr)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(in) :: bl, bc, br
|
||||
real(kind=8), dimension(n), intent(out) :: wl, wc, wr
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(1:n) :: tt
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate the factor τ
|
||||
!
|
||||
tt(1:n) = abs(bl(1:n) - br(1:n))
|
||||
|
||||
! calculate the weights
|
||||
!
|
||||
wl(1:n) = 1.0d+00 + (tt(1:n) / (bl(1:n) + eps))**2
|
||||
wc(1:n) = 1.0d+00 + (tt(1:n) / (bc(1:n) + eps))**2
|
||||
wr(1:n) = 1.0d+00 + (tt(1:n) / (br(1:n) + eps))**2
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine stencil_weights_z
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine STENCIL_WEIGHTS_YC:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine calculate the stencil weights using Yalamleev-Carpenter method.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! n - the length of the input vector;
|
||||
! f - the input vector of cell averaged values;
|
||||
! bl - the smoothness indicators for the left stencil;
|
||||
! bc - the smoothness indicators for the central stencil;
|
||||
! br - the smoothness indicators for the right stencil;
|
||||
! wl - the weights the left stencil;
|
||||
! wc - the weights for the central stencil;
|
||||
! wr - the weights for the right stencil;
|
||||
!
|
||||
! 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 stencil_weights_yc(n, f, bl, bc, br, wl, wc, wr)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(in) :: bl, bc, br
|
||||
real(kind=8), dimension(n), intent(out) :: wl, wc, wr
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: im2, im1, i , ip1, ip2
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(1:n) :: tt
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate the factor τ
|
||||
!
|
||||
do i = 1, n
|
||||
im2 = max(1, i - 2)
|
||||
im1 = max(1, i - 1)
|
||||
ip1 = min(n, i + 1)
|
||||
ip2 = min(n, i + 2)
|
||||
|
||||
tt(i) = (6.0d+00 * f(i) + (f(im2) + f(ip2)) &
|
||||
- 4.0d+00 * (f(im1) + f(ip1)))**2
|
||||
|
||||
end do
|
||||
|
||||
! calculate the weights
|
||||
!
|
||||
wl(1:n) = 1.0d+00 + (tt(1:n) / (bl(1:n) + eps))**2
|
||||
wc(1:n) = 1.0d+00 + (tt(1:n) / (bc(1:n) + eps))**2
|
||||
wr(1:n) = 1.0d+00 + (tt(1:n) / (br(1:n) + eps))**2
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine stencil_weights_yc
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! function LIMITER_ZERO:
|
||||
! ---------------------
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user