From 76c7059a372c1d2f8b94e2f56c57e86b6410ae67 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 25 Feb 2023 15:59:17 -0300 Subject: [PATCH 1/6] INTERPOLATIONS: Update coefficients for OCMP5 method. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index ab77d7a..7d208fd 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -4283,7 +4283,7 @@ module interpolations ! The method uses an implicit tri-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 6.4e-12. +! discretization E to be 6.017e-12 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4316,8 +4316,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 5.0148583220385203994203557d-01 - real(kind=8), parameter :: a2 = 2.5452055302783787784904423d-01 + real(kind=8), parameter :: a1 = 5.01241562476900742404893381352899d-01 + real(kind=8), parameter :: a2 = 2.54639102475156530564823003088612d-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, & From 89a2a2b3c56500d56983bbeb10fcf51f205cd0ba Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 25 Feb 2023 16:01:41 -0300 Subject: [PATCH 2/6] INTERPOLATIONS: Update coefficients for OCMP7P method. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 7d208fd..d55d1b3 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -4679,7 +4679,7 @@ module interpolations ! The method uses an implicit penta-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.38e-11. +! discretization E to be 1.2984e-11 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4703,8 +4703,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 6.7100477315625665203995935d-01 - real(kind=8), parameter :: a2 = 3.4097048071278782662422942d-01 + real(kind=8), parameter :: a1 = 6.70862177813208397304414380169977d-01 + real(kind=8), parameter :: a2 = 3.41040434946989987878686778898846d-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, & From 42b40c5238239e9b6a0186218a801deaa1010aad Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 25 Feb 2023 16:03:37 -0300 Subject: [PATCH 3/6] INTERPOLATIONS: Update coefficients for OCMP9P method. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index d55d1b3..ef65413 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -4815,7 +4815,7 @@ module interpolations ! The method uses an implicit penta-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.54e-11. +! discretization E to be 1.4776e-11 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4839,8 +4839,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 6.7028117788580347100541103d-01 - real(kind=8), parameter :: a2 = 4.1002518667894543684869155d-01 + real(kind=8), parameter :: a1 = 6.69938726201951413971659371944416d-01 + real(kind=8), parameter :: a2 = 4.10221751253870603082971042746238d-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, & From d1a386fcf7913ea1f82fa768a6c576ba40cec11a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 25 Feb 2023 17:05:50 -0300 Subject: [PATCH 4/6] INTERPOLATIONS: Update coefficients for OCMP7T method. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index ef65413..8695d36 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -4415,7 +4415,7 @@ module interpolations ! The method uses an implicit tri-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.0e-11. +! discretization E to be 9.4859e-12 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4439,8 +4439,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 5.0119544102147968726471782d-01 - real(kind=8), parameter :: a2 = 3.0649545768561641522180678d-01 + real(kind=8), parameter :: a1 = 5.00765434358710549355767793738926d-01 + real(kind=8), parameter :: a2 = 3.06739643875254943361647028328465d-01 real(kind=8), dimension(3), parameter :: di7 = [ a1, 1.0d+00, a2 ] real(kind=8), dimension(7), parameter :: & ci7 = [ 4.00d+00 * a1 + 4.00d+00 * a2 - 3.00d+00, & From 5b26a7f7266badd492f48ce6ff4745e61b78f215 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 25 Feb 2023 18:01:53 -0300 Subject: [PATCH 5/6] INTERPOLATIONS: Update coefficients for OCMP9T method. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 8695d36..4393c46 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -4544,7 +4544,7 @@ module interpolations ! The method uses an implicit tri-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.32e-11. +! discretization E to be 1.2513e-11 (at local minimum of E/r²). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4568,8 +4568,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 5.0087697436175480833840357d-01 - real(kind=8), parameter :: a2 = 3.4091832657841234002695365d-01 + real(kind=8), parameter :: a1 = 5.00415305771646324666320357024935d-01 + real(kind=8), parameter :: a2 = 3.41203006456869647299948174549092d-01 real(kind=8), dimension(3), parameter :: di9 = [ a1, 1.0d+00, a2 ] real(kind=8), dimension(9), parameter :: & ci9 = [ - 5.000d+00 * a1 - 5.000d+00 * a2 + 4.000d+00, & From 99d5196736b8cbc71fb1cc7b40460b78f8d979d7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 27 Feb 2023 10:04:06 -0300 Subject: [PATCH 6/6] FORCING: Allow to specify the period of turbulence injection. The period is controlled by two parameters: 'injection_time_start' - the beginning of the injection period, 'injection_time_stop' - the end of the injection period. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 2 +- sources/forcing.F90 | 51 +++++++++++++++++++++++++++++++++---------- 2 files changed, 41 insertions(+), 12 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index c953f32..8b7fccc 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1109,7 +1109,7 @@ module evolution ! add forcing contribution ! - if (forcing_enabled) call update_forcing(dt) + if (forcing_enabled) call update_forcing(time, dt) ! check if we need to perform the refinement step ! diff --git a/sources/forcing.F90 b/sources/forcing.F90 index 5a5858b..af248b5 100644 --- a/sources/forcing.F90 +++ b/sources/forcing.F90 @@ -35,8 +35,8 @@ module forcing ! interfaces for procedure pointers ! abstract interface - subroutine update_forcing_iface(dt) - real(kind=8), intent(in) :: dt + subroutine update_forcing_iface(t, dt) + real(kind=8), intent(in) :: t, dt end subroutine end interface @@ -89,6 +89,8 @@ module forcing real(kind=8), save :: zth = 0.0d+00 real(kind=8), save :: zfc = 1.0d+00 #endif /* NDIMS == 3 */ + real(kind=8) :: tinj_start = 0.00000d+00 + real(kind=8) :: tinj_stop = 9.99999d+99 ! Fourier profile parameters ! @@ -244,6 +246,21 @@ module forcing end if #endif /* NDIMS == 3 */ +! get the injection period +! + call get_parameter("injection_time_start", tinj_start) + call get_parameter("injection_time_stop" , tinj_stop ) + if (tinj_start >= tinj_stop) then + if (verbose) then + write(*,*) + write(*,wfmt) "Injection cannot stop before it starts!" + write(*,wfmt) "Please, modify parameter 'injection_time_start'" // & + " or 'injection_time_stop'." + write(*,wfmt) "Resetting 'injection_time_start' to 0.0!" + tinj_start = 0.0d+00 + end if + end if + ! select the energy injection method ! select case(trim(injection)) @@ -929,6 +946,9 @@ module forcing call print_parameter(verbose, "amplitude threshold", fmin) end if + call print_parameter(verbose, "injection starts at", tinj_start) + call print_parameter(verbose, "injection stops at", tinj_stop) + if (.not. periodic(1)) & call print_parameter(verbose, "profile X transition between", xlo, xup) if (.not. periodic(2)) & @@ -953,18 +973,19 @@ module forcing ! ! Arguments: ! -! dt - time increment; +! t - time; +! dt - time step; ! !=============================================================================== ! - subroutine update_forcing_eddies(dt) + subroutine update_forcing_eddies(t, dt) use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen use random , only : randuni, randsym implicit none - real(kind=8), intent(in) :: dt + real(kind=8), intent(in) :: t, dt integer :: ni, n #if NDIMS == 3 @@ -974,6 +995,8 @@ module forcing !------------------------------------------------------------------------------- ! + if (t < tinj_start .or. t > tinj_stop) return + ! calculate the number of eddies to be injected during this interval ! dinj = dinj + rate * dt @@ -1029,11 +1052,12 @@ module forcing ! ! Arguments: ! -! dt - time increment; +! t - time; +! dt - time step; ! !=============================================================================== ! - subroutine update_forcing_oh(dt) + subroutine update_forcing_oh(t, dt) #if NDIMS == 3 use constants, only : pi2 @@ -1042,7 +1066,7 @@ module forcing implicit none - real(kind=8), intent(in) :: dt + real(kind=8), intent(in) :: t, dt integer :: l real(kind=8) :: acoeff, dcoeff @@ -1055,6 +1079,8 @@ module forcing !------------------------------------------------------------------------------- ! + if (t < tinj_start .or. t > tinj_stop) return + ! calculate drift and diffusion coefficients ! acoeff = exp(- dt / tscale) @@ -1113,18 +1139,19 @@ module forcing ! ! Arguments: ! -! dt - time increment; +! t - time; +! dt - time step; ! !=============================================================================== ! - subroutine update_forcing_alvelius(dt) + subroutine update_forcing_alvelius(t, dt) use constants, only : pi2 use random , only : randuni implicit none - real(kind=8), intent(in) :: dt + real(kind=8), intent(in) :: t, dt integer :: l real(kind=8) :: th1, dinj, sqdt @@ -1139,6 +1166,8 @@ module forcing !------------------------------------------------------------------------------- ! + if (t < tinj_start .or. t > tinj_stop) return + ! determine velocify coefficients in Fourier space ! call get_vcoefs()