Merge branch 'master' into reconnection

This commit is contained in:
Grzegorz Kowal 2023-02-08 12:34:37 -03:00
commit bccd0501cc

View File

@ -138,6 +138,15 @@ module interpolations
real(kind=8), dimension(2), parameter :: & real(kind=8), dimension(2), parameter :: &
ce2 = [ 5.0d-01, 5.0d-01 ] ce2 = [ 5.0d-01, 5.0d-01 ]
! coefficients for the optimized compact schemes
!
logical :: oc_scheme = .false.
logical :: oc_stable = .false.
real(kind=8), dimension(2) :: oc_a = 0.0d+00
real(kind=8), dimension(3) :: di5
real(kind=8), dimension(5) :: di7, di9, ci5, ci7
real(kind=8), dimension(7) :: ci9
! by default everything is private ! by default everything is private
! !
private private
@ -222,6 +231,15 @@ module interpolations
call get_parameter("central_weight" , cweight ) call get_parameter("central_weight" , cweight )
call get_parameter("cfl" , cfl ) call get_parameter("cfl" , cfl )
stmp = 'minimum'
call get_parameter("ocmp_scheme_mode", stmp)
select case(trim(stmp))
case ("stable", "STABLE", "stability", "STABILITY")
oc_stable = .true.
case default
oc_stable = .false.
end select
! calculate κ = (1 - ν) / ν ! calculate κ = (1 - ν) / ν
! !
kappa = min(kappa, (1.0d+00 - cfl) / cfl) kappa = min(kappa, (1.0d+00 - cfl) / cfl)
@ -347,18 +365,60 @@ module interpolations
reconstruct_states => reconstruct_ocmp5 reconstruct_states => reconstruct_ocmp5
order = 5 order = 5
nghosts = max(nghosts, 4) nghosts = max(nghosts, 4)
oc_scheme = .true.
if (oc_stable) then
oc_a = [ 4.98486328125000044409d-01, 2.50756835937499977796d-01 ]
else
oc_a = [ 5.01779599466318559919d-01, 2.54377898581479855444d-01 ]
end if
di5 = [ oc_a(1), 1.0d+00, oc_a(2) ]
ci5 = [- 3.0d+00 * oc_a(1) - 3.0d+00 * oc_a(2) + 2.0d+00, &
2.7d+01 * oc_a(1) + 1.7d+01 * oc_a(2) - 1.3d+01, &
4.7d+01 * oc_a(1) - 4.3d+01 * oc_a(2) + 4.7d+01, &
- 1.3d+01 * oc_a(1) + 7.7d+01 * oc_a(2) + 2.7d+01, &
2.0d+00 * oc_a(1) + 1.2d+01 * oc_a(2) - 3.0d+00 ] / 6.0d+01
case ("ocmp7", "OCMP7") case ("ocmp7", "OCMP7")
name_rec = "7th order Optimized Compact Monotonicity Preserving" name_rec = "7th order Optimized Compact Monotonicity Preserving"
interfaces => interfaces_dir interfaces => interfaces_dir
reconstruct_states => reconstruct_ocmp7 reconstruct_states => reconstruct_ocmp7
order = 7 order = 7
nghosts = max(nghosts, 6) nghosts = max(nghosts, 6)
oc_scheme = .true.
if (oc_stable) then
oc_a = [ 6.64390055338541563046d-01, 3.34471638997395903647d-01 ]
else
oc_a = [ 6.71157762607103580699d-01, 3.40895359923885754583d-01 ]
end if
di7 = [ (3.0d+00 * oc_a(1) + 2.0d+00 * oc_a(2) - 2.0d+00) / 8.0d+00, &
oc_a(1), 1.0d+00, oc_a(2), &
( oc_a(1) + 6.0d+00 * oc_a(2) - 2.0d+00) / 4.0d+01 ]
ci7 = [ 1.80d+01 * oc_a(1) + 1.80d+01 * oc_a(2) - 1.60d+01, &
5.43d+02 * oc_a(1) + 2.68d+02 * oc_a(2) - 2.91d+02, &
3.43d+02 * oc_a(1) - 3.32d+02 * oc_a(2) + 5.09d+02, &
-1.07d+02 * oc_a(1) + 5.68d+02 * oc_a(2) + 3.09d+02, &
4.30d+01 * oc_a(1) + 3.18d+02 * oc_a(2) - 9.10d+01 ] / 6.0d+02
case ("ocmp9", "OCMP9") case ("ocmp9", "OCMP9")
name_rec = "9th order Optimized Compact Monotonicity Preserving" name_rec = "9th order Optimized Compact Monotonicity Preserving"
interfaces => interfaces_dir interfaces => interfaces_dir
reconstruct_states => reconstruct_ocmp9 reconstruct_states => reconstruct_ocmp9
order = 9 order = 9
nghosts = max(nghosts, 6) nghosts = max(nghosts, 6)
oc_scheme = .true.
if (oc_stable) then
oc_a = [ 6.64322884877522779057d-01, 4.01406269073486310361d-01 ]
else
oc_a = [ 6.70794426721399439373d-01, 4.09731006175920509094d-01 ]
end if
di9 = [ (9.0d+00 * oc_a(1) + 5.0d+00 * oc_a(2) - 6.0d+00) / 2.0d+01, &
oc_a(1), 1.0d+00, oc_a(2), &
( oc_a(1) + 5.0d+00 * oc_a(2) - 2.0d+00) / 2.0d+01 ]
ci9 = [-1.000d+01 * oc_a(1) -1.000d+01 * oc_a(2) +1.000d+01, &
2.420d+02 * oc_a(1) +2.000d+02 * oc_a(2) -2.140d+02, &
4.085d+03 * oc_a(1) +1.635d+03 * oc_a(2) -2.146d+03, &
2.335d+03 * oc_a(1) -1.865d+03 * oc_a(2) +3.454d+03, &
-8.150d+02 * oc_a(1) +3.385d+03 * oc_a(2) +2.404d+03, &
4.450d+02 * oc_a(1) +2.895d+03 * oc_a(2) -9.560d+02, &
1.800d+01 * oc_a(1) +6.000d+01 * oc_a(2) -3.200d+01 ] / 4.2d+03
case ("gp", "GP") case ("gp", "GP")
write(stmp, '(f16.1)') sgp write(stmp, '(f16.1)') sgp
write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp & write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp &
@ -649,6 +709,13 @@ module interpolations
call print_section(verbose, "Interpolations") call print_section(verbose, "Interpolations")
call print_parameter(verbose, "reconstruction" , name_rec ) call print_parameter(verbose, "reconstruction" , name_rec )
if (oc_scheme) then
if (oc_stable) then
call print_parameter(verbose, "scheme optimized for", "stability")
else
call print_parameter(verbose, "scheme optimized for", "minimum error")
end if
end if
call print_parameter(verbose, "TVD limiter" , name_tlim) call print_parameter(verbose, "TVD limiter" , name_tlim)
call print_parameter(verbose, "prolongation limiter", name_plim) call print_parameter(verbose, "prolongation limiter", name_plim)
if (mlp) then if (mlp) then
@ -4295,17 +4362,6 @@ module interpolations
real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u 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) n = size(fc)
@ -4423,19 +4479,6 @@ module interpolations
real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u 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) n = size(fc)
@ -4565,21 +4608,6 @@ module interpolations
real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u 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) n = size(fc)