MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
nonlininteg_vecconvection_pa_diag.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11#pragma once
12
13#include "../../config/config.hpp"
14#include "../../general/forall.hpp"
15#include "../../linalg/dtensor.hpp"
16#include "../kernels.hpp"
17#include "../nonlininteg.hpp"
18
19namespace mfem
20{
21
22/// \cond DO_NOT_DOCUMENT
23
24namespace internal
25{
26
27template<int T_D1D = 0, int T_Q1D = 0>
28inline void SmemPAConvectionNLGradDiagonal2D(const int NE,
29 const real_t *b,
30 const real_t *g,
31 const real_t *a,
32 const real_t *u,
33 real_t *de,
34 const int d1d,
35 const int q1d)
36{
37 static constexpr int VDIM = 2, DIM = 2;
38 const int D1D = T_D1D ? T_D1D : d1d;
39 const int Q1D = T_Q1D ? T_Q1D : q1d;
40
41 const auto A = Reshape(a, VDIM, DIM, Q1D, Q1D, NE);
42 const auto U = Reshape(u, D1D, D1D, VDIM, NE);
43 auto D = Reshape(de, D1D, D1D, VDIM, NE);
44
45 mfem::forall_2D<T_Q1D * T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
46 {
47 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
48 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
49
50 MFEM_SHARED real_t sM[3][MQ1][MQ1], sQ[3][MQ1][MQ1];
51 MFEM_SHARED real_t sB[MD1][MQ1], sG[MD1][MQ1];
52
53 kernels::internal::v_regs2d_t<VDIM, MQ1> r0, r1;
54 kernels::internal::vd_regs2d_t<VDIM, DIM, MQ1> g0, g1;
55
56 kernels::internal::LoadMatrix(D1D, Q1D, b, sB);
57 kernels::internal::LoadMatrix(D1D, Q1D, g, sG);
58
59 kernels::internal::LoadDofs2d(e, D1D, U, r0);
60 kernels::internal::Eval2d(D1D, Q1D, sM[0], sB, r0, r1);
61
62 kernels::internal::LoadDofs2d(e, D1D, U, g0);
63 kernels::internal::Grad2d(D1D, Q1D, sM[0], sB, sG, g0, g1);
64
65 for (int v = 0; v < VDIM; ++v)
66 {
67 future::tensor<real_t, VDIM> e_v = {};
68 e_v[v] = real_t(1);
69 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
70 {
71 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
72 {
73 const future::tensor<real_t, VDIM> u_val =
74 {
75 r1[0][qy][qx], r1[1][qy][qx]
76 };
77 const future::tensor<real_t, VDIM, DIM> Q_adj =
78 {
79 { { A(0, 0, qx, qy, e), A(1, 0, qx, qy, e) },
80 { A(0, 1, qx, qy, e), A(1, 1, qx, qy, e) }
81 }
82 };
83 const future::tensor<real_t, VDIM, DIM> grad_U =
84 {
85 { { g1[0][0][qy][qx], g1[1][0][qy][qx] },
86 { g1[0][1][qy][qx], g1[1][1][qy][qx] }
87 }
88 };
89 const auto one = Q_adj * u_val;
90 const auto two = transpose(grad_U) * (Q_adj * e_v);
91 sQ[0][qx][qy] = one[0];
92 sQ[1][qx][qy] = one[1];
93 sQ[2][qx][qy] = two[v];
94 }
95 }
96 MFEM_SYNC_THREAD;
97
98 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
99 {
100 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
101 {
102 real_t s[3] = {};
103 for (int qy = 0; qy < Q1D; ++qy)
104 {
105 const real_t By = sB[dy][qy], Gy = sG[dy][qy];
106 s[0] += By * By * sQ[0][qx][qy];
107 s[1] += Gy * By * sQ[1][qx][qy];
108 s[2] += By * By * sQ[2][qx][qy];
109 }
110 sM[0][qx][dy] = s[0];
111 sM[1][qx][dy] = s[1];
112 sM[2][qx][dy] = s[2];
113 }
114 }
115 MFEM_SYNC_THREAD;
116
117 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
118 {
119 MFEM_FOREACH_THREAD_DIRECT(dx, x, D1D)
120 {
121 real_t d = 0.0;
122 for (int qx = 0; qx < Q1D; ++qx)
123 {
124 const real_t Bx = sB[dx][qx], Gx = sG[dx][qx];
125 d += Gx * Bx * sM[0][qx][dy] +
126 Bx * Bx * sM[1][qx][dy] +
127 Bx * Bx * sM[2][qx][dy];
128 }
129 D(dx, dy, v, e) += d;
130 }
131 }
132 MFEM_SYNC_THREAD;
133 }
134 });
135}
136
137template<int T_D1D = 0, int T_Q1D = 0>
138inline void SmemPAConvectionNLGradDiagonal3D(const int NE,
139 const real_t *b,
140 const real_t *g,
141 const real_t *a,
142 const real_t *u,
143 real_t *de,
144 const int d1d,
145 const int q1d)
146{
147 static constexpr int VDIM = 3, DIM = 3;
148 const int D1D = T_D1D ? T_D1D : d1d;
149 const int Q1D = T_Q1D ? T_Q1D : q1d;
150
151 const auto A = Reshape(a, VDIM, DIM, Q1D, Q1D, Q1D, NE);
152 const auto U = Reshape(u, D1D, D1D, D1D, VDIM, NE);
153 auto D = Reshape(de, D1D, D1D, D1D, VDIM, NE);
154
155 mfem::forall_2D<T_Q1D * T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
156 {
157 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
158 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
159
160 MFEM_SHARED real_t sM[4][MQ1][MQ1], sQ[4][MQ1][MQ1];
161 MFEM_SHARED real_t sB[MD1][MQ1], sG[MD1][MQ1];
162
163 kernels::internal::v_regs3d_t<VDIM, MQ1> r0, r1;
164 kernels::internal::vd_regs3d_t<VDIM, DIM, MQ1> g0, g1;
165
166 kernels::internal::LoadMatrix(D1D, Q1D, b, sB);
167 kernels::internal::LoadMatrix(D1D, Q1D, g, sG);
168
169 kernels::internal::LoadDofs3d(e, D1D, U, r0);
170 kernels::internal::Eval3d(D1D, Q1D, sM[0], sB, r0, r1);
171
172 kernels::internal::LoadDofs3d(e, D1D, U, g0);
173 kernels::internal::Grad3d(D1D, Q1D, sM[0], sB, sG, g0, g1);
174
175 for (int v = 0; v < VDIM; ++v)
176 {
177 future::tensor<real_t, VDIM> e_v = {};
178 e_v[v] = real_t(1);
179 for (int dz = 0; dz < D1D; ++dz)
180 {
181 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
182 {
183 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
184 {
185 real_t s[4] = {};
186 for (int qz = 0; qz < Q1D; ++qz)
187 {
188 const future::tensor<real_t, VDIM> u_val =
189 {
190 r1[0][qz][qy][qx], r1[1][qz][qy][qx], r1[2][qz][qy][qx]
191 };
192 const future::tensor<real_t, VDIM, DIM> Q_adj = {{
193 {A(0,0,qx,qy,qz,e), A(1,0,qx,qy,qz,e), A(2,0,qx,qy,qz,e)},
194 {A(0,1,qx,qy,qz,e), A(1,1,qx,qy,qz,e), A(2,1,qx,qy,qz,e)},
195 {A(0,2,qx,qy,qz,e), A(1,2,qx,qy,qz,e), A(2,2,qx,qy,qz,e)}
196 }
197 };
198 const future::tensor<real_t, VDIM, DIM> grad_U = {{
199 {g1[0][0][qz][qy][qx], g1[1][0][qz][qy][qx], g1[2][0][qz][qy][qx]},
200 {g1[0][1][qz][qy][qx], g1[1][1][qz][qy][qx], g1[2][1][qz][qy][qx]},
201 {g1[0][2][qz][qy][qx], g1[1][2][qz][qy][qx], g1[2][2][qz][qy][qx]}
202 }
203 };
204 const auto one = Q_adj * u_val;
205 const auto two = transpose(grad_U) * (Q_adj * e_v);
206
207 const real_t Bz = sB[dz][qz], Gz = sG[dz][qz];
208 s[0] += one[0] * Bz * Bz;
209 s[1] += one[1] * Bz * Bz;
210 s[2] += one[2] * Bz * Gz;
211 s[3] += two[v] * Bz * Bz;
212 }
213 sQ[0][qx][qy] = s[0];
214 sQ[1][qx][qy] = s[1];
215 sQ[2][qx][qy] = s[2];
216 sQ[3][qx][qy] = s[3];
217 }
218 }
219 MFEM_SYNC_THREAD;
220
221 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
222 {
223 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
224 {
225 real_t s[4] = {};
226 for (int qy = 0; qy < Q1D; ++qy)
227 {
228 const real_t By = sB[dy][qy], Gy = sG[dy][qy];
229 s[0] += By * By * sQ[0][qx][qy];
230 s[1] += Gy * By * sQ[1][qx][qy];
231 s[2] += By * By * sQ[2][qx][qy];
232 s[3] += By * By * sQ[3][qx][qy];
233 }
234 sM[0][dy][qx] = s[0];
235 sM[1][dy][qx] = s[1];
236 sM[2][dy][qx] = s[2];
237 sM[3][dy][qx] = s[3];
238 }
239 }
240 MFEM_SYNC_THREAD;
241
242 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
243 {
244 MFEM_FOREACH_THREAD_DIRECT(dx, x, D1D)
245 {
246 real_t d = 0.0;
247 for (int qx = 0; qx < Q1D; ++qx)
248 {
249 const real_t Bx = sB[dx][qx], Gx = sG[dx][qx];
250 d += Gx * Bx * sM[0][dy][qx];
251 d += Bx * Bx * sM[1][dy][qx];
252 d += Bx * Bx * sM[2][dy][qx];
253 d += Bx * Bx * sM[3][dy][qx];
254 }
255 D(dx, dy, dz, v, e) += d;
256 }
257 }
258 MFEM_SYNC_THREAD;
259 }
260 }
261 });
262}
263
264} // namespace internal
265
266template<int T_D1D, int T_Q1D>
267VectorConvectionNLFIntegrator::GradDiagPAType
268VectorConvectionNLFIntegrator::GradDiagPA2D::Kernel()
269{
270 static_assert(T_D1D <= T_Q1D, "d1d > q1d is not supported");
271 return internal::SmemPAConvectionNLGradDiagonal2D<T_D1D, T_Q1D>;
272}
273
274inline VectorConvectionNLFIntegrator::GradDiagPAType
275VectorConvectionNLFIntegrator::GradDiagPA2D::Fallback(int d1d, int q1d)
276{
277 MFEM_VERIFY(d1d <= q1d, "d1d > q1d is not supported");
278 MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
279 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
280 return internal::SmemPAConvectionNLGradDiagonal2D<>;
281}
282
283template<int T_D1D, int T_Q1D>
284VectorConvectionNLFIntegrator::GradDiagPAType
285VectorConvectionNLFIntegrator::GradDiagPA3D::Kernel()
286{
287 static_assert(T_D1D <= T_Q1D, "d1d > q1d is not supported");
288 return internal::SmemPAConvectionNLGradDiagonal3D<T_D1D, T_Q1D>;
289}
290
291inline VectorConvectionNLFIntegrator::GradDiagPAType
292VectorConvectionNLFIntegrator::GradDiagPA3D::Fallback(int d1d, int q1d)
293{
294 MFEM_VERIFY(d1d <= q1d, "d1d > q1d is not supported");
295 MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
296 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
297 return internal::SmemPAConvectionNLGradDiagonal3D<>;
298}
299
300/// \endcond DO_NOT_DOCUMENT
301
302} // namespace mfem
mfem::VectorConvectionNLFIntegrator::GradDiagPAType
void(*)(const int ne, const real_t *B, const real_t *G, const real_t *A, const real_t *u, real_t *y, const int d1d, const int q1d) GradDiagPAType
Definition nonlininteg.hpp:439
b
real_t b
Definition lissajous.cpp:42
a
real_t a
Definition lissajous.cpp:41
DIM
constexpr int DIM
Definition minimal-surface.cpp:71
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
mfem::future::transpose
MFEM_HOST_DEVICE tensor< T, n, m > transpose(const tensor< T, m, n > &A)
Returns the transpose of the matrix.
Definition tensor.hpp:1388
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::Reshape
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:138
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:1220
mfem::real_t
float real_t
Definition config.hpp:46
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0