Merge pull request idaholab#15684 from WilkAndy/geochem_transport_156…

…68_15683

Transport kernels added to geochem module

framework/include/kernels/ConservativeAdvection.h

@@ -35,7 +35,7 @@ class ConservativeAdvection : public Kernel
   virtual void computeJacobian() override;
 
   /// advection velocity
-  RealVectorValue _velocity;
+  const VectorVariableValue & _velocity;
 
   /// enum to make the code clearer
   enum class JacRes
@@ -62,4 +62,3 @@ class ConservativeAdvection : public Kernel
   /// Calculates the fully-upwind Residual and Jacobian (depending on res_or_jac)
   void fullUpwind(JacRes res_or_jac);
 };
-

framework/src/kernels/ConservativeAdvection.C

@@ -20,7 +20,7 @@ ConservativeAdvection::validParams()
   InputParameters params = Kernel::validParams();
   params.addClassDescription("Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak "
                              "form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.");
-  params.addRequiredParam<RealVectorValue>("velocity", "Velocity vector");
+  params.addRequiredCoupledVar("velocity", "Velocity vector");
   MooseEnum upwinding_type("none full", "none");
   params.addParam<MooseEnum>("upwinding_type",
                              upwinding_type,
@@ -32,7 +32,7 @@ ConservativeAdvection::validParams()
 
 ConservativeAdvection::ConservativeAdvection(const InputParameters & parameters)
   : Kernel(parameters),
-    _velocity(getParam<RealVectorValue>("velocity")),
+    _velocity(coupledVectorValue("velocity")),
     _upwinding(getParam<MooseEnum>("upwinding_type").getEnum<UpwindingType>()),
     _u_nodal(_var.dofValues()),
     _upwind_node(0),
@@ -43,7 +43,7 @@ ConservativeAdvection::ConservativeAdvection(const InputParameters & parameters)
 Real
 ConservativeAdvection::negSpeedQp() const
 {
-  return -_grad_test[_i][_qp] * _velocity;
+  return -_grad_test[_i][_qp] * _velocity[_qp];
 }
 
 Real

modules/geochemistry/doc/content/source/kernels/GeochemistryDispersion.md

@@ -0,0 +1,19 @@
+# GeochemistryDispersion
+
+This Kernel implements the PDE fragment
+\begin{equation}
+\nabla(\phi D\nabla c) \ ,
+\end{equation}
+which is part of the [transport](transport.md) equations.  Here
+
+- $\phi$ is the porosity, which may be a fixed real number or may be an `AuxVariable` that is spatially-dependent
+- $D$ is the hydrodynamic dispersion tensor (called `tensor_coeff` in the input file)
+- $c$ is the concentration (mol/m$^{3}$(aqueous solution)) of an aqueous species
+- $\nabla$ denotes the vector of spatial derivatives
+
+!syntax parameters /Kernels/GeochemistryDispersion
+
+!syntax inputs /Kernels/GeochemistryDispersion
+
+!syntax children /Kernels/GeochemistryDispersion
+

modules/geochemistry/doc/content/source/kernels/GeochemistryTimeDerivative.md

@@ -0,0 +1,23 @@
+# GeochemistryTimeDerivative
+
+This Kernel implements the PDE fragment
+\begin{equation}
+\phi \frac{\partial c}{\partial t} \ ,
+\end{equation}
+which is part of the [transport](transport.md) equations.  Here
+
+- $\phi$ is the porosity, which may be a fixed real number or may be an `AuxVariable` that is spatially-dependent
+- $c$ is the concentration (mol/m$^{3}$(aqueous solution)) of an aqueous species
+- $t$ is time
+
+Two notable features of this Kernel are:
+
+- it should not be used for cases were porosity is time-dependent (in this case $\phi$ should be inside the time derivative);
+- mass lumping to the nodes is used in order to reduce the possibility that concentration becomes negative.
+
+!syntax parameters /Kernels/GeochemistryTimeDerivative
+
+!syntax inputs /Kernels/GeochemistryTimeDerivative
+
+!syntax children /Kernels/GeochemistryTimeDerivative
+

modules/geochemistry/include/kernels/GeochemistryDispersion.h

@@ -0,0 +1,32 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "AnisotropicDiffusion.h"
+
+/**
+ * Kernel describing grad(porosity * dispersion * grad(concentration)), where porosity is an
+ * AuxVariable, dispersion is the hydrodynamic dispersion tensor (input by user as tensor_coeff),
+ * and concentration is the variable for this Kernel
+ */
+class GeochemistryDispersion : public AnisotropicDiffusion
+{
+public:
+  static InputParameters validParams();
+
+  GeochemistryDispersion(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+
+private:
+  const VariableValue & _porosity;
+};

modules/geochemistry/include/kernels/GeochemistryTimeDerivative.h

@@ -0,0 +1,34 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "TimeKernel.h"
+
+/**
+ * Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable.
+ * This Kernel should not be used if porosity is time-dependent.
+ * Mass lumping is employed for numerical stability
+ */
+class GeochemistryTimeDerivative : public TimeKernel
+{
+public:
+  static InputParameters validParams();
+
+  GeochemistryTimeDerivative(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+
+private:
+  const VariableValue & _nodal_u_dot;
+  const VariableValue & _nodal_du_dot_du;
+  const VariableValue & _porosity;
+};

modules/geochemistry/src/kernels/GeochemistryDispersion.C

@@ -0,0 +1,41 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "GeochemistryDispersion.h"
+
+registerMooseObject("GeochemistryApp", GeochemistryDispersion);
+
+InputParameters
+GeochemistryDispersion::validParams()
+{
+  InputParameters params = AnisotropicDiffusion::validParams();
+  params.addCoupledVar("porosity", 1.0, "Porosity");
+  params.addClassDescription("Kernel describing grad(porosity * tensor_coeff * "
+                             "grad(concentration)), where porosity is an AuxVariable (or just "
+                             "a real number), tensor_coeff is the hydrodynamic dispersion tensor "
+                             "and concentration is the 'variable' for this Kernel");
+  return params;
+}
+
+GeochemistryDispersion::GeochemistryDispersion(const InputParameters & parameters)
+  : AnisotropicDiffusion(parameters), _porosity(coupledValue("porosity"))
+{
+}
+
+Real
+GeochemistryDispersion::computeQpResidual()
+{
+  return _porosity[_qp] * AnisotropicDiffusion::computeQpResidual();
+}
+
+Real
+GeochemistryDispersion::computeQpJacobian()
+{
+  return _porosity[_qp] * AnisotropicDiffusion::computeQpJacobian();
+}

modules/geochemistry/src/kernels/GeochemistryTimeDerivative.C

@@ -0,0 +1,47 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "GeochemistryTimeDerivative.h"
+
+registerMooseObject("GeochemistryApp", GeochemistryTimeDerivative);
+
+InputParameters
+GeochemistryTimeDerivative::validParams()
+{
+  InputParameters params = TimeKernel::validParams();
+  params.addCoupledVar("porosity", 1.0, "Porosity");
+  params.addClassDescription(
+      "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just "
+      "a real number) and concentration is the 'variable' for this Kernel.  This Kernel should not "
+      "be used if porosity is time-dependent.  Mass lumping is employed for numerical stability");
+  return params;
+}
+
+GeochemistryTimeDerivative::GeochemistryTimeDerivative(const InputParameters & parameters)
+  : TimeKernel(parameters),
+    _nodal_u_dot(_var.dofValuesDot()),
+    _nodal_du_dot_du(_var.dofValuesDuDotDu()),
+    _porosity(coupledValue("porosity"))
+{
+}
+
+Real
+GeochemistryTimeDerivative::computeQpResidual()
+{
+  return _test[_i][_qp] * _porosity[_qp] * _nodal_u_dot[_i];
+}
+
+Real
+GeochemistryTimeDerivative::computeQpJacobian()
+{
+  if (_i == _j)
+    return _test[_i][_qp] * _porosity[_qp] * _nodal_du_dot_du[_j];
+  else
+    return 0.0;
+}

modules/geochemistry/test/tests/kernels/advection_1.i

@@ -0,0 +1,87 @@
+# A step-like initial concentration is advected to the right using a constant velocity.
+# Because of the Dirichlet BC on the left, the step-like concentration profile is maintained (up to the usual numerical diffusion)
+# Because upwinding_type=full in the ConservativeAdvection Kernel, there are no overshoots and undershoots
+# The total amount of "conc" should increase by dt * velocity every timestep, as recorded by the front_position Postprocessor
+[Mesh]
+  type = GeneratedMesh
+  dim = 1
+  nx = 100
+[]
+
+[Variables]
+  [./conc]
+  [../]
+[]
+
+[ICs]
+  [./conc]
+    type = FunctionIC
+    function = 'if(x<=0.25, 1, 0)'
+    variable = conc
+  [../]
+[]
+
+[BCs]
+  [./left]
+    type = DirichletBC
+    boundary = left
+    value = 1.0
+    variable = conc
+  [../]
+[]
+
+[Kernels]
+  [./dot]
+    type = GeochemistryTimeDerivative
+    variable = conc
+  [../]
+  [./adv]
+    type = ConservativeAdvection
+    velocity = velocity
+    upwinding_type = full
+    variable = conc
+  [../]
+[]
+
+[AuxVariables]
+  [./velocity]
+    family = MONOMIAL_VEC
+    order = CONSTANT
+  [../]
+[]
+
+[AuxKernels]
+  [./velocity]
+    type = VectorFunctionAux
+    function = vel_fcn
+    variable = velocity
+  [../]
+[]
+
+[Functions]
+  [./vel_fcn]
+    type = ParsedVectorFunction
+    value_x = 1
+    value_y = 0
+    value_z = 0
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = Newton
+  dt = 0.01
+  end_time = 0.1
+[]
+
+[Postprocessors]
+  [./front_position]
+    type = ElementIntegralVariablePostprocessor
+    variable = conc
+  [../]
+[]
+
+[Outputs]
+  csv = true
+[]
+