pAdaptationRL Gauss-Lobatto nodes test case and AdaptiveTimeStep for MixedRK scheme

.github/workflows/CI_parallel.yml

@@ -1015,6 +1015,23 @@ jobs:
           # mpiexec -n 64 ./horses3d.ns ActuatorLineController.control
         # if: '!cancelled()'
 
+#
+# 34) Cylinder with Reinforcement Learning p-adaptation and Gauss-Lobatto nodes
+# --------------------------------------------
+      - name: Build Cylinder_pAdaptationRL_GaussLobatto
+        working-directory: ./Solver/test/NavierStokes/Cylinder_pAdaptationRL_GaussLobatto/SETUP
+        run: |
+          source /opt/intel/oneapi/setvars.sh || true
+          make MODE=${{matrix.mode}} COMPILER=${{matrix.compiler}} COMM=${{matrix.comm}} ENABLE_THREADS=${{matrix.enable_threads}} WITH_METIS=YES
+        if: '!cancelled()'
+
+      - name: Run Cylinder_pAdaptationRL_GaussLobatto
+        working-directory: ./Solver/test/NavierStokes/Cylinder_pAdaptationRL_GaussLobatto
+        run: |
+          source /opt/intel/oneapi/setvars.sh || true
+          mpiexec -n 64 ./horses3d.ns Cylinder_pAdaptationRL_GaussLobatto.control
+        if: '!cancelled()'
+
 ########################################################################
 #######                       PARTICLES                         ########
 ########################################################################

.github/workflows/CI_sequential 2.yml

@@ -804,3 +804,21 @@ jobs:
           ./horses3d.ns Cylinder_AdaptiveTimeStep_pAdaptationRL.control
         if: '!cancelled()'
 
+#
+# 34) Cylinder with Reinforcement Learning p-adaptation and Gauss-Lobatto nodes
+# --------------------------------------------
+
+      - name: Build Cylinder_pAdaptationRL_GaussLobatto
+        working-directory: ./Solver/test/NavierStokes/Cylinder_pAdaptationRL_GaussLobatto/SETUP
+        run: |
+          source /opt/intel/oneapi/setvars.sh || true
+          make MODE=${{matrix.mode}} COMPILER=${{matrix.compiler}} COMM=${{matrix.comm}} ENABLE_THREADS=${{matrix.enable_threads}} WITH_MKL=${{matrix.mkl}} WITH_HDF5=${{matrix.hdf5}}
+        if: '!cancelled()'
+
+      - name: Run Cylinder_pAdaptationRL_GaussLobatto
+        working-directory: ./Solver/test/NavierStokes/Cylinder_pAdaptationRL_GaussLobatto
+        run: |
+          source /opt/intel/oneapi/setvars.sh || true
+          ./horses3d.ns Cylinder_pAdaptationRL_GaussLobatto.control
+        if: '!cancelled()'
+

Solver/configure

@@ -80,6 +80,7 @@ TEST_CASES="./Euler/BoxAroundCircle \
             ./NavierStokes/ActuatorLineInterpolation \
             ./NavierStokes/ActuatorLineProjection \
             ./NavierStokes/Cylinder_AdaptiveTimeStep_pAdaptationRL \
+            ./NavierStokes/Cylinder_pAdaptationRL_GaussLobatto \
             ./NavierStokes/ABLBoundaryMapping \
             ./NavierStokes/ActuatorLineController \
             ./IncompressibleNS/Convergence \

Solver/src/libs/timeintegrator/ExplicitMethods.f90

@@ -276,6 +276,14 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
       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]
       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]
 
+      logical :: updateQLowRK
+
+      if (present(dtAdaptation)) then
+         updateQLowRK = dtAdaptation
+      else
+         updateQLowRK = .false.
+      end if
+
       allocate(phase1_mask(size(mesh % elements)))
       allocate(phase2_mask(size(mesh % elements)))
 
@@ -287,10 +295,12 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
 !$omp parallel do schedule(runtime) reduction(+:phase1_count,phase2_count) private(id)
       do id = 1, size(mesh % elements)
          if ( all(mesh % elements(id) % storage % Q(1,:,:,:) > 1.0_RP - 1e-8_RP) ) then
+             mesh % elements(id) % MLevel = 1 !required for adaptive time step
              phase1_mask(id) = .true.
              phase2_mask(id) = .false.
              phase1_count = phase1_count + 1
          else
+            mesh % elements(id) % MLevel = 4 !required for adaptive time step
             phase1_mask(id) = .false.
             phase2_mask(id) = .true.
             phase2_count = phase2_count + 1
@@ -374,6 +384,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
                if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
             end if
 
+            if (k==1 .and. updateQLowRK) then !For adaptive time step only, update QLowRK
+!$omp parallel do schedule(runtime)
+               do id = 1, SIZE( mesh % elements )
+                  if (phase1_mask(id)) then
+#ifdef FLOW 
+                     mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + dt_vec(id)*mesh % elements(id) % storage % QDot
+#endif
+                  end if
+               end do ! id
+!$omp end parallel do
+            end if
+
 !$omp parallel do schedule(runtime)
             DO id = 1, SIZE( mesh % elements )
                if (phase1_mask(id)) then
@@ -433,6 +455,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
                if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
             end if
 
+            if (k==1 .and. updateQLowRK) then !For adaptive time step only, update QLowRK
+!$omp parallel do schedule(runtime)
+               do id = 1, SIZE( mesh % elements )
+                  if (phase2_mask(id)) then
+#ifdef FLOW 
+                     mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + dt_vec(id)*mesh % elements(id) % storage % QDot
+#endif
+                  end if
+               end do ! id
+!$omp end parallel do
+            end if
+
 !$omp parallel do schedule(runtime)
             DO id = 1, SIZE( mesh % elements )
                if (phase2_mask(id)) then
@@ -479,6 +513,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
                if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
             end if
 
+            if (k==1 .and. updateQLowRK) then !For adaptive time step only, update QLowRK
+!$omp parallel do schedule(runtime)
+               do id = 1, SIZE( mesh % elements )
+                  if (phase1_mask(id)) then
+#ifdef FLOW 
+                     mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + deltaT*mesh % elements(id) % storage % QDot
+#endif
+                  end if
+               end do ! id
+!$omp end parallel do
+            end if
+
 !$omp parallel do schedule(runtime)
             DO id = 1, SIZE( mesh % elements )
                if (phase1_mask(id)) then
@@ -538,6 +584,18 @@ SUBROUTINE TakeMixedRKStep( mesh, particles, t, deltaT, ComputeTimeDerivative ,
                if (dts) call ComputePseudoTimeDerivative(mesh, tk, global_dt)
             end if
 
+            if (k==1 .and. updateQLowRK) then !For adaptive time step only, update QLowRK
+!$omp parallel do schedule(runtime)
+               do id = 1, SIZE( mesh % elements )
+                  if (phase2_mask(id)) then
+#ifdef FLOW 
+                     mesh % elements(id) % storage % QLowRK = mesh % elements(id) % storage % Q + deltaT*mesh % elements(id) % storage % QDot
+#endif
+                  end if
+               end do ! id
+!$omp end parallel do
+            end if
+
 !$omp parallel do schedule(runtime)
             DO id = 1, SIZE( mesh % elements )
                if (phase2_mask(id)) then

Solver/src/libs/timeintegrator/TimeIntegrator.f90

@@ -156,8 +156,12 @@ SUBROUTINE constructTimeIntegrator(self,controlVariables, sem, initial_iter, ini
                         if (MPI_Process % isRoot ) then
                            write(STD_OUT,*) "Using 'ML-RK3' method with adaptive dt."
                         end if
+                     case(MIXED_RK_NAME)
+                        if (MPI_Process % isRoot ) then
+                           write(STD_OUT,*) "Using 'MIXED-RK' method with adaptive dt."
+                        end if
                      case default
-                        error stop "Error, only 'RK3' method is implemented for adaptive dt"
+                        error stop "Error, only 'RK3'-based methods are implemented for adaptive dt"
                      end select
                   else
                      write(STD_OUT,*) "Using 'RK3' method with adaptive dt by default."

Solver/test/NavierStokes/Cylinder_pAdaptationRL_GaussLobatto/Cylinder_pAdaptationRL_GaussLobatto.control

@@ -0,0 +1,74 @@
+Flow equations        = "NS"
+mesh file name        = "MESH/cyl_circ.msh"
+Polynomial order      = 2
+Number of time steps  = 75
+Output Interval       = 10
+Convergence tolerance = 1.d-10
+cfl                   = 0.1
+dcfl                  = 0.1
+mach number           = 0.2
+Reynolds number       = 40.0
+AOA theta             = 0.0
+AOA phi               = 0.0
+solution file name    = "RESULTS/Cylinder.hsol"
+save gradients with solution = .false.
+restart               = .false.
+restart file name     = "RESULTS/Cylinder.hsol"
+riemann solver        = roe
+simulation type       = steady-state
+time integration      = explicit
+explicit method       = RK3
+Discretization nodes  = Gauss-Lobatto
+
+#define boundary cylinder
+  type = NoSlipWall
+#end
+
+#define boundary left__right
+  type = FreeSlipWall
+#end
+
+#define boundary inlet
+  type = Inflow
+#end
+
+#define boundary outlet
+  type = Outflow
+#end
+
+!-Adaptation-----------------------------------
+#define p-adaptation
+   adaptation type       = RL
+   agent file            = policy_padaptation_gauss_lobatto/p_adaptation_policy
+   tolerance             = 1d-2
+   Nmax                  = [6, 6, 2]
+   Nmin                  = [2, 2, 2]
+   adjust nz             = .FALSE.
+   restart files         = .FALSE.
+   order across faces    = N*2/3
+   mode                  = iteration
+   interval              = 50
+#end
+!---------------------------------------------
+
+#define surface monitor 1
+   Name = cyl-drag
+   Marker = cylinder
+   Variable = drag
+   Direction = [1.0,0.0,0.0]
+   Reference surface = 1.0	
+#end
+
+#define surface monitor 2
+   Name = cyl-lift
+   Marker = cylinder
+   Variable = lift
+   Direction = [0.0,1.0,0.0]
+   Reference surface = 1.0
+#end
+
+#define probe 1
+   Name = wake_u
+   Position = [2.0,0.0,0.5]
+   Variable = u
+#end