PACMAN 0.1.0
Portable Algorithms for Coupling, Mapping, and Adaptive iNterpolation
Loading...
Searching...
No Matches
NewtonKokkos.hpp
Go to the documentation of this file.
1//
2// This file is subject to the terms and conditions defined in
3// file 'LICENSE', which is part of this source code package.
4//
5
6#pragma once
7
8#include <KokkosBatched_Gemm_Decl.hpp>
9#include <KokkosBatched_LU_Decl.hpp>
10#include <KokkosBatched_Trsv_Decl.hpp>
11#include <KokkosBlas1_axpby.hpp>
12#include <KokkosBlas2_serial_gemv.hpp>
13
14#include "finite_elements/utils/FETools.hpp"
15
16namespace PACMAN {
17namespace FiniteElements {
18
20#define D_XI_ETA_PHI_TOL_SQUARED 1e-10
22#define MAX_NEWTON_ITER 15
23
24#define CaseSwitch(CellType_STRUCT) \
25 case CellType_STRUCT: { \
26 return ApplyNewton<ExecSpace, \
27 LagrangeSpaceGeo<ExecSpace, CellType_STRUCT>, Dim>( \
28 Xcoor, targetPoint, weights, forceEvaluation); \
29 break; \
30 }
31
48template <typename ExecSpace, typename FEspace, int_t Dim>
49KOKKOS_FUNCTION bool
50ApplyNewton(const Kokkos::View<coordinates_t **, ExecSpace> Xcoor,
51 const Kokkos::View<coordinates_t[Dim], ExecSpace> targetPoint,
52 Kokkos::View<fp_t *, ExecSpace> weights,
53 const bool forceEvaluation = false) {
54
55 // Shape function values at the current xi,eta,phi stored in weights
56 // KOKKOS_ASSERT(Xcoor.extent_int(0) == FEspace::numberOfShapeFunctions);
57
58 Kokkos::Array<fp_t, 3> xietaphi_arr = {0.5, 0.5, 0.5};
59 Kokkos::View<fp_t *, ExecSpace> xietaphi(xietaphi_arr.data(), Dim);
60
61 // Coordinates in the global system [DIM]
62 Kokkos::Array<fp_t, Dim> f_arr;
63 Kokkos::View<fp_t *, ExecSpace> f(f_arr.data(), Dim);
64 // Derivatives of shape functions at the current xi,eta,phi [FE::DIM;
65 // FE::NSHAPES]
66 Kokkos::Array<fp_t, FEspace::dimensionality * FEspace::numberOfShapeFunctions>
67 dN_arr;
68 Kokkos::View<fp_t **, ExecSpace> dN(dN_arr.data(), FEspace::dimensionality,
69 FEspace::numberOfShapeFunctions);
70 // Derivative of f w.r.t xi,eta,phi [FE::DIM; DIM]
71 Kokkos::Array<fp_t, FEspace::dimensionality * Dim> dNXcoor_arr;
72 Kokkos::View<fp_t **, ExecSpace> dNXcoor(dNXcoor_arr.data(),
73 FEspace::dimensionality, Dim);
74 // Numerical derivative of f w.r.t xi,eta,phi [FE::DIM]
75 Kokkos::Array<fp_t, FEspace::dimensionality> dfNum_arr;
76 Kokkos::View<fp_t *, ExecSpace> dfNum(dfNum_arr.data(),
77 FEspace::dimensionality);
78 // Newton function to solve h * dxietaphi = dfNum [FE::DIM, FE::DIM]
79 Kokkos::Array<fp_t, FEspace::dimensionality * FEspace::dimensionality> h_arr;
80 Kokkos::View<fp_t **, ExecSpace> h(h_arr.data(), FEspace::dimensionality,
81 FEspace::dimensionality);
82
83 FEspace::UpdateShapeFunctionsValues(xietaphi_arr[0], xietaphi_arr[1],
84 xietaphi_arr[2], weights);
85 // f= weights.transpose() * xcoor;
86 KokkosBlas::Experimental::serial_gemv('T', fp_consts::one(), Xcoor, weights,
87 fp_consts::zero(), f);
88 // f -= targetPoint;
89 KokkosBlas::serial_axpy(-fp_consts::one(), targetPoint, f);
90
91 int_t cnt = 0;
92 double error = 0.0;
93
94 using namespace KokkosBatched;
95
96 for (; cnt < MAX_NEWTON_ITER; ++cnt) {
97
98 FEspace::UpdateShapeFunctionsDerValues(xietaphi_arr[0], xietaphi_arr[1],
99 xietaphi_arr[2], dN);
100
101 // dfNum = (f * ((dN * xcoor).transpose())).transpose();
102 // here h = (dN * xcoor)
103 SerialGemm<Trans::NoTranspose, Trans::NoTranspose,
104 Algo::Gemm::Unblocked>::invoke(fp_consts::one(), dN, Xcoor,
105 fp_consts::zero(), dNXcoor);
106 KokkosBlas::Experimental::serial_gemv('N', fp_consts::one(), dNXcoor, f,
107 fp_consts::zero(), dfNum);
108
109 // h = (dN * xcoor) * (dN * xcoor).transpose();
110 SerialGemm<Trans::NoTranspose, Trans::Transpose,
111 Algo::Gemm::Unblocked>::invoke(fp_consts::one(), dNXcoor,
112 dNXcoor, fp_consts::zero(), h);
113
114 if (cnt > 9) {
115 Kokkos::Array<fp_t, FEspace::dimensionality * FEspace::dimensionality>
116 ddNi_arr;
117 Kokkos::View<fp_t **, ExecSpace> ddNi(
118 ddNi_arr.data(), FEspace::dimensionality, FEspace::dimensionality);
119 for (int_t j = 0; j < Dim; j++) { // cols
120 for (int_t i = 0; i < FEspace::numberOfShapeFunctions; i++) { // rows
121 for (int_t k = 0;
122 k < FEspace::dimensionality * FEspace::dimensionality; k++) {
123 ddNi_arr[k] = 0.0;
124 }
125 FEspace::UpdateShapeFunctionsDerDerValues(
126 i, xietaphi_arr[0], xietaphi_arr[1], xietaphi_arr[2], ddNi);
127 KokkosBlas::serial_axpy(Xcoor(i, j), ddNi, h);
128 }
129 }
130 }
131
132 SerialLU<Algo::LU::Unblocked>::invoke(h);
133 SerialTrsv<Uplo::Lower, Trans::NoTranspose, Diag::Unit,
134 Algo::Trsv::Unblocked>::invoke(fp_consts::one(), h, dfNum);
135 SerialTrsv<Uplo::Upper, Trans::NoTranspose, Diag::NonUnit,
136 Algo::Trsv::Unblocked>::invoke(fp_consts::one(), h, dfNum);
137
138 for (int_t j = 0; j < FEspace::dimensionality; j++) {
139 xietaphi(j) -= dfNum(j);
140 }
141
142 FEspace::UpdateShapeFunctionsValues(xietaphi_arr[0], xietaphi_arr[1],
143 xietaphi_arr[2], weights);
144 KokkosBlas::Experimental::serial_gemv('T', fp_consts::one(), Xcoor, weights,
145 fp_consts::zero(),
146 f); // f= weights.transpose() * xcoor;
147 KokkosBlas::serial_axpy(-fp_consts::one(), targetPoint,
148 f); // f -= targetPoint;
149
150 // const fp_t error = KokkosBlas::dot(f, f) / KokkosBlas::dot(targetPoint,
151 // targetPoint);
152 fp_t norm2 = fp_consts::zero();
153 fp_t scale = fp_consts::zero();
154 for (int_t j = 0; j < Dim; j++) {
155 norm2 += f(j) * f(j);
156 scale += targetPoint(j) * targetPoint(j);
157 }
158 error = norm2 / scale;
159 if (error < D_XI_ETA_PHI_TOL_SQUARED) {
160 break;
161 }
162 }
163
164 if (forceEvaluation) {
165 return true;
166 }
167
168 if ((cnt < MAX_NEWTON_ITER) &&
169 FEspace::isInside(xietaphi_arr[0], xietaphi_arr[1], xietaphi_arr[2])) {
170 return true;
171 }
172
173 return false;
174}
175
189template <typename ExecSpace, int_t Dim>
190KOKKOS_FUNCTION bool ApplyNewtonOnElement(
191 const cell_t type, const Kokkos::View<coordinates_t **, ExecSpace> Xcoor,
192 const Kokkos::View<coordinates_t[Dim], ExecSpace> targetPoint,
193 Kokkos::View<fp_t *, ExecSpace> weights,
194 const bool forceEvaluation = false) {
195 // clang-format off
196 switch (static_cast<CellType>(type)) {
197 CaseSwitch(CellType::VTK_LINE)
198 CaseSwitch(CellType::VTK_QUADRATIC_HEXAHEDRON)
199 CaseSwitch(CellType::VTK_QUADRATIC_EDGE)
200 CaseSwitch(CellType::VTK_HEXAHEDRON)
201 CaseSwitch(CellType::VTK_QUADRATIC_PYRAMID)
202 CaseSwitch(CellType::VTK_PYRAMID)
203 CaseSwitch(CellType::VTK_QUAD)
204 CaseSwitch(CellType::VTK_QUADRATIC_QUAD)
205 CaseSwitch(CellType::VTK_QUADRATIC_TETRA)
206 CaseSwitch(CellType::VTK_TETRA)
207 CaseSwitch(CellType::VTK_TRIANGLE)
208 CaseSwitch(CellType::VTK_QUADRATIC_TRIANGLE)
209 CaseSwitch(CellType::VTK_QUADRATIC_WEDGE)
210 CaseSwitch(CellType::VTK_WEDGE)
211// LCOV_EXCL_START
212 default:
213 char errorMsg[] = "Unexpected element type XXX";
214 errorMsg[sizeof(errorMsg) - 4] = '0' + ((int)type / 100);
215 errorMsg[sizeof(errorMsg) - 3] = '0' + (((int)type / 10) % 10);
216 errorMsg[sizeof(errorMsg) - 2] = '0' + ((int)type % 10);
217 Kokkos::abort(errorMsg);
218 return false;
219// LCOV_EXCL_STOP
220 }
221 // clang-format on
222 return false;
223}
224} // namespace FiniteElements
225} // namespace PACMAN
FETools.hpp
Finite-element geometric helpers and Lagrange-space specializations.
CaseSwitch
#define CaseSwitch(CellType_STRUCT)
Definition NewtonKokkos.hpp:24
MAX_NEWTON_ITER
#define MAX_NEWTON_ITER
Maximal number of newton iters.
Definition NewtonKokkos.hpp:22
D_XI_ETA_PHI_TOL_SQUARED
#define D_XI_ETA_PHI_TOL_SQUARED
Tolerance to stop the newton, dxietaphi.squaredNorm().
Definition NewtonKokkos.hpp:20
PACMAN::FiniteElements::ApplyNewtonOnElement
KOKKOS_FUNCTION bool ApplyNewtonOnElement(const cell_t type, const Kokkos::View< coordinates_t **, ExecSpace > Xcoor, const Kokkos::View< coordinates_t[Dim], ExecSpace > targetPoint, Kokkos::View< fp_t *, ExecSpace > weights, const bool forceEvaluation=false)
Dispatch Newton-based FE inversion to the correct element type.
Definition NewtonKokkos.hpp:190
PACMAN::FiniteElements::ApplyNewton
KOKKOS_FUNCTION bool ApplyNewton(const Kokkos::View< coordinates_t **, ExecSpace > Xcoor, const Kokkos::View< coordinates_t[Dim], ExecSpace > targetPoint, Kokkos::View< fp_t *, ExecSpace > weights, const bool forceEvaluation=false)
Solve FE inverse mapping with Newton iterations and compute interpolation weights.
Definition NewtonKokkos.hpp:50
PACMAN::fp_consts::zero
KOKKOS_INLINE_FUNCTION consteval fp_t zero(void)
Definition types.hpp:92
PACMAN::fp_consts::one
KOKKOS_INLINE_FUNCTION consteval fp_t one(void)
Definition types.hpp:95
PACMAN::fp_t
double fp_t
Definition types.hpp:15
PACMAN::coordinates_t
double coordinates_t
Definition types.hpp:16
PACMAN::cell_t
int32_t cell_t
Definition types.hpp:25
PACMAN::int_t
int32_t int_t
Definition types.hpp:18
PACMAN::CellType
CellType
Definition types.hpp:60
PACMAN::CellType::VTK_QUADRATIC_HEXAHEDRON
@ VTK_QUADRATIC_HEXAHEDRON
Definition types.hpp:79
PACMAN::CellType::VTK_QUADRATIC_PYRAMID
@ VTK_QUADRATIC_PYRAMID
Definition types.hpp:81
PACMAN::CellType::VTK_QUADRATIC_TRIANGLE
@ VTK_QUADRATIC_TRIANGLE
Definition types.hpp:71