Skip to content

Commit 03bba8b

Browse files
DhueperDhueper
committed
Adaptive time step
1 parent be6d065 commit 03bba8b

10 files changed

Lines changed: 595 additions & 88 deletions

‎Solver/src/libs/mesh/StorageClass.f90‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,6 +79,7 @@ module StorageClass
7979
real(kind=RP), allocatable :: S_NSP(:,:,:,:) ! NSE Particles source term
8080
real(kind=RP), allocatable :: QbaseSponge(:,:,:,:) ! Base Flow State vector for sponges
8181
real(kind=RP), allocatable :: intensitySponge(:,:,:) ! Intensity for sponges
82+
real(kind=RP), allocatable :: QLowRK(:,:,:,:) ! NSE State vector solved with a lower order RK method (required for adaptive time-stepping)
8283
#ifndef ACOUSTIC
8384
real(kind=RP), allocatable :: mu_NS(:,:,:,:) ! (mu, beta, kappa) artificial
8485
real(kind=RP), allocatable :: mu_turb_NS(:,:,:) ! mu of LES
@@ -791,6 +792,7 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
791792
!
792793
#ifdef FLOW
793794
allocate ( self % QNS (1:NCONS,0:Nx,0:Ny,0:Nz) )
795+
allocate ( self % QLowRK(1:NCONS,0:Nx,0:Ny,0:Nz) )
794796
allocate (self % QbaseSponge(1:NCONS,0:Nx,0:Ny,0:Nz) )
795797
allocate (self % intensitySponge(0:Nx,0:Ny,0:Nz) )
796798
allocate ( self % QdotNS(1:NCONS,0:Nx,0:Ny,0:Nz) )
@@ -877,6 +879,7 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
877879
self % S_NS = 0.0_RP
878880
self % S_NSP = 0.0_RP
879881
self % QNS = 0.0_RP
882+
self % QLowRK = 0.0_RP
880883
self % QbaseSponge = 0.0_RP
881884
self % intensitySponge = 0.0_RP
882885
self % QDotNS = 0.0_RP
@@ -980,6 +983,7 @@ elemental subroutine ElementStorage_Assign(to, from)
980983

981984
#ifdef FLOW
982985
to % QNS = from % QNS
986+
to % QLowRK = from % QLowRK
983987
to % QbaseSponge = from % QbaseSponge
984988
to % intensitySponge = from % intensitySponge
985989

@@ -1076,6 +1080,7 @@ elemental subroutine ElementStorage_Destruct(self)
10761080
safedeallocate(self % QDotNS)
10771081
safedeallocate(self % QbaseSponge)
10781082
safedeallocate(self % intensitySponge)
1083+
safedeallocate(self % QLowRK)
10791084

10801085
if ( allocated(self % PrevQ) ) then
10811086
num_prevSol = size(self % PrevQ)
@@ -1109,6 +1114,9 @@ elemental subroutine ElementStorage_Destruct(self)
11091114
!if (self % anJacobian) then ! Not needed since there's only one variable (= one if)
11101115
safedeallocate(self % dF_dgradQ)
11111116
!end if
1117+
1118+
! Destruct statistics
1119+
call self % stats % destruct()
11121120
#endif
11131121

11141122
#ifdef CAHNHILLIARD
@@ -1139,6 +1147,18 @@ elemental subroutine ElementStorage_Destruct(self)
11391147

11401148
safedeallocate(self % PrevQ)
11411149

1150+
! Deallocate RKSteps
1151+
! ------------------
1152+
#ifdef MULTIPHASE
1153+
if ( allocated(self % RKSteps) ) then
1154+
do k = 1, size(self % RKSteps)
1155+
safedeallocate(self % RKSteps(k) % K)
1156+
safedeallocate(self % RKSteps(k) % hatK)
1157+
end do
1158+
deallocate(self % RKSteps)
1159+
end if
1160+
#endif
1161+
11421162
end subroutine ElementStorage_Destruct
11431163
#ifdef FLOW
11441164
pure subroutine ElementStorage_SetStorageToNS(self)
@@ -1271,6 +1291,7 @@ impure elemental subroutine ElementStorage_InterpolateSolution(this,other,nodes,
12711291
if (all(this % Nxyz == other % Nxyz)) then
12721292
other % Q = this % Q
12731293
other % QbaseSponge = this % QbaseSponge
1294+
other % QLowRK = this % QLowRK
12741295
else
12751296
!$omp critical
12761297
call NodalStorage(this % Nxyz(1)) % construct(nodes,this % Nxyz(1))
@@ -1319,6 +1340,15 @@ impure elemental subroutine ElementStorage_InterpolateSolution(this,other,nodes,
13191340
this % Q(1:,0:,0:,0:) => this % QbaseSponge
13201341
other % Q(1:,0:,0:,0:) => other % QbaseSponge
13211342

1343+
call Interp3DArrays (Nvars = NCONS , &
1344+
Nin = this % Nxyz , &
1345+
inArray = this % Q , &
1346+
Nout = other % Nxyz , &
1347+
outArray = other % Q )
1348+
1349+
this % Q(1:,0:,0:,0:) => this % QLowRK
1350+
other % Q(1:,0:,0:,0:) => other % QLowRK
1351+
13221352
call Interp3DArrays (Nvars = NCONS , &
13231353
Nin = this % Nxyz , &
13241354
inArray = this % Q , &

‎Solver/src/libs/timeintegrator/AdaptiveTimeStepClass.f90‎

Lines changed: 280 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,280 @@
1+
!
2+
!//////////////////////////////////////////////////////
3+
!
4+
! Class cotaining routines for adapting time step based on Reinforcement Learning.
5+
! -> The adaptation procedure is performed with a RL agent trined with a Value Iteration algorithm.
6+
! -> The current implementation is compatible with OpenMP and MPI.
7+
!
8+
!////////////////////////////////////////////////////////////////////////
9+
!
10+
module AdaptiveTimeStepClass
11+
use SMConstants
12+
#ifdef NAVIERSTOKES
13+
use PhysicsStorage , only: CTD_IGNORE_MODE, flowIsNavierStokes
14+
use FluidData , only: thermodynamics
15+
#elif defined(MULTIPHASE)
16+
use PhysicsStorage , only: CTD_IGNORE_MODE, IMP
17+
#else
18+
use PhysicsStorage , only: CTD_IGNORE_MODE
19+
#endif
20+
use ElementClass
21+
use HexMeshClass
22+
USE DGSEMClass
23+
use FTValueDictionaryClass , only: FTValueDictionary
24+
use StorageClass
25+
use MPI_Process_Info
26+
27+
#ifdef _HAS_MPI_
28+
use mpi
29+
#endif
30+
implicit none
31+
32+
#include "Includes.h"
33+
private
34+
public adaptiveTimeStep_t
35+
36+
type :: adaptiveTimeStep_t
37+
real(kind=RP) :: min_dt = 1e-8_RP ! Minimum time step
38+
real(kind=RP) :: max_dt = 1e-2_RP ! Maximum time step
39+
real(kind=RP) :: b1 = 0.0 !PID 1st parameter
40+
real(kind=RP) :: b2 = -0.2 !PID 2nd parameter
41+
real(kind=RP) :: b3 = 0.2 !PID 3rd parameter
42+
real(kind=RP) :: k = 3.0 !Order of the RK method (only RK3 is implemented)
43+
real(kind=RP) :: dt_eps(3) !Epsilon for the PID controller
44+
logical :: constructed = .FALSE. !If the adaptive time step is constructed
45+
real(kind=RP) :: dtAdaptationStep = 1e-3_RP !Adaptation step
46+
real(kind=RP) :: lastAdaptationTime = 0.0_RP !Last adaptation time
47+
integer :: error_variable !1:u, 2:v, 3:w, 4:rho*u, 5:rho*v, 6:rho*w, 7:p (only for Navier-Stokes)
48+
logical :: pAdapted = .FALSE.
49+
50+
contains
51+
procedure :: construct => adaptiveTimeStep_Construct
52+
procedure :: destruct => adaptiveTimeStep_Destruct
53+
procedure :: hasToAdapt => adaptiveTimeStep_HasToAdapt
54+
procedure :: update => adaptiveTimeStep_Update
55+
end type adaptiveTimeStep_t
56+
57+
contains
58+
59+
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
60+
!
61+
! -----------------------------------------------
62+
! Routine for constructing the adaptive time step
63+
! -----------------------------------------------
64+
subroutine adaptiveTimeStep_Construct(this, controlVariables, t0)
65+
implicit none
66+
class(adaptiveTimeStep_t), intent(inout) :: this
67+
type(FTValueDictionary) , intent(in) :: controlVariables
68+
real(kind=RP) , intent(in) :: t0
69+
!! ---------------
70+
this % dt_eps(1) = 1.0_RP
71+
this % dt_eps(2) = 1.0_RP
72+
this % dt_eps(3) = 1.0_RP
73+
this % constructed = .TRUE.
74+
this % lastAdaptationTime = t0
75+
76+
if (controlVariables % containsKey("dt adaptation step")) then
77+
this % dtAdaptationStep = controlVariables % doublePrecisionValueForKey("dt adaptation step")
78+
end if
79+
80+
if (controlVariables % containsKey("minimum dt")) then
81+
this % min_dt = controlVariables % doublePrecisionValueForKey("minimum dt")
82+
end if
83+
84+
if (controlVariables % containsKey("maximum dt")) then
85+
this % max_dt = controlVariables % doublePrecisionValueForKey("maximum dt")
86+
end if
87+
end subroutine adaptiveTimeStep_Construct
88+
!
89+
!
90+
! -----------------------------------------------
91+
! Subroutine for destructing the adaptive time step
92+
! -----------------------------------------------
93+
subroutine adaptiveTimeStep_Destruct(this)
94+
implicit none
95+
class(adaptiveTimeStep_t), intent(inout) :: this
96+
!! ---------------
97+
end subroutine adaptiveTimeStep_Destruct
98+
!
99+
! -----------------------------------------------
100+
! Subroutine for updating the adaptive time step
101+
! -----------------------------------------------
102+
subroutine adaptiveTimeStep_Update(this, mesh, t, dt)
103+
implicit none
104+
class(adaptiveTimeStep_t) , intent(inout) :: this
105+
class(HexMesh) , intent(in) :: mesh
106+
real(kind=RP) , intent(in) :: t
107+
real(kind=RP) , intent(inout) :: dt
108+
!! ---------------
109+
! Local variables
110+
!! ---------------
111+
integer :: eID, ierr, i
112+
real(kind=RP) :: dt_weight, sum_dt_weight, dt_lim, min_dt_lim, max_dt_lim
113+
114+
min_dt_lim = 0.5_RP ! Minimum limit for the time step
115+
max_dt_lim = 1.3_RP ! Maximum limit for the time step
116+
117+
this % lastAdaptationTime = t
118+
dt_weight = 0.0_RP
119+
sum_dt_weight = 0.0_RP
120+
!$omp parallel shared(dt_weight, mesh, this)
121+
!$omp do reduction(+:dt_weight) schedule(runtime)
122+
do eID = 1, mesh % no_of_elements
123+
dt_weight = dt_weight + adaptiveTimeStep_ComputeWeights(this, mesh % elements(eID))
124+
end do
125+
!$omp end do
126+
!$omp end parallel
127+
128+
if ( MPI_Process % doMPIAction ) then
129+
#ifdef _HAS_MPI_
130+
call mpi_allreduce ( dt_weight, sum_dt_weight, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr )
131+
#endif
132+
else
133+
sum_dt_weight = dt_weight
134+
end if
135+
136+
if (isnan(sum_dt_weight)) then
137+
this % dt_eps(3) = 1e-10_RP
138+
else
139+
this % dt_eps(3) = 1.0_RP / max(sqrt(sum_dt_weight / mesh % no_of_allElements), 1e-10_RP)
140+
end if
141+
142+
if ((this % dt_eps(1) == this % dt_eps(2)) .and. (this % dt_eps(2) == 1.0_RP)) then
143+
this % dt_eps(1) = this % dt_eps(3) ! Initialize the first epsilon value
144+
this % dt_eps(2) = this % dt_eps(3) ! Initialize the second epsilon value
145+
else
146+
this % dt_eps(3) = min(this % dt_eps(3), this % dt_eps(2) * 10.0_RP) ! Ensure that the third epsilon is not too large
147+
end if
148+
149+
if (this % dt_eps(2) > this % dt_eps(3)) then
150+
if (this % pAdapted) then !Right after p-Adaptation
151+
this % b3 = 0.15_RP
152+
this % b2 = -0.15_RP
153+
this % pAdapted = .FALSE.
154+
else
155+
this % b3 = 0.2_RP
156+
this % b2 = -0.2_RP
157+
end if
158+
else
159+
this % b3 = 0.2_RP
160+
if (this % dt_eps(2) < 1.0_RP) then
161+
this % b2 = -0.1_RP
162+
else if (this % dt_eps(2) / this % dt_eps(3) < 0.5_RP) then
163+
if (this % pAdapted) then
164+
this % b3 = 0.15_RP
165+
this % b2 = max(-0.15_RP * (1.0_RP + 0.01_RP * this % dt_eps(3) / this % dt_eps(2)), -0.2_RP)
166+
this % pAdapted = .FALSE.
167+
else
168+
this % b2 = max(-0.2_RP * (1.0_RP + 0.01_RP * this % dt_eps(3) / this % dt_eps(2)), -0.25_RP)
169+
end if
170+
else
171+
this % b2 = -0.2_RP
172+
end if
173+
end if
174+
175+
dt_lim = dt_limiter(this % dt_eps(3) ** (this % b3/this % k) * &
176+
this % dt_eps(2) ** (this % b2/this % k) * &
177+
this % dt_eps(1) ** (this % b1/this % k))
178+
179+
if (dt_lim > 1.05_RP .or. dt_lim < 0.95_RP) then
180+
181+
this % dt_eps(1) = this % dt_eps(2)
182+
this % dt_eps(2) = this % dt_eps(3)
183+
184+
if (dt_lim < 0.85_RP) then
185+
dt_lim = dt_lim * 0.9_RP
186+
else if (dt_lim < 0.95_RP) then
187+
dt_lim = dt_lim ** 1.2_RP
188+
else if (dt_lim < 1.0_RP) then
189+
dt_lim = dt_lim
190+
else
191+
if (this % dt_eps(3) > 100.0_RP) then
192+
dt_lim = dt_lim ** (1.1_RP + this % dt_eps(3) / 10000.0_RP)
193+
else
194+
dt_lim = dt_lim * (1.0_RP - 0.2_RP / this % dt_eps(3))
195+
end if
196+
end if
197+
198+
dt_lim = min(max(dt_lim, min_dt_lim), max_dt_lim) ! Limit the time step
199+
dt = dt * dt_lim
200+
201+
dt = min(max(dt, this % min_dt), this % max_dt) ! Ensure dt is not below the minimum value or above the maximum value
202+
203+
this % dtAdaptationStep = dt
204+
205+
end if !end if dt_lim > 1.05_RP .or. dt_lim < 0.95_RP
206+
207+
end subroutine AdaptiveTimeStep_Update
208+
209+
function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight)
210+
implicit none
211+
class(adaptiveTimeStep_t), intent(in) :: this
212+
type(Element) , intent(in) :: e
213+
real(kind=RP) :: dt_weight
214+
!! ---------------
215+
integer :: i, j, k, dir, Ndir
216+
integer :: Pxyz(3)
217+
real(kind=RP) :: Q_error(0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3)), QLowRK_error(0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3))
218+
219+
#ifdef FLOW
220+
! Initialization of P
221+
Pxyz = e % Nxyz
222+
223+
Ndir = 7 !Number of available error variables
224+
225+
dt_weight = 0.0_RP
226+
227+
!Select the error variable
228+
do i = 0, Pxyz(3)
229+
do j = 0, Pxyz(2)
230+
do k = 0, Pxyz(1)
231+
do dir = 1, Ndir
232+
if(mod(this % error_variable, 3) == mod(dir, 3) .and. dir <= 3) then
233+
if (this % error_variable <= 6) then
234+
Q_error(k, j, i) = e % storage % Q(dir+1, k, j, i)
235+
QLowRK_error(k, j, i) = e % storage % QLowRK(dir+1, k, j, i)
236+
if (this % error_variable <= 3) then
237+
Q_error(k, j, i) = Q_error(k, j, i) / e % storage % Q(1, k, j, i)
238+
QLowRK_error(k, j, i) = QLowRK_error(k, j, i) / e % storage % QLowRK(1, k, j, i)
239+
end if
240+
else if (this % error_variable == 7) then
241+
#ifdef NAVIERSTOKES
242+
!Pressure
243+
Q_error(k, j, i) = thermodynamics % gammaMinus1*(e % storage % Q(5,k,j,i) - 0.5_RP*(e % storage % Q(2,k,j,i)**2 + e % storage % Q(3,k,j,i)**2 + e % storage % Q(4,k,j,i)**2)/e % storage % Q(1,k,j,i))
244+
QLowRK_error(k, j, i) = thermodynamics % gammaMinus1*(e % storage % QLowRK(5,k,j,i) - 0.5_RP*(e % storage % QLowRK(2,k,j,i)**2 + e % storage % QLowRK(3,k,j,i)**2 + e % storage % QLowRK(4,k,j,i)**2)/e % storage % QLowRK(1,k,j,i))
245+
#elif defined(MULTIPHASE)
246+
!Pressure
247+
Q_error(k, j, i) = e % storage % QNS(IMP,k,j,i)
248+
QLowRK_error(k, j, i) = e % storage % QLowRK(IMP,k,j,i)
249+
#endif
250+
end if
251+
end if
252+
end do
253+
dt_weight = dt_weight + (Q_error(k, j, i) - QLowRK_error(k, j, i))**2.0_RP
254+
end do
255+
end do
256+
end do
257+
258+
dt_weight = dt_weight / (e % storage % sensor)**2.0_RP
259+
dt_weight = dt_weight / ((Pxyz(1)+1) * (Pxyz(2)+1) * (Pxyz(3)+1)) !Average over all Gauss points
260+
#endif
261+
end function adaptiveTimeStep_ComputeWeights
262+
263+
subroutine adaptiveTimeStep_HasToAdapt(this, t, dt, dtHasToAdapt)
264+
implicit none
265+
class(adaptiveTimeStep_t), intent(in) :: this
266+
real(kind=RP), intent(in) :: t, dt
267+
logical, intent(out) :: dtHasToAdapt
268+
!! ---------------
269+
dtHasToAdapt = (t - this % lastAdaptationTime + dt) >= this % dtAdaptationStep
270+
271+
end subroutine adaptiveTimeStep_HasToAdapt
272+
273+
function dt_limiter(a) result(dt)
274+
implicit none
275+
real(kind=RP), intent(in) :: a
276+
real(kind=RP) :: dt
277+
!! ---------------
278+
dt = 1.0_RP + atan(a-1.0_RP)
279+
end function dt_limiter
280+
end module AdaptiveTimeStepClass

0 commit comments

Comments
 (0)