Enhancements for Euler-Lagrange bubble/particle model

docs/documentation/case.md

@@ -482,9 +482,9 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
 
 - `mp_weno` activates monotonicity preservation in the WENO reconstruction (MPWENO) such that the values of reconstructed variables do not reside outside the range spanned by WENO stencil ([Balsara and Shu, 2000](references.md); [Suresh and Huynh, 1997](references.md)).
 
-- `muscl_order` specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using `muscl_order = 2`, `muscl_lim` must be defined. 
+- `muscl_order` specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using `muscl_order = 2`, `muscl_lim` must be defined.
 
-- `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5. 
+- `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5.
 `muscl_lim = 1`, `2`, `3`, `4`, and `5` correspond to minmod, monotonized central, Van Albada, Van Leer, and SUPERBEE, respectively.
 
 - `int_comp` activates interface compression using THINC used in MUSCL Reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
@@ -599,6 +599,24 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running.
 | `output_partial_domain` | Logical | Output part of the domain |
 | `[x,y,z]_output%beg` | Real    | Beginning of the output domain in the [x,y,z]-direction |
 | `[x,y,z]_output%end` | Real    | End of the output domain in the [x,y,z]-direction |
+| `lag_txt_wrt`        | Logical | Write Lagrangian bubble data to `.dat` files |
+| `lag_header`         | Logical | Write header to Lagrangian bubble `.dat` files |
+| `lag_db_wrt`         | Logical | Write Lagrangian bubble data to silo/hdf5 database files |
+| `lag_id_wrt`         | Logical | Add the global bubble idea to the database file |
+| `lag_pos_wrt`        | Logical | Add the bubble position to the database file |
+| `lag_pos_prev_wrt`   | Logical | Add the previous bubble position to the database file |
+| `lag_vel_wrt`        | Logical | Add the bubble translational velocity to the database file |
+| `lag_rad_wrt`        | Logical | Add the bubble radius to the database file |
+| `lag_rvel_wrt`       | Logical | Add the bubble radial velocity to the database file |
+| `lag_r0_wrt`         | Logical | Add the bubble initial radius to the database file |
+| `lag_rmax_wrt`       | Logical | Add the bubble maximum radius to the database file |
+| `lag_rmin_wrt`       | Logical | Add the bubble minimum radius to the database file |
+| `lag_dphidt_wrt`     | Logical | Add the bubble subgrid velocity potential to the database file |
+| `lag_pres_wrt`       | Logical | Add the bubble pressure to the database file |
+| `lag_mv_wrt`         | Logical | Add the bubble vapor mass to the database file |
+| `lag_mg_wrt`         | Logical | Add the bubble gas mass to the database file |
+| `lag_betaT_wrt`      | Logical | Add the bubble heat flux model coefficient to the database file |
+| `lag_betaC_wrt`      | Logical | Add the bubble mass flux model coefficient to the database file |
 
 The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing.
 
@@ -628,7 +646,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 
 - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%beg` and `[x,y,z]_output%end`.
 This is useful for large domains where only a portion of the domain is of interest.
-It is not supported when `precision = 1` and `format = 1`. 
+It is not supported when `precision = 1` and `format = 1`.
 It also cannot be enabled with `flux_wrt`, `heat_ratio_wrt`, `pres_inf_wrt`, `c_wrt`, `omega_wrt`, `ib`, `schlieren_wrt`, `qm_wrt`, or 'liutex_wrt'.
 
 ### 8. Acoustic Source {#acoustic-source}

examples/2D_lagrange_bubblescreen/case.py

@@ -111,6 +111,7 @@
             "precision": 2,
             "prim_vars_wrt": "T",
             "parallel_io": "T",
+            "lag_db_wrt": "T",
             # Patch 1: Water (left)
             "patch_icpp(1)%geometry": 3,
             "patch_icpp(1)%x_centroid": 0.0,

examples/3D_lagrange_bubblescreen/case.py

@@ -120,6 +120,7 @@
             "precision": 2,
             "prim_vars_wrt": "T",
             "parallel_io": "T",
+            "lag_db_wrt": "T",
             # Patch 1: Water (left)
             "patch_icpp(1)%geometry": 9,
             "patch_icpp(1)%x_centroid": 0.0,

src/common/include/3dHardcodedIC.fpp

@@ -7,6 +7,67 @@
     ! Case 302 - IGR Jet
     real(wp) :: r, ux_th, ux_am, p_th, p_am, rho_th, rho_am, y_th, z_th, r_th, eps_smooth
 
+    real(wp), allocatable, dimension(:, :) :: ih
+    integer :: pos, start, end
+    logical :: file_exist
+    character(len=10000) :: line
+    character(len=25) :: value
+
+    if (patch_icpp(patch_id)%hcid == 303) then
+        allocate(ih(0:n_glb, 0:p_glb))
+
+        if (interface_file == '.') then
+            call s_mpi_abort("Error: interface_file must be specified for hcid=303")
+        else
+            inquire (file=trim(interface_file), exist=file_exist)
+            if (file_exist) then
+                open(unit=10, file=trim(interface_file), status="old", action="read")
+                do i = 0, m_glb
+                    read(10, '(A)') line  ! Read a full line as a string
+                    start = 1
+
+                    do j = 0, p_glb
+                        end = index(line(start:), ',')  ! Find the next comma
+                        if (end == 0) then
+                            value = trim(adjustl(line(start:)))  ! Last value in the line
+                        else
+                            value = trim(adjustl(line(start:start+end-2)))  ! Extract substring
+                            start = start + end  ! Move to next value
+                        end if
+                        read(value, *) ih(i, j)  ! Convert string to numeric value
+                        if (.not. f_is_default(normMag)) ih(i, j) = ih(i, j) * normMag
+                        if (.not. f_is_default(normFac)) ih(i, j) = ih(i, j) + normFac
+                    end do
+                end do
+                close(10)
+            else
+                call s_mpi_abort("Error: interface_file specified for hcid=303 does not exist")
+            end if
+        end if
+    end if
+
+    if (patch_icpp(patch_id)%hcid == 304) then
+        allocate(ih(0:n_glb, 0:0))
+        if (interface_file == '.') then
+            call s_mpi_abort("Error: interface_file must be specified for hcid=304")
+        else
+            inquire (file=trim(interface_file), exist=file_exist)
+            if (file_exist) then
+                open(unit=10, file=trim(interface_file), status="old", action="read")
+                do i = 0, n_glb
+                    read(10, '(A)') line  ! Read a full line as a string
+                    value = trim(line)
+                    read(value, *) ih(i, 0)  ! Convert string to numeric value
+                    if (.not. f_is_default(normMag)) ih(i, 0) = ih(i, 0) * normMag
+                    if (.not. f_is_default(normFac)) ih(i, 0) = ih(i, 0) + normFac
+                end do
+                close(10)
+            else
+                call s_mpi_abort("Error: interface_file specified for hcid=304 does not exist")
+            end if
+        end if
+    end if
+
     eps = 1e-9_wp
 
 #:enddef
@@ -86,6 +147,43 @@
 
         q_prim_vf(E_idx)%sf(i, j, k) = p_th*f_cut_on(r - r_th, eps_smooth)*f_cut_on(x_cc(i), eps_smooth) + p_am
 
+    case (303) ! 3D Interface from file cartesian
+
+        alph = 0.5_wp * (1 + (1._wp - 2._wp * eps) * &
+                    tanh((ih(start_idx(2) + j,start_idx(3) + k)  - x_cc(i))*(0.5_wp / dx)))
+
+        q_prim_vf(advxb)%sf(i,j,k) = 1._wp - alph
+        q_prim_vf(advxe)%sf(i,j,k) = alph
+
+        q_prim_vf(contxb)%sf(i,j,k) = q_prim_vf(advxb)%sf(i,j,k) * (950._wp / 1000._wp)
+        q_prim_vf(contxe)%sf(i,j,k) = q_prim_vf(advxe)%sf(i,j,k) * (1._wp / 1000._wp)
+
+        h = x_cc(i) - ih(start_idx(2) + j, start_idx(3) + k)
+        q_prim_vf(momxb)%sf(i,j,k) = -1._wp * (ih(start_idx(2) + j, start_idx(3) + k) - normFac) * exp(-h * h / 1000) / 100._wp
+
+        q_prim_vf(E_idx)%sf(i,j,k) = p0 + &
+            (q_prim_vf(contxb)%sf(i,j,k) + q_prim_vf(contxe)%sf(i,j,k)) * g0 * &
+            (ih(start_idx(2) + j, start_idx(3) + k) - x_cc(i))
+
+        if (surface_tension) q_prim_vf(c_idx)%sf(i,j,k) = alph
+
+    case (304) ! 3D Interface from file axisymmetric
+
+        alph = 0.5_wp * (1 + (1._wp - 2._wp * eps) * &
+                    tanh((ih(start_idx(2) + j,0)  - x_cc(i))*(0.01_wp / dx)))
+
+        q_prim_vf(advxb)%sf(i,j,k) = alph
+        q_prim_vf(advxe)%sf(i,j,k) = 1._wp - alph
+
+        q_prim_vf(contxb)%sf(i,j,k) = q_prim_vf(advxb)%sf(i,j,k) * 1._wp
+        q_prim_vf(contxe)%sf(i,j,k) = q_prim_vf(advxe)%sf(i,j,k) * (1._wp / 950._wp)
+
+        q_prim_vf(E_idx)%sf(i,j,k) = p0 + &
+            (q_prim_vf(contxb)%sf(i,j,k) + q_prim_vf(contxe)%sf(i,j,k)) * g0 * &
+            (ih(start_idx(1) + i, start_idx(3) + k) - y_cc(j))
+
+        if (surface_tension) q_prim_vf(c_idx)%sf(i,j,k) = alph
+
     case (370)
         ! This hardcoded case extrudes a 2D profile to initialize a 3D simulation domain
         @: HardcodedReadValues()

src/common/m_boundary_common.fpp

@@ -1942,6 +1942,21 @@ contains
         offset_x%beg = buff_size; offset_x%end = buff_size
         offset_y%beg = buff_size; offset_y%end = buff_size
         offset_z%beg = buff_size; offset_z%end = buff_size
+
+#ifdef MFC_MPI
+        ! Populate global domain boundaries with stretched grids
+        call s_mpi_allreduce_min(x_cb(-1), glb_bounds(1)%beg)
+        call s_mpi_allreduce_max(x_cb(m), glb_bounds(1)%end)
+
+        if (n > 0) then
+            call s_mpi_allreduce_min(y_cb(-1), glb_bounds(2)%beg)
+            call s_mpi_allreduce_max(y_cb(n), glb_bounds(2)%end)
+            if (p > 0) then
+                call s_mpi_allreduce_min(z_cb(-1), glb_bounds(3)%beg)
+                call s_mpi_allreduce_max(z_cb(p), glb_bounds(3)%end)
+            end if
+        end if
+#endif
 #endif
 
 #ifndef MFC_PRE_PROCESS

src/common/m_constants.fpp

@@ -13,6 +13,7 @@ module m_constants
     real(wp), parameter :: small_alf = 1.e-11_wp                !< Small alf tolerance
     real(wp), parameter :: pi = 3.141592653589793_wp !< Pi
     real(wp), parameter :: verysmall = 1.e-12_wp              !< Very small number
+    real(wp), parameter :: Re_b_min = 1.e-6_wp                !< minimum bubble reynolds number for drag coeff calc.
 
     integer, parameter :: num_stcls_min = 5                       !< Minimum # of stencils
     integer, parameter :: path_len = 400                          !< Maximum path length
@@ -21,7 +22,7 @@ module m_constants
     integer, parameter :: fourier_rings = 5                       !< Fourier filter ring limit
     integer, parameter :: num_fluids_max = 10                     !< Maximum number of fluids in the simulation
     integer, parameter :: num_probes_max = 10                     !< Maximum number of flow probes in the simulation
-    integer, parameter :: num_patches_max = 10
+    integer, parameter :: num_patches_max = 20
     integer, parameter :: num_bc_patches_max = 10
     integer, parameter :: pathlen_max = 400
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
@@ -61,6 +62,7 @@ module m_constants
     ! Lagrange bubbles 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
+    integer, parameter :: lag_io_vars = 21 ! Number of variables per particle for MPI_IO
 
     ! Strang Splitting constants
     real(wp), parameter :: dflt_adap_dt_tol = 1.e-4_wp !< Default tolerance for adaptive step size

src/common/m_derived_types.fpp

@@ -447,6 +447,12 @@ module m_derived_types
         logical :: write_bubbles            !< Write files to track the bubble evolution each time step
         logical :: write_bubbles_stats      !< Write the maximum and minimum radius of each bubble
         integer :: nBubs_glb                !< Global number of bubbles
+        integer :: vel_model                !< Particle velocity model
+        integer :: drag_model               !< Particle drag model
+        logical :: pressure_force           !< Include pressure force translational motion
+        logical :: gravity_force            !< Include gravity force in translational motion
+        logical :: momentum_transfer_force  !< Include momentum transfer from radial dynamics in translational motion
+        real(wp) :: c_d                     !< Drag coefficient
         real(wp) :: epsilonb         !< Standard deviation scaling for the gaussian function
         real(wp) :: charwidth        !< Domain virtual depth (z direction, for 2D simulations)
         real(wp) :: valmaxvoid       !< Maximum void fraction permitted

src/common/m_helper_basic.fpp

@@ -112,10 +112,11 @@ contains
 
     pure subroutine s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, &
                                                   igr_order, buff_size, idwint, idwbuff, &
-                                                  viscous, bubbles_lagrange, m, n, p, num_dims, igr, ib)
+                                                  viscous, bubbles_lagrange, m, n, p, num_dims, &
+                                                  igr, ib, fd_number)
 
         integer, intent(in) :: recon_type, weno_polyn, muscl_polyn
-        integer, intent(in) :: m, n, p, num_dims, igr_order
+        integer, intent(in) :: m, n, p, num_dims, igr_order, fd_number
         integer, intent(inout) :: buff_size
         type(int_bounds_info), dimension(3), intent(inout) :: idwint, idwbuff
         logical, intent(in) :: viscous, bubbles_lagrange
@@ -140,7 +141,7 @@ contains
 
         ! Correction for smearing function in the lagrangian subgrid bubble model
         if (bubbles_lagrange) then
-            buff_size = max(buff_size, 6)
+            buff_size = max(buff_size + fd_number, 6 + fd_number)
         end if
 
         if (ib) then