trixi-framework
diff --git a/Diff for: ‎.typos.toml
+1 b/Diff for: ‎.typos.toml
+1
diff --git a/Diff for: ‎NEWS.md
+23 b/Diff for: ‎NEWS.md
+23
diff --git a/Diff for: ‎Project.toml
+1-1 b/Diff for: ‎Project.toml
+1-1
diff --git a/Diff for: ‎docs/src/refs.bib
+48 b/Diff for: ‎docs/src/refs.bib
+48
diff --git a/Diff for: ‎docs/src/systems/fluid.md
+180-32 b/Diff for: ‎docs/src/systems/fluid.md
+180-32
diff --git a/Diff for: ‎examples/fluid/dam_break_oil_film_2d.jl
+10-2 b/Diff for: ‎examples/fluid/dam_break_oil_film_2d.jl
+10-2
@@ -9,3 +9,4 @@ Strack = "Strack"
 Lok = "Lok"
 PN = "PN"
 Typ = "Typ"
+Yiu = "Yiu"
@@ -3,6 +3,29 @@
 TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.2.7
+
+### Features
+
+- Adds the classic **Continuum Surface Force (CSF)** model based on Morris 2000 (#584), which computes
+  surface tension as a **body force** proportional to curvature and directed along the interface normal.
+  This method is efficient and accurate for capillary effects but does not explicitly conserve momentum.
+
+- Added the classic **Continuum Surface Stress (CSS)** model based on Morris 2000 (#584), which is
+  a **momentum-conserving** approach that formulates surface tension as the **divergence of a stress tensor**.
+  However, it requires additional computation and stabilization to handle **high-density interfaces** and reduce numerical instabilities.
+
+- Implement `BoundaryZone` to allow for bidirectional flow (#623)
+
+
+### Testing
+
+- Run CI tests on GPUs via Buildkite CI (#723)
+
+### Bugs
+
+- Fix GPU computations (#689)
+
 ## Version 0.2.6
 
 ### Features
 
@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.2.7-dev"
+version = "0.2.7"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 
@@ -491,6 +491,19 @@ @Article{Morris1997
   publisher = {Elsevier BV},
 }
 
+@article{Morris2000,
+  author = {Morris, Joseph P.},
+  title = {Simulating surface tension with smoothed particle hydrodynamics},
+  journal = {International Journal for Numerical Methods in Fluids},
+  volume = {33},
+  number = {3},
+  pages = {333-353},
+  keywords = {interfacial flow, meshless methods, surface tension},
+  doi = {https://doi.org/10.1002/1097-0363(20000615)33:3<333::AID-FLD11>3.0.CO;2-7},
+  year = {2000}
+}
+
+
 @InProceedings{Mueller2003,
   author    = {M{\"u}ller, Matthias and Charypar, David and Gross, Markus},
   booktitle = {Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation},
@@ -622,3 +635,38 @@ @Article{Wendland1995
 }
 
 @Comment{jabref-meta: databaseType:bibtex;}
+
+@book{Lange2005,
+  author = {{Lange's Handbook of Chemistry}},
+  title = {Lange's Handbook of Chemistry},
+  edition = {16th},
+  publisher = {McGraw-Hill},
+  year = {2005}
+}
+
+@article{MeloEspinosa2014,
+  author = {Melo-Espinosa, Eliezer A. and S{\'a}nchez-Borroto, Yisel and Errasti, Michel and Hansen, Alan C.},
+  title = {Surface Tension Prediction of Vegetable Oils Using Artificial Neural Networks and Multiple Linear Regression},
+  journal = {Energy Conversion and Management},
+  volume = {84},
+  pages = {50--60},
+  doi = {https://doi.org/10.1016/j.enconman.2014.08.004},
+  year = {2014}
+}
+
+@book{Hui1992,
+  author = {Hui, Yiu H.},
+  title = {Encyclopedia of Food Science and Technology},
+  year = {1992},
+  publisher = {Academic Press},
+  address = {London},
+  isbn = {9780122774308}
+}
+
+@book{Poling2001,
+  title = {The Properties of Gases and Liquids},
+  author = {Poling, Bruce E. and Prausnitz, John M. and O'Connell, John P.},
+  year = {2001},
+  publisher = {McGraw-Hill},
+  address = {New York}
+}
@@ -1,10 +1,22 @@
+
 # [Fluid Models](@id fluid_models)
-Currently available fluid methods are the [weakly compressible SPH method](@ref wcsph) and the [entropically damped artificial compressibility for SPH](@ref edac).  
+
+Currently available fluid methods are the [weakly compressible SPH method](@ref wcsph) and the
+[entropically damped artificial compressibility for SPH](@ref edac).  
 This page lists models and techniques that apply to both of these methods.  
 
 ## [Viscosity](@id viscosity_wcsph)
 
-TODO: Explain viscosity.
+Viscosity is a critical physical property governing momentum diffusion within a fluid.
+In the context of SPH, viscosity determines how rapidly velocity gradients are smoothed out,
+influencing key flow characteristics such as boundary layer formation, vorticity diffusion,
+and dissipation of kinetic energy. It also helps determine whether a flow is laminar or turbulent
+under a given set of conditions.
+
+Implementing viscosity correctly in SPH is essential for producing physically accurate results,
+and different methods exist to capture both numerical stabilization and true viscous effects.
+
+### API
 
 ```@autodocs
 Modules = [TrixiParticles]
@@ -18,71 +30,207 @@ Modules = [TrixiParticles]
 Pages = [joinpath("general", "corrections.jl")]
 ```
 
-
+---
 
 ## [Surface Normals](@id surface_normal)
+
+### Overview of surface normal calculation in SPH
+
+Surface normals are essential for modeling surface tension as they provide the directionality
+of forces acting at the fluid interface. They are calculated based on the particle properties and
+their spatial distribution.
+
+#### Color field and gradient-based surface normals
+
+The surface normal at a particle is derived from the color field, a scalar field assigned to particles
+to distinguish between different fluid phases or between fluid and air. The color field gradients point
+towards the interface, and the normalized gradient defines the surface normal direction.
+
+The simplest SPH formulation for a surface normal, ``n_a`` is given as
+
+```math
+n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab},
+```
+
+where:
+
+- ``c_b`` is the color field value for particle ``b``,
+- ``m_b`` is the mass of particle ``b``,
+- ``\rho_b`` is the density of particle ``b``,
+- ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel ``W_{ab}`` with respect to particle ``a``.
+
+#### Normalization of surface normals
+
+The calculated normals are normalized to unit vectors:
+
+```math
+\hat{n}_a = \frac{n_a}{\Vert n_a \Vert}.
+```
+
+Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces.
+
+#### Handling noise and errors in normal calculation
+
+In regions distant from the interface, the calculated normals may be small or inaccurate due to the
+smoothing kernel's support radius. To mitigate this:
+
+1. Normals below a threshold are excluded from further calculations.
+2. Curvature calculations use a corrected formulation to reduce errors near interface fringes.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "surface_normal_sph.jl")]
 ```
 
-## [Surface Tension](@id surface_tension)
+---
 
-### Akinci-based intra-particle force surface tension and wall adhesion model
-The work by Akinci proposes three forces:
-- a cohesion force
-- a surface area minimization force
-- a wall adhesion force
+## [Surface Tension](@id surface_tension)
 
-The classical model is composed of the curvature minimization and cohesion force.
+Surface tension is a key phenomenon in fluid dynamics, influencing the behavior of droplets, bubbles, and fluid interfaces.
+In SPH, surface tension is modeled as forces arising due to surface curvature and relative particle movement, ensuring realistic
+simulation of capillary effects, droplet coalescence, and fragmentation.
+
+The surface tension coefficient ``\sigma`` is a physical parameter that quantifies the energy required to increase the surface area
+of a fluid by a unit amount. A higher value of ``\sigma`` indicates that the fluid resists changes to its surface area more strongly,
+causing droplets or bubbles to assume shapes (often spherical) that minimize their surface. In practice, ``\sigma`` can be measured
+experimentally through techniques such as the pendant drop method, the Wilhelmy plate method, or the du Noüy ring method,
+each of which relates a measurable force or change in shape to the fluid’s surface tension. For pure substances,
+tabulated reference values of ``\sigma`` at given temperatures are commonly used, while for mixtures or complex fluids,
+direct experimental measurements or values can be estimated from empirical equation (see [Poling](@cite Poling2001) or [Lange](@cite Lange2005)).
+In the following table some values are shown for reference. The values marked with a '~' are complex mixtures that are estimated by an empirical equation (see [Poling](@cite Poling2001)).
+
+| **Fluid**    | **Surface Tension (``\sigma``) [N/m at 20°C]** |
+|--------------|----------------------------------------------:|
+| **Gasoline**    | ~0.022   [Poling](@cite Poling2001)             |
+| **Ethanol**     | 0.022386 [Lange](@cite Lange2005)               |
+| **Acetone**     | 0.02402  [Lange](@cite Lange2005)               |
+| **Mineral Oil** | ~0.030   [Poling](@cite Poling2001)             |
+| **Olive Oil**   | 0.03303  [Hui](@cite Hui1992), [MeloEspinosa](@cite MeloEspinosa2014) |
+| **Glycerol**    | 0.06314  [Lange](@cite Lange2005)               |
+| **Water**       | 0.07288  [Lange](@cite Lange2005)               |
+| **Mercury**     | 0.486502 [Lange](@cite Lange2005)               |
+
+### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf)
+
+The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components:
 
 #### Cohesion force
-The model calculates the cohesion force based on the support radius ``h_c`` and the distance between particles.
-This force is determined using two distinct regimes within the support radius:
-- For particles closer than half the support radius,
-  a repulsive force is calculated to prevent particles from clustering too tightly,
-  enhancing the simulation's stability and realism.
-- Beyond half the support radius and within the full support radius,
-  an attractive force is computed, simulating the effects of surface tension that draw particles together.
-The cohesion force, ``F_{\text{cohesion}}``, for a pair of particles is given by:
+
+The cohesion force captures the attraction between particles at the fluid interface, creating the effect of surface tension.
+It is defined by the distance between particles and the support radius ``h_c``, using a kernel-based formulation.
+
+**Key features:**
+
+- Particles within half the support radius experience a repulsive force to prevent clustering.
+- Particles beyond half the radius but within the support radius experience an attractive force to simulate cohesion.
+
+Mathematically:
+
 ```math
 F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert},
 ```
-where:
-- ``\sigma`` represents the surface tension coefficient, adjusting the overall strength of the cohesion effect.
-- ``C`` is a scalar function of the distance between particles.
 
-The cohesion kernel ``C`` is defined as
+where ``C(r)``, the cohesion kernel, is defined as:
+
 ```math
 C(r)=\frac{32}{\pi h_c^9}
 \begin{cases}
-(h_c-r)^3 r^3, & \text{if } 2r > h_c \\
-2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c \\
-0, & \text{otherwise}
+(h_c-r)^3 r^3, & \text{if } 2r > h_c, \\
+2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\
+0, & \text{otherwise.}
 \end{cases}
 ```
 
 #### Surface area minimization force
-To model the minimization of the surface area and curvature of the fluid, a curvature force is used, which is calculated as
+
+The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area.
+It acts based on the difference in surface normals:
+
 ```math
-F_{\text{curvature}} = -\sigma (n_a - n_b)
+F_{\text{curvature}} = -\sigma (n_a - n_b),
 ```
 
+where ``n_a`` and ``n_b`` are the surface normals of the interacting particles.
+
 #### Wall adhesion force
-The wall adhesion model proposed by Akinci et al. is based on a kernel function which is 0 from 0.0 to 0.5 support radiia with a maximum at 0.75.
-With the force calculated with an adhesion coefficient ``\beta`` as
+
+This force models the interaction between fluid and solid boundaries, simulating adhesion effects at walls.
+It uses a custom kernel with a peak at 0.75 times the support radius:
+
 ```math
 F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert},
 ```
-with ``A`` being the adhesion kernel defined as
+
+where ``A(r)`` is the adhesion kernel:
+
 ```math
-A(r)= \frac{0.007}{h_c^{3.25}}
+A(r) = \frac{0.007}{h_c^{3.25}}
 \begin{cases}
-\sqrt[4]{- \frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c \\
+\sqrt[4]{-\frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c, \\
 0, & \text{otherwise.}
 \end{cases}
 ```
 
+---
+
+### [Morris surface tension model](@id morris_csf)
+
+The method described by [Morris](@cite Morris2000) estimates curvature by combining particle color gradients (see [`surface_normal`](@ref)) and smoothing functions to derive surface normals.
+The computed curvature is then used to determine forces acting perpendicular to the interface.
+While this method provides accurate surface tension forces, it does not explicitly conserve momentum.
+
+In the Morris model, surface tension is computed based on local interface curvature ``\kappa`` and the unit surface normal ``\hat{n}.``
+By estimating ``\hat{n}`` and ``\kappa`` at each particle near the interface, the surface tension force for particle a can be written as:
+
+```math
+F_{\text{surface tension}} = - \sigma \frac{\kappa_a}{\rho_a}\hat{n}_a
+```
+
+This formulation focuses directly on geometric properties of the interface, making it relatively straightforward to implement when a reliable interface detection
+(e.g., a color function) is available. However, accurately estimating ``\kappa`` and ``n`` may require fine resolutions.
+
+---
+
+### [Morris-based momentum-conserving surface tension model](@id moriss_css)
+
+In addition to the simpler curvature-based formulation, [Morris](@cite Morris2000) introduced a momentum-conserving approach.
+This method treats surface tension forces as arising from the divergence of a stress tensor, ensuring exact conservation
+of linear momentum and offering more robust behavior for high-resolution or long-duration simulations
+where accumulated numerical error can be significant.
+
+#### Stress tensor formulation
+
+The surface tension force can be seen as a divergence of a stress tensor ``S``
+
+```math
+F_{\text{surface tension}} = \nabla \cdot S,
+```
+
+with ``S`` defined as
+
+```math
+S = \sigma \delta_s (I - \hat{n} \otimes \hat{n}),
+```
+
+with:
+
+- ``\delta_s``: Surface delta function,
+- ``\hat{n}``: Unit normal vector,
+- ``I``: Identity matrix.
+
+This divergence can be computed numerically in the SPH framework as
+
+```math
+\sum_b \frac{m_b}{\rho_a \rho_b} (S_a + S_b) \nabla W_{ab} 
+```
+
+#### Advantages and limitations
+
+While momentum conservation makes this model attractive, it requires additional computational effort and stabilization
+techniques to address instabilities in high-density regions.
+
+### API
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "surface_tension.jl")]
 
@@ -18,15 +18,16 @@ smoothing_length = 3.5 * fluid_particle_spacing
 sound_speed = 20 * sqrt(gravity * H)
 
 # Physical values
-nu_water = 8.9E-7
-nu_oil = 6E-5
+nu_water = 8.9e-7
+nu_oil = 6e-5
 nu_ratio = nu_water / nu_oil
 
 nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
 nu_sim_water = nu_ratio * nu_sim_oil
 
 oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 
+# TODO: broken if both systems use surface tension
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
               viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
@@ -60,6 +61,13 @@ oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
                                          correction=AkinciFreeSurfaceCorrection(oil_density),
                                          reference_particle_spacing=fluid_particle_spacing)
 
+# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
+#                                          oil_eos, smoothing_kernel,
+#                                          smoothing_length, viscosity=oil_viscosity,
+#                                          acceleration=(0.0, -gravity),
+#                                          surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03),
+#                                          reference_particle_spacing=fluid_particle_spacing)
+
 # ==========================================================================================
 # ==== Simulation
 semi = Semidiscretization(fluid_system, oil_system, boundary_system,