@@ -276,6 +276,14 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
276276 real (kind= RP), parameter :: b_RK14_4(N_STAGES_RK14_4) = [0.0000000000000000_RP , 0.0367762454319673_RP , 0.1249685262725025_RP , 0.2446177702277698_RP , 0.2476149531070420_RP , 0.2969311120382472_RP , 0.3978149645802642_RP , 0.5270854589440328_RP , 0.6981269994175695_RP , 0.8190890835352128_RP , 0.8527059887098624_RP , 0.8604711817462826_RP , 0.8627060376969976_RP , 0.8734213127600976_RP ]
277277 real (kind= RP), parameter :: c_RK14_4(N_STAGES_RK14_4) = [0.0367762454319673_RP , 0.3136296607553959_RP , 0.1531848691869027_RP , 0.0030097086818182_RP , 0.3326293790646110_RP , 0.2440251405350864_RP , 0.3718879239592277_RP , 0.6204126221582444_RP , 0.1524043173028741_RP , 0.0760894927419266_RP , 0.0077604214040978_RP , 0.0024647284755382_RP , 0.0780348340049386_RP , 5.5059777270269628_RP ]
278278
279+ logical :: updateQLowRK
280+
281+ if (present (dtAdaptation)) then
282+ updateQLowRK = dtAdaptation
283+ else
284+ updateQLowRK = .false.
285+ end if
286+
279287 allocate (phase1_mask(size (mesh % elements)))
280288 allocate (phase2_mask(size (mesh % elements)))
281289
@@ -287,10 +295,12 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
287295! $omp parallel do schedule(runtime) reduction(+:phase1_count,phase2_count) private(id)
288296 do id = 1 , size (mesh % elements)
289297 if ( all (mesh % elements(id) % storage % Q(1 ,:,:,:) > 1.0_RP - 1e-8_RP ) ) then
298+ mesh % elements(id) % MLevel = 1 ! required for adaptive time step
290299 phase1_mask(id) = .true.
291300 phase2_mask(id) = .false.
292301 phase1_count = phase1_count + 1
293302 else
303+ mesh % elements(id) % MLevel = 4 ! required for adaptive time step
294304 phase1_mask(id) = .false.
295305 phase2_mask(id) = .true.
296306 phase2_count = phase2_count + 1
@@ -374,6 +384,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
374384 if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
375385 end if
376386
387+ if (k== 1 .and. updateQLowRK) then ! For adaptive time step only, update QLowRK
388+ ! $omp parallel do schedule(runtime)
389+ do id = 1 , SIZE ( mesh % elements )
390+ if (phase1_mask(id)) then
391+ #ifdef FLOW
392+ mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + dt_vec(id)* mesh % elements(id) % storage % QDot
393+ #endif
394+ end if
395+ end do ! id
396+ ! $omp end parallel do
397+ end if
398+
377399! $omp parallel do schedule(runtime)
378400 DO id = 1 , SIZE ( mesh % elements )
379401 if (phase1_mask(id)) then
@@ -433,6 +455,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
433455 if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
434456 end if
435457
458+ if (k== 1 .and. updateQLowRK) then ! For adaptive time step only, update QLowRK
459+ ! $omp parallel do schedule(runtime)
460+ do id = 1 , SIZE ( mesh % elements )
461+ if (phase2_mask(id)) then
462+ #ifdef FLOW
463+ mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + dt_vec(id)* mesh % elements(id) % storage % QDot
464+ #endif
465+ end if
466+ end do ! id
467+ ! $omp end parallel do
468+ end if
469+
436470! $omp parallel do schedule(runtime)
437471 DO id = 1 , SIZE ( mesh % elements )
438472 if (phase2_mask(id)) then
@@ -479,6 +513,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
479513 if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
480514 end if
481515
516+ if (k== 1 .and. updateQLowRK) then ! For adaptive time step only, update QLowRK
517+ ! $omp parallel do schedule(runtime)
518+ do id = 1 , SIZE ( mesh % elements )
519+ if (phase1_mask(id)) then
520+ #ifdef FLOW
521+ mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + deltaT* mesh % elements(id) % storage % QDot
522+ #endif
523+ end if
524+ end do ! id
525+ ! $omp end parallel do
526+ end if
527+
482528! $omp parallel do schedule(runtime)
483529 DO id = 1 , SIZE ( mesh % elements )
484530 if (phase1_mask(id)) then
@@ -538,6 +584,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
538584 if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
539585 end if
540586
587+ if (k== 1 .and. updateQLowRK) then ! For adaptive time step only, update QLowRK
588+ ! $omp parallel do schedule(runtime)
589+ do id = 1 , SIZE ( mesh % elements )
590+ if (phase2_mask(id)) then
591+ #ifdef FLOW
592+ mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + deltaT* mesh % elements(id) % storage % QDot
593+ #endif
594+ end if
595+ end do ! id
596+ ! $omp end parallel do
597+ end if
598+
541599! $omp parallel do schedule(runtime)
542600 DO id = 1 , SIZE ( mesh % elements )
543601 if (phase2_mask(id)) then
0 commit comments