Permalink
Browse files

Merge pull request #12599 from lindsayad/ad-tensor-mech-ex

Add TensorMechanics automatic differentiation example
  • Loading branch information...
dschwen committed Dec 14, 2018
2 parents 7ae2ab3 + 14a5d40 commit 743232c179835f77e5220541bc3e140234e10d77
Showing with 2,403 additions and 940 deletions.
  1. +3 −1 framework/include/functions/PiecewiseBilinear.h
  2. +38 −2 framework/include/interfaces/Coupleable.h
  3. +5 −1 framework/include/materials/ADMaterial.h
  4. +4 −4 framework/include/materials/MaterialData.h
  5. +31 −31 framework/include/materials/MaterialProperty.h
  6. +1 −1 framework/include/materials/MaterialPropertyStorage.h
  7. +59 −13 framework/include/restart/DataIO.h
  8. +67 −0 framework/include/utils/ADReal.h
  9. +192 −145 framework/include/utils/ColumnMajorMatrix.h
  10. +182 −50 framework/include/utils/MooseADWrapper.h
  11. +62 −53 framework/include/utils/MooseTypes.h
  12. +23 −6 framework/include/utils/MooseUtils.h
  13. +138 −54 framework/include/utils/RankFourTensor.h
  14. +6 −2 framework/include/utils/RankThreeTensor.h
  15. +133 −102 framework/include/utils/RankTwoTensor.h
  16. +14 −0 framework/src/interfaces/Coupleable.C
  17. +1 −13 framework/src/materials/Material.C
  18. +12 −1 framework/src/materials/MaterialPropertyStorage.C
  19. +3 −3 framework/src/problems/FEProblemBase.C
  20. +0 −30 framework/src/restart/DataIO.C
  21. +6 −6 framework/src/utils/BicubicInterpolation.C
  22. +85 −40 framework/src/utils/ColumnMajorMatrix.C
  23. +273 −0 framework/src/utils/MooseADWrapper.C
  24. +0 −12 framework/src/utils/MooseUtils.C
  25. +216 −88 framework/src/utils/RankFourTensor.C
  26. +325 −225 framework/src/utils/RankTwoTensor.C
  27. +6 −2 modules/phase_field/include/kernels/ACGrGrElasticDrivingForce.h
  28. +3 −1 modules/phase_field/include/kernels/ACInterfaceStress.h
  29. +3 −1 modules/phase_field/include/kernels/AllenCahnElasticEnergyOffDiag.h
  30. +6 −2 modules/phase_field/include/materials/ElasticEnergyMaterial.h
  31. +3 −1 modules/phase_field/include/materials/StrainGradDispDerivatives.h
  32. +0 −1 modules/solid_mechanics/include/kernels/HomogenizationKernel.h
  33. +0 −1 modules/solid_mechanics/include/kernels/StressDivergence.h
  34. +1 −2 modules/solid_mechanics/include/postprocessors/HomogenizedElasticConstants.h
  35. +3 −1 modules/solid_mechanics/include/postprocessors/InteractionIntegralSM.h
  36. +3 −1 modules/tensor_mechanics/include/auxkernels/GlobalDisplacementAux.h
  37. +6 −2 modules/tensor_mechanics/include/kernels/GeneralizedPlaneStrainOffDiag.h
  38. +6 −3 modules/tensor_mechanics/include/kernels/MomentBalancing.h
  39. +0 −1 modules/tensor_mechanics/include/kernels/PhaseFieldFractureMechanicsOffDiag.h
  40. +3 −1 modules/tensor_mechanics/include/kernels/StressDivergenceBeam.h
  41. +0 −2 modules/tensor_mechanics/include/kernels/StressDivergenceTensors.h
  42. +6 −2 modules/tensor_mechanics/include/kernels/WeakPlaneStress.h
  43. +3 −1 modules/tensor_mechanics/include/materials/ComputeEigenstrainBase.h
  44. +3 −1 modules/tensor_mechanics/include/materials/ComputeGlobalStrain.h
  45. +3 −1 modules/tensor_mechanics/include/materials/ComputeInterfaceStress.h
  46. +4 −2 modules/tensor_mechanics/include/materials/ComputeThermalExpansionEigenstrainBase.h
  47. +3 −1 modules/tensor_mechanics/include/materials/EshelbyTensor.h
  48. +6 −2 modules/tensor_mechanics/include/materials/MultiPhaseStressMaterial.h
  49. +3 −1 modules/tensor_mechanics/include/materials/StrainEnergyDensity.h
  50. +3 −1 modules/tensor_mechanics/include/materials/ThermalFractureIntegral.h
  51. +6 −2 modules/tensor_mechanics/include/materials/TwoPhaseStressMaterial.h
  52. +3 −1 modules/tensor_mechanics/include/postprocessors/InteractionIntegral.h
  53. +3 −1 modules/tensor_mechanics/include/postprocessors/JIntegral.h
  54. +6 −2 modules/tensor_mechanics/include/scalarkernels/GlobalStrain.h
  55. +3 −1 modules/tensor_mechanics/include/userobjects/CrackFrontDefinition.h
  56. +6 −2 modules/tensor_mechanics/include/userobjects/GeneralizedPlaneStrainUserObject.h
  57. +10 −1 modules/tensor_mechanics/include/utils/ElasticityTensorTools.h
  58. +8 −6 modules/tensor_mechanics/include/utils/RankTwoScalarTools.h
  59. +2 −2 modules/tensor_mechanics/src/userobjects/CrackFrontDefinition.C
  60. +35 −0 modules/tensor_mechanics/test/include/kernels/ADStressDivergenceTest.h
  61. +49 −0 modules/tensor_mechanics/test/include/materials/TensorMechADMatTest.h
  62. +31 −0 modules/tensor_mechanics/test/src/kernels/ADStressDivergenceTest.C
  63. +59 −0 modules/tensor_mechanics/test/src/materials/TensorMechADMatTest.C
  64. BIN modules/tensor_mechanics/test/tests/ad_simple_linear/gold/linear-out.e
  65. +96 −0 modules/tensor_mechanics/test/tests/ad_simple_linear/linear-ad.i
  66. +82 −0 modules/tensor_mechanics/test/tests/ad_simple_linear/linear-hand-coded.i
  67. +39 −0 modules/tensor_mechanics/test/tests/ad_simple_linear/tests
  68. +0 −2 modules/xfem/include/kernels/CrackTipEnrichmentStressDivergenceTensors.h
  69. +3 −1 modules/xfem/include/userobjects/XFEMRankTwoTensorMarkerUserObject.h
  70. +6 −2 test/include/materials/OutputTestMaterial.h
@@ -13,7 +13,9 @@
#include "Function.h"

class PiecewiseBilinear;
class ColumnMajorMatrix;
template <typename>
class ColumnMajorMatrixTempl;
typedef ColumnMajorMatrixTempl<Real> ColumnMajorMatrix;
class BilinearInterpolation;

template <>
@@ -26,8 +26,10 @@ template <typename T>
class DenseVector;
}

#define adCoupledValue(Name) this->template adCoupledValueTemplate<compute_stage>(Name)
#define adCoupledGradient(Name) this->template adCoupledGradientTemplate<compute_stage>(Name)
#define adCoupledValue this->template adCoupledValueTemplate<compute_stage>
#define adCoupledGradient this->template adCoupledGradientTemplate<compute_stage>
#define adZero this->template adZeroTemplate<compute_stage>
#define adGradZero this->template adGradZeroTemplate<compute_stage>

/**
* Interface for objects that needs coupling capabilities
@@ -596,6 +598,20 @@ class Coupleable
virtual const DenseVector<Number> & coupledSolutionDoFsOlder(const std::string & var_name,
unsigned int comp = 0);

/**
* Template method that returns _zero to RESIDUAL computing objects and _ad_zero to JACOBIAN
* computing objects
*/
template <ComputeStage compute_stage>
const typename VariableValueType<compute_stage>::type & adZeroTemplate();

/**
* Template method that returns _grad_zero to RESIDUAL computing objects and _ad_grad_zero to
* JACOBIAN computing objects
*/
template <ComputeStage compute_stage>
const typename VariableGradientType<compute_stage>::type & adGradZeroTemplate();

protected:
// Reference to the interface's input parameters
const InputParameters & _c_parameters;
@@ -868,4 +884,24 @@ Coupleable::getADDefaultGradient()
template <>
VariableGradient & Coupleable::getADDefaultGradient<RESIDUAL>();

template <ComputeStage compute_stage>
const typename VariableValueType<compute_stage>::type &
Coupleable::adZeroTemplate()
{
return _zero;
}

template <>
const typename VariableValueType<JACOBIAN>::type & Coupleable::adZeroTemplate<JACOBIAN>();

template <ComputeStage compute_stage>
const typename VariableGradientType<compute_stage>::type &
Coupleable::adGradZeroTemplate()
{
return _grad_zero;
}

template <>
const typename VariableGradientType<JACOBIAN>::type & Coupleable::adGradZeroTemplate<JACOBIAN>();

#endif /* COUPLEABLE_H */
@@ -13,7 +13,11 @@
#include "Material.h"
#include "MooseTypes.h"

#define usingMaterialMembers using ADMaterial<compute_stage>::_qp
#define usingMaterialMembers \
using ADMaterial<compute_stage>::_qp; \
using ADMaterial<compute_stage>::_ad_grad_zero; \
using ADMaterial<compute_stage>::_ad_zero; \
using ADMaterial<compute_stage>::coupledComponents

// forward declarations
template <ComputeStage>
@@ -165,7 +165,7 @@ class MaterialData
* Calls resizeProps helper function for regular material properties
*/
template <typename T>
void resizeProps(unsigned int size);
void resizeProps(unsigned int size, bool declared_ad = false);

/// Status of storage swapping (calling swap sets this to true; swapBack sets it to false)
bool _swapped;
@@ -211,7 +211,7 @@ MaterialData::haveProperty(const std::string & prop_name) const

template <typename T>
void
MaterialData::resizeProps(unsigned int size)
MaterialData::resizeProps(unsigned int size, bool declared_ad)
{
auto n = size + 1;
if (_props.size() < n)
@@ -222,7 +222,7 @@ MaterialData::resizeProps(unsigned int size)
_props_older.resize(n, nullptr);

if (_props[size] == nullptr)
_props[size] = new ADMaterialPropertyObject<T>;
_props[size] = new ADMaterialPropertyObject<T>(declared_ad);
if (_props_old[size] == nullptr)
_props_old[size] = new ADMaterialPropertyObject<T>;
if (_props_older[size] == nullptr)
@@ -277,7 +277,7 @@ MaterialData::declareADHelper(MaterialProperties & props,
const std::string & libmesh_dbg_var(prop_name),
unsigned int prop_id)
{
resizeProps<T>(prop_id);
resizeProps<T>(prop_id, true);
auto prop = dynamic_cast<ADMaterialPropertyObject<T> *>(props[prop_id]);
mooseAssert(prop != nullptr, "Internal error in declaring material property: " + prop_name);
return *prop;
@@ -27,7 +27,7 @@ class PropertyValue;
* types
*/
template <typename P>
PropertyValue * _init_helper(int size, PropertyValue * prop, const P * the_type);
PropertyValue * _init_helper(int size, PropertyValue * prop, const P * the_type, bool use_ad);

/**
* Abstract definition of a property value.
@@ -70,11 +70,8 @@ class PropertyValue
virtual void store(std::ostream & stream) = 0;
virtual void load(std::istream & stream) = 0;

/**
* copy the Real version to the DualNumber<Real> version of the material property for the
* specified quadrature point
*/
virtual void copyValueToDualNumber(const unsigned int i) = 0;
/// Mark whether this property is in AD mode. This method is necessary for switching the state after swapping material properties during stateful material calculations
virtual void markAD(bool use_ad) = 0;

/**
* copy the value portion (not the derivatives) of the DualNumber<Real> version of the material
@@ -106,11 +103,11 @@ class MaterialProperty : public PropertyValue
{
public:
/// Explicitly declare a public constructor because we made the copy constructor private
MaterialProperty() : PropertyValue()
MaterialProperty(bool use_ad = false) : PropertyValue(), _use_ad(use_ad)
{ /* */
}

virtual ~MaterialProperty() { _value.release(); }
virtual ~MaterialProperty() {}

/**
* @returns a read-only reference to the parameter value.
@@ -174,28 +171,14 @@ class MaterialProperty : public PropertyValue
*/
virtual void load(std::istream & stream) override;

/**
* Friend helper function to handle scalar material property initializations
* @param size - the size corresponding to the quadrature rule
* @param prop - The Material property that we will resize since this is not a member
* @param the_type - This is just a template parameter used to identify the
* difference between the scalar and vector template functions
*/
template <typename P>
friend PropertyValue * _init_helper(int size, PropertyValue * prop, const P * the_type);

/**
* copy the Real version to the DualNumber<Real> version of the material property for the
* specified quadrature point
*/
void copyValueToDualNumber(const unsigned int i) override { _value[i].copyValueToDualNumber(); }

/**
* copy the value portion (not the derivatives) of the DualNumber<Real> version of the material
* property to the Real version for the specified quadrature point
*/
void copyDualNumberToValue(const unsigned int i) override { _value[i].copyDualNumberToValue(); }

void markAD(bool use_ad) override;

private:
/// private copy constructor to avoid shallow copying of material properties
MaterialProperty(const MaterialProperty<T> & /*src*/)
@@ -210,12 +193,24 @@ class MaterialProperty : public PropertyValue
}

protected:
/// Whether this property was declared as AD
bool _use_ad;
/// Stored parameter value.
MooseArray<MooseADWrapper<T>> _value;
std::vector<MooseADWrapper<T>> _value;
};

// ------------------------------------------------------------
// Material::Property<> class inline methods

template <typename T>
inline void
MaterialProperty<T>::markAD(bool use_ad)
{
_use_ad = use_ad;
for (auto && moose_wrapper : _value)
moose_wrapper.markAD(use_ad);
}

template <typename T>
inline std::string
MaterialProperty<T>::type()
@@ -227,14 +222,19 @@ template <typename T>
inline PropertyValue *
MaterialProperty<T>::init(int size)
{
return _init_helper(size, this, static_cast<T *>(0));
return _init_helper(size, this, static_cast<T *>(0), _use_ad);
}

template <typename T>
inline void
MaterialProperty<T>::resize(int n)
{
_value.resize(n);
auto diff = n - static_cast<int>(_value.size());
if (diff < 0)
_value.erase(_value.end() + diff, _value.end());
else
for (decltype(diff) i = 0; i < diff; ++i)
_value.emplace_back(MooseADWrapper<T>(_use_ad));
}

template <typename T>
@@ -275,7 +275,7 @@ template <typename T>
class ADMaterialPropertyObject : public MaterialProperty<T>
{
public:
ADMaterialPropertyObject() : MaterialProperty<T>() {}
ADMaterialPropertyObject(bool use_ad = false) : MaterialProperty<T>(use_ad) {}

/**
* Get element i out of the array as a writeable reference.
@@ -358,10 +358,10 @@ dataLoad(std::istream & stream, MaterialProperties & v, void * context)
// Scalar Init Helper Function
template <typename P>
PropertyValue *
_init_helper(int size, PropertyValue * /*prop*/, const P *)
_init_helper(int size, PropertyValue * /*prop*/, const P *, bool use_ad)
{
MaterialProperty<P> * copy = new MaterialProperty<P>;
copy->_value.resize(size, MooseADWrapper<P>{});
ADMaterialPropertyObject<P> * copy = new ADMaterialPropertyObject<P>(use_ad);
copy->resize(size);
return copy;
}

@@ -125,7 +125,7 @@ class MaterialPropertyStorage
* properties are
* reused for computing current properties. This is called when solve succeeded.
*/
void shift();
void shift(const FEProblemBase & fe_problem);

/**
* Copy material properties from elem_from to elem_to. Thread safe.
@@ -16,6 +16,8 @@
#include "HashMap.h"
#include "MooseError.h"
#include "Backup.h"
#include "RankTwoTensor.h"
#include "RankFourTensor.h"

#include "libmesh/vector_value.h"
#include "libmesh/tensor_value.h"
@@ -36,9 +38,9 @@
#include <memory>

// Forward declarations
class ColumnMajorMatrix;
class RankTwoTensor;
class RankFourTensor;
template <typename>
class ColumnMajorMatrixTempl;
typedef ColumnMajorMatrixTempl<Real> ColumnMajorMatrix;
namespace libMesh
{
template <typename T>
@@ -334,10 +336,6 @@ template <>
void dataStore(std::ostream & stream, std::stringstream & s, void * context);
template <>
void dataStore(std::ostream & stream, std::stringstream *& s, void * context);
template <>
void dataStore(std::ostream & stream, RankTwoTensor &, void *);
template <>
void dataStore(std::ostream & stream, RankFourTensor &, void *);

inline void
dataStore(std::ostream & stream, ADReal & dn, void * context)
@@ -347,6 +345,14 @@ dataStore(std::ostream & stream, ADReal & dn, void * context)
dataStore(stream, dn.derivatives()[i], context);
}

template <std::size_t N>
inline void
dataStore(std::ostream & stream, ADReal (&dn)[N], void * context)
{
for (std::size_t i = 0; i < N; ++i)
dataStore(stream, dn[i], context);
}

template <typename T>
void
dataStore(std::ostream & stream, NumericVector<T> & v, void * context)
@@ -403,12 +409,36 @@ dataStore(std::ostream & stream, VectorValue<T> & v, void * context)
}
}

template <typename T>
void
dataStore(std::ostream & stream, RankTwoTensorTempl<T> & rtt, void * context)
{
dataStore(stream, rtt._coords, context);
}

template <typename T>
void
dataStore(std::ostream & stream, RankFourTensorTempl<T> & rft, void * context)
{
dataStore(stream, rft._vals, context);
}

template <typename T>
inline void
dataStore(std::ostream & stream, MooseADWrapper<T> & dn_wrapper, void * context)
{
dataStore(stream, dn_wrapper.value(), context);
dataStore(stream, dn_wrapper.dn(false), context);
if (dn_wrapper._dual_number)
{
unsigned int m = 1;
stream.write((char *)&m, sizeof(m));
dataStore(stream, *dn_wrapper._dual_number, context);
}
else
{
unsigned int m = 0;
stream.write((char *)&m, sizeof(m));
}
}

// DO NOT MODIFY THE NEXT LINE - It is used by MOOSEDocs
@@ -578,10 +608,6 @@ template <>
void dataLoad(std::istream & stream, std::stringstream & s, void * context);
template <>
void dataLoad(std::istream & stream, std::stringstream *& s, void * context);
template <>
void dataLoad(std::istream & stream, RankTwoTensor &, void *);
template <>
void dataLoad(std::istream & stream, RankFourTensor &, void *);

inline void
dataLoad(std::istream & stream, ADReal & dn, void * context)
@@ -653,12 +679,32 @@ dataLoad(std::istream & stream, VectorValue<T> & v, void * context)
}
}

template <typename T>
void
dataLoad(std::istream & stream, RankTwoTensorTempl<T> & rtt, void * context)
{
dataLoad(stream, rtt._coords, context);
}

template <typename T>
void
dataLoad(std::istream & stream, RankFourTensorTempl<T> & rft, void * context)
{
dataLoad(stream, rft._vals, context);
}

template <typename T>
inline void
dataLoad(std::istream & stream, MooseADWrapper<T> & dn_wrapper, void * context)
{
dataLoad(stream, dn_wrapper.value(), context);
dataLoad(stream, dn_wrapper.dn(false), context);
unsigned int n = 0;
stream.read((char *)&n, sizeof(n));
if (n)
{
dn_wrapper._dual_number = libmesh_make_unique<typename MooseADWrapper<T>::DNType>();
dataLoad(stream, *dn_wrapper._dual_number, context);
}
}

// Scalar Helper Function
Oops, something went wrong.

0 comments on commit 743232c

Please sign in to comment.