From ee81564d48fd1dd56bfa035e5eef76a749db73f1 Mon Sep 17 00:00:00 2001
From: Thomas Roy <roy27@llnl.gov>
Date: Tue, 14 May 2024 17:04:56 -0700
Subject: [PATCH] Add equations for fluid solvers in docs

---
 docs/source/user_guide/fluid_solver.rst | 73 ++++++++++++++++++++-----
 echemfem/flow_solver.py                 |  2 +
 2 files changed, 62 insertions(+), 13 deletions(-)

diff --git a/docs/source/user_guide/fluid_solver.rst b/docs/source/user_guide/fluid_solver.rst
index 84205d3..cd0ec53 100644
--- a/docs/source/user_guide/fluid_solver.rst
+++ b/docs/source/user_guide/fluid_solver.rst
@@ -9,13 +9,33 @@ Navier-Stokes Solver
 --------------------
 
 The :class:`echemfem.NavierStokesFlowSolver` class contains an implementation
-for the incompressible Navier-Stokes equations. Both nondimensional and
-dimensional equations are available. Physical parameters are passed as a
-:py:class:`dict` in the ``fluid_params`` argument.
+for the steady-state incompressible Navier-Stokes equations. Both
+nondimensional and dimensional equations are available.
 
-To use the nondimensional version, the user must pass the following key:
+We are solving for velocity :math:`\mathbf u` and pressure :math:`p`. The
+dimensional version of the momentum equation is
 
-* ``"Reynolds number"``
+.. math::
+
+   -\nu \nabla^2 \mathbf u + \mathbf u \cdot \nabla \mathbf u + \frac{1}{\rho} \nabla p = 0,
+
+where :math:`\nu` is the kinematic viscosity and :math:`\rho` is the density.
+It is also common to use the dynamic viscosity :math:`\mu = \nu \rho`
+
+The nondimensional version is
+
+.. math::
+
+   - \frac{1}{\mathrm{Re}}\nabla^2 \mathbf u + \mathbf u \cdot \nabla \mathbf u + \nabla p = 0,
+
+where all quantities have been nondimensionalized, and :math:`\mathrm{Re}` is the Reynolds number. In both cases, we also have the incompressibility condition
+
+.. math::
+
+   \nabla \cdot \mathbf{u} = 0.
+
+Physical parameters are passed as a
+:py:class:`dict` in the ``fluid_params`` argument.
 
 For the dimensional version, the user must pass the following keys:
 
@@ -23,31 +43,58 @@ For the dimensional version, the user must pass the following keys:
 
 * ``"dynamic viscosity"`` or ``"kinematic viscosity"``
 
+To use the nondimensional version, the user must pass the following key:
+
+* ``"Reynolds number"``
+
 Navier-Stokes-Brinkman Solver
 -----------------------------
 
 The :class:`echemfem.NavierStokesBrinkmanFlowSolver` class contains an
-implementation for the incompressible Navier-Stokes-Brinkman equations. Both
-nondimensional and dimensional equations are available. In addition to the
-parameters provided for Navier-Stokes, the physical parameters below must be
-provided.
+implementation for the steady-state incompressible Navier-Stokes-Brinkman
+equations. Both nondimensional and dimensional equations are available.
 
-To use the nondimensional version, the user must also pass the following key:
+The dimensional version of the momentum equation is
 
-* ``"Darcy number"``
+.. math::
+
+   -\nabla\left(\nu_\mathrm{eff} \mathbf u \right)+ \mathbf u \cdot \nabla \mathbf u + \frac{1}{\rho} \nabla p + \nu K^{-1} \mathbf u = 0,
+
+where :math:`\nu_\mathrm{eff}` is the effective viscosity in the porous medium,
+and :math:`K`, its permeability. The inverse permeability can be provided
+directly in cases where it is zero in some regions, i.e. liquid-only regions.
+It is common to take :math:`\nu_\mathrm{eff}=\nu` for simplicity (the default
+here).
+
+The nondimensional implementation currently assume :math:`\nu_\mathrm{eff}=\nu`
+and :math:`K>0`, such that
+
+
+.. math::
+
+   - \frac{1}{\mathrm{Re}}\nabla^2 \mathbf u + \mathbf u \cdot \nabla \mathbf u + \nabla p + \frac{1}{\mathrm{Re}\mathrm{Da}} \mathbf u = 0,
+
+where all quantities have been nondimensionalized and :math:`\mathrm{Da}` is the Darcy number.
+
+In addition to the parameters provided for Navier-Stokes, the physical
+parameters below must be provided.
 
 For the dimensional version, the user must also pass the following keys:
 
 * ``"permeability"`` or ``"inverse permeability"``
 
-* Optional: ``"effective kinematic viscosity"``
+* Optional: ``"effective kinematic viscosity"`` or ``"effective dynamic viscosity"``
+
+To use the nondimensional version, the user must also pass the following key:
+
+* ``"Darcy number"``
 
 Boundary conditions
 -------------------
 
 A :py:class:`dict` containing boundary markers is passed when creating a
 ``FlowSolver`` object. Below are the different options for ``boundary_markers``
-keys. The velocity is denoted as :math:`\mathbf u` and pressure as :math:`p`.
+keys.
 
 * ``"no slip"``: Sets velocity to zero, :math:`\mathbf u = 0`.
 
diff --git a/echemfem/flow_solver.py b/echemfem/flow_solver.py
index f62e51b..24ccd78 100644
--- a/echemfem/flow_solver.py
+++ b/echemfem/flow_solver.py
@@ -271,6 +271,8 @@ def setup_problem(self):
                 inv_K = params["inverse permeability"]
             if params.get("effective kinematic viscosity"):
                 nu_eff = params.get("effective kinematic viscosity")
+            elif params.get("effective dynamic viscosity"):
+                nu_eff = params.get("effective dynamic viscosity") / rho
             else:
                 nu_eff = nu
             F = nu_eff * inner(grad(u), grad(v)) * dx \