Permalink
Browse files

Merge pull request #12657 from dschwen/adtm_12650

Cartesian AD stress divergence kernel and small strain linear elasticity
  • Loading branch information...
lindsayad committed Jan 11, 2019
2 parents 846968f + 7079e69 commit 1abf649307c970e593233366e909053126cb5402
Showing with 794 additions and 157 deletions.
  1. +7 −4 framework/include/base/MooseObject.h
  2. +15 −2 framework/include/kernels/ADKernel.h
  3. +10 −5 framework/include/materials/ADMaterial.h
  4. +4 −2 framework/src/kernels/ADKernel.C
  5. +3 −3 framework/src/materials/ADMaterial.C
  6. +1 −1 modules/combined/doc/content/index.md
  7. +25 −25 modules/porous_flow/doc/content/source/actions/PorousFlowBasicTHM.md
  8. +24 −0 modules/tensor_mechanics/doc/content/source/kernels/ADStressDivergenceTensors.md
  9. +17 −0 modules/tensor_mechanics/doc/content/source/materials/ADComputeLinearElasticStress.md
  10. +17 −0 modules/tensor_mechanics/doc/content/source/materials/ADComputeSmallStrain.md
  11. +58 −0 modules/tensor_mechanics/include/kernels/ADStressDivergenceTensors.h
  12. +37 −0 modules/tensor_mechanics/include/materials/ADComputeLinearElasticStress.h
  13. +35 −0 modules/tensor_mechanics/include/materials/ADComputeSmallStrain.h
  14. +71 −0 modules/tensor_mechanics/include/materials/ADComputeStrainBase.h
  15. +68 −0 modules/tensor_mechanics/include/materials/ADComputeStressBase.h
  16. +93 −0 modules/tensor_mechanics/src/kernels/ADStressDivergenceTensors.C
  17. +1 −1 modules/tensor_mechanics/src/kernels/StressDivergenceTensors.C
  18. +45 −0 modules/tensor_mechanics/src/materials/ADComputeLinearElasticStress.C
  19. +66 −0 modules/tensor_mechanics/src/materials/ADComputeSmallStrain.C
  20. +105 −0 modules/tensor_mechanics/src/materials/ADComputeStrainBase.C
  21. +56 −0 modules/tensor_mechanics/src/materials/ADComputeStressBase.C
  22. +1 −1 modules/tensor_mechanics/src/materials/ComputeElasticityBeam.C
  23. +1 −1 modules/tensor_mechanics/src/materials/ComputeElasticityTensorBase.C
  24. +1 −1 modules/tensor_mechanics/src/materials/ComputeElasticityTensorCP.C
  25. +6 −5 modules/tensor_mechanics/src/materials/ComputeMultiPlasticityStress.C
  26. +1 −1 modules/tensor_mechanics/src/materials/ComputeStrainBase.C
  27. +2 −2 modules/tensor_mechanics/src/materials/FluxBasedStrainIncrement.C
  28. +1 −1 modules/tensor_mechanics/src/materials/GBRelaxationStrainIncrement.C
  29. +3 −2 modules/tensor_mechanics/src/materials/StrainEnergyDensity.C
  30. +0 −31 modules/tensor_mechanics/test/src/kernels/ADStressDivergenceTest.C
  31. +0 −59 modules/tensor_mechanics/test/src/materials/TensorMechADMatTest.C
  32. +19 −9 modules/tensor_mechanics/test/tests/ad_simple_linear/linear-ad.i
  33. +1 −1 modules/tensor_mechanics/test/tests/ad_simple_linear/tests
@@ -26,6 +26,11 @@ InputParameters validParams<MooseObject>();
// needed to avoid #include cycle with MooseApp and MooseObject
[[noreturn]] void callMooseErrorRaw(std::string & msg, MooseApp * app);

// helper macro to explicitly instantiate AD classes
#define adBaseClass(X) \
template class X<RESIDUAL>; \
template class X<JACOBIAN>

#define adGetParam this->template getParam

/**
@@ -97,8 +102,7 @@ class MooseObject : public ConsoleStreamInterface, public libMesh::ParallelObjec
* back to the normal behavior of mooseError - only printing a message using the given args.
*/
template <typename... Args>
[[noreturn]] void paramError(const std::string & param, Args... args)
{
[[noreturn]] void paramError(const std::string & param, Args... args) {
auto prefix = param + ": ";
if (!_pars.inputLocation(param).empty())
prefix = _pars.inputLocation(param) + ": (" + _pars.paramFullpath(param) + "):\n";
@@ -137,8 +141,7 @@ class MooseObject : public ConsoleStreamInterface, public libMesh::ParallelObjec
}

template <typename... Args>
[[noreturn]] void mooseError(Args &&... args) const
{
[[noreturn]] void mooseError(Args &&... args) const {
std::ostringstream oss;
moose::internal::mooseStreamAll(oss, std::forward<Args>(args)...);
std::string msg = oss.str();
@@ -19,8 +19,21 @@
using ADKernel<compute_stage>::_u; \
using ADKernel<compute_stage>::_var; \
using ADKernel<compute_stage>::_grad_test; \
using ADKernel<compute_stage>::_grad_u

using ADKernel<compute_stage>::_grad_u; \
using ADKernel<compute_stage>::_JxW; \
using ADKernel<compute_stage>::_coord; \
using ADKernel<compute_stage>::_local_re; \
using ADKernel<compute_stage>::_qrule; \
using ADKernel<compute_stage>::_save_in; \
using ADKernel<compute_stage>::_has_save_in; \
using ADKernel<compute_stage>::_current_elem_volume; \
using ADKernel<compute_stage>::coupled; \
using ADKernel<compute_stage>::coupledComponents; \
using ADKernel<compute_stage>::getBlockCoordSystem; \
using ADKernel<compute_stage>::paramError; \
using ADKernel<compute_stage>::isParamValid

// forward declarations
template <ComputeStage compute_stage>
class ADKernel;

@@ -17,16 +17,21 @@
using ADMaterial<compute_stage>::_qp; \
using ADMaterial<compute_stage>::_ad_grad_zero; \
using ADMaterial<compute_stage>::_ad_zero; \
using ADMaterial<compute_stage>::coupledComponents
using ADMaterial<compute_stage>::_qrule; \
using ADMaterial<compute_stage>::_JxW; \
using ADMaterial<compute_stage>::_coord; \
using ADMaterial<compute_stage>::_assembly; \
using ADMaterial<compute_stage>::_mesh; \
using ADMaterial<compute_stage>::coupled; \
using ADMaterial<compute_stage>::coupledComponents; \
using ADMaterial<compute_stage>::isParamValid; \
using ADMaterial<compute_stage>::paramError

// forward declarations
template <ComputeStage>
class ADMaterial;

template <>
InputParameters validParams<ADMaterial<RESIDUAL>>();
template <>
InputParameters validParams<ADMaterial<JACOBIAN>>();
declareADValidParams(ADMaterial);

/**
* ADMaterials compute ADMaterialProperties.
@@ -89,6 +89,7 @@ ADKernel<compute_stage>::computeResidual()
{
prepareVectorTag(_assembly, _var.number());

precalculateResidual();
for (_i = 0; _i < _test.size(); _i++)
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
_local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
@@ -117,6 +118,7 @@ ADKernel<compute_stage>::computeJacobian()

size_t ad_offset = _var.number() * _sys.getMaxVarNDofsPerElem();

precalculateResidual();
for (_i = 0; _i < _test.size(); _i++)
{
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
@@ -203,5 +205,5 @@ ADKernel<compute_stage>::computeOffDiagJacobianScalar(unsigned int /*jvar*/)
*/
}

template class ADKernel<RESIDUAL>;
template class ADKernel<JACOBIAN>;
// explicit instantiation is required for AD base classes
adBaseClass(ADKernel);
@@ -22,13 +22,13 @@ template <>
InputParameters
validParams<ADMaterial<JACOBIAN>>()
{
return validParams<ADMaterial<JACOBIAN>>();
return validParams<ADMaterial<RESIDUAL>>();
}

template <ComputeStage compute_stage>
ADMaterial<compute_stage>::ADMaterial(const InputParameters & parameters) : Material(parameters)
{
}

template class ADMaterial<RESIDUAL>;
template class ADMaterial<JACOBIAN>;
// explicit instantiation is required for AD base classes
adBaseClass(ADMaterial);
@@ -1,4 +1,4 @@
# Combined Module

The combined module is for testing only, it should not be used as a module in the
traditional since.
traditional sense.
@@ -4,31 +4,31 @@ This action allows simple simulation of fully-saturated, single-phase,
single-component fluid flow. The fluid flow may be optionally coupled
to mechanics and/or heat flow using the `coupling_type` flag.

The fluid equation is a *simplified* form of the full [PorousFlow fluid equation](governing_equations.md) (see [PorousFlowFullySaturatedMassTimeDerivative](PorousFlowFullySaturatedMassTimeDerivative.md) and [PorousFlowFullySaturatedDarcyBase](PorousFlowFullySaturatedDarcyBase.md) for the derivation):
The fluid equation is a *simplified* form of the full [PorousFlow fluid equation](/porous_flow/governing_equations.md) (see [PorousFlowFullySaturatedMassTimeDerivative](/PorousFlowFullySaturatedMassTimeDerivative.md) and [PorousFlowFullySaturatedDarcyBase](/PorousFlowFullySaturatedDarcyBase.md) for the derivation):
\begin{equation}
\label{eq:basicthm}
0 = (\rho)\left(\frac{\dot{P}}{M} + \alpha_{B}\dot{\epsilon}_{v} - A\dot{T}\right) -
\nabla_{i}\left((\rho) k_{ij}\left(\nabla_{j}P - \rho g_{j}\right)/\mu\right)
\ .
\end{equation}
Note that the fluid-mass time derivative is close to linear, and is perfectly linear if `multiply_by_density=false`, and this also almost linearises the flow term. Extremely good nonlinear convergence should therefore be expected, but there are some knock-on effects that are documented in [PorousFlowFullySaturatedMassTimeDerivative](PorousFlowFullySaturatedMassTimeDerivative.md).
Note that the fluid-mass time derivative is close to linear, and is perfectly linear if `multiply_by_density=false`, and this also almost linearises the flow term. Extremely good nonlinear convergence should therefore be expected, but there are some knock-on effects that are documented in [PorousFlowFullySaturatedMassTimeDerivative](/PorousFlowFullySaturatedMassTimeDerivative.md).

In fully-saturated, single-phase simulations [upwinding](upwinding.md)
In fully-saturated, single-phase simulations [upwinding](/porous_flow/upwinding.md)
is typically unnecessary. Moreover, the standard PorousFlow Kernels
are somewhat inefficient in the fully-saturated case since Material
Properties such as relative permeabilities and saturations actually do
not need to be computed.

In fully-saturated, single-phase, single-component simulations, the
[mass lumping](lumping.md) is also typically unncessary. More
[mass lumping](/porous_flow/mass_lumping.md) is also typically unncessary. More
importantly, in many real-life situations, as may be seen in
[eq:basicthm] the time derivative of the fluid mass may be linearised,
which leads to improved convergence.

To simulate [eq:basicthm] the `PorousFlowBasicTHM` Action employs the following Kernels:

- [PorousFlowFullySaturatedMassTimeDerivative](PorousFlowFullySaturatedMassTimeDerivative.md), if the simulation is of Transient type.
- [PorousFlowFullySaturatedDarcyBase](PorousFlowFullySaturatedDarcyBase.md)
- [PorousFlowFullySaturatedMassTimeDerivative](/PorousFlowFullySaturatedMassTimeDerivative.md), if the simulation is of Transient type.
- [PorousFlowFullySaturatedDarcyBase](/PorousFlowFullySaturatedDarcyBase.md)

For isothermal simulations, the fluid properties may still depend on temperature, so the `temperature` input parameter may be set to any real number, or to an AuxVariable if desired.

@@ -38,42 +38,42 @@ For anisothermal simulations, the energy equation reads
\end{equation}
where the final term is only used if coupling with mechanics is also desired. To simulate this DE, `PorousFlowBasicTHM` uses the following kernels:

- [PorousFlowEnergyTimeDerivative](PorousFlowEnergyTimeDerivative.md) if the simulation is of Transient type
- [PorousFlowFullySaturatedHeatAdvection](PorousFlowFullySaturatedHeatAdvection.md) with `multiply_by_density=true` (irrespective of the setting of this flag in the `PorousflowBasicTHM` Action)
- [PorousFlowHeatConduction](PorousFlowHeatConduction.md)
- [PorousFlowHeatVolumetricExpansion](PorousFlowHeatVolumetricExpansion.md) if coupling with temperature and mechanics is desired, and if the simulation is of Transient type
- [PorousFlowEnergyTimeDerivative](/PorousFlowEnergyTimeDerivative.md) if the simulation is of Transient type
- [PorousFlowFullySaturatedHeatAdvection](/PorousFlowFullySaturatedHeatAdvection.md) with `multiply_by_density=true` (irrespective of the setting of this flag in the `PorousflowBasicTHM` Action)
- [PorousFlowHeatConduction](/PorousFlowHeatConduction.md)
- [PorousFlowHeatVolumetricExpansion](/PorousFlowHeatVolumetricExpansion.md) if coupling with temperature and mechanics is desired, and if the simulation is of Transient type

For simulations that couple fluid flow to mechanics, the equations are already written in [governing equations](governing_equations.md), and `PorousFlowBasicTHM` implements these by using the following kernels:
For simulations that couple fluid flow to mechanics, the equations are already written in [governing equations](/porous_flow/governing_equations.md), and `PorousFlowBasicTHM` implements these by using the following kernels:

- [StressDivergenceTensors](/StressDivergenceTensors.md)
- [Gravity](/Gravity.md)
- [PorousFlowEffectiveStressCoupling](PorousFlowEffectiveStressCoupling.md)
- [PorousFlowEffectiveStressCoupling](/PorousFlowEffectiveStressCoupling.md)

`PorousFlowBasicTHM` adds many Materials automatically, however, to run a simulation you will need to provide more Materials for each mesh block, depending on your simulation type, viz:

- One of the `PorousFlowPermeability` classes, eg [PorousFlowPermeabilityConst](PorousFlowPermeabilityConst.md)
- [PorousFlowConstantBiotModulus](PorousFlowConstantBiotModulus.md)
- [PorousFlowConstantThermalExpansionCoefficient](PorousFlowConstantThermalExpansionCoefficient.md)
- [PorousFlowPorosity](PorousFlowPorosity.md)
- [PorousFlowMatrixInternalEnergy](PorousFlowMatrixInternalEnergy.md)
- A thermal conductivity calculator, eg [PorousFlowThermalConductivityIdeal](PorousFlowThermalConductivityIdeal.md)
- An elasticity tensor, eg, [ComputeIsotropicElasticityTensor](ComputeIsotropicElasticityTensor.md)
- A strain calculator, eg, [ComputeSmallStrain](ComputeSmallStrain.md)
- A thermal expansion eigenstrain calculator, eg, [ComputeThermalExpansionEigenstrain](ComputeThermalExpansionEigenstrain.md)
- A stress calculator, eg [ComputeLinearElasticStress](ComputeLinearElasticStress.md)
- [PorousFlowConstantBiotModulus](/PorousFlowConstantBiotModulus.md)
- [PorousFlowConstantThermalExpansionCoefficient](/PorousFlowConstantThermalExpansionCoefficient.md)
- [PorousFlowPorosity](/PorousFlowPorosity.md)
- [PorousFlowMatrixInternalEnergy](/PorousFlowMatrixInternalEnergy.md)
- A thermal conductivity calculator, eg [PorousFlowThermalConductivityIdeal](/PorousFlowThermalConductivityIdeal.md)
- An elasticity tensor, eg, [ComputeIsotropicElasticityTensor](/ComputeIsotropicElasticityTensor.md)
- A strain calculator, eg, [ComputeSmallStrain](/ComputeSmallStrain.md)
- A thermal expansion eigenstrain calculator, eg, [ComputeThermalExpansionEigenstrain](/ComputeThermalExpansionEigenstrain.md)
- A stress calculator, eg [ComputeLinearElasticStress](/ComputeLinearElasticStress.md)

!alert note
Since upwinding and no mass lumping of the fluid mass are used (for simplicity, efficiency and to reduce [numerical diffusion](numerical_diffusion.md)), the results may be slightly different to simulations that employ upwinding and mass lumping.
Since upwinding and no mass lumping of the fluid mass are used (for simplicity, efficiency and to reduce [numerical diffusion](/porous_flow/numerical_diffusion.md)), the results may be slightly different to simulations that employ upwinding and mass lumping.

A simple example of using `PorousFlowBasicTHM` is documented in the [PorousFlow tutorial](tutorial_01.md) with input file:
A simple example of using `PorousFlowBasicTHM` is documented in the [PorousFlow tutorial](/porous_flow/tutorial_01.md) with input file:

!listing modules/porous_flow/examples/tutorial/01.i

A TH example that uses `PorousFlowBasicTHM` is also documented in the [PorousFlow tutorial](tutorial_03.md) with input file:
A TH example that uses `PorousFlowBasicTHM` is also documented in the [PorousFlow tutorial](/porous_flow/tutorial_03.md) with input file:

!listing modules/porous_flow/examples/tutorial/03.i

A THM example that uses `PorousFlowBasicTHM` is also documented in the [PorousFlow tutorial](tutorial_04.md) with input file:
A THM example that uses `PorousFlowBasicTHM` is also documented in the [PorousFlow tutorial](/porous_flow/tutorial_04.md) with input file:

!listing modules/porous_flow/examples/tutorial/04.i

@@ -0,0 +1,24 @@
# Stress Divergence Tensors with Automatic Differentiation

!syntax description /ADKernels/ADStressDivergenceTensors<RESIDUAL>

## Description

The `ADStressDivergenceTensors` kernel calculates the residual of the stress
divergence for 1D, 2D, and 3D problems in the Cartesian coordinate system.

The Jacobian in `ADStressDivergenceTensors` is computed using forward automatic
differentiation.

## Residual Calculation

!include modules/tensor_mechanics/common/supplementalStressDivergenceKernels.md

Either 1, 2, or 3 displacement variables can be used in the stress divergence
calculator for the Cartesian system.

!syntax parameters /ADKernels/ADStressDivergenceTensors<RESIDUAL>

!syntax inputs /ADKernels/ADStressDivergenceTensors<RESIDUAL>

!syntax children /ADKernels/ADStressDivergenceTensors<RESIDUAL>
@@ -0,0 +1,17 @@
# ADComputeLinearElasticStress

!syntax description /ADMaterials/ADComputeLinearElasticStress<RESIDUAL>

This material is a version of
[ComputeLinearElasticStress](/materials/ComputeLinearElasticStress.md) which
uses forward mode automatic differentiation to carry along derivatives of its
material properties for use with the automatic differentiation tensor mechanics
kernels (such as [ADStressDivergence](/kernels/ADStressDivergenceTensors.md)).

!syntax parameters /ADMaterials/ADComputeLinearElasticStress<RESIDUAL>

!syntax inputs /ADMaterials/ADComputeLinearElasticStress<RESIDUAL>

!syntax children /ADMaterials/ADComputeLinearElasticStress<RESIDUAL>

!bibtex bibliography
@@ -0,0 +1,17 @@
# ADComputeSmallStrain

!syntax description /ADMaterials/ADComputeSmallStrain<RESIDUAL>

This material is a version of
[ComputeSmallStrain](/materials/ComputeSmallStrain.md) which
uses forward mode automatic differentiation to carry along derivatives of its
material properties for use with the automatic differentiation tensor mechanics
kernels (such as [ADStressDivergence](/kernels/ADStressDivergenceTensors.md)).

!syntax parameters /ADMaterials/ADComputeSmallStrain<RESIDUAL>

!syntax inputs /ADMaterials/ADComputeSmallStrain<RESIDUAL>

!syntax children /ADMaterials/ADComputeSmallStrain<RESIDUAL>

!bibtex bibliography
@@ -0,0 +1,58 @@
//* 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

#ifndef ADSTRESSDIVERGENCETENSORS_H
#define ADSTRESSDIVERGENCETENSORS_H

#include "ADKernel.h"

// Forward Declarations
template <ComputeStage>
class ADStressDivergenceTensors;
template <typename>
class RankTwoTensorTempl;
typedef RankTwoTensorTempl<Real> RankTwoTensor;
typedef RankTwoTensorTempl<ADReal> ADRankTwoTensor;

declareADValidParams(ADStressDivergenceTensors);

/**
* ADStressDivergenceTensors is the automatic differentiation version of StressDivergenceTensors
*/
template <ComputeStage compute_stage>
class ADStressDivergenceTensors : public ADKernel<compute_stage>
{
public:
ADStressDivergenceTensors(const InputParameters & parameters);

protected:
virtual void initialSetup() override;

virtual ADResidual computeQpResidual() override;
virtual void precalculateResidual() override;

const std::string _base_name;

const ADMaterialProperty(RankTwoTensor) & _stress;
const unsigned int _component;

/// Coupled displacement variables
const unsigned int _ndisp;
std::vector<unsigned int> _disp_var;

/// Gradient of test function averaged over the element. Used in volumetric locking correction calculation.
std::vector<typename RealVectorValueType<compute_stage>::type> _avg_grad_test;

/// Flag for volumetric locking correction
const bool _volumetric_locking_correction;

usingKernelMembers;
};

#endif // ADSTRESSDIVERGENCETENSORS_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

#ifndef ADCOMPUTELINEARELASTICSTRESS_H
#define ADCOMPUTELINEARELASTICSTRESS_H

#include "ADComputeStressBase.h"

template <ComputeStage>
class ADComputeLinearElasticStress;

declareADValidParams(ADComputeLinearElasticStress);

/**
* ADComputeLinearElasticStress computes the stress following linear elasticity theory (small
* strains)
*/
template <ComputeStage compute_stage>
class ADComputeLinearElasticStress : public ADComputeStressBase<compute_stage>
{
public:
ADComputeLinearElasticStress(const InputParameters & parameters);
virtual void initialSetup();

protected:
virtual void computeQpStress();

usingComputeStressBaseMembers;
};

#endif // ADCOMPUTELINEARELASTICSTRESS_H
Oops, something went wrong.

0 comments on commit 1abf649

Please sign in to comment.