diff --git a/include/kernels/ADStokesIncompressibility.h b/include/kernels/ADStokesIncompressibility.h index a17b405..6f78643 100644 --- a/include/kernels/ADStokesIncompressibility.h +++ b/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; }; \ No newline at end of file diff --git a/include/kernels/ADStokesStressDivergence.h b/include/kernels/ADStokesStressDivergence.h index 3d56455..964a80e 100644 --- a/include/kernels/ADStokesStressDivergence.h +++ b/include/kernels/ADStokesStressDivergence.h @@ -22,4 +22,5 @@ class ADStokesStressDivergence : public ADVectorKernel virtual ADReal computeQpResidual() override; const ADMaterialProperty & _stress; + const ADVariableValue & _pressure; }; \ No newline at end of file diff --git a/include/materials/StokesLinearViscous.h b/include/materials/StokesLinearViscous.h new file mode 100644 index 0000000..d28e5e1 --- /dev/null +++ b/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 +{ +public: + static InputParameters validParams(); + + StokesLinearViscous(const InputParameters & parameters); + +protected: + /// @brief Calculate the strain + void computeQpProperties() override; + + /// The stress + ADMaterialProperty & _stress; + + /// The strain rate + const ADMaterialProperty & _strain_rate; + + /// Power law prefactor + const ADMaterialProperty & _mu; +}; \ No newline at end of file diff --git a/include/materials/StokesPowerLawCreep.h b/include/materials/StokesPowerLawCreep.h new file mode 100644 index 0000000..32699ec --- /dev/null +++ b/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 +{ +public: + static InputParameters validParams(); + + StokesPowerLawCreep(const InputParameters & parameters); + +protected: + /// @brief Calculate the strain + void computeQpProperties() override; + + /// The stress + ADMaterialProperty & _stress; + + /// The strain rate + const ADMaterialProperty & _strain_rate; + + /// Power law prefactor + const ADMaterialProperty & _A; + + /// Power law exponent + const ADMaterialProperty & _n; +}; \ No newline at end of file diff --git a/include/materials/StokesStrainRate.h b/include/materials/StokesStrainRate.h new file mode 100644 index 0000000..1114602 --- /dev/null +++ b/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 +{ +public: + static InputParameters validParams(); + + StokesStrainRate(const InputParameters & parameters); + +protected: + /// @brief Calculate the strain + void computeQpProperties() override; + + /// The strain rate + ADMaterialProperty & _strain_rate; + + /// The coupled velocity gradient + const ADVectorVariableGradient & _grad_vel; +}; \ No newline at end of file diff --git a/include/postprocessors/ScalarPrincipalValues.h b/include/postprocessors/ScalarPrincipalValues.h index c9ab500..39710a8 100644 --- a/include/postprocessors/ScalarPrincipalValues.h +++ b/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; }; diff --git a/include/postprocessors/SideExtremePostprocessor.h b/include/postprocessors/SideExtremePostprocessor.h index 20254c3..f6e0c93 100644 --- a/include/postprocessors/SideExtremePostprocessor.h +++ b/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; diff --git a/include/postprocessors/TimeDerivativePostprocessor.h b/include/postprocessors/TimeDerivativePostprocessor.h index 0aa51ab..b0291d5 100644 --- a/include/postprocessors/TimeDerivativePostprocessor.h +++ b/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 diff --git a/include/utils/miconoss/miconosstype.h b/include/utils/miconoss/miconosstype.h index 048f5fc..f90c990 100644 --- a/include/utils/miconoss/miconosstype.h +++ b/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 diff --git a/src/kernels/ADStokesIncompressibility.C b/src/kernels/ADStokesIncompressibility.C index f50db5b..ee0aeb3 100644 --- a/src/kernels/ADStokesIncompressibility.C +++ b/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]; } \ No newline at end of file diff --git a/src/kernels/ADStokesStressDivergence.C b/src/kernels/ADStokesStressDivergence.C index d701a4c..a566bcb 100644 --- a/src/kernels/ADStokesStressDivergence.C +++ b/src/kernels/ADStokesStressDivergence.C @@ -9,7 +9,7 @@ #include "ADStokesStressDivergence.h" -registerMooseObject("MooseApp", ADStokesStressDivergence); +registerMooseObject("DeerApp", ADStokesStressDivergence); InputParameters ADStokesStressDivergence::validParams() @@ -17,16 +17,22 @@ 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("stress")) + : ADVectorKernel(parameters), + _stress(getADMaterialPropertyByName("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]; } \ No newline at end of file diff --git a/src/materials/StokesLinearViscous.C b/src/materials/StokesLinearViscous.C new file mode 100644 index 0000000..4534f93 --- /dev/null +++ b/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::validParams(); + params.addClassDescription("Calculate the stress for Stokes flow using the simple linear model."); + + params.addRequiredParam("mu", "Viscosity"); + + params.suppressParameter("use_displaced_mesh"); + return params; +} + +StokesLinearViscous::StokesLinearViscous(const InputParameters & parameters) + : DerivativeMaterialInterface(parameters), + _stress(declareADProperty("stress")), + _strain_rate(getADMaterialPropertyByName("strain_rate")), + _mu(getADMaterialProperty("mu")) +{ +} + +void +StokesLinearViscous::computeQpProperties() +{ + _stress[_qp] = _strain_rate[_qp] * _mu[_qp]; +} \ No newline at end of file diff --git a/src/materials/StokesPowerLawCreep.C b/src/materials/StokesPowerLawCreep.C new file mode 100644 index 0000000..3536e63 --- /dev/null +++ b/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::validParams(); + params.addClassDescription("Calculate the stress for Stokes flow using a power law creep model"); + + params.addRequiredParam("A", "Creep prefactor"); + params.addRequiredParam("n", "Creep exponent"); + + params.suppressParameter("use_displaced_mesh"); + return params; +} + +StokesPowerLawCreep::StokesPowerLawCreep(const InputParameters & parameters) + : DerivativeMaterialInterface(parameters), + _stress(declareADProperty("stress")), + _strain_rate(getADMaterialPropertyByName("strain_rate")), + _A(getADMaterialProperty("A")), + _n(getADMaterialProperty("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)); +} \ No newline at end of file diff --git a/src/materials/StokesStrainRate.C b/src/materials/StokesStrainRate.C new file mode 100644 index 0000000..d095a32 --- /dev/null +++ b/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::validParams(); + params.addClassDescription("Calculate the strain rate as the symmetric gradient of the velocity"); + params.addRequiredCoupledVar("velocity", "The vector-valued velocity"); + params.suppressParameter("use_displaced_mesh"); + return params; +} + +StokesStrainRate::StokesStrainRate(const InputParameters & parameters) + : DerivativeMaterialInterface(parameters), + _strain_rate(declareADProperty("strain_rate")), + _grad_vel(adCoupledVectorGradient("velocity")) + +{ +} + +void +StokesStrainRate::computeQpProperties() +{ + _strain_rate[_qp] = 0.5 * (_grad_vel[_qp] + _grad_vel[_qp].transpose()); +} \ No newline at end of file diff --git a/src/postprocessors/ScalarPrincipalValues.C b/src/postprocessors/ScalarPrincipalValues.C index 1e65683..2869cb0 100644 --- a/src/postprocessors/ScalarPrincipalValues.C +++ b/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("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) diff --git a/src/utils/miconoss/miconosstype.C b/src/utils/miconoss/miconosstype.C index 66bb1e8..1c88c8f 100644 --- a/src/utils/miconoss/miconosstype.C +++ b/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"; diff --git a/tests/stokes/flow.i b/tests/stokes/flow.i new file mode 100644 index 0000000..cf4a2f5 --- /dev/null +++ b/tests/stokes/flow.i @@ -0,0 +1,155 @@ +[Mesh] + [mesh] + type = GeneratedMeshGenerator + dim = 3 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + nx = 5 + ny = 5 + nz = 5 + elem_type = HEX20 + [] +[] + +[Variables] + [p] + order = FIRST + family = LAGRANGE + [] + [u] + order = SECOND + family = LAGRANGE_VEC + [] +[] + +[AuxVariables] + [stress_xx] + family = MONOMIAL + order = CONSTANT + [] + [stress_yy] + family = MONOMIAL + order = CONSTANT + [] + [stress_zz] + family = MONOMIAL + order = CONSTANT + [] +[] + +[AuxKernels] + [stress_xx] + type = ADMaterialRankTwoTensorAux + variable = stress_xx + property = 'stress' + i = 0 + j = 0 + [] + [stress_yy] + type = ADMaterialRankTwoTensorAux + variable = stress_yy + property = 'stress' + i = 1 + j = 1 + [] + [stress_zz] + type = ADMaterialRankTwoTensorAux + variable = stress_zz + property = 'stress' + i = 2 + j = 2 + [] +[] + +[ICs] + [u] + type = VectorConstantIC + variable = u + x_value = 1e-15 + y_value = 1e-15 + z_value = 1e-15 + [] +[] + +[Kernels] + [equil] + type = ADStokesStressDivergence + variable = u + pressure = p + [] + [incompressible] + type = ADStokesIncompressibility + variable = p + velocity = u + [] +[] + +[Materials] + [strain] + type = StokesStrainRate + velocity = u + [] + [stress] + type = StokesLinearViscous + mu = 1.0 + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = u + boundary = 'bottom top' + [] + [inflow] + type = ADVectorFunctionDirichletBC + variable = u + boundary = 'right' + function_x = '-sin(pi * y)' + [] +[] + +[Preconditioning] + [FSP] + type = FSP + topsplit = 'up' + [up] + splitting = 'u p' + splitting_type = schur + petsc_options_iname = '-pc_fieldsplit_schur_fact_type -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side' + petsc_options_value = 'full self 300 fgmres right' + [] + [u] + vars = 'u' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + [] + [p] + vars = 'p' + petsc_options = '-pc_lsc_scale_diag' + petsc_options_iname = '-ksp_type -ksp_gmres_restart -pc_type -ksp_pc_side -lsc_pc_type -lsc_ksp_type -lsc_ksp_pc_side -lsc_ksp_gmres_restart' + petsc_options_value = 'fgmres 300 lsc right lu gmres right 300' + [] + [] +[] + +[Executioner] + type = Steady + + solve_type = 'newton' + line_search = none + + l_max_its = 10 + l_tol = 1e-14 + nl_max_its = 15 + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-6 +[] + +[Outputs] + exodus = true +[] diff --git a/tests/stokes/gold/flow_out.e b/tests/stokes/gold/flow_out.e new file mode 100644 index 0000000..9a08444 Binary files /dev/null and b/tests/stokes/gold/flow_out.e differ diff --git a/tests/stokes/tests b/tests/stokes/tests new file mode 100644 index 0000000..3ca9200 --- /dev/null +++ b/tests/stokes/tests @@ -0,0 +1,8 @@ +[Tests] + [./flow] + type = Exodiff + input = 'flow.i' + exodiff = 'flow_out.e' + requirement = "Stokes flow gives the right result for simple creeping flow" + [../] +[]