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