MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_vecdiv_pa.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/array.hpp"
15#include "../../general/forall.hpp"
16#include "../../linalg/dtensor.hpp"
17#include "../../linalg/vector.hpp"
18#include "../bilininteg.hpp"
19#include "../kernels.hpp"
20
21namespace mfem
22{
23
24/// \cond DO_NOT_DOCUMENT
25
26namespace internal
27{
28
29// Shared memory PA Divergence Apply 2D kernel
30template<int T_TR_D1D = 0, int T_TE_D1D = 0, int T_Q1D = 0>
31inline void SmemPADivergenceApply2D(const int NE,
32 const Array<real_t> &b_,
33 const Array<real_t> &g_,
34 const Array<real_t> &bt_,
35 const Vector &q_,
36 const Vector &x_,
37 Vector &y_,
38 const int tr_d1d = 0,
39 const int te_d1d = 0,
40 const int q1d = 0)
41{
42 const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
43 const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
44 const int Q1D = T_Q1D ? T_Q1D : q1d;
45
46 MFEM_VERIFY(TR_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
47 MFEM_VERIFY(TE_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
48 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
49
50 const auto B = b_.Read(), G = g_.Read(), Bt = bt_.Read();
51 const auto Q = Reshape(q_.Read(), Q1D, Q1D, 2, 2, NE);
52 const auto X = Reshape(x_.Read(), TR_D1D, TR_D1D, 2, NE);
53 auto Y = Reshape(y_.ReadWrite(), TE_D1D, TE_D1D, 1, NE);
54
55 mfem::forall_2D<T_Q1D * T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
56 {
57 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
58
59 MFEM_SHARED real_t smem[MQ1][MQ1];
60 MFEM_SHARED real_t sB[MQ1][MQ1], sG[MQ1][MQ1];
61
62 kernels::internal::vd_regs2d_t<2, 2, MQ1> g0, g1;
63 kernels::internal::v_regs2d_t<1, MQ1> r0, r1;
64
65 kernels::internal::LoadMatrix(TR_D1D, Q1D, B, sB);
66 kernels::internal::LoadMatrix(TR_D1D, Q1D, G, sG);
67
68 kernels::internal::LoadDofs2d(e, TR_D1D, X, g0);
69 kernels::internal::Grad2d(TR_D1D, Q1D, smem, sB, sG, g0, g1);
70
71 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
72 {
73 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
74 {
75 r0[0][qy][qx] =
76 g1[0][0][qy][qx] * Q(qx, qy, 0, 0, e) +
77 g1[0][1][qy][qx] * Q(qx, qy, 1, 0, e) +
78 g1[1][0][qy][qx] * Q(qx, qy, 0, 1, e) +
79 g1[1][1][qy][qx] * Q(qx, qy, 1, 1, e);
80 }
81 }
82 MFEM_SYNC_THREAD;
83
84 kernels::internal::LoadMatrix<MQ1,true>(TE_D1D, Q1D, Bt, sB);
85 kernels::internal::EvalTranspose2d(TE_D1D, Q1D, smem, sB, r0, r1);
86 kernels::internal::WriteDofs2d(e, TE_D1D, r1, Y);
87 });
88}
89
90// Shared memory PA Divergence Apply 2D kernel transpose
91template<int T_TR_D1D = 0, int T_TE_D1D = 0, int T_Q1D = 0>
92inline void SmemPADivergenceApplyTranspose2D(const int NE,
93 const Array<real_t> &bt,
94 const Array<real_t> &gt,
95 const Array<real_t> &b,
96 const Vector &q_,
97 const Vector &x_,
98 Vector &y_,
99 const int tr_d1d = 0,
100 const int te_d1d = 0,
101 const int q1d = 0)
102{
103 const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
104 const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
105 const int Q1D = T_Q1D ? T_Q1D : q1d;
106
107 MFEM_VERIFY(TR_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
108 MFEM_VERIFY(TE_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
109 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
110
111 const auto Bt = bt.Read(), Gt = gt.Read(), B = b.Read();
112 const auto Q = Reshape(q_.Read(), Q1D, Q1D, 2, 2, NE);
113 const auto X = Reshape(x_.Read(), TE_D1D, TE_D1D, 1, NE);
114 auto Y = Reshape(y_.ReadWrite(), TR_D1D, TR_D1D, 2, NE);
115
116 mfem::forall_2D<T_Q1D * T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
117 {
118 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
119
120 MFEM_SHARED real_t smem[MQ1][MQ1];
121 MFEM_SHARED real_t sB[MQ1][MQ1], sG[MQ1][MQ1];
122
123 kernels::internal::v_regs2d_t<1, MQ1> r0, r1;
124 kernels::internal::vd_regs2d_t<2, 2, MQ1> g0, g1;
125
126 kernels::internal::LoadMatrix(TE_D1D, Q1D, B, sB);
127 kernels::internal::LoadDofs2d(e, TE_D1D, X, r0);
128 kernels::internal::Eval2d(TE_D1D, Q1D, smem, sB, r0, r1);
129
130 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
131 {
132 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
133 {
134 g0[0][0][qy][qx] = r1[0][qy][qx] * Q(qx, qy, 0, 0, e);
135 g0[0][1][qy][qx] = r1[0][qy][qx] * Q(qx, qy, 1, 0, e);
136 g0[1][0][qy][qx] = r1[0][qy][qx] * Q(qx, qy, 0, 1, e);
137 g0[1][1][qy][qx] = r1[0][qy][qx] * Q(qx, qy, 1, 1, e);
138 }
139 }
140 MFEM_SYNC_THREAD;
141
142 kernels::internal::LoadMatrix<MQ1,true>(TR_D1D, Q1D, Bt, sB);
143 kernels::internal::LoadMatrix<MQ1,true>(TR_D1D, Q1D, Gt, sG);
144 kernels::internal::GradTranspose2d(TR_D1D, Q1D, smem, sB, sG, g0, g1);
145 kernels::internal::WriteDofs2d(e, TR_D1D, g1, Y);
146 });
147}
148
149// Shared memory PA Divergence Apply 3D kernel transpose
150template<int T_TR_D1D = 0, int T_TE_D1D = 0, int T_Q1D = 0>
151inline void SmemPADivergenceApplyTranspose3D(const int NE,
152 const Array<real_t> &bt,
153 const Array<real_t> &gt,
154 const Array<real_t> &b,
155 const Vector &q_,
156 const Vector &x_,
157 Vector &y_,
158 int tr_d1d = 0,
159 int te_d1d = 0,
160 int q1d = 0)
161{
162 const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
163 const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
164 const int Q1D = T_Q1D ? T_Q1D : q1d;
165
166 MFEM_VERIFY(TR_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
167 MFEM_VERIFY(TE_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
168 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
169
170 const auto Bt = bt.Read(), Gt = gt.Read(), B = b.Read();
171 const auto Q = Reshape(q_.Read(), Q1D, Q1D, Q1D, 3, 3, NE);
172 const auto X = Reshape(x_.Read(), TE_D1D, TE_D1D, TE_D1D, 1, NE);
173 auto Y = Reshape(y_.ReadWrite(), TR_D1D, TR_D1D, TR_D1D, 3, NE);
174
175 mfem::forall_2D<T_Q1D * T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
176 {
177 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
178
179 MFEM_SHARED real_t smem[MQ1][MQ1];
180 MFEM_SHARED real_t sB[MQ1][MQ1], sG[MQ1][MQ1];
181
182 kernels::internal::v_regs3d_t<1, MQ1> r0, r1;
183 kernels::internal::vd_regs3d_t<3, 3, MQ1> g0, g1;
184
185 kernels::internal::LoadMatrix(TE_D1D, Q1D, B, sB);
186 kernels::internal::LoadDofs3d(e, TE_D1D, X, r0);
187 kernels::internal::Eval3d(TE_D1D, Q1D, smem, sB, r0, r1);
188
189 for (int qz = 0; qz < Q1D; qz++)
190 {
191 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
192 {
193 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
194 {
195 const auto r = r1[0][qz][qy][qx];
196 g0[0][0][qz][qy][qx] = r * Q(qx, qy, qz, 0, 0, e);
197 g0[0][1][qz][qy][qx] = r * Q(qx, qy, qz, 1, 0, e);
198 g0[0][2][qz][qy][qx] = r * Q(qx, qy, qz, 2, 0, e);
199
200 g0[1][0][qz][qy][qx] = r * Q(qx, qy, qz, 0, 1, e);
201 g0[1][1][qz][qy][qx] = r * Q(qx, qy, qz, 1, 1, e);
202 g0[1][2][qz][qy][qx] = r * Q(qx, qy, qz, 2, 1, e);
203
204 g0[2][0][qz][qy][qx] = r * Q(qx, qy, qz, 0, 2, e);
205 g0[2][1][qz][qy][qx] = r * Q(qx, qy, qz, 1, 2, e);
206 g0[2][2][qz][qy][qx] = r * Q(qx, qy, qz, 2, 2, e);
207 }
208 }
209 }
210 MFEM_SYNC_THREAD;
211
212 kernels::internal::LoadMatrix<MQ1,true>(TR_D1D, Q1D, Bt, sB);
213 kernels::internal::LoadMatrix<MQ1,true>(TR_D1D, Q1D, Gt, sG);
214 kernels::internal::GradTranspose3d(TR_D1D, Q1D, smem, sB, sG, g0, g1);
215 kernels::internal::WriteDofs3d(e, TR_D1D, g1, Y);
216 });
217}
218
219// Shared memory PA Divergence Apply 3D kernel
220template<int T_TR_D1D = 0, int T_TE_D1D = 0, int T_Q1D = 0>
221inline void SmemPADivergenceApply3D(const int NE,
222 const Array<real_t> &b_,
223 const Array<real_t> &g_,
224 const Array<real_t> &bt_,
225 const Vector &q_,
226 const Vector &x_,
227 Vector &y_,
228 const int tr_d1d = 0,
229 const int te_d1d = 0,
230 const int q1d = 0)
231{
232 const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
233 const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
234 const int Q1D = T_Q1D ? T_Q1D : q1d;
235
236 MFEM_VERIFY(TR_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
237 MFEM_VERIFY(TE_D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
238 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
239
240 const auto B = b_.Read(), G = g_.Read(), Bt = bt_.Read();
241 const auto Q = Reshape(q_.Read(), Q1D, Q1D, Q1D, 3,3, NE);
242 const auto X = Reshape(x_.Read(), TR_D1D, TR_D1D, TR_D1D, 3, NE);
243 auto Y = Reshape(y_.ReadWrite(), TE_D1D, TE_D1D, TE_D1D, 1, NE);
244
245 mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
246 {
247 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
248
249 MFEM_SHARED real_t smem[MQ1][MQ1];
250 MFEM_SHARED real_t sB[MQ1][MQ1], sG[MQ1][MQ1];
251
252 kernels::internal::vd_regs3d_t<3, 3, MQ1> g0, g1;
253 kernels::internal::v_regs3d_t<1, MQ1> r0, r1;
254
255 kernels::internal::LoadMatrix(TR_D1D, Q1D, B, sB);
256 kernels::internal::LoadMatrix(TR_D1D, Q1D, G, sG);
257
258 kernels::internal::LoadDofs3d(e, TR_D1D, X, g0);
259 kernels::internal::Grad3d(TR_D1D, Q1D, smem, sB, sG, g0, g1);
260
261 for (int qz = 0; qz < Q1D; qz++)
262 {
263 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
264 {
265 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
266 {
267 r0[0][qz][qy][qx] =
268 // c = 0
269 g1[0][0][qz][qy][qx] * Q(qx, qy, qz, 0, 0, e) +
270 g1[0][1][qz][qy][qx] * Q(qx, qy, qz, 1, 0, e) +
271 g1[0][2][qz][qy][qx] * Q(qx, qy, qz, 2, 0, e) +
272 // c = 1
273 g1[1][0][qz][qy][qx] * Q(qx, qy, qz, 0, 1, e) +
274 g1[1][1][qz][qy][qx] * Q(qx, qy, qz, 1, 1, e) +
275 g1[1][2][qz][qy][qx] * Q(qx, qy, qz, 2, 1, e) +
276 // c = 2
277 g1[2][0][qz][qy][qx] * Q(qx, qy, qz, 0, 2, e) +
278 g1[2][1][qz][qy][qx] * Q(qx, qy, qz, 1, 2, e) +
279 g1[2][2][qz][qy][qx] * Q(qx, qy, qz, 2, 2, e);
280 }
281 }
282 }
283 MFEM_SYNC_THREAD;
284
285 kernels::internal::LoadMatrix<MQ1, true>(TE_D1D, Q1D, Bt, sB);
286 kernels::internal::EvalTranspose3d(TE_D1D, Q1D, smem, sB, r0, r1);
287 kernels::internal::WriteDofs3d(e, TE_D1D, r1, Y);
288 });
289}
290
291} // namespace internal
292
293template<int DIM, int T_TR_D1D, int T_TE_D1D, int T_Q1D>
294VectorDivergenceIntegrator::VectorDivergenceAddMultPAType
295VectorDivergenceIntegrator::VectorDivergenceAddMultPA::Kernel()
296{
297 static_assert(T_TR_D1D <= T_Q1D && T_TE_D1D <= T_Q1D);
298 if constexpr (DIM == 2)
299 {
300 return internal::SmemPADivergenceApply2D<T_TR_D1D, T_TE_D1D, T_Q1D>;
301 }
302 else if constexpr (DIM == 3)
303 {
304 return internal::SmemPADivergenceApply3D<T_TR_D1D, T_TE_D1D, T_Q1D>;
305 }
306 MFEM_ABORT("Unsupported kernel");
307}
308
309inline VectorDivergenceIntegrator::VectorDivergenceAddMultPAType
310VectorDivergenceIntegrator::VectorDivergenceAddMultPA::Fallback
311(int dim, int tr_d1d, int te_d1d, int q1d)
312{
313 MFEM_VERIFY(tr_d1d <= q1d && te_d1d <= q1d, "");
314 MFEM_VERIFY(tr_d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
315 MFEM_VERIFY(te_d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
316 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
317 if (dim == 2)
318 {
319 return internal::SmemPADivergenceApply2D;
320 }
321 else if (dim == 3)
322 {
323 return internal::SmemPADivergenceApply3D;
324 }
325 MFEM_ABORT("Unsupported kernel");
326}
327
328template<int DIM, int T_TR_D1D, int T_TE_D1D, int T_Q1D>
329VectorDivergenceIntegrator::VectorDivergenceAddMultTransposePAType
330VectorDivergenceIntegrator::VectorDivergenceAddMultTransposePA::Kernel()
331{
332 static_assert(T_TR_D1D <= T_Q1D && T_TE_D1D <= T_Q1D);
333 if constexpr (DIM == 2)
334 {
335 return internal::SmemPADivergenceApplyTranspose2D<T_TR_D1D, T_TE_D1D, T_Q1D>;
336 }
337 else if constexpr (DIM == 3)
338 {
339 return internal::SmemPADivergenceApplyTranspose3D<T_TR_D1D, T_TE_D1D, T_Q1D>;
340 }
341 MFEM_ABORT("Unsupported kernel");
342}
343
344inline VectorDivergenceIntegrator::VectorDivergenceAddMultTransposePAType
345VectorDivergenceIntegrator::VectorDivergenceAddMultTransposePA::Fallback
346(int dim, int tr_d1d, int te_d1d, int q1d)
347{
348 MFEM_VERIFY(tr_d1d <= q1d && te_d1d <= q1d, "");
349 MFEM_VERIFY(tr_d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
350 MFEM_VERIFY(te_d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
351 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
352 if (dim == 2)
353 {
354 return internal::SmemPADivergenceApplyTranspose2D;
355 }
356 else if (dim == 3)
357 {
358 return internal::SmemPADivergenceApplyTranspose3D;
359 }
360 MFEM_ABORT("Unsupported kernel");
361}
362
363/// \endcond DO_NOT_DOCUMENT
364
365} // namespace mfem
mfem::VectorDivergenceIntegrator::VectorDivergenceAddMultTransposePAType
void(*)(const int ne, const Array< real_t > &bt, const Array< real_t > &gt, const Array< real_t > &b, const Vector &q, const Vector &x, Vector &y, const int tr_d1d, const int te_d1d, const int q1d) VectorDivergenceAddMultTransposePAType
Definition bilininteg.hpp:3206
mfem::VectorDivergenceIntegrator::VectorDivergenceAddMultPAType
void(*)(const int ne, const Array< real_t > &b, const Array< real_t > &g, const Array< real_t > &bt, const Vector &op, const Vector &x, Vector &y, const int tr_d1d, const int te_d1d, const int q1d) VectorDivergenceAddMultPAType
Definition bilininteg.hpp:3197
dim
int dim
Definition ex24.cpp:53
b
real_t b
Definition lissajous.cpp:42
DIM
constexpr int DIM
Definition minimal-surface.cpp:71
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
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::DeviceDofQuadLimits::Get
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:138
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0