MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
nonlininteg_vecconvection_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/forall.hpp"
15#include "../../linalg/dtensor.hpp"
16#include "../kernels.hpp"
17#include "../nonlininteg.hpp"
18
19namespace mfem
20{
21
22/// \cond DO_NOT_DOCUMENT
23
24namespace internal
25{
26
27// PA Convection NL 2D kernel
28template<int T_D1D = 0, int T_Q1D = 0>
29inline void SmemPAConvectionNLApply2D(const int NE,
30 const real_t *b,
31 const real_t *g,
32 const real_t *a,
33 const real_t *x,
34 real_t *y,
35 const int d1d = 0,
36 const int q1d = 0)
37{
38 static constexpr int VDIM = 2, DIM = 2;
39 const int D1D = T_D1D ? T_D1D : d1d;
40 const int Q1D = T_Q1D ? T_Q1D : q1d;
41
42 const auto B = Reshape(b, Q1D, D1D);
43 const auto G = Reshape(g, Q1D, D1D);
44 const auto A = Reshape(a, VDIM, DIM, Q1D, Q1D, NE);
45 const auto X = Reshape(x, D1D, D1D, VDIM, NE);
46 auto Y = Reshape(y, D1D, D1D, VDIM, NE);
47
48 mfem::forall_2D<T_Q1D * T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
49 {
50 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
51 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
52
53 MFEM_SHARED real_t smem[MQ1][MQ1], sB[MD1][MQ1], sG[MD1][MQ1];
54
55 kernels::internal::vd_regs2d_t<VDIM, DIM, MQ1> g0, g1;
56 kernels::internal::v_regs2d_t<VDIM, MQ1> r0, r1;
57 kernels::internal::v_regs2d_t<VDIM, MQ1> s0, s1;
58
59 kernels::internal::LoadMatrix(D1D, Q1D, B, sB);
60 kernels::internal::LoadMatrix(D1D, Q1D, G, sG);
61
62 kernels::internal::LoadDofs2d(e, D1D, X, r0);
63 kernels::internal::Eval2d(D1D, Q1D, smem, sB, r0, r1); // u vector-value
64 kernels::internal::LoadDofs2d(e, D1D, X, g0);
65 kernels::internal::Grad2d(D1D, Q1D, smem, sB, sG, g0, g1); // u vector-gradient
66
67 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
68 {
69 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
70 {
71 const future::tensor<real_t, 2> U =
72 {
73 r1[0][qy][qx], r1[1][qy][qx]
74 };
75 const future::tensor<real_t, 2,2> gradU = {{
76 {g1[0][0][qy][qx], g1[1][0][qy][qx]},
77 {g1[0][1][qy][qx], g1[1][1][qy][qx]},
78 }
79 };
80 const future::tensor<real_t, 2,2> Q = {{
81 {A(0,0,qx,qy,e), A(1,0,qx,qy,e)},
82 {A(0,1,qx,qy,e), A(1,1,qx,qy,e)},
83 }
84 };
85 const future::tensor<real_t, 2> conv = transpose(gradU) * (Q * U);
86 s0[0][qy][qx] = conv[0];
87 s0[1][qy][qx] = conv[1];
88 }
89 }
90 MFEM_SYNC_THREAD;
91 kernels::internal::EvalTranspose2d(D1D, Q1D, smem, sB, s0, s1);
92 kernels::internal::WriteDofs2d(e, D1D, s1, Y);
93 });
94}
95
96// PA Convection NL 3D kernel
97template<int T_D1D = 0, int T_Q1D = 0>
98inline void SmemPAConvectionNLApply3D(const int NE,
99 const real_t *b,
100 const real_t *g,
101 const real_t *a,
102 const real_t *x,
103 real_t *y,
104 const int d1d = 0,
105 const int q1d = 0)
106{
107 static constexpr int VDIM = 3, DIM = 3;
108 const int D1D = T_D1D ? T_D1D : d1d;
109 const int Q1D = T_Q1D ? T_Q1D : q1d;
110
111 const auto B = Reshape(b, Q1D, D1D);
112 const auto G = Reshape(g, Q1D, D1D);
113 const auto A = Reshape(a, VDIM, DIM, Q1D, Q1D, Q1D, NE);
114 const auto X = Reshape(x, D1D, D1D, D1D, VDIM, NE);
115 auto Y = Reshape(y, D1D, D1D, D1D, VDIM, NE);
116
117 mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
118 {
119 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
120 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
121
122 MFEM_SHARED real_t smem[MQ1][MQ1], sB[MD1][MQ1], sG[MD1][MQ1];
123
124 kernels::internal::vd_regs3d_t<VDIM, DIM, MQ1> g0, g1;
125 kernels::internal::v_regs3d_t<VDIM, MQ1> r0, r1;
126 kernels::internal::v_regs3d_t<VDIM, MQ1> s0, s1;
127
128 kernels::internal::LoadMatrix(D1D, Q1D, B, sB);
129 kernels::internal::LoadMatrix(D1D, Q1D, G, sG);
130
131 kernels::internal::LoadDofs3d(e, D1D, X, r0);
132 kernels::internal::Eval3d(D1D, Q1D, smem, sB, r0, r1); // u vector-value
133 kernels::internal::LoadDofs3d(e, D1D, X, g0);
134 kernels::internal::Grad3d(D1D, Q1D, smem, sB, sG, g0, g1); // u vector-gradient
135
136 for (int qz = 0; qz < Q1D; qz++)
137 {
138 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
139 {
140 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
141 {
142 const future::tensor<real_t, 3> U =
143 {
144 r1[0][qz][qy][qx], r1[1][qz][qy][qx], r1[2][qz][qy][qx]
145 };
146 const future::tensor<real_t, 3,3> gradU = {{
147 {g1[0][0][qz][qy][qx], g1[1][0][qz][qy][qx], g1[2][0][qz][qy][qx]},
148 {g1[0][1][qz][qy][qx], g1[1][1][qz][qy][qx], g1[2][1][qz][qy][qx]},
149 {g1[0][2][qz][qy][qx], g1[1][2][qz][qy][qx], g1[2][2][qz][qy][qx]}
150 }
151 };
152 const future::tensor<real_t, 3,3> Q = {{
153 {A(0,0,qx,qy,qz,e), A(1,0,qx,qy,qz,e), A(2,0,qx,qy,qz,e)},
154 {A(0,1,qx,qy,qz,e), A(1,1,qx,qy,qz,e), A(2,1,qx,qy,qz,e)},
155 {A(0,2,qx,qy,qz,e), A(1,2,qx,qy,qz,e), A(2,2,qx,qy,qz,e)}
156 }
157 };
158 const future::tensor<real_t, 3> conv = transpose(gradU) * (Q * U);
159 s0[0][qz][qy][qx] = conv[0];
160 s0[1][qz][qy][qx] = conv[1];
161 s0[2][qz][qy][qx] = conv[2];
162 }
163 }
164 }
165 MFEM_SYNC_THREAD;
166 kernels::internal::EvalTranspose3d(D1D, Q1D, smem, sB, s0, s1);
167 kernels::internal::WriteDofs3d(e, D1D, s1, Y);
168 });
169}
170
171} // namespace internal
172
173template<int DIM, int T_D1D, int T_Q1D>
174VectorConvectionNLFIntegrator::AddMultPAType
175VectorConvectionNLFIntegrator::AddMultPAKernels::Kernel()
176{
177 static_assert(T_D1D <= T_Q1D, "d1d > q1d is not supported");
178 if constexpr (DIM == 2)
179 {
180 return internal::SmemPAConvectionNLApply2D<T_D1D, T_Q1D>;
181 }
182 else if constexpr (DIM == 3)
183 {
184 return internal::SmemPAConvectionNLApply3D<T_D1D, T_Q1D>;
185 }
186 MFEM_ABORT("Unsupported kernel");
187}
188
189inline VectorConvectionNLFIntegrator::AddMultPAType
190VectorConvectionNLFIntegrator::AddMultPAKernels::Fallback
191(int dim, int d1d, int q1d)
192{
193 MFEM_VERIFY(d1d <= q1d, "d1d > q1d is not supported");
194 MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
195 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
196 if (dim == 2)
197 {
198 return internal::SmemPAConvectionNLApply2D<>;
199 }
200 else if (dim == 3)
201 {
202 return internal::SmemPAConvectionNLApply3D<>;
203 }
204 MFEM_ABORT("Unsupported kernel");
205}
206
207/// \endcond DO_NOT_DOCUMENT
208
209} // namespace mfem
mfem::VectorConvectionNLFIntegrator::AddMultPAType
void(*)(const int ne, const real_t *B, const real_t *G, const real_t *A, const real_t *x, real_t *y, const int d1d, const int q1d) AddMultPAType
Definition nonlininteg.hpp:421
dim
int dim
Definition ex24.cpp:53
b
real_t b
Definition lissajous.cpp:42
a
real_t a
Definition lissajous.cpp:41
DIM
constexpr int DIM
Definition minimal-surface.cpp:71
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
mfem::future::transpose
MFEM_HOST_DEVICE tensor< T, n, m > transpose(const tensor< T, m, n > &A)
Returns the transpose of the matrix.
Definition tensor.hpp:1388
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
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0