Strang Splitting for lagrange bubbles (#827)

Co-authored-by: Diego Vaca <diegovaca@autoreg-6682246.dyn.wpi.edu>
Co-authored-by: Diego Vaca <diegovaca@Mac.lan>
Co-authored-by: Diego Vaca <diegovaca@Diegos-MacBook-Pro-128.local>

Author: 4 people

docs/documentation/case.md

@@ -436,7 +436,7 @@ The effect and use of the source term are assessed by [Schmidmayer et al., 2019]
 - `time_stepper` specifies the order of the Runge-Kutta (RK) time integration scheme that is used for temporal integration in simulation, from the 1st to 5th order by corresponding integer.
 Note that `time_stepper = 3` specifies the total variation diminishing (TVD), third order RK scheme ([Gottlieb and Shu, 1998](references.md)).
 
-- `adap_dt` activates the Strang operator splitting scheme which splits flux and source terms in time marching, and an adaptive time stepping strategy is implemented for the source term. It requires ``bubbles = 'T'``, ``polytropic = 'T'``, ``adv_n = 'T'`` and `time_stepper = 3`.
+- `adap_dt` activates the Strang operator splitting scheme which splits flux and source terms in time marching, and an adaptive time stepping strategy is implemented for the source term. It requires ``bubbles_euler = 'T'``, ``polytropic = 'T'``, ``adv_n = 'T'`` and `time_stepper = 3`. Additionally, it can be used with ``bubbles_lagrange = 'T'`` and `time_stepper = 3`
 
 - `weno_order` specifies the order of WENO scheme that is used for spatial reconstruction of variables by an integer of 1, 3, 5, and 7, that correspond to the 1st, 3rd, 5th, and 7th order, respectively.
 
@@ -790,8 +790,6 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 | `x0`                  | Real    | Reference length                                          |
 | `Thost`               | Real    | Temperature of the surrounding liquid (host)              |
 | `diffcoefvap`         | Real    | Vapor diffusivity in the gas                              |
-| `rkck_adap_dt`        | Logical | Activates the adaptive rkck time stepping algorithm       |
-| `rkck_tolerance`      | Real    | Admissible error truncation tolerance in the rkck stepper  |
 
 - `nBubs_glb` Total number of bubbles. Their initial conditions need to be specified in the ./input/lag_bubbles.dat file. See the example cases for additional information.
 
@@ -805,8 +803,6 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 
 - `massTransfer_model` Activates the mass transfer model at the bubble's interface based on ([Preston et al., 2007](references.md)).
 
-- `rkck_adap_dt` Activates the adaptive 4th/5th order Runge—Kutta–Cash–Karp (RKCK) time-stepping algorithm (requires `time_stepper ==4`). A maximum error between the 4th and 5th order Runge-Kutta-Cash-Karp solutions for the same time step size is calculated. If the error is smaller than a tolerance (`rkck_tolerance`), then the algorithm employs the 5th order solution, while if not, both eulerian/lagrangian variables are re-calculated with a smaller time step size.
-
 ### 10. Velocity Field Setup
 
 | Parameter              | Type    | Description |

examples/2D_lagrange_bubblescreen/case.py

@@ -0,0 +1,170 @@
+#!/usr/bin/env python3
+import math
+import json
+
+# Bubble screen
+# Description: A planar acoustic wave interacts with a bubble cloud
+# in water. The background field is modeled in using an Eulerian framework,
+# while the bubbles are tracked using a Lagrangian framework.
+
+# Reference values for nondimensionalization
+x0 = 1.0e-03  # length - m
+rho0 = 1.0e03  # density - kg/m3
+c0 = 1475.0  # speed of sound - m/s
+p0 = rho0 * c0 * c0  # pressure - Pa
+T0 = 298  # temperature - K
+
+# Host properties (water)
+gamma_host = 2.7466  # Specific heat ratio
+pi_inf_host = 792.02e06  # Stiffness - Pa
+mu_host = 1e-3  # Dynamic viscosity - Pa.s
+c_host = 1475.0  # speed of sound - m/s
+rho_host = 1000  # density kg/m3
+T_host = 298  # temperature K
+
+# Lagrangian bubbles' properties
+R_uni = 8314  # Universal gas constant - J/kmol/K
+MW_g = 28.0  # Molar weight of the gas - kg/kmol
+MW_v = 18.0  # Molar weight of the vapor - kg/kmol
+gamma_g = 1.4  # Specific heat ratio of the gas
+gamma_v = 1.333  # Specific heat ratio of the vapor
+pv = 2350  # Vapor pressure of the host - Pa
+cp_g = 1.0e3  # Specific heat of the gas - J/kg/K
+cp_v = 2.1e3  # Specific heat of the vapor - J/kg/K
+k_g = 0.025  # Thermal conductivity of the gas - W/m/K
+k_v = 0.02  # Thermal conductivity of the vapor - W/m/K
+diffVapor = 2.5e-5  # Diffusivity coefficient of the vapor - m2/s
+sigBubble = 0.069  # Surface tension of the bubble - N/m
+mu_g = 1.48e-5
+
+# Acoustic source properties
+patm = 101325.0  # Atmospheric pressure - Pa
+pamp = 1.0e5  # Amplitude of the acoustic source - Pa
+freq = 300e03  # Source frequency - Hz
+wlen = c_host / freq  # Wavelength - m
+
+# Domain and time set up
+
+xb = -12.0e-3  # Domain boundaries - m (x direction)
+xe = 12.0e-3
+yb = -2.5e-3  # Domain boundaries - m (y direction)
+ye = 2.5e-3
+z_virtual = 5.0e-3  # Virtual depth (z direction)
+
+Nx = 240  # number of elements into x direction
+Ny = 50  # number of elements into y direction
+
+dt = 7.5e-9  # constant time-step - sec
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": xb / x0,
+            "x_domain%end": xe / x0,
+            "y_domain%beg": yb / x0,
+            "y_domain%end": ye / x0,
+            "stretch_y": "F",
+            "stretch_x": "F",
+            "m": Nx,
+            "n": Ny,
+            "p": 0,
+            "dt": dt * (c0 / x0),
+            "t_step_start": 0,
+            "t_step_stop": 3000,
+            "t_step_save": 500,
+            # Simulation Algorithm Parameters
+            "model_eqns": 2,
+            "time_stepper": 3,
+            "num_fluids": 2,
+            "num_patches": 1,
+            "viscous": "T",
+            "mpp_lim": "F",
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -6,
+            "bc_x%end": -6,
+            "bc_y%beg": -1,
+            "bc_y%end": -1,
+            # Acoustic source
+            "acoustic_source": "T",
+            "num_source": 1,
+            "acoustic(1)%support": 2,
+            "acoustic(1)%pulse": 1,
+            "acoustic(1)%npulse": 1,
+            "acoustic(1)%mag": pamp / p0,
+            "acoustic(1)%wavelength": wlen / x0,
+            "acoustic(1)%length": 2 * (ye - yb) / x0,
+            "acoustic(1)%loc(1)": -7.0e-03 / x0,
+            "acoustic(1)%loc(2)": 0.0,
+            "acoustic(1)%dir": 0.0,
+            "acoustic(1)%delay": 0.0,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Patch 1: Water (left)
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.0,
+            "patch_icpp(1)%y_centroid": 0.0,
+            "patch_icpp(1)%length_x": 2 * (xe - xb) / x0,
+            "patch_icpp(1)%length_y": 2 * (ye - yb) / x0,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": patm / p0,
+            "patch_icpp(1)%alpha_rho(1)": rho_host / rho0,
+            "patch_icpp(1)%alpha_rho(2)": 0.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "patch_icpp(1)%alpha(2)": 0.0,
+            # Lagrangian Bubbles
+            "bubbles_lagrange": "T",
+            "bubble_model": 2,  # Keller-Miksis model
+            "lag_params%nBubs_glb": 1194,  # Number of bubbles
+            "lag_params%solver_approach": 2,
+            "lag_params%cluster_type": 2,
+            "lag_params%pressure_corrector": "T",
+            "lag_params%smooth_type": 1,
+            "lag_params%heatTransfer_model": "T",
+            "lag_params%massTransfer_model": "T",
+            "lag_params%epsilonb": 1.0,
+            "lag_params%valmaxvoid": 0.9,
+            "lag_params%write_bubbles": "F",
+            "lag_params%write_bubbles_stats": "F",
+            "lag_params%c0": c0,
+            "lag_params%rho0": rho0,
+            "lag_params%T0": T0,
+            "lag_params%x0": x0,
+            "lag_params%diffcoefvap": diffVapor,
+            "lag_params%Thost": T_host,
+            "lag_params%charwidth": z_virtual / x0,
+            # Fluids Physical Parameters
+            # Host medium
+            "fluid_pp(1)%gamma": 1.0 / (gamma_host - 1.0),
+            "fluid_pp(1)%pi_inf": gamma_host * (pi_inf_host / p0) / (gamma_host - 1.0),
+            "fluid_pp(1)%Re(1)": 1.0 / (mu_host / (rho0 * c0 * x0)),
+            "fluid_pp(1)%mul0": mu_host,
+            "fluid_pp(1)%ss": sigBubble,
+            "fluid_pp(1)%pv": pv,
+            "fluid_pp(1)%gamma_v": gamma_v,
+            "fluid_pp(1)%M_v": MW_v,
+            "fluid_pp(1)%k_v": k_v,
+            "fluid_pp(1)%cp_v": cp_v,
+            # Bubble gas state
+            "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0e00,
+            "fluid_pp(2)%Re(1)": 1.0 / (mu_g / (rho0 * c0 * x0)),
+            "fluid_pp(2)%gamma_v": gamma_g,
+            "fluid_pp(2)%M_v": MW_g,
+            "fluid_pp(2)%k_v": k_g,
+            "fluid_pp(2)%cp_v": cp_g,
+        }
+    )
+)

examples/2D_lagrange_bubblescreen/input/README.txt

@@ -0,0 +1,5 @@
+
+The user input file 'input/lag_bubbles.dat' contains the initial conditions of the lagrangian bubbles.
+Each row represents the initial state of one specific bubble, which are:
+
+xPosition/x0    yPosition/x0    zPosition/x0    xVel/c0     yVel/c0     zVel/c0     radius/x0       interfaceVelocity/c0

examples/3D_lagrange_shbubcollapse/case.py

@@ -80,6 +80,7 @@
             "n": Ny,
             "p": Nz,
             "dt": round(dt * c0 / x0, 6),
+            "adap_dt": "T",
             "n_start": 0,
             "t_save": saveTime * (c0 / x0),
             "t_stop": stopTime * (c0 / x0),
@@ -89,7 +90,7 @@
             "num_patches": 1,
             "mpp_lim": "F",
             "viscous": "T",
-            "time_stepper": 4,  # 4th/5th RKCK
+            "time_stepper": 3,
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "mapped_weno": "T",
@@ -141,8 +142,6 @@
             # Lagrangian Bubbles
             "bubbles_lagrange": "T",
             "bubble_model": 2,  # Keller-Miksis model
-            "rkck_adap_dt": "T",  # Activate adaptive time stepper
-            "rkck_tolerance": 1.0e-05,
             "lag_params%nBubs_glb": 1,
             "lag_params%solver_approach": 2,  # Two-way coupled
             "lag_params%cluster_type": 2,

src/common/m_checker_common.fpp

@@ -62,7 +62,7 @@ contains
         !! Called by s_check_inputs_common for simulation and post-processing
     subroutine s_check_inputs_time_stepping
         if (cfl_dt) then
-            @:PROHIBIT((cfl_target < 0 .or. cfl_target > 1._wp) .and. .not. rkck_adap_dt)
+            @:PROHIBIT(cfl_target < 0 .or. cfl_target > 1._wp)
             @:PROHIBIT(t_stop <= 0)
             @:PROHIBIT(t_save <= 0)
             @:PROHIBIT(t_save > t_stop)

src/common/m_constants.fpp

@@ -50,34 +50,16 @@ module m_constants
     integer, parameter :: mapCells = 3 !< Number of cells around the bubble where the smoothening function will have effect
     real(wp), parameter :: R_uni = 8314._wp ! Universal gas constant - J/kmol/K
 
-    ! RKCK constants
-    integer, parameter :: num_ts_rkck = 6 !< Number of time-stages in the RKCK stepper
-    ! RKCK 4th/5th time stepper coefficients based on Cash J. and Karp A. (1990)
-    real(wp), parameter :: rkck_c1 = 0._wp, rkck_c2 = 0.2_wp, rkck_c3 = 0.3_wp, rkck_c4 = 0.6_wp, &     ! c1 c2 c3 c4 c5 c6
-                           rkck_c5 = 1._wp, rkck_c6 = 0.875_wp
-    real(wp), dimension(6), parameter :: rkck_coef1 = (/0.2_wp, 0._wp, 0._wp, 0._wp, 0._wp, 0._wp/)     ! a21
-    real(wp), dimension(6), parameter :: rkck_coef2 = (/3._wp/40._wp, 9._wp/40._wp, 0._wp, 0._wp, &     ! a31 a32
-                                                        0._wp, 0._wp/)
-    real(wp), dimension(6), parameter :: rkck_coef3 = (/0.3_wp, -0.9_wp, 1.2_wp, 0._wp, 0._wp, 0._wp/)  ! a41 a42 a43
-    real(wp), dimension(6), parameter :: rkck_coef4 = (/-11._wp/54._wp, 2.5_wp, -70._wp/27._wp, &       ! a51 a52 a53 a54
-                                                        35._wp/27._wp, 0._wp, 0._wp/)
-    real(wp), dimension(6), parameter :: rkck_coef5 = (/1631._wp/55296._wp, 175._wp/512._wp, &          ! a61 a62 a63 a64 a65
-                                                        575._wp/13824._wp, 44275._wp/110592._wp, &
-                                                        253._wp/4096._wp, 0._wp/)
-    real(wp), dimension(6), parameter :: rkck_coef6 = (/37._wp/378._wp, 0._wp, 250._wp/621._wp, &       ! b1 b2 b3 b4 b5 b6
-                                                        125._wp/594._wp, 0._wp, 512._wp/1771._wp/)
-    real(wp), dimension(6), parameter :: rkck_coefE = (/37._wp/378._wp - 2825._wp/27648._wp, 0._wp, &   ! er1 er2 er3 er4 er5 er6 (4th/5th error)
-                                                        250._wp/621._wp - 18575._wp/48384._wp, &
-                                                        125._wp/594._wp - 13525._wp/55296._wp, &
-                                                        -277._wp/14336._wp, 512._wp/1771._wp - 0.25_wp/)
-    ! Adaptive rkck constants
-    real(wp), parameter :: verysmall_dt = 1e-14_wp !< Very small dt, stop stepper
-    real(wp), parameter :: SAFETY = 0.9_wp !< Next dt will be maximum 0.9*dt if truncation error is above tolerance.
-    real(wp), parameter :: RNDDEC = 1e8_wp !< Round calculated dt (16th digit) to avoid the inclusion of random decimals
-    real(wp), parameter :: PSHRNK = -0.25_wp !< Factor to reduce dt when truncation error above tolerance
-    real(wp), parameter :: SHRNKDT = 0.5_wp !< Factor to reduce dt due to negative bubble radius
-    real(wp), parameter :: ERRCON = 1.89e-4_wp !< Limit to slightly increase dt when truncation error is between ERRCON and 1
-    real(wp), parameter :: PGROW = -0.2_wp !< Factor to increase dt when truncation error is between ERRCON and 1
+    ! Strang Splitting constants
+    real(wp), parameter :: dflt_adap_dt_tol = 1e-4_wp !< Default tolerance for adaptive step size
+    integer, parameter :: adap_dt_max_iters = 100 !< Maximum number of iterations
+    ! Constants of the algorithm described by Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary Differential Equations I, Chapter II.4
+    ! to choose the initial time step size for the adaptive time stepping routine
+    real(wp), parameter :: threshold_first_guess = 1e-5_wp
+    real(wp), parameter :: threshold_second_guess = 1e-15_wp
+    real(wp), parameter :: scale_first_guess = 1e-3_wp
+    real(wp), parameter :: scale_guess = 1e-2_wp
+    real(wp), parameter :: small_guess = 1e-6_wp
 
     ! Relativity
     integer, parameter :: relativity_cons_to_prim_max_iter = 100

src/post_process/m_global_parameters.fpp

@@ -314,7 +314,7 @@ module m_global_parameters
 
     !> @name Lagrangian bubbles
     !> @{
-    logical :: bubbles_lagrange, rkck_adap_dt
+    logical :: bubbles_lagrange
     !> @}
 
     real(wp) :: Bx0 !< Constant magnetic field in the x-direction (1D)
@@ -443,7 +443,6 @@ contains
 
         ! Lagrangian bubbles modeling
         bubbles_lagrange = .false.
-        rkck_adap_dt = .false.
 
         ! IBM
         num_ibs = dflt_int

src/post_process/m_mpi_proxy.fpp

@@ -170,8 +170,7 @@ contains
             & 'file_per_process', 'relax', 'cf_wrt',                           &
             & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt',          &
             & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange',        &
-            & 'rkck_adap_dt', 'output_partial_domain', 'relativity',           &
-            & 'cont_damage' ]
+            & 'output_partial_domain', 'relativity', 'cont_damage' ]
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

src/post_process/m_start_up.f90

@@ -83,7 +83,7 @@ subroutine s_read_input_file
             polydisperse, poly_sigma, file_per_process, relax, &
             relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, &
             cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
-            cfl_target, surface_tension, bubbles_lagrange, rkck_adap_dt, &
+            cfl_target, surface_tension, bubbles_lagrange, &
             sim_data, hyperelasticity, Bx0, relativity, cont_damage
 
         ! Inquiring the status of the post_process.inp file
@@ -113,7 +113,7 @@ subroutine s_read_input_file
 
             nGlobal = (m_glb + 1)*(n_glb + 1)*(p_glb + 1)
 
-            if (cfl_adap_dt .or. cfl_const_dt .or. rkck_adap_dt) cfl_dt = .true.
+            if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true.
 
         else
             call s_mpi_abort('File post_process.inp is missing. Exiting.')

src/pre_process/m_global_parameters.fpp

@@ -267,10 +267,6 @@ module m_global_parameters
     integer :: chemxb, chemxe
     !> @}
 
-    !> @ lagrangian solver parameters
-    logical :: rkck_adap_dt
-    !> @}
-
     integer, allocatable, dimension(:, :, :) :: logic_grid
 
     type(pres_field) :: pb
@@ -545,9 +541,6 @@ contains
             fluid_pp(i)%G = 0._wp
         end do
 
-        ! Lagrangian solver
-        rkck_adap_dt = .false.
-
         Bx0 = dflt_real
 
     end subroutine s_assign_default_values_to_user_inputs