MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
nonlininteg_vecconvection_pa_grad.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
27template<int T_D1D = 0, int T_Q1D = 0>
28inline void SmemPAConvectionNLGradApply2D(const int ne,
29 const real_t *b,
30 const real_t *g,
31 const real_t *a,
32 const real_t *u,
33 const real_t *du,
34 real_t *y,
35 const int d1d,
36 const int q1d)
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 A = Reshape(a, VDIM, DIM, Q1D, Q1D, ne);
43 const auto U = Reshape(u, D1D, D1D, VDIM, ne);
44 const auto dU = Reshape(du, D1D, D1D, VDIM, ne);
45 auto Y = Reshape(y, D1D, D1D, VDIM, ne);
46
47 mfem::forall_2D<T_Q1D * T_Q1D>(ne, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
48 {
49 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
50 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
51
52 MFEM_SHARED real_t smem[MQ1][MQ1];
53 MFEM_SHARED real_t sB[MD1][MQ1], sG[MD1][MQ1];
54
55 kernels::internal::vd_regs2d_t<VDIM, DIM, MQ1> g0, g1, g2;
56 kernels::internal::v_regs2d_t<DIM, MQ1> r0, r1, r2;
57
58 kernels::internal::LoadMatrix(D1D, Q1D, b, sB);
59 kernels::internal::LoadMatrix(D1D, Q1D, g, sG);
60
61 kernels::internal::LoadDofs2d(e, D1D, dU, g0);
62 kernels::internal::Grad2d(D1D, Q1D, smem, sB, sG, g0, g1); // δu gradient
63
64 kernels::internal::LoadDofs2d(e, D1D, U, r0);
65 kernels::internal::Eval2d(D1D, Q1D, smem, sB, r0, r2); // u value
66
67 kernels::internal::LoadDofs2d(e, D1D, dU, r0);
68 kernels::internal::Eval2d(D1D, Q1D, smem, sB, r0, r1); // δu value
69
70 kernels::internal::LoadDofs2d(e, D1D, U, g0);
71 kernels::internal::Grad2d(D1D, Q1D, smem, sB, sG, g0, g2); // u gradient
72
73 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
74 {
75 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
76 {
77 // First part of the Jacobian: u·∇δu
78 const future::tensor<real_t, DIM> u_val =
79 {
80 r2[0][qy][qx], r2[1][qy][qx]
81 };
82 const future::tensor<real_t, VDIM, DIM> Q_adj =
83 {
84 { { A(0, 0, qx, qy, e), A(1, 0, qx, qy, e) },
85 { A(0, 1, qx, qy, e), A(1, 1, qx, qy, e) }
86 }
87 };
88 const future::tensor<real_t, VDIM, DIM> grad_dU =
89 {
90 { { g1[0][0][qy][qx], g1[1][0][qy][qx] },
91 { g1[0][1][qy][qx], g1[1][1][qy][qx] }
92 }
93 };
94 const auto one = transpose(grad_dU) * (Q_adj * u_val);
95
96 // Second part of the Jacobian: δu·∇u
97 const future::tensor<real_t, DIM> du_val =
98 {
99 r1[0][qy][qx], r1[1][qy][qx]
100 };
101 const future::tensor<real_t, VDIM, DIM> grad_U =
102 {
103 { { g2[0][0][qy][qx], g2[1][0][qy][qx] },
104 { g2[0][1][qy][qx], g2[1][1][qy][qx] }
105 }
106 };
107 const auto two = transpose(grad_U) * (Q_adj * du_val);
108
109 // u⋅∇δu + δu⋅∇u
110 r0[0][qy][qx] = one[0] + two[0];
111 r0[1][qy][qx] = one[1] + two[1];
112 }
113 }
114 MFEM_SYNC_THREAD;
115 kernels::internal::EvalTranspose2d(D1D, Q1D, smem, sB, r0, r1);
116 kernels::internal::WriteDofs2d(e, D1D, r1, Y);
117 });
118}
119
120template<int T_D1D = 0, int T_Q1D = 0>
121inline void SmemPAConvectionNLGradApply3D(const int ne,
122 const real_t *b,
123 const real_t *g,
124 const real_t *a,
125 const real_t *u,
126 const real_t *du,
127 real_t *y,
128 const int d1d,
129 const int q1d)
130{
131 static constexpr int VDIM = 3, DIM = 3;
132 const int D1D = T_D1D ? T_D1D : d1d;
133 const int Q1D = T_Q1D ? T_Q1D : q1d;
134
135 const auto A = Reshape(a, VDIM, DIM, Q1D, Q1D, Q1D, ne);
136 const auto U = Reshape(u, D1D, D1D, D1D, VDIM, ne);
137 const auto dU = Reshape(du, D1D, D1D, D1D, VDIM, ne);
138 auto Y = Reshape(y, D1D, D1D, D1D, VDIM, ne);
139
140 mfem::forall_2D<T_Q1D * T_Q1D>(ne, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
141 {
142 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
143 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
144
145 MFEM_SHARED real_t smem[MQ1][MQ1];
146 MFEM_SHARED real_t sB[MD1][MQ1], sG[MD1][MQ1];
147
148 kernels::internal::v_regs3d_t<VDIM, MQ1> r0, r1, r2;
149 kernels::internal::vd_regs3d_t<VDIM, DIM, MQ1> g0, g1, g2;
150
151 kernels::internal::LoadMatrix(D1D, Q1D, b, sB);
152 kernels::internal::LoadMatrix(D1D, Q1D, g, sG);
153
154 kernels::internal::LoadDofs3d(e, D1D, dU, g0);
155 kernels::internal::Grad3d(D1D, Q1D, smem, sB, sG, g0, g1); // δu gradient
156
157 kernels::internal::LoadDofs3d(e, D1D, U, r0);
158 kernels::internal::Eval3d(D1D, Q1D, smem, sB, r0, r2); // u value
159
160 kernels::internal::LoadDofs3d(e, D1D, dU, r0);
161 kernels::internal::Eval3d(D1D, Q1D, smem, sB, r0, r1); // δu value
162
163 kernels::internal::LoadDofs3d(e, D1D, U, g0);
164 kernels::internal::Grad3d(D1D, Q1D, smem, sB, sG, g0, g2); // u gradient
165
166 for (int qz = 0; qz < Q1D; qz++)
167 {
168 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
169 {
170 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
171 {
172 // First part of the Jacobian: u·∇δu
173 const future::tensor<real_t, DIM> u_val =
174 {
175 r2[0][qz][qy][qx],
176 r2[1][qz][qy][qx],
177 r2[2][qz][qy][qx]
178 };
179 const future::tensor<real_t, VDIM, DIM> Q_adj = {{
180 {A(0,0,qx,qy,qz,e), A(1,0,qx,qy,qz,e), A(2,0,qx,qy,qz,e)},
181 {A(0,1,qx,qy,qz,e), A(1,1,qx,qy,qz,e), A(2,1,qx,qy,qz,e)},
182 {A(0,2,qx,qy,qz,e), A(1,2,qx,qy,qz,e), A(2,2,qx,qy,qz,e)}
183 }
184 };
185 const future::tensor<real_t, DIM, DIM> grad_dU = {{
186 {g1[0][0][qz][qy][qx], g1[1][0][qz][qy][qx], g1[2][0][qz][qy][qx]},
187 {g1[0][1][qz][qy][qx], g1[1][1][qz][qy][qx], g1[2][1][qz][qy][qx]},
188 {g1[0][2][qz][qy][qx], g1[1][2][qz][qy][qx], g1[2][2][qz][qy][qx]}
189 }
190 };
191 const auto one = transpose(grad_dU) * (Q_adj * u_val);
192
193 // Second part of the Jacobian: δu·∇u
194 const future::tensor<real_t, DIM> du_val =
195 {
196 r1[0][qz][qy][qx], r1[1][qz][qy][qx], r1[2][qz][qy][qx]
197 };
198 const future::tensor<real_t, VDIM, DIM> grad_U = {{
199 {g2[0][0][qz][qy][qx], g2[1][0][qz][qy][qx], g2[2][0][qz][qy][qx]},
200 {g2[0][1][qz][qy][qx], g2[1][1][qz][qy][qx], g2[2][1][qz][qy][qx]},
201 {g2[0][2][qz][qy][qx], g2[1][2][qz][qy][qx], g2[2][2][qz][qy][qx]}
202 }
203 };
204 const auto two = transpose(grad_U) * (Q_adj * du_val);
205
206 // u⋅∇δu + δu⋅∇u
207 r0[0][qz][qy][qx] = one[0] + two[0];
208 r0[1][qz][qy][qx] = one[1] + two[1];
209 r0[2][qz][qy][qx] = one[2] + two[2];
210 }
211 }
212 }
213 MFEM_SYNC_THREAD;
214 kernels::internal::EvalTranspose3d(D1D, Q1D, smem, sB, r0, r1);
215 kernels::internal::WriteDofs3d(e, D1D, r1, Y);
216 });
217}
218
219} // namespace internal
220
221template<int T_D1D, int T_Q1D>
222VectorConvectionNLFIntegrator::AddMultGradPAType
223VectorConvectionNLFIntegrator::AddMultGradPA2D::Kernel()
224{
225 static_assert(T_D1D <= T_Q1D, "d1d > q1d is not supported");
226 return internal::SmemPAConvectionNLGradApply2D<T_D1D, T_Q1D>;
227}
228
229inline VectorConvectionNLFIntegrator::AddMultGradPAType
230VectorConvectionNLFIntegrator::AddMultGradPA2D::Fallback(int d1d, int q1d)
231{
232 MFEM_VERIFY(d1d <= q1d, "d1d > q1d is not supported");
233 MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
234 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
235 return internal::SmemPAConvectionNLGradApply2D<>;
236}
237
238template<int T_D1D, int T_Q1D>
239VectorConvectionNLFIntegrator::AddMultGradPAType
240VectorConvectionNLFIntegrator::AddMultGradPA3D::Kernel()
241{
242 static_assert(T_D1D <= T_Q1D, "d1d > q1d is not supported");
243 return internal::SmemPAConvectionNLGradApply3D<T_D1D, T_Q1D>;
244}
245
246inline VectorConvectionNLFIntegrator::AddMultGradPAType
247VectorConvectionNLFIntegrator::AddMultGradPA3D::Fallback(int d1d, int q1d)
248{
249 MFEM_VERIFY(d1d <= q1d, "d1d > q1d is not supported");
250 MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
251 MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
252 return internal::SmemPAConvectionNLGradApply3D<>;
253}
254
255/// \endcond DO_NOT_DOCUMENT
256
257} // namespace mfem
mfem::VectorConvectionNLFIntegrator::AddMultGradPAType
void(*)(const int ne, const real_t *B, const real_t *G, const real_t *A, const real_t *u, const real_t *x, real_t *y, const int d1d, const int q1d) AddMultGradPAType
Definition nonlininteg.hpp:429
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::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::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