Theseus
Compressible flow solver
Loading...
Searching...
No Matches
DGSEMIntegrator.hpp
Go to the documentation of this file.
1// Copyright (c) 2025-2026 Board of Trustees of the University of Illinois
2//
3// This file is part of Theseus.
4//
5// SPDX-License-Identifier: BSD-3-Clause
6#pragma once
7
8#include "mfem.hpp"
9#include "AxisymmetricSource.hpp"
10#include "dgsem_cache_utilities.hpp"
11#include "bc_kernels.hpp"
12
13namespace Theseus
14{
15
16 namespace DGSEMIntegrator
17 {
18
19 template<typename ContextType>
20 MFEM_HOST_DEVICE inline
21 static mfem::real_t AssembleElementVolumeKernel(const ContextType &ctx,
22 const mfem::real_t *el_u, const mfem::real_t *elJac_d,
23 const mfem::real_t *elMetric_d, mfem::real_t *el_dudt)
24 {
25
26 const int Np_x = ctx.Np_x;
27 const int Np_y = ctx.Np_y;
28 const int Np_z = ctx.Np_z;
29 const int dim = ctx.dim;
30 const int neq = ctx.num_equations;
31 const int dof = Np_x * Np_y * Np_z;
32 const mfem::real_t *Dhat2_d = ctx.Dhat2_d;
33
34 mfem::real_t f[Theseus::MAXEQ] = {0.,0.,0.,0.,0.};
35 mfem::real_t state1[Theseus::MAXEQ];
36 mfem::real_t state2[Theseus::MAXEQ];
37 mfem::real_t max_char_speed = 0.0;
38
39 { // X-direction (metric row 0)
40 // Zero'ing probably unnecessary: Chandrashekar flux overwrites it every time
41 // for(int q = 0;q < neq;q++) f[q] = 0.0;
42 for (int k = 0; k < Np_z; k++)
43 for (int j = 0; j < Np_y; j++)
44 for (int i = 0; i < Np_x; i++)
45 {
46 int id1 = k * Np_y * Np_x + j * Np_x + i;
47 Kernels::el_gather_state(el_u, dof, neq, id1, state1);
48 const mfem::real_t *met1 = elMetric_d+id1*dim*dim;
49 for (int m = i + 1; m < Np_x; m++)
50 {
51 int id2 = k * Np_y * Np_x + j * Np_x + m;
52 Kernels::el_gather_state(el_u, dof, neq, id2, state2);
53 const mfem::real_t *met2 = elMetric_d + id2*dim*dim;
54
55 const mfem::real_t cs = ctx.iflux.ComputeVolumeFlux(ctx.gas, state1, state2, met1, met2, f);
56 max_char_speed = Kernels::rmax(cs, max_char_speed);
57
58 const mfem::real_t c1 = Dhat2_d[m + Np_x*i];
59 const mfem::real_t c2 = Dhat2_d[i + Np_x*m];
60 Kernels::el_scatter_add(f, dof, neq, id1, c1, el_dudt);
61 Kernels::el_scatter_add(f, dof, neq, id2, c2, el_dudt);
62
63 }
64 }
65 } // X-direction block
66
67 // Y-direction (metric row 1)
68 if(dim > 1) {
69 for (int k = 0; k < Np_z; ++k)
70 for (int j = 0; j < Np_y; ++j)
71 for (int i = 0; i < Np_x; ++i)
72 {
73 const int id1 = k*Np_y*Np_x + j*Np_x + i;
74 Kernels::el_gather_state(el_u, dof, neq, id1, state1);
75 const mfem::real_t *met1 = elMetric_d + id1*dim*dim + 1*dim;
76
77 for (int m = j+1; m < Np_y; ++m)
78 {
79 const int id2 = k*Np_y*Np_x + m*Np_x + i;
80 Kernels::el_gather_state(el_u, dof, neq, id2, state2);
81 const mfem::real_t *met2 = elMetric_d + id2*dim*dim + dim;
82 // ComputeVolumeFlux *overwrites* f, so don't worry about reuse
83 const mfem::real_t cs = ctx.iflux.ComputeVolumeFlux(ctx.gas, state1, state2, met1, met2, f);
84 max_char_speed = Kernels::rmax(max_char_speed, cs);
85
86 const mfem::real_t c1 = Dhat2_d[m + Np_y*j]; // column j, entry m
87 const mfem::real_t c2 = Dhat2_d[j + Np_y*m]; // column m, entry j
88 Kernels::el_scatter_add(f, dof, neq, id1, c1, el_dudt);
89 Kernels::el_scatter_add(f, dof, neq, id2, c2, el_dudt);
90
91 }
92 }
93 } // Y-direction block
94
95 if (dim > 2) { // Z-direction (metric row 2)
96 for (int k = 0; k < Np_z; ++k)
97 for (int j = 0; j < Np_y; ++j)
98 for (int i = 0; i < Np_x; ++i)
99 {
100 const int id1 = k*Np_y*Np_x + j*Np_x + i;
101 Kernels::el_gather_state(el_u, dof, neq, id1, state1);
102 const mfem::real_t *met1 = elMetric_d + id1*dim*dim + 2*dim;
103
104 for (int m = k+1; m < Np_z; ++m)
105 {
106 const int id2 = m*Np_y*Np_x + j*Np_x + i;
107 Kernels::el_gather_state(el_u, dof, neq, id2, state2);
108 const mfem::real_t *met2 = elMetric_d + id2*dim*dim + 2*dim;
109
110 const mfem::real_t cs = ctx.iflux.ComputeVolumeFlux(ctx.gas, state1, state2, met1, met2, f);
111 max_char_speed = Kernels::rmax(max_char_speed, cs);
112
113 const mfem::real_t c1 = Dhat2_d[m + Np_z*k];
114 const mfem::real_t c2 = Dhat2_d[k + Np_z*m];
115 Kernels::el_scatter_add(f, dof, neq, id1, c1, el_dudt);
116 Kernels::el_scatter_add(f, dof, neq, id2, c2, el_dudt);
117 }
118 }
119 } // Z-direction block
120 // const int NPtot = Np_x * Np_y * Np_z; // = Np_x * Np_x * Np_x (!)
121 Kernels::el_scale(elJac_d, -1.0, dof, neq, el_dudt);
122
123 return max_char_speed;
124 }
125
126 template<typename ContextT>
127 MFEM_HOST_DEVICE static mfem::real_t AssembleElementFaceKernel(const ContextT &ctx, const mfem::real_t *u_face,
128 const mfem::real_t *nor_face,const mfem::real_t *w_minus,
129 const mfem::real_t *w_plus, mfem::real_t *rhs_face)
130 {
131 mfem::real_t max_char_speed = 0.0;
132 mfem::real_t point_flux[Theseus::MAXEQ];
133 mfem::real_t qMinus[Theseus::MAXEQ];
134 mfem::real_t qPlus[Theseus::MAXEQ];
135 const int nfp = ctx.num_face_points;
136 const int neq = ctx.num_equations;
137 const int dim = ctx.dim;
138 // auto idx = [=](int side, int fp, int eq) -> int
139 // {
140 // return (((side)*neq + eq)*nfp + fp);
141 // };
142 for (int i = 0; i < nfp; i++)
143 {
144 const mfem::real_t *nor_d = nor_face + i*dim;
145 const mfem::real_t wminus = -w_minus[i];
146 const mfem::real_t wplus = w_plus[i];
147 // Could avoid these copy-in,out
148 for(int j = 0;j < neq;j++){
149 qMinus[j] = u_face[ctx.iface_idx(0,i,j)];
150 qPlus[j] = u_face[ctx.iface_idx(1,i,j)];
151 }
152 max_char_speed = \
153 Kernels::rmax(max_char_speed, ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus,
154 nor_d, point_flux));
155 for(int j = 0;j < neq;j++){
156 rhs_face[ctx.iface_idx(0, i, j)] = wminus * point_flux[j];
157 rhs_face[ctx.iface_idx(1, i, j)] = wplus * point_flux[j];
158 }
159 }
160 return max_char_speed;
161 }
162
163 template<typename ContextT>
164 MFEM_HOST_DEVICE static mfem::real_t AssembleViscousElementFaceKernel(const ContextT &ctx, const mfem::real_t *u_face,
165 const mfem::real_t *nor_face,const mfem::real_t *w_minus,
166 const mfem::real_t *w_plus, const mfem::real_t *dprim_face_x,
167 const mfem::real_t *dprim_face_y, const mfem::real_t*dprim_face_z,
168 const mfem::real_t *face_radius,
169 mfem::real_t *rhs_face)
170 {
171 mfem::real_t max_char_speed = 0.0;
172 mfem::real_t point_flux[Theseus::MAXEQ];
173 mfem::real_t vflux_minus[Theseus::MAXEQ][Theseus::MAXDIM];
174 mfem::real_t vflux_plus[Theseus::MAXEQ][Theseus::MAXDIM];
175 mfem::real_t qMinus[Theseus::MAXEQ];
176 mfem::real_t qPlus[Theseus::MAXEQ];
177 mfem::real_t gradPrim_plus[Theseus::MAXDIM][Theseus::MAXEQ];
178 mfem::real_t gradPrim_minus[Theseus::MAXDIM][Theseus::MAXEQ];
179 const mfem::real_t *dprim_face[Theseus::MAXDIM] = {dprim_face_x, dprim_face_y, dprim_face_z};
180 const int nfp = ctx.num_face_points;
181 const int neq = ctx.num_equations;
182 const int dim = ctx.dim;
183 // auto idx = [=](int side, int fp, int eq) -> int
184 // {
185 // return (((side)*neq + eq)*nfp + fp);
186 // };
187 for (int i = 0; i < nfp; i++)
188 {
189 const mfem::real_t *nor_d = nor_face + i*dim;
190 const mfem::real_t wminus = -w_minus[i];
191 const mfem::real_t wplus = w_plus[i];
192 // Could avoid these copy-in,out
193 for(int j = 0;j < neq;j++){
194 int minus_index = ctx.iface_idx(0, i, j);
195 int plus_index = ctx.iface_idx(1, i, j);
196 qMinus[j] = u_face[minus_index];
197 qPlus[j] = u_face[plus_index];
198 for(int idim = 0;idim < dim;idim++){
199 gradPrim_minus[idim][j] = dprim_face[idim][minus_index];
200 gradPrim_plus[idim][j] = dprim_face[idim][plus_index];
201 }
202 }
203 max_char_speed = \
204 Kernels::rmax(max_char_speed, ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus,
205 nor_d, point_flux));
206
207 // Here, point_flux is +(F_inv * Normal)
208
209 // Grab the viscous flux
210 NavierStokesFlux::ComputeViscousFluxKernel(ctx.gas, qMinus,
211 gradPrim_minus[0],
212 gradPrim_minus[1],
213 gradPrim_minus[2], vflux_minus,
214 ctx.axisymmetric,
215 ctx.axisymmetric ? face_radius[i] : 0.0);
216 NavierStokesFlux::ComputeViscousFluxKernel(ctx.gas, qPlus,
217 gradPrim_plus[0],
218 gradPrim_plus[1],
219 gradPrim_plus[2], vflux_plus,
220 ctx.axisymmetric,
221 ctx.axisymmetric ? face_radius[i] : 0.0);
222
223 // Now we have vflux(+) and vflux(-)
224 // In this loop:
225 // - average vflux
226 // - dot avg vflux with nor
227 // - accumulate dotted (avg*n) into point_flux
228 for(int j = 0;j < neq;j++){
229 for(int idim = 0;idim < dim;idim++){
230 mfem::real_t avg = 0.5*(vflux_minus[j][idim] + vflux_plus[j][idim]);
231 point_flux[j] -= nor_d[idim]*avg;
232 }
233 }
234 // So now: point_flux = +(F_inv * Normal) -(F^bar_visc * Normal)
235 // in this loop:
236 // - SET/Overwrite rhs_face
237 // - NEGATE the (-) face point_flux to properly orient
238 for(int j = 0;j < neq;j++){
239 rhs_face[ctx.iface_idx(0, i, j)] = wminus * point_flux[j];
240 rhs_face[ctx.iface_idx(1, i, j)] = wplus * point_flux[j];
241 }
242 }
243 return max_char_speed;
244 }
245
246 template<typename ContextT>
247 MFEM_HOST_DEVICE inline static mfem::real_t ComputeFVFluxesKernel(const ContextT &ctx,
248 const mfem::real_t *el_u,
249 const mfem::real_t *elJac,
250 const mfem::real_t *el_metric_xi,
251 const mfem::real_t *el_metric_eta,
252 const mfem::real_t *el_metric_zeta,
253 mfem::real_t *el_dudt)
254 {
255 const int dim = ctx.dim;
256 const int Np_x = ctx.Np_x;
257 const int Np_y = ctx.Np_y;
258 const int Np_z = ctx.Np_z;
259 const int neq = ctx.num_equations;
260 const int npe = Np_x * Np_y * Np_z;
261 const mfem::real_t *qWgt = ctx.subcell_weights_d;
262
263 mfem::real_t max_char_speed = 0.0;
264 mfem::real_t flux_num[Theseus::MAXEQ];
265 mfem::real_t du_subcell[Theseus::MAXEQ];
266 mfem::real_t state1_local[Theseus::MAXEQ];
267 mfem::real_t state2_local[Theseus::MAXEQ];
268
269 for(int i = 0;i < npe*neq;i++)
270 el_dudt[i] = 0.0;
271
272 for (int k = 0; k < Np_z; k++)
273 {
274 for (int j = 0; j < Np_y; j++)
275 {
276 for(int q = 0; q < neq;q++){
277 du_subcell[q] = 0.0;
278 }
279 int id1 = k * Np_y * Np_x + j * Np_x;
280 Kernels::el_gather_state(el_u, npe, neq, id1, state1_local);
281 for (int i = 0; i < Np_x - 1; i++)
282 {
283 int id2 = id1 + 1;
284 Kernels::el_gather_state(el_u, npe, neq, id2, state2_local);
285 const mfem::real_t *nor = el_metric_xi + id2*dim;
286
287 max_char_speed = \
288 Kernels::rmax(max_char_speed,
289 ctx.iflux.ComputeFaceFlux(ctx.gas, state1_local,
290 state2_local, nor, flux_num));
291 for(int q = 0; q < neq;q++){
292 du_subcell[q] -= flux_num[q];
293 }
294 for(int q = 0; q < neq;q++){
295 du_subcell[q] /= (elJac[id1] * qWgt[i]);
296 }
297 Kernels::el_scatter_assign(du_subcell, npe, neq, id1, 1.0, el_dudt);
298 for(int q = 0; q < neq;q++){
299 du_subcell[q] = flux_num[q];
300 }
301 for(int q = 0;q < neq;q++){
302 state1_local[q] = state2_local[q];
303 }
304 id1 = id2;
305 }
306 for(int q = 0;q < neq;q++){
307 du_subcell[q] /= (elJac[id1] * qWgt[Np_x-1]);
308 }
309 Kernels::el_scatter_assign(du_subcell, npe, neq, id1, 1.0, el_dudt);
310 }
311 }
312
313 if (dim > 1)
314 {
315 for (int k = 0; k < Np_z; k++)
316 {
317 for (int i = 0; i < Np_x; i++)
318 {
319 for(int q = 0; q < neq;q++){
320 du_subcell[q] = 0.0;
321 }
322 int id1 = k * Np_y * Np_x + i;
323 Kernels::el_gather_state(el_u, npe, neq, id1,
324 state1_local);
325 for (int j = 0; j < Np_y - 1; j++)
326 {
327 int id2 = k * Np_y * Np_x + (j + 1) * Np_x + i;
328 Kernels::el_gather_state(el_u, npe, neq, id2,
329 state2_local);
330 const mfem::real_t *nor = el_metric_eta + id2*dim;
331 max_char_speed = \
332 Kernels::rmax(max_char_speed,
333 ctx.iflux.ComputeFaceFlux(ctx.gas,
334 state1_local,
335 state2_local,
336 nor, flux_num));
337 for(int q = 0;q < neq;q++){
338 du_subcell[q] -= flux_num[q];
339 }
340 for(int q = 0;q < neq;q++){
341 du_subcell[q] /= (elJac[id1] * qWgt[j]);
342 }
343 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
344 for(int q = 0;q < neq;q++){
345 du_subcell[q] = flux_num[q];
346 state1_local[q] = state2_local[q];
347 }
348 id1 = id2;
349 }
350 for(int q = 0;q < neq;q++){
351 du_subcell[q] /= (elJac[id1] * qWgt[Np_y - 1]);
352 }
353 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
354 }
355 }
356 if (dim > 2)
357 {
358 for (int j = 0; j < Np_y; j++)
359 {
360 for (int i = 0; i < Np_x; i++)
361 {
362 for(int q = 0; q < neq;q++){
363 du_subcell[q] = 0.0;
364 }
365 int id1 = j * Np_x + i;
366 Kernels::el_gather_state(el_u, npe, neq, id1,
367 state1_local);
368 for (int k = 0; k < Np_z - 1; k++)
369 {
370 int id2 = (k + 1) * Np_y * Np_x + j * Np_x + i;
371 Kernels::el_gather_state(el_u, npe, neq, id2,
372 state2_local);
373 const mfem::real_t *nor = el_metric_zeta + id2*dim;
374 max_char_speed = \
375 Kernels::rmax(max_char_speed,
376 ctx.iflux.ComputeFaceFlux(ctx.gas, state1_local,
377 state2_local, nor, flux_num));
378 for(int q = 0;q < neq;q++){
379 du_subcell[q] -= flux_num[q];
380 }
381 for(int q = 0;q < neq;q++){
382 du_subcell[q] /= (elJac[id1] * qWgt[k]);
383 }
384 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
385
386 for(int q = 0;q < neq;q++){
387 du_subcell[q] = flux_num[q];
388 state1_local[q] = state2_local[q];
389 }
390 id1 = id2;
391 }
392 for(int q = 0;q < neq;q++){
393 du_subcell[q] /= (elJac[id1] * qWgt[Np_z - 1]);
394 }
395 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
396 }
397 }
398 }
399 }
400 return max_char_speed;
401 }
402
403
404 template<typename ContextType>
405 MFEM_HOST_DEVICE inline
406 static void AssembleViscousElementVolumeKernel(const ContextType &ctx,
407 const mfem::real_t *el_u,
408 const mfem::real_t *elJac_d,
409 const mfem::real_t *elMetric_d,
410 const mfem::real_t *elRadius_d,
411 const mfem::real_t *el_gradprim_x,
412 const mfem::real_t *el_gradprim_y,
413 const mfem::real_t *el_gradprim_z,
414 mfem::real_t *el_dudt)
415 {
416 const int Np_x = ctx.Np_x;
417 const int Np_y = ctx.Np_y;
418 const int Np_z = ctx.Np_z;
419 const int dim = ctx.dim;
420 const int neq = ctx.num_equations;
421 const int dof = Np_x * Np_y * Np_z;
422 const mfem::real_t *Dhat_d = ctx.Dhat_d;
423
424 // One source-point scratch
425 mfem::real_t state[Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
426 mfem::real_t dqx [Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
427 mfem::real_t dqy [Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
428 mfem::real_t dqz [Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
429
430 // flux(eq,dir)
431 // one transformed reference-direction flux vector
432 mfem::real_t f_ref[Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
433
434 for (int k = 0; k < Np_z; ++k)
435 {
436 for (int j = 0; j < Np_y; ++j)
437 {
438 for (int i = 0; i < Np_x; ++i)
439 {
440 const int id1 = k * Np_y * Np_x + j * Np_x + i;
441 const mfem::real_t J = elJac_d[id1];
442 const mfem::real_t jInv = 1.0/J;
443
444 mfem::real_t dU_viscous[Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
445
446 // xi contribution
447 for (int l = 0; l < Np_x; ++l)
448 {
449 const int idl = k * Np_y * Np_x + j * Np_x + l;
450 const mfem::real_t c = Dhat_d[l + Np_x * i];
451
452 Kernels::el_gather_state(el_u, dof, neq, idl, state);
453 Kernels::el_gather_grad_state(el_gradprim_x, el_gradprim_y,
454 el_gradprim_z, dim, dof, neq, idl,
455 dqx, dqy, dqz);
456
457 const mfem::real_t *adj_row = elMetric_d + idl * dim * dim + 0 * dim;
458 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
459 ctx.gas, dim, neq, state, dqx, dqy, dqz, adj_row,
460 f_ref, ctx.axisymmetric,
461 ctx.axisymmetric ? elRadius_d[idl] : 0.0);
462 for (int q = 0; q < neq; ++q)
463 {
464 dU_viscous[q] += c * f_ref[q];
465 }
466 }
467
468 // eta contribution
469 if (dim > 1)
470 {
471 for (int l = 0; l < Np_y; ++l)
472 {
473 const int idl = k * Np_y * Np_x + l * Np_x + i;
474 const mfem::real_t c = Dhat_d[l + Np_y * j];
475
476 Kernels::el_gather_state(el_u, dof, neq, idl, state);
477 Kernels::el_gather_grad_state(el_gradprim_x, el_gradprim_y,
478 el_gradprim_z, dim, dof, neq, idl,
479 dqx, dqy, dqz);
480
481 const mfem::real_t *adj_row = elMetric_d + idl * dim * dim + 1 * dim;
482 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
483 ctx.gas, dim, neq, state, dqx, dqy, dqz, adj_row,
484 f_ref, ctx.axisymmetric,
485 ctx.axisymmetric ? elRadius_d[idl] : 0.0);
486
487 for (int q = 0; q < neq; ++q)
488 {
489 dU_viscous[q] += c * f_ref[q];
490 }
491 }
492 }
493
494 // zeta contribution
495 if (dim > 2)
496 {
497 for (int l = 0; l < Np_z; ++l)
498 {
499 const int idl = l * Np_y * Np_x + j * Np_x + i;
500 const mfem::real_t c = Dhat_d[l + Np_z * k];
501
502 Kernels::el_gather_state(el_u, dof, neq, idl, state);
503 Kernels::el_gather_grad_state(el_gradprim_x, el_gradprim_y,
504 el_gradprim_z, dim, dof, neq, idl,
505 dqx, dqy, dqz);
506
507 const mfem::real_t *adj_row = elMetric_d + idl * dim * dim + 2 * dim;
508 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
509 ctx.gas, dim, neq, state, dqx, dqy, dqz, adj_row,
510 f_ref, ctx.axisymmetric,
511 ctx.axisymmetric ? elRadius_d[idl] : 0.0);
512
513 for (int q = 0; q < neq; ++q)
514 {
515 dU_viscous[q] += c * f_ref[q];
516 }
517 }
518 }
519 Kernels::el_scatter_add(dU_viscous, dof, neq, id1, jInv, el_dudt);
520 if (ctx.axisymmetric)
521 {
522 Kernels::el_gather_state(el_u, dof, neq, id1, state);
523 Kernels::el_gather_grad_state(
524 el_gradprim_x, el_gradprim_y, el_gradprim_z, dim,
525 dof, neq, id1, dqx, dqy, dqz);
526 mfem::real_t source[Theseus::MAXEQ] = {0.0};
527 if (!AddAxisymmetricViscousSourceAwayFromAxis(
528 ctx.gas, state, dqx, dqy, dqz,
529 elRadius_d[id1], source))
530 {
531 AddAxisymmetricViscousSourceAtAxis(
532 ctx, el_u, el_gradprim_x, el_gradprim_y,
533 el_gradprim_z, elRadius_d, elJac_d, elMetric_d,
534 id1, source);
535 }
536 Kernels::el_scatter_add(source, dof, neq, id1, 1.0,
537 el_dudt);
538 }
539 }
540 }
541 }
542 }
543
544 template <typename ContextType>
545 MFEM_HOST_DEVICE inline
546 static void AssembleGradElementVolumeKernel(const ContextType &ctx,
547 const mfem::real_t *el_u,
548 const mfem::real_t *elJac_d,
549 const mfem::real_t *elMetric_d,
550 mfem::real_t *el_grad_u[Theseus::MAXDIM])
551 {
552 const int Np_x = ctx.Np_x;
553 const int Np_y = ctx.Np_y;
554 const int Np_z = ctx.Np_z;
555 const int neq = ctx.num_equations;
556 const int dim = ctx.dim;
557 const int dof = Np_x * Np_y * Np_z;
558 const mfem::real_t *D_d = ctx.D_d;
559
560 if(dim == 1){
561
562 // Keep MAX_EQ in mind later if neq can exceed 5.
563 mfem::real_t dudxi[Theseus::MAXEQ];
564
565 for (int i = 0; i < Np_x; ++i)
566 {
567 const int id = i;
568
569 for (int q = 0; q < neq; ++q)
570 {
571 dudxi[q] = 0.0;
572 }
573 // Reference-space derivatives.
574 for (int l = 0; l < Np_x; ++l)
575 {
576 const int id_x = l;
577 const mfem::real_t c_xi = D_d[i*Np_x + l]; // legacy D_T(l, i)
578
579 for (int q = 0; q < neq; ++q)
580 {
581 dudxi[q] += el_u[id_x + q * dof] * c_xi;
582 }
583 }
584
585 const mfem::real_t invJ = 1.0 / elJac_d[id];
586 const mfem::real_t *adj = elMetric_d + id * dim * dim;
587
588 for (int q = 0; q < neq; ++q)
589 {
590 el_grad_u[0][id + q * dof] = invJ * (dudxi[q] * adj[0]);
591 }
592 }
593 } else if(dim == 2){
594 mfem::real_t dudxi[Theseus::MAXEQ];
595 mfem::real_t dudeta[Theseus::MAXEQ];
596
597 for (int j = 0; j < Np_y; ++j)
598 {
599 for (int i = 0; i < Np_x; ++i)
600 {
601 const int id = j * Np_x + i;
602
603 for (int q = 0; q < neq; ++q)
604 {
605 dudxi[q] = 0.0;
606 dudeta[q] = 0.0;
607 }
608
609 // Reference-space derivatives.
610 for (int l = 0; l < Np_x; ++l)
611 {
612 const int id_x = j * Np_x + l;
613 const int id_y = l * Np_x + i;
614
615 const mfem::real_t c_xi = D_d[l + Np_x * i]; // legacy D_T(l,i)
616 const mfem::real_t c_eta = D_d[l + Np_x * j]; // legacy D_T(l,j)
617
618 for (int q = 0; q < neq; ++q)
619 {
620 dudxi[q] += el_u[id_x + q * dof] * c_xi;
621 dudeta[q] += el_u[id_y + q * dof] * c_eta;
622 }
623 }
624
625 const mfem::real_t invJ = 1.0 / elJac_d[id];
626 const mfem::real_t *adj = elMetric_d + id * dim * dim;
627
628 // adj stored row-major per point:
629 // [ adj[0] adj[1] ]
630 // [ adj[2] adj[3] ]
631 for (int q = 0; q < neq; ++q)
632 {
633 el_grad_u[0][id + q * dof] = invJ * (dudxi[q] * adj[0] + dudeta[q] * adj[2]);
634 el_grad_u[1][id + q * dof] = invJ * (dudxi[q] * adj[1] + dudeta[q] * adj[3]);
635 }
636 }
637 }
638 } else if (dim == 3) {
639
640 mfem::real_t dudxi[Theseus::MAXEQ];
641 mfem::real_t dudeta[Theseus::MAXEQ];
642 mfem::real_t dudzeta[Theseus::MAXEQ];
643
644 for (int k = 0; k < Np_z; ++k)
645 {
646 for (int j = 0; j < Np_y; ++j)
647 {
648 for (int i = 0; i < Np_x; ++i)
649 {
650 const int id = k * Np_x * Np_y + j * Np_x + i;
651
652 for (int q = 0; q < neq; ++q)
653 {
654 dudxi[q] = 0.0;
655 dudeta[q] = 0.0;
656 dudzeta[q] = 0.0;
657 }
658
659 // Reference-space derivatives.
660 for (int l = 0; l < Np_x; ++l)
661 {
662 const int id_x = k * Np_x * Np_y + j * Np_x + l;
663 const int id_y = k * Np_x * Np_y + l * Np_x + i;
664 const int id_z = l * Np_x * Np_y + j * Np_x + i;
665 const mfem::real_t c_xi = D_d[l + Np_x * i]; // legacy D_T(l,i)
666 const mfem::real_t c_eta = D_d[l + Np_x * j]; // legacy D_T(l,j)
667 const mfem::real_t c_zeta = D_d[l + Np_x * k]; // legacy D_T(l,k)
668 for (int q = 0; q < neq; ++q)
669 {
670 dudxi[q] += el_u[id_x + q * dof] * c_xi;
671 dudeta[q] += el_u[id_y + q * dof] * c_eta;
672 dudzeta[q] += el_u[id_z + q * dof] * c_zeta;
673 }
674 }
675
676 const mfem::real_t invJ = 1.0 / elJac_d[id];
677 const mfem::real_t *adj = elMetric_d + id * dim * dim;
678
679 // adj stored row-major per point:
680 // [ adj[0] adj[1] adj[2] ]
681 // [ adj[3] adj[4] adj[5] ]
682 // [ adj[6] adj[7] adj[8] ]
683 for (int q = 0; q < neq; ++q)
684 {
685 el_grad_u[0][id + q * dof] = invJ * (dudxi[q] * adj[0] +
686 dudeta[q] * adj[3] +
687 dudzeta[q] * adj[6]);
688
689 el_grad_u[1][id + q * dof] = invJ * (dudxi[q] * adj[1] +
690 dudeta[q] * adj[4] +
691 dudzeta[q] * adj[7]);
692
693 el_grad_u[2][id + q * dof] = invJ * (dudxi[q] * adj[2] +
694 dudeta[q] * adj[5] +
695 dudzeta[q] * adj[8]);
696 }
697 }
698 }
699 }
700 }
701 }
702
703 template <typename ContextT>
704 MFEM_HOST_DEVICE inline
705 static void AssembleGradInteriorFaceKernel(const ContextT &ctx,
706 const mfem::real_t *u_face,
707 const mfem::real_t *nor_face,
708 const mfem::real_t *w_minus,
709 const mfem::real_t *w_plus,
710 mfem::real_t *rhs_face[Theseus::MAXDIM])
711 {
712 const int nfp = ctx.num_face_points;
713 const int neq = ctx.num_equations;
714 const int dim = ctx.dim;
715
716 mfem::real_t qMinus[Theseus::MAXEQ];
717 mfem::real_t qPlus[Theseus::MAXEQ];
718 mfem::real_t jump[Theseus::MAXEQ];
719
720 for (int i = 0; i < nfp; ++i)
721 {
722 const mfem::real_t *nor_d = nor_face + i * dim;
723
724 const mfem::real_t wminus = w_minus[i];
725 const mfem::real_t wplus = w_plus[i];
726
727 for (int q = 0; q < neq; ++q)
728 {
729 qMinus[q] = u_face[ctx.iface_idx(0, i, q)];
730 qPlus[q] = u_face[ctx.iface_idx(1, i, q)];
731 jump[q] = mfem::real_t(0.5) * (qPlus[q] - qMinus[q]);
732 }
733
734 for ( int idim = 0;idim < dim;idim++){
735 mfem::real_t *rhs_d = rhs_face[idim];
736 const mfem::real_t n_d = nor_d[idim];
737 for (int q = 0; q < neq; ++q)
738 {
739 const mfem::real_t f_d = jump[q]*n_d;
740 rhs_d[ctx.iface_idx(0, i, q)] = wminus * f_d;
741 rhs_d[ctx.iface_idx(1, i, q)] = wplus * f_d;
742 }
743 }
744 }
745 }
746
747 template <typename DeviceCacheT>
748 MFEM_HOST_DEVICE inline
749 static void AssembleGradBoundaryPointKernel(const DeviceCacheT &dc,
750 const Theseus::BCDescriptor &bc,
751 const mfem::real_t *u_face,
752 const mfem::real_t *nor_point,
753 const mfem::real_t scale,
754 const int fp,
755 mfem::real_t *rhs_face[Theseus::MAXDIM])
756 {
757 const int dim = dc.dim;
758 const int nfp = dc.num_face_points;
759 const int neq = dc.num_equations;
760
761 mfem::real_t state1[Theseus::MAXEQ];
762 mfem::real_t fluxN[Theseus::MAXEQ];
763 mfem::real_t flux_dir[Theseus::MAXEQ];
764
765 Theseus::Kernels::el_gather_state(u_face, nfp, neq, fp, state1);
766
767 Theseus::BC::ComputeBdrFaceGradFlux(dc, bc, state1, fluxN);
768
769 for(int idim = 0;idim < dim;idim++){
770 for(int q = 0;q < neq;q++){
771 flux_dir[q] = fluxN[q]*nor_point[idim];
772 }
773 Theseus::Kernels::el_scatter_add(flux_dir, nfp, neq, fp, scale, rhs_face[idim]);
774 }
775 }
776
777 };
778}
bc_kernels.hpp
dgsem_cache_utilities.hpp
Theseus::BC::ComputeBdrFaceGradFlux
MFEM_HOST_DEVICE void ComputeBdrFaceGradFlux(const DeviceCacheT &dc, const Theseus::BCDescriptor &bc, const mfem::real_t *state1, mfem::real_t *fluxN)
Definition bc_kernels.hpp:34
Theseus::Kernels::el_scale
MFEM_HOST_DEVICE void el_scale(const mfem::real_t *scale_d, const mfem::real_t fac, const int dof, const int neq, mfem::real_t *el_soln)
Definition theseus_kernels.hpp:187
Theseus::Kernels::el_scatter_add
MFEM_HOST_DEVICE void el_scatter_add(const mfem::real_t *f, const int dof, const int num_eq, const int id, const mfem::real_t scale, mfem::real_t *du)
Definition theseus_kernels.hpp:153
Theseus::Kernels::el_gather_state
MFEM_HOST_DEVICE void el_gather_state(const mfem::real_t *u, const int dof, const int num_eq, const int id, mfem::real_t *dst)
Definition theseus_kernels.hpp:135
Theseus::Kernels::el_scatter_assign
MFEM_HOST_DEVICE void el_scatter_assign(const mfem::real_t *f, const int dof, const int num_eq, const int id, const mfem::real_t scale, mfem::real_t *du)
Definition theseus_kernels.hpp:170
Theseus::Kernels::el_gather_grad_state
MFEM_HOST_DEVICE void el_gather_grad_state(const mfem::real_t *grad_state_x, const mfem::real_t *grad_state_y, const mfem::real_t *grad_state_z, const int dim, const int dof, const int neq, const int id, mfem::real_t *dqx, mfem::real_t *dqy, mfem::real_t *dqz)
Definition theseus_kernels.hpp:143
Theseus::Kernels::rmax
MFEM_HOST_DEVICE mfem::real_t rmax(mfem::real_t a, mfem::real_t b)
Definition theseus_kernels.hpp:17
Theseus
Definition AxisymmetricGeometry.hpp:15
Theseus::AddAxisymmetricViscousSourceAwayFromAxis
MFEM_HOST_DEVICE bool AddAxisymmetricViscousSourceAwayFromAxis(const GasT &gas, const mfem::real_t *state, const mfem::real_t *dprim_x, const mfem::real_t *dprim_y, const mfem::real_t *dprim_z, mfem::real_t radius, mfem::real_t *state_rate)
Definition AxisymmetricSource.hpp:68
Theseus::MAXDIM
constexpr const int MAXDIM
Definition theseus_kernels.hpp:14
Theseus::MAXEQ
constexpr const int MAXEQ
Definition theseus_kernels.hpp:13
Theseus::AddAxisymmetricViscousSourceAtAxis
MFEM_HOST_DEVICE void AddAxisymmetricViscousSourceAtAxis(const ContextT &ctx, const mfem::real_t *element_state, const mfem::real_t *element_gradprim_x, const mfem::real_t *element_gradprim_y, const mfem::real_t *element_gradprim_z, const mfem::real_t *element_radius, const mfem::real_t *element_jacobian, const mfem::real_t *element_metric, int point, mfem::real_t *state_rate)
Definition AxisymmetricSource.hpp:161
Theseus::BCDescriptor
Definition bc_cache_utilities.hpp:35