Merge pull request #13441 from dschwen/adlinearinterpolation_13432

Enable LinearInterpolation for AD

framework/include/utils/LinearInterpolation.h

@@ -9,25 +9,28 @@
 
 #pragma once
 
+#include "Moose.h"
+#include "MooseTypes.h"
+#include "DualReal.h"
+
 #include <vector>
 #include <string>
 
-#include "Moose.h"
-
 /**
  * This class interpolates values given a set of data pairs and an abscissa.
  */
-class LinearInterpolation
+template <typename T>
+class LinearInterpolationTempl
 {
 public:
   /* Constructor, Takes two vectors of points for which to apply the fit.  One should be of the
    * independent variable while the other should be of the dependent variable.  These values should
    * correspond to one and other in the same position.
    */
-  LinearInterpolation(const std::vector<Real> & X, const std::vector<Real> & Y);
-  LinearInterpolation() : _x(std::vector<Real>()), _y(std::vector<Real>()) {}
+  LinearInterpolationTempl(const std::vector<Real> & X, const std::vector<Real> & Y);
+  LinearInterpolationTempl() : _x(std::vector<Real>()), _y(std::vector<Real>()) {}
 
-  virtual ~LinearInterpolation() = default;
+  virtual ~LinearInterpolationTempl() = default;
 
   /**
    * Set the x and y values.
@@ -45,14 +48,14 @@ class LinearInterpolation
    * This function will take an independent variable input and will return the dependent variable
    * based on the generated fit
    */
-  Real sample(Real x) const;
+  T sample(const T & x) const;
 
   /**
    * This function will take an independent variable input and will return the derivative of the
    * dependent variable
    * with respect to the independent variable based on the generated fit
    */
-  Real sampleDerivative(Real x) const;
+  T sampleDerivative(const T & x) const;
 
   /**
    * This function returns the size of the array holding the points, i.e. the number of sample
@@ -75,3 +78,18 @@ class LinearInterpolation
   static int _file_number;
 };
 
+#define ADLinearInterpolation typename LinearInterpolationType<compute_stage>::type
+
+typedef LinearInterpolationTempl<Real> LinearInterpolation;
+typedef LinearInterpolationTempl<DualReal> DualLinearInterpolation;
+
+template <ComputeStage compute_stage>
+struct LinearInterpolationType
+{
+  typedef LinearInterpolation type;
+};
+template <>
+struct LinearInterpolationType<JACOBIAN>
+{
+  typedef DualLinearInterpolation type;
+};

framework/src/utils/LinearInterpolation.C

@@ -9,20 +9,27 @@
 
 #include "LinearInterpolation.h"
 
+#include "metaphysicl/numberarray.h"
+#include "metaphysicl/dualnumber.h"
+
 #include <cassert>
 #include <fstream>
 #include <stdexcept>
 
-int LinearInterpolation::_file_number = 0;
+template <typename T>
+int LinearInterpolationTempl<T>::_file_number = 0;
 
-LinearInterpolation::LinearInterpolation(const std::vector<Real> & x, const std::vector<Real> & y)
+template <typename T>
+LinearInterpolationTempl<T>::LinearInterpolationTempl(const std::vector<Real> & x,
+                                                      const std::vector<Real> & y)
   : _x(x), _y(y)
 {
   errorCheck();
 }
 
+template <typename T>
 void
-LinearInterpolation::errorCheck()
+LinearInterpolationTempl<T>::errorCheck()
 {
   if (_x.size() != _y.size())
     throw std::domain_error("Vectors are not the same length");
@@ -37,10 +44,11 @@ LinearInterpolation::errorCheck()
     }
 }
 
-Real
-LinearInterpolation::sample(Real x) const
+template <typename T>
+T
+LinearInterpolationTempl<T>::sample(const T & x) const
 {
-  // sanity check (empty LinearInterpolations get constructed in many places
+  // sanity check (empty LinearInterpolationTempls get constructed in many places
   // so we cannot put this into the errorCheck)
   assert(_x.size() > 0);
 
@@ -58,8 +66,9 @@ LinearInterpolation::sample(Real x) const
   return 0;
 }
 
-Real
-LinearInterpolation::sampleDerivative(Real x) const
+template <typename T>
+T
+LinearInterpolationTempl<T>::sampleDerivative(const T & x) const
 {
   // endpoint cases
   if (x < _x[0])
@@ -75,8 +84,9 @@ LinearInterpolation::sampleDerivative(Real x) const
   return 0;
 }
 
+template <typename T>
 Real
-LinearInterpolation::integrate()
+LinearInterpolationTempl<T>::integrate()
 {
   Real answer = 0;
   for (unsigned int i = 1; i < _x.size(); ++i)
@@ -85,20 +95,26 @@ LinearInterpolation::integrate()
   return answer;
 }
 
+template <typename T>
 Real
-LinearInterpolation::domain(int i) const
+LinearInterpolationTempl<T>::domain(int i) const
 {
   return _x[i];
 }
 
+template <typename T>
 Real
-LinearInterpolation::range(int i) const
+LinearInterpolationTempl<T>::range(int i) const
 {
   return _y[i];
 }
 
+template <typename T>
 unsigned int
-LinearInterpolation::getSampleSize()
+LinearInterpolationTempl<T>::getSampleSize()
 {
   return _x.size();
 }
+
+template class LinearInterpolationTempl<Real>;
+template class LinearInterpolationTempl<DualReal>;

modules/solid_mechanics/include/materials/IsotropicTempDepHardening.h

@@ -10,9 +10,9 @@
 #pragma once
 
 #include "IsotropicPlasticity.h"
+#include "LinearInterpolation.h"
 
 class PiecewiseLinear;
-class LinearInterpolation;
 
 class IsotropicTempDepHardening;
 
@@ -42,4 +42,3 @@ class IsotropicTempDepHardening : public IsotropicPlasticity
   unsigned int _hf_index_hi;
   Real _hf_fraction;
 };
-

modules/tensor_mechanics/include/materials/TemperatureDependentHardeningStressUpdate.h

@@ -10,9 +10,9 @@
 #pragma once
 
 #include "IsotropicPlasticityStressUpdate.h"
+#include "LinearInterpolation.h"
 
 class PiecewiseLinear;
-class LinearInterpolation;
 
 class TemperatureDependentHardeningStressUpdate;
 
@@ -67,4 +67,3 @@ class TemperatureDependentHardeningStressUpdate : public IsotropicPlasticityStre
    */
   Real _hf_fraction;
 };
-

test/include/materials/LinearInterpolationMaterial.h

@@ -10,9 +10,9 @@
 #pragma once
 
 #include "Material.h"
+#include "LinearInterpolation.h"
 
 class LinearInterpolationMaterial;
-class LinearInterpolation;
 class PolynomialFit;
 
 template <>
@@ -31,4 +31,3 @@ class LinearInterpolationMaterial : public Material
   std::unique_ptr<PolynomialFit> _poly_fit;
   MaterialProperty<Real> & _property;
 };
-

unit/src/LinearInterpolationTest.C

@@ -10,21 +10,25 @@
 #include "gtest/gtest.h"
 
 #include "LinearInterpolation.h"
+#include "DualReal.h"
+
+#include "metaphysicl/numberarray.h"
+#include "metaphysicl/dualnumber.h"
 
 #include <cmath>
 
 TEST(LinearInterpolationTest, getSampleSize)
 {
-  std::vector<double> x = {1, 2, 3, 5};
-  std::vector<double> y = {0, 5, 6, 8};
+  std::vector<Real> x = {1, 2, 3, 5};
+  std::vector<Real> y = {0, 5, 6, 8};
   LinearInterpolation interp(x, y);
   ASSERT_EQ(interp.getSampleSize(), x.size());
 }
 
 TEST(LinearInterpolationTest, sample)
 {
-  std::vector<double> x = {1, 2, 3, 5};
-  std::vector<double> y = {0, 5, 6, 8};
+  std::vector<Real> x = {1, 2, 3, 5};
+  std::vector<Real> y = {0, 5, 6, 8};
   LinearInterpolation interp(x, y);
 
   EXPECT_DOUBLE_EQ(interp.sample(0.), 0.);
@@ -42,3 +46,17 @@ TEST(LinearInterpolationTest, sample)
   EXPECT_DOUBLE_EQ(interp.sampleDerivative(2.), 1.);
   EXPECT_DOUBLE_EQ(interp.sampleDerivative(2.1), 1.);
 }
+
+TEST(LinearInterpolationTest, automatic_differentiation_sample)
+{
+  std::vector<Real> x = {1, 2};
+  std::vector<Real> y = {0, 5};
+  DualLinearInterpolation interp(x, y);
+
+  DualReal xx = 1.5;
+  xx.derivatives()[0] = 1;
+  auto yy = interp.sample(xx);
+
+  EXPECT_DOUBLE_EQ(yy.value(), 2.5);
+  EXPECT_DOUBLE_EQ(yy.derivatives()[0], 5.0);
+}