diff --git a/Solver/src/libs/mesh/HexElementClass.f90 b/Solver/src/libs/mesh/HexElementClass.f90 index 1c856a0d39..6abe67b41b 100644 --- a/Solver/src/libs/mesh/HexElementClass.f90 +++ b/Solver/src/libs/mesh/HexElementClass.f90 @@ -66,8 +66,9 @@ Module ElementClass INTEGER :: nodeIDs(8) integer :: faceIDs(6) integer :: faceSide(6) - integer :: MLevel ! RK Level - real(kind=RP) :: ML_CFL ! CFL storage for Multi Level RK + integer :: MLevel ! RK Level + real(kind=RP) :: ML_CFL ! CFL storage for Multi Level RK + real(kind=RP) :: ML_error_ratio ! Ratio between temporal and spatial error relative to the global ratio INTEGER, DIMENSION(3) :: Nxyz ! Polynomial orders in every direction (Nx,Ny,Nz) real(kind=RP) :: hn ! Ratio of size and polynomial order TYPE(MappedGeometry) :: geom @@ -124,7 +125,8 @@ SUBROUTINE HexElement_Construct( self, Nx, Ny, Nz, nodeIDs, eID, globID) self % boundaryName = emptyBCName self % hasSharedFaces = .false. self % NumberOfConnections = 0 - self % MLevel = 1 + self % MLevel = 1 + self % ML_error_ratio = 1.0_RP ! ! ---------------------------------------- ! Solution Storage is allocated separately diff --git a/Solver/src/libs/timeintegrator/AdaptiveTimeStepClass.f90 b/Solver/src/libs/timeintegrator/AdaptiveTimeStepClass.f90 index c6d6ddff25..4ad319a512 100644 --- a/Solver/src/libs/timeintegrator/AdaptiveTimeStepClass.f90 +++ b/Solver/src/libs/timeintegrator/AdaptiveTimeStepClass.f90 @@ -101,15 +101,15 @@ end subroutine adaptiveTimeStep_Destruct ! ----------------------------------------------- subroutine adaptiveTimeStep_Update(this, mesh, t, dt) implicit none - class(adaptiveTimeStep_t) , intent(inout) :: this - class(HexMesh) , intent(in) :: mesh - real(kind=RP) , intent(in) :: t - real(kind=RP) , intent(inout) :: dt + class(adaptiveTimeStep_t) , intent(inout) :: this + class(HexMesh) , intent(inout) :: mesh + real(kind=RP) , intent(in) :: t + real(kind=RP) , intent(inout) :: dt !! --------------- ! Local variables !! --------------- integer :: eID, ierr, i - real(kind=RP) :: dt_weight, sum_dt_weight, dt_lim, min_dt_lim, max_dt_lim + real(kind=RP) :: dt_weight, sum_dt_weight, avg_sum_dt_weight, dt_lim, min_dt_lim, max_dt_lim min_dt_lim = 0.5_RP ! Minimum limit for the time step max_dt_lim = 1.3_RP ! Maximum limit for the time step @@ -117,6 +117,7 @@ subroutine adaptiveTimeStep_Update(this, mesh, t, dt) this % lastAdaptationTime = t dt_weight = 0.0_RP sum_dt_weight = 0.0_RP + avg_sum_dt_weight = 0.0_RP !$omp parallel shared(dt_weight, mesh, this) !$omp do reduction(+:dt_weight) schedule(runtime) do eID = 1, mesh % no_of_elements @@ -133,10 +134,18 @@ subroutine adaptiveTimeStep_Update(this, mesh, t, dt) sum_dt_weight = dt_weight end if + avg_sum_dt_weight = sum_dt_weight / mesh % no_of_allElements + +!$omp parallel do schedule(runtime) + do eID = 1, mesh % no_of_elements + mesh % elements(eID) % ML_error_ratio = mesh % elements(eID) % ML_error_ratio / avg_sum_dt_weight + end do +!$omp end parallel do + if (isnan(sum_dt_weight)) then this % dt_eps(3) = 1e-10_RP else - this % dt_eps(3) = 1.0_RP / max(sqrt(sum_dt_weight / mesh % no_of_allElements), 1e-10_RP) + this % dt_eps(3) = 1.0_RP / max(sqrt(avg_sum_dt_weight), 1e-10_RP) end if 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 function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight) implicit none - class(adaptiveTimeStep_t), intent(in) :: this - type(Element) , intent(in) :: e - real(kind=RP) :: dt_weight + class(adaptiveTimeStep_t), intent(in) :: this + type(Element) , intent(inout) :: e + real(kind=RP) :: dt_weight !! --------------- integer :: i, j, k, dir, Ndir integer :: Pxyz(3) @@ -220,7 +229,7 @@ function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight) ! Initialization of P Pxyz = e % Nxyz - Ndir = 7 !Number of available error variables + Ndir = 3 dt_weight = 0.0_RP @@ -247,6 +256,10 @@ function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight) Q_error(k, j, i) = e % storage % QNS(IMP,k,j,i) QLowRK_error(k, j, i) = e % storage % QLowRK(IMP,k,j,i) #endif + else if (this % error_variable == 8) then + !Density + Q_error(k, j, i) = e % storage % Q(1,k,j,i) + QLowRK_error(k, j, i) = e % storage % QLowRK(1,k,j,i) end if end if end do @@ -257,6 +270,9 @@ function adaptiveTimeStep_ComputeWeights(this, e) result(dt_weight) dt_weight = dt_weight / (e % storage % sensor)**2.0_RP dt_weight = dt_weight / ((Pxyz(1)+1) * (Pxyz(2)+1) * (Pxyz(3)+1)) !Average over all Gauss points + ! MLRK correction + dt_weight = dt_weight / (3.0_RP ** (e % MLevel - 1))**2.0_RP + e % ML_error_ratio = dt_weight #endif end function adaptiveTimeStep_ComputeWeights diff --git a/Solver/src/libs/timeintegrator/ExplicitMethods.f90 b/Solver/src/libs/timeintegrator/ExplicitMethods.f90 index 7f236b0723..6a17972bae 100644 --- a/Solver/src/libs/timeintegrator/ExplicitMethods.f90 +++ b/Solver/src/libs/timeintegrator/ExplicitMethods.f90 @@ -1529,10 +1529,18 @@ SUBROUTINE TakeMLRK3Step( mesh, particles, t, deltaT, ComputeTimeDerivative, dt_ - INTEGER :: i, j, id, lID, locLevel, k2, k3, k1, nLevel + INTEGER :: i, j, id, lID, locLevel, k2, k3, k1, nLevel INTEGER, ALLOCATABLE :: k(:) REAL(KIND=RP) :: deltaStep(3), corrector, deltaTLF - REAL(KIND=RP), ALLOCATABLE :: cL(:) , tk(:), deltaTL(:) ! + REAL(KIND=RP), ALLOCATABLE :: cL(:) , tk(:), deltaTL(:) + + logical :: updateQLowRK + + if (present(dtAdaptation)) then + updateQLowRK = dtAdaptation + else + updateQLowRK = .false. + end if nLevel = mesh % MLRK % maxLevel allocate(k(nLevel), cL(nLevel), tk(nLevel), deltaTL(nLevel)) @@ -1543,8 +1551,8 @@ SUBROUTINE TakeMLRK3Step( mesh, particles, t, deltaT, ComputeTimeDerivative, dt_ d(3)=1.0_RP/4.0_RP deltaStep(1) = b(2) - deltaStep(2) = b(3)-b(2) - deltaStep(3) = 1.0_RP-b(3) + deltaStep(2) = b(3)-b(2) + deltaStep(3) = 1.0_RP-b(3) deltaTL(:) = deltaT tk(:) = t @@ -1553,8 +1561,19 @@ SUBROUTINE TakeMLRK3Step( mesh, particles, t, deltaT, ComputeTimeDerivative, dt_ k(:) = 3 do k1 = 1,3 - tk(:) = t + b(k1)*deltaT - call ComputeTimeDerivative( mesh, particles, tk(1), CTD_IGNORE_MODE) + tk(:) = t + b(k1)*deltaT + call ComputeTimeDerivative( mesh, particles, tk(1), CTD_IGNORE_MODE) + + if (k1==1 .and. updateQLowRK) then !For adaptive time step only, update QLowRK +!$omp parallel do schedule(runtime) + do id = 1, SIZE( mesh % elements ) +#ifdef FLOW + mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + deltaT*mesh % elements(id) % storage % QDot +#endif + end do ! id +!$omp end parallel do + end if + locLevel = 1 ! ------------------------------------------------------------------------------------------------------------------------------- ! LEVEL 2-LEVEL N-1 diff --git a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 index 61fe680193..deddabe258 100644 --- a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 +++ b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 @@ -140,14 +140,22 @@ SUBROUTINE constructTimeIntegrator(self,controlVariables, sem, initial_iter, ini self % dt = controlVariables % doublePrecisionValueForKey("dt") else self % dt = 1e-8_RP - write(STD_OUT,*) "Warning: 'dt' keyword not specified for adaptive dt. Using initial minimum value of ", self % dt + if (MPI_Process % isRoot ) then + write(STD_OUT,*) "Warning: 'dt' keyword not specified for adaptive dt. Using initial minimum value of ", self % dt + end if end if if ( controlVariables % ContainsKey("explicit method") ) then keyword = controlVariables % StringValueForKey("explicit method",LINE_LENGTH) call toLower(keyword) select case (keyword) case(RK3_NAME) - write(STD_OUT,*) "Using 'RK3' method with adaptive dt." + if (MPI_Process % isRoot ) then + write(STD_OUT,*) "Using 'RK3' method with adaptive dt." + end if + case(ML_RK3_NAME) + if (MPI_Process % isRoot ) then + write(STD_OUT,*) "Using 'ML-RK3' method with adaptive dt." + end if case default error stop "Error, only 'RK3' method is implemented for adaptive dt" end select diff --git a/Solver/src/libs/timeintegrator/pAdaptationClass.f90 b/Solver/src/libs/timeintegrator/pAdaptationClass.f90 index 23fe0318cf..23457cccf9 100644 --- a/Solver/src/libs/timeintegrator/pAdaptationClass.f90 +++ b/Solver/src/libs/timeintegrator/pAdaptationClass.f90 @@ -482,6 +482,7 @@ subroutine OverEnrichRegions(overenriching,mesh,NNew,Nmax,Nmin) integer :: cornerID ! Corner counter logical :: enriched(mesh % no_of_elements) logical :: is_inside + integer :: i !--------------------------------------- if (.not. allocated(overenriching) ) return