MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_diffusion_kernels.cpp
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#include "bilininteg_diffusion_kernels.hpp"
13#include "bilininteg_diffusion_pa_simplices.hpp" // IWYU pragma: keep
14
15namespace mfem
16{
17
18// PA Diffusion Integrator
19
20DiffusionIntegrator::Kernels::Kernels()
21{
22 // 2D
23 // Q = P, only for simplex
24 DiffusionIntegrator::AddSimplexSpecialization<2,2,1>();
25 DiffusionIntegrator::AddSimplexSpecialization<2,3,2>();
26 DiffusionIntegrator::AddSimplexSpecialization<2,4,3>();
27 DiffusionIntegrator::AddSimplexSpecialization<2,5,4>();
28 DiffusionIntegrator::AddSimplexSpecialization<2,6,5>();
29 DiffusionIntegrator::AddSimplexSpecialization<2,7,6>();
30 // Q = P+1
31 DiffusionIntegrator::AddSpecialization<2,1,1>();
32 DiffusionIntegrator::AddSpecialization<2,2,2>();
33 DiffusionIntegrator::AddSpecialization<2,3,3>();
34 DiffusionIntegrator::AddSpecialization<2,4,4>();
35 DiffusionIntegrator::AddSpecialization<2,5,5>();
36 DiffusionIntegrator::AddSpecialization<2,6,6>();
37 DiffusionIntegrator::AddSpecialization<2,7,7>();
38 DiffusionIntegrator::AddSpecialization<2,8,8>();
39 DiffusionIntegrator::AddSpecialization<2,9,9>();
40 // Q = P+2
41 DiffusionIntegrator::AddSpecialization<2,1,2>();
42 DiffusionIntegrator::AddSpecialization<2,2,3>();
43 DiffusionIntegrator::AddSpecialization<2,3,4>();
44 DiffusionIntegrator::AddSpecialization<2,4,5>();
45 DiffusionIntegrator::AddSpecialization<2,5,6>();
46 DiffusionIntegrator::AddSpecialization<2,6,7>();
47 DiffusionIntegrator::AddSpecialization<2,7,8>();
48 DiffusionIntegrator::AddSpecialization<2,8,9>();
49 DiffusionIntegrator::AddSpecialization<2,9,10>();
50 // others
51 DiffusionIntegrator::AddSimplexSpecialization<2,2,5>();
52 DiffusionIntegrator::AddSimplexSpecialization<2,3,6>();
53
54 // 3D
55 // Q = P, only for simplex
56 DiffusionIntegrator::AddSimplexSpecialization<3,2,1>();
57 DiffusionIntegrator::AddSimplexSpecialization<3,3,2>();
58 DiffusionIntegrator::AddSimplexSpecialization<3,4,3>();
59 DiffusionIntegrator::AddSimplexSpecialization<3,5,4>();
60 DiffusionIntegrator::AddSimplexSpecialization<3,6,5>();
61 DiffusionIntegrator::AddSimplexSpecialization<3,7,6>();
62 DiffusionIntegrator::AddSimplexSpecialization<3,8,7>();
63 // Q = P+1
64 DiffusionIntegrator::AddSpecialization<3,1,1>();
65 DiffusionIntegrator::AddSpecialization<3,2,2>();
66 DiffusionIntegrator::AddSpecialization<3,3,3>();
67 DiffusionIntegrator::AddSpecialization<3,4,4>();
68 DiffusionIntegrator::AddSpecialization<3,5,5>();
69 DiffusionIntegrator::AddSpecialization<3,6,6>();
70 DiffusionIntegrator::AddSpecialization<3,7,7>();
71 DiffusionIntegrator::AddSpecialization<3,8,8>();
72 // Q = P+2
73 DiffusionIntegrator::AddSpecialization<3,1,2>();
74 DiffusionIntegrator::AddSpecialization<3,2,3>();
75 DiffusionIntegrator::AddSpecialization<3,3,4>();
76 DiffusionIntegrator::AddSpecialization<3,4,5>();
77 DiffusionIntegrator::AddSpecialization<3,5,6>();
78 DiffusionIntegrator::AddSpecialization<3,6,7>();
79 DiffusionIntegrator::AddSpecialization<3,7,8>();
80 DiffusionIntegrator::AddSpecialization<3,8,9>();
81 // others
82 DiffusionIntegrator::AddSpecialization<3,4,6>();
83 DiffusionIntegrator::AddSpecialization<3,5,8>();
84}
85
86namespace internal
87{
88
89template<>
90void PADiffusionSetup2D<2>(const int Q1D,
91 const int coeffDim,
92 const int NE,
93 const Array<real_t> &w,
94 const Vector &j,
95 const Vector &c,
96 Vector &d);
97
98template<>
99void PADiffusionSetup2D<3>(const int Q1D,
100 const int coeffDim,
101 const int NE,
102 const Array<real_t> &w,
103 const Vector &j,
104 const Vector &c,
105 Vector &d);
106
107void PADiffusionSetup(const int dim,
108 const int sdim,
109 const int D1D,
110 const int Q1D,
111 const int coeffDim,
112 const int NE,
113 const Array<real_t> &W,
114 const Vector &J,
115 const Vector &C,
116 Vector &D)
117{
118 if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADiffusionSetup"); }
119 if (dim == 2)
120 {
121#ifdef MFEM_USE_OCCA
122 if (DeviceCanUseOcca())
123 {
124 OccaPADiffusionSetup2D(D1D, Q1D, NE, W, J, C, D);
125 return;
126 }
127#else
128 MFEM_CONTRACT_VAR(D1D);
129#endif // MFEM_USE_OCCA
130 if (sdim == 2) { PADiffusionSetup2D<2>(Q1D, coeffDim, NE, W, J, C, D); }
131 if (sdim == 3) { PADiffusionSetup2D<3>(Q1D, coeffDim, NE, W, J, C, D); }
132 }
133 if (dim == 3)
134 {
135#ifdef MFEM_USE_OCCA
136 if (DeviceCanUseOcca())
137 {
138 OccaPADiffusionSetup3D(D1D, Q1D, NE, W, J, C, D);
139 return;
140 }
141#endif // MFEM_USE_OCCA
142 PADiffusionSetup3D(Q1D, coeffDim, NE, W, J, C, D);
143 }
144}
145
146template<>
147void PADiffusionSetup2D<2>(const int Q1D,
148 const int coeffDim,
149 const int NE,
150 const Array<real_t> &w,
151 const Vector &j,
152 const Vector &c,
153 Vector &d)
154{
155 const bool symmetric = (coeffDim != 4);
156 const bool const_c = c.Size() == coeffDim;
157 const auto W = Reshape(w.Read(), Q1D,Q1D);
158 const auto J = Reshape(j.Read(), Q1D,Q1D,2,2,NE);
159 const auto C = const_c ? Reshape(c.Read(), coeffDim,1,1,1) :
160 Reshape(c.Read(), coeffDim,Q1D,Q1D,NE);
161 auto D = Reshape(d.Write(), Q1D,Q1D, symmetric ? 3 : 4, NE);
162
163 auto get_coeff = [const_c] MFEM_HOST_DEVICE
164 (const decltype(C) &C, int i, int qx, int qy, int e)
165 {
166 return const_c ? C(i,0,0,0) : C(i,qx,qy,e);
167 };
168
169 mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
170 {
171 MFEM_FOREACH_THREAD(qx,x,Q1D)
172 {
173 MFEM_FOREACH_THREAD(qy,y,Q1D)
174 {
175 const real_t J11 = J(qx,qy,0,0,e);
176 const real_t J21 = J(qx,qy,1,0,e);
177 const real_t J12 = J(qx,qy,0,1,e);
178 const real_t J22 = J(qx,qy,1,1,e);
179 const real_t w_detJ = W(qx,qy) / ((J11*J22)-(J21*J12));
180 if (coeffDim == 3 || coeffDim == 4) // Matrix coefficient
181 {
182 // First compute entries of R = MJ^{-T}, without det J factor.
183 const real_t M11 = get_coeff(C,0,qx,qy,e);
184 const real_t M12 = get_coeff(C,1,qx,qy,e);
185 const real_t M21 = symmetric ? M12 : get_coeff(C,2,qx,qy,e);
186 const real_t M22 = symmetric ? get_coeff(C,2,qx,qy,e)
187 : get_coeff(C,3,qx,qy,e);
188 const real_t R11 = M11*J22 - M12*J12;
189 const real_t R21 = M21*J22 - M22*J12;
190 const real_t R12 = -M11*J21 + M12*J11;
191 const real_t R22 = -M21*J21 + M22*J11;
192
193 // Now set y to J^{-1}R.
194 D(qx,qy,0,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
195 D(qx,qy,1,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1
196 D(qx,qy,2,e) = w_detJ * (symmetric ? (-J21*R12 + J11*R22) :
197 (J22*R12 - J12*R22)); // 2,2 or 1,2
198 if (!symmetric)
199 {
200 D(qx,qy,3,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
201 }
202 }
203 else // Vector or scalar coefficient
204 {
205 const real_t C1 = get_coeff(C,0,qx,qy,e);
206 const real_t C2 = get_coeff(C,coeffDim==2?1:0,qx,qy,e);
207
208 D(qx,qy,0,e) = w_detJ * (C2*J12*J12 + C1*J22*J22); // 1,1
209 D(qx,qy,1,e) = -w_detJ * (C2*J12*J11 + C1*J22*J21); // 1,2
210 D(qx,qy,2,e) = w_detJ * (C2*J11*J11 + C1*J21*J21); // 2,2
211 }
212 }
213 }
214 });
215}
216
217template<>
218void PADiffusionSetup2D<3>(const int Q1D,
219 const int coeffDim,
220 const int NE,
221 const Array<real_t> &w,
222 const Vector &j,
223 const Vector &c,
224 Vector &d)
225{
226 MFEM_VERIFY(coeffDim == 1, "Matrix and vector coefficients not supported");
227 constexpr int DIM = 2;
228 constexpr int SDIM = 3;
229 const bool const_c = c.Size() == 1;
230 const auto W = Reshape(w.Read(), Q1D,Q1D);
231 const auto J = Reshape(j.Read(), Q1D,Q1D,SDIM,DIM,NE);
232 const auto C = const_c ? Reshape(c.Read(), 1,1,1) :
233 Reshape(c.Read(), Q1D,Q1D,NE);
234 auto D = Reshape(d.Write(), Q1D,Q1D, 3, NE);
235 mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
236 {
237 MFEM_FOREACH_THREAD(qx,x,Q1D)
238 {
239 MFEM_FOREACH_THREAD(qy,y,Q1D)
240 {
241 const real_t wq = W(qx,qy);
242 const real_t J11 = J(qx,qy,0,0,e);
243 const real_t J21 = J(qx,qy,1,0,e);
244 const real_t J31 = J(qx,qy,2,0,e);
245 const real_t J12 = J(qx,qy,0,1,e);
246 const real_t J22 = J(qx,qy,1,1,e);
247 const real_t J32 = J(qx,qy,2,1,e);
248 const real_t E = J11*J11 + J21*J21 + J31*J31;
249 const real_t G = J12*J12 + J22*J22 + J32*J32;
250 const real_t F = J11*J12 + J21*J22 + J31*J32;
251 const real_t iw = 1.0 / std::sqrt(E*G - F*F);
252 const real_t coeff = const_c ? C(0,0,0) : C(qx,qy,e);
253 const real_t alpha = wq * coeff * iw;
254 D(qx,qy,0,e) = alpha * G; // 1,1
255 D(qx,qy,1,e) = -alpha * F; // 1,2
256 D(qx,qy,2,e) = alpha * E; // 2,2
257 }
258 }
259 });
260}
261
262void PADiffusionSetup3D(const int Q1D,
263 const int coeffDim,
264 const int NE,
265 const Array<real_t> &w,
266 const Vector &j,
267 const Vector &c,
268 Vector &d)
269{
270 const bool symmetric = (coeffDim != 9);
271 const bool const_c = c.Size() == coeffDim;
272 const auto W = Reshape(w.Read(), Q1D,Q1D,Q1D);
273 const auto J = Reshape(j.Read(), Q1D,Q1D,Q1D,3,3,NE);
274 const auto C = const_c ? Reshape(c.Read(), coeffDim,1,1,1,1) :
275 Reshape(c.Read(), coeffDim,Q1D,Q1D,Q1D,NE);
276 auto D = Reshape(d.Write(), Q1D,Q1D,Q1D, symmetric ? 6 : 9, NE);
277
278 auto get_coeff = [const_c] MFEM_HOST_DEVICE
279 (const decltype(C) &C, int i, int qx, int qy, int qz, int e)
280 {
281 return const_c ? C(i,0,0,0,0) : C(i,qx,qy,qz,e);
282 };
283
284 mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
285 {
286 MFEM_FOREACH_THREAD(qx,x,Q1D)
287 {
288 MFEM_FOREACH_THREAD(qy,y,Q1D)
289 {
290 MFEM_FOREACH_THREAD(qz,z,Q1D)
291 {
292 const real_t J11 = J(qx,qy,qz,0,0,e);
293 const real_t J21 = J(qx,qy,qz,1,0,e);
294 const real_t J31 = J(qx,qy,qz,2,0,e);
295 const real_t J12 = J(qx,qy,qz,0,1,e);
296 const real_t J22 = J(qx,qy,qz,1,1,e);
297 const real_t J32 = J(qx,qy,qz,2,1,e);
298 const real_t J13 = J(qx,qy,qz,0,2,e);
299 const real_t J23 = J(qx,qy,qz,1,2,e);
300 const real_t J33 = J(qx,qy,qz,2,2,e);
301 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
302 J21 * (J12 * J33 - J32 * J13) +
303 J31 * (J12 * J23 - J22 * J13);
304 const real_t w_detJ = W(qx,qy,qz) / detJ;
305 // adj(J)
306 const real_t A11 = (J22 * J33) - (J23 * J32);
307 const real_t A12 = (J32 * J13) - (J12 * J33);
308 const real_t A13 = (J12 * J23) - (J22 * J13);
309 const real_t A21 = (J31 * J23) - (J21 * J33);
310 const real_t A22 = (J11 * J33) - (J13 * J31);
311 const real_t A23 = (J21 * J13) - (J11 * J23);
312 const real_t A31 = (J21 * J32) - (J31 * J22);
313 const real_t A32 = (J31 * J12) - (J11 * J32);
314 const real_t A33 = (J11 * J22) - (J12 * J21);
315
316 if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
317 {
318 // Compute entries of R = MJ^{-T} = M adj(J)^T, without det J.
319 const real_t M11 = get_coeff(C, 0, qx,qy,qz, e);
320 const real_t M12 = get_coeff(C, 1, qx,qy,qz, e);
321 const real_t M13 = get_coeff(C, 2, qx,qy,qz, e);
322 const real_t M21 = (!symmetric) ? get_coeff(C, 3, qx,qy,qz, e) : M12;
323 const real_t M22 = (!symmetric) ? get_coeff(C, 4, qx,qy,qz, e)
324 : get_coeff(C, 3, qx,qy,qz, e);
325 const real_t M23 = (!symmetric) ? get_coeff(C, 5, qx,qy,qz, e)
326 : get_coeff(C, 4, qx,qy,qz, e);
327 const real_t M31 = (!symmetric) ? get_coeff(C, 6, qx,qy,qz, e) : M13;
328 const real_t M32 = (!symmetric) ? get_coeff(C, 7, qx,qy,qz, e) : M23;
329 const real_t M33 = (!symmetric) ? get_coeff(C, 8, qx,qy,qz, e)
330 : get_coeff(C, 5, qx,qy,qz, e);
331
332 const real_t R11 = M11*A11 + M12*A12 + M13*A13;
333 const real_t R12 = M11*A21 + M12*A22 + M13*A23;
334 const real_t R13 = M11*A31 + M12*A32 + M13*A33;
335 const real_t R21 = M21*A11 + M22*A12 + M23*A13;
336 const real_t R22 = M21*A21 + M22*A22 + M23*A23;
337 const real_t R23 = M21*A31 + M22*A32 + M23*A33;
338 const real_t R31 = M31*A11 + M32*A12 + M33*A13;
339 const real_t R32 = M31*A21 + M32*A22 + M33*A23;
340 const real_t R33 = M31*A31 + M32*A32 + M33*A33;
341
342 // Now set D to J^{-1} R = adj(J) R
343 D(qx,qy,qz,0,e) = w_detJ * (A11*R11 + A12*R21 + A13*R31); // 1,1
344 const real_t D12 = w_detJ * (A11*R12 + A12*R22 + A13*R32);
345 D(qx,qy,qz,1,e) = D12; // 1,2
346 D(qx,qy,qz,2,e) = w_detJ * (A11*R13 + A12*R23 + A13*R33); // 1,3
347
348 const real_t D22 = w_detJ * (A21*R12 + A22*R22 + A23*R32);
349 const real_t D23 = w_detJ * (A21*R13 + A22*R23 + A23*R33);
350
351 const real_t D33 = w_detJ * (A31*R13 + A32*R23 + A33*R33);
352
353 D(qx,qy,qz,4,e) = symmetric ? D23 : D22; // 2,3 or 2,2
354 D(qx,qy,qz,5,e) = symmetric ? D33 : D23; // 3,3 or 2,3
355
356 if (symmetric)
357 {
358 D(qx,qy,qz,3,e) = D22; // 2,2
359 }
360 else
361 {
362 D(qx,qy,qz,3,e) = w_detJ * (A21*R11 + A22*R21 + A23*R31); // 2,1
363 D(qx,qy,qz,6,e) = w_detJ * (A31*R11 + A32*R21 + A33*R31); // 3,1
364 D(qx,qy,qz,7,e) = w_detJ * (A31*R12 + A32*R22 + A33*R32); // 3,2
365 D(qx,qy,qz,8,e) = D33; // 3,3
366 }
367 }
368 else // Vector or scalar coefficient version
369 {
370 const real_t C1 = get_coeff(C,0,qx,qy,qz,e);
371 const real_t C2 = get_coeff(C,coeffDim==3?1:0,qx,qy,qz,e);
372 const real_t C3 = get_coeff(C,coeffDim==3?2:0,qx,qy,qz,e);
373
374 // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
375 D(qx,qy,qz,0,e) = w_detJ * (C1*A11*A11 + C2*A12*A12 + C3*A13*A13); // 1,1
376 D(qx,qy,qz,1,e) = w_detJ * (C1*A11*A21 + C2*A12*A22 + C3*A13*A23); // 2,1
377 D(qx,qy,qz,2,e) = w_detJ * (C1*A11*A31 + C2*A12*A32 + C3*A13*A33); // 3,1
378 D(qx,qy,qz,3,e) = w_detJ * (C1*A21*A21 + C2*A22*A22 + C3*A23*A23); // 2,2
379 D(qx,qy,qz,4,e) = w_detJ * (C1*A21*A31 + C2*A22*A32 + C3*A23*A33); // 3,2
380 D(qx,qy,qz,5,e) = w_detJ * (C1*A31*A31 + C2*A32*A32 + C3*A33*A33); // 3,3
381 }
382 }
383 }
384 }
385 });
386}
387
388#ifdef MFEM_USE_OCCA
389void OccaPADiffusionSetup2D(const int D1D,
390 const int Q1D,
391 const int NE,
392 const Array<real_t> &W,
393 const Vector &J,
394 const Vector &C,
395 Vector &op)
396{
397 occa::properties props;
398 props["defines/D1D"] = D1D;
399 props["defines/Q1D"] = Q1D;
400 const occa::memory o_W = OccaMemoryRead(W.GetMemory(), W.Size());
401 const occa::memory o_J = OccaMemoryRead(J.GetMemory(), J.Size());
402 const occa::memory o_C = OccaMemoryRead(C.GetMemory(), C.Size());
403 occa::memory o_op = OccaMemoryWrite(op.GetMemory(), op.Size());
404 const bool const_c = C.Size() == 1;
405 const occa_id_t id = std::make_pair(D1D,Q1D);
406 static occa_kernel_t OccaDiffSetup2D_ker;
407 if (OccaDiffSetup2D_ker.find(id) == OccaDiffSetup2D_ker.end())
408 {
409 const occa::kernel DiffusionSetup2D =
410 mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl",
411 "DiffusionSetup2D", props);
412 OccaDiffSetup2D_ker.emplace(id, DiffusionSetup2D);
413 }
414 OccaDiffSetup2D_ker.at(id)(NE, o_W, o_J, o_C, o_op, const_c);
415}
416
417void OccaPADiffusionSetup3D(const int D1D,
418 const int Q1D,
419 const int NE,
420 const Array<real_t> &W,
421 const Vector &J,
422 const Vector &C,
423 Vector &op)
424{
425 occa::properties props;
426 props["defines/D1D"] = D1D;
427 props["defines/Q1D"] = Q1D;
428 const occa::memory o_W = OccaMemoryRead(W.GetMemory(), W.Size());
429 const occa::memory o_J = OccaMemoryRead(J.GetMemory(), J.Size());
430 const occa::memory o_C = OccaMemoryRead(C.GetMemory(), C.Size());
431 occa::memory o_op = OccaMemoryWrite(op.GetMemory(), op.Size());
432 const bool const_c = C.Size() == 1;
433 const occa_id_t id = std::make_pair(D1D,Q1D);
434 static occa_kernel_t OccaDiffSetup3D_ker;
435 if (OccaDiffSetup3D_ker.find(id) == OccaDiffSetup3D_ker.end())
436 {
437 const occa::kernel DiffusionSetup3D =
438 mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl",
439 "DiffusionSetup3D", props);
440 OccaDiffSetup3D_ker.emplace(id, DiffusionSetup3D);
441 }
442 OccaDiffSetup3D_ker.at(id)(NE, o_W, o_J, o_C, o_op, const_c);
443}
444
445void OccaPADiffusionApply2D(const int D1D,
446 const int Q1D,
447 const int NE,
448 const Array<real_t> &B,
449 const Array<real_t> &G,
450 const Array<real_t> &Bt,
451 const Array<real_t> &Gt,
452 const Vector &D,
453 const Vector &X,
454 Vector &Y)
455{
456 occa::properties props;
457 props["defines/D1D"] = D1D;
458 props["defines/Q1D"] = Q1D;
459 const occa::memory o_B = OccaMemoryRead(B.GetMemory(), B.Size());
460 const occa::memory o_G = OccaMemoryRead(G.GetMemory(), G.Size());
461 const occa::memory o_Bt = OccaMemoryRead(Bt.GetMemory(), Bt.Size());
462 const occa::memory o_Gt = OccaMemoryRead(Gt.GetMemory(), Gt.Size());
463 const occa::memory o_D = OccaMemoryRead(D.GetMemory(), D.Size());
464 const occa::memory o_X = OccaMemoryRead(X.GetMemory(), X.Size());
465 occa::memory o_Y = OccaMemoryReadWrite(Y.GetMemory(), Y.Size());
466 const occa_id_t id = std::make_pair(D1D,Q1D);
467 if (!Device::Allows(Backend::OCCA_CUDA))
468 {
469 static occa_kernel_t OccaDiffApply2D_cpu;
470 if (OccaDiffApply2D_cpu.find(id) == OccaDiffApply2D_cpu.end())
471 {
472 const occa::kernel DiffusionApply2D_CPU =
473 mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl",
474 "DiffusionApply2D_CPU", props);
475 OccaDiffApply2D_cpu.emplace(id, DiffusionApply2D_CPU);
476 }
477 OccaDiffApply2D_cpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_D, o_X, o_Y);
478 }
479 else
480 {
481 static occa_kernel_t OccaDiffApply2D_gpu;
482 if (OccaDiffApply2D_gpu.find(id) == OccaDiffApply2D_gpu.end())
483 {
484 const occa::kernel DiffusionApply2D_GPU =
485 mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl",
486 "DiffusionApply2D_GPU", props);
487 OccaDiffApply2D_gpu.emplace(id, DiffusionApply2D_GPU);
488 }
489 OccaDiffApply2D_gpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_D, o_X, o_Y);
490 }
491}
492
493void OccaPADiffusionApply3D(const int D1D,
494 const int Q1D,
495 const int NE,
496 const Array<real_t> &B,
497 const Array<real_t> &G,
498 const Array<real_t> &Bt,
499 const Array<real_t> &Gt,
500 const Vector &D,
501 const Vector &X,
502 Vector &Y)
503{
504 occa::properties props;
505 props["defines/D1D"] = D1D;
506 props["defines/Q1D"] = Q1D;
507 const occa::memory o_B = OccaMemoryRead(B.GetMemory(), B.Size());
508 const occa::memory o_G = OccaMemoryRead(G.GetMemory(), G.Size());
509 const occa::memory o_Bt = OccaMemoryRead(Bt.GetMemory(), Bt.Size());
510 const occa::memory o_Gt = OccaMemoryRead(Gt.GetMemory(), Gt.Size());
511 const occa::memory o_D = OccaMemoryRead(D.GetMemory(), D.Size());
512 const occa::memory o_X = OccaMemoryRead(X.GetMemory(), X.Size());
513 occa::memory o_Y = OccaMemoryReadWrite(Y.GetMemory(), Y.Size());
514 const occa_id_t id = std::make_pair(D1D,Q1D);
515 if (!Device::Allows(Backend::OCCA_CUDA))
516 {
517 static occa_kernel_t OccaDiffApply3D_cpu;
518 if (OccaDiffApply3D_cpu.find(id) == OccaDiffApply3D_cpu.end())
519 {
520 const occa::kernel DiffusionApply3D_CPU =
521 mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl",
522 "DiffusionApply3D_CPU", props);
523 OccaDiffApply3D_cpu.emplace(id, DiffusionApply3D_CPU);
524 }
525 OccaDiffApply3D_cpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_D, o_X, o_Y);
526 }
527 else
528 {
529 static occa_kernel_t OccaDiffApply3D_gpu;
530 if (OccaDiffApply3D_gpu.find(id) == OccaDiffApply3D_gpu.end())
531 {
532 const occa::kernel DiffusionApply3D_GPU =
533 mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl",
534 "DiffusionApply3D_GPU", props);
535 OccaDiffApply3D_gpu.emplace(id, DiffusionApply3D_GPU);
536 }
537 OccaDiffApply3D_gpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_D, o_X, o_Y);
538 }
539}
540#endif // MFEM_USE_OCCA
541
542} // namespace internal
543
544} // namespace mfem
bilininteg_diffusion_kernels.hpp
bilininteg_diffusion_pa_simplices.hpp
mfem::DiffusionIntegrator::AddSimplexSpecialization
static void AddSimplexSpecialization()
Definition bilininteg.hpp:2436
mfem::DiffusionIntegrator::AddSpecialization
static void AddSpecialization()
Definition bilininteg.hpp:2428
mfem::Vector
Vector data type.
Definition vector.hpp:82
alpha
const real_t alpha
Definition ex15.cpp:369
dim
int dim
Definition ex24.cpp:53
SDIM
constexpr int SDIM
Definition minimal-surface.cpp:72
DIM
constexpr int DIM
Definition minimal-surface.cpp:71
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
mfem::OccaMemoryReadWrite
occa::memory OccaMemoryReadWrite(Memory< T > &mem, size_t size)
Wrap a Memory object as occa::memory for read-write access with the mfem::Device MemoryClass....
Definition occa.hpp:59
mfem::Read
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition device.hpp:369
mfem::OccaMemoryWrite
occa::memory OccaMemoryWrite(Memory< T > &mem, size_t size)
Wrap a Memory object as occa::memory for write only access with the mfem::Device MemoryClass....
Definition occa.hpp:48
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::occa_kernel_t
std::map< occa_id_t, occa::kernel > occa_kernel_t
Definition occa.hpp:79
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:1244
mfem::DeviceCanUseOcca
bool DeviceCanUseOcca()
Function that determines if an OCCA kernel should be used, based on the current mfem::Device configur...
Definition occa.hpp:69
mfem::OccaMemoryRead
const occa::memory OccaMemoryRead(const Memory< T > &mem, size_t size)
Wrap a Memory object as occa::memory for read only access with the mfem::Device MemoryClass....
Definition occa.hpp:37
mfem::OccaDev
occa::device & OccaDev()
Return the default occa::device used by MFEM.
Definition occa.cpp:27
mfem::occa_id_t
std::pair< int, int > occa_id_t
Definition occa.hpp:78
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0