diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 97f01ac..c04bad5 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -138,6 +138,15 @@ module interpolations real(kind=8), dimension(2), parameter :: & 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 ! private @@ -222,6 +231,15 @@ module interpolations call get_parameter("central_weight" , cweight ) 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 - ν) / ν ! kappa = min(kappa, (1.0d+00 - cfl) / cfl) @@ -347,18 +365,60 @@ module interpolations reconstruct_states => reconstruct_ocmp5 order = 5 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") name_rec = "7th order Optimized Compact Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_ocmp7 order = 7 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") name_rec = "9th order Optimized Compact Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_ocmp9 order = 9 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") write(stmp, '(f16.1)') sgp write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp & @@ -649,6 +709,13 @@ module interpolations call print_section(verbose, "Interpolations") 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, "prolongation limiter", name_plim) if (mlp) then @@ -4295,17 +4362,6 @@ module interpolations 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) @@ -4423,19 +4479,6 @@ module interpolations 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) @@ -4565,21 +4608,6 @@ module interpolations 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)