Skip to content

Commit f0599d7

Browse files
DhueperDhueper
authored
Merge pull request #269 from loganoz/david_development
Adaptive Time Step - Multi Level RK compatibility
2 parents ba26917 + 81a62d4 commit f0599d7

File tree

5 files changed

+67
-21
lines changed

5 files changed

+67
-21
lines changed

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

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -66,8 +66,9 @@ Module ElementClass
6666
INTEGER :: nodeIDs(8)
6767
integer :: faceIDs(6)
6868
integer :: faceSide(6)
69-
integer :: MLevel ! RK Level
70-
real(kind=RP) :: ML_CFL ! CFL storage for Multi Level RK
69+
integer :: MLevel ! RK Level
70+
real(kind=RP) :: ML_CFL ! CFL storage for Multi Level RK
71+
real(kind=RP) :: ML_error_ratio ! Ratio between temporal and spatial error relative to the global ratio
7172
INTEGER, DIMENSION(3) :: Nxyz ! Polynomial orders in every direction (Nx,Ny,Nz)
7273
real(kind=RP) :: hn ! Ratio of size and polynomial order
7374
TYPE(MappedGeometry) :: geom
@@ -124,7 +125,8 @@ SUBROUTINE HexElement_Construct( self, Nx, Ny, Nz, nodeIDs, eID, globID)
124125
self % boundaryName = emptyBCName
125126
self % hasSharedFaces = .false.
126127
self % NumberOfConnections = 0
127-
self % MLevel = 1
128+
self % MLevel = 1
129+
self % ML_error_ratio = 1.0_RP
128130
!
129131
! ----------------------------------------
130132
! Solution Storage is allocated separately

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

Lines changed: 26 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -101,22 +101,23 @@ end subroutine adaptiveTimeStep_Destruct
101101
! -----------------------------------------------
102102
subroutine adaptiveTimeStep_Update(this, mesh, t, dt)
103103
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
104+
class(adaptiveTimeStep_t) , intent(inout) :: this
105+
class(HexMesh) , intent(inout) :: mesh
106+
real(kind=RP) , intent(in) :: t
107+
real(kind=RP) , intent(inout) :: dt
108108
!! ---------------
109109
! Local variables
110110
!! ---------------
111111
integer :: eID, ierr, i
112-
real(kind=RP) :: dt_weight, sum_dt_weight, dt_lim, min_dt_lim, max_dt_lim
112+
real(kind=RP) :: dt_weight, sum_dt_weight, avg_sum_dt_weight, dt_lim, min_dt_lim, max_dt_lim
113113

114114
min_dt_lim = 0.5_RP ! Minimum limit for the time step
115115
max_dt_lim = 1.3_RP ! Maximum limit for the time step
116116

117117
this % lastAdaptationTime = t
118118
dt_weight = 0.0_RP
119119
sum_dt_weight = 0.0_RP
120+
avg_sum_dt_weight = 0.0_RP
120121
!$omp parallel shared(dt_weight, mesh, this)
121122
!$omp do reduction(+:dt_weight) schedule(runtime)
122123
do eID = 1, mesh % no_of_elements
@@ -133,10 +134,18 @@ subroutine adaptiveTimeStep_Update(this, mesh, t, dt)
133134
sum_dt_weight = dt_weight
134135
end if
135136

137+
avg_sum_dt_weight = sum_dt_weight / mesh % no_of_allElements
138+
139+
!$omp parallel do schedule(runtime)
140+
do eID = 1, mesh % no_of_elements
141+
mesh % elements(eID) % ML_error_ratio = mesh % elements(eID) % ML_error_ratio / avg_sum_dt_weight
142+
end do
143+
!$omp end parallel do
144+
136145
if (isnan(sum_dt_weight)) then
137146
this % dt_eps(3) = 1e-10_RP
138147
else
139-
this % dt_eps(3) = 1.0_RP / max(sqrt(sum_dt_weight / mesh % no_of_allElements), 1e-10_RP)
148+
this % dt_eps(3) = 1.0_RP / max(sqrt(avg_sum_dt_weight), 1e-10_RP)
140149
end if
141150

142151
if ((this % dt_eps(1) == this % dt_eps(2)) .and. (this % dt_eps(2) == 1.0_RP)) then
@@ -208,9 +217,9 @@ end subroutine AdaptiveTimeStep_Update
208217

209218
function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight)
210219
implicit none
211-
class(adaptiveTimeStep_t), intent(in) :: this
212-
type(Element) , intent(in) :: e
213-
real(kind=RP) :: dt_weight
220+
class(adaptiveTimeStep_t), intent(in) :: this
221+
type(Element) , intent(inout) :: e
222+
real(kind=RP) :: dt_weight
214223
!! ---------------
215224
integer :: i, j, k, dir, Ndir
216225
integer :: Pxyz(3)
@@ -220,7 +229,7 @@ function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight)
220229
! Initialization of P
221230
Pxyz = e % Nxyz
222231

223-
Ndir = 7 !Number of available error variables
232+
Ndir = 3
224233

225234
dt_weight = 0.0_RP
226235

@@ -247,6 +256,10 @@ function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight)
247256
Q_error(k, j, i) = e % storage % QNS(IMP,k,j,i)
248257
QLowRK_error(k, j, i) = e % storage % QLowRK(IMP,k,j,i)
249258
#endif
259+
else if (this % error_variable == 8) then
260+
!Density
261+
Q_error(k, j, i) = e % storage % Q(1,k,j,i)
262+
QLowRK_error(k, j, i) = e % storage % QLowRK(1,k,j,i)
250263
end if
251264
end if
252265
end do
@@ -257,6 +270,9 @@ function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight)
257270

258271
dt_weight = dt_weight / (e % storage % sensor)**2.0_RP
259272
dt_weight = dt_weight / ((Pxyz(1)+1) * (Pxyz(2)+1) * (Pxyz(3)+1)) !Average over all Gauss points
273+
! MLRK correction
274+
dt_weight = dt_weight / (3.0_RP ** (e % MLevel - 1))**2.0_RP
275+
e % ML_error_ratio = dt_weight
260276
#endif
261277
end function adaptiveTimeStep_ComputeWeights
262278

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

Lines changed: 25 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1529,10 +1529,18 @@ SUBROUTINE TakeMLRK3Step( mesh, particles, t, deltaT, ComputeTimeDerivative, dt_
15291529

15301530

15311531

1532-
INTEGER :: i, j, id, lID, locLevel, k2, k3, k1, nLevel
1532+
INTEGER :: i, j, id, lID, locLevel, k2, k3, k1, nLevel
15331533
INTEGER, ALLOCATABLE :: k(:)
15341534
REAL(KIND=RP) :: deltaStep(3), corrector, deltaTLF
1535-
REAL(KIND=RP), ALLOCATABLE :: cL(:) , tk(:), deltaTL(:) !
1535+
REAL(KIND=RP), ALLOCATABLE :: cL(:) , tk(:), deltaTL(:)
1536+
1537+
logical :: updateQLowRK
1538+
1539+
if (present(dtAdaptation)) then
1540+
updateQLowRK = dtAdaptation
1541+
else
1542+
updateQLowRK = .false.
1543+
end if
15361544

15371545
nLevel = mesh % MLRK % maxLevel
15381546
allocate(k(nLevel), cL(nLevel), tk(nLevel), deltaTL(nLevel))
@@ -1543,8 +1551,8 @@ SUBROUTINE TakeMLRK3Step( mesh, particles, t, deltaT, ComputeTimeDerivative, dt_
15431551
d(3)=1.0_RP/4.0_RP
15441552

15451553
deltaStep(1) = b(2)
1546-
deltaStep(2) = b(3)-b(2)
1547-
deltaStep(3) = 1.0_RP-b(3)
1554+
deltaStep(2) = b(3)-b(2)
1555+
deltaStep(3) = 1.0_RP-b(3)
15481556

15491557
deltaTL(:) = deltaT
15501558
tk(:) = t
@@ -1553,8 +1561,19 @@ SUBROUTINE TakeMLRK3Step( mesh, particles, t, deltaT, ComputeTimeDerivative, dt_
15531561

15541562
k(:) = 3
15551563
do k1 = 1,3
1556-
tk(:) = t + b(k1)*deltaT
1557-
call ComputeTimeDerivative( mesh, particles, tk(1), CTD_IGNORE_MODE)
1564+
tk(:) = t + b(k1)*deltaT
1565+
call ComputeTimeDerivative( mesh, particles, tk(1), CTD_IGNORE_MODE)
1566+
1567+
if (k1==1 .and. updateQLowRK) then !For adaptive time step only, update QLowRK
1568+
!$omp parallel do schedule(runtime)
1569+
do id = 1, SIZE( mesh % elements )
1570+
#ifdef FLOW
1571+
mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + deltaT*mesh % elements(id) % storage % QDot
1572+
#endif
1573+
end do ! id
1574+
!$omp end parallel do
1575+
end if
1576+
15581577
locLevel = 1
15591578
! -------------------------------------------------------------------------------------------------------------------------------
15601579
! LEVEL 2-LEVEL N-1

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

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -140,14 +140,22 @@ SUBROUTINE constructTimeIntegrator(self,controlVariables, sem, initial_iter, ini
140140
self % dt = controlVariables % doublePrecisionValueForKey("dt")
141141
else
142142
self % dt = 1e-8_RP
143-
write(STD_OUT,*) "Warning: 'dt' keyword not specified for adaptive dt. Using initial minimum value of ", self % dt
143+
if (MPI_Process % isRoot ) then
144+
write(STD_OUT,*) "Warning: 'dt' keyword not specified for adaptive dt. Using initial minimum value of ", self % dt
145+
end if
144146
end if
145147
if ( controlVariables % ContainsKey("explicit method") ) then
146148
keyword = controlVariables % StringValueForKey("explicit method",LINE_LENGTH)
147149
call toLower(keyword)
148150
select case (keyword)
149151
case(RK3_NAME)
150-
write(STD_OUT,*) "Using 'RK3' method with adaptive dt."
152+
if (MPI_Process % isRoot ) then
153+
write(STD_OUT,*) "Using 'RK3' method with adaptive dt."
154+
end if
155+
case(ML_RK3_NAME)
156+
if (MPI_Process % isRoot ) then
157+
write(STD_OUT,*) "Using 'ML-RK3' method with adaptive dt."
158+
end if
151159
case default
152160
error stop "Error, only 'RK3' method is implemented for adaptive dt"
153161
end select

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -482,6 +482,7 @@ subroutine OverEnrichRegions(overenriching,mesh,NNew,Nmax,Nmin)
482482
integer :: cornerID ! Corner counter
483483
logical :: enriched(mesh % no_of_elements)
484484
logical :: is_inside
485+
integer :: i
485486
!---------------------------------------
486487

487488
if (.not. allocated(overenriching) ) return

0 commit comments

Comments
 (0)