Merge branch 'muscl' of https://github.com/okBrian/MFC into muscl

okBrian · okBrian · commit 6695ef541872 · 2025-07-28T17:30:16.000-04:00
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -380,6 +380,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `time_stepper`         | Integer | Runge--Kutta order [1-3] |
 | `adap_dt`              | Logical | Strang splitting scheme with adaptive time stepping |
 | `recon_type`           | Integer | Reconstruction Type: [1] WENO; [2] MUSCL |
+| `adap_dt_tol`          | Real    | Tolerance for adaptive time stepping in Strang splitting scheme|
+| `adap_dt_max_iters`    | Integer | Max iteration for adaptive time stepping in Strang splitting scheme |
 | `weno_order`	         | Integer | WENO order [1,3,5] |
 | `weno_eps`	         | Real    | WENO perturbation (avoid division by zero) |
 | `mapped_weno`	         | Logical | WENO-M (WENO with mapping of nonlinear weights) |
@@ -459,7 +461,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_euler = 'T'``, ``polytropic = 'T'``, ``adv_n = 'T'`` and `time_stepper = 3`. Additionally, it can be used with ``bubbles_lagrange = '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`. `adap_dt_tol` and `adap_dt_max_iters` are 1e-4 and 100, respectively, by default.
 
 - `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.
 
@@ -756,7 +758,6 @@ Implementation of the parameters into the model follow [Ando (2010)](references.
 | `polydisperse`   | Logical	| Polydispersity in equilibrium bubble radius R0 |
 | `nb` 			     | Integer	| Number of bins: [1] Monodisperse; [$>1$] Polydisperse |
 | `poly_sigma` 	       | Real 		|	Standard deviation for probability density function of polydisperse bubble populations |
-| `R0_type` 	       | Integer 		|	Quadrature rule for probability density function of polydisperse bubble populations |
 | `Ca` 			     | Real		  | Cavitation number |
 | `Web` 			   | Real		  | Weber number |
 | `Re_inv` 		   | Real		  | Inverse Reynolds number |
@@ -774,15 +775,12 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 - `thermal` specifies a model for heat transfer across the bubble interface by an integer from 1 through 3.
 `thermal = 1`, `2`, and `3` correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on [Preston et al., 2007](references.md), respectively.
 
-- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius.
+- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius. Simpson's rule is used for integrating the log-normal PDF of equilibrium bubble radius for polydisperse populations.
 
 - `R0ref` specifies the reference bubble radius.
 
 - `nb` specifies the number of discrete bins that define the probability density function (PDF) of the equilibrium bubble radius.
 
-- `R0_type` specifies the quadrature rule for integrating the log-normal PDF of equilibrium bubble radius for polydisperse populations.
-`R0_type = 1` corresponds to Simpson's rule.
-
 - `poly_sigma` specifies the standard deviation of the log-normal PDF of equilibrium bubble radius for polydisperse populations.
 
 - `Ca`, `Web`, and `Re_inv` respectively specify the Cavitation number, Weber number, and the inverse Reynolds number that characterize the offset of the gas pressure from the vapor pressure, surface tension, and liquid viscosity when the polytropic gas compression model is used.
diff --git a/examples/0D_bubblecollapse_adap/case.py b/examples/0D_bubblecollapse_adap/case.py
@@ -118,6 +118,7 @@
             "adv_n": "T",
             # adap_dt
             "adap_dt": "T",
+            "adap_dt_max_iters": 200,
             # Gas compression model
             "polytropic": "T",
             "thermal": 1,
diff --git a/examples/1D_poly_bubscreen/case.py b/examples/1D_poly_bubscreen/case.py
@@ -157,7 +157,6 @@
             "bubbles_euler": "T",
             "bubble_model": 2,
             "polytropic": "T",
-            "R0_type": 1,
             "thermal": 3,
             "R0ref": myr0,
             "nb": 1,
diff --git a/examples/1D_qbmm/case.py b/examples/1D_qbmm/case.py
@@ -171,7 +171,6 @@
             "bubble_model": 2,
             "polytropic": "F",
             "polydisperse": "T",
-            "R0_type": 1,
             "poly_sigma": 0.3,
             "thermal": 3,
             "R0ref": myr0,
diff --git a/examples/2D_bubbly_steady_shock/case.py b/examples/2D_bubbly_steady_shock/case.py
@@ -191,7 +191,6 @@
             "bubble_model": 2,
             "polytropic": "T",
             "polydisperse": "F",
-            "R0_type": 1,
             "poly_sigma": 0.3,
             "thermal": 3,
             "R0ref": x0,
diff --git a/examples/2D_whale_bubble_annulus/case.py b/examples/2D_whale_bubble_annulus/case.py
@@ -152,7 +152,6 @@
             "bubbles_euler": "T",
             "bubble_model": 3,
             "polytropic": "T",
-            "R0_type": 1,
             "thermal": 3,
             "R0ref": myr0,
             "nb": 1,
diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
@@ -64,7 +64,8 @@ module m_constants
 
     ! Strang Splitting constants
     real(wp), parameter :: dflt_adap_dt_tol = 1.e-4_wp !< Default tolerance for adaptive step size
-    integer, parameter :: adap_dt_max_iters = 100 !< Maximum number of iterations
+    integer, parameter :: dflt_adap_dt_max_iters = 100 !< Default max iteration for adaptive step size
+
     ! 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 = 1.e-5_wp
diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp
@@ -284,7 +284,7 @@ module m_global_parameters
     integer :: nb
     real(wp) :: R0ref
     real(wp) :: Ca, Web, Re_inv
-    real(wp), dimension(:), allocatable :: weight, R0, V0
+    real(wp), dimension(:), allocatable :: weight, R0
     logical :: bubbles_euler
     logical :: qbmm
     logical :: polytropic
@@ -560,7 +560,7 @@ contains
 
                 allocate (bub_idx%rs(nb), bub_idx%vs(nb))
                 allocate (bub_idx%ps(nb), bub_idx%ms(nb))
-                allocate (weight(nb), R0(nb), V0(nb))
+                allocate (weight(nb), R0(nb))
 
                 if (qbmm) then
                     allocate (bub_idx%moms(nb, nmom))
@@ -592,11 +592,7 @@ contains
                 if (nb == 1) then
                     weight(:) = 1._wp
                     R0(:) = 1._wp
-                    V0(:) = 0._wp
-                else if (nb > 1) then
-                    !call s_simpson
-                    V0(:) = 0._wp
-                else
+                else if (nb < 1) then
                     stop 'Invalid value of nb'
                 end if
 
@@ -663,7 +659,7 @@ contains
 
                 allocate (bub_idx%rs(nb), bub_idx%vs(nb))
                 allocate (bub_idx%ps(nb), bub_idx%ms(nb))
-                allocate (weight(nb), R0(nb), V0(nb))
+                allocate (weight(nb), R0(nb))
 
                 do i = 1, nb
                     if (polytropic .neqv. .true.) then
@@ -684,10 +680,7 @@ contains
                 if (nb == 1) then
                     weight(:) = 1._wp
                     R0(:) = 1._wp
-                    V0(:) = 0._wp
-                else if (nb > 1) then
-                    V0(:) = 0._wp
-                else
+                else if (nb < 1) then
                     stop 'Invalid value of nb'
                 end if
 
diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp
@@ -405,8 +405,8 @@ contains
         ! Bubbles euler variables
         if (bubbles_euler) then
             do i = 1, nb
-                muR = R0(i)*patch_icpp(smooth_patch_id)%r0 ! = R0(i)
-                muV = V0(i)*patch_icpp(smooth_patch_id)%v0 ! = 0
+                muR = R0(i)*patch_icpp(smooth_patch_id)%r0
+                muV = patch_icpp(smooth_patch_id)%v0
                 if (qbmm) then
                     ! Initialize the moment set
                     if (dist_type == 1) then
@@ -610,8 +610,8 @@ contains
         ! Smoothed bubble variables
         if (bubbles_euler) then
             do i = 1, nb
-                muR = R0(i)*patch_icpp(patch_id)%r0 ! = 1*R0(i)
-                muV = V0(i)*patch_icpp(patch_id)%v0 ! = 1*V0(i)
+                muR = R0(i)*patch_icpp(patch_id)%r0
+                muV = patch_icpp(patch_id)%v0
                 if (qbmm) then
                     ! Initialize the moment set
                     if (dist_type == 1) then
@@ -630,12 +630,6 @@ contains
                         q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
                     end if
                 else
-                    ! q_prim_vf(bub_idx%rs(i))%sf(j,k,l) = &
-                    !     (eta * R0(i)*patch_icpp(patch_id)%r0 &
-                    !     + (1._wp-eta)*orig_prim_vf(bub_idx%rs(i)))
-                    ! q_prim_vf(bub_idx%vs(i))%sf(j,k,l) = &
-                    !     (eta * V0(i)*patch_icpp(patch_id)%v0 &
-                    !     + (1._wp-eta)*orig_prim_vf(bub_idx%vs(i)))
                     q_prim_vf(bub_idx%rs(i))%sf(j, k, l) = muR
                     q_prim_vf(bub_idx%vs(i))%sf(j, k, l) = muV
 
diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp
@@ -137,11 +137,10 @@ contains
     end subroutine s_check_inputs_grid_stretching
 
     !> Checks constraints on the QBMM and polydisperse bubble parameters
-        !! (qbmm, polydisperse, dist_type, rhoRV, and R0_type)
+        !! (qbmm, polydisperse, dist_type and rhoRV)
     impure subroutine s_check_inputs_qbmm_and_polydisperse
         @:PROHIBIT(qbmm .and. dist_type == dflt_int, "dist_type must be set if using QBMM")
         @:PROHIBIT(qbmm .and. dist_type /= 1 .and. rhoRV > 0._wp, "rhoRV cannot be used with dist_type != 1")
-        @:PROHIBIT(polydisperse .and. R0_type == dflt_int, "R0 type must be set if using Polydisperse")
     end subroutine s_check_inputs_qbmm_and_polydisperse
 
     !> Checks constraints on initial partial density perturbation
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
@@ -219,7 +219,7 @@ module m_global_parameters
     integer :: nb
     real(wp) :: R0ref
     real(wp) :: Ca, Web, Re_inv
-    real(wp), dimension(:), allocatable :: weight, R0, V0
+    real(wp), dimension(:), allocatable :: weight, R0
     logical :: bubbles_euler
     logical :: qbmm      !< Quadrature moment method
     integer :: nmom  !< Number of carried moments
@@ -255,7 +255,6 @@ module m_global_parameters
     real(wp) :: gamma_m, gamma_n, mu_n
     real(wp) :: poly_sigma
     integer :: dist_type !1 = binormal, 2 = lognormal-normal
-    integer :: R0_type   !1 = simpson
     !> @}
 
     !> @name Surface Tension Modeling
@@ -498,7 +497,6 @@ contains
         sigV = dflt_real
         rhoRV = 0._wp
         dist_type = dflt_int
-        R0_type = dflt_int
 
         R_n = dflt_real
         R_v = dflt_real
@@ -648,7 +646,7 @@ contains
                     sys_size = n_idx
                 end if
 
-                allocate (weight(nb), R0(nb), V0(nb))
+                allocate (weight(nb), R0(nb))
                 allocate (bub_idx%rs(nb), bub_idx%vs(nb))
                 allocate (bub_idx%ps(nb), bub_idx%ms(nb))
 
@@ -689,11 +687,7 @@ contains
                 if (nb == 1) then
                     weight(:) = 1._wp
                     R0(:) = 1._wp
-                    V0(:) = 1._wp
-                else if (nb > 1) then
-                    V0(:) = 1._wp
-                    !R0 and weight initialized in s_simpson
-                else
+                else if (nb < 1) then
                     stop 'Invalid value of nb'
                 end if
 
@@ -768,7 +762,7 @@ contains
 
                 allocate (bub_idx%rs(nb), bub_idx%vs(nb))
                 allocate (bub_idx%ps(nb), bub_idx%ms(nb))
-                allocate (weight(nb), R0(nb), V0(nb))
+                allocate (weight(nb), R0(nb))
 
                 do i = 1, nb
                     if (.not. polytropic) then
@@ -789,10 +783,7 @@ contains
                 if (nb == 1) then
                     weight(:) = 1._wp
                     R0(:) = 1._wp
-                    V0(:) = 0._wp
-                else if (nb > 1) then
-                    V0(:) = 1._wp
-                else
+                else if (nb < 1) then
                     stop 'Invalid value of nb'
                 end if
 
diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp
@@ -45,7 +45,7 @@ contains
             & 'loops_x', 'loops_y', 'loops_z', 'model_eqns', 'num_fluids',     &
             & 'weno_order', 'precision', 'perturb_flow_fluid', &
             & 'perturb_sph_fluid', 'num_patches', 'thermal', 'nb', 'dist_type',&
-            & 'R0_type', 'relax_model', 'num_ibs', 'n_start', 'elliptic_smoothing_iters', &
+            & 'relax_model', 'num_ibs', 'n_start', 'elliptic_smoothing_iters', &
             & 'num_bc_patches', 'mixlayer_perturb_nk', 'recon_type', 'muscl_order']
             call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
         #:endfor
diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp
@@ -142,7 +142,7 @@ contains
             rhoref, pref, bubbles_euler, R0ref, nb, &
             polytropic, thermal, Ca, Web, Re_inv, &
             polydisperse, poly_sigma, qbmm, &
-            sigR, sigV, dist_type, rhoRV, R0_type, &
+            sigR, sigV, dist_type, rhoRV, &
             file_per_process, relax, relax_model, &
             palpha_eps, ptgalpha_eps, ib, num_ibs, patch_ib, &
             sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, &
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
@@ -410,8 +410,7 @@ module m_global_parameters
 
     real(wp), dimension(:), allocatable :: weight !< Simpson quadrature weights
     real(wp), dimension(:), allocatable :: R0     !< Bubble sizes
-    real(wp), dimension(:), allocatable :: V0     !< Bubble velocities
-    $:GPU_DECLARE(create='[weight,R0,V0]')
+    $:GPU_DECLARE(create='[weight,R0]')
 
     logical :: bubbles_euler      !< Bubbles euler on/off
     logical :: polytropic   !< Polytropic  switch
@@ -421,7 +420,8 @@ module m_global_parameters
     logical :: adv_n        !< Solve the number density equation and compute alpha from number density
     logical :: adap_dt      !< Adaptive step size control
     real(wp) :: adap_dt_tol !< Tolerance to control adaptive step size
-    $:GPU_DECLARE(create='[adv_n,adap_dt,adap_dt_tol]')
+    integer :: adap_dt_max_iters !< Maximum number of iterations
+    $:GPU_DECLARE(create='[adv_n,adap_dt,adap_dt_tol,adap_dt_max_iters]')
 
     integer :: bubble_model !< Gilmore or Keller--Miksis bubble model
     integer :: thermal      !< Thermal behavior. 1 = adiabatic, 2 = isotherm, 3 = transfer
@@ -436,10 +436,9 @@ module m_global_parameters
     integer, parameter :: nmom = 6 !< Number of carried moments per R0 location
     integer :: nmomsp    !< Number of moments required by ensemble-averaging
     integer :: nmomtot   !< Total number of carried moments moments/transport equations
-    integer :: R0_type
 
     real(wp) :: pi_fac   !< Factor for artificial pi_inf
-    $:GPU_DECLARE(create='[qbmm, nmomsp,nmomtot,R0_type,pi_fac]')
+    $:GPU_DECLARE(create='[qbmm, nmomsp,nmomtot,pi_fac]')
 
     #:if not MFC_CASE_OPTIMIZATION
         $:GPU_DECLARE(create='[nb]')
@@ -689,11 +688,10 @@ contains
             num_fluids = dflt_int
         #:endif
 
-        R0_type = dflt_int
-
         adv_n = .false.
         adap_dt = .false.
         adap_dt_tol = dflt_real
+        adap_dt_max_iters = dflt_adap_dt_max_iters
 
         pi_fac = 1._wp
 
@@ -924,7 +922,7 @@ contains
                         sys_size = n_idx
                     end if
 
-                    @:ALLOCATE(weight(nb), R0(nb), V0(nb))
+                    @:ALLOCATE(weight(nb), R0(nb))
                     @:ALLOCATE(bub_idx%rs(nb), bub_idx%vs(nb))
                     @:ALLOCATE(bub_idx%ps(nb), bub_idx%ms(nb))
 
@@ -965,11 +963,7 @@ contains
                     if (nb == 1) then
                         weight(:) = 1._wp
                         R0(:) = 1._wp
-                        V0(:) = 1._wp
-                    else if (nb > 1) then
-                        V0(:) = 1._wp
-                        !R0 and weight initialized in s_simpson
-                    else
+                    else if (nb < 1) then
                         stop 'Invalid value of nb'
                     end if
 
@@ -1041,7 +1035,7 @@ contains
 
                     @:ALLOCATE(bub_idx%rs(nb), bub_idx%vs(nb))
                     @:ALLOCATE(bub_idx%ps(nb), bub_idx%ms(nb))
-                    @:ALLOCATE(weight(nb), R0(nb), V0(nb))
+                    @:ALLOCATE(weight(nb), R0(nb))
 
                     do i = 1, nb
                         if (polytropic) then
@@ -1061,10 +1055,7 @@ contains
                     if (nb == 1) then
                         weight(:) = 1._wp
                         R0(:) = 1._wp
-                        V0(:) = 0._wp
-                    else if (nb > 1) then
-                        V0(:) = 1._wp
-                    else
+                    else if (nb < 1) then
                         stop 'Invalid value of nb'
                     end if
 
@@ -1276,7 +1267,7 @@ contains
         $:GPU_UPDATE(device='[momxb,momxe,advxb,advxe,contxb,contxe, &
             & bubxb,bubxe,intxb,intxe,sys_size,buff_size,E_idx, &
             & alf_idx,n_idx,adv_n,adap_dt,pi_fac,strxb,strxe, &
-            & chemxb,chemxe,c_idx]')
+            & chemxb,chemxe,c_idx,adap_dt_tol,adap_dt_max_iters]')
         $:GPU_UPDATE(device='[b_size,xibeg,xiend,tensor_size]')
 
         $:GPU_UPDATE(device='[species_idx]')
@@ -1303,7 +1294,7 @@ contains
         #:endif
 
         $:GPU_ENTER_DATA(copyin='[nb,R0ref,Ca,Web,Re_inv,weight,R0, &
-            & V0,bubbles_euler,polytropic,polydisperse,qbmm,R0_type, &
+            & bubbles_euler,polytropic,polydisperse,qbmm, &
             & ptil,bubble_model,thermal,poly_sigma]')
         $:GPU_ENTER_DATA(copyin='[R_n,R_v,phi_vn,phi_nv,Pe_c,Tw,pv, &
             & M_n,M_v,k_n,k_v,pb0,mass_n0,mass_v0,Pe_T, &
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
@@ -98,8 +98,9 @@ contains
             & 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', &
             & 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end',  'fd_order',     &
             & 'num_probes', 'num_integrals', 'bubble_model', 'thermal',        &
-            & 'R0_type', 'num_source', 'relax_model', 'num_ibs', 'n_start',    &
-            & 'num_bc_patches', 'num_igr_iters', 'num_igr_warm_start_iters' ]
+            & 'num_source', 'relax_model', 'num_ibs', 'n_start',    &
+            & 'num_bc_patches', 'num_igr_iters', 'num_igr_warm_start_iters', &
+            & 'adap_dt_max_iters' ]
             call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py
diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py
