fix some WARNING, add bib ref

.github/workflows/test_IAMReX.yml

@@ -4,6 +4,9 @@ on:
   # schedule:
   #   - cron: "*/5 * * * *"
   push:
+    paths:
+      - 'Source/**'
+
 jobs:
   test-IAMReX:
     runs-on: ubuntu-latest
@@ -31,4 +34,4 @@ jobs:
           python Test_IAMReX/test_IAMReX.py
       - name: Finished.
         run: |
-          echo "Tests completed!"
+          echo "Tests completed!"

Docs/IAMReX_documentation/source/IAMRandIAMReX.rst

@@ -7,7 +7,7 @@ IAMR (Base)
 - Standard incompressible Navier-Stokes
 
 IAMReX (Extended)
-----------------
+-----------------
 - ``NUM_STATE_MAX = AMREX_SPACEDIM+5`` (adds level set field)
 - Additional advance methods:
   - ``advance_semistaggered_twophase_ls()``: Two-phase flow with level sets
@@ -72,4 +72,4 @@ Understanding the Flow
 ----------------------
 1. Variable setup happens in ``NS_setup.cpp``
 2. Time advancement occurs in ``NavierStokes::advance()`` methods
-3. Boundary conditions are applied via functions in NS_BC.H
+3. Boundary conditions are applied via functions in NS_BC.H

Docs/IAMReX_documentation/source/LevelSet.rst

@@ -1,7 +1,7 @@
 .. _LevelSetMethod:
 
 Level Set Method
-===============
+================
 
 The level set method is a powerful tool for capturing the interface between two immiscible fluids. It is based on the idea of a signed distance function, which is a function that assigns a signed distance to each point in the domain. The level set function is defined as follows:
 
@@ -12,7 +12,7 @@ The level set method is a powerful tool for capturing the interface between two
 where :math:`\phi` is the level set function, :math:`\mathbf{x}` is the position vector, and :math:`t` is the time.
 
 Material Properties
-------------------
+-------------------
 
 The density and viscosity are defined as functions of the level set field:
 
@@ -36,7 +36,7 @@ where :math:`H(\phi)` is the Heaviside function:
    \end{cases}
 
 Time Discretization
-------------------
+-------------------
 
 For a single level, the momentum equation :eq:`eq:ns` is advanced by a fractional step method with the approximate projection to enforce the incompressibility condition (equation :eq:`eq:div`). The LS advection equation :eq:`eq:phi` is updated using the Godunov scheme.
 
@@ -56,9 +56,9 @@ At the beginning of each time advancement of level :math:`l`, the velocity :math
    .. math::
       :label: eq:viscsolve
 
-      \begin{aligned} 
-      &\boldsymbol{u^{*,n+1}}-\frac{\Delta t}{2\rho(\phi^{n+1/2})Re}{\nabla} \cdot {\mu(\phi^{n+1})}{\nabla}\boldsymbol{u^{*,n+1}} = 
-      \boldsymbol{u^n}-\Delta t \left[\nabla \cdot (\boldsymbol{uu})\right]^{n+1/2} + \\ &\frac{\Delta t}{\rho(\phi^{n+1/2})}\bigg[-\nabla p^{n-1/2}+ 
+      \begin{aligned}
+      &\boldsymbol{u^{*,n+1}}-\frac{\Delta t}{2\rho(\phi^{n+1/2})Re}{\nabla} \cdot {\mu(\phi^{n+1})}{\nabla}\boldsymbol{u^{*,n+1}} =
+      \boldsymbol{u^n}-\Delta t \left[\nabla \cdot (\boldsymbol{uu})\right]^{n+1/2} + \\ &\frac{\Delta t}{\rho(\phi^{n+1/2})}\bigg[-\nabla p^{n-1/2}+
       \frac{1}{2Re}{\nabla}\cdot{\mu(\phi^{n})}{\nabla}\boldsymbol{u^n} + \rho(\phi^{n+1/2}) \frac{z}{Fr^2} -  \frac{1}{We}\kappa(\phi^{n+1/2})\delta(x^{n+1/2})\boldsymbol{n}\bigg].
       \end{aligned}
 
@@ -74,25 +74,25 @@ At the beginning of each time advancement of level :math:`l`, the velocity :math
 3. **Apply the projection method** to obtain the pressure and a solenoidal velocity field. To conduct the level projection, a temporary variable :math:`\boldsymbol{V}` is defined as
 
    .. math::
-      :label: eq:ns_lp1
+      :label: eq:ns_lp_ls1
 
       \boldsymbol{V} =  \frac{\boldsymbol{{u^{*,n+1}}}}{\Delta t} + \frac{1}{\rho(\phi^{n+1/2})} \nabla p^{n-1/2}
 
    Then the updated pressure :math:`p^{n+1/2}` is calculated by
 
    .. math::
-      :label: eq:ns_lp2
+      :label: eq:ns_lp_ls2
 
       L^{cc,\mathrm{level}}_{\rho^{n+1/2}} p^{n+1/2} =  \nabla \cdot \boldsymbol{V}
 
    where :math:`L^{cc,\mathrm{level}}_{\rho^{n+1/2}}p^{n+1/2}` is a density-weighted approximation to :math:`\nabla \cdot (1/\rho^{n+1/2} \nabla p^{n+1/2})`. Finally, the velocity can be calculated as
 
    .. math::
-      :label: eq:ns_lp3
+      :label: eq:ns_lp_ls3
 
       \boldsymbol{{u^{n+1}}} = \Delta t \left(\boldsymbol{V} - \frac{1}{\rho^{n+1/2}} \nabla p^{n+1/2}\right)
 
-   As defined in the AMReX framework, :math:`\nabla \cdot` and :math:`\nabla` are the cell-centered level divergence operator :math:`D^{cc,\mathrm{level}}` and level gradient operator :math:`G^{cc,\mathrm{level}}`, respectively. The level gradient operator :math:`G^{cc,\mathrm{level}}` is not the minus transpose of the level divergence operator :math:`D^{cc,\mathrm{level}}`, i.e., :math:`G^{cc,\mathrm{level}} \neq -(D^{cc,\mathrm{level}})^T`. As a result, the idempotency of the approximate projection :math:`\boldsymbol{P} = I - G^{cc,\mathrm{level}}(L^{cc,\mathrm{level}})^{-1}D^{cc,\mathrm{level}}` is not ensured, i.e., :math:`\boldsymbol{P}^{2} \neq \boldsymbol{P}`. Yet, this nonidempotent approximate projection is stable and appears to be well-behaved in various numerical tests and practical applications. Notably, for a uniform single grid with periodic boundary conditions, Lai theoretically proved that this approximate projection method is stable, in that :math:`\|\boldsymbol{P}\| \leq 1`. It should be noted that the approximate projection is applied to the intermediate velocity :math:`\boldsymbol{{u^{*,n+1}}}` (equation :eq:`eq:ns_lp1`). Compared with the form that projects the increment velocity :math:`\boldsymbol{u^{*,n+1}}-\boldsymbol{u^n}`, e.g. as that used in other methods, the projection method used here can reduce the accumulation of pressure errors and lead to a more stable algorithm. The effectiveness and stability of this approximate projection has been validated through various numerical tests.
+   As defined in the AMReX framework, :math:`\nabla \cdot` and :math:`\nabla` are the cell-centered level divergence operator :math:`D^{cc,\mathrm{level}}` and level gradient operator :math:`G^{cc,\mathrm{level}}`, respectively. The level gradient operator :math:`G^{cc,\mathrm{level}}` is not the minus transpose of the level divergence operator :math:`D^{cc,\mathrm{level}}`, i.e., :math:`G^{cc,\mathrm{level}} \neq -(D^{cc,\mathrm{level}})^T`. As a result, the idempotency of the approximate projection :math:`\boldsymbol{P} = I - G^{cc,\mathrm{level}}(L^{cc,\mathrm{level}})^{-1}D^{cc,\mathrm{level}}` is not ensured, i.e., :math:`\boldsymbol{P}^{2} \neq \boldsymbol{P}`. Yet, this nonidempotent approximate projection is stable and appears to be well-behaved in various numerical tests and practical applications. Notably, for a uniform single grid with periodic boundary conditions, Lai theoretically proved that this approximate projection method is stable, in that :math:`\|\boldsymbol{P}\| \leq 1`. It should be noted that the approximate projection is applied to the intermediate velocity :math:`\boldsymbol{{u^{*,n+1}}}` (equation :eq:`eq:ns_lp_ls1`). Compared with the form that projects the increment velocity :math:`\boldsymbol{u^{*,n+1}}-\boldsymbol{u^n}`, e.g. as that used in other methods, the projection method used here can reduce the accumulation of pressure errors and lead to a more stable algorithm. The effectiveness and stability of this approximate projection has been validated through various numerical tests.
 
 4. **Reinitialize the LS function** :math:`\phi` to maintain :math:`\phi` as a signed distance function of the interface and guarantee the conservation of the mass of the two phases. In this step, a temporary LS function :math:`d(\boldsymbol{x},\tau)` is updated iteratively using the following pseudo evolution equation:
 
@@ -121,6 +121,6 @@ At last, we give a summary of the single-level advancement algorithm as follows.
 
 1. Advance the LS function using equation :eq:`eq:s0phin1`
 2. Solve the intermediate velocity using equation :eq:`eq:viscsolve`
-3. Apply the projection method to update the pressure and velocity field following equations :eq:`eq:ns_lp1`--:eq:`eq:ns_lp3`
+3. Apply the projection method to update the pressure and velocity field following equations :eq:`eq:ns_lp_ls1`--:eq:`eq:ns_lp_ls3`
 4. Re-initialize the LS function on the single level using equations :eq:`eq:ns_reinit1`--:eq:`eq:ns_reinit3`
 

Docs/IAMReX_documentation/source/Results.rst

@@ -29,7 +29,7 @@ where :math:`H` is the Heaviside function, defined by:
 
 .. math::
 
-   H(\phi) = 
+   H(\phi) =
    \begin{cases}
    0, & \phi \le 0\\
    1, & \phi > 0
@@ -101,7 +101,7 @@ for any Eulerian cell :math:`(i,j,k)`. The computational domain is :math:`L_x \t
    | 128   | 0.40209828964  | 6.359·10⁻³          | 2.0678      |
    +-------+----------------+---------------------+-------------+
 
-The numerical errors decrease with the increase of the :math:`d/h`, where :math:`h` is the Cartesian grid spacing on level :math:`0`. If the resolution on the finest level keeps unchanged, we validated that the results of a three-level grid are the same as those of the corresponding single-level grid. In addition, our results show the second-order convergence and agree well with the results in kempe. It also matches the overall second-order accuracy of the basic fluid solver.
+The numerical errors decrease with the increase of the :math:`d/h`, where :math:`h` is the Cartesian grid spacing on level :math:`0`. If the resolution on the finest level keeps unchanged, we validated that the results of a three-level grid are the same as those of the corresponding single-level grid. In addition, our results show the second-order convergence and agree well with the results in Kempe :cite:`kempe2012improved`. It also matches the overall second-order accuracy of the basic fluid solver.
 
 Lastly, it is noted that this method is also applicable when multiple particles are close to each other or their surfaces are in direct contact. Because the PVF calculation is a separate operation for each particle, the total volume fraction is not needed as long as the Eulerian force considers the effects of all particles .
 
@@ -114,7 +114,7 @@ Flow Past Sphere
 
     The schematic of the flow passing through the spherical particles
 
-We validate the accuracy and efficacy of our adaptive solver by simulating a spherical particle in uniform flow with different particle Reynolds numbers. The diameter of the particle is :math:`D_p = 1`, the computational domain is :math:`L_x \times L_y \times L_z = 20D_p \times 10D_p \times 10D_p`, the distance of the particle from the inlet is :math:`d = 5D_p` and located in the center of the yz plane. The inlet and outlet boundaries are applied in the x direction and the inlet velocity :math:`U` is :math:`1m/s`. Both y and z directions are periodic boundaries. 
+We validate the accuracy and efficacy of our adaptive solver by simulating a spherical particle in uniform flow with different particle Reynolds numbers :cite:`schiller1933uber,zhu2022particle`. The diameter of the particle is :math:`D_p = 1`, the computational domain is :math:`L_x \times L_y \times L_z = 20D_p \times 10D_p \times 10D_p`, the distance of the particle from the inlet is :math:`d = 5D_p` and located in the center of the yz plane. The inlet and outlet boundaries are applied in the x direction and the inlet velocity :math:`U` is :math:`1m/s`. Both y and z directions are periodic boundaries.
 
 The influence of AMR on the simulation results is investigated by using the subcycling method with different levels. As shown in Fig.
 
@@ -134,7 +134,7 @@ The theoretical S-N law for calculating the drag coefficient of the shaped parti
 .. math::
     C_D = (24/Re_p)(1+0.15Re_p^{0.687}),
 
-which is proposed by schiller et al., and :math:`Re_p = UD_p/\nu` represents the particle Reynolds number. It can be seen from above Fig that the present results under different particle Reynolds numbers are in good agreement with S-N law. The fact that different levels of grid produce the nearly identical results validated the accuracy of our solver on the adaptive grid.
+which is proposed by Schiller :cite:`schiller1933uber`., and :math:`Re_p = UD_p/\nu` represents the particle Reynolds number. It can be seen from above Fig that the present results under different particle Reynolds numbers are in good agreement with S-N law. The fact that different levels of grid produce the nearly identical results validated the accuracy of our solver on the adaptive grid.
 
 
 Cluster of monodisperse particles
@@ -148,7 +148,7 @@ We demonstrate the accuracy and efficacy of our codes for simulating clusters of
 
     Monodisperse particles on a three-level AMR grid
 
-80 particles of diameter :math:`D = 1` are randomly distributed in a channel of size :math:`L_x\times L_y \times L_z = 10\times 20 \times 10`. To choose an optimal interaction number :math:`N_s` in this complex configuration, the maximum error of the no-slip boundary condition among 80 particles is tested with a unit flow field :math:`u=(1,0,0)`.the maximum error of no-slip condition decreases as :math:`N_s` increases and it is strongly reduced for :math:`N_s=2`. According to the selection suggestions provided by Breugem et al., :math:`N_s=2` is the optimal value for balancing the accuracy of the no-slip boundary and the computational efficiency. After determining :math:`N_s`, the fluid flow is driven by applying a pressure gradient of 1.0 in the z direction. This case can represent a porous medium with a volume fraction of 0.02. Three levels of the AMR grid is applied. The grid resolution on the finest level is :math:`d/h=16`. Since the multi-direct forcing immersed boundary method and fictitious domain method require cube grid cells, the grid cell requirement is equals to case 1 in Table.
+80 particles of diameter :math:`D = 1` are randomly distributed in a channel of size :math:`L_x\times L_y \times L_z = 10\times 20 \times 10`. To choose an optimal interaction number :math:`N_s` in this complex configuration, the maximum error of the no-slip boundary condition among 80 particles is tested with a unit flow field :math:`u=(1,0,0)`.the maximum error of no-slip condition decreases as :math:`N_s` increases and it is strongly reduced for :math:`N_s=2`. According to the selection suggestions provided by Breugem :cite:`breugem2012second`. :math:`N_s=2` is the optimal value for balancing the accuracy of the no-slip boundary and the computational efficiency. After determining :math:`N_s`, the fluid flow is driven by applying a pressure gradient of 1.0 in the z direction. This case can represent a porous medium with a volume fraction of 0.02. Three levels of the AMR grid is applied. The grid resolution on the finest level is :math:`d/h=16`. Since the multi-direct forcing immersed boundary method :cite:`kidanemariam2022open,yousefi2023role` and fictitious domain method :cite:`xia2020effects,fan2023three` require cube grid cells, the grid cell requirement is equals to case 1 in Table.
 
     +------+---------+---------+---------+-------------+
     | case | level 0 | level 1 | level 2 | Total cells |
@@ -158,14 +158,14 @@ We demonstrate the accuracy and efficacy of our codes for simulating clusters of
     |   2  |  128000 |  534656 | 1593664 |   2256320   |
     +------+---------+---------+---------+-------------+
 
-Compared with them, our algorithm has a 72.5% grid reduction with :math:`d/h=16`. And it has a 62.5% Lagrangian markers reduction compared with the DLM method with :math:`d/h=16`.
+Compared with them, our algorithm has a 72.5% grid reduction with :math:`d/h=16`. And it has a 62.5% Lagrangian markers reduction compared with the DLM method :cite:`sharma2022coupled,zeng2022subcycling` with :math:`d/h=16`.
 
 When the simulation reaches the steady state, the total pressure drop balances the IB force generated by all particles in the streamwise z direction. Following the equation in akiki et al., the theoretical drag force is given by
 
 .. math::
     F_{theory} = (\frac{\Delta p}{\Delta z}L_z)L_xL_y
 
-Fig. shows the time series of total IB force for all particles. The resistance gradually reaches a steady state after 40000 steps. In this case, the theoretical value of drag force given by :math:`F_{theory}` is 2000, while the present average values at steady state with :math:`N_s=2` and 4 are all converged around 2000. It indicates that :math:`N_s=2` is sufficient for this case. And the agreement between theory and present results validates the accuracy of our proposed framework in dealing with large amounts of particles in the fluid system. 
+Fig. shows the time series of total IB force for all particles. The resistance gradually reaches a steady state after 40000 steps. In this case, the theoretical value of drag force given by :math:`F_{theory}` is 2000, while the present average values at steady state with :math:`N_s=2` and 4 are all converged around 2000. It indicates that :math:`N_s=2` is sufficient for this case. And the agreement between theory and present results validates the accuracy of our proposed framework in dealing with large amounts of particles in the fluid system.
 
 .. figure:: ./Results/MonodispersPDrop.png
     :align: center

Docs/IAMReX_documentation/source/Software_Chapter.rst

@@ -49,8 +49,8 @@ Key Computational Functions
 - ``NavierStokesBase::rk_second_reinit()``: Second RK step during reinitialization
 - ``NavierStokesBase::mass_fix()``: Mass fix subroutine during reinitialization
 
-Features
-^^^^^^^^
+Features(**Level Set**)
+^^^^^^^^^^^^^^^^^^^^^^^
 - **Sharp Interface**: Maintains interface thickness of 1.5~2 grid cells
 - **Mass Conservation**: Conservative advection schemes
 - **Reinitialization**: Periodic distance function correction
@@ -82,8 +82,8 @@ Particle Data Structure (**Kernel**)
         Vector<Real> phiK, thetaK;              // Spherical marker distribution
     }
 
-Key Computational functions
-^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Key Computational functions(**IBM**)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 1. Fluid-Solid Coupling
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -109,14 +109,14 @@ Key Computational functions
 4. Particle Dynamics
 ~~~~~~~~~~~~~~~~~~~~
 - ``mParticle::UpdateParticles()``: 6-DOF motion integration
-- **Constraint Handling**: 
+- **Constraint Handling**:
   - Translation locks (``TL[i]``): 0=fixed, 1=prescribed, 2=free
   - Rotation locks (``RL[i]``): Similar constraint system
 - **Collision Integration**: Seamless coupling with collision forces
 - ``nodal_phi_to_pvf()``: Particle volume fraction calculation
 
-Features
-^^^^^^^^
+Features(**IBM**)
+^^^^^^^^^^^^^^^^^
 
 Marker Distribution
 ~~~~~~~~~~~~~~~~~~~
@@ -135,10 +135,10 @@ Core Architecture of Particle Collision
 
 Primary Classes (**Collision**)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-- **``ParticleCollision``**: Main collision management
-- **``CollisionParticle``**: Particle representation for collisions
-- **``CollisionPair``**: Collision pair structure
-- **``CollisionCell``**: Spatial hashing cell
+- ``ParticleCollision``: Main collision management
+- ``CollisionParticle``: Particle representation for collisions
+- ``CollisionPair``: Collision pair structure
+- ``CollisionCell``: Spatial hashing cell
 
 Collision Detection Algorithm
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -196,4 +196,4 @@ V. Summary
 
 IAMeX represents a substantial advancement in computational multi-physics, transforming the single-phase IAMR code into a comprehensive platform for complex fluid-solid interaction simulations. The integration of diffused immersed boundary methods, particle collision dynamics, and multiple interface tracking approaches provides researchers with unprecedented capabilities for studying real-world multi-physics phenomena while maintaining computational efficiency and scalability.
 
-The modular design and robust implementation ensure that IAMeX serves as both a production simulation tool and a research platform for developing next-generation multi-physics algorithms. Its contributions to the computational fluid dynamics community extend beyond mere feature additions, representing fundamental advances in the numerical treatment of complex interfacial and particulate flows.
+The modular design and robust implementation ensure that IAMeX serves as both a production simulation tool and a research platform for developing next-generation multi-physics algorithms. Its contributions to the computational fluid dynamics community extend beyond mere feature additions, representing fundamental advances in the numerical treatment of complex interfacial and particulate flows.