Potentially working stokes flow, need a manufactured solution

include/kernels/ADStokesIncompressibility.h

@@ -21,5 +21,5 @@ class ADStokesIncompressibility : public ADKernel
 protected:
   virtual ADReal computeQpResidual() override;
 
-  const ADVectorVariableGradient & _grad_disp;
+  const ADVectorVariableGradient & _grad_vel;
 };

include/kernels/ADStokesStressDivergence.h

@@ -22,4 +22,5 @@ class ADStokesStressDivergence : public ADVectorKernel
   virtual ADReal computeQpResidual() override;
 
   const ADMaterialProperty<RankTwoTensor> & _stress;
+  const ADVariableValue & _pressure;
 };

include/materials/StokesLinearViscous.h

@@ -0,0 +1,37 @@
+//* 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 "Material.h"
+#include "DerivativeMaterialInterface.h"
+
+/**
+ * Calculate the simple small strain rate
+ */
+class StokesLinearViscous : public DerivativeMaterialInterface<Material>
+{
+public:
+  static InputParameters validParams();
+
+  StokesLinearViscous(const InputParameters & parameters);
+
+protected:
+  /// @brief  Calculate the strain
+  void computeQpProperties() override;
+
+  /// The stress
+  ADMaterialProperty<RankTwoTensor> & _stress;
+
+  /// The strain rate
+  const ADMaterialProperty<RankTwoTensor> & _strain_rate;
+
+  /// Power law prefactor
+  const ADMaterialProperty<Real> & _mu;
+};

include/materials/StokesPowerLawCreep.h

@@ -0,0 +1,40 @@
+//* 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 "Material.h"
+#include "DerivativeMaterialInterface.h"
+
+/**
+ * Calculate the simple small strain rate
+ */
+class StokesPowerLawCreep : public DerivativeMaterialInterface<Material>
+{
+public:
+  static InputParameters validParams();
+
+  StokesPowerLawCreep(const InputParameters & parameters);
+
+protected:
+  /// @brief  Calculate the strain
+  void computeQpProperties() override;
+
+  /// The stress
+  ADMaterialProperty<RankTwoTensor> & _stress;
+
+  /// The strain rate
+  const ADMaterialProperty<RankTwoTensor> & _strain_rate;
+
+  /// Power law prefactor
+  const ADMaterialProperty<Real> & _A;
+
+  /// Power law exponent
+  const ADMaterialProperty<Real> & _n;
+};

include/materials/StokesStrainRate.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 "Material.h"
+#include "DerivativeMaterialInterface.h"
+
+/**
+ * Calculate the simple small strain rate
+ */
+class StokesStrainRate : public DerivativeMaterialInterface<Material>
+{
+public:
+  static InputParameters validParams();
+
+  StokesStrainRate(const InputParameters & parameters);
+
+protected:
+  /// @brief  Calculate the strain
+  void computeQpProperties() override;
+
+  /// The strain rate
+  ADMaterialProperty<RankTwoTensor> & _strain_rate;
+
+  /// The coupled velocity gradient
+  const ADVectorVariableGradient & _grad_vel;
+};

include/postprocessors/ScalarPrincipalValues.h

@@ -20,11 +20,10 @@ class ScalarPrincipalValues : public GeneralPostprocessor
 
   virtual void initialize() override;
   virtual void execute() override;
-  virtual PostprocessorValue getValue() const;
+  virtual PostprocessorValue getValue() const override;
 
 protected:
   const VariableValue & _scalar_var;
-  const Order & _scalar_order;
   const size_t _rank;
   PostprocessorValue _value;
 };

include/postprocessors/SideExtremePostprocessor.h

@@ -17,7 +17,7 @@ class SideExtremePostprocessor : public SidePostprocessor
 
   virtual void initialize() override;
   virtual void execute() override;
-  virtual PostprocessorValue getValue() const;
+  virtual PostprocessorValue getValue() const override;
   virtual void threadJoin(const UserObject & y) override;
   virtual void finalize() override;
 

include/postprocessors/TimeDerivativePostprocessor.h

@@ -26,7 +26,7 @@ class TimeDerivativePostprocessor : public GeneralPostprocessor
 
   virtual void initialize() override;
   virtual void execute() override;
-  virtual PostprocessorValue getValue() const;
+  virtual PostprocessorValue getValue() const override;
 
 protected:
   /// cumulative sum of the post-processor value

include/utils/miconoss/miconosstype.h

@@ -26,7 +26,7 @@ void printMatrix(const matrixD & M, const std::string & Mname);
 namespace miconossmath
 {
 
-extern "C" void dgesv_(int *, int *, double *, int *, int *, double *, int *, int *);
+extern "C"  void dgesv_(const int & n, const int & nrhs, double * A, const int & lda, int * ipiv, double * b, const int & ldb, int & info);
 
 /// the allowed norm types
 enum normtype

src/kernels/ADStokesIncompressibility.C

@@ -9,26 +9,26 @@
 
 #include "ADStokesIncompressibility.h"
 
-registerMooseObject("MooseApp", ADStokesIncompressibility);
+registerMooseObject("DeerApp", ADStokesIncompressibility);
 
 InputParameters
 ADStokesIncompressibility::validParams()
 {
-  InputParameters params = ADVectorKernel::validParams();
+  InputParameters params = ADKernel::validParams();
   params.addClassDescription("The incompressibility equation for Stokes flow");
 
-  params.addRequiredCoupledVar("displacement", "The vector-valued displacement");
+  params.addRequiredCoupledVar("velocity", "The vector-valued velocity");
 
   return params;
 }
 
 ADStokesIncompressibility::ADStokesIncompressibility(const InputParameters & parameters)
-  : ADVectorKernel(parameters), _grad_disp(adCoupledVectorGradient("displacement"))
+  : ADKernel(parameters), _grad_vel(adCoupledVectorGradient("velocity"))
 {
 }
 
 ADReal
 ADStokesIncompressibility::computeQpResidual()
 {
-  return (_grad_disp[_qp](0, 0) + _grad_disp[_qp](1, 1) + _grad_disp[_qp](2, 2)) * _phi[_qp];
+  return (_grad_vel[_qp](0, 0) + _grad_vel[_qp](1, 1) + _grad_vel[_qp](2, 2)) * _test[_i][_qp];
 }

src/kernels/ADStokesStressDivergence.C

@@ -9,24 +9,30 @@
 
 #include "ADStokesStressDivergence.h"
 
-registerMooseObject("MooseApp", ADStokesStressDivergence);
+registerMooseObject("DeerApp", ADStokesStressDivergence);
 
 InputParameters
 ADStokesStressDivergence::validParams()
 {
   InputParameters params = ADVectorKernel::validParams();
   params.addClassDescription("The unstabilized stress diveregence kernel for Stokes flow");
 
+  params.addRequiredCoupledVar("pressure", "The pressure");
+
   return params;
 }
 
 ADStokesStressDivergence::ADStokesStressDivergence(const InputParameters & parameters)
-  : ADVectorKernel(parameters), _stress(getADMaterialPropertyByName<RankTwoTensor>("stress"))
+  : ADVectorKernel(parameters),
+    _stress(getADMaterialPropertyByName<RankTwoTensor>("stress")),
+    _pressure(adCoupledValue("pressure"))
 {
 }
 
 ADReal
 ADStokesStressDivergence::computeQpResidual()
 {
-  return _stress[_qp].contract(_grad_phi[_qp]);
+  return _stress[_qp].contract(_grad_test[_i][_qp]) -
+         (_grad_test[_i][_qp](0, 0) + _grad_test[_i][_qp](1, 1) + _grad_test[_i][_qp](2, 2)) *
+             _pressure[_qp];
 }

src/materials/StokesLinearViscous.C

@@ -0,0 +1,37 @@
+//* 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 "StokesLinearViscous.h"
+
+registerMooseObject("DeerApp", StokesLinearViscous);
+
+InputParameters
+StokesLinearViscous::validParams()
+{
+  InputParameters params = DerivativeMaterialInterface<Material>::validParams();
+  params.addClassDescription("Calculate the stress for Stokes flow using the simple linear model.");
+
+  params.addRequiredParam<MaterialPropertyName>("mu", "Viscosity");
+
+  params.suppressParameter<bool>("use_displaced_mesh");
+  return params;
+}
+
+StokesLinearViscous::StokesLinearViscous(const InputParameters & parameters)
+  : DerivativeMaterialInterface<Material>(parameters),
+    _stress(declareADProperty<RankTwoTensor>("stress")),
+    _strain_rate(getADMaterialPropertyByName<RankTwoTensor>("strain_rate")),
+    _mu(getADMaterialProperty<Real>("mu"))
+{
+}
+
+void
+StokesLinearViscous::computeQpProperties()
+{
+  _stress[_qp] = _strain_rate[_qp] * _mu[_qp];
+}

src/materials/StokesPowerLawCreep.C

@@ -0,0 +1,46 @@
+//* 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 "StokesPowerLawCreep.h"
+
+registerMooseObject("DeerApp", StokesPowerLawCreep);
+
+InputParameters
+StokesPowerLawCreep::validParams()
+{
+  InputParameters params = DerivativeMaterialInterface<Material>::validParams();
+  params.addClassDescription("Calculate the stress for Stokes flow using a power law creep model");
+
+  params.addRequiredParam<MaterialPropertyName>("A", "Creep prefactor");
+  params.addRequiredParam<MaterialPropertyName>("n", "Creep exponent");
+
+  params.suppressParameter<bool>("use_displaced_mesh");
+  return params;
+}
+
+StokesPowerLawCreep::StokesPowerLawCreep(const InputParameters & parameters)
+  : DerivativeMaterialInterface<Material>(parameters),
+    _stress(declareADProperty<RankTwoTensor>("stress")),
+    _strain_rate(getADMaterialPropertyByName<RankTwoTensor>("strain_rate")),
+    _A(getADMaterialProperty<Real>("A")),
+    _n(getADMaterialProperty<Real>("n"))
+{
+}
+
+void
+StokesPowerLawCreep::computeQpProperties()
+{
+  auto A = _A[_qp];
+  auto n = _n[_qp];
+  auto dev_strain_rate = _strain_rate[_qp].deviatoric();
+
+  auto effective_strain_rate = std::sqrt(2.0 / 3.0 * dev_strain_rate.contract(dev_strain_rate));
+  auto effective_stress = std::pow(effective_strain_rate / A, 1.0 / n);
+
+  _stress[_qp] = 2.0 / 3.0 * dev_strain_rate / (A * std::pow(effective_stress, n - 1.0));
+}

src/materials/StokesStrainRate.C

@@ -0,0 +1,35 @@
+//* 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 "StokesStrainRate.h"
+
+registerMooseObject("DeerApp", StokesStrainRate);
+
+InputParameters
+StokesStrainRate::validParams()
+{
+  InputParameters params = DerivativeMaterialInterface<Material>::validParams();
+  params.addClassDescription("Calculate the strain rate as the symmetric gradient of the velocity");
+  params.addRequiredCoupledVar("velocity", "The vector-valued velocity");
+  params.suppressParameter<bool>("use_displaced_mesh");
+  return params;
+}
+
+StokesStrainRate::StokesStrainRate(const InputParameters & parameters)
+  : DerivativeMaterialInterface<Material>(parameters),
+    _strain_rate(declareADProperty<RankTwoTensor>("strain_rate")),
+    _grad_vel(adCoupledVectorGradient("velocity"))
+
+{
+}
+
+void
+StokesStrainRate::computeQpProperties()
+{
+  _strain_rate[_qp] = 0.5 * (_grad_vel[_qp] + _grad_vel[_qp].transpose());
+}

src/postprocessors/ScalarPrincipalValues.C

@@ -25,11 +25,10 @@ ScalarPrincipalValues::validParams()
 ScalarPrincipalValues::ScalarPrincipalValues(const InputParameters & parameters)
   : GeneralPostprocessor(parameters),
     _scalar_var(coupledScalarValue("scalar_variable")),
-    _scalar_order(coupledScalarOrder("scalar_variable")),
     _rank(getParam<size_t>("rank"))
 {
   // Could add general case later
-  if (_scalar_order != 6)
+  if (coupledScalarOrder("scalar_variable") != 6)
     mooseError("Coupled scalar variable needs to have order 6");
 
   if (_rank < 0)

src/utils/miconoss/miconosstype.C

@@ -115,7 +115,7 @@ solveAxb(const matrixD & A, const vecD & b, const uint & syssize, vecD & b_out)
       k += 1;
     }
 
-  dgesv_(&neq, &nrhs, &jac[0], &neq, &IPIV[0], &b_out[0], &neq, &info);
+  dgesv_(neq, nrhs, &jac[0], neq, &IPIV[0], &b_out[0], neq, info);
 
   if (info < 0)
     std::cerr << "solveAxb the " + std::to_string(info) + "the argument has an illegal value";
@@ -154,7 +154,7 @@ solveAxNb(const matrixD & A, const matrixD & b, const uint & syssize, matrixD &
       k += 1;
     }
 
-  dgesv_(&neq, &nrhs, &A_to_use[0], &neq, &IPIV[0], &b_to_use[0], &neq, &info);
+  dgesv_(neq, nrhs, &A_to_use[0], neq, &IPIV[0], &b_to_use[0], neq, info);
 
   if (info < 0)
     std::cerr << "solveAxNb the " + std::to_string(-info) + "th argument has an illegal value";