MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_vecmass_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
21using mfem::kernels::internal::SetMaxOf;
22
23namespace mfem
24{
25
26/// \cond DO_NOT_DOCUMENT
27
28namespace internal
29{
30
31template <int T_D1D = 0, int T_Q1D = 0>
32void SmemPAVectorMassApply2D(const int NE,
33 const int coeff_vdim,
34 const Array<real_t> &b,
35 const Vector &d,
36 const Vector &x,
37 Vector &y,
38 const int d1d = 0,
39 const int q1d = 0)
40{
41 static constexpr int DIM = 2, VDIM = 2;
42 const int D1D = T_D1D ? T_D1D : d1d;
43 const int Q1D = T_Q1D ? T_Q1D : q1d;
44
45 const bool const_coeff = coeff_vdim == 1;
46 const bool vector_coeff = coeff_vdim == DIM;
47 const bool matrix_coeff = coeff_vdim == DIM*DIM;
48
49 const auto B = b.Read();
50 const auto D = Reshape(d.Read(), Q1D, Q1D, coeff_vdim, NE);
51 const auto X = Reshape(x.Read(), D1D, D1D, VDIM, NE);
52 auto Y = Reshape(y.ReadWrite(), D1D, D1D, VDIM, NE);
53
54 mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
55 {
56 constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
57 constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
58
59 MFEM_SHARED real_t sB[MD1][MQ1], smem[MQ1][MQ1];
60 kernels::internal::v_regs2d_t<VDIM, MQ1> r0, r1;
61 kernels::internal::LoadMatrix(D1D, Q1D, B, sB);
62 kernels::internal::LoadDofs2d(e, D1D, X, r0);
63 kernels::internal::Eval2d(D1D, Q1D, smem, sB, r0, r1);
64
65 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
66 {
67 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
68 {
69 const real_t Qx = r1[0][qy][qx];
70 const real_t Qy = r1[1][qy][qx];
71 const real_t D0 = D(qx, qy, 0, e);
72
73 if (const_coeff)
74 {
75 r0[0][qy][qx] = D0 * Qx;
76 r0[1][qy][qx] = D0 * Qy;
77 }
78 if (vector_coeff)
79 {
80 const real_t D1 = D(qx, qy, 1, e);
81 r0[0][qy][qx] = D0 * Qx;
82 r0[1][qy][qx] = D1 * Qy;
83 }
84 if (matrix_coeff)
85 {
86 const real_t D1 = D(qx, qy, 1, e);
87 const real_t D2 = D(qx, qy, 2, e);
88 const real_t D3 = D(qx, qy, 3, e);
89 r0[0][qy][qx] = D0 * Qx + D1 * Qy;
90 r0[1][qy][qx] = D2 * Qx + D3 * Qy;
91 }
92 }
93 }
94 kernels::internal::EvalTranspose2d(D1D, Q1D, smem, sB, r0, r1);
95 kernels::internal::WriteDofs2d(e, D1D, r1, Y);
96 });
97}
98
99template <int T_D1D = 0, int T_Q1D = 0>
100void SmemPAVectorMassApply3D(const int NE,
101 const int coeff_vdim,
102 const Array<real_t> &b,
103 const Vector &d,
104 const Vector &x,
105 Vector &y,
106 const int d1d = 0,
107 const int q1d = 0)
108{
109 static constexpr int VDIM = 3;
110 const int D1D = T_D1D ? T_D1D : d1d;
111 const int Q1D = T_Q1D ? T_Q1D : q1d;
112
113 const bool const_coeff = coeff_vdim == 1;
114 const bool vector_coeff = coeff_vdim == VDIM;
115 const bool matrix_coeff = coeff_vdim == VDIM*VDIM;
116
117 const auto B = b.Read();
118 const auto D = Reshape(d.Read(), Q1D, Q1D, Q1D, coeff_vdim, NE);
119 const auto X = Reshape(x.Read(), D1D, D1D, D1D, VDIM, NE);
120 auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
121
122 mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
123 {
124 constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
125 constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
126
127 MFEM_SHARED real_t sB[MD1][MQ1], smem[MQ1][MQ1];
128 kernels::internal::v_regs3d_t<VDIM, MQ1> r0, r1;
129 kernels::internal::LoadMatrix(D1D, Q1D, B, sB);
130 kernels::internal::LoadDofs3d(e, D1D, X, r0);
131 kernels::internal::Eval3d(D1D, Q1D, smem, sB, r0, r1);
132
133 for (int qz = 0; qz < Q1D; qz++)
134 {
135 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
136 {
137 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
138 {
139 const real_t Qx = r1[0][qz][qy][qx];
140 const real_t Qy = r1[1][qz][qy][qx];
141 const real_t Qz = r1[2][qz][qy][qx];
142 const real_t D0 = D(qx, qy, qz, 0, e);
143 if (const_coeff)
144 {
145 r0[0][qz][qy][qx] = D0 * Qx;
146 r0[1][qz][qy][qx] = D0 * Qy;
147 r0[2][qz][qy][qx] = D0 * Qz;
148 }
149 if (vector_coeff)
150 {
151 const real_t D1 = D(qx, qy, qz, 1, e);
152 const real_t D2 = D(qx, qy, qz, 2, e);
153 r0[0][qz][qy][qx] = D0 * Qx;
154 r0[1][qz][qy][qx] = D1 * Qy;
155 r0[2][qz][qy][qx] = D2 * Qz;
156 }
157 if (matrix_coeff)
158 {
159 const real_t D1 = D(qx, qy, qz, 1, e);
160 const real_t D2 = D(qx, qy, qz, 2, e);
161 const real_t D3 = D(qx, qy, qz, 3, e);
162 const real_t D4 = D(qx, qy, qz, 4, e);
163 const real_t D5 = D(qx, qy, qz, 5, e);
164 const real_t D6 = D(qx, qy, qz, 6, e);
165 const real_t D7 = D(qx, qy, qz, 7, e);
166 const real_t D8 = D(qx, qy, qz, 8, e);
167 r0[0][qz][qy][qx] = D0 * Qx + D1 * Qy + D2 * Qz;
168 r0[1][qz][qy][qx] = D3 * Qx + D4 * Qy + D5 * Qz;
169 r0[2][qz][qy][qx] = D6 * Qx + D7 * Qy + D8 * Qz;
170 }
171 }
172 }
173 }
174 kernels::internal::EvalTranspose3d(D1D, Q1D, smem, sB, r0, r1);
175 kernels::internal::WriteDofs3d(e, D1D, r1, Y);
176 });
177}
178
179template <int T_Q1D = 0, int T_MDQ = 16>
180void SmemPAVectorMassAssembleDiagonal2D(const int ne, const int d1d,
181 const int q1d, const real_t *b_r,
182 const real_t *d_r, real_t *y_rw)
183{
184 constexpr int VDIM = 2;
185
186 const int D1D = d1d;
187 const int Q1D = T_Q1D ? T_Q1D : q1d;
188
189 MFEM_VERIFY(Q1D <= T_MDQ && D1D <= Q1D, "");
190
191 const auto B = Reshape(b_r, Q1D, D1D);
192 const auto D = Reshape(d_r, Q1D, Q1D, ne);
193 auto Y = Reshape(y_rw, D1D, D1D, VDIM, ne);
194
195 mfem::forall_2D<T_Q1D*T_Q1D>(
196 ne, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
197 {
198 constexpr int MQ1 = T_Q1D ? T_Q1D : T_MDQ;
199
200 MFEM_SHARED real_t sm[MQ1][MQ1];
201
202 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
203 {
204 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
205 {
206 real_t u = 0.0;
207 for (int qy = 0; qy < Q1D; ++qy)
208 {
209 u += B(qy, dy) * B(qy, dy) * D(qx, qy, e);
210 }
211 sm[qx][dy] = u;
212 }
213 }
214 MFEM_SYNC_THREAD;
215
216 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
217 {
218 MFEM_FOREACH_THREAD_DIRECT(dx, x, D1D)
219 {
220 real_t u = 0.0;
221 for (int qx = 0; qx < Q1D; ++qx)
222 {
223 u += B(qx, dx) * B(qx, dx) * sm[qx][dy];
224 }
225 Y(dx, dy, 0, e) += u;
226 Y(dx, dy, 1, e) += u;
227 }
228 }
229 });
230}
231
232// T_MDQ <= 10 so the Q1D^3 thread block stays within the 1024/block GPU limit
233template <int T_Q1D = 0, int T_MDQ = 10>
234void SmemPAVectorMassAssembleDiagonal3D(const int ne, const int d1d,
235 const int q1d, const real_t *b_r,
236 const real_t *d_r, real_t *y_rw)
237{
238 constexpr int VDIM = 3;
239
240 const int D1D = d1d;
241 const int Q1D = T_Q1D ? T_Q1D : q1d;
242
243 MFEM_VERIFY(Q1D <= T_MDQ && D1D <= Q1D, "");
244
245 const auto B = Reshape(b_r, Q1D, D1D);
246 const auto D = Reshape(d_r, Q1D, Q1D, Q1D, ne);
247 auto Y = Reshape(y_rw, D1D, D1D, D1D, VDIM, ne);
248
249 mfem::forall_3D<T_Q1D*T_Q1D*T_Q1D>(
250 ne, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
251 {
252 constexpr int MQ1 = T_Q1D ? T_Q1D : T_MDQ;
253
254 MFEM_SHARED real_t sm[2][MQ1][MQ1][MQ1];
255
256 MFEM_FOREACH_THREAD_DIRECT(dz, z, D1D)
257 {
258 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
259 {
260 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
261 {
262 real_t u = 0.0;
263 for (int qz = 0; qz < Q1D; ++qz)
264 {
265 u += B(qz, dz) * B(qz, dz) * D(qx, qy, qz, e);
266 }
267 sm[0][dz][qy][qx] = u;
268 }
269 }
270 }
271 MFEM_SYNC_THREAD;
272
273 MFEM_FOREACH_THREAD_DIRECT(dz, z, D1D)
274 {
275 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
276 {
277 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
278 {
279 real_t u = 0.0;
280 for (int qy = 0; qy < Q1D; ++qy)
281 {
282 u += B(qy, dy) * B(qy, dy) * sm[0][dz][qy][qx];
283 }
284 sm[1][dz][dy][qx] = u;
285 }
286 }
287 }
288 MFEM_SYNC_THREAD;
289
290 MFEM_FOREACH_THREAD_DIRECT(dz, z, D1D)
291 {
292 MFEM_FOREACH_THREAD_DIRECT(dy, y, D1D)
293 {
294 MFEM_FOREACH_THREAD_DIRECT(dx, x, D1D)
295 {
296 real_t u = 0.0;
297 for (int qx = 0; qx < Q1D; ++qx)
298 {
299 u += B(qx, dx) * B(qx, dx) * sm[1][dz][dy][qx];
300 }
301 Y(dx, dy, dz, 0, e) += u;
302 Y(dx, dy, dz, 1, e) += u;
303 Y(dx, dy, dz, 2, e) += u;
304 }
305 }
306 }
307 });
308}
309
310} // namespace internal
311
312// AddMultPA kernels
313template<int DIM, int T_D1D, int T_Q1D>
314VectorMassIntegrator::VectorMassAddMultPAType
315VectorMassIntegrator::VectorMassAddMultPA::Kernel()
316{
317 if constexpr (DIM == 2)
318 {
319 return internal::SmemPAVectorMassApply2D<T_D1D,T_Q1D>;
320 }
321 else if constexpr (DIM == 3)
322 {
323 return internal::SmemPAVectorMassApply3D<T_D1D, T_Q1D>;
324 }
325 MFEM_ABORT("Unsupported kernel");
326}
327
328inline VectorMassIntegrator::VectorMassAddMultPAType
329VectorMassIntegrator::VectorMassAddMultPA::Fallback(int dim, int, int)
330{
331 if (dim == 2)
332 {
333 return internal::SmemPAVectorMassApply2D;
334 }
335 else if (dim == 3)
336 {
337 return internal::SmemPAVectorMassApply3D;
338 }
339 MFEM_ABORT("Unsupported kernel");
340}
341
342// DiagonalPA kernels
343template<int DIM, int T_Q1D>
344VectorMassIntegrator::VectorMassAssembleDiagonalPAType
345VectorMassIntegrator::VectorMassAssembleDiagonalPA::Kernel()
346{
347 if constexpr (DIM == 2)
348 {
349 return internal::SmemPAVectorMassAssembleDiagonal2D<T_Q1D>;
350 }
351 else if constexpr (DIM == 3)
352 {
353 return internal::SmemPAVectorMassAssembleDiagonal3D<T_Q1D>;
354 }
355 MFEM_ABORT("Unsupported kernel");
356}
357
358inline VectorMassIntegrator::VectorMassAssembleDiagonalPAType
359VectorMassIntegrator::VectorMassAssembleDiagonalPA::Fallback(int dim, int)
360{
361 if (dim == 2)
362 {
363 return internal::SmemPAVectorMassAssembleDiagonal2D;
364 }
365 else if (dim == 3)
366 {
367 return internal::SmemPAVectorMassAssembleDiagonal3D;
368 }
369 MFEM_ABORT("Unsupported kernel");
370}
371
372/// \endcond DO_NOT_DOCUMENT
373
374} // namespace mfem
mfem::VectorMassIntegrator::VectorMassAssembleDiagonalPAType
void(*)(const int, const int, const int, const real_t *, const real_t *, real_t *) VectorMassAssembleDiagonalPAType
Definition bilininteg.hpp:2767
mfem::VectorMassIntegrator::VectorMassAddMultPAType
void(*)(const int, const int, const Array< real_t > &, const Vector &, const Vector &, Vector &, const int, const int) VectorMassAddMultPAType
Definition bilininteg.hpp:2758
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::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::DofQuadLimits
internal::DofQuadLimits_CUDA DofQuadLimits
Maximum number of 1D DOFs or quadrature points for the architecture currently being compiled for (use...
Definition forall.hpp:108
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:1220
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:1244
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0