Fix to Courant-based timestepping

Author: mbahlali

ICFERST/src/Extract_From_State.F90

@@ -2209,9 +2209,22 @@ subroutine Adaptive_NonLinear(Mdims, packed_state, reference_field, its, itime,&
       call real_from_python(pyfunc, acctim, min_ts)
     end if
     !Now check the CFL ones just in case that approach is being used in combination with the non-linear case
-    call get_option('/timestepping/adaptive_timestep/minimum_timestep', auxR, stat= auxI)
+    if (have_option('/timestepping/adaptive_timestep/minimum_timestep/constant')) then
+      call get_option('/timestepping/adaptive_timestep/minimum_timestep/constant', auxR, stat= auxI)
+    else if (have_option('/timestepping/adaptive_timestep/minimum_timestep/python')) then
+      call get_option('/timestepping/adaptive_timestep/minimum_timestep/python', pyfunc)
+      call real_from_python(pyfunc, acctim, auxR)
+    end if
+    
     if (auxI == 0) min_ts = min(min_ts, auxR)
-    call get_option('/timestepping/adaptive_timestep/maximum_timestep', auxR, stat= auxI)
+
+    if (have_option('/timestepping/adaptive_timestep/maximum_timestep/constant')) then
+      call get_option('/timestepping/adaptive_timestep/maximum_timestep/constant', auxR, stat= auxI)
+    else if (have_option('/timestepping/adaptive_timestep/maximum_timestep/python')) then
+      call get_option('/timestepping/adaptive_timestep/maximum_timestep/python', pyfunc)
+      call real_from_python(pyfunc, acctim, auxR)
+    end if
+
     if (auxI == 0) max_ts = min(max_ts, auxR)
 
     !Ensure that even adapting the time, the final time is matched

ICFERST/src/Multiphase_TimeLoop.F90

@@ -254,6 +254,8 @@ subroutine MultiFluids_SolveTimeLoop( state, &
         double precision, ALLOCATABLE, dimension(:,:) :: concetration_phreeqc
         real :: total_mass_metal_before_adapt, total_mass_metal_after_adapt, total_mass_metal_after_correction, total_mass_metal_after_bound
         logical :: viscosity_EOS
+        character(len=PYTHON_FUNC_LEN) :: pyfunc
+        
 #ifdef HAVE_ZOLTAN
       real(zoltan_float) :: ver
       integer(zoltan_int) :: ierr
@@ -1021,8 +1023,19 @@ subroutine MultiFluids_SolveTimeLoop( state, &
                 nonlinear_dt = dt!To use if also the nonlinear adapt time-step is on
                 c = -66.6 ; minc = 0. ; maxc = 66.e6 ; ic = 1.1!66.e6
                 call get_option( '/timestepping/adaptive_timestep/requested_cfl', rc )
-                call get_option( '/timestepping/adaptive_timestep/minimum_timestep', minc, stat, default = 0.)
-                call get_option( '/timestepping/adaptive_timestep/maximum_timestep', maxc, stat, default = 66.e6)
+
+                call get_option('/timestepping/adaptive_timestep/minimum_timestep/constant', minc, stat, default = 0.)
+                if (have_option('/timestepping/adaptive_timestep/minimum_timestep/python')) then
+                    call get_option('/timestepping/adaptive_timestep/minimum_timestep/python', pyfunc)
+                    call real_from_python(pyfunc, acctim, minc)
+                end if
+
+                call get_option('/timestepping/adaptive_timestep/maximum_timestep/constant', maxc, stat, default = 66.e6)
+                if (have_option('/timestepping/adaptive_timestep/maximum_timestep/python')) then
+                    call get_option('/timestepping/adaptive_timestep/maximum_timestep/python', pyfunc)
+                    call real_from_python(pyfunc, acctim, maxc)
+                end if
+
                 call get_option( '/timestepping/adaptive_timestep/increase_tolerance', ic, stat, default = 1.1)
                 !For porous media we need to use the Courant number obtained in cv_assemb
                 if (is_porous_media) then

femtools/Adaptive_Timestepping.F90

@@ -31,7 +31,7 @@ module adaptive_timestepping
   !!< Contains new style adaptive timestepping routines
 
   use fldebug
-  use global_parameters, only : OPTION_PATH_LEN, FIELD_NAME_LEN
+  use global_parameters, only : OPTION_PATH_LEN, FIELD_NAME_LEN, PYTHON_FUNC_LEN
   use spud
   use unittest_tools
   use fields
@@ -70,6 +70,8 @@ subroutine calc_cflnumber_field_based_dt(state, dt, force_calculation)
     type(scalar_field), pointer :: cflnumber_field
     type(vector_field), pointer :: velocity_field
     type(mesh_type), pointer :: mesh
+    character(len=PYTHON_FUNC_LEN) :: pyfunc
+    real :: current_time
 
     ewrite(1, *) "In calc_cflnumber_field_based_dt"
 
@@ -79,9 +81,19 @@ subroutine calc_cflnumber_field_based_dt(state, dt, force_calculation)
       lforce_calculation = .false.
     end if
 
+    call get_option("/timestepping/current_time", current_time)
+
     call get_option(base_path // "/requested_cfl", max_cfl)
-    call get_option(base_path // "/minimum_timestep", min_dt, default = tiny(0.0))
-    call get_option(base_path // "/maximum_timestep", max_dt, default = huge(0.0))
+    call get_option(base_path // "/minimum_timestep/constant", min_dt, default = tiny(0.0))
+    if (have_option(base_path // '/minimum_timestep/python')) then
+      call get_option(base_path // '/minimum_timestep/python', pyfunc)
+      call real_from_python(pyfunc, current_time, min_dt)
+    end if
+    call get_option(base_path // "/maximum_timestep/constant", max_dt, default = huge(0.0))
+    if (have_option(base_path // '/maximum_timestep/python')) then
+      call get_option(base_path // '/maximum_timestep/python', pyfunc)
+      call real_from_python(pyfunc, current_time, max_dt)
+    end if
     call get_option(base_path // "/increase_tolerance", increase_tolerance, default = huge(0.0) * epsilon(0.0))
     ! what type of cfl number are we using?
     call get_option(base_path//"/courant_number[0]/name", cfl_type)
@@ -187,6 +199,8 @@ subroutine adaptive_timestepping_check_options
     character(len = OPTION_PATH_LEN) :: base_path
     integer :: stat
     real :: max_cfl, max_dt, min_dt, increase_tolerance
+    character(len=PYTHON_FUNC_LEN) :: pyfunc
+    real :: current_time
 
     base_path = "/timestepping/adaptive_timestep"
 
@@ -205,7 +219,15 @@ subroutine adaptive_timestepping_check_options
       FLExit("Maximum adaptive timestepping CFL number must be positive")
     end if
 
-    call get_option(trim(base_path) // "/minimum_timestep", min_dt, stat)
+    call get_option("/timestepping/current_time", current_time)
+
+    if (have_option(trim(base_path) // '/minimum_timestep/constant')) then
+      call get_option(trim(base_path) // "/minimum_timestep/constant", min_dt, stat)
+    else if (have_option(trim(base_path) // '/minimum_timestep/python')) then
+      call get_option(trim(base_path) // '/minimum_timestep/python', pyfunc)
+      call real_from_python(pyfunc, current_time, min_dt)
+    end if
+
     if(stat == SPUD_NO_ERROR) then
       if(min_dt < 0.0) then
         FLExit("Minimum timestep size cannot be negative")

main/Fluids.F90

@@ -36,7 +36,7 @@ module fluids_module
        simulation_start_time, &
        simulation_start_cpu_time, &
        simulation_start_wall_time, &
-       topology_mesh_name, FIELD_NAME_LEN, is_porous_media
+       topology_mesh_name, FIELD_NAME_LEN, is_porous_media, PYTHON_FUNC_LEN
   use futils, only: int2str
   use reference_counting, only: print_references
   use parallel_tools
@@ -912,6 +912,7 @@ subroutine update_state_post_adapt(state, metric_tensor, dt, sub_state, nonlinea
     real, intent(inout) :: dt
     integer, intent(inout) :: nonlinear_iterations, nonlinear_iterations_adapt
     type(state_type), dimension(:), pointer :: sub_state
+    character(len=PYTHON_FUNC_LEN) :: pyfunc
 
     ! Overwrite the number of nonlinear iterations if the option is switched on
     if(have_option("/timestepping/nonlinear_iterations/nonlinear_iterations_at_adapt")) then
@@ -939,8 +940,12 @@ subroutine update_state_post_adapt(state, metric_tensor, dt, sub_state, nonlinea
     call enforce_discrete_properties(state)
 
     if (have_option("/mesh_adaptivity/hr_adaptivity/adaptive_timestep_at_adapt")) then
-        if (have_option("/timestepping/adaptive_timestep/minimum_timestep")) then
-            call get_option("/timestepping/adaptive_timestep/minimum_timestep", dt)
+        if (have_option('/timestepping/adaptive_timestep/minimum_timestep/constant')) then
+            call get_option('/timestepping/adaptive_timestep/minimum_timestep/constant', dt)
+            call set_option("/timestepping/timestep", dt)
+        else if (have_option('/timestepping/adaptive_timestep/minimum_timestep/python')) then
+            call get_option('/timestepping/adaptive_timestep/minimum_timestep/python', pyfunc)
+            call real_from_python(pyfunc, current_time, dt)
             call set_option("/timestepping/timestep", dt)
         else
             ewrite(-1,*) "Warning: you have adaptive timestep adjustment after &&