Doco and tests of ConservativeAdvection

Fixes idaholab#11852

framework/doc/content/bib/moose.bib

@@ -173,3 +173,25 @@ @article{Nguyen2014
 volume = {74},
 year = {2014}
 }
+
+@article{dalen1979,
+abstract = {This paper summarizes some research that was conducted to construct finite-element models for reservoir flow problems. The models are based on Galerkin's method, but the method is applied in an unorthodox manner to simplify calculation of coefficients and to improve stability. Specifically, techniques of compatibility relaxation, capacity lumping, and upstream mobility weighting are used, and schemes are obtained that seem to combine the simplicity and high stability of conventional finite-difference models with the generality and modeling flexibility of finite-element methods. The development of a model for single-phase gas flow and a two-phase oil/water model is described. Numerical examples are included to illustrate the usefulness of finite elements. In particular, the triangular element with linear interpolation is shown to be an attractive alternative to the standard five-point, finite-difference approximation for two-dimensional analysis.},
+author = {Dalen,~V},
+doi = {10.2118/7196-PA},
+journal = {Society of Petroleum Engineers Journal},
+pages = {1-11},
+title = {{Simplified Finite-Element Models for Reservoir Flow Problems}},
+volume = {19},
+year = {1979}
+}
+
+@article{adhikary2011,
+abstract = {We demonstrate that care must be taken in a finite element discretisation of multi-phase compressible Darcy flow, otherwise constraints of non-negativity of fluid mass can become violated. Generalising a technique pioneered by Dalen, a lumped, finite-difference-inspired approximation with upstream weighting is described. This is numerically cheap, physically elegant, and the fluid mass at all nodes remains non-negative.},
+author = {Adhikary~D. and Wilkins~A.},
+doi = {10.4028/www.scientific.net/DDF.318.33},
+journal = {Defect and Diffusion Forum},
+pages = {33--40},
+title = {{Analysis of finite element discretisation schemes for multi-phase Darcy flow}},
+volume = {318},
+year = {2011}
+}

framework/doc/content/source/kernels/ConservativeAdvection.md

@@ -2,7 +2,7 @@
 
 ## Description
 
-The `ConservativeAdvection` kernel implements an advection term given for the domain ($\Omega$) define as
+The `ConservativeAdvection` kernel implements an advection term given for the domain ($\Omega$) defined as
 
 \begin{equation}
 \underbrace{\nabla \cdot \vec{v} u}_{\textrm{ConservativeAdvection}} + \sum_{i=1}^n \beta_i = 0 \in \Omega,
@@ -11,7 +11,7 @@ where $v$ is the advecting velocity and the second term on the left hand side
 represents the strong forms of other kernels. `ConservativeAdvection` does not assume
 that the velocity is divergence free and instead applies $\nabla$ to the test
 function $\psi_i$ in the weak variational form after integrating by parts,
-as in
+which results in the following (without numerical stabilization)
 
 \begin{equation}
 R_i(u_h) = \underbrace{-(\nabla \psi_i, \vec{v} u)}_{\textrm{ConservativeAdvection}} + \langle\psi_i, \vec{v} u
@@ -22,24 +22,96 @@ element solution of the weak formulation. The first term is the volumetric term
 is a surface term describing the advective flux out of the
 volume. `ConservativeAdvection` corresponds to the former volumetric term.
 
-The corresponding Jacobian is given by
+Without numerical stabilization the corresponding Jacobian is given by
 
 \begin{equation}
 \frac{\partial R_i(u_h)}{\partial u_j} = -(\nabla \psi_i, \vec{v} \phi_j).
 \end{equation}
 
+## Full upwinding
+
+Advective flow is notoriously prone to physically-incorrect overshoots
+and undershoots, so in many simulations some numerical stabilization
+is used to reduce or eliminate this spurious behaviour.
+Full-upwinding [cite:dalen1979,adhikary2011] is an example of
+numerical stabilization and this essentially adds numerical diffusion
+to completely eliminate overshoots and undershoots.  Full-upwinding is
+available in `ConservativeAdvection` by setting the `upwinding_type`
+appropriately.
+
+For each element in the mesh, full upwinding is implemented in the
+following way.  For each $i$, the quantity
+\begin{equation}
+\tilde{R}_{i} = -(\nabla\psi_i, \vec{v}) \ ,
+\end{equation}
+is evaluated.  If $\tilde{R}_{i}>0$ is positive then mass (or heat, or whatever $u$
+represents) is flowing *out* of node $i$, and this is called an
+"upwind" node.  If $\tilde{R}_{i}\geq 0$, the residual for node $i$
+is set to
+\begin{equation}
+R_{i} = u_{i}\tilde{R}_{i} \ \ \texttt{ for } \ \ \tilde{R}_{i}\geq 0.
+\end{equation}
+Here $u_{i}$ is the value of $u$ at the node $i$.  Define the total
+mass flowing from the upwind nodes:
+\begin{equation}
+M_{\mathrm{out}} = \sum_{i \mathrm{ with } \tilde{R}_{i}\geq 0}
+u_{i}\tilde{R}_{i} \ .
+\end{equation}
+Similarly, define
+\begin{equation}
+\tilde{M}_{\mathrm{in}} = - \sum_{i \mathrm{ with } \tilde{R}_{i}< 0}
+\tilde{R}_{i} \ .
+\end{equation}
+Mass is conserved if the residual for the downwind nodes is defined to
+be
+\begin{equation}
+R_{i} = \tilde{R}_{i} M_{\mathrm{out}} / M_{\mathrm{in}}  \ \ \texttt{ for } \ \ \tilde{R}_{i}< 0.
+\end{equation}
+
+## Dirichlet boundary conditions for pure advection
+
+In the purely advective case, Dirichlet boundary conditions may be
+used to inject mass (or heat, or whatever $u$ corresponds to
+physically) into the domain, or extract it from the domain.
+Specifically:
+
+- if $\vec{v}\cdot\vec{n}<0$ (with $\vec{n}$ being the outward normal) then the velocity field is "blowing" mass into the domain.  Fixing $u=u_{B}$ on this boundary means the mass flux (units kg.s$^{-1}$.m$^{-2}$) entering the domain is $-\vec{v}\cdot\vec{n}u_{B}$.
+- if $\vec{v}\cdot\vec{n}<0$ (with $\vec{n}$ being the outward normal) then the velocity field is "blowing" mass out of the domain.  If no Dirichlet boundary condition is used then mass is prevented from leaving the domain via this boundary: mass arriving here will "build up", causing $u$ to increase.  Alternatively, if $u=0$ is fixed at the boundary all mass arriving here will exit the domain immediately.
+
+
 ## Example Syntax
 
 `ConservativeAdvection` can be used in a variety of problems, including
 advection-diffusion-reaction. The syntax for `ConservativeAdvection` is
-demonstrated in this `Kernels` block from an advection-diffusion-reaction test
-case:
+demonstrated in this `Kernels` block from a pure advection test case:
 
 !listing
-test/tests/dgkernels/advection_diffusion_mixed_bcs_test_resid_jac/dg_advection_diffusion_test.i
+test/tests/kernels/conservative_advection/no_upwinding_1D.i
 block=Kernels label=false
 
-The velocity is supplied as a three component vector with order $v_x$, $v_y$, and  $v_z$.
+The velocity is supplied as a three component vector with order $v_x$,
+$v_y$, and  $v_z$.
+
+## Comparison of no upwinding and full upwinding
+
+The tests
+[no_upwinding_1D](test/tests/kernels/conservative_advection/no_upwinding_1D.i)
+and
+[full_upwinding_1D](test/tests/kernels/conservative_advection/full_upwinding_1D.i)
+describe the same physical situation: advection from left to right
+along a line.  A source at the left boundary introduces $u$ into the
+domain (via a Dirichlet boundary condition).  The right boundary is
+impermeable (no Dirichlet boundary condition), so when $u$ arrives
+there is starts building up at the boundary.  It is clear from the
+figures below that no upwinding leads to unphysical overshoots and
+undershoots, while full upwinding results in none of that oscillatory
+behaviour but instead produces more numerical diffusion.
+
+!media media/framework/kernels/conservative_advection_1d_1.png
+ caption=$u$ after 0.1 seconds of advection
+
+!media media/framework/kernels/conservative_advection_1d_5.png
+ caption=$u$ after 0.5 seconds of advection
 
 !syntax parameters /Kernels/ConservativeAdvection
 

test/tests/kernels/conservative_advection/full_upwinding_1D.i

@@ -0,0 +1,48 @@
+# ConservativeAdvection with upwinding_type = full
+# Apply a velocity = (1, 0, 0) and see a pulse advect to the right
+# Note that the pulse diffuses more than with no upwinding,
+# but there are no overshoots and undershoots and that the
+# center of the pulse at u=0.5 advects with the correct velocity
+[Mesh]
+  type = GeneratedMesh
+  dim = 1
+  nx = 10
+[]
+
+[Variables]
+  [./u]
+  [../]
+[]
+
+[BCs]
+  [./u_fixed_left]
+    type = DirichletBC
+    boundary = left
+    variable = u
+    value = 1
+  [../]
+[]
+
+[Kernels]
+  [./udot]
+    type = MassLumpedTimeDerivative
+    variable = u
+  [../]
+  [./advection]
+    type = ConservativeAdvection
+    variable = u
+    velocity = '1 0 0'
+    upwinding_type = full
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = LINEAR
+  dt = 0.1
+  end_time = 1
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/kernels/conservative_advection/full_upwinding_2D.i

@@ -0,0 +1,47 @@
+# 2D test of advection with full upwinding
+# Note there are no overshoots or undershoots
+# but there is numerical diffusion.
+# The center of the blob advects with the correct velocity
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 40
+  ny = 40
+[]
+
+[Variables]
+  [./u]
+  [../]
+[]
+
+[ICs]
+  [./u_blob]
+    type = FunctionIC
+    variable = u
+    function = 'if(x<0.2,if(y<0.2,1,0),0)'
+  [../]
+[]
+
+[Kernels]
+  [./udot]
+    type = MassLumpedTimeDerivative
+    variable = u
+  [../]
+  [./advection]
+    type = ConservativeAdvection
+    variable = u
+    upwinding_type = full
+    velocity = '2 1 0'
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = LINEAR
+  dt = 0.01
+  end_time = 0.1
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/kernels/conservative_advection/full_upwinding_jacobian.i

@@ -0,0 +1,53 @@
+# Test of advection with full upwinding
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  nx = 3
+  ny = 2
+  nz = 1
+[]
+
+[Variables]
+  [./u]
+  [../]
+[]
+
+[ICs]
+  [./u]
+    type = RandomIC
+    variable = u
+  [../]
+[]
+
+[BCs]
+  [./u_fixed_left]
+    type = DirichletBC
+    boundary = left
+    variable = u
+    value = 1
+  [../]
+[]
+
+[Kernels]
+  [./advection]
+    type = ConservativeAdvection
+    variable = u
+    upwinding_type = full
+    velocity = '2 -1.1 1.23'
+  [../]
+[]
+
+[Preconditioning]
+  [./andy]
+    type = SMP
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = NEWTON
+  petsc_options_iname = '-snes_type'
+  petsc_options_value = 'test'
+  dt = 2
+  end_time = 2
+[]

test/tests/kernels/conservative_advection/gold/none_in_none_out_out.csv

@@ -0,0 +1,12 @@
+time,total_mass
+0,0
+1,25
+2,25
+3,25
+4,25
+5,25
+6,25
+7,25
+8,25
+9,25
+10,25

test/tests/kernels/conservative_advection/no_upwinding_1D.i

@@ -0,0 +1,45 @@
+# ConservativeAdvection with upwinding_type = None
+# Apply a velocity = (1, 0, 0) and see a pulse advect to the right
+# Note there are overshoots and undershoots
+[Mesh]
+  type = GeneratedMesh
+  dim = 1
+  nx = 10
+[]
+
+[Variables]
+  [./u]
+  [../]
+[]
+
+[BCs]
+  [./u_fixed_left]
+    type = DirichletBC
+    boundary = left
+    variable = u
+    value = 1
+  [../]
+[]
+
+[Kernels]
+  [./udot]
+    type = TimeDerivative
+    variable = u
+  [../]
+  [./advection]
+    type = ConservativeAdvection
+    variable = u
+    velocity = '1 0 0'
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = LINEAR
+  dt = 0.1
+  end_time = 1
+[]
+
+[Outputs]
+  exodus = true
+[]