Nalu
Nalu: a generalized unstructured massively parallel low Mach flow code designed to support a variety of energy applications of interest (most notably Wind ECP) built on the Sierra Toolkit and Trilinos solver Tpetra/Epetra stack. The open source BSD, clause 3 license model has been chosen for the code base. See LICENSE for more information. http://NaluCFD.org
MasterElementFunctions.h
Go to the documentation of this file.
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2014 Sandia Corporation. */
3 /* This software is released under the license detailed */
4 /* in the file, LICENSE, which is located in the top-level Nalu */
5 /* directory structure */
6 /*------------------------------------------------------------------------*/
7 
8 #ifndef MasterElementFunctions_h
9 #define MasterElementFunctions_h
10 
11 #include <AlgTraits.h>
12 
13 #include <master_element/MasterElement.h>
14 #include <master_element/TensorOps.h>
15 
16 #include <SimdInterface.h>
17 #include <Kokkos_Core.hpp>
18 
19 #include <stk_util/environment/ReportHandler.hpp>
20 
21 #include <vector>
22 #include <cstdlib>
23 #include <stdexcept>
24 #include <string>
25 #include <array>
26 #include <type_traits>
27 
28 namespace sierra {
29 namespace nalu {
30 
31  template <typename AlgTraits, typename GradViewType, typename CoordViewType, typename OutputViewType>
32  void generic_grad_op_3d(GradViewType referenceGradWeights, CoordViewType coords, OutputViewType weights)
33  {
34  using ftype = typename CoordViewType::value_type;
35  static_assert(std::is_same<ftype, typename GradViewType::value_type>::value, "Incompatiable value type for views");
36  static_assert(std::is_same<ftype, typename OutputViewType::value_type>::value, "Incompatiable value type for views");
37  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
38  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
39  static_assert(OutputViewType::Rank == 3, "Weight view assumed to be 3D");
40  static_assert(AlgTraits::nDim_ == 3, "3D method");
41 
42  ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
43  ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
44  ThrowAssert(weights.extent(0) == referenceGradWeights.extent(0));
45  ThrowAssert(weights.extent(1) == referenceGradWeights.extent(1));
46  ThrowAssert(weights.extent(2) == referenceGradWeights.extent(2));
47 
48  for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
49  ftype jact[3][3] = {
50  {0.0, 0.0, 0.0},
51  {0.0, 0.0, 0.0},
52  {0.0, 0.0, 0.0}
53  };
54 
55  ftype refGrad[AlgTraits::nodesPerElement_][3];
56  for (int n = 0; n < AlgTraits::nodesPerElement_; ++n) {
57  refGrad[n][0] = referenceGradWeights(ip, n, 0);
58  refGrad[n][1] = referenceGradWeights(ip, n, 1);
59  refGrad[n][2] = referenceGradWeights(ip, n, 2);
60 
61  jact[0][0] += refGrad[n][0] * coords(n, 0);
62  jact[0][1] += refGrad[n][1] * coords(n, 0);
63  jact[0][2] += refGrad[n][2] * coords(n, 0);
64 
65  jact[1][0] += refGrad[n][0] * coords(n, 1);
66  jact[1][1] += refGrad[n][1] * coords(n, 1);
67  jact[1][2] += refGrad[n][2] * coords(n, 1);
68 
69  jact[2][0] += refGrad[n][0] * coords(n, 2);
70  jact[2][1] += refGrad[n][1] * coords(n, 2);
71  jact[2][2] += refGrad[n][2] * coords(n, 2);
72  }
73 
74  ftype adjJac[3][3];
75  adjJac[0][0] = jact[1][1] * jact[2][2] - jact[2][1] * jact[1][2];
76  adjJac[0][1] = jact[1][2] * jact[2][0] - jact[2][2] * jact[1][0];
77  adjJac[0][2] = jact[1][0] * jact[2][1] - jact[2][0] * jact[1][1];
78 
79  adjJac[1][0] = jact[0][2] * jact[2][1] - jact[2][2] * jact[0][1];
80  adjJac[1][1] = jact[0][0] * jact[2][2] - jact[2][0] * jact[0][2];
81  adjJac[1][2] = jact[0][1] * jact[2][0] - jact[2][1] * jact[0][0];
82 
83  adjJac[2][0] = jact[0][1] * jact[1][2] - jact[1][1] * jact[0][2];
84  adjJac[2][1] = jact[0][2] * jact[1][0] - jact[1][2] * jact[0][0];
85  adjJac[2][2] = jact[0][0] * jact[1][1] - jact[1][0] * jact[0][1];
86 
87  ThrowAssertMsg(
88  stk::simd::are_any(
89  jact[0][0] * adjJac[0][0] + jact[1][0] * adjJac[1][0] + jact[2][0] * adjJac[2][0]
90  > tiny_positive_value()
91  ),
92  "Problem with Jacobian determinant"
93  );
94 
95  const ftype inv_detj = ftype(1.0) /
96  (jact[0][0] * adjJac[0][0] + jact[1][0] * adjJac[1][0] + jact[2][0] * adjJac[2][0]);
97 
98  for (int n = 0; n < AlgTraits::nodesPerElement_; ++n) {
99  weights(ip, n, 0) = inv_detj *
100  (adjJac[0][0] * refGrad[n][0] + adjJac[0][1] * refGrad[n][1] + adjJac[0][2] * refGrad[n][2]);
101 
102  weights(ip, n, 1) = inv_detj *
103  (adjJac[1][0] * refGrad[n][0] + adjJac[1][1] * refGrad[n][1] + adjJac[1][2] * refGrad[n][2]);
104 
105  weights(ip, n, 2) = inv_detj *
106  (adjJac[2][0] * refGrad[n][0] + adjJac[2][1] * refGrad[n][1] + adjJac[2][2] * refGrad[n][2]);
107  }
108  }
109  }
110 
111  template <typename AlgTraits, typename GradViewType, typename CoordViewType, typename OutputViewType>
112  void generic_gij_3d(
113  GradViewType referenceGradWeights,
114  CoordViewType coords,
115  OutputViewType gup,
116  OutputViewType glo)
117  {
118  using ftype = typename CoordViewType::value_type;
119  static_assert(std::is_same<ftype, typename GradViewType::value_type>::value,
120  "Incompatiable value type for views");
121  static_assert(std::is_same<ftype, typename OutputViewType::value_type>::value,
122  "Incompatiable value type for views");
123  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
124  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
125  static_assert(OutputViewType::Rank == 3, "gij view assumed to be 3D");
126  static_assert(AlgTraits::nDim_ == 3, "3D method");
127 
128  for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
129 
130  ftype jac[3][3] = { {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0} };
131  for (int n = 0; n < AlgTraits::nodesPerElement_; ++n) {
132  jac[0][0] += referenceGradWeights(ip, n, 0) * coords(n, 0);
133  jac[0][1] += referenceGradWeights(ip, n, 1) * coords(n, 0);
134  jac[0][2] += referenceGradWeights(ip, n, 2) * coords(n, 0);
135 
136  jac[1][0] += referenceGradWeights(ip, n, 0) * coords(n, 1);
137  jac[1][1] += referenceGradWeights(ip, n, 1) * coords(n, 1);
138  jac[1][2] += referenceGradWeights(ip, n, 2) * coords(n, 1);
139 
140  jac[2][0] += referenceGradWeights(ip, n, 0) * coords(n, 2);
141  jac[2][1] += referenceGradWeights(ip, n, 1) * coords(n, 2);
142  jac[2][2] += referenceGradWeights(ip, n, 2) * coords(n, 2);
143  }
144 
145  gup(ip, 0, 0) = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1] + jac[0][2] * jac[0][2];
146  gup(ip, 0, 1) = jac[0][0] * jac[1][0] + jac[0][1] * jac[1][1] + jac[0][2] * jac[1][2];
147  gup(ip, 0, 2) = jac[0][0] * jac[2][0] + jac[0][1] * jac[2][1] + jac[0][2] * jac[2][2];
148 
149  gup(ip, 1, 0) = gup(ip, 0, 1);
150  gup(ip, 1, 1) = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1] + jac[1][2] * jac[1][2];
151  gup(ip, 1, 2) = jac[1][0] * jac[2][0] + jac[1][1] * jac[2][1] + jac[1][2] * jac[2][2];
152 
153  gup(ip, 2, 0) = gup(ip, 0, 2);
154  gup(ip, 2, 1) = gup(ip, 1, 2);
155  gup(ip, 2, 2) = jac[2][0] * jac[2][0] + jac[2][1] * jac[2][1] + jac[2][2] * jac[2][2];
156 
157  // the covariant is the inverse of the contravariant by definition
158  // gUpper is symmetric
159  const ftype inv_detj = ftype(1.0) / (
160  gup(ip, 0, 0) * ( gup(ip, 1, 1) * gup(ip, 2, 2) - gup(ip, 1, 2) * gup(ip, 1, 2) )
161  - gup(ip, 0, 1) * ( gup(ip, 0, 1) * gup(ip, 2, 2) - gup(ip, 1, 2) * gup(ip, 0, 2) )
162  + gup(ip, 0, 2) * ( gup(ip, 0, 1) * gup(ip, 1, 2) - gup(ip, 1, 1) * gup(ip, 0, 2) )
163  );
164 
165  glo(ip, 0, 0) = inv_detj * (gup(ip, 1, 1) * gup(ip, 2, 2) - gup(ip, 1, 2) * gup(ip, 1, 2));
166  glo(ip, 0, 1) = inv_detj * (gup(ip, 0, 2) * gup(ip, 1, 2) - gup(ip, 0, 1) * gup(ip, 2, 2));
167  glo(ip, 0, 2) = inv_detj * (gup(ip, 0, 1) * gup(ip, 1, 2) - gup(ip, 0, 2) * gup(ip, 1, 1));
168 
169  glo(ip, 1, 0) = glo(ip, 0, 1);
170  glo(ip, 1, 1) = inv_detj * (gup(ip, 0, 0) * gup(ip, 2, 2) - gup(ip, 0, 2) * gup(ip, 0, 2));
171  glo(ip, 1, 2) = inv_detj * (gup(ip, 0, 2) * gup(ip, 0, 1) - gup(ip, 0, 0) * gup(ip, 1, 2));
172 
173 
174  glo(ip, 2, 0) = glo(ip, 0, 2);
175  glo(ip, 2, 1) = glo(ip, 1, 2);
176  glo(ip, 2, 2) = inv_detj * (gup(ip, 0, 0) * gup(ip, 1, 1) - gup(ip, 0, 1) * gup(ip, 0, 1));
177  }
178  }
179 
180  template <typename AlgTraits, typename GradViewType, typename CoordViewType, typename OutputViewType>
181  void generic_determinant_3d(GradViewType referenceGradWeights, CoordViewType coords, OutputViewType detj)
182  {
183  using ftype = typename CoordViewType::value_type;
184  static_assert(std::is_same<ftype, typename GradViewType::value_type>::value, "Incompatiable value type for views");
185  static_assert(std::is_same<ftype, typename OutputViewType::value_type>::value, "Incompatiable value type for views");
186  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
187  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
188  static_assert(OutputViewType::Rank == 1, "Weight view assumed to be 1D");
189  static_assert(AlgTraits::nDim_ == 3, "3D method");
190 
191  ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
192  ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
193 
194  ThrowAssert(detj.extent(0) == referenceGradWeights.extent(0));
195 
196  for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
197  ftype jac[3][3] = { {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0} };
198  for (int n = 0; n < AlgTraits::nodesPerElement_; ++n) {
199  jac[0][0] += referenceGradWeights(ip, n, 0) * coords(n, 0);
200  jac[0][1] += referenceGradWeights(ip, n, 1) * coords(n, 0);
201  jac[0][2] += referenceGradWeights(ip, n, 2) * coords(n, 0);
202 
203  jac[1][0] += referenceGradWeights(ip, n, 0) * coords(n, 1);
204  jac[1][1] += referenceGradWeights(ip, n, 1) * coords(n, 1);
205  jac[1][2] += referenceGradWeights(ip, n, 2) * coords(n, 1);
206 
207  jac[2][0] += referenceGradWeights(ip, n, 0) * coords(n, 2);
208  jac[2][1] += referenceGradWeights(ip, n, 1) * coords(n, 2);
209  jac[2][2] += referenceGradWeights(ip, n, 2) * coords(n, 2);
210  }
211  detj(ip) = determinant33(&jac[0][0]);
212  }
213  }
214 
215 } // namespace nalu
216 } // namespace Sierra
217 
218 #endif
sierra::nalu::generic_gij_3d
void generic_gij_3d(GradViewType referenceGradWeights, CoordViewType coords, OutputViewType gup, OutputViewType glo)
Definition: MasterElementFunctions.h:112
sierra
Definition: ABLForcingAlgorithm.C:26
sierra::nalu::generic_determinant_3d
void generic_determinant_3d(GradViewType referenceGradWeights, CoordViewType coords, OutputViewType detj)
Definition: MasterElementFunctions.h:181
SimdInterface.h
STK SIMD Interface.
sierra::nalu::determinant33
KOKKOS_FORCEINLINE_FUNCTION ScalarType determinant33(const ScalarType *mat)
Definition: TensorOps.h:69
sierra::nalu::tiny_positive_value
constexpr double tiny_positive_value()
Definition: TensorOps.h:29
sierra::nalu::generic_grad_op_3d
void generic_grad_op_3d(GradViewType referenceGradWeights, CoordViewType coords, OutputViewType weights)
Definition: MasterElementFunctions.h:32
AlgTraits.h
MasterElement.h
TensorOps.h
Generated by   doxygen 1.8.11