Add CFL Based Adaptive Time-Stepping (#515)

Co-authored-by: wilfonba <[email protected]>
Co-authored-by: Benjamin Alexander Wilfong <[email protected]>

Author: wilfonba

docs/documentation/case.md

@@ -108,7 +108,7 @@ When used, `--omni` will output profiling information for all subroutines, inclu
 Adding this argument will moderately slow down the simulation and run the MFC executable several times.
 For this reason, it should only be used with case files with few timesteps.
 
-
+<a name="restarting-cases"></a>
 ### Restarting Cases
 
 When running a simulation, MFC generates a `./restart_data` folder in the case directory that contains `lustre_*.dat` files that can be used to restart a simulation from saved timesteps.
@@ -117,11 +117,15 @@ The user can also choose to add new patches at the intermediate timestep.
 
 If you want to restart a simulation, 
 
-- Set up the initial simulation, with:
+- For a simulation that uses a constant time step set up the initial case file with: 
     - `t_step_start` : $t_i$
     - `t_step_stop`  : $t_f$
-    - `t_step_save`  : $SF$
-in which $t_i$ is the starting time, $t_f$ is the final time, and $SF$ is the saving frequency time.
+    - `t_step_save`  : $SF$ in which $t_i$ is the starting time, $t_f$ is the final time, and $SF$ is the saving frequency time.
+    For a simulation that uses adaptive time-stepping, set up the initial case file with:
+    - `n_start` : $t_i$
+    - `t_stop`  : $t_f$
+    - `t_save`  : $SF$ in which $t_i$ is the starting time, $t_f$ is the final time, and $SF$ is the saving frequency time.
+
 - Run `pre_process` and `simulation` on the case.
     - `./mfc.sh run case.py -t pre_process simulation `
 - As the simulation runs, Lustre files will be created for each saved timestep in `./restart_data`.
@@ -137,10 +141,16 @@ in which $t_i$ is the starting time, $t_f$ is the final time, and $SF$ is the sa
 			- `a_(xyz)`
 			- `(xyz)_a`
 			- `(xyz)_b`
-	- Alter the following:
+	- When using a constant time-step, alter the following:
 		- `t_step_start` : $t_s$ (the point at which the simulation will restart)
 		- `t_step_stop`  : $t_{f2}$ (new final simulation time, which can be the same as $t_f$)
 		- `t_step_save`  : ${SF}_2$ (if interested in changing the saving frequency)
+
+        If using a CFL-based time-step, alter the following:
+		- `n_start` : $t_s$ (the save file at which the simulation will restart)
+		- `t_stop`  : $t_{f2}$ (new final simulation time, which can be the same as $t_f$)
+		- `t_save`  : ${SF}_2$ (if interested in changing the saving frequency)
+
 	- Add the following:
 		- `old_ic` : 'T' (to specify that we have initial conditions from previous simulations)
 		- `old_grid` : 'T' (to specify that we have a grid from previous simulations)

docs/documentation/running.md

@@ -0,0 +1,5 @@
+# 2D IBM CFL dt (2D)
+
+## Result
+
+![Result](result.png)

examples/2D_ibm_cfl_dt/README.md

@@ -0,0 +1,107 @@
+import json
+import math
+
+Mu = 1.84E-05
+gam_a = 1.4
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    # For these computations, the cylinder is placed at the (0,0,0)
+    # domain origin. 
+    # axial direction
+    'x_domain%beg'                 : 0.0E+00,
+    'x_domain%end'                 : 8.0E-03,
+    # r direction
+    'y_domain%beg'                 : 0.0E+00,
+    'y_domain%end'                 : 6.0E-03,
+    'cyl_coord'                    : 'F',
+    'm'                            : 199,
+    'n'                            : 139,
+    'p'                            : 0,
+    'cfl_adap_dt'                  : 'T',
+    'cfl_target'                   : 0.3,
+    'n_start'                      : 0,
+    't_save'                       : 8e-5,
+    't_stop'                       : 8e-3,
+    # ==========================================================================
+    
+    # Simulation Algorithm Parameters ==========================================
+    # Only one patches are necessary, the air tube
+    'num_patches'                  : 1,
+    # Use the 5 equation model
+    'model_eqns'                   : 2,
+    'alt_soundspeed'               : 'F',
+    # One fluids: air
+    'num_fluids'                   : 1,
+    # No need to ensure the volume fractions sum to unity at the end of each
+    # time step
+    'mpp_lim'                      : 'F',
+    # Correct errors when computing speed of sound
+    'mixture_err'                  : 'T',
+    # Use TVD RK3 for time marching
+    'time_stepper'                 : 3,
+    # Use WENO5
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1.E-16,
+    # 'weno_Re_flux'                 : 'T',
+    'weno_avg'                     : 'T',
+    'avg_state'                    : 2,
+    'mapped_weno'                  : 'T',
+    'null_weights'                 : 'F',
+    'mp_weno'                      : 'T',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    # We use ghost-cell 
+    'bc_x%beg'                     : -3,
+    'bc_x%end'                     : -3,
+    'bc_y%beg'                     : -2,
+    'bc_y%end'                     : -2,
+    # Set IB to True and add 1 patch
+    'ib'                           : 'T',
+    'num_ibs'                      : 1,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'E_wrt'                        :'T',
+    'parallel_io'                  :'T',
+    'c_wrt'                        :'T',
+    # ==========================================================================
+
+    # Patch: Constant Tube filled with air =====================================
+    # Specify the cylindrical air tube grid geometry
+    'patch_icpp(1)%geometry'       : 3,
+    'patch_icpp(1)%x_centroid'     : 4.0E-03,
+    # Uniform medium density, centroid is at the center of the domain
+    'patch_icpp(1)%y_centroid'     : 3.0E-03,
+    'patch_icpp(1)%length_x'       : 8.0E-03,
+    'patch_icpp(1)%length_y'       : 6.0E-03,  
+    # Specify the patch primitive variables 
+    'patch_icpp(1)%vel(1)'         : 5E+00,
+    'patch_icpp(1)%vel(2)'         : 0.0E+00,
+    'patch_icpp(1)%pres'           : 1.E+00,
+    'patch_icpp(1)%alpha_rho(1)'   : 1.E+00,
+    'patch_icpp(1)%alpha(1)'       : 1.E+00,
+    # # ========================================================================
+
+    # Patch: Cylinder Immersed Boundary ========================================
+    'patch_ib(1)%geometry'       : 2,
+    'patch_ib(1)%x_centroid'     : 1.5E-03,
+    'patch_ib(1)%y_centroid'     : 3E-03,
+    'patch_ib(1)%radius'         : 0.4E-03,
+    'patch_ib(1)%slip'           : 'F',
+    # # ========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.E+00/(gam_a-1.E+00),  # 2.50(Not 1.40)
+    'fluid_pp(1)%pi_inf'           : 0,
+    'fluid_pp(1)%Re(1)'            : 250000,
+    # ==========================================================================
+}))

examples/2D_ibm_cfl_dt/case.py

@@ -60,9 +60,17 @@ contains
     !> Checks constraints on the time-stepping parameters.
         !! Called by s_check_inputs_common for simulation and post-processing
     subroutine s_check_inputs_time_stepping
-        @:PROHIBIT(t_step_start < 0)
-        @:PROHIBIT(t_step_stop <= t_step_start)
-        @:PROHIBIT(t_step_save > t_step_stop - t_step_start)
+        if (cfl_dt) then
+            @:PROHIBIT(cfl_target < 0 .or. cfl_target > 1d0)
+            @:PROHIBIT(t_stop <= 0)
+            @:PROHIBIT(t_save <= 0)
+            @:PROHIBIT(t_save > t_stop)
+            @:PROHIBIT(n_start < 0)
+        else
+            @:PROHIBIT(t_step_start < 0)
+            @:PROHIBIT(t_step_stop <= t_step_start)
+            @:PROHIBIT(t_step_save > t_step_stop - t_step_start)
+        end if
     end subroutine s_check_inputs_time_stepping
 
     !> Checks constraints on the finite difference parameters.

examples/2D_ibm_cfl_dt/result.png

@@ -24,5 +24,6 @@ module m_constants
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
     real(kind(0d0)), parameter :: capillary_cutoff = 1e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
     real(kind(0d0)), parameter :: acoustic_spatial_support_width = 2.5d0 !< Spatial support width of acoustic source, used in s_source_spatial
+    real(kind(0d0)), parameter :: dflt_vcfl_dt = 100d0 !< value of vcfl_dt when viscosity is off for computing adaptive timestep size
 
 end module m_constants

src/common/m_checker_common.fpp

@@ -76,6 +76,16 @@ module m_global_parameters
     integer :: t_step_stop   !< Last time-step directory
     integer :: t_step_save   !< Interval between consecutive time-step directory
 
+    !> @name IO options for adaptive time-stepping
+    !> @{
+    logical :: cfl_adap_dt, cfl_const_dt, cfl_dt
+    real(kind(0d0)) :: t_save
+    real(kind(0d0)) :: t_stop
+    real(kind(0d0)) :: cfl_target
+    integer :: n_save
+    integer :: n_start
+    !> @}
+
     ! NOTE: The variables m_root, x_root_cb and x_root_cc contain the grid data
     ! of the defragmented computational domain. They are only used in 1D. For
     ! serial simulations, they are equal to m, x_cb and x_cc, respectively.
@@ -278,6 +288,14 @@ contains
         t_step_stop = dflt_int
         t_step_save = dflt_int
 
+        cfl_adap_dt = .false.
+        cfl_const_dt = .false.
+        cfl_dt = .false.
+        cfl_target = dflt_real
+        t_save = dflt_real
+        n_start = dflt_int
+        t_stop = dflt_real
+
         ! Simulation algorithm parameters
         model_eqns = dflt_int
         num_fluids = dflt_int

src/common/m_constants.fpp

@@ -158,7 +158,8 @@ contains
             & 't_step_start', 't_step_stop', 't_step_save', 'weno_order',      &
             & 'model_eqns', 'num_fluids', 'bc_x%beg', 'bc_x%end', 'bc_y%beg',  &
             & 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'flux_lim', 'format',        &
-            & 'precision', 'fd_order', 'thermal', 'nb', 'relax_model' ]
+            & 'precision', 'fd_order', 'thermal', 'nb', 'relax_model',         &
+            & 'n_start' ]
             call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
@@ -168,7 +169,7 @@ contains
             & 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', &
             & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles', 'qbmm',   &
             & 'polytropic', 'polydisperse', 'file_per_process', 'relax', 'cf_wrt',     &
-            & 'adv_n', 'ib' ]
+            & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt' ]
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
@@ -189,7 +190,7 @@ contains
         end do
 
         #:for VAR in [ 'pref', 'rhoref', 'R0ref', 'poly_sigma', 'Web', 'Ca', &
-            & 'Re_inv', 'sigma' ]
+            & 'Re_inv', 'sigma', 't_save', 't_stop' ]
             call MPI_BCAST(${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
         #:endfor
         call MPI_BCAST(schlieren_alpha(1), num_fluids_max, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)

src/post_process/m_global_parameters.fpp

@@ -74,7 +74,9 @@ subroutine s_read_input_file
             parallel_io, rhoref, pref, bubbles, qbmm, sigR, &
             R0ref, nb, polytropic, thermal, Ca, Web, Re_inv, &
             polydisperse, poly_sigma, file_per_process, relax, &
-            relax_model, cf_wrt, sigma, adv_n, ib
+            relax_model, cf_wrt, sigma, adv_n, ib, &
+            cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
+            cfl_target
 
         ! Inquiring the status of the post_process.inp file
         file_loc = 'post_process.inp'
@@ -102,6 +104,9 @@ subroutine s_read_input_file
             p_glb = p
 
             nGlobal = (m_glb + 1)*(n_glb + 1)*(p_glb + 1)
+
+            if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true.
+
         else
             call s_mpi_abort('File post_process.inp is missing. Exiting ...')
         end if
@@ -142,11 +147,17 @@ subroutine s_perform_time_step(t_step)
 
         integer, intent(inout) :: t_step
         if (proc_rank == 0) then
-            print '(" ["I3"%]  Saving "I8" of "I0" @ t_step = "I0"")', &
-                int(ceiling(100d0*(real(t_step - t_step_start)/(t_step_stop - t_step_start + 1)))), &
-                (t_step - t_step_start)/t_step_save + 1, &
-                (t_step_stop - t_step_start)/t_step_save + 1, &
-                t_step
+            if (cfl_dt) then
+                print '(" ["I3"%]  Saving "I8" of "I0"")', &
+                    int(ceiling(100d0*(real(t_step - n_start)/(n_save)))), &
+                    t_step, n_save
+            else
+                print '(" ["I3"%]  Saving "I8" of "I0" @ t_step = "I0"")', &
+                    int(ceiling(100d0*(real(t_step - t_step_start)/(t_step_stop - t_step_start + 1)))), &
+                    (t_step - t_step_start)/t_step_save + 1, &
+                    (t_step_stop - t_step_start)/t_step_save + 1, &
+                    t_step
+            end if
         end if
 
         ! Populating the grid and conservative variables