MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_dgtrace_kernels.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
12#ifndef BILININTEG_DGTRACE_KERNELS_HPP
13#define BILININTEG_DGTRACE_KERNELS_HPP
14
15#include "../../general/forall.hpp"
16#include "../bilininteg.hpp"
17#include "../gridfunc.hpp"
18#include "../qfunction.hpp"
19#include "../restriction.hpp"
20
21/// \cond DO_NOT_DOCUMENT
22namespace mfem
23{
24
25namespace internal
26{
27
28// PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
29template <int T_D1D = 0, int T_Q1D = 0>
30void PADGTraceApply2D(const int NF, const Array<real_t> &b,
31 const Array<real_t> &bt, const Vector &op_,
32 const Vector &x_, Vector &y_, const int d1d = 0,
33 const int q1d = 0)
34{
35 const int VDIM = 1;
36 const int D1D = T_D1D ? T_D1D : d1d;
37 const int Q1D = T_Q1D ? T_Q1D : q1d;
38 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
39 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
40 auto B = Reshape(b.Read(), Q1D, D1D);
41 auto Bt = Reshape(bt.Read(), D1D, Q1D);
42 auto op = Reshape(op_.Read(), Q1D, 2, 2, NF);
43 auto x = Reshape(x_.Read(), D1D, VDIM, 2, NF);
44 auto y = Reshape(y_.ReadWrite(), D1D, VDIM, 2, NF);
45
46 mfem::forall(NF, [=] MFEM_HOST_DEVICE(int f)
47 {
48 const int VDIM = 1;
49 const int D1D = T_D1D ? T_D1D : d1d;
50 const int Q1D = T_Q1D ? T_Q1D : q1d;
51 // the following variables are evaluated at compile time
52 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
53 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
54 real_t u0[max_D1D][VDIM];
55 real_t u1[max_D1D][VDIM];
56 for (int d = 0; d < D1D; d++)
57 {
58 for (int c = 0; c < VDIM; c++)
59 {
60 u0[d][c] = x(d, c, 0, f);
61 u1[d][c] = x(d, c, 1, f);
62 }
63 }
64 real_t Bu0[max_Q1D][VDIM];
65 real_t Bu1[max_Q1D][VDIM];
66 for (int q = 0; q < Q1D; ++q)
67 {
68 for (int c = 0; c < VDIM; c++)
69 {
70 Bu0[q][c] = 0.0;
71 Bu1[q][c] = 0.0;
72 }
73 for (int d = 0; d < D1D; ++d)
74 {
75 const real_t b = B(q, d);
76 for (int c = 0; c < VDIM; c++)
77 {
78 Bu0[q][c] += b * u0[d][c];
79 Bu1[q][c] += b * u1[d][c];
80 }
81 }
82 }
83 real_t DBu[max_Q1D][VDIM];
84 for (int q = 0; q < Q1D; ++q)
85 {
86 for (int c = 0; c < VDIM; c++)
87 {
88 DBu[q][c] = op(q, 0, 0, f) * Bu0[q][c] + op(q, 1, 0, f) * Bu1[q][c];
89 }
90 }
91 real_t BDBu[max_D1D][VDIM];
92 for (int d = 0; d < D1D; ++d)
93 {
94 for (int c = 0; c < VDIM; c++)
95 {
96 BDBu[d][c] = 0.0;
97 }
98 for (int q = 0; q < Q1D; ++q)
99 {
100 const real_t b = Bt(d, q);
101 for (int c = 0; c < VDIM; c++)
102 {
103 BDBu[d][c] += b * DBu[q][c];
104 }
105 }
106 for (int c = 0; c < VDIM; c++)
107 {
108 y(d, c, 0, f) += BDBu[d][c];
109 y(d, c, 1, f) += -BDBu[d][c];
110 }
111 }
112 });
113}
114
115// PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
116template <int T_D1D = 0, int T_Q1D = 0>
117void PADGTraceApply3D(const int NF, const Array<real_t> &b,
118 const Array<real_t> &bt, const Vector &op_,
119 const Vector &x_, Vector &y_, const int d1d = 0,
120 const int q1d = 0)
121{
122 const int VDIM = 1;
123 const int D1D = T_D1D ? T_D1D : d1d;
124 const int Q1D = T_Q1D ? T_Q1D : q1d;
125 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
126 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
127 auto B = Reshape(b.Read(), Q1D, D1D);
128 auto Bt = Reshape(bt.Read(), D1D, Q1D);
129 auto op = Reshape(op_.Read(), Q1D, Q1D, 2, 2, NF);
130 auto x = Reshape(x_.Read(), D1D, D1D, VDIM, 2, NF);
131 auto y = Reshape(y_.ReadWrite(), D1D, D1D, VDIM, 2, NF);
132
133 mfem::forall(NF, [=] MFEM_HOST_DEVICE(int f)
134 {
135 const int VDIM = 1;
136 const int D1D = T_D1D ? T_D1D : d1d;
137 const int Q1D = T_Q1D ? T_Q1D : q1d;
138 // the following variables are evaluated at compile time
139 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
140 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
141 real_t u0[max_D1D][max_D1D][VDIM];
142 real_t u1[max_D1D][max_D1D][VDIM];
143 for (int d1 = 0; d1 < D1D; d1++)
144 {
145 for (int d2 = 0; d2 < D1D; d2++)
146 {
147 for (int c = 0; c < VDIM; c++)
148 {
149 u0[d1][d2][c] = x(d1, d2, c, 0, f);
150 u1[d1][d2][c] = x(d1, d2, c, 1, f);
151 }
152 }
153 }
154 real_t Bu0[max_Q1D][max_D1D][VDIM];
155 real_t Bu1[max_Q1D][max_D1D][VDIM];
156 for (int q = 0; q < Q1D; ++q)
157 {
158 for (int d2 = 0; d2 < D1D; d2++)
159 {
160 for (int c = 0; c < VDIM; c++)
161 {
162 Bu0[q][d2][c] = 0.0;
163 Bu1[q][d2][c] = 0.0;
164 }
165 for (int d1 = 0; d1 < D1D; ++d1)
166 {
167 const real_t b = B(q, d1);
168 for (int c = 0; c < VDIM; c++)
169 {
170 Bu0[q][d2][c] += b * u0[d1][d2][c];
171 Bu1[q][d2][c] += b * u1[d1][d2][c];
172 }
173 }
174 }
175 }
176 real_t BBu0[max_Q1D][max_Q1D][VDIM];
177 real_t BBu1[max_Q1D][max_Q1D][VDIM];
178 for (int q1 = 0; q1 < Q1D; ++q1)
179 {
180 for (int q2 = 0; q2 < Q1D; q2++)
181 {
182 for (int c = 0; c < VDIM; c++)
183 {
184 BBu0[q1][q2][c] = 0.0;
185 BBu1[q1][q2][c] = 0.0;
186 }
187 for (int d2 = 0; d2 < D1D; ++d2)
188 {
189 const real_t b = B(q2, d2);
190 for (int c = 0; c < VDIM; c++)
191 {
192 BBu0[q1][q2][c] += b * Bu0[q1][d2][c];
193 BBu1[q1][q2][c] += b * Bu1[q1][d2][c];
194 }
195 }
196 }
197 }
198 real_t DBBu[max_Q1D][max_Q1D][VDIM];
199 for (int q1 = 0; q1 < Q1D; ++q1)
200 {
201 for (int q2 = 0; q2 < Q1D; q2++)
202 {
203 for (int c = 0; c < VDIM; c++)
204 {
205 DBBu[q1][q2][c] = op(q1, q2, 0, 0, f) * BBu0[q1][q2][c] +
206 op(q1, q2, 1, 0, f) * BBu1[q1][q2][c];
207 }
208 }
209 }
210 real_t BDBBu[max_Q1D][max_D1D][VDIM];
211 for (int q1 = 0; q1 < Q1D; ++q1)
212 {
213 for (int d2 = 0; d2 < D1D; d2++)
214 {
215 for (int c = 0; c < VDIM; c++)
216 {
217 BDBBu[q1][d2][c] = 0.0;
218 }
219 for (int q2 = 0; q2 < Q1D; ++q2)
220 {
221 const real_t b = Bt(d2, q2);
222 for (int c = 0; c < VDIM; c++)
223 {
224 BDBBu[q1][d2][c] += b * DBBu[q1][q2][c];
225 }
226 }
227 }
228 }
229 real_t BBDBBu[max_D1D][max_D1D][VDIM];
230 for (int d1 = 0; d1 < D1D; ++d1)
231 {
232 for (int d2 = 0; d2 < D1D; d2++)
233 {
234 for (int c = 0; c < VDIM; c++)
235 {
236 BBDBBu[d1][d2][c] = 0.0;
237 }
238 for (int q1 = 0; q1 < Q1D; ++q1)
239 {
240 const real_t b = Bt(d1, q1);
241 for (int c = 0; c < VDIM; c++)
242 {
243 BBDBBu[d1][d2][c] += b * BDBBu[q1][d2][c];
244 }
245 }
246 for (int c = 0; c < VDIM; c++)
247 {
248 y(d1, d2, c, 0, f) += BBDBBu[d1][d2][c];
249 y(d1, d2, c, 1, f) += -BBDBBu[d1][d2][c];
250 }
251 }
252 }
253 });
254}
255
256// Optimized PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
257template <int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
258void SmemPADGTraceApply3D(const int NF, const Array<real_t> &b,
259 const Array<real_t> &bt, const Vector &op_,
260 const Vector &x_, Vector &y_, const int d1d = 0,
261 const int q1d = 0)
262{
263 const int D1D = T_D1D ? T_D1D : d1d;
264 const int Q1D = T_Q1D ? T_Q1D : q1d;
265 constexpr int NBZ = T_NBZ ? T_NBZ : 1;
266 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
267 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
268 auto B = Reshape(b.Read(), Q1D, D1D);
269 auto Bt = Reshape(bt.Read(), D1D, Q1D);
270 auto op = Reshape(op_.Read(), Q1D, Q1D, 2, 2, NF);
271 auto x = Reshape(x_.Read(), D1D, D1D, 2, NF);
272 auto y = Reshape(y_.ReadWrite(), D1D, D1D, 2, NF);
273
274 mfem::forall_2D_batch(NF, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE(int f)
275 {
276 const int tidz = MFEM_THREAD_ID(z);
277 const int D1D = T_D1D ? T_D1D : d1d;
278 const int Q1D = T_Q1D ? T_Q1D : q1d;
279 // the following variables are evaluated at compile time
280 constexpr int NBZ = T_NBZ ? T_NBZ : 1;
281 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
282 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
283 MFEM_SHARED real_t u0[NBZ][max_D1D][max_D1D];
284 MFEM_SHARED real_t u1[NBZ][max_D1D][max_D1D];
285 MFEM_FOREACH_THREAD(d1, x, D1D)
286 {
287 MFEM_FOREACH_THREAD(d2, y, D1D)
288 {
289 u0[tidz][d1][d2] = x(d1, d2, 0, f);
290 u1[tidz][d1][d2] = x(d1, d2, 1, f);
291 }
292 }
293 MFEM_SYNC_THREAD;
294 MFEM_SHARED real_t Bu0[NBZ][max_Q1D][max_D1D];
295 MFEM_SHARED real_t Bu1[NBZ][max_Q1D][max_D1D];
296 MFEM_FOREACH_THREAD(q1, x, Q1D)
297 {
298 MFEM_FOREACH_THREAD(d2, y, D1D)
299 {
300 real_t Bu0_ = 0.0;
301 real_t Bu1_ = 0.0;
302 for (int d1 = 0; d1 < D1D; ++d1)
303 {
304 const real_t b = B(q1, d1);
305 Bu0_ += b * u0[tidz][d1][d2];
306 Bu1_ += b * u1[tidz][d1][d2];
307 }
308 Bu0[tidz][q1][d2] = Bu0_;
309 Bu1[tidz][q1][d2] = Bu1_;
310 }
311 }
312 MFEM_SYNC_THREAD;
313 MFEM_SHARED real_t BBu0[NBZ][max_Q1D][max_Q1D];
314 MFEM_SHARED real_t BBu1[NBZ][max_Q1D][max_Q1D];
315 MFEM_FOREACH_THREAD(q1, x, Q1D)
316 {
317 MFEM_FOREACH_THREAD(q2, y, Q1D)
318 {
319 real_t BBu0_ = 0.0;
320 real_t BBu1_ = 0.0;
321 for (int d2 = 0; d2 < D1D; ++d2)
322 {
323 const real_t b = B(q2, d2);
324 BBu0_ += b * Bu0[tidz][q1][d2];
325 BBu1_ += b * Bu1[tidz][q1][d2];
326 }
327 BBu0[tidz][q1][q2] = BBu0_;
328 BBu1[tidz][q1][q2] = BBu1_;
329 }
330 }
331 MFEM_SYNC_THREAD;
332 MFEM_SHARED real_t DBBu[NBZ][max_Q1D][max_Q1D];
333 MFEM_FOREACH_THREAD(q1, x, Q1D)
334 {
335 MFEM_FOREACH_THREAD(q2, y, Q1D)
336 {
337 DBBu[tidz][q1][q2] = op(q1, q2, 0, 0, f) * BBu0[tidz][q1][q2] +
338 op(q1, q2, 1, 0, f) * BBu1[tidz][q1][q2];
339 }
340 }
341 MFEM_SYNC_THREAD;
342 MFEM_SHARED real_t BDBBu[NBZ][max_Q1D][max_D1D];
343 MFEM_FOREACH_THREAD(q1, x, Q1D)
344 {
345 MFEM_FOREACH_THREAD(d2, y, D1D)
346 {
347 real_t BDBBu_ = 0.0;
348 for (int q2 = 0; q2 < Q1D; ++q2)
349 {
350 const real_t b = Bt(d2, q2);
351 BDBBu_ += b * DBBu[tidz][q1][q2];
352 }
353 BDBBu[tidz][q1][d2] = BDBBu_;
354 }
355 }
356 MFEM_SYNC_THREAD;
357 MFEM_FOREACH_THREAD(d1, x, D1D)
358 {
359 MFEM_FOREACH_THREAD(d2, y, D1D)
360 {
361 real_t BBDBBu_ = 0.0;
362 for (int q1 = 0; q1 < Q1D; ++q1)
363 {
364 const real_t b = Bt(d1, q1);
365 BBDBBu_ += b * BDBBu[tidz][q1][d2];
366 }
367 y(d1, d2, 0, f) += BBDBBu_;
368 y(d1, d2, 1, f) += -BBDBBu_;
369 }
370 }
371 });
372}
373
374// PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
375template <int T_D1D = 0, int T_Q1D = 0>
376void PADGTraceApplyTranspose2D(const int NF, const Array<real_t> &b,
377 const Array<real_t> &bt, const Vector &op_,
378 const Vector &x_, Vector &y_, const int d1d = 0,
379 const int q1d = 0)
380{
381 const int VDIM = 1;
382 const int D1D = T_D1D ? T_D1D : d1d;
383 const int Q1D = T_Q1D ? T_Q1D : q1d;
384 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
385 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
386 auto B = Reshape(b.Read(), Q1D, D1D);
387 auto Bt = Reshape(bt.Read(), D1D, Q1D);
388 auto op = Reshape(op_.Read(), Q1D, 2, 2, NF);
389 auto x = Reshape(x_.Read(), D1D, VDIM, 2, NF);
390 auto y = Reshape(y_.ReadWrite(), D1D, VDIM, 2, NF);
391
392 mfem::forall(NF, [=] MFEM_HOST_DEVICE(int f)
393 {
394 const int VDIM = 1;
395 const int D1D = T_D1D ? T_D1D : d1d;
396 const int Q1D = T_Q1D ? T_Q1D : q1d;
397 // the following variables are evaluated at compile time
398 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
399 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
400 real_t u0[max_D1D][VDIM];
401 real_t u1[max_D1D][VDIM];
402 for (int d = 0; d < D1D; d++)
403 {
404 for (int c = 0; c < VDIM; c++)
405 {
406 u0[d][c] = x(d, c, 0, f);
407 u1[d][c] = x(d, c, 1, f);
408 }
409 }
410 real_t Bu0[max_Q1D][VDIM];
411 real_t Bu1[max_Q1D][VDIM];
412 for (int q = 0; q < Q1D; ++q)
413 {
414 for (int c = 0; c < VDIM; c++)
415 {
416 Bu0[q][c] = 0.0;
417 Bu1[q][c] = 0.0;
418 }
419 for (int d = 0; d < D1D; ++d)
420 {
421 const real_t b = B(q, d);
422 for (int c = 0; c < VDIM; c++)
423 {
424 Bu0[q][c] += b * u0[d][c];
425 Bu1[q][c] += b * u1[d][c];
426 }
427 }
428 }
429 real_t DBu0[max_Q1D][VDIM];
430 real_t DBu1[max_Q1D][VDIM];
431 for (int q = 0; q < Q1D; ++q)
432 {
433 for (int c = 0; c < VDIM; c++)
434 {
435 DBu0[q][c] =
436 op(q, 0, 0, f) * Bu0[q][c] + op(q, 0, 1, f) * Bu1[q][c];
437 DBu1[q][c] =
438 op(q, 1, 0, f) * Bu0[q][c] + op(q, 1, 1, f) * Bu1[q][c];
439 }
440 }
441 real_t BDBu0[max_D1D][VDIM];
442 real_t BDBu1[max_D1D][VDIM];
443 for (int d = 0; d < D1D; ++d)
444 {
445 for (int c = 0; c < VDIM; c++)
446 {
447 BDBu0[d][c] = 0.0;
448 BDBu1[d][c] = 0.0;
449 }
450 for (int q = 0; q < Q1D; ++q)
451 {
452 const real_t b = Bt(d, q);
453 for (int c = 0; c < VDIM; c++)
454 {
455 BDBu0[d][c] += b * DBu0[q][c];
456 BDBu1[d][c] += b * DBu1[q][c];
457 }
458 }
459 for (int c = 0; c < VDIM; c++)
460 {
461 y(d, c, 0, f) += BDBu0[d][c];
462 y(d, c, 1, f) += BDBu1[d][c];
463 }
464 }
465 });
466}
467
468// PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
469template <int T_D1D = 0, int T_Q1D = 0>
470void PADGTraceApplyTranspose3D(const int NF, const Array<real_t> &b,
471 const Array<real_t> &bt, const Vector &op_,
472 const Vector &x_, Vector &y_, const int d1d = 0,
473 const int q1d = 0)
474{
475 const int VDIM = 1;
476 const int D1D = T_D1D ? T_D1D : d1d;
477 const int Q1D = T_Q1D ? T_Q1D : q1d;
478 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
479 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
480 auto B = Reshape(b.Read(), Q1D, D1D);
481 auto Bt = Reshape(bt.Read(), D1D, Q1D);
482 auto op = Reshape(op_.Read(), Q1D, Q1D, 2, 2, NF);
483 auto x = Reshape(x_.Read(), D1D, D1D, VDIM, 2, NF);
484 auto y = Reshape(y_.ReadWrite(), D1D, D1D, VDIM, 2, NF);
485
486 mfem::forall(NF, [=] MFEM_HOST_DEVICE(int f)
487 {
488 const int VDIM = 1;
489 const int D1D = T_D1D ? T_D1D : d1d;
490 const int Q1D = T_Q1D ? T_Q1D : q1d;
491 // the following variables are evaluated at compile time
492 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
493 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
494 real_t u0[max_D1D][max_D1D][VDIM];
495 real_t u1[max_D1D][max_D1D][VDIM];
496 for (int d1 = 0; d1 < D1D; d1++)
497 {
498 for (int d2 = 0; d2 < D1D; d2++)
499 {
500 for (int c = 0; c < VDIM; c++)
501 {
502 u0[d1][d2][c] = x(d1, d2, c, 0, f);
503 u1[d1][d2][c] = x(d1, d2, c, 1, f);
504 }
505 }
506 }
507 real_t Bu0[max_Q1D][max_D1D][VDIM];
508 real_t Bu1[max_Q1D][max_D1D][VDIM];
509 for (int q1 = 0; q1 < Q1D; ++q1)
510 {
511 for (int d2 = 0; d2 < D1D; ++d2)
512 {
513 for (int c = 0; c < VDIM; c++)
514 {
515 Bu0[q1][d2][c] = 0.0;
516 Bu1[q1][d2][c] = 0.0;
517 }
518 for (int d1 = 0; d1 < D1D; ++d1)
519 {
520 const real_t b = B(q1, d1);
521 for (int c = 0; c < VDIM; c++)
522 {
523 Bu0[q1][d2][c] += b * u0[d1][d2][c];
524 Bu1[q1][d2][c] += b * u1[d1][d2][c];
525 }
526 }
527 }
528 }
529 real_t BBu0[max_Q1D][max_Q1D][VDIM];
530 real_t BBu1[max_Q1D][max_Q1D][VDIM];
531 for (int q1 = 0; q1 < Q1D; ++q1)
532 {
533 for (int q2 = 0; q2 < Q1D; ++q2)
534 {
535 for (int c = 0; c < VDIM; c++)
536 {
537 BBu0[q1][q2][c] = 0.0;
538 BBu1[q1][q2][c] = 0.0;
539 }
540 for (int d2 = 0; d2 < D1D; ++d2)
541 {
542 const real_t b = B(q2, d2);
543 for (int c = 0; c < VDIM; c++)
544 {
545 BBu0[q1][q2][c] += b * Bu0[q1][d2][c];
546 BBu1[q1][q2][c] += b * Bu1[q1][d2][c];
547 }
548 }
549 }
550 }
551 real_t DBu0[max_Q1D][max_Q1D][VDIM];
552 real_t DBu1[max_Q1D][max_Q1D][VDIM];
553 for (int q1 = 0; q1 < Q1D; ++q1)
554 {
555 for (int q2 = 0; q2 < Q1D; ++q2)
556 {
557 const real_t D00 = op(q1, q2, 0, 0, f);
558 const real_t D01 = op(q1, q2, 0, 1, f);
559 const real_t D10 = op(q1, q2, 1, 0, f);
560 const real_t D11 = op(q1, q2, 1, 1, f);
561 for (int c = 0; c < VDIM; c++)
562 {
563 DBu0[q1][q2][c] = D00 * BBu0[q1][q2][c] + D01 * BBu1[q1][q2][c];
564 DBu1[q1][q2][c] = D10 * BBu0[q1][q2][c] + D11 * BBu1[q1][q2][c];
565 }
566 }
567 }
568 real_t BDBu0[max_D1D][max_Q1D][VDIM];
569 real_t BDBu1[max_D1D][max_Q1D][VDIM];
570 for (int d1 = 0; d1 < D1D; ++d1)
571 {
572 for (int q2 = 0; q2 < Q1D; ++q2)
573 {
574 for (int c = 0; c < VDIM; c++)
575 {
576 BDBu0[d1][q2][c] = 0.0;
577 BDBu1[d1][q2][c] = 0.0;
578 }
579 for (int q1 = 0; q1 < Q1D; ++q1)
580 {
581 const real_t b = Bt(d1, q1);
582 for (int c = 0; c < VDIM; c++)
583 {
584 BDBu0[d1][q2][c] += b * DBu0[q1][q2][c];
585 BDBu1[d1][q2][c] += b * DBu1[q1][q2][c];
586 }
587 }
588 }
589 }
590 real_t BBDBu0[max_D1D][max_D1D][VDIM];
591 real_t BBDBu1[max_D1D][max_D1D][VDIM];
592 for (int d1 = 0; d1 < D1D; ++d1)
593 {
594 for (int d2 = 0; d2 < D1D; ++d2)
595 {
596 for (int c = 0; c < VDIM; c++)
597 {
598 BBDBu0[d1][d2][c] = 0.0;
599 BBDBu1[d1][d2][c] = 0.0;
600 }
601 for (int q2 = 0; q2 < Q1D; ++q2)
602 {
603 const real_t b = Bt(d2, q2);
604 for (int c = 0; c < VDIM; c++)
605 {
606 BBDBu0[d1][d2][c] += b * BDBu0[d1][q2][c];
607 BBDBu1[d1][d2][c] += b * BDBu1[d1][q2][c];
608 }
609 }
610 for (int c = 0; c < VDIM; c++)
611 {
612 y(d1, d2, c, 0, f) += BBDBu0[d1][d2][c];
613 y(d1, d2, c, 1, f) += BBDBu1[d1][d2][c];
614 }
615 }
616 }
617 });
618}
619
620// Optimized PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
621template <int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
622void SmemPADGTraceApplyTranspose3D(const int NF, const Array<real_t> &b,
623 const Array<real_t> &bt, const Vector &op_,
624 const Vector &x_, Vector &y_,
625 const int d1d = 0, const int q1d = 0)
626{
627 const int D1D = T_D1D ? T_D1D : d1d;
628 const int Q1D = T_Q1D ? T_Q1D : q1d;
629 constexpr int NBZ = T_NBZ ? T_NBZ : 1;
630 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
631 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
632 auto B = Reshape(b.Read(), Q1D, D1D);
633 auto Bt = Reshape(bt.Read(), D1D, Q1D);
634 auto op = Reshape(op_.Read(), Q1D, Q1D, 2, 2, NF);
635 auto x = Reshape(x_.Read(), D1D, D1D, 2, NF);
636 auto y = Reshape(y_.ReadWrite(), D1D, D1D, 2, NF);
637
638 mfem::forall_2D_batch(NF, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE(int f)
639 {
640 const int tidz = MFEM_THREAD_ID(z);
641 const int D1D = T_D1D ? T_D1D : d1d;
642 const int Q1D = T_Q1D ? T_Q1D : q1d;
643 // the following variables are evaluated at compile time
644 constexpr int NBZ = T_NBZ ? T_NBZ : 1;
645 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
646 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
647 MFEM_SHARED real_t u0[NBZ][max_D1D][max_D1D];
648 MFEM_SHARED real_t u1[NBZ][max_D1D][max_D1D];
649 MFEM_FOREACH_THREAD(d1, x, D1D)
650 {
651 MFEM_FOREACH_THREAD(d2, y, D1D)
652 {
653 u0[tidz][d1][d2] = x(d1, d2, 0, f);
654 u1[tidz][d1][d2] = x(d1, d2, 1, f);
655 }
656 }
657 MFEM_SYNC_THREAD;
658 MFEM_SHARED real_t Bu0[NBZ][max_Q1D][max_D1D];
659 MFEM_SHARED real_t Bu1[NBZ][max_Q1D][max_D1D];
660 MFEM_FOREACH_THREAD(q1, x, Q1D)
661 {
662 MFEM_FOREACH_THREAD(d2, y, D1D)
663 {
664 real_t Bu0_ = 0.0;
665 real_t Bu1_ = 0.0;
666 for (int d1 = 0; d1 < D1D; ++d1)
667 {
668 const real_t b = B(q1, d1);
669 Bu0_ += b * u0[tidz][d1][d2];
670 Bu1_ += b * u1[tidz][d1][d2];
671 }
672 Bu0[tidz][q1][d2] = Bu0_;
673 Bu1[tidz][q1][d2] = Bu1_;
674 }
675 }
676 MFEM_SYNC_THREAD;
677 MFEM_SHARED real_t BBu0[NBZ][max_Q1D][max_Q1D];
678 MFEM_SHARED real_t BBu1[NBZ][max_Q1D][max_Q1D];
679 MFEM_FOREACH_THREAD(q1, x, Q1D)
680 {
681 MFEM_FOREACH_THREAD(q2, y, Q1D)
682 {
683 real_t BBu0_ = 0.0;
684 real_t BBu1_ = 0.0;
685 for (int d2 = 0; d2 < D1D; ++d2)
686 {
687 const real_t b = B(q2, d2);
688 BBu0_ += b * Bu0[tidz][q1][d2];
689 BBu1_ += b * Bu1[tidz][q1][d2];
690 }
691 BBu0[tidz][q1][q2] = BBu0_;
692 BBu1[tidz][q1][q2] = BBu1_;
693 }
694 }
695 MFEM_SYNC_THREAD;
696 MFEM_SHARED real_t DBBu0[NBZ][max_Q1D][max_Q1D];
697 MFEM_SHARED real_t DBBu1[NBZ][max_Q1D][max_Q1D];
698 MFEM_FOREACH_THREAD(q1, x, Q1D)
699 {
700 MFEM_FOREACH_THREAD(q2, y, Q1D)
701 {
702 const real_t D00 = op(q1, q2, 0, 0, f);
703 const real_t D01 = op(q1, q2, 0, 1, f);
704 const real_t D10 = op(q1, q2, 1, 0, f);
705 const real_t D11 = op(q1, q2, 1, 1, f);
706 const real_t u0q = BBu0[tidz][q1][q2];
707 const real_t u1q = BBu1[tidz][q1][q2];
708 DBBu0[tidz][q1][q2] = D00 * u0q + D01 * u1q;
709 DBBu1[tidz][q1][q2] = D10 * u0q + D11 * u1q;
710 }
711 }
712 MFEM_SYNC_THREAD;
713 MFEM_SHARED real_t BDBBu0[NBZ][max_Q1D][max_D1D];
714 MFEM_SHARED real_t BDBBu1[NBZ][max_Q1D][max_D1D];
715 MFEM_FOREACH_THREAD(q1, x, Q1D)
716 {
717 MFEM_FOREACH_THREAD(d2, y, D1D)
718 {
719 real_t BDBBu0_ = 0.0;
720 real_t BDBBu1_ = 0.0;
721 for (int q2 = 0; q2 < Q1D; ++q2)
722 {
723 const real_t b = Bt(d2, q2);
724 BDBBu0_ += b * DBBu0[tidz][q1][q2];
725 BDBBu1_ += b * DBBu1[tidz][q1][q2];
726 }
727 BDBBu0[tidz][q1][d2] = BDBBu0_;
728 BDBBu1[tidz][q1][d2] = BDBBu1_;
729 }
730 }
731 MFEM_SYNC_THREAD;
732 MFEM_FOREACH_THREAD(d1, x, D1D)
733 {
734 MFEM_FOREACH_THREAD(d2, y, D1D)
735 {
736 real_t BBDBBu0_ = 0.0;
737 real_t BBDBBu1_ = 0.0;
738 for (int q1 = 0; q1 < Q1D; ++q1)
739 {
740 const real_t b = Bt(d1, q1);
741 BBDBBu0_ += b * BDBBu0[tidz][q1][d2];
742 BBDBBu1_ += b * BDBBu1[tidz][q1][d2];
743 }
744 y(d1, d2, 0, f) += BBDBBu0_;
745 y(d1, d2, 1, f) += BBDBBu1_;
746 }
747 }
748 });
749}
750
751} // namespace internal
752
753template <int DIM, int D1D, int Q1D>
754DGTraceIntegrator::ApplyKernelType DGTraceIntegrator::ApplyPAKernels::Kernel()
755{
756 if constexpr (DIM == 2)
757 {
758 return internal::PADGTraceApply2D<D1D, Q1D>;
759 }
760 else if constexpr (DIM == 3)
761 {
762 if constexpr (D1D == 3 || D1D == 4)
763 {
764 return internal::SmemPADGTraceApply3D<D1D, Q1D, 2>;
765 }
766 else
767 {
768 return internal::SmemPADGTraceApply3D<D1D, Q1D>;
769 }
770 }
771 MFEM_ABORT("");
772}
773
774template <int DIM, int D1D, int Q1D>
775DGTraceIntegrator::ApplyKernelType DGTraceIntegrator::ApplyPATKernels::Kernel()
776{
777 if constexpr (DIM == 2)
778 {
779 return internal::PADGTraceApplyTranspose2D<D1D, Q1D>;
780 }
781 else if constexpr (DIM == 3)
782 {
783 return internal::SmemPADGTraceApplyTranspose3D<D1D, Q1D>;
784 }
785 MFEM_ABORT("");
786}
787} // namespace mfem
788
789/// \endcond DO_NOT_DOCUMENT
790#endif
mfem::DGTraceIntegrator::ApplyKernelType
void(*)(const int, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, Vector &, const int, const int) ApplyKernelType
arguments: nf, B, Bt, pa_data, x, y, dofs1D, quad1D
Definition bilininteg.hpp:3593
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_batch
void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
Definition forall.hpp:1232
mfem::f
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:1134
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