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 :
15 : namespace Eigen
16 : {
17 : namespace internal
18 : {
19 : template <typename V, typename D, bool asd>
20 : inline bool
21 168 : isinf_impl(const MetaPhysicL::DualNumber<V, D, asd> & a)
22 : {
23 168 : return std::isinf(a);
24 : }
25 :
26 : template <typename V, typename D, bool asd>
27 : inline bool
28 168 : isnan_impl(const MetaPhysicL::DualNumber<V, D, asd> & a)
29 : {
30 168 : return std::isnan(a);
31 : }
32 :
33 : template <typename V, typename D, bool asd>
34 : inline MetaPhysicL::DualNumber<V, D, asd>
35 : sqrt(const MetaPhysicL::DualNumber<V, D, asd> & a)
36 : {
37 : return std::sqrt(a);
38 : }
39 :
40 : template <typename V, typename D, bool asd>
41 : inline MetaPhysicL::DualNumber<V, D, asd>
42 : abs(const MetaPhysicL::DualNumber<V, D, asd> & a)
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 1 : static value_type value(const Eigen::Matrix<T, M, N, O, M2, N2> & in)
65 : {
66 11 : 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 : {
86 : typedef MetaPhysicL::DualNumber<V, D, asd> Real;
87 : typedef MetaPhysicL::DualNumber<V, D, asd> NonInteger;
88 : typedef MetaPhysicL::DualNumber<V, D, asd> Nested;
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 : {
105 : typedef MetaPhysicL::DualNumber<V, D, asd> ReturnType;
106 : };
107 : template <typename BinaryOp, typename V, typename D, bool asd>
108 : struct ScalarBinaryOpTraits<MetaPhysicL::DualNumber<V, D, asd>, Real, BinaryOp>
109 : {
110 : typedef MetaPhysicL::DualNumber<V, D, asd> ReturnType;
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
131 : cast<ADReal, long double>(const ADReal & x)
132 : {
133 : return MetaPhysicL::raw_value(x);
134 : }
135 :
136 : template <>
137 : inline double
138 : cast<ADReal, double>(const ADReal & x)
139 : {
140 : return MetaPhysicL::raw_value(x);
141 : }
142 :
143 : template <>
144 : inline long
145 : cast<ADReal, long>(const ADReal & x)
146 : {
147 : return MetaPhysicL::raw_value(x);
148 : }
149 :
150 : template <>
151 : inline int
152 : cast<ADReal, int>(const ADReal & x)
153 : {
154 : return MetaPhysicL::raw_value(x);
155 : }
156 :
157 : /**
158 : * Override number traits for the ADReal type used in libEigen calculations.
159 : * nr and mr are used to determine block matrix (panel) sizes. We keep them
160 : * as small as possible to avoid panel sizes that exceed the allowable stack
161 : * size. Those numbers could be made a function of the Eigen stack limit and
162 : * the number of derivative entries.
163 : */
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 : };
181 : typedef ResScalar LhsPacket;
182 : typedef ResScalar RhsPacket;
183 : typedef ResScalar ResPacket;
184 : typedef ResScalar AccPacket;
185 : typedef ResScalar LhsPacket4Packing;
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 27 : 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 27 : if (rows == 0 || cols == 0 || depth == 0)
205 0 : return;
206 :
207 27 : ADReal acc1;
208 :
209 27 : if (strideA == -1)
210 1 : strideA = depth;
211 27 : if (strideB == -1)
212 1 : strideB = depth;
213 :
214 234 : for (Index i = 0; i < rows; ++i)
215 420 : for (Index j = 0; j < cols; ++j)
216 : {
217 213 : const ADReal * A = blockA + i * strideA + offsetA;
218 213 : const ADReal * B = blockB + j * strideB + offsetB;
219 :
220 213 : acc1 = 0;
221 1050 : for (Index k = 0; k < depth; k++)
222 837 : acc1 += A[k] * B[k];
223 417 : res(i, j) += acc1 * alpha;
224 : }
225 27 : }
226 : };
227 :
228 : }
|
