MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilininteg_convection_pa.cpp
Go to the documentation of this file.
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/convection.hpp"
16 #include "quadinterpolator.hpp"
17 
18 namespace mfem
19 {
20 
21 // PA Convection Integrator
22 
23 // PA Convection Assemble 2D kernel
24 static void PAConvectionSetup2D(const int NQ,
25  const int NE,
26  const Array<double> &w,
27  const Vector &j,
28  const Vector &vel,
29  const double alpha,
30  Vector &op)
31 {
32  constexpr int DIM = 2;
33 
34  const bool const_v = vel.Size() == DIM;
35 
36  const auto W = w.Read();
37  const auto J = Reshape(j.Read(), NQ,DIM,DIM,NE);
38  const auto V = const_v ?
39  Reshape(vel.Read(), DIM,1,1) :
40  Reshape(vel.Read(), DIM,NQ,NE);
41  auto y = Reshape(op.Write(), NQ,DIM,NE);
42 
43  MFEM_FORALL(q_global, NE*NQ,
44  {
45  const int e = q_global / NQ;
46  const int q = q_global % NQ;
47  const double J11 = J(q,0,0,e);
48  const double J21 = J(q,1,0,e);
49  const double J12 = J(q,0,1,e);
50  const double J22 = J(q,1,1,e);
51  const double w = alpha * W[q];
52  const double v0 = const_v ? V(0,0,0) : V(0,q,e);
53  const double v1 = const_v ? V(1,0,0) : V(1,q,e);
54  const double wx = w * v0;
55  const double wy = w * v1;
56  // y = alpha * W * det(J) * J^{-1} . v = adj(J) . { wx, wy }
57  y(q,0,e) = wx * J22 - wy * J12; // 1
58  y(q,1,e) = -wx * J21 + wy * J11; // 2
59  });
60 }
61 
62 // PA Convection Assemble 3D kernel
63 static void PAConvectionSetup3D(const int NQ,
64  const int NE,
65  const Array<double> &w,
66  const Vector &j,
67  const Vector &vel,
68  const double alpha,
69  Vector &op)
70 {
71  constexpr int DIM = 3;
72  constexpr int SDIM = DIM;
73  const auto W = Reshape(w.Read(), NQ);
74  const auto J = Reshape(j.Read(), NQ,SDIM,DIM,NE);
75  const bool const_v = vel.Size() == DIM;
76  const auto V = const_v ?
77  Reshape(vel.Read(), 3,1,1) :
78  Reshape(vel.Read(), 3,NQ,NE);
79  auto y = Reshape(op.Write(), NQ,3,NE);
80  MFEM_FORALL(q_global, NE*NQ,
81  {
82  const int e = q_global / NQ;
83  const int q = q_global % NQ;
84  const double J11 = J(q,0,0,e);
85  const double J12 = J(q,0,1,e);
86  const double J13 = J(q,0,2,e);
87  const double J21 = J(q,1,0,e);
88  const double J22 = J(q,1,1,e);
89  const double J23 = J(q,1,2,e);
90  const double J31 = J(q,2,0,e);
91  const double J32 = J(q,2,1,e);
92  const double J33 = J(q,2,2,e);
93  const double w = alpha * W(q);
94  const double v0 = const_v ? V(0,0,0) : V(0,q,e);
95  const double v1 = const_v ? V(1,0,0) : V(1,q,e);
96  const double v2 = const_v ? V(2,0,0) : V(2,q,e);
97  const double wx = w * v0;
98  const double wy = w * v1;
99  const double wz = w * v2;
100  // A = adj(J)
101  const double A11 = (J22 * J33) - (J23 * J32);
102  const double A12 = (J32 * J13) - (J12 * J33);
103  const double A13 = (J12 * J23) - (J22 * J13);
104  const double A21 = (J31 * J23) - (J21 * J33);
105  const double A22 = (J11 * J33) - (J13 * J31);
106  const double A23 = (J21 * J13) - (J11 * J23);
107  const double A31 = (J21 * J32) - (J31 * J22);
108  const double A32 = (J31 * J12) - (J11 * J32);
109  const double A33 = (J11 * J22) - (J12 * J21);
110  // y = alpha * W * det(J) * J^{-1} . v = adj(J) . { wx, wy, wz }
111  y(q,0,e) = wx * A11 + wy * A12 + wz * A13;
112  y(q,1,e) = wx * A21 + wy * A22 + wz * A23;
113  y(q,2,e) = wx * A31 + wy * A32 + wz * A33;
114  });
115 }
116 
117 static void PAConvectionSetup(const int dim,
118  const int NQ,
119  const int NE,
120  const Array<double> &W,
121  const Vector &J,
122  const Vector &coeff,
123  const double alpha,
124  Vector &op)
125 {
126  if (dim == 1) { MFEM_ABORT("dim==1 not supported in PAConvectionSetup"); }
127  if (dim == 2)
128  {
129  PAConvectionSetup2D(NQ, NE, W, J, coeff, alpha, op);
130  }
131  if (dim == 3)
132  {
133  PAConvectionSetup3D(NQ, NE, W, J, coeff, alpha, op);
134  }
135 }
136 
137 // PA Convection Apply 2D kernel
138 template<int T_D1D = 0, int T_Q1D = 0> static
139 void PAConvectionApply2D(const int ne,
140  const Array<double> &b,
141  const Array<double> &g,
142  const Array<double> &bt,
143  const Array<double> &gt,
144  const Vector &op_,
145  const Vector &x_,
146  Vector &y_,
147  const int d1d = 0,
148  const int q1d = 0)
149 {
150  const int NE = ne;
151  const int D1D = T_D1D ? T_D1D : d1d;
152  const int Q1D = T_Q1D ? T_Q1D : q1d;
153  MFEM_VERIFY(D1D <= MAX_D1D, "");
154  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
155  auto B = Reshape(b.Read(), Q1D, D1D);
156  auto G = Reshape(g.Read(), Q1D, D1D);
157  auto Bt = Reshape(bt.Read(), D1D, Q1D);
158  auto op = Reshape(op_.Read(), Q1D, Q1D, 2, NE);
159  auto x = Reshape(x_.Read(), D1D, D1D, NE);
160  auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
161  MFEM_FORALL(e, NE,
162  {
163  const int D1D = T_D1D ? T_D1D : d1d;
164  const int Q1D = T_Q1D ? T_Q1D : q1d;
165  // the following variables are evaluated at compile time
166  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
167  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
168 
169  double u[max_D1D][max_D1D];
170  for (int dy = 0; dy < D1D; ++dy)
171  {
172  for (int dx = 0; dx < D1D; ++dx)
173  {
174  u[dy][dx] = x(dx,dy,e);
175  }
176  }
177  double Bu[max_D1D][max_Q1D];
178  double Gu[max_D1D][max_Q1D];
179  for (int dy = 0; dy < D1D; ++dy)
180  {
181  for (int qx = 0; qx < Q1D; ++qx)
182  {
183  Bu[dy][qx] = 0.0;
184  Gu[dy][qx] = 0.0;
185  for (int dx = 0; dx < D1D; ++dx)
186  {
187  const double bx = B(qx,dx);
188  const double gx = G(qx,dx);
189  const double x = u[dy][dx];
190  Bu[dy][qx] += bx * x;
191  Gu[dy][qx] += gx * x;
192  }
193  }
194  }
195  double GBu[max_Q1D][max_Q1D];
196  double BGu[max_Q1D][max_Q1D];
197  for (int qx = 0; qx < Q1D; ++qx)
198  {
199  for (int qy = 0; qy < Q1D; ++qy)
200  {
201  GBu[qy][qx] = 0.0;
202  BGu[qy][qx] = 0.0;
203  for (int dy = 0; dy < D1D; ++dy)
204  {
205  const double bx = B(qy,dy);
206  const double gx = G(qy,dy);
207  GBu[qy][qx] += gx * Bu[dy][qx];
208  BGu[qy][qx] += bx * Gu[dy][qx];
209  }
210  }
211  }
212  // Calculate Dxy, xDy in plane
213  double DGu[max_Q1D][max_Q1D];
214  for (int qy = 0; qy < Q1D; ++qy)
215  {
216  for (int qx = 0; qx < Q1D; ++qx)
217  {
218  const double O1 = op(qx,qy,0,e);
219  const double O2 = op(qx,qy,1,e);
220 
221  const double gradX = BGu[qy][qx];
222  const double gradY = GBu[qy][qx];
223 
224  DGu[qy][qx] = (O1 * gradX) + (O2 * gradY);
225  }
226  }
227  double BDGu[max_D1D][max_Q1D];
228  for (int qx = 0; qx < Q1D; ++qx)
229  {
230  for (int dy = 0; dy < D1D; ++dy)
231  {
232  BDGu[dy][qx] = 0.0;
233  for (int qy = 0; qy < Q1D; ++qy)
234  {
235  const double w = Bt(dy,qy);
236  BDGu[dy][qx] += w * DGu[qy][qx];
237  }
238  }
239  }
240  for (int dx = 0; dx < D1D; ++dx)
241  {
242  for (int dy = 0; dy < D1D; ++dy)
243  {
244  double BBDGu = 0.0;
245  for (int qx = 0; qx < Q1D; ++qx)
246  {
247  const double w = Bt(dx,qx);
248  BBDGu += w * BDGu[dy][qx];
249  }
250  y(dx,dy,e) += BBDGu;
251  }
252  }
253  });
254 }
255 
256 // Optimized PA Convection Apply 2D kernel
257 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
258 void SmemPAConvectionApply2D(const int ne,
259  const Array<double> &b,
260  const Array<double> &g,
261  const Array<double> &bt,
262  const Array<double> &gt,
263  const Vector &op_,
264  const Vector &x_,
265  Vector &y_,
266  const int d1d = 0,
267  const int q1d = 0)
268 {
269  const int NE = ne;
270  const int D1D = T_D1D ? T_D1D : d1d;
271  const int Q1D = T_Q1D ? T_Q1D : q1d;
272  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
273  MFEM_VERIFY(D1D <= MAX_D1D, "");
274  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
275  auto B = Reshape(b.Read(), Q1D, D1D);
276  auto G = Reshape(g.Read(), Q1D, D1D);
277  auto Bt = Reshape(bt.Read(), D1D, Q1D);
278  auto op = Reshape(op_.Read(), Q1D, Q1D, 2, NE);
279  auto x = Reshape(x_.Read(), D1D, D1D, NE);
280  auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
281  MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
282  {
283  const int tidz = MFEM_THREAD_ID(z);
284  const int D1D = T_D1D ? T_D1D : d1d;
285  const int Q1D = T_Q1D ? T_Q1D : q1d;
286  // the following variables are evaluated at compile time
287  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
288  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
289  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
290  // constexpr int MDQ = (max_Q1D > max_D1D) ? max_Q1D : max_D1D;
291  MFEM_SHARED double u[NBZ][max_D1D][max_D1D];
292  MFEM_FOREACH_THREAD(dy,y,D1D)
293  {
294  MFEM_FOREACH_THREAD(dx,x,D1D)
295  {
296  // e is really equal to e+tidz
297  u[tidz][dy][dx] = x(dx,dy,e);
298  }
299  }
300  MFEM_SYNC_THREAD;
301  MFEM_SHARED double Bu[NBZ][max_D1D][max_Q1D];
302  MFEM_SHARED double Gu[NBZ][max_D1D][max_Q1D];
303  MFEM_FOREACH_THREAD(dy,y,D1D)
304  {
305  MFEM_FOREACH_THREAD(qx,x,Q1D)
306  {
307  Bu[tidz][dy][qx] = 0.0;
308  Gu[tidz][dy][qx] = 0.0;
309  for (int dx = 0; dx < D1D; ++dx)
310  {
311  const double bx = B(qx,dx);
312  const double gx = G(qx,dx);
313  const double x = u[tidz][dy][dx];
314  Bu[tidz][dy][qx] += bx * x;
315  Gu[tidz][dy][qx] += gx * x;
316  }
317  }
318  }
319  MFEM_SYNC_THREAD;
320  MFEM_SHARED double GBu[NBZ][max_Q1D][max_Q1D];
321  MFEM_SHARED double BGu[NBZ][max_Q1D][max_Q1D];
322  MFEM_FOREACH_THREAD(qx,x,Q1D)
323  {
324  MFEM_FOREACH_THREAD(qy,y,Q1D)
325  {
326  GBu[tidz][qy][qx] = 0.0;
327  BGu[tidz][qy][qx] = 0.0;
328  for (int dy = 0; dy < D1D; ++dy)
329  {
330  const double bx = B(qy,dy);
331  const double gx = G(qy,dy);
332  GBu[tidz][qy][qx] += gx * Bu[tidz][dy][qx];
333  BGu[tidz][qy][qx] += bx * Gu[tidz][dy][qx];
334  }
335  }
336  }
337  MFEM_SYNC_THREAD;
338  // Calculate Dxy, xDy in plane
339  MFEM_SHARED double DGu[NBZ][max_Q1D][max_Q1D];
340  MFEM_FOREACH_THREAD(qy,y,Q1D)
341  {
342  MFEM_FOREACH_THREAD(qx,x,Q1D)
343  {
344  const double O1 = op(qx,qy,0,e);
345  const double O2 = op(qx,qy,1,e);
346 
347  const double gradX = BGu[tidz][qy][qx];
348  const double gradY = GBu[tidz][qy][qx];
349 
350  DGu[tidz][qy][qx] = (O1 * gradX) + (O2 * gradY);
351  }
352  }
353  MFEM_SYNC_THREAD;
354  MFEM_SHARED double BDGu[NBZ][max_D1D][max_Q1D];
355  MFEM_FOREACH_THREAD(qx,x,Q1D)
356  {
357  MFEM_FOREACH_THREAD(dy,y,D1D)
358  {
359  BDGu[tidz][dy][qx] = 0.0;
360  for (int qy = 0; qy < Q1D; ++qy)
361  {
362  const double w = Bt(dy,qy);
363  BDGu[tidz][dy][qx] += w * DGu[tidz][qy][qx];
364  }
365  }
366  }
367  MFEM_SYNC_THREAD;
368  MFEM_FOREACH_THREAD(dx,x,D1D)
369  {
370  MFEM_FOREACH_THREAD(dy,y,D1D)
371  {
372  double BBDGu = 0.0;
373  for (int qx = 0; qx < Q1D; ++qx)
374  {
375  const double w = Bt(dx,qx);
376  BBDGu += w * BDGu[tidz][dy][qx];
377  }
378  y(dx,dy,e) += BBDGu;
379  }
380  }
381  });
382 }
383 
384 // PA Convection Apply 3D kernel
385 template<int T_D1D = 0, int T_Q1D = 0> static
386 void PAConvectionApply3D(const int ne,
387  const Array<double> &b,
388  const Array<double> &g,
389  const Array<double> &bt,
390  const Array<double> &gt,
391  const Vector &op_,
392  const Vector &x_,
393  Vector &y_,
394  const int d1d = 0,
395  const int q1d = 0)
396 {
397  const int NE = ne;
398  const int D1D = T_D1D ? T_D1D : d1d;
399  const int Q1D = T_Q1D ? T_Q1D : q1d;
400  MFEM_VERIFY(D1D <= MAX_D1D, "");
401  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
402  auto B = Reshape(b.Read(), Q1D, D1D);
403  auto G = Reshape(g.Read(), Q1D, D1D);
404  auto Bt = Reshape(bt.Read(), D1D, Q1D);
405  auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, 3, NE);
406  auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
407  auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
408  MFEM_FORALL(e, NE,
409  {
410  const int D1D = T_D1D ? T_D1D : d1d;
411  const int Q1D = T_Q1D ? T_Q1D : q1d;
412  // the following variables are evaluated at compile time
413  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
414  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
415 
416  double u[max_D1D][max_D1D][max_D1D];
417  for (int dz = 0; dz < D1D; ++dz)
418  {
419  for (int dy = 0; dy < D1D; ++dy)
420  {
421  for (int dx = 0; dx < D1D; ++dx)
422  {
423  u[dz][dy][dx] = x(dx,dy,dz,e);
424  }
425  }
426  }
427  double Bu[max_D1D][max_D1D][max_Q1D];
428  double Gu[max_D1D][max_D1D][max_Q1D];
429  for (int dz = 0; dz < D1D; ++dz)
430  {
431  for (int dy = 0; dy < D1D; ++dy)
432  {
433  for (int qx = 0; qx < Q1D; ++qx)
434  {
435  Bu[dz][dy][qx] = 0.0;
436  Gu[dz][dy][qx] = 0.0;
437  for (int dx = 0; dx < D1D; ++dx)
438  {
439  const double bx = B(qx,dx);
440  const double gx = G(qx,dx);
441  const double x = u[dz][dy][dx];
442  Bu[dz][dy][qx] += bx * x;
443  Gu[dz][dy][qx] += gx * x;
444  }
445  }
446  }
447  }
448  double BBu[max_D1D][max_Q1D][max_Q1D];
449  double GBu[max_D1D][max_Q1D][max_Q1D];
450  double BGu[max_D1D][max_Q1D][max_Q1D];
451  for (int dz = 0; dz < D1D; ++dz)
452  {
453  for (int qx = 0; qx < Q1D; ++qx)
454  {
455  for (int qy = 0; qy < Q1D; ++qy)
456  {
457  BBu[dz][qy][qx] = 0.0;
458  GBu[dz][qy][qx] = 0.0;
459  BGu[dz][qy][qx] = 0.0;
460  for (int dy = 0; dy < D1D; ++dy)
461  {
462  const double bx = B(qy,dy);
463  const double gx = G(qy,dy);
464  BBu[dz][qy][qx] += bx * Bu[dz][dy][qx];
465  GBu[dz][qy][qx] += gx * Bu[dz][dy][qx];
466  BGu[dz][qy][qx] += bx * Gu[dz][dy][qx];
467  }
468  }
469  }
470  }
471  double GBBu[max_Q1D][max_Q1D][max_Q1D];
472  double BGBu[max_Q1D][max_Q1D][max_Q1D];
473  double BBGu[max_Q1D][max_Q1D][max_Q1D];
474  for (int qx = 0; qx < Q1D; ++qx)
475  {
476  for (int qy = 0; qy < Q1D; ++qy)
477  {
478  for (int qz = 0; qz < Q1D; ++qz)
479  {
480  GBBu[qz][qy][qx] = 0.0;
481  BGBu[qz][qy][qx] = 0.0;
482  BBGu[qz][qy][qx] = 0.0;
483  for (int dz = 0; dz < D1D; ++dz)
484  {
485  const double bx = B(qz,dz);
486  const double gx = G(qz,dz);
487  GBBu[qz][qy][qx] += gx * BBu[dz][qy][qx];
488  BGBu[qz][qy][qx] += bx * GBu[dz][qy][qx];
489  BBGu[qz][qy][qx] += bx * BGu[dz][qy][qx];
490  }
491  }
492  }
493  }
494  // Calculate Dxy, xDy in plane
495  double DGu[max_Q1D][max_Q1D][max_Q1D];
496  for (int qz = 0; qz < Q1D; ++qz)
497  {
498  for (int qy = 0; qy < Q1D; ++qy)
499  {
500  for (int qx = 0; qx < Q1D; ++qx)
501  {
502  const double O1 = op(qx,qy,qz,0,e);
503  const double O2 = op(qx,qy,qz,1,e);
504  const double O3 = op(qx,qy,qz,2,e);
505 
506  const double gradX = BBGu[qz][qy][qx];
507  const double gradY = BGBu[qz][qy][qx];
508  const double gradZ = GBBu[qz][qy][qx];
509 
510  DGu[qz][qy][qx] = (O1 * gradX) + (O2 * gradY) + (O3 * gradZ);
511  }
512  }
513  }
514  double BDGu[max_D1D][max_Q1D][max_Q1D];
515  for (int qx = 0; qx < Q1D; ++qx)
516  {
517  for (int qy = 0; qy < Q1D; ++qy)
518  {
519  for (int dz = 0; dz < D1D; ++dz)
520  {
521  BDGu[dz][qy][qx] = 0.0;
522  for (int qz = 0; qz < Q1D; ++qz)
523  {
524  const double w = Bt(dz,qz);
525  BDGu[dz][qy][qx] += w * DGu[qz][qy][qx];
526  }
527  }
528  }
529  }
530  double BBDGu[max_D1D][max_D1D][max_Q1D];
531  for (int dz = 0; dz < D1D; ++dz)
532  {
533  for (int qx = 0; qx < Q1D; ++qx)
534  {
535  for (int dy = 0; dy < D1D; ++dy)
536  {
537  BBDGu[dz][dy][qx] = 0.0;
538  for (int qy = 0; qy < Q1D; ++qy)
539  {
540  const double w = Bt(dy,qy);
541  BBDGu[dz][dy][qx] += w * BDGu[dz][qy][qx];
542  }
543  }
544  }
545  }
546  for (int dz = 0; dz < D1D; ++dz)
547  {
548  for (int dy = 0; dy < D1D; ++dy)
549  {
550  for (int dx = 0; dx < D1D; ++dx)
551  {
552  double BBBDGu = 0.0;
553  for (int qx = 0; qx < Q1D; ++qx)
554  {
555  const double w = Bt(dx,qx);
556  BBBDGu += w * BBDGu[dz][dy][qx];
557  }
558  y(dx,dy,dz,e) += BBBDGu;
559  }
560  }
561  }
562  });
563 }
564 
565 // Optimized PA Convection Apply 3D kernel
566 template<int T_D1D = 0, int T_Q1D = 0> static
567 void SmemPAConvectionApply3D(const int ne,
568  const Array<double> &b,
569  const Array<double> &g,
570  const Array<double> &bt,
571  const Array<double> &gt,
572  const Vector &op_,
573  const Vector &x_,
574  Vector &y_,
575  const int d1d = 0,
576  const int q1d = 0)
577 {
578  const int NE = ne;
579  const int D1D = T_D1D ? T_D1D : d1d;
580  const int Q1D = T_Q1D ? T_Q1D : q1d;
581  MFEM_VERIFY(D1D <= MAX_D1D, "");
582  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
583  auto B = Reshape(b.Read(), Q1D, D1D);
584  auto G = Reshape(g.Read(), Q1D, D1D);
585  auto Bt = Reshape(bt.Read(), D1D, Q1D);
586  auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, 3, NE);
587  auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
588  auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
589  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
590  {
591  const int D1D = T_D1D ? T_D1D : d1d;
592  const int Q1D = T_Q1D ? T_Q1D : q1d;
593  // the following variables are evaluated at compile time
594  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
595  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
596  constexpr int max_DQ = (max_Q1D > max_D1D) ? max_Q1D : max_D1D;
597  MFEM_SHARED double sm0[max_DQ*max_DQ*max_DQ];
598  MFEM_SHARED double sm1[max_DQ*max_DQ*max_DQ];
599  MFEM_SHARED double sm2[max_DQ*max_DQ*max_DQ];
600  MFEM_SHARED double sm3[max_DQ*max_DQ*max_DQ];
601  MFEM_SHARED double sm4[max_DQ*max_DQ*max_DQ];
602  MFEM_SHARED double sm5[max_DQ*max_DQ*max_DQ];
603 
604  double (*u)[max_D1D][max_D1D] = (double (*)[max_D1D][max_D1D]) sm0;
605  MFEM_FOREACH_THREAD(dz,z,D1D)
606  {
607  MFEM_FOREACH_THREAD(dy,y,D1D)
608  {
609  MFEM_FOREACH_THREAD(dx,x,D1D)
610  {
611  u[dz][dy][dx] = x(dx,dy,dz,e);
612  }
613  }
614  }
615  MFEM_SYNC_THREAD;
616  double (*Bu)[max_D1D][max_Q1D] = (double (*)[max_D1D][max_Q1D])sm1;
617  double (*Gu)[max_D1D][max_Q1D] = (double (*)[max_D1D][max_Q1D])sm2;
618  MFEM_FOREACH_THREAD(dz,z,D1D)
619  {
620  MFEM_FOREACH_THREAD(dy,y,D1D)
621  {
622  MFEM_FOREACH_THREAD(qx,x,Q1D)
623  {
624  double Bu_ = 0.0;
625  double Gu_ = 0.0;
626  for (int dx = 0; dx < D1D; ++dx)
627  {
628  const double bx = B(qx,dx);
629  const double gx = G(qx,dx);
630  const double x = u[dz][dy][dx];
631  Bu_ += bx * x;
632  Gu_ += gx * x;
633  }
634  Bu[dz][dy][qx] = Bu_;
635  Gu[dz][dy][qx] = Gu_;
636  }
637  }
638  }
639  MFEM_SYNC_THREAD;
640  double (*BBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm3;
641  double (*GBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm4;
642  double (*BGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm5;
643  MFEM_FOREACH_THREAD(dz,z,D1D)
644  {
645  MFEM_FOREACH_THREAD(qx,x,Q1D)
646  {
647  MFEM_FOREACH_THREAD(qy,y,Q1D)
648  {
649  double BBu_ = 0.0;
650  double GBu_ = 0.0;
651  double BGu_ = 0.0;
652  for (int dy = 0; dy < D1D; ++dy)
653  {
654  const double bx = B(qy,dy);
655  const double gx = G(qy,dy);
656  BBu_ += bx * Bu[dz][dy][qx];
657  GBu_ += gx * Bu[dz][dy][qx];
658  BGu_ += bx * Gu[dz][dy][qx];
659  }
660  BBu[dz][qy][qx] = BBu_;
661  GBu[dz][qy][qx] = GBu_;
662  BGu[dz][qy][qx] = BGu_;
663  }
664  }
665  }
666  MFEM_SYNC_THREAD;
667  double (*GBBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm0;
668  double (*BGBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm1;
669  double (*BBGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm2;
670  MFEM_FOREACH_THREAD(qx,x,Q1D)
671  {
672  MFEM_FOREACH_THREAD(qy,y,Q1D)
673  {
674  MFEM_FOREACH_THREAD(qz,z,Q1D)
675  {
676  double GBBu_ = 0.0;
677  double BGBu_ = 0.0;
678  double BBGu_ = 0.0;
679  for (int dz = 0; dz < D1D; ++dz)
680  {
681  const double bx = B(qz,dz);
682  const double gx = G(qz,dz);
683  GBBu_ += gx * BBu[dz][qy][qx];
684  BGBu_ += bx * GBu[dz][qy][qx];
685  BBGu_ += bx * BGu[dz][qy][qx];
686  }
687  GBBu[qz][qy][qx] = GBBu_;
688  BGBu[qz][qy][qx] = BGBu_;
689  BBGu[qz][qy][qx] = BBGu_;
690  }
691  }
692  }
693  MFEM_SYNC_THREAD;
694  double (*DGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm3;
695  MFEM_FOREACH_THREAD(qz,z,Q1D)
696  {
697  MFEM_FOREACH_THREAD(qy,y,Q1D)
698  {
699  MFEM_FOREACH_THREAD(qx,x,Q1D)
700  {
701  const double O1 = op(qx,qy,qz,0,e);
702  const double O2 = op(qx,qy,qz,1,e);
703  const double O3 = op(qx,qy,qz,2,e);
704 
705  const double gradX = BBGu[qz][qy][qx];
706  const double gradY = BGBu[qz][qy][qx];
707  const double gradZ = GBBu[qz][qy][qx];
708 
709  DGu[qz][qy][qx] = (O1 * gradX) + (O2 * gradY) + (O3 * gradZ);
710  }
711  }
712  }
713  MFEM_SYNC_THREAD;
714  double (*BDGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm4;
715  MFEM_FOREACH_THREAD(qx,x,Q1D)
716  {
717  MFEM_FOREACH_THREAD(qy,y,Q1D)
718  {
719  MFEM_FOREACH_THREAD(dz,z,D1D)
720  {
721  double BDGu_ = 0.0;
722  for (int qz = 0; qz < Q1D; ++qz)
723  {
724  const double w = Bt(dz,qz);
725  BDGu_ += w * DGu[qz][qy][qx];
726  }
727  BDGu[dz][qy][qx] = BDGu_;
728  }
729  }
730  }
731  MFEM_SYNC_THREAD;
732  double (*BBDGu)[max_D1D][max_Q1D] = (double (*)[max_D1D][max_Q1D])sm5;
733  MFEM_FOREACH_THREAD(dz,z,D1D)
734  {
735  MFEM_FOREACH_THREAD(qx,x,Q1D)
736  {
737  MFEM_FOREACH_THREAD(dy,y,D1D)
738  {
739  double BBDGu_ = 0.0;
740  for (int qy = 0; qy < Q1D; ++qy)
741  {
742  const double w = Bt(dy,qy);
743  BBDGu_ += w * BDGu[dz][qy][qx];
744  }
745  BBDGu[dz][dy][qx] = BBDGu_;
746  }
747  }
748  }
749  MFEM_SYNC_THREAD;
750  MFEM_FOREACH_THREAD(dz,z,D1D)
751  {
752  MFEM_FOREACH_THREAD(dy,y,D1D)
753  {
754  MFEM_FOREACH_THREAD(dx,x,D1D)
755  {
756  double BBBDGu = 0.0;
757  for (int qx = 0; qx < Q1D; ++qx)
758  {
759  const double w = Bt(dx,qx);
760  BBBDGu += w * BBDGu[dz][dy][qx];
761  }
762  y(dx,dy,dz,e) = BBBDGu;
763  }
764  }
765  }
766  });
767 }
768 
770 {
771  const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
773  // Assumes tensor-product elements
774  Mesh *mesh = fes.GetMesh();
775  const FiniteElement &el = *fes.GetFE(0);
777  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, Trans);
778  if (DeviceCanUseCeed())
779  {
780  delete ceedOp;
781  ceedOp = new ceed::PAConvectionIntegrator(fes, *ir, Q, alpha);
782  return;
783  }
784  const int dims = el.GetDim();
785  const int symmDims = dims;
786  const int nq = ir->GetNPoints();
787  dim = mesh->Dimension();
788  ne = fes.GetNE();
790  maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
791  dofs1D = maps->ndof;
792  quad1D = maps->nqpt;
793  pa_data.SetSize(symmDims * nq * ne, mt);
794  Vector vel;
795  if (VectorConstantCoefficient *cQ =
796  dynamic_cast<VectorConstantCoefficient*>(Q))
797  {
798  vel = cQ->GetVec();
799  }
800  else if (VectorGridFunctionCoefficient *vgfQ =
801  dynamic_cast<VectorGridFunctionCoefficient*>(Q))
802  {
803  vel.SetSize(dim * nq * ne, mt);
804 
805  const GridFunction *gf = vgfQ->GetGridFunction();
806  const FiniteElementSpace &gf_fes = *gf->FESpace();
807  const QuadratureInterpolator *qi(gf_fes.GetQuadratureInterpolator(*ir));
808  const bool use_tensor_products = UsesTensorBasis(gf_fes);
809  const ElementDofOrdering ordering = use_tensor_products ?
812  const Operator *R = gf_fes.GetElementRestriction(ordering);
813 
814  Vector xe(R->Height(), mt);
815  xe.UseDevice(true);
816 
817  R->Mult(*gf, xe);
818  qi->SetOutputLayout(QVectorLayout::byVDIM);
819  qi->DisableTensorProducts(!use_tensor_products);
820  qi->Values(xe,vel);
821  }
823  dynamic_cast<VectorQuadratureFunctionCoefficient*>(Q))
824  {
825  const QuadratureFunction &qFun = cQ->GetQuadFunction();
826  MFEM_VERIFY(qFun.Size() == dim * nq * ne,
827  "Incompatible QuadratureFunction dimension \n");
828 
829  MFEM_VERIFY(ir == &qFun.GetSpace()->GetElementIntRule(0),
830  "IntegrationRule used within integrator and in"
831  " QuadratureFunction appear to be different");
832 
833  qFun.Read();
834  vel.MakeRef(const_cast<QuadratureFunction &>(qFun),0);
835  }
836  else
837  {
838  vel.SetSize(dim * nq * ne);
839  auto C = Reshape(vel.HostWrite(), dim, nq, ne);
840  DenseMatrix Q_ir;
841  for (int e = 0; e < ne; ++e)
842  {
844  Q->Eval(Q_ir, T, *ir);
845  for (int q = 0; q < nq; ++q)
846  {
847  for (int i = 0; i < dim; ++i)
848  {
849  C(i,q,e) = Q_ir(i,q);
850  }
851  }
852  }
853  }
854  PAConvectionSetup(dim, nq, ne, ir->GetWeights(), geom->J,
855  vel, alpha, pa_data);
856 }
857 
858 static void PAConvectionApply(const int dim,
859  const int D1D,
860  const int Q1D,
861  const int NE,
862  const Array<double> &B,
863  const Array<double> &G,
864  const Array<double> &Bt,
865  const Array<double> &Gt,
866  const Vector &op,
867  const Vector &x,
868  Vector &y)
869 {
870  if (dim == 2)
871  {
872  switch ((D1D << 4 ) | Q1D)
873  {
874  case 0x22: return SmemPAConvectionApply2D<2,2,8>(NE,B,G,Bt,Gt,op,x,y);
875  case 0x33: return SmemPAConvectionApply2D<3,3,4>(NE,B,G,Bt,Gt,op,x,y);
876  case 0x34: return SmemPAConvectionApply2D<3,4,4>(NE,B,G,Bt,Gt,op,x,y);
877  case 0x44: return SmemPAConvectionApply2D<4,4,4>(NE,B,G,Bt,Gt,op,x,y);
878  case 0x46: return SmemPAConvectionApply2D<4,6,4>(NE,B,G,Bt,Gt,op,x,y);
879  case 0x55: return SmemPAConvectionApply2D<5,5,2>(NE,B,G,Bt,Gt,op,x,y);
880  case 0x58: return SmemPAConvectionApply2D<5,8,2>(NE,B,G,Bt,Gt,op,x,y);
881  case 0x66: return SmemPAConvectionApply2D<6,6,1>(NE,B,G,Bt,Gt,op,x,y);
882  case 0x77: return SmemPAConvectionApply2D<7,7,1>(NE,B,G,Bt,Gt,op,x,y);
883  case 0x88: return SmemPAConvectionApply2D<8,8,1>(NE,B,G,Bt,Gt,op,x,y);
884  case 0x99: return SmemPAConvectionApply2D<9,9,1>(NE,B,G,Bt,Gt,op,x,y);
885  default: return PAConvectionApply2D(NE,B,G,Bt,Gt,op,x,y,D1D,Q1D);
886  }
887  }
888  else if (dim == 3)
889  {
890  switch ((D1D << 4 ) | Q1D)
891  {
892  case 0x23: return SmemPAConvectionApply3D<2,3>(NE,B,G,Bt,Gt,op,x,y);
893  case 0x24: return SmemPAConvectionApply3D<2,4>(NE,B,G,Bt,Gt,op,x,y);
894  case 0x26: return SmemPAConvectionApply3D<2,6>(NE,B,G,Bt,Gt,op,x,y);
895  case 0x34: return SmemPAConvectionApply3D<3,4>(NE,B,G,Bt,Gt,op,x,y);
896  case 0x35: return SmemPAConvectionApply3D<3,5>(NE,B,G,Bt,Gt,op,x,y);
897  case 0x45: return SmemPAConvectionApply3D<4,5>(NE,B,G,Bt,Gt,op,x,y);
898  case 0x48: return SmemPAConvectionApply3D<4,8>(NE,B,G,Bt,Gt,op,x,y);
899  case 0x56: return SmemPAConvectionApply3D<5,6>(NE,B,G,Bt,Gt,op,x,y);
900  case 0x67: return SmemPAConvectionApply3D<6,7>(NE,B,G,Bt,Gt,op,x,y);
901  case 0x78: return SmemPAConvectionApply3D<7,8>(NE,B,G,Bt,Gt,op,x,y);
902  case 0x89: return SmemPAConvectionApply3D<8,9>(NE,B,G,Bt,Gt,op,x,y);
903  default: return PAConvectionApply3D(NE,B,G,Bt,Gt,op,x,y,D1D,Q1D);
904  }
905  }
906  MFEM_ABORT("Unknown kernel.");
907 }
908 
909 // PA Convection Apply kernel
911 {
912  if (DeviceCanUseCeed())
913  {
914  ceedOp->AddMult(x, y);
915  }
916  else
917  {
918  PAConvectionApply(dim, dofs1D, quad1D, ne,
919  maps->B, maps->G, maps->Bt, maps->Gt,
920  pa_data, x, y);
921  }
922 }
923 
925 {
926  if (DeviceCanUseCeed())
927  {
928  ceedOp->GetDiagonal(diag);
929  }
930  else
931  {
932  MFEM_ABORT("AssembleDiagonalPA not yet implemented for"
933  " ConvectionIntegrator.");
934  }
935 }
936 
937 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe.hpp:243
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:317
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe.hpp:165
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:784
const IntegrationRule & GetElementIntRule(int idx) const
Get the IntegrationRule associated with mesh element idx.
Definition: fespace.hpp:950
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
const GeometricFactors * geom
Not owned.
Vector quadrature function coefficient which requires that the quadrature rules used for this vector ...
Vector coefficient that is constant in space and time.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe.hpp:177
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:438
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
Array< double > Gt
Transpose of G.
Definition: fe.hpp:216
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe.hpp:173
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:83
constexpr int DIM
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:957
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:108
virtual void AssembleDiagonalPA(Vector &diag)
Assemble diagonal and add it to Vector diag.
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1195
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:136
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:567
Vector J
Jacobians of the element transformations at all quadrature points.
Definition: mesh.hpp:1596
Native ordering as defined by the FiniteElement.
const int MAX_Q1D
Definition: forall.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
double b
Definition: lissajous.cpp:42
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
QuadratureSpace * GetSpace() const
Get the associated QuadratureSpace.
Definition: gridfunc.hpp:775
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
int Dimension() const
Definition: mesh.hpp:911
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:600
Array< double > Bt
Transpose of B.
Definition: fe.hpp:194
virtual void AssemblePA(const FiniteElementSpace &)
Method defining partial assembly.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
static const IntegrationRule & GetRule(const FiniteElement &el, ElementTransformation &Trans)
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:30
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe.hpp:188
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe.cpp:367
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition: fespace.cpp:1256
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:65
int dim
Definition: ex24.cpp:53
const int MAX_D1D
Definition: forall.hpp:27
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition: fe.hpp:209
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2388
constexpr int SDIM
Lexicographic ordering for tensor-product FiniteElements.
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:577
Vector coefficient defined by a vector GridFunction.
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Abstract operator.
Definition: operator.hpp:24
Class representing a function through its values (scalar or vector) at quadrature points...
Definition: gridfunc.hpp:734
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:426
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
const DofToQuad * maps
Not owned.
VectorCoefficient * Q