Skip to content

Commit 593cae7

Browse files
DhueperDhueper
committed
Adaptive Time Step - Multi Level RK compatibility
1 parent ba26917 commit 593cae7

5 files changed

Lines changed: 168 additions & 66 deletions

File tree

‎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 & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -140,14 +140,21 @@ 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
144-
end if
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
145146
if ( controlVariables % ContainsKey("explicit method") ) then
146147
keyword = controlVariables % StringValueForKey("explicit method",LINE_LENGTH)
147148
call toLower(keyword)
148149
select case (keyword)
149150
case(RK3_NAME)
150-
write(STD_OUT,*) "Using 'RK3' method with adaptive dt."
151+
if (MPI_Process % isRoot ) then
152+
write(STD_OUT,*) "Using 'RK3' method with adaptive dt."
153+
end if
154+
case(ML_RK3_NAME)
155+
if (MPI_Process % isRoot ) then
156+
write(STD_OUT,*) "Using 'ML-RK3' method with adaptive dt."
157+
end if
151158
case default
152159
error stop "Error, only 'RK3' method is implemented for adaptive dt"
153160
end select

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

Lines changed: 102 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -482,64 +482,122 @@ 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
488489

489490
enriched = .FALSE.
490491

491-
do oID = 1, size(overenriching)
492-
associate (box => overenriching(oID) )
493-
494-
element_loop: do eID=1, mesh % no_of_elements
492+
if (.not. allocated(mesh % elements_aerodynamics)) then
493+
do oID = 1, size(overenriching)
494+
associate (box => overenriching(oID) )
495495

496-
if ( enriched(eID) ) cycle element_loop
496+
element_loop: do eID=1, mesh % no_of_elements
497+
498+
if ( enriched(eID) ) cycle element_loop
499+
500+
associate ( corners => mesh % elements(eID) % hexMap % corners )
501+
502+
!
503+
! Enrich element if any of the corners is inside the region
504+
! ---------------------------------------------------------
505+
corner_loop: do cornerID=1, 8
506+
507+
is_inside = ( (corners(1,cornerID) > box % x_span(1) .and. corners(1,cornerID) < box % x_span(2)) .and. &
508+
(corners(2,cornerID) > box % y_span(1) .and. corners(2,cornerID) < box % y_span(2)) .and. &
509+
(corners(3,cornerID) > box % z_span(1) .and. corners(3,cornerID) < box % z_span(2)) )
510+
511+
if( (box % region == 1 .and. is_inside) .or. (box % region == 2 .and. (.not. is_inside))) then
512+
513+
if (box % mode == 1) then !Increase mode
514+
NNew(1,eID) = max(min(NNew(1,eID) + box % order, Nmax(1)), Nmin(1))
515+
NNew(2,eID) = max(min(NNew(2,eID) + box % order, Nmax(2)), Nmin(2))
516+
NNew(3,eID) = max(min(NNew(3,eID) + box % order, Nmax(3)), Nmin(3))
517+
else if (box % mode == 2) then !Set mode
518+
NNew(1,eID) = max(min(box % polynomial(1), Nmax(1)), Nmin(1))
519+
NNew(2,eID) = max(min(box % polynomial(2), Nmax(2)), Nmin(2))
520+
NNew(3,eID) = max(min(box % polynomial(3), Nmax(3)), Nmin(3))
521+
else if (box % mode == 3) then !Max mode
522+
NNew(1,eID) = min(max(NNew(1,eID), box % polynomial(1)), Nmax(1))
523+
NNew(2,eID) = min(max(NNew(2,eID), box % polynomial(2)), Nmax(2))
524+
NNew(3,eID) = min(max(NNew(3,eID), box % polynomial(3)), Nmax(3))
525+
else if (box % mode == 4) then !Min mode
526+
NNew(1,eID) = max(min(NNew(1,eID), box % polynomial(1)), Nmin(1))
527+
NNew(2,eID) = max(min(NNew(2,eID), box % polynomial(2)), Nmin(2))
528+
NNew(3,eID) = max(min(NNew(3,eID), box % polynomial(3)), Nmin(3))
529+
else if (box % mode == 5) then !Freeze mode
530+
NNew(1,eID) = mesh % elements(eID) % Nxyz(1)
531+
NNew(2,eID) = mesh % elements(eID) % Nxyz(2)
532+
NNew(3,eID) = mesh % elements(eID) % Nxyz(3)
533+
end if
534+
535+
enriched(eID) = .TRUE.
536+
exit corner_loop
537+
end if
538+
end do corner_loop
539+
540+
end associate
541+
end do element_loop
497542

498-
associate ( corners => mesh % elements(eID) % hexMap % corners )
543+
end associate
544+
end do
545+
else
546+
do oID = 1, size(overenriching)
547+
associate (box => overenriching(oID) )
499548

500-
!
501-
! Enrich element if any of the corners is inside the region
502-
! ---------------------------------------------------------
503-
corner_loop: do cornerID=1, 8
549+
element_loop: do i = 1, size(mesh % elements_aerodynamics)
550+
551+
eID = mesh % elements_aerodynamics(i)
552+
553+
if ( enriched(eID) ) cycle element_loop
554+
555+
associate ( corners => mesh % elements(eID) % hexMap % corners )
556+
557+
!
558+
! Enrich element if any of the corners is inside the region
559+
! ---------------------------------------------------------
560+
corner_loop: do cornerID=1, 8
504561

505-
is_inside = ( (corners(1,cornerID) > box % x_span(1) .and. corners(1,cornerID) < box % x_span(2)) .and. &
506-
(corners(2,cornerID) > box % y_span(1) .and. corners(2,cornerID) < box % y_span(2)) .and. &
507-
(corners(3,cornerID) > box % z_span(1) .and. corners(3,cornerID) < box % z_span(2)) )
562+
is_inside = ( (corners(1,cornerID) > box % x_span(1) .and. corners(1,cornerID) < box % x_span(2)) .and. &
563+
(corners(2,cornerID) > box % y_span(1) .and. corners(2,cornerID) < box % y_span(2)) .and. &
564+
(corners(3,cornerID) > box % z_span(1) .and. corners(3,cornerID) < box % z_span(2)) )
508565

509-
if( (box % region == 1 .and. is_inside) .or. (box % region == 2 .and. (.not. is_inside))) then
566+
if( (box % region == 1 .and. is_inside) .or. (box % region == 2 .and. (.not. is_inside))) then
510567

511-
if (box % mode == 1) then !Increase mode
512-
NNew(1,eID) = max(min(NNew(1,eID) + box % order, Nmax(1)), Nmin(1))
513-
NNew(2,eID) = max(min(NNew(2,eID) + box % order, Nmax(2)), Nmin(2))
514-
NNew(3,eID) = max(min(NNew(3,eID) + box % order, Nmax(3)), Nmin(3))
515-
else if (box % mode == 2) then !Set mode
516-
NNew(1,eID) = max(min(box % polynomial(1), Nmax(1)), Nmin(1))
517-
NNew(2,eID) = max(min(box % polynomial(2), Nmax(2)), Nmin(2))
518-
NNew(3,eID) = max(min(box % polynomial(3), Nmax(3)), Nmin(3))
519-
else if (box % mode == 3) then !Max mode
520-
NNew(1,eID) = min(max(NNew(1,eID), box % polynomial(1)), Nmax(1))
521-
NNew(2,eID) = min(max(NNew(2,eID), box % polynomial(2)), Nmax(2))
522-
NNew(3,eID) = min(max(NNew(3,eID), box % polynomial(3)), Nmax(3))
523-
else if (box % mode == 4) then !Min mode
524-
NNew(1,eID) = max(min(NNew(1,eID), box % polynomial(1)), Nmin(1))
525-
NNew(2,eID) = max(min(NNew(2,eID), box % polynomial(2)), Nmin(2))
526-
NNew(3,eID) = max(min(NNew(3,eID), box % polynomial(3)), Nmin(3))
527-
else if (box % mode == 5) then !Freeze mode
528-
NNew(1,eID) = mesh % elements(eID) % Nxyz(1)
529-
NNew(2,eID) = mesh % elements(eID) % Nxyz(2)
530-
NNew(3,eID) = mesh % elements(eID) % Nxyz(3)
531-
end if
532-
533-
enriched(eID) = .TRUE.
534-
exit corner_loop
535-
end if
536-
end do corner_loop
568+
if (box % mode == 1) then !Increase mode
569+
NNew(1,eID) = max(min(NNew(1,eID) + box % order, Nmax(1)), Nmin(1))
570+
NNew(2,eID) = max(min(NNew(2,eID) + box % order, Nmax(2)), Nmin(2))
571+
NNew(3,eID) = max(min(NNew(3,eID) + box % order, Nmax(3)), Nmin(3))
572+
else if (box % mode == 2) then !Set mode
573+
NNew(1,eID) = max(min(box % polynomial(1), Nmax(1)), Nmin(1))
574+
NNew(2,eID) = max(min(box % polynomial(2), Nmax(2)), Nmin(2))
575+
NNew(3,eID) = max(min(box % polynomial(3), Nmax(3)), Nmin(3))
576+
else if (box % mode == 3) then !Max mode
577+
NNew(1,eID) = min(max(NNew(1,eID), box % polynomial(1)), Nmax(1))
578+
NNew(2,eID) = min(max(NNew(2,eID), box % polynomial(2)), Nmax(2))
579+
NNew(3,eID) = min(max(NNew(3,eID), box % polynomial(3)), Nmax(3))
580+
else if (box % mode == 4) then !Min mode
581+
NNew(1,eID) = max(min(NNew(1,eID), box % polynomial(1)), Nmin(1))
582+
NNew(2,eID) = max(min(NNew(2,eID), box % polynomial(2)), Nmin(2))
583+
NNew(3,eID) = max(min(NNew(3,eID), box % polynomial(3)), Nmin(3))
584+
else if (box % mode == 5) then !Freeze mode
585+
NNew(1,eID) = mesh % elements(eID) % Nxyz(1)
586+
NNew(2,eID) = mesh % elements(eID) % Nxyz(2)
587+
NNew(3,eID) = mesh % elements(eID) % Nxyz(3)
588+
end if
589+
590+
enriched(eID) = .TRUE.
591+
exit corner_loop
592+
end if
593+
end do corner_loop
594+
595+
end associate
596+
end do element_loop
537597

538598
end associate
539-
end do element_loop
540-
541-
end associate
542-
end do
599+
end do
600+
end if
543601

544602
end subroutine OverEnrichRegions
545603
!

0 commit comments

Comments
 (0)