https://mooseframework.inl.gov
EigenADReal.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #pragma once
11 
12 #include "ADReal.h"
13 #include "metaphysicl/raw_type.h"
14 
15 namespace Eigen
16 {
17 namespace internal
18 {
19 template <typename V, typename D, bool asd>
20 inline bool
22 {
23  return std::isinf(a);
24 }
25 
26 template <typename V, typename D, bool asd>
27 inline bool
29 {
30  return std::isnan(a);
31 }
32 
33 template <typename V, typename D, bool asd>
36 {
37  return std::sqrt(a);
38 }
39 
40 template <typename V, typename D, bool asd>
43 {
44  return std::abs(a);
45 }
46 }
47 } // namespace Eigen
48 
49 // this include _must_ come after the Eigen::internal overloads above. We also ignore a warning
50 // about an Eigen internal use of a potentially uninitialized variable
51 #include "libmesh/ignore_warnings.h"
52 #include <Eigen/Core>
53 #include "libmesh/restore_warnings.h"
54 
55 // Eigen needs this
56 namespace MetaPhysicL
57 {
58 // raw_value AD->non-AD conversion for ADReal valued Eigen::Matrix objects
59 template <typename T, int M, int N, int O, int M2, int N2>
60 struct RawType<Eigen::Matrix<T, M, N, O, M2, N2>>
61 {
62  typedef Eigen::Matrix<typename RawType<T>::value_type, M, N, O, M2, N2> value_type;
63 
64  static value_type value(const Eigen::Matrix<T, M, N, O, M2, N2> & in)
65  {
66  return value_type::NullaryExpr([&in](Eigen::Index i) { return raw_value(in(i)); });
67  }
68 };
69 
70 // raw_value overload for Map type objects that forces evaluation
71 template <typename T>
72 auto
73 raw_value(const Eigen::Map<T> & in)
74 {
75  return raw_value(in.eval());
76 }
77 } // namespace MetaPhysicL
78 
79 namespace Eigen
80 {
81 // libEigen support for dual number types
82 template <typename V, typename D, bool asd>
83 struct NumTraits<MetaPhysicL::DualNumber<V, D, asd>>
84  : NumTraits<V> // permits to get the epsilon, dummy_precision, lowest, highest functions
85 {
89 
90  enum
91  {
92  IsComplex = 0,
93  IsInteger = 0,
94  IsSigned = 1,
95  RequireInitialization = 1,
96  ReadCost = HugeCost,
97  AddCost = HugeCost,
98  MulCost = HugeCost
99  };
100 };
101 
102 template <typename BinaryOp, typename V, typename D, bool asd>
103 struct ScalarBinaryOpTraits<Real, MetaPhysicL::DualNumber<V, D, asd>, BinaryOp>
104 {
106 };
107 template <typename BinaryOp, typename V, typename D, bool asd>
108 struct ScalarBinaryOpTraits<MetaPhysicL::DualNumber<V, D, asd>, Real, BinaryOp>
109 {
111 };
112 } // namespace Eigen
113 
114 namespace Moose
115 {
116 template <typename T>
117 struct ADType;
118 
119 template <typename T, int M, int N, int O, int M2, int N2>
120 struct ADType<Eigen::Matrix<T, M, N, O, M2, N2>>
121 {
122  typedef typename Eigen::Matrix<typename ADType<T>::type, M, N, O, M2, N2> type;
123 };
124 }
125 
126 namespace Eigen::internal
127 {
128 
129 template <>
130 inline long double
132 {
133  return MetaPhysicL::raw_value(x);
134 }
135 
136 template <>
137 inline double
139 {
140  return MetaPhysicL::raw_value(x);
141 }
142 
143 template <>
144 inline long
146 {
147  return MetaPhysicL::raw_value(x);
148 }
149 
150 template <>
151 inline int
153 {
154  return MetaPhysicL::raw_value(x);
155 }
156 
164 template <>
165 class gebp_traits<ADReal, ADReal, false, false>
166 {
167 public:
168  typedef ADReal ResScalar;
169  enum
170  {
171  Vectorizable = false,
172  LhsPacketSize = 1,
173  RhsPacketSize = 1,
174  ResPacketSize = 1,
175  NumberOfRegisters = 1,
176  nr = 1,
177  mr = 1,
178  LhsProgress = 1,
179  RhsProgress = 1
180  };
186 };
187 
188 template <typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs>
189 struct gebp_kernel<ADReal, ADReal, Index, DataMapper, 1, 1, ConjugateLhs, ConjugateRhs>
190 {
191  EIGEN_DONT_INLINE
192  void operator()(const DataMapper & res,
193  const ADReal * blockA,
194  const ADReal * blockB,
195  Index rows,
196  Index depth,
197  Index cols,
198  const ADReal & alpha,
199  Index strideA = -1,
200  Index strideB = -1,
201  Index offsetA = 0,
202  Index offsetB = 0)
203  {
204  if (rows == 0 || cols == 0 || depth == 0)
205  return;
206 
207  ADReal acc1;
208 
209  if (strideA == -1)
210  strideA = depth;
211  if (strideB == -1)
212  strideB = depth;
213 
214  for (Index i = 0; i < rows; ++i)
215  for (Index j = 0; j < cols; ++j)
216  {
217  const ADReal * A = blockA + i * strideA + offsetA;
218  const ADReal * B = blockB + j * strideB + offsetB;
219 
220  acc1 = 0;
221  for (Index k = 0; k < depth; k++)
222  acc1 += A[k] * B[k];
223  res(i, j) += acc1 * alpha;
224  }
225  }
226 };
227 
228 }
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:42
Eigen::Matrix< typename RawType< T >::value_type, M, N, O, M2, N2 > value_type
Definition: EigenADReal.h:62
MetaPhysicL::DualNumber< V, D, asd > Nested
Definition: EigenADReal.h:88
int cast< ADReal, int >(const ADReal &x)
Definition: EigenADReal.h:152
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
We need to instantiate the following CompareTypes to tell the compiler that ADReal is a subtype of Ch...
DualNumber< Real, DNDerivativeType, true > ADReal
Definition: ADRealForward.h:46
double cast< ADReal, double >(const ADReal &x)
Definition: EigenADReal.h:138
static value_type value(const Eigen::Matrix< T, M, N, O, M2, N2 > &in)
Definition: EigenADReal.h:64
MetaPhysicL::DualNumber< V, D, asd > sqrt(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:35
EIGEN_DONT_INLINE void operator()(const DataMapper &res, const ADReal *blockA, const ADReal *blockB, Index rows, Index depth, Index cols, const ADReal &alpha, Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
Definition: EigenADReal.h:192
long double cast< ADReal, long double >(const ADReal &x)
Definition: EigenADReal.h:131
bool isinf_impl(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:21
MetaPhysicL::DualNumber< V, D, asd > Real
Definition: EigenADReal.h:86
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool isnan_impl(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:28
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
long cast< ADReal, long >(const ADReal &x)
Definition: EigenADReal.h:145
MetaPhysicL::DualNumber< V, D, asd > NonInteger
Definition: EigenADReal.h:87
Eigen::Matrix< typename ADType< T >::type, M, N, O, M2, N2 > type
Definition: EigenADReal.h:122