Line data Source code
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 "libmesh/int_range.h"
14 : #include "metaphysicl/raw_type.h"
15 : #include "metaphysicl/metaphysicl_version.h"
16 :
17 : namespace Eigen
18 : {
19 : namespace internal
20 : {
21 : template <typename V, typename D, bool asd>
22 : inline bool
23 336 : isinf_impl(const MetaPhysicL::DualNumber<V, D, asd> & a)
24 : {
25 : using std::isinf;
26 336 : return isinf(a);
27 : }
28 :
29 : template <typename V, typename D, bool asd>
30 : inline bool
31 336 : isnan_impl(const MetaPhysicL::DualNumber<V, D, asd> & a)
32 : {
33 : using std::isnan;
34 336 : return isnan(a);
35 : }
36 :
37 : template <typename V, typename D, bool asd>
38 : inline MetaPhysicL::DualNumber<V, D, asd>
39 : sqrt(const MetaPhysicL::DualNumber<V, D, asd> & a)
40 : {
41 : #if METAPHYSICL_MAJOR_VERSION < 2
42 : return std::sqrt(a);
43 : #else
44 : return MetaPhysicL::sqrt(a);
45 : #endif
46 : }
47 :
48 : template <typename V, typename D, bool asd>
49 : inline MetaPhysicL::DualNumber<V, D, asd>
50 : abs(const MetaPhysicL::DualNumber<V, D, asd> & a)
51 : {
52 : #if METAPHYSICL_MAJOR_VERSION < 2
53 : return std::abs(a);
54 : #else
55 : return MetaPhysicL::abs(a);
56 : #endif
57 : }
58 : }
59 : } // namespace Eigen
60 :
61 : // this include _must_ come after the Eigen::internal overloads above. We also ignore a warning
62 : // about an Eigen internal use of a potentially uninitialized variable
63 : #include "libmesh/ignore_warnings.h"
64 : #include <Eigen/Core>
65 : #include "libmesh/restore_warnings.h"
66 :
67 : // Eigen needs this
68 : namespace MetaPhysicL
69 : {
70 : // raw_value AD->non-AD conversion for ADReal valued Eigen::Matrix objects
71 : template <typename T, int M, int N, int O, int M2, int N2>
72 : struct RawType<Eigen::Matrix<T, M, N, O, M2, N2>>
73 : {
74 : typedef Eigen::Matrix<typename RawType<T>::value_type, M, N, O, M2, N2> value_type;
75 :
76 4 : static value_type value(const Eigen::Matrix<T, M, N, O, M2, N2> & in)
77 : {
78 40 : return value_type::NullaryExpr([&in](Eigen::Index i) { return raw_value(in(i)); });
79 : }
80 : };
81 :
82 : // specialized for RealEigenVector
83 : template <typename T, int Options, int MaxSize>
84 : struct RawType<Eigen::Matrix<T, -1, 1, Options, MaxSize, 1>>
85 : {
86 : typedef Eigen::Matrix<typename RawType<T>::value_type, -1, 1, Options, MaxSize, 1> value_type;
87 :
88 525444 : static value_type value(const Eigen::Matrix<T, -1, 1, Options, MaxSize, 1> & in)
89 : {
90 525444 : value_type ret(in.size());
91 1576332 : for (const auto i : libMesh::make_range(in.size()))
92 1050888 : ret(i) = raw_value(in(i));
93 525444 : return ret;
94 0 : }
95 : };
96 :
97 : // raw_value overload for Map type objects that forces evaluation
98 : template <typename T>
99 : auto
100 : raw_value(const Eigen::Map<T> & in)
101 : {
102 : return raw_value(in.eval());
103 : }
104 : } // namespace MetaPhysicL
105 :
106 : namespace Eigen
107 : {
108 : // libEigen support for dual number types
109 : template <typename V, typename D, bool asd>
110 : struct NumTraits<MetaPhysicL::DualNumber<V, D, asd>>
111 : : NumTraits<V> // permits to get the epsilon, dummy_precision, lowest, highest functions
112 : {
113 : typedef MetaPhysicL::DualNumber<V, D, asd> Real;
114 : typedef MetaPhysicL::DualNumber<V, D, asd> NonInteger;
115 : typedef MetaPhysicL::DualNumber<V, D, asd> Nested;
116 :
117 : enum
118 : {
119 : IsComplex = 0,
120 : IsInteger = 0,
121 : IsSigned = 1,
122 : RequireInitialization = 1,
123 : ReadCost = HugeCost,
124 : AddCost = HugeCost,
125 : MulCost = HugeCost
126 : };
127 : };
128 :
129 : template <typename BinaryOp, typename V, typename D, bool asd>
130 : struct ScalarBinaryOpTraits<Real, MetaPhysicL::DualNumber<V, D, asd>, BinaryOp>
131 : {
132 : typedef MetaPhysicL::DualNumber<V, D, asd> ReturnType;
133 : };
134 : template <typename BinaryOp, typename V, typename D, bool asd>
135 : struct ScalarBinaryOpTraits<MetaPhysicL::DualNumber<V, D, asd>, Real, BinaryOp>
136 : {
137 : typedef MetaPhysicL::DualNumber<V, D, asd> ReturnType;
138 : };
139 : } // namespace Eigen
140 :
141 : namespace Moose
142 : {
143 : template <typename T>
144 : struct ADType;
145 :
146 : template <typename T, int M, int N, int O, int M2, int N2>
147 : struct ADType<Eigen::Matrix<T, M, N, O, M2, N2>>
148 : {
149 : typedef typename Eigen::Matrix<typename ADType<T>::type, M, N, O, M2, N2> type;
150 : };
151 : }
152 :
153 : namespace Eigen::internal
154 : {
155 :
156 : template <>
157 : inline long double
158 : cast<ADReal, long double>(const ADReal & x)
159 : {
160 : return MetaPhysicL::raw_value(x);
161 : }
162 :
163 : template <>
164 : inline double
165 : cast<ADReal, double>(const ADReal & x)
166 : {
167 : return MetaPhysicL::raw_value(x);
168 : }
169 :
170 : template <>
171 : inline long
172 : cast<ADReal, long>(const ADReal & x)
173 : {
174 : return MetaPhysicL::raw_value(x);
175 : }
176 :
177 : template <>
178 : inline int
179 : cast<ADReal, int>(const ADReal & x)
180 : {
181 : return MetaPhysicL::raw_value(x);
182 : }
183 :
184 : /**
185 : * Override number traits for the ADReal type used in libEigen calculations.
186 : * nr and mr are used to determine block matrix (panel) sizes. We keep them
187 : * as small as possible to avoid panel sizes that exceed the allowable stack
188 : * size. Those numbers could be made a function of the Eigen stack limit and
189 : * the number of derivative entries.
190 : */
191 : template <>
192 : class gebp_traits<ADReal, ADReal, false, false>
193 : {
194 : public:
195 : typedef ADReal ResScalar;
196 : enum
197 : {
198 : Vectorizable = false,
199 : LhsPacketSize = 1,
200 : RhsPacketSize = 1,
201 : ResPacketSize = 1,
202 : NumberOfRegisters = 1,
203 : nr = 1,
204 : mr = 1,
205 : LhsProgress = 1,
206 : RhsProgress = 1
207 : };
208 : typedef ResScalar LhsPacket;
209 : typedef ResScalar RhsPacket;
210 : typedef ResScalar ResPacket;
211 : typedef ResScalar AccPacket;
212 : typedef ResScalar LhsPacket4Packing;
213 : };
214 :
215 : template <typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs>
216 : struct gebp_kernel<ADReal, ADReal, Index, DataMapper, 1, 1, ConjugateLhs, ConjugateRhs>
217 : {
218 : EIGEN_DONT_INLINE
219 54 : void operator()(const DataMapper & res,
220 : const ADReal * blockA,
221 : const ADReal * blockB,
222 : Index rows,
223 : Index depth,
224 : Index cols,
225 : const ADReal & alpha,
226 : Index strideA = -1,
227 : Index strideB = -1,
228 : Index offsetA = 0,
229 : Index offsetB = 0)
230 : {
231 54 : if (rows == 0 || cols == 0 || depth == 0)
232 0 : return;
233 :
234 54 : ADReal acc1;
235 :
236 54 : if (strideA == -1)
237 2 : strideA = depth;
238 54 : if (strideB == -1)
239 2 : strideB = depth;
240 :
241 468 : for (Index i = 0; i < rows; ++i)
242 840 : for (Index j = 0; j < cols; ++j)
243 : {
244 426 : const ADReal * A = blockA + i * strideA + offsetA;
245 426 : const ADReal * B = blockB + j * strideB + offsetB;
246 :
247 426 : acc1 = 0;
248 2100 : for (Index k = 0; k < depth; k++)
249 1674 : acc1 += A[k] * B[k];
250 834 : res(i, j) += acc1 * alpha;
251 : }
252 54 : }
253 : };
254 :
255 : }
|
