Merge pull request #24204 from lindsayad/vector-ins-convected-mesh

Some FSI additions

.coverage

@@ -36,7 +36,7 @@ require_total = 85
 require_total = 84
 
 [fsi]
-require_total = 90
+require_total = 85
 
 [functional_expansion_tools]
 require_total = 82

framework/include/interfacekernels/ADInterfaceKernel.h

@@ -35,28 +35,27 @@ class ADInterfaceKernelTempl : public InterfaceKernelBase, public NeighborMooseV
    * Using the passed DGResidual type, selects the correct test function space and residual block,
    * and then calls computeQpResidual
    */
-  void computeElemNeighResidual(Moose::DGResidualType type) override final;
+  void computeElemNeighResidual(Moose::DGResidualType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpJacobian
    */
-  void computeElemNeighJacobian(Moose::DGJacobianType type) override final;
+  void computeElemNeighJacobian(Moose::DGJacobianType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpOffDiagJacobian with the passed jvar
    */
-  void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
-                                       unsigned int jvar) override final;
+  void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar);
 
   /// Selects the correct Jacobian type and routine to call for the primary variable jacobian
-  void computeElementOffDiagJacobian(unsigned int jvar) override final;
+  virtual void computeElementOffDiagJacobian(unsigned int jvar) override final;
 
   /// Selects the correct Jacobian type and routine to call for the secondary variable jacobian
-  void computeNeighborOffDiagJacobian(unsigned int jvar) override final;
+  virtual void computeNeighborOffDiagJacobian(unsigned int jvar) override final;
 
   /// Computes the residual for the current side.
   void computeResidual() override final;

framework/include/interfacekernels/InterfaceKernel.h

@@ -48,22 +48,21 @@ class InterfaceKernelTempl : public InterfaceKernelBase, public NeighborMooseVar
    * Using the passed DGResidual type, selects the correct test function space and residual block,
    * and then calls computeQpResidual
    */
-  virtual void computeElemNeighResidual(Moose::DGResidualType type) override;
+  virtual void computeElemNeighResidual(Moose::DGResidualType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpJacobian
    */
-  virtual void computeElemNeighJacobian(Moose::DGJacobianType type) override;
+  virtual void computeElemNeighJacobian(Moose::DGJacobianType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpOffDiagJacobian with the passed jvar
    */
-  virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
-                                               unsigned int jvar) override;
+  virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar);
 
   /// Selects the correct Jacobian type and routine to call for the primary variable jacobian
   virtual void computeElementOffDiagJacobian(unsigned int jvar) override;

framework/include/interfacekernels/InterfaceKernelBase.h

@@ -35,26 +35,6 @@ class InterfaceKernelBase : public NeighborResidualObject,
   /// The neighbor variable number that this interface kernel operates on
   virtual const MooseVariableFEBase & neighborVariable() const = 0;
 
-  /**
-   * Using the passed DGResidual type, selects the correct test function space and residual block,
-   * and then calls computeQpResidual
-   */
-  virtual void computeElemNeighResidual(Moose::DGResidualType type) = 0;
-
-  /**
-   * Using the passed DGJacobian type, selects the correct test function and trial function spaces
-   * and
-   * jacobian block, and then calls computeQpJacobian
-   */
-  virtual void computeElemNeighJacobian(Moose::DGJacobianType type) = 0;
-
-  /**
-   * Using the passed DGJacobian type, selects the correct test function and trial function spaces
-   * and
-   * jacobian block, and then calls computeQpOffDiagJacobian with the passed jvar
-   */
-  virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar) = 0;
-
   /// Selects the correct Jacobian type and routine to call for the primary variable jacobian
   virtual void computeElementOffDiagJacobian(unsigned int jvar) = 0;
 

framework/src/kernels/ADKernelSUPG.C

@@ -21,15 +21,20 @@ ADKernelSUPGTempl<T>::validParams()
   InputParameters params = ADKernelStabilizedTempl<T>::validParams();
   params.addParam<MaterialPropertyName>(
       "tau_name", "tau", "The name of the stabilization parameter tau.");
-  params.addRequiredCoupledVar("velocity", "The velocity variable.");
+  params.addCoupledVar("velocity", "The velocity variable.");
+  params.addParam<MaterialPropertyName>("material_velocity",
+                                        "A material property describing the velocity");
   return params;
 }
 
 template <typename T>
 ADKernelSUPGTempl<T>::ADKernelSUPGTempl(const InputParameters & parameters)
   : ADKernelStabilizedTempl<T>(parameters),
     _tau(this->template getADMaterialProperty<Real>("tau_name")),
-    _velocity(this->adCoupledVectorValue("velocity"))
+    _velocity(
+        this->isParamValid("velocity")
+            ? this->adCoupledVectorValue("velocity")
+            : this->template getADMaterialProperty<RealVectorValue>("material_velocity").get())
 {
 }
 

modules/fsi/doc/content/source/interfacekernels/ADPenaltyVelocityContinuity.md

@@ -0,0 +1,15 @@
+# ADPenaltyVelocityContinuity
+
+!syntax description /InterfaceKernels/ADPenaltyVelocityContinuity
+
+This object is meant for use in coupling the Navier-Stokes fluid equations with
+solid mechanics equations at a fluid-structure interface. In that context this
+object imposes both continuity of stress and continuity of material velocity via
+a penalty condition. For more information on the error of penalty methods,
+please see [CoupledPenaltyInterfaceDiffusion.md].
+
+!syntax parameters /InterfaceKernels/ADPenaltyVelocityContinuity
+
+!syntax inputs /InterfaceKernels/ADPenaltyVelocityContinuity
+
+!syntax children /InterfaceKernels/ADPenaltyVelocityContinuity

modules/fsi/doc/content/source/kernels/ConvectedMesh.md

@@ -1,15 +1,22 @@
 # ConvectedMesh
 
-## Description
+!syntax description /Kernels/ConvectedMesh
 
 `ConvectedMesh` implements the corresponding weak form for the components of
 the term:
 
 \begin{equation}
+\label{convection}
 -\rho \left(\frac{\partial\vec{d_m}}{\partial t} \cdot \nabla\right) \vec{u}
 \end{equation}
 
 where $\rho$ is the density, $\vec{d_m}$ is the fluid mesh displacements, and
 $\vec{u}$ is the fluid velocity. This term is essential for obtaining the
 correct convective derivative of the fluid in cases where the fluid mesh is
 dynamic, e.g. in simulations of fluid-structure interaction.
+
+!syntax parameters /Kernels/ConvectedMesh
+
+!syntax inputs /Kernels/ConvectedMesh
+
+!syntax children /Kernels/ConvectedMesh

modules/fsi/doc/content/source/kernels/ConvectedMeshPSPG.md

@@ -0,0 +1,31 @@
+# ConvectedMeshPSPG
+
+!syntax description /Kernels/ConvectedMeshPSPG
+
+This object adds the pressure-stabilized Petrov-Galerkin term to the pressure
+equation corresponding to the [ConvectedMesh.md] object. It implements the
+weak form
+
+\begin{equation}
+\label{pspg}
+\int \nabla \psi \cdot \left(\tau\frac{\partial\vec{d_m}}{\partial t}\cdot\nabla\vec{u}\right) dV
+\end{equation}
+
+where $\nabla \psi$ is the gradient of the pressure test function, $\tau$ is a
+stabilization parameter computed automatically by the `navier_stokes` base class
+`INSBase`, $\vec{d_m}$ is the fluid mesh displacements, and
+$\vec{u}$ is the fluid velocity. Note that when comparing [pspg] with
+[!eqref](ConvectedMesh.md#convection) that the minus sign and density $\rho$ have
+disappeared. This is because the form of PSPG stabilization is
+
+\begin{equation}
+\int \nabla \psi \cdot \left(-\frac{\tau}{\rho}\vec{R}\right)
+\end{equation}
+
+where $\vec{R}$ denotes the strong form of the momentum residuals.
+
+!syntax parameters /Kernels/ConvectedMeshPSPG
+
+!syntax inputs /Kernels/ConvectedMeshPSPG
+
+!syntax children /Kernels/ConvectedMeshPSPG

modules/fsi/include/interfacekernels/ADPenaltyVelocityContinuity.h

@@ -0,0 +1,101 @@
+///* 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 "InterfaceKernelBase.h"
+#include "MooseVariableFE.h"
+
+/**
+ * Interface kernel for enforcing continuity of stress and velocity
+ */
+class ADPenaltyVelocityContinuity : public InterfaceKernelBase
+{
+public:
+  static InputParameters validParams();
+
+  ADPenaltyVelocityContinuity(const InputParameters & parameters);
+
+  virtual void computeResidual() override;
+  virtual void computeResidualAndJacobian() override;
+  virtual void computeJacobian() override;
+  virtual void computeElementOffDiagJacobian(unsigned int jvar) override;
+  virtual void computeNeighborOffDiagJacobian(unsigned int jvar) override;
+  virtual const MooseVariableFieldBase & variable() const override;
+  virtual const MooseVariableFieldBase & neighborVariable() const override;
+
+protected:
+  /// The penalty factor
+  const Real _penalty;
+
+  /// Fluid velocity variable
+  const VectorMooseVariable * const _velocity_var;
+
+  /// Fluid velocity values
+  const ADVectorVariableValue & _velocity;
+
+  /// Solid velocity values
+  std::vector<const ADVariableValue *> _solid_velocities;
+
+  /// Displacement variables
+  std::vector<const MooseVariable *> _displacements;
+
+  /// JxW with displacement derivatives
+  const MooseArray<ADReal> & _ad_JxW;
+
+  /// Coordinate transformation with displacement derivatives
+  const MooseArray<ADReal> & _ad_coord;
+
+  /// Residuals data member to avoid constant heap allocation
+  std::vector<ADReal> _residuals;
+
+  /// Jump data member to avoid constant heap allocations
+  std::vector<ADRealVectorValue> _qp_jumps;
+};
+
+inline void
+ADPenaltyVelocityContinuity::computeJacobian()
+{
+  computeResidual();
+}
+
+inline void
+ADPenaltyVelocityContinuity::computeResidualAndJacobian()
+{
+  computeResidual();
+}
+
+inline void
+ADPenaltyVelocityContinuity::computeElementOffDiagJacobian(const unsigned int jvar)
+{
+  if (jvar == _velocity_var->number())
+    // Only need to do this once because AD does everything all at once
+    computeResidual();
+}
+
+inline void
+ADPenaltyVelocityContinuity::computeNeighborOffDiagJacobian(unsigned int)
+{
+}
+
+inline const MooseVariableFieldBase &
+ADPenaltyVelocityContinuity::variable() const
+{
+  return *_velocity_var;
+}
+
+inline const MooseVariableFieldBase &
+ADPenaltyVelocityContinuity::neighborVariable() const
+{
+  if (_displacements.empty() || !_displacements.front())
+    mooseError("The 'neighborVariable' method was called which requires that displacements be "
+               "actual variables.");
+
+  return *_displacements.front();
+}

modules/fsi/include/kernels/ConvectedMesh.h

@@ -9,12 +9,12 @@
 
 #pragma once
 
-#include "Kernel.h"
+#include "INSBase.h"
 
 /**
  * This calculates the time derivative for a coupled variable
  **/
-class ConvectedMesh : public Kernel
+class ConvectedMesh : public INSBase
 {
 public:
   static InputParameters validParams();
@@ -26,6 +26,17 @@ class ConvectedMesh : public Kernel
   virtual Real computeQpJacobian() override;
   virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
 
+  /**
+   * Compute the Jacobian of the Petrov-Galerkin stabilization addition with respect to the provided
+   * velocity component
+   */
+  Real computePGVelocityJacobian(unsigned short component);
+
+  /**
+   * Compute the strong residual, e.g. -rho * ddisp/dt * grad(u)
+   */
+  Real strongResidual();
+
   const VariableValue & _disp_x_dot;
   const VariableValue & _d_disp_x_dot;
   const unsigned int _disp_x_id;
@@ -37,4 +48,10 @@ class ConvectedMesh : public Kernel
   const unsigned int _disp_z_id;
 
   const MaterialProperty<Real> & _rho;
+
+  /// Whether this kernel should include Streamline-Upwind Petrov-Galerkin stabilization
+  const bool _supg;
+
+  /// The velocity component this kernel is acting on
+  unsigned short _component;
 };

modules/fsi/include/kernels/ConvectedMeshPSPG.h

@@ -0,0 +1,48 @@
+//* 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 "INSBase.h"
+
+/**
+ *  Computes residual and Jacobian contributions for the PSPG stabilization term for mesh advection
+ **/
+class ConvectedMeshPSPG : public INSBase
+{
+public:
+  static InputParameters validParams();
+
+  ConvectedMeshPSPG(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+  RealVectorValue dStrongResidualDDisp(unsigned short component);
+  RealVectorValue dStrongResidualDVel(unsigned short component);
+
+  /**
+   * Compute the strong residual, e.g. -rho * ddisp/dt * grad(u)
+   */
+  RealVectorValue strongResidual();
+
+  const VariableValue & _disp_x_dot;
+  const VariableValue & _d_disp_x_dot;
+  const unsigned int _disp_x_id;
+  const VariableValue & _disp_y_dot;
+  const VariableValue & _d_disp_y_dot;
+  const unsigned int _disp_y_id;
+  const VariableValue & _disp_z_dot;
+  const VariableValue & _d_disp_z_dot;
+  const unsigned int _disp_z_id;
+
+  const MaterialProperty<Real> & _rho;
+};