1 // Copyright (c) 2010-2021, 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 "../general/forall.hpp"
13 #include "bilininteg.hpp"
14 #include "gridfunc.hpp"
15 #include "ceed/mass.hpp"
16
17 using namespace std;
18
19 namespace mfem
20 {
21
22 // PA Mass Integrator
23
24 // PA Mass Assemble kernel
AssemblePA(const FiniteElementSpace & fes)25 void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
26 {
27 // Assuming the same element type
28 Mesh *mesh = fes.GetMesh();
29 if (mesh->GetNE() == 0) { return; }
30 const FiniteElement &el = *fes.GetFE(0);
31 ElementTransformation *T = mesh->GetElementTransformation(0);
32 const IntegrationRule *ir
33 = IntRule ? IntRule : &MassIntegrator::GetRule(el, el, *T);
34 if (DeviceCanUseCeed())
35 {
36 delete ceedOp;
37 ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
38 return;
39 }
40 dim = mesh->Dimension();
41 ne = fes.GetMesh()->GetNE();
42 nq = ir->GetNPoints();
43 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::COORDINATES |
44 GeometricFactors::JACOBIANS);
45 maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
46 dofs1D = maps->ndof;
47 quad1D = maps->nqpt;
48 pa_data.SetSize(ne*nq, Device::GetDeviceMemoryType());
49 double coeff = 1.0;
50 if (Q)
51 {
52 ConstantCoefficient *cQ = dynamic_cast<ConstantCoefficient*>(Q);
53 MFEM_VERIFY(cQ != NULL, "Only ConstantCoefficient is supported.");
54 coeff = cQ->constant;
55 }
56 if (!(dim == 2 || dim == 3))
57 {
58 MFEM_ABORT("Dimension not supported.");
59 }
60 if (dim == 2)
61 {
62 const double constant = coeff;
63 const int NE = ne;
64 const int NQ = nq;
65 auto w = ir->GetWeights().Read();
66 auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
67 auto v = Reshape(pa_data.Write(), NQ, NE);
68 MFEM_FORALL(e, NE,
69 {
70 for (int q = 0; q < NQ; ++q)
71 {
72 const double J11 = J(q,0,0,e);
73 const double J12 = J(q,1,0,e);
74 const double J21 = J(q,0,1,e);
75 const double J22 = J(q,1,1,e);
76 const double detJ = (J11*J22)-(J21*J12);
77 v(q,e) = w[q] * constant * detJ;
78 }
79 });
80 }
81 if (dim == 3)
82 {
83 const double constant = coeff;
84 const int NE = ne;
85 const int NQ = nq;
86 auto W = ir->GetWeights().Read();
87 auto J = Reshape(geom->J.Read(), NQ,3,3,NE);
88 auto v = Reshape(pa_data.Write(), NQ,NE);
89 MFEM_FORALL(e, NE,
90 {
91 for (int q = 0; q < NQ; ++q)
92 {
93 const double J11 = J(q,0,0,e), J12 = J(q,0,1,e), J13 = J(q,0,2,e);
94 const double J21 = J(q,1,0,e), J22 = J(q,1,1,e), J23 = J(q,1,2,e);
95 const double J31 = J(q,2,0,e), J32 = J(q,2,1,e), J33 = J(q,2,2,e);
96 const double detJ = J11 * (J22 * J33 - J32 * J23) -
97 /* */ J21 * (J12 * J33 - J32 * J13) +
98 /* */ J31 * (J12 * J23 - J22 * J13);
99 v(q,e) = W[q] * constant * detJ;
100 }
101 });
102 }
103 }
104
105 template<const int T_D1D = 0,
106 const int T_Q1D = 0>
PAVectorMassApply2D(const int NE,const Array<double> & B_,const Array<double> & Bt_,const Vector & op_,const Vector & x_,Vector & y_,const int d1d=0,const int q1d=0)107 static void PAVectorMassApply2D(const int NE,
108 const Array<double> &B_,
109 const Array<double> &Bt_,
110 const Vector &op_,
111 const Vector &x_,
112 Vector &y_,
113 const int d1d = 0,
114 const int q1d = 0)
115 {
116 const int D1D = T_D1D ? T_D1D : d1d;
117 const int Q1D = T_Q1D ? T_Q1D : q1d;
118 constexpr int VDIM = 2;
119 MFEM_VERIFY(D1D <= MAX_D1D, "");
120 MFEM_VERIFY(Q1D <= MAX_Q1D, "");
121 auto B = Reshape(B_.Read(), Q1D, D1D);
122 auto Bt = Reshape(Bt_.Read(), D1D, Q1D);
123 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
124 auto x = Reshape(x_.Read(), D1D, D1D, VDIM, NE);
125 auto y = Reshape(y_.ReadWrite(), D1D, D1D, VDIM, NE);
126 MFEM_FORALL(e, NE,
127 {
128 const int D1D = T_D1D ? T_D1D : d1d; // nvcc workaround
129 const int Q1D = T_Q1D ? T_Q1D : q1d;
130 // the following variables are evaluated at compile time
131 constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
132 constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
133 double sol_xy[max_Q1D][max_Q1D];
134 for (int c = 0; c < VDIM; ++c)
135 {
136 for (int qy = 0; qy < Q1D; ++qy)
137 {
138 for (int qx = 0; qx < Q1D; ++qx)
139 {
140 sol_xy[qy][qx] = 0.0;
141 }
142 }
143 for (int dy = 0; dy < D1D; ++dy)
144 {
145 double sol_x[max_Q1D];
146 for (int qy = 0; qy < Q1D; ++qy)
147 {
148 sol_x[qy] = 0.0;
149 }
150 for (int dx = 0; dx < D1D; ++dx)
151 {
152 const double s = x(dx,dy,c,e);
153 for (int qx = 0; qx < Q1D; ++qx)
154 {
155 sol_x[qx] += B(qx,dx)* s;
156 }
157 }
158 for (int qy = 0; qy < Q1D; ++qy)
159 {
160 const double d2q = B(qy,dy);
161 for (int qx = 0; qx < Q1D; ++qx)
162 {
163 sol_xy[qy][qx] += d2q * sol_x[qx];
164 }
165 }
166 }
167 for (int qy = 0; qy < Q1D; ++qy)
168 {
169 for (int qx = 0; qx < Q1D; ++qx)
170 {
171 sol_xy[qy][qx] *= op(qx,qy,e);
172 }
173 }
174 for (int qy = 0; qy < Q1D; ++qy)
175 {
176 double sol_x[max_D1D];
177 for (int dx = 0; dx < D1D; ++dx)
178 {
179 sol_x[dx] = 0.0;
180 }
181 for (int qx = 0; qx < Q1D; ++qx)
182 {
183 const double s = sol_xy[qy][qx];
184 for (int dx = 0; dx < D1D; ++dx)
185 {
186 sol_x[dx] += Bt(dx,qx) * s;
187 }
188 }
189 for (int dy = 0; dy < D1D; ++dy)
190 {
191 const double q2d = Bt(dy,qy);
192 for (int dx = 0; dx < D1D; ++dx)
193 {
194 y(dx,dy,c,e) += q2d * sol_x[dx];
195 }
196 }
197 }
198 }
199 });
200 }
201
202 template<const int T_D1D = 0,
203 const int T_Q1D = 0>
PAVectorMassApply3D(const int NE,const Array<double> & B_,const Array<double> & Bt_,const Vector & op_,const Vector & x_,Vector & y_,const int d1d=0,const int q1d=0)204 static void PAVectorMassApply3D(const int NE,
205 const Array<double> &B_,
206 const Array<double> &Bt_,
207 const Vector &op_,
208 const Vector &x_,
209 Vector &y_,
210 const int d1d = 0,
211 const int q1d = 0)
212 {
213 const int D1D = T_D1D ? T_D1D : d1d;
214 const int Q1D = T_Q1D ? T_Q1D : q1d;
215 constexpr int VDIM = 3;
216 MFEM_VERIFY(D1D <= MAX_D1D, "");
217 MFEM_VERIFY(Q1D <= MAX_Q1D, "");
218 auto B = Reshape(B_.Read(), Q1D, D1D);
219 auto Bt = Reshape(Bt_.Read(), D1D, Q1D);
220 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
221 auto x = Reshape(x_.Read(), D1D, D1D, D1D, VDIM, NE);
222 auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
223 MFEM_FORALL(e, NE,
224 {
225 const int D1D = T_D1D ? T_D1D : d1d;
226 const int Q1D = T_Q1D ? T_Q1D : q1d;
227 constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
228 constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
229 double sol_xyz[max_Q1D][max_Q1D][max_Q1D];
230 for (int c = 0; c < VDIM; ++ c)
231 {
232 for (int qz = 0; qz < Q1D; ++qz)
233 {
234 for (int qy = 0; qy < Q1D; ++qy)
235 {
236 for (int qx = 0; qx < Q1D; ++qx)
237 {
238 sol_xyz[qz][qy][qx] = 0.0;
239 }
240 }
241 }
242 for (int dz = 0; dz < D1D; ++dz)
243 {
244 double sol_xy[max_Q1D][max_Q1D];
245 for (int qy = 0; qy < Q1D; ++qy)
246 {
247 for (int qx = 0; qx < Q1D; ++qx)
248 {
249 sol_xy[qy][qx] = 0.0;
250 }
251 }
252 for (int dy = 0; dy < D1D; ++dy)
253 {
254 double sol_x[max_Q1D];
255 for (int qx = 0; qx < Q1D; ++qx)
256 {
257 sol_x[qx] = 0;
258 }
259 for (int dx = 0; dx < D1D; ++dx)
260 {
261 const double s = x(dx,dy,dz,c,e);
262 for (int qx = 0; qx < Q1D; ++qx)
263 {
264 sol_x[qx] += B(qx,dx) * s;
265 }
266 }
267 for (int qy = 0; qy < Q1D; ++qy)
268 {
269 const double wy = B(qy,dy);
270 for (int qx = 0; qx < Q1D; ++qx)
271 {
272 sol_xy[qy][qx] += wy * sol_x[qx];
273 }
274 }
275 }
276 for (int qz = 0; qz < Q1D; ++qz)
277 {
278 const double wz = B(qz,dz);
279 for (int qy = 0; qy < Q1D; ++qy)
280 {
281 for (int qx = 0; qx < Q1D; ++qx)
282 {
283 sol_xyz[qz][qy][qx] += wz * sol_xy[qy][qx];
284 }
285 }
286 }
287 }
288 for (int qz = 0; qz < Q1D; ++qz)
289 {
290 for (int qy = 0; qy < Q1D; ++qy)
291 {
292 for (int qx = 0; qx < Q1D; ++qx)
293 {
294 sol_xyz[qz][qy][qx] *= op(qx,qy,qz,e);
295 }
296 }
297 }
298 for (int qz = 0; qz < Q1D; ++qz)
299 {
300 double sol_xy[max_D1D][max_D1D];
301 for (int dy = 0; dy < D1D; ++dy)
302 {
303 for (int dx = 0; dx < D1D; ++dx)
304 {
305 sol_xy[dy][dx] = 0;
306 }
307 }
308 for (int qy = 0; qy < Q1D; ++qy)
309 {
310 double sol_x[max_D1D];
311 for (int dx = 0; dx < D1D; ++dx)
312 {
313 sol_x[dx] = 0;
314 }
315 for (int qx = 0; qx < Q1D; ++qx)
316 {
317 const double s = sol_xyz[qz][qy][qx];
318 for (int dx = 0; dx < D1D; ++dx)
319 {
320 sol_x[dx] += Bt(dx,qx) * s;
321 }
322 }
323 for (int dy = 0; dy < D1D; ++dy)
324 {
325 const double wy = Bt(dy,qy);
326 for (int dx = 0; dx < D1D; ++dx)
327 {
328 sol_xy[dy][dx] += wy * sol_x[dx];
329 }
330 }
331 }
332 for (int dz = 0; dz < D1D; ++dz)
333 {
334 const double wz = Bt(dz,qz);
335 for (int dy = 0; dy < D1D; ++dy)
336 {
337 for (int dx = 0; dx < D1D; ++dx)
338 {
339 y(dx,dy,dz,c,e) += wz * sol_xy[dy][dx];
340 }
341 }
342 }
343 }
344 }
345 });
346 }
347
PAVectorMassApply(const int dim,const int D1D,const int Q1D,const int NE,const Array<double> & B,const Array<double> & Bt,const Vector & op,const Vector & x,Vector & y)348 static void PAVectorMassApply(const int dim,
349 const int D1D,
350 const int Q1D,
351 const int NE,
352 const Array<double> &B,
353 const Array<double> &Bt,
354 const Vector &op,
355 const Vector &x,
356 Vector &y)
357 {
358 if (dim == 2)
359 {
360 return PAVectorMassApply2D(NE, B, Bt, op, x, y, D1D, Q1D);
361 }
362 if (dim == 3)
363 {
364 return PAVectorMassApply3D(NE, B, Bt, op, x, y, D1D, Q1D);
365 }
366 MFEM_ABORT("Unknown kernel.");
367 }
368
AddMultPA(const Vector & x,Vector & y) const369 void VectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
370 {
371 if (DeviceCanUseCeed())
372 {
373 ceedOp->AddMult(x, y);
374 }
375 else
376 {
377 PAVectorMassApply(dim, dofs1D, quad1D, ne, maps->B, maps->Bt, pa_data, x, y);
378 }
379 }
380
381 template<const int T_D1D = 0, const int T_Q1D = 0>
PAVectorMassAssembleDiagonal2D(const int NE,const Array<double> & B_,const Array<double> & Bt_,const Vector & op_,Vector & diag_,const int d1d=0,const int q1d=0)382 static void PAVectorMassAssembleDiagonal2D(const int NE,
383 const Array<double> &B_,
384 const Array<double> &Bt_,
385 const Vector &op_,
386 Vector &diag_,
387 const int d1d = 0,
388 const int q1d = 0)
389 {
390 const int D1D = T_D1D ? T_D1D : d1d;
391 const int Q1D = T_Q1D ? T_Q1D : q1d;
392 constexpr int VDIM = 2;
393 MFEM_VERIFY(D1D <= MAX_D1D, "");
394 MFEM_VERIFY(Q1D <= MAX_Q1D, "");
395 auto B = Reshape(B_.Read(), Q1D, D1D);
396 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
397 auto y = Reshape(diag_.ReadWrite(), D1D, D1D, VDIM, NE);
398 MFEM_FORALL(e, NE,
399 {
400 const int D1D = T_D1D ? T_D1D : d1d;
401 const int Q1D = T_Q1D ? T_Q1D : q1d;
402 constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
403 constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
404
405 double temp[max_Q1D][max_D1D];
406 for (int qx = 0; qx < Q1D; ++qx)
407 {
408 for (int dy = 0; dy < D1D; ++dy)
409 {
410 temp[qx][dy] = 0.0;
411 for (int qy = 0; qy < Q1D; ++qy)
412 {
413 temp[qx][dy] += B(qy, dy) * B(qy, dy) * op(qx, qy, e);
414 }
415 }
416 }
417 for (int dy = 0; dy < D1D; ++dy)
418 {
419 for (int dx = 0; dx < D1D; ++dx)
420 {
421 double temp1 = 0.0;
422 for (int qx = 0; qx < Q1D; ++qx)
423 {
424 temp1 += B(qx, dx) * B(qx, dx) * temp[qx][dy];
425 }
426 y(dx, dy, 0, e) = temp1;
427 y(dx, dy, 1, e) = temp1;
428 }
429 }
430 });
431 }
432
433 template<const int T_D1D = 0, const int T_Q1D = 0>
PAVectorMassAssembleDiagonal3D(const int NE,const Array<double> & B_,const Array<double> & Bt_,const Vector & op_,Vector & diag_,const int d1d=0,const int q1d=0)434 static void PAVectorMassAssembleDiagonal3D(const int NE,
435 const Array<double> &B_,
436 const Array<double> &Bt_,
437 const Vector &op_,
438 Vector &diag_,
439 const int d1d = 0,
440 const int q1d = 0)
441 {
442 const int D1D = T_D1D ? T_D1D : d1d;
443 const int Q1D = T_Q1D ? T_Q1D : q1d;
444 constexpr int VDIM = 3;
445 MFEM_VERIFY(D1D <= MAX_D1D, "");
446 MFEM_VERIFY(Q1D <= MAX_Q1D, "");
447 auto B = Reshape(B_.Read(), Q1D, D1D);
448 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
449 auto y = Reshape(diag_.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
450 MFEM_FORALL(e, NE,
451 {
452 const int D1D = T_D1D ? T_D1D : d1d; // nvcc workaround
453 const int Q1D = T_Q1D ? T_Q1D : q1d;
454 // the following variables are evaluated at compile time
455 constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
456 constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
457
458 double temp[max_Q1D][max_Q1D][max_D1D];
459 for (int qx = 0; qx < Q1D; ++qx)
460 {
461 for (int qy = 0; qy < Q1D; ++qy)
462 {
463 for (int dz = 0; dz < D1D; ++dz)
464 {
465 temp[qx][qy][dz] = 0.0;
466 for (int qz = 0; qz < Q1D; ++qz)
467 {
468 temp[qx][qy][dz] += B(qz, dz) * B(qz, dz) * op(qx, qy, qz, e);
469 }
470 }
471 }
472 }
473 double temp2[max_Q1D][max_D1D][max_D1D];
474 for (int qx = 0; qx < Q1D; ++qx)
475 {
476 for (int dz = 0; dz < D1D; ++dz)
477 {
478 for (int dy = 0; dy < D1D; ++dy)
479 {
480 temp2[qx][dy][dz] = 0.0;
481 for (int qy = 0; qy < Q1D; ++qy)
482 {
483 temp2[qx][dy][dz] += B(qy, dy) * B(qy, dy) * temp[qx][qy][dz];
484 }
485 }
486 }
487 }
488 for (int dz = 0; dz < D1D; ++dz)
489 {
490 for (int dy = 0; dy < D1D; ++dy)
491 {
492 for (int dx = 0; dx < D1D; ++dx)
493 {
494 double temp3 = 0.0;
495 for (int qx = 0; qx < Q1D; ++qx)
496 {
497 temp3 += B(qx, dx) * B(qx, dx)
498 * temp2[qx][dy][dz];
499 }
500 y(dx, dy, dz, 0, e) = temp3;
501 y(dx, dy, dz, 1, e) = temp3;
502 y(dx, dy, dz, 2, e) = temp3;
503 }
504 }
505 }
506 });
507 }
508
PAVectorMassAssembleDiagonal(const int dim,const int D1D,const int Q1D,const int NE,const Array<double> & B,const Array<double> & Bt,const Vector & op,Vector & y)509 static void PAVectorMassAssembleDiagonal(const int dim,
510 const int D1D,
511 const int Q1D,
512 const int NE,
513 const Array<double> &B,
514 const Array<double> &Bt,
515 const Vector &op,
516 Vector &y)
517 {
518 if (dim == 2)
519 {
520 return PAVectorMassAssembleDiagonal2D(NE, B, Bt, op, y, D1D, Q1D);
521 }
522 else if (dim == 3)
523 {
524 return PAVectorMassAssembleDiagonal3D(NE, B, Bt, op, y, D1D, Q1D);
525 }
526 MFEM_ABORT("Dimension not implemented.");
527 }
528
AssembleDiagonalPA(Vector & diag)529 void VectorMassIntegrator::AssembleDiagonalPA(Vector &diag)
530 {
531 if (DeviceCanUseCeed())
532 {
533 ceedOp->GetDiagonal(diag);
534 }
535 else
536 {
537 PAVectorMassAssembleDiagonal(dim,
538 dofs1D,
539 quad1D,
540 ne,
541 maps->B,
542 maps->Bt,
543 pa_data,
544 diag);
545 }
546 }
547
548 } // namespace mfem
549