From 222b08e4a12efd5c79d425c6ece79f5d16f01fbb Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 27 Oct 2020 10:49:06 +0100
Subject: [PATCH] kr=0 symmetric initialization and enforcement at each time
step
---
src/compute_Sapj.F90 | 77 ++++++++---------------------
src/diagnose.F90 | 1 +
src/inital.F90 | 20 +++++---
src/poisson.F90 | 75 ++++++++++------------------
src/stepon.F90 | 23 +++++++++
src/time_integration_mod.F90 | 95 ++++++++++++++++++++++++++++++------
6 files changed, 162 insertions(+), 129 deletions(-)
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index be2e7f4f..6618ce0c 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -24,32 +24,32 @@ SUBROUTINE compute_Sapj
sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln argument
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
- ! ploope: DO ip = 1,1
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
- ! jloope: DO ij = 1,1
Sepj(ip,ij,:,:) = 0._dp
factn = 1
nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
- ! nloope: DO in = 1,1
krloope1: DO ikr = 1,Nkr ! Loop over kr
kzloope1: DO ikz = 1,Nkz ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
- IF ( DK ) THEN
- be_2 = 0._dp
- ELSE
- be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
- ENDIF
+ be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
kerneln = be_2**(in-1)/factn * EXP(-be_2)
-
! First convolution term
IF ( NON_LIN ) THEN
F_(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
ELSE
F_(ikr,ikz) = 0._dp
ENDIF
+ ! Background sinusoidal electrostatic potential phi_0 for KH inst.
+ IF ( q_e .NE. 0._dp ) THEN ! If electron are not removed
+ IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
+ IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0
+ F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz)
+ ENDIF
+ ENDIF
+ ENDIF
! Second convolution term
G_(ikr,ikz) = 0._dp ! initialization of the sum
DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
@@ -57,14 +57,6 @@ SUBROUTINE compute_Sapj
dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
ENDDO
G_(ikr,ikz) = imagu*kz*G_(ikr,ikz)
- ! G_(ikr,ikz) = imagu*kz* moments_e(ip,ij,ikr,ikz,updatetlevel)
-
- ! Background sinusoidal electrostatic potential phi_0 for KH inst.
- IF ( (kz .EQ. 0._dp) .AND. (q_e .NE. 0._dp) ) THEN ! Kronecker kz=0
- IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0
- F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz)
- ENDIF
- ENDIF
ENDDO kzloope1
ENDDO krloope1
@@ -72,23 +64,15 @@ SUBROUTINE compute_Sapj
CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve Fourier to Fourier
Sepj(ip,ij,:,:) = Sepj(ip,ij,:,:) + CONV ! Add it to Sepj
- IF ( NON_LIN ) THEN
+ IF ( NON_LIN ) THEN ! Fully non linear term ikz phi * ikr Napj
krloope2: DO ikr = 1,Nkr ! Loop over kr
kzloope2: DO ikz = 1,Nkz ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
-
- IF ( DK ) THEN
- be_2 = 0._dp
- ELSE
- be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
- ENDIF
-
+ be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
kerneln = be_2**(in-1)/factn * EXP(-be_2)
-
! First convolution term
- F_(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
-
+ F_(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
! Second convolution term
G_(ikr,ikz) = 0._dp ! initialization of the sum
DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
@@ -117,26 +101,17 @@ SUBROUTINE compute_Sapj
sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
- ! ploopi: DO ip = 1,1 ! Loop over Hermite moments
jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
- ! jloopi: DO ij = 1,1 ! Loop over Laguerre moments
Sipj(ip,ij,:,:) = 0._dp
factn = 1
nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
- ! nloopi: DO in = 1,1 ! Loop over laguerre for the sum
krloopi1: DO ikr = 1,Nkr ! Loop over kr
kzloopi1: DO ikz = 1,Nkz ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
-
- IF ( DK ) THEN
- bi_2 = 0._dp
- ELSE
- bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
- ENDIF
-
+ bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
! First convolution term
@@ -145,6 +120,12 @@ SUBROUTINE compute_Sapj
ELSE
F_(ikr,ikz) = 0._dp
ENDIF
+ ! Background sinusoidal electrostatic potential phi_0 for KH inst.
+ IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
+ IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0
+ F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz)
+ ENDIF
+ ENDIF
! Second convolution term
G_(ikr,ikz) = 0._dp ! initialization of the sum
DO is = 1, MIN( in+ij-1, jmaxi+1 )
@@ -152,14 +133,6 @@ SUBROUTINE compute_Sapj
dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
ENDDO
G_(ikr,ikz) = imagu*kz*G_(ikr,ikz)
- ! G_(ikr,ikz) = imagu*kz* moments_i(ip,ij,ikr,ikz,updatetlevel)
-
- ! Background sinusoidal electrostatic potential phi_0 for KH inst.
- IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
- IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0
- F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz)
- ENDIF
- ENDIF
ENDDO kzloopi1
ENDDO krloopi1
@@ -171,18 +144,10 @@ SUBROUTINE compute_Sapj
kzloopi2: DO ikz = 1,Nkz ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
-
- IF ( DK ) THEN
- bi_2 = 0._dp
- ELSE
- bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
- ENDIF
-
+ bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
-
! First convolution term
F_(ikr,ikz) = imagu*kz*phi(ikr,ikz) * kerneln
-
! Second convolution term
G_(ikr,ikz) = 0._dp ! initialization of the sum
IF ( NON_LIN ) THEN
@@ -192,7 +157,6 @@ SUBROUTINE compute_Sapj
ENDDO
G_(ikr,ikz) = imagu*kr*G_(ikr,ikz)
ENDIF
-
! Background sinusoidal ion denstiy n_i0 for KH inst.
IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=+-kr0
@@ -209,6 +173,7 @@ SUBROUTINE compute_Sapj
IF ( in + 1 .LE. jmaxi+1 ) THEN
factn = real(in,dp) * factn ! factorial(n+1)
ENDIF
+
ENDDO nloopi
ENDDO jloopi
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index bba9f425..82fa6c8c 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -219,6 +219,7 @@ SUBROUTINE diagnose(kstep)
ELSEIF (kstep .EQ. -1) THEN
CALL cpu_time(finish)
+ CALL attach(fidres, "/data/input","cpu_time",finish-start)
! Close all diagnostic files
CALL closef(fidres)
END IF
diff --git a/src/inital.F90 b/src/inital.F90
index 822b5dbb..b3d9fc9b 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -89,12 +89,14 @@ SUBROUTINE init_moments
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
- moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- ! moments_e( ip,ij, Nkz-ikz, ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
- ! moments_e( ip,ij, ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
- ! moments_e( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) &
+ * AA_r(ikr) * AA_z(ikz)
END DO
END DO
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ CALL RANDOM_NUMBER(noise)
+ moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :)
+ END DO
END DO
END DO
DO ip=ips_i,ipe_i
@@ -102,12 +104,14 @@ SUBROUTINE init_moments
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
- moments_i( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- ! moments_i( ip,ij, Nkz-ikz, ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
- ! moments_i( ip,ij, ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
- ! moments_i( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) &
+ * AA_r(ikr) * AA_z(ikz)
END DO
END DO
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ CALL RANDOM_NUMBER(noise)
+ moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :)
+ END DO
END DO
END DO
diff --git a/src/poisson.F90 b/src/poisson.F90
index 305a6b32..9af24e7a 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -40,46 +40,35 @@ SUBROUTINE poisson
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
!! Sum(kernel * Moments_0n)
- ! IF ( DK ) THEN
- ! sum_kernel_mom_e = moments_e(1, 1, ikr, ikz, updatetlevel)
- ! ELSE
- ! Initialization for n = 0 (ine = 1)
- Kne = exp(-be_2)
- sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
- ! loop over n only if the max polynomial degree is 1 or more
- if (jmaxe .GT. 0) then
- DO ine=2,jmaxe+1 ! ine = n+1
- ine_dp = REAL(ine-1,dp)
- ! We update iteratively the kernel functions (to spare factorial computations)
- Kne = Kne * be_2/ine_dp
- sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
- END DO
- endif
- ! ENDIF
- !! Sum(kernel**2)
- ! IF ( DK ) THEN
- ! sum_kernel2_e = 1._dp - 2._dp*be_2 + be_2**2!(1._dp-be_2)**2
- ! ELSE
- ! Initialization for n = 0 (ine = 1)
- Kne = exp(-be_2)
- sum_kernel2_e = Kne**2
- ! loop over n only without caring of jmax since no moment dependency
- ! DO ine=2,10
- if (jmaxe .GT. 0) then
- DO ine=2,jmaxe+1 ! ine = n+1
- ine_dp = REAL(ine-1,dp) ! Real index (0 to jmax)
- Kne = Kne * be_2/ine_dp ! update kernel_n
- sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
- END DO
- ENDIF
+ ! Initialization for n = 0 (ine = 1)
+ Kne = exp(-be_2)
+ sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
+ ! loop over n only if the max polynomial degree is 1 or more
+ if (jmaxe .GT. 0) then
+ DO ine=2,jmaxe+1 ! ine = n+1
+ ine_dp = REAL(ine-1,dp)
+ ! We update iteratively the kernel functions (to spare factorial computations)
+ Kne = Kne * be_2/ine_dp
+ sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
+ END DO
+ endif
+ ! Initialization for n = 0 (ine = 1)
+ Kne = exp(-be_2)
+ sum_kernel2_e = Kne**2
+ ! loop over n only without caring of jmax since no moment dependency
+ ! DO ine=2,10
+ if (jmaxe .GT. 0) then
+ DO ine=2,jmaxe+1 ! ine = n+1
+ ine_dp = REAL(ine-1,dp) ! Real index (0 to jmax)
+ Kne = Kne * be_2/ine_dp ! update kernel_n
+ sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
+ END DO
+ ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
!! Sum(kernel * Moments_0n)
! Initialization for n = 0 (ini = 1)
- ! IF ( DK ) THEN
- ! sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel)
- ! ELSE
Kni = exp(-bi_2)
sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel)
! loop over n only if the max polynomial degree is 1 or more
@@ -91,11 +80,7 @@ SUBROUTINE poisson
sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
END DO
endif
- ! ENDIF
- !! Sum(kernel**2)
- ! IF ( DK ) THEN
- ! sum_kernel2_i = 1._dp - 2._dp*bi_2 + bi_2**2 !(1._dp-bi_2)**2
- ! ELSE
+
! Initialization for n = 0 (ini = 1)
Kni = exp(-bi_2)
sum_kernel2_i = Kni**2
@@ -113,10 +98,6 @@ SUBROUTINE poisson
alphaD = kperp2 * lambdaD**2
gammaD = alphaD + qe2_taue * (1._dp - sum_kernel2_e) & ! Called Poisson_ in MOLI
+ qi2_taui * (1._dp - sum_kernel2_i)
- ! ! Adiabatic electron response
- ! IF ( q_e .EQ. 0 ) THEN
- ! gammaD = gammaD + 1._dp
- ! ENDIF
gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
@@ -125,12 +106,6 @@ SUBROUTINE poisson
ENDIF
phi(ikr, ikz) = gammaD_phi/gammaD
- ! write(*,*) -2._dp*gammaD,' ',-kperp2
- ! phi(ikr, ikz) = -gammaD_phi/kperp2
- ! write(*,*) 'gD = ', gammaD,', kperp = ',-kperp2
- ! write(*,*) 'sum_kernel2_i = ', sum_kernel2_i, ', sum_kernel2_e', sum_kernel2_e
- ! phi(ikr, ikz) = gammaD_phi/(qi2_taui * (sum_kernel2_i))
- ! phi(ikr, ikz) = -moments_i(1,1,ikr,ikz,1)/kperp2
END DO
END DO
diff --git a/src/stepon.F90 b/src/stepon.F90
index eb55a104..a900af8f 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -40,6 +40,8 @@ SUBROUTINE stepon
ENDIF
!(The two routines above are called in inital for t=0)
+ CALL enforce_symetry() ! Enforcing symmetry on kr = 0
+
CALL checkfield_all()
END DO
@@ -62,4 +64,25 @@ SUBROUTINE stepon
ENDIF
END SUBROUTINE checkfield_all
+ SUBROUTINE enforce_symetry
+
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ moments_e( ip,ij,1,ikz, :) = CONJG(moments_e( ip,ij,1,Nkz+2-ikz, :))
+ END DO
+ END DO
+ END DO
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ moments_i( ip,ij,1,ikz, :) = CONJG(moments_i( ip,ij,1,Nkz+2-ikz, :))
+ END DO
+ END DO
+ END DO
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ phi(1,ikz) = phi(1,Nkz+2-ikz)
+ END DO
+ END SUBROUTINE enforce_symetry
+
END SUBROUTINE stepon
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 4936e6f9..3f473707 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -8,8 +8,8 @@ MODULE time_integration
INTEGER, PUBLIC, PROTECTED :: updatetlevel ! Current time level to be updated
real(dp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E,A_I
- real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_I
- real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
+ real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_Es,b_I
+ real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
character(len=10),PUBLIC,PROTECTED :: numerical_scheme='RK4'
@@ -47,9 +47,9 @@ CONTAINS
IMPLICIT NONE
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
-
+
CALL attach(fidres, TRIM(str), "numerical_scheme", numerical_scheme)
-
+
END SUBROUTINE time_integration_outputinputs
@@ -57,7 +57,7 @@ CONTAINS
SUBROUTINE set_numerical_scheme
! Initialize Butcher coefficient of numerical_scheme
- use basic
+ use basic
IMPLICIT NONE
SELECT CASE (numerical_scheme)
@@ -65,7 +65,9 @@ CONTAINS
CALL RK4
CASE ('DOPRI5')
CALL DOPRI5
- CASE DEFAULT
+ CASE ('DOPRI5_ADAPT')
+ CALL DOPRI5_ADAPT
+ CASE DEFAULT
WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
END SELECT
@@ -82,20 +84,20 @@ CONTAINS
INTEGER,PARAMETER :: nbstep = 4
CALL allocate_array_dp1(c_E,1,nbstep)
- CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 4
- c_E(1) = 0.0_dp
+ c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp
c_E(3) = 1.0_dp/2.0_dp
c_E(4) = 1.0_dp
b_E(1) = 1.0_dp/6.0_dp
- b_E(2) = 1.0_dp/3.0_dp
+ b_E(2) = 1.0_dp/3.0_dp
b_E(3) = 1.0_dp/3.0_dp
- b_E(4) = 1.0_dp/6.0_dp
+ b_E(4) = 1.0_dp/6.0_dp
A_E(2,1) = 1.0_dp/2.0_dp
A_E(3,2) = 1.0_dp/2.0_dp
@@ -108,22 +110,22 @@ CONTAINS
SUBROUTINE DOPRI5
! Butcher coeff for DOPRI5 --> Stiffness detection
! DOPRI5 used for stiffness detection.
- ! 5 order method/7 stages
+ ! 5 order method/7 stages
USE basic
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep =7
CALL allocate_array_dp1(c_E,1,nbstep)
- CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 7
!
c_E(1) = 0._dp
- c_E(2) = 1.0_dp/5.0_dp
- c_E(3) = 3.0_dp /10.0_dp
- c_E(4) = 4.0_dp/5.0_dp
+ c_E(2) = 1.0_dp/5.0_dp
+ c_E(3) = 3.0_dp /10.0_dp
+ c_E(4) = 4.0_dp/5.0_dp
c_E(5) = 8.0_dp/9.0_dp
c_E(6) = 1.0_dp
c_E(7) = 1.0_dp
@@ -159,4 +161,67 @@ CONTAINS
!
END SUBROUTINE DOPRI5
+ SUBROUTINE DOPRI5_ADAPT
+ ! Butcher coeff for DOPRI5 --> Stiffness detection
+ ! DOPRI5 used for stiffness detection.
+ ! 5 order method/7 stages
+
+ USE basic
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep =7
+
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp1(b_Es,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+
+ ntimelevel = 7
+ !
+ c_E(1) = 0._dp
+ c_E(2) = 1.0_dp/5.0_dp
+ c_E(3) = 3.0_dp /10.0_dp
+ c_E(4) = 4.0_dp/5.0_dp
+ c_E(5) = 8.0_dp/9.0_dp
+ c_E(6) = 1.0_dp
+ c_E(7) = 1.0_dp
+ !
+ A_E(2,1) = 1.0_dp/5.0_dp
+ A_E(3,1) = 3.0_dp/40.0_dp
+ A_E(3,2) = 9.0_dp/40.0_dp
+ A_E(4,1) = 44.0_dp/45.0_dp
+ A_E(4,2) = -56.0_dp/15.0_dp
+ A_E(4,3) = 32.0_dp/9.0_dp
+ A_E(5,1 ) = 19372.0_dp/6561.0_dp
+ A_E(5,2) = -25360.0_dp/2187.0_dp
+ A_E(5,3) = 64448.0_dp/6561.0_dp
+ A_E(5,4) = -212.0_dp/729.0_dp
+ A_E(6,1) = 9017.0_dp/3168.0_dp
+ A_E(6,2)= -355.0_dp/33.0_dp
+ A_E(6,3) = 46732.0_dp/5247.0_dp
+ A_E(6,4) = 49.0_dp/176.0_dp
+ A_E(6,5) = -5103.0_dp/18656.0_dp
+ A_E(7,1) = 35.0_dp/384.0_dp
+ A_E(7,3) = 500.0_dp/1113.0_dp
+ A_E(7,4) = 125.0_dp/192.0_dp
+ A_E(7,5) = -2187.0_dp/6784.0_dp
+ A_E(7,6) = 11.0_dp/84.0_dp
+ !
+ b_E(1) = 35.0_dp/384.0_dp
+ b_E(2) = 0._dp
+ b_E(3) = 500.0_dp/1113.0_dp
+ b_E(4) = 125.0_dp/192.0_dp
+ b_E(5) = -2187.0_dp/6784.0_dp
+ b_E(6) = 11.0_dp/84.0_dp
+ b_E(7) = 0._dp
+ !
+ b_Es(1) = 5179.0_dp/57600.0_dp
+ b_Es(2) = 0._dp
+ b_Es(3) = 7571.0_dp/16695.0_dp
+ b_Es(4) = 393.0_dp/640.0_dp
+ b_Es(5) = -92097.0_dp/339200.0_dp
+ b_Es(6) = 187.0_dp/2100.0_dp
+ b_Es(7) = 1._dp/40._dp
+ !
+ END SUBROUTINE DOPRI5_ADAPT
+
END MODULE time_integration
--
GitLab