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 void 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
38 { // X-direction (metric row 0)
39 // Zero'ing probably unnecessary: Chandrashekar flux overwrites it every time
40 // for(int q = 0;q < neq;q++) f[q] = 0.0;
41 for (int k = 0; k < Np_z; k++)
42 for (int j = 0; j < Np_y; j++)
43 for (int i = 0; i < Np_x; i++)
44 {
45 int id1 = k * Np_y * Np_x + j * Np_x + i;
46 Kernels::el_gather_state(el_u, dof, neq, id1, state1);
47 const mfem::real_t *met1 = elMetric_d+id1*dim*dim;
48 for (int m = i + 1; m < Np_x; m++)
49 {
50 int id2 = k * Np_y * Np_x + j * Np_x + m;
51 Kernels::el_gather_state(el_u, dof, neq, id2, state2);
52 const mfem::real_t *met2 = elMetric_d + id2*dim*dim;
53
54 ctx.iflux.ComputeVolumeFlux(ctx.gas, state1, state2,
55 met1, met2, f);
56
57 const mfem::real_t c1 = Dhat2_d[m + Np_x*i];
58 const mfem::real_t c2 = Dhat2_d[i + Np_x*m];
59 Kernels::el_scatter_add(f, dof, neq, id1, c1, el_dudt);
60 Kernels::el_scatter_add(f, dof, neq, id2, c2, el_dudt);
61
62 }
63 }
64 } // X-direction block
65
66 // Y-direction (metric row 1)
67 if(dim > 1) {
68 for (int k = 0; k < Np_z; ++k)
69 for (int j = 0; j < Np_y; ++j)
70 for (int i = 0; i < Np_x; ++i)
71 {
72 const int id1 = k*Np_y*Np_x + j*Np_x + i;
73 Kernels::el_gather_state(el_u, dof, neq, id1, state1);
74 const mfem::real_t *met1 = elMetric_d + id1*dim*dim + 1*dim;
75
76 for (int m = j+1; m < Np_y; ++m)
77 {
78 const int id2 = k*Np_y*Np_x + m*Np_x + i;
79 Kernels::el_gather_state(el_u, dof, neq, id2, state2);
80 const mfem::real_t *met2 = elMetric_d + id2*dim*dim + dim;
81 // ComputeVolumeFlux *overwrites* f, so don't worry about reuse
82 ctx.iflux.ComputeVolumeFlux(ctx.gas, state1, state2,
83 met1, met2, f);
84
85 const mfem::real_t c1 = Dhat2_d[m + Np_y*j]; // column j, entry m
86 const mfem::real_t c2 = Dhat2_d[j + Np_y*m]; // column m, entry j
87 Kernels::el_scatter_add(f, dof, neq, id1, c1, el_dudt);
88 Kernels::el_scatter_add(f, dof, neq, id2, c2, el_dudt);
89
90 }
91 }
92 } // Y-direction block
93
94 if (dim > 2) { // Z-direction (metric row 2)
95 for (int k = 0; k < Np_z; ++k)
96 for (int j = 0; j < Np_y; ++j)
97 for (int i = 0; i < Np_x; ++i)
98 {
99 const int id1 = k*Np_y*Np_x + j*Np_x + i;
100 Kernels::el_gather_state(el_u, dof, neq, id1, state1);
101 const mfem::real_t *met1 = elMetric_d + id1*dim*dim + 2*dim;
102
103 for (int m = k+1; m < Np_z; ++m)
104 {
105 const int id2 = m*Np_y*Np_x + j*Np_x + i;
106 Kernels::el_gather_state(el_u, dof, neq, id2, state2);
107 const mfem::real_t *met2 = elMetric_d + id2*dim*dim + 2*dim;
108
109 ctx.iflux.ComputeVolumeFlux(ctx.gas, state1, state2,
110 met1, met2, f);
111
112 const mfem::real_t c1 = Dhat2_d[m + Np_z*k];
113 const mfem::real_t c2 = Dhat2_d[k + Np_z*m];
114 Kernels::el_scatter_add(f, dof, neq, id1, c1, el_dudt);
115 Kernels::el_scatter_add(f, dof, neq, id2, c2, el_dudt);
116 }
117 }
118 } // Z-direction block
119 // const int NPtot = Np_x * Np_y * Np_z; // = Np_x * Np_x * Np_x (!)
120 Kernels::el_scale(elJac_d, -1.0, dof, neq, el_dudt);
121
122 }
123
124 template<typename ContextType>
125 MFEM_HOST_DEVICE inline
126 static void AssembleVolumePointKernel(
127 const ContextType &ctx, const mfem::real_t *el_u,
128 const mfem::real_t *elJac_d, const mfem::real_t *elMetric_d,
129 const int point, mfem::real_t *el_dudt)
130 {
131 const int Np_x = ctx.Np_x;
132 const int Np_y = ctx.Np_y;
133 const int Np_z = ctx.Np_z;
134 const int dim = ctx.dim;
135 const int neq = ctx.num_equations;
136 const int dof = ctx.ndof_scalar_el;
137 const mfem::real_t *Dhat2_d = ctx.Dhat2_d;
138 const int i = point % Np_x;
139 const int j = (point / Np_x) % Np_y;
140 const int k = point / (Np_x*Np_y);
141
142 mfem::real_t state_lower[Theseus::MAXEQ];
143 mfem::real_t state_upper[Theseus::MAXEQ];
144 mfem::real_t flux[Theseus::MAXEQ];
145 mfem::real_t point_rate[Theseus::MAXEQ] = {0.0};
146
147 for (int m = 0; m < Np_x; ++m)
148 {
149 if (m == i) { continue; }
150 const int lower = m < i ? m : i;
151 const int upper = m < i ? i : m;
152 const int lower_point = k*Np_y*Np_x + j*Np_x + lower;
153 const int upper_point = k*Np_y*Np_x + j*Np_x + upper;
154 Kernels::el_gather_state(el_u, dof, neq, lower_point, state_lower);
155 Kernels::el_gather_state(el_u, dof, neq, upper_point, state_upper);
156 ctx.iflux.ComputeVolumeFlux(ctx.gas, state_lower, state_upper,
157 elMetric_d + lower_point*dim*dim,
158 elMetric_d + upper_point*dim*dim, flux);
159 const mfem::real_t coefficient = Dhat2_d[m + Np_x*i];
160 for (int q = 0; q < neq; ++q)
161 {
162 point_rate[q] += coefficient*flux[q];
163 }
164 }
165
166 if (dim > 1)
167 {
168 for (int m = 0; m < Np_y; ++m)
169 {
170 if (m == j) { continue; }
171 const int lower = m < j ? m : j;
172 const int upper = m < j ? j : m;
173 const int lower_point = k*Np_y*Np_x + lower*Np_x + i;
174 const int upper_point = k*Np_y*Np_x + upper*Np_x + i;
175 Kernels::el_gather_state(
176 el_u, dof, neq, lower_point, state_lower);
177 Kernels::el_gather_state(
178 el_u, dof, neq, upper_point, state_upper);
179 ctx.iflux.ComputeVolumeFlux(
180 ctx.gas, state_lower, state_upper,
181 elMetric_d + lower_point*dim*dim + dim,
182 elMetric_d + upper_point*dim*dim + dim, flux);
183 const mfem::real_t coefficient = Dhat2_d[m + Np_y*j];
184 for (int q = 0; q < neq; ++q)
185 {
186 point_rate[q] += coefficient*flux[q];
187 }
188 }
189 }
190
191 if (dim > 2)
192 {
193 for (int m = 0; m < Np_z; ++m)
194 {
195 if (m == k) { continue; }
196 const int lower = m < k ? m : k;
197 const int upper = m < k ? k : m;
198 const int lower_point = lower*Np_y*Np_x + j*Np_x + i;
199 const int upper_point = upper*Np_y*Np_x + j*Np_x + i;
200 Kernels::el_gather_state(
201 el_u, dof, neq, lower_point, state_lower);
202 Kernels::el_gather_state(
203 el_u, dof, neq, upper_point, state_upper);
204 ctx.iflux.ComputeVolumeFlux(
205 ctx.gas, state_lower, state_upper,
206 elMetric_d + lower_point*dim*dim + 2*dim,
207 elMetric_d + upper_point*dim*dim + 2*dim, flux);
208 const mfem::real_t coefficient = Dhat2_d[m + Np_z*k];
209 for (int q = 0; q < neq; ++q)
210 {
211 point_rate[q] += coefficient*flux[q];
212 }
213 }
214 }
215
216 Kernels::el_scatter_assign(
217 point_rate, dof, neq, point, -1.0/elJac_d[point], el_dudt);
218 }
219
220 template<typename ContextT>
221 MFEM_HOST_DEVICE static void AssembleFacePointKernel(const ContextT &ctx,
222 const mfem::real_t *u_face,
223 const mfem::real_t *nor_point,
224 const mfem::real_t w_minus,
225 const mfem::real_t w_plus,
226 const int fp,
227 mfem::real_t *rhs_face)
228 {
229 mfem::real_t point_flux[Theseus::MAXEQ];
230 mfem::real_t qMinus[Theseus::MAXEQ];
231 mfem::real_t qPlus[Theseus::MAXEQ];
232 const int neq = ctx.num_equations;
233 for(int q = 0; q < neq; ++q){
234 qMinus[q] = u_face[ctx.iface_idx(0, fp, q)];
235 qPlus[q] = u_face[ctx.iface_idx(1, fp, q)];
236 }
237
238 ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus, nor_point, point_flux);
239
240 for(int q = 0; q < neq; ++q){
241 rhs_face[ctx.iface_idx(0, fp, q)] = -w_minus * point_flux[q];
242 rhs_face[ctx.iface_idx(1, fp, q)] = w_plus * point_flux[q];
243 }
244
245 }
246
247 template<typename ContextT>
248 MFEM_HOST_DEVICE static void AssembleElementFaceKernel(const ContextT &ctx, const mfem::real_t *u_face,
249 const mfem::real_t *nor_face,const mfem::real_t *w_minus,
250 const mfem::real_t *w_plus, mfem::real_t *rhs_face)
251 {
252 const int nfp = ctx.num_face_points;
253 const int dim = ctx.dim;
254 for (int fp = 0; fp < nfp; ++fp)
255 {
256 AssembleFacePointKernel(ctx, u_face, nor_face + fp*dim,
257 w_minus[fp], w_plus[fp], fp, rhs_face);
258 }
259 }
260
261 template<typename ContextT>
262 MFEM_HOST_DEVICE static void AssembleViscousFacePointKernel(
263 const ContextT &ctx, const mfem::real_t *u_face,
264 const mfem::real_t *nor_point, const mfem::real_t w_minus,
265 const mfem::real_t w_plus, const mfem::real_t *dprim_face_x,
266 const mfem::real_t *dprim_face_y, const mfem::real_t *dprim_face_z,
267 const mfem::real_t radius, const int fp, mfem::real_t *rhs_face)
268 {
269 mfem::real_t point_flux[Theseus::MAXEQ];
270 mfem::real_t vflux_minus[Theseus::MAXEQ][Theseus::MAXDIM];
271 mfem::real_t vflux_plus[Theseus::MAXEQ][Theseus::MAXDIM];
272 mfem::real_t qMinus[Theseus::MAXEQ];
273 mfem::real_t qPlus[Theseus::MAXEQ];
274 mfem::real_t gradPrim_plus[Theseus::MAXDIM][Theseus::MAXEQ];
275 mfem::real_t gradPrim_minus[Theseus::MAXDIM][Theseus::MAXEQ];
276 const mfem::real_t *dprim_face[Theseus::MAXDIM] = {
277 dprim_face_x, dprim_face_y, dprim_face_z};
278 const int neq = ctx.num_equations;
279 const int dim = ctx.dim;
280
281 for(int q = 0; q < neq; ++q){
282 const int minus_index = ctx.iface_idx(0, fp, q);
283 const int plus_index = ctx.iface_idx(1, fp, q);
284 qMinus[q] = u_face[minus_index];
285 qPlus[q] = u_face[plus_index];
286 for(int idim = 0; idim < dim; ++idim){
287 gradPrim_minus[idim][q] = dprim_face[idim][minus_index];
288 gradPrim_plus[idim][q] = dprim_face[idim][plus_index];
289 }
290 }
291
292 ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus, nor_point, point_flux);
293
294 NavierStokesFlux::ComputeViscousFluxKernel(
295 ctx.gas, qMinus, gradPrim_minus[0], gradPrim_minus[1],
296 gradPrim_minus[2], vflux_minus, ctx.axisymmetric,
297 ctx.axisymmetric ? radius : 0.0);
298 NavierStokesFlux::ComputeViscousFluxKernel(
299 ctx.gas, qPlus, gradPrim_plus[0], gradPrim_plus[1],
300 gradPrim_plus[2], vflux_plus, ctx.axisymmetric,
301 ctx.axisymmetric ? radius : 0.0);
302
303 for(int q = 0; q < neq; ++q){
304 for(int idim = 0; idim < dim; ++idim){
305 const mfem::real_t avg =
306 0.5*(vflux_minus[q][idim] + vflux_plus[q][idim]);
307 point_flux[q] -= nor_point[idim]*avg;
308 }
309 }
310
311 for(int q = 0; q < neq; ++q){
312 rhs_face[ctx.iface_idx(0, fp, q)] = -w_minus * point_flux[q];
313 rhs_face[ctx.iface_idx(1, fp, q)] = w_plus * point_flux[q];
314 }
315
316 }
317
318 template<typename ContextT>
319 MFEM_HOST_DEVICE static void AssembleViscousElementFaceKernel(const ContextT &ctx, const mfem::real_t *u_face,
320 const mfem::real_t *nor_face,const mfem::real_t *w_minus,
321 const mfem::real_t *w_plus, const mfem::real_t *dprim_face_x,
322 const mfem::real_t *dprim_face_y, const mfem::real_t*dprim_face_z,
323 const mfem::real_t *face_radius,
324 mfem::real_t *rhs_face)
325 {
326 const int nfp = ctx.num_face_points;
327 const int dim = ctx.dim;
328 for (int fp = 0; fp < nfp; ++fp)
329 {
330 const mfem::real_t radius =
331 ctx.axisymmetric ? face_radius[fp] : 0.0;
332 AssembleViscousFacePointKernel(
333 ctx, u_face, nor_face + fp*dim, w_minus[fp], w_plus[fp],
334 dprim_face_x, dprim_face_y, dprim_face_z, radius, fp, rhs_face);
335 }
336 }
337
338 template<typename ContextT>
339 MFEM_HOST_DEVICE inline static void ComputeFVFluxesKernel(const ContextT &ctx,
340 const mfem::real_t *el_u,
341 const mfem::real_t *elJac,
342 const mfem::real_t *el_metric_xi,
343 const mfem::real_t *el_metric_eta,
344 const mfem::real_t *el_metric_zeta,
345 mfem::real_t *el_dudt)
346 {
347 const int dim = ctx.dim;
348 const int Np_x = ctx.Np_x;
349 const int Np_y = ctx.Np_y;
350 const int Np_z = ctx.Np_z;
351 const int neq = ctx.num_equations;
352 const int npe = Np_x * Np_y * Np_z;
353 const mfem::real_t *qWgt = ctx.subcell_weights_d;
354
355 mfem::real_t flux_num[Theseus::MAXEQ];
356 mfem::real_t du_subcell[Theseus::MAXEQ];
357 mfem::real_t state1_local[Theseus::MAXEQ];
358 mfem::real_t state2_local[Theseus::MAXEQ];
359
360 for(int i = 0;i < npe*neq;i++)
361 el_dudt[i] = 0.0;
362
363 for (int k = 0; k < Np_z; k++)
364 {
365 for (int j = 0; j < Np_y; j++)
366 {
367 for(int q = 0; q < neq;q++){
368 du_subcell[q] = 0.0;
369 }
370 int id1 = k * Np_y * Np_x + j * Np_x;
371 Kernels::el_gather_state(el_u, npe, neq, id1, state1_local);
372 for (int i = 0; i < Np_x - 1; i++)
373 {
374 int id2 = id1 + 1;
375 Kernels::el_gather_state(el_u, npe, neq, id2, state2_local);
376 const mfem::real_t *nor = el_metric_xi + id2*dim;
377
378 ctx.iflux.ComputeFaceFlux(ctx.gas, state1_local,
379 state2_local, nor, flux_num);
380 for(int q = 0; q < neq;q++){
381 du_subcell[q] -= flux_num[q];
382 }
383 for(int q = 0; q < neq;q++){
384 du_subcell[q] /= (elJac[id1] * qWgt[i]);
385 }
386 Kernels::el_scatter_assign(du_subcell, npe, neq, id1, 1.0, el_dudt);
387 for(int q = 0; q < neq;q++){
388 du_subcell[q] = flux_num[q];
389 }
390 for(int q = 0;q < neq;q++){
391 state1_local[q] = state2_local[q];
392 }
393 id1 = id2;
394 }
395 for(int q = 0;q < neq;q++){
396 du_subcell[q] /= (elJac[id1] * qWgt[Np_x-1]);
397 }
398 Kernels::el_scatter_assign(du_subcell, npe, neq, id1, 1.0, el_dudt);
399 }
400 }
401
402 if (dim > 1)
403 {
404 for (int k = 0; k < Np_z; k++)
405 {
406 for (int i = 0; i < Np_x; i++)
407 {
408 for(int q = 0; q < neq;q++){
409 du_subcell[q] = 0.0;
410 }
411 int id1 = k * Np_y * Np_x + i;
412 Kernels::el_gather_state(el_u, npe, neq, id1,
413 state1_local);
414 for (int j = 0; j < Np_y - 1; j++)
415 {
416 int id2 = k * Np_y * Np_x + (j + 1) * Np_x + i;
417 Kernels::el_gather_state(el_u, npe, neq, id2,
418 state2_local);
419 const mfem::real_t *nor = el_metric_eta + id2*dim;
420 ctx.iflux.ComputeFaceFlux(ctx.gas, state1_local,
421 state2_local, nor, flux_num);
422 for(int q = 0;q < neq;q++){
423 du_subcell[q] -= flux_num[q];
424 }
425 for(int q = 0;q < neq;q++){
426 du_subcell[q] /= (elJac[id1] * qWgt[j]);
427 }
428 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
429 for(int q = 0;q < neq;q++){
430 du_subcell[q] = flux_num[q];
431 state1_local[q] = state2_local[q];
432 }
433 id1 = id2;
434 }
435 for(int q = 0;q < neq;q++){
436 du_subcell[q] /= (elJac[id1] * qWgt[Np_y - 1]);
437 }
438 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
439 }
440 }
441 if (dim > 2)
442 {
443 for (int j = 0; j < Np_y; j++)
444 {
445 for (int i = 0; i < Np_x; i++)
446 {
447 for(int q = 0; q < neq;q++){
448 du_subcell[q] = 0.0;
449 }
450 int id1 = j * Np_x + i;
451 Kernels::el_gather_state(el_u, npe, neq, id1,
452 state1_local);
453 for (int k = 0; k < Np_z - 1; k++)
454 {
455 int id2 = (k + 1) * Np_y * Np_x + j * Np_x + i;
456 Kernels::el_gather_state(el_u, npe, neq, id2,
457 state2_local);
458 const mfem::real_t *nor = el_metric_zeta + id2*dim;
459 ctx.iflux.ComputeFaceFlux(ctx.gas, state1_local,
460 state2_local, nor, flux_num);
461 for(int q = 0;q < neq;q++){
462 du_subcell[q] -= flux_num[q];
463 }
464 for(int q = 0;q < neq;q++){
465 du_subcell[q] /= (elJac[id1] * qWgt[k]);
466 }
467 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
468
469 for(int q = 0;q < neq;q++){
470 du_subcell[q] = flux_num[q];
471 state1_local[q] = state2_local[q];
472 }
473 id1 = id2;
474 }
475 for(int q = 0;q < neq;q++){
476 du_subcell[q] /= (elJac[id1] * qWgt[Np_z - 1]);
477 }
478 Kernels::el_scatter_add(du_subcell, npe, neq, id1, 1.0, el_dudt);
479 }
480 }
481 }
482 }
483 }
484
485
486 template<typename ContextType>
487 MFEM_HOST_DEVICE inline
488 static void AssembleViscousVolumePointKernel(
489 const ContextType &ctx, const mfem::real_t *el_u,
490 const mfem::real_t *elJac_d, const mfem::real_t *elMetric_d,
491 const mfem::real_t *elRadius_d,
492 const mfem::real_t *el_gradprim_x,
493 const mfem::real_t *el_gradprim_y,
494 const mfem::real_t *el_gradprim_z,
495 const int point, mfem::real_t *el_dudt)
496 {
497 const int Np_x = ctx.Np_x;
498 const int Np_y = ctx.Np_y;
499 const int Np_z = ctx.Np_z;
500 const int dim = ctx.dim;
501 const int neq = ctx.num_equations;
502 const int dof = ctx.ndof_scalar_el;
503 const mfem::real_t *Dhat_d = ctx.Dhat_d;
504 const int i = point % Np_x;
505 const int j = (point / Np_x) % Np_y;
506 const int k = point / (Np_x*Np_y);
507
508 mfem::real_t state[Theseus::MAXEQ] = {0.0};
509 mfem::real_t dqx[Theseus::MAXEQ] = {0.0};
510 mfem::real_t dqy[Theseus::MAXEQ] = {0.0};
511 mfem::real_t dqz[Theseus::MAXEQ] = {0.0};
512 mfem::real_t f_ref[Theseus::MAXEQ] = {0.0};
513 mfem::real_t dU_viscous[Theseus::MAXEQ] = {0.0};
514
515 for (int l = 0; l < Np_x; ++l)
516 {
517 const int sample = k*Np_y*Np_x + j*Np_x + l;
518 const mfem::real_t coefficient = Dhat_d[l + Np_x*i];
519 Kernels::el_gather_state(el_u, dof, neq, sample, state);
520 Kernels::el_gather_grad_state(
521 el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq,
522 sample, dqx, dqy, dqz);
523 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
524 ctx.gas, dim, neq, state, dqx, dqy, dqz,
525 elMetric_d + sample*dim*dim, f_ref, ctx.axisymmetric,
526 ctx.axisymmetric ? elRadius_d[sample] : 0.0);
527 for (int q = 0; q < neq; ++q)
528 {
529 dU_viscous[q] += coefficient*f_ref[q];
530 }
531 }
532
533 if (dim > 1)
534 {
535 for (int l = 0; l < Np_y; ++l)
536 {
537 const int sample = k*Np_y*Np_x + l*Np_x + i;
538 const mfem::real_t coefficient = Dhat_d[l + Np_y*j];
539 Kernels::el_gather_state(el_u, dof, neq, sample, state);
540 Kernels::el_gather_grad_state(
541 el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq,
542 sample, dqx, dqy, dqz);
543 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
544 ctx.gas, dim, neq, state, dqx, dqy, dqz,
545 elMetric_d + sample*dim*dim + dim, f_ref,
546 ctx.axisymmetric,
547 ctx.axisymmetric ? elRadius_d[sample] : 0.0);
548 for (int q = 0; q < neq; ++q)
549 {
550 dU_viscous[q] += coefficient*f_ref[q];
551 }
552 }
553 }
554
555 if (dim > 2)
556 {
557 for (int l = 0; l < Np_z; ++l)
558 {
559 const int sample = l*Np_y*Np_x + j*Np_x + i;
560 const mfem::real_t coefficient = Dhat_d[l + Np_z*k];
561 Kernels::el_gather_state(el_u, dof, neq, sample, state);
562 Kernels::el_gather_grad_state(
563 el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq,
564 sample, dqx, dqy, dqz);
565 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
566 ctx.gas, dim, neq, state, dqx, dqy, dqz,
567 elMetric_d + sample*dim*dim + 2*dim, f_ref,
568 ctx.axisymmetric,
569 ctx.axisymmetric ? elRadius_d[sample] : 0.0);
570 for (int q = 0; q < neq; ++q)
571 {
572 dU_viscous[q] += coefficient*f_ref[q];
573 }
574 }
575 }
576
577 Kernels::el_scatter_add(
578 dU_viscous, dof, neq, point, 1.0/elJac_d[point], el_dudt);
579 if (ctx.axisymmetric)
580 {
581 Kernels::el_gather_state(el_u, dof, neq, point, state);
582 Kernels::el_gather_grad_state(
583 el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq,
584 point, dqx, dqy, dqz);
585 mfem::real_t source[Theseus::MAXEQ] = {0.0};
586 if (!AddAxisymmetricViscousSourceAwayFromAxis(
587 ctx.gas, state, dqx, dqy, dqz, elRadius_d[point], source))
588 {
589 AddAxisymmetricViscousSourceAtAxis(
590 ctx, el_u, el_gradprim_x, el_gradprim_y, el_gradprim_z,
591 elRadius_d, elJac_d, elMetric_d, point, source);
592 }
593 Kernels::el_scatter_add(source, dof, neq, point, 1.0, el_dudt);
594 }
595 }
596
597 template<typename ContextType>
598 MFEM_HOST_DEVICE inline
599 static void AssembleViscousElementVolumeKernel(const ContextType &ctx,
600 const mfem::real_t *el_u,
601 const mfem::real_t *elJac_d,
602 const mfem::real_t *elMetric_d,
603 const mfem::real_t *elRadius_d,
604 const mfem::real_t *el_gradprim_x,
605 const mfem::real_t *el_gradprim_y,
606 const mfem::real_t *el_gradprim_z,
607 mfem::real_t *el_dudt)
608 {
609 const int Np_x = ctx.Np_x;
610 const int Np_y = ctx.Np_y;
611 const int Np_z = ctx.Np_z;
612 const int dim = ctx.dim;
613 const int neq = ctx.num_equations;
614 const int dof = Np_x * Np_y * Np_z;
615 const mfem::real_t *Dhat_d = ctx.Dhat_d;
616
617 // One source-point scratch
618 mfem::real_t state[Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
619 mfem::real_t dqx [Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
620 mfem::real_t dqy [Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
621 mfem::real_t dqz [Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
622
623 // flux(eq,dir)
624 // one transformed reference-direction flux vector
625 mfem::real_t f_ref[Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
626
627 for (int k = 0; k < Np_z; ++k)
628 {
629 for (int j = 0; j < Np_y; ++j)
630 {
631 for (int i = 0; i < Np_x; ++i)
632 {
633 const int id1 = k * Np_y * Np_x + j * Np_x + i;
634 const mfem::real_t J = elJac_d[id1];
635 const mfem::real_t jInv = 1.0/J;
636
637 mfem::real_t dU_viscous[Theseus::MAXEQ] = {0., 0., 0., 0., 0.};
638
639 // xi contribution
640 for (int l = 0; l < Np_x; ++l)
641 {
642 const int idl = k * Np_y * Np_x + j * Np_x + l;
643 const mfem::real_t c = Dhat_d[l + Np_x * i];
644
645 Kernels::el_gather_state(el_u, dof, neq, idl, state);
646 Kernels::el_gather_grad_state(el_gradprim_x, el_gradprim_y,
647 el_gradprim_z, dim, dof, neq, idl,
648 dqx, dqy, dqz);
649
650 const mfem::real_t *adj_row = elMetric_d + idl * dim * dim + 0 * dim;
651 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
652 ctx.gas, dim, neq, state, dqx, dqy, dqz, adj_row,
653 f_ref, ctx.axisymmetric,
654 ctx.axisymmetric ? elRadius_d[idl] : 0.0);
655 for (int q = 0; q < neq; ++q)
656 {
657 dU_viscous[q] += c * f_ref[q];
658 }
659 }
660
661 // eta contribution
662 if (dim > 1)
663 {
664 for (int l = 0; l < Np_y; ++l)
665 {
666 const int idl = k * Np_y * Np_x + l * Np_x + i;
667 const mfem::real_t c = Dhat_d[l + Np_y * j];
668
669 Kernels::el_gather_state(el_u, dof, neq, idl, state);
670 Kernels::el_gather_grad_state(el_gradprim_x, el_gradprim_y,
671 el_gradprim_z, dim, dof, neq, idl,
672 dqx, dqy, dqz);
673
674 const mfem::real_t *adj_row = elMetric_d + idl * dim * dim + 1 * dim;
675 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
676 ctx.gas, dim, neq, state, dqx, dqy, dqz, adj_row,
677 f_ref, ctx.axisymmetric,
678 ctx.axisymmetric ? elRadius_d[idl] : 0.0);
679
680 for (int q = 0; q < neq; ++q)
681 {
682 dU_viscous[q] += c * f_ref[q];
683 }
684 }
685 }
686
687 // zeta contribution
688 if (dim > 2)
689 {
690 for (int l = 0; l < Np_z; ++l)
691 {
692 const int idl = l * Np_y * Np_x + j * Np_x + i;
693 const mfem::real_t c = Dhat_d[l + Np_z * k];
694
695 Kernels::el_gather_state(el_u, dof, neq, idl, state);
696 Kernels::el_gather_grad_state(el_gradprim_x, el_gradprim_y,
697 el_gradprim_z, dim, dof, neq, idl,
698 dqx, dqy, dqz);
699
700 const mfem::real_t *adj_row = elMetric_d + idl * dim * dim + 2 * dim;
701 Theseus::NavierStokesFlux::compute_ref_viscous_flux(
702 ctx.gas, dim, neq, state, dqx, dqy, dqz, adj_row,
703 f_ref, ctx.axisymmetric,
704 ctx.axisymmetric ? elRadius_d[idl] : 0.0);
705
706 for (int q = 0; q < neq; ++q)
707 {
708 dU_viscous[q] += c * f_ref[q];
709 }
710 }
711 }
712 Kernels::el_scatter_add(dU_viscous, dof, neq, id1, jInv, el_dudt);
713 if (ctx.axisymmetric)
714 {
715 Kernels::el_gather_state(el_u, dof, neq, id1, state);
716 Kernels::el_gather_grad_state(
717 el_gradprim_x, el_gradprim_y, el_gradprim_z, dim,
718 dof, neq, id1, dqx, dqy, dqz);
719 mfem::real_t source[Theseus::MAXEQ] = {0.0};
720 if (!AddAxisymmetricViscousSourceAwayFromAxis(
721 ctx.gas, state, dqx, dqy, dqz,
722 elRadius_d[id1], source))
723 {
724 AddAxisymmetricViscousSourceAtAxis(
725 ctx, el_u, el_gradprim_x, el_gradprim_y,
726 el_gradprim_z, elRadius_d, elJac_d, elMetric_d,
727 id1, source);
728 }
729 Kernels::el_scatter_add(source, dof, neq, id1, 1.0,
730 el_dudt);
731 }
732 }
733 }
734 }
735 }
736
737 template <typename ContextType>
738 MFEM_HOST_DEVICE inline
739 static void AssembleGradVolumePointKernel(
740 const ContextType &ctx, const mfem::real_t *el_u,
741 const mfem::real_t *elJac_d, const mfem::real_t *elMetric_d,
742 const int point, mfem::real_t *el_grad_u[Theseus::MAXDIM])
743 {
744 const int Np_x = ctx.Np_x;
745 const int Np_y = ctx.Np_y;
746 const int neq = ctx.num_equations;
747 const int dim = ctx.dim;
748 const int dof = ctx.ndof_scalar_el;
749 const mfem::real_t *D_d = ctx.D_d;
750
751 const int i = point % Np_x;
752 const int j = (point / Np_x) % Np_y;
753 const int k = point / (Np_x * Np_y);
754
755 mfem::real_t dudxi[Theseus::MAXEQ] = {0.0};
756 mfem::real_t dudeta[Theseus::MAXEQ] = {0.0};
757 mfem::real_t dudzeta[Theseus::MAXEQ] = {0.0};
758
759 for (int l = 0; l < Np_x; ++l)
760 {
761 const int sample = k*Np_y*Np_x + j*Np_x + l;
762 const mfem::real_t coefficient = D_d[l + Np_x*i];
763 for (int q = 0; q < neq; ++q)
764 {
765 dudxi[q] += el_u[sample + q*dof] * coefficient;
766 }
767 }
768
769 if (dim > 1)
770 {
771 for (int l = 0; l < Np_y; ++l)
772 {
773 const int sample = k*Np_y*Np_x + l*Np_x + i;
774 const mfem::real_t coefficient = D_d[l + Np_y*j];
775 for (int q = 0; q < neq; ++q)
776 {
777 dudeta[q] += el_u[sample + q*dof] * coefficient;
778 }
779 }
780 }
781
782 if (dim > 2)
783 {
784 for (int l = 0; l < ctx.Np_z; ++l)
785 {
786 const int sample = l*Np_y*Np_x + j*Np_x + i;
787 const mfem::real_t coefficient = D_d[l + ctx.Np_z*k];
788 for (int q = 0; q < neq; ++q)
789 {
790 dudzeta[q] += el_u[sample + q*dof] * coefficient;
791 }
792 }
793 }
794
795 const mfem::real_t invJ = 1.0 / elJac_d[point];
796 const mfem::real_t *adj = elMetric_d + point*dim*dim;
797 for (int q = 0; q < neq; ++q)
798 {
799 if (dim == 1)
800 {
801 el_grad_u[0][point + q*dof] = invJ*dudxi[q]*adj[0];
802 }
803 else if (dim == 2)
804 {
805 el_grad_u[0][point + q*dof] =
806 invJ*(dudxi[q]*adj[0] + dudeta[q]*adj[2]);
807 el_grad_u[1][point + q*dof] =
808 invJ*(dudxi[q]*adj[1] + dudeta[q]*adj[3]);
809 }
810 else
811 {
812 el_grad_u[0][point + q*dof] =
813 invJ*(dudxi[q]*adj[0] + dudeta[q]*adj[3] +
814 dudzeta[q]*adj[6]);
815 el_grad_u[1][point + q*dof] =
816 invJ*(dudxi[q]*adj[1] + dudeta[q]*adj[4] +
817 dudzeta[q]*adj[7]);
818 el_grad_u[2][point + q*dof] =
819 invJ*(dudxi[q]*adj[2] + dudeta[q]*adj[5] +
820 dudzeta[q]*adj[8]);
821 }
822 }
823 }
824
825 template <typename ContextType>
826 MFEM_HOST_DEVICE inline
827 static void AssembleGradElementVolumeKernel(const ContextType &ctx,
828 const mfem::real_t *el_u,
829 const mfem::real_t *elJac_d,
830 const mfem::real_t *elMetric_d,
831 mfem::real_t *el_grad_u[Theseus::MAXDIM])
832 {
833 const int Np_x = ctx.Np_x;
834 const int Np_y = ctx.Np_y;
835 const int Np_z = ctx.Np_z;
836 const int neq = ctx.num_equations;
837 const int dim = ctx.dim;
838 const int dof = Np_x * Np_y * Np_z;
839 const mfem::real_t *D_d = ctx.D_d;
840
841 if(dim == 1){
842
843 // Keep MAX_EQ in mind later if neq can exceed 5.
844 mfem::real_t dudxi[Theseus::MAXEQ];
845
846 for (int i = 0; i < Np_x; ++i)
847 {
848 const int id = i;
849
850 for (int q = 0; q < neq; ++q)
851 {
852 dudxi[q] = 0.0;
853 }
854 // Reference-space derivatives.
855 for (int l = 0; l < Np_x; ++l)
856 {
857 const int id_x = l;
858 const mfem::real_t c_xi = D_d[i*Np_x + l]; // legacy D_T(l, i)
859
860 for (int q = 0; q < neq; ++q)
861 {
862 dudxi[q] += el_u[id_x + q * dof] * c_xi;
863 }
864 }
865
866 const mfem::real_t invJ = 1.0 / elJac_d[id];
867 const mfem::real_t *adj = elMetric_d + id * dim * dim;
868
869 for (int q = 0; q < neq; ++q)
870 {
871 el_grad_u[0][id + q * dof] = invJ * (dudxi[q] * adj[0]);
872 }
873 }
874 } else if(dim == 2){
875 mfem::real_t dudxi[Theseus::MAXEQ];
876 mfem::real_t dudeta[Theseus::MAXEQ];
877
878 for (int j = 0; j < Np_y; ++j)
879 {
880 for (int i = 0; i < Np_x; ++i)
881 {
882 const int id = j * Np_x + i;
883
884 for (int q = 0; q < neq; ++q)
885 {
886 dudxi[q] = 0.0;
887 dudeta[q] = 0.0;
888 }
889
890 // Reference-space derivatives.
891 for (int l = 0; l < Np_x; ++l)
892 {
893 const int id_x = j * Np_x + l;
894 const int id_y = l * Np_x + i;
895
896 const mfem::real_t c_xi = D_d[l + Np_x * i]; // legacy D_T(l,i)
897 const mfem::real_t c_eta = D_d[l + Np_x * j]; // legacy D_T(l,j)
898
899 for (int q = 0; q < neq; ++q)
900 {
901 dudxi[q] += el_u[id_x + q * dof] * c_xi;
902 dudeta[q] += el_u[id_y + q * dof] * c_eta;
903 }
904 }
905
906 const mfem::real_t invJ = 1.0 / elJac_d[id];
907 const mfem::real_t *adj = elMetric_d + id * dim * dim;
908
909 // adj stored row-major per point:
910 // [ adj[0] adj[1] ]
911 // [ adj[2] adj[3] ]
912 for (int q = 0; q < neq; ++q)
913 {
914 el_grad_u[0][id + q * dof] = invJ * (dudxi[q] * adj[0] + dudeta[q] * adj[2]);
915 el_grad_u[1][id + q * dof] = invJ * (dudxi[q] * adj[1] + dudeta[q] * adj[3]);
916 }
917 }
918 }
919 } else if (dim == 3) {
920
921 mfem::real_t dudxi[Theseus::MAXEQ];
922 mfem::real_t dudeta[Theseus::MAXEQ];
923 mfem::real_t dudzeta[Theseus::MAXEQ];
924
925 for (int k = 0; k < Np_z; ++k)
926 {
927 for (int j = 0; j < Np_y; ++j)
928 {
929 for (int i = 0; i < Np_x; ++i)
930 {
931 const int id = k * Np_x * Np_y + j * Np_x + i;
932
933 for (int q = 0; q < neq; ++q)
934 {
935 dudxi[q] = 0.0;
936 dudeta[q] = 0.0;
937 dudzeta[q] = 0.0;
938 }
939
940 // Reference-space derivatives.
941 for (int l = 0; l < Np_x; ++l)
942 {
943 const int id_x = k * Np_x * Np_y + j * Np_x + l;
944 const int id_y = k * Np_x * Np_y + l * Np_x + i;
945 const int id_z = l * Np_x * Np_y + j * Np_x + i;
946 const mfem::real_t c_xi = D_d[l + Np_x * i]; // legacy D_T(l,i)
947 const mfem::real_t c_eta = D_d[l + Np_x * j]; // legacy D_T(l,j)
948 const mfem::real_t c_zeta = D_d[l + Np_x * k]; // legacy D_T(l,k)
949 for (int q = 0; q < neq; ++q)
950 {
951 dudxi[q] += el_u[id_x + q * dof] * c_xi;
952 dudeta[q] += el_u[id_y + q * dof] * c_eta;
953 dudzeta[q] += el_u[id_z + q * dof] * c_zeta;
954 }
955 }
956
957 const mfem::real_t invJ = 1.0 / elJac_d[id];
958 const mfem::real_t *adj = elMetric_d + id * dim * dim;
959
960 // adj stored row-major per point:
961 // [ adj[0] adj[1] adj[2] ]
962 // [ adj[3] adj[4] adj[5] ]
963 // [ adj[6] adj[7] adj[8] ]
964 for (int q = 0; q < neq; ++q)
965 {
966 el_grad_u[0][id + q * dof] = invJ * (dudxi[q] * adj[0] +
967 dudeta[q] * adj[3] +
968 dudzeta[q] * adj[6]);
969
970 el_grad_u[1][id + q * dof] = invJ * (dudxi[q] * adj[1] +
971 dudeta[q] * adj[4] +
972 dudzeta[q] * adj[7]);
973
974 el_grad_u[2][id + q * dof] = invJ * (dudxi[q] * adj[2] +
975 dudeta[q] * adj[5] +
976 dudzeta[q] * adj[8]);
977 }
978 }
979 }
980 }
981 }
982 }
983
984 template <typename ContextT>
985 MFEM_HOST_DEVICE inline
986 static void AssembleGradInteriorFacePointKernel(
987 const ContextT &ctx,
988 const mfem::real_t *u_face,
989 const mfem::real_t *nor_point,
990 const mfem::real_t w_minus,
991 const mfem::real_t w_plus,
992 const int fp,
993 mfem::real_t *rhs_face[Theseus::MAXDIM])
994 {
995 const int neq = ctx.num_equations;
996 const int dim = ctx.dim;
997
998 mfem::real_t jump[Theseus::MAXEQ];
999
1000 for (int q = 0; q < neq; ++q)
1001 {
1002 jump[q] = mfem::real_t(0.5) *
1003 (u_face[ctx.iface_idx(1, fp, q)] -
1004 u_face[ctx.iface_idx(0, fp, q)]);
1005 }
1006
1007 for (int idim = 0; idim < dim; ++idim){
1008 mfem::real_t *rhs_d = rhs_face[idim];
1009 const mfem::real_t n_d = nor_point[idim];
1010 for (int q = 0; q < neq; ++q)
1011 {
1012 const mfem::real_t f_d = jump[q]*n_d;
1013 rhs_d[ctx.iface_idx(0, fp, q)] = w_minus * f_d;
1014 rhs_d[ctx.iface_idx(1, fp, q)] = w_plus * f_d;
1015 }
1016 }
1017 }
1018
1019 template <typename ContextT>
1020 MFEM_HOST_DEVICE inline
1021 static void AssembleGradInteriorFaceKernel(const ContextT &ctx,
1022 const mfem::real_t *u_face,
1023 const mfem::real_t *nor_face,
1024 const mfem::real_t *w_minus,
1025 const mfem::real_t *w_plus,
1026 mfem::real_t *rhs_face[Theseus::MAXDIM])
1027 {
1028 const int nfp = ctx.num_face_points;
1029 const int dim = ctx.dim;
1030
1031 for (int fp = 0; fp < nfp; ++fp)
1032 {
1033 AssembleGradInteriorFacePointKernel(
1034 ctx, u_face, nor_face + fp*dim, w_minus[fp], w_plus[fp], fp,
1035 rhs_face);
1036 }
1037 }
1038
1039 template <typename DeviceCacheT>
1040 MFEM_HOST_DEVICE inline
1041 static void AssembleGradBoundaryPointKernel(const DeviceCacheT &dc,
1042 const Theseus::BCDescriptor &bc,
1043 const mfem::real_t *u_face,
1044 const mfem::real_t *nor_point,
1045 const mfem::real_t scale,
1046 const int fp,
1047 mfem::real_t *rhs_face[Theseus::MAXDIM])
1048 {
1049 const int dim = dc.dim;
1050 const int nfp = dc.num_face_points;
1051 const int neq = dc.num_equations;
1052
1053 mfem::real_t state1[Theseus::MAXEQ];
1054 mfem::real_t fluxN[Theseus::MAXEQ];
1055 mfem::real_t flux_dir[Theseus::MAXEQ];
1056
1057 Theseus::Kernels::el_gather_state(u_face, nfp, neq, fp, state1);
1058
1059 Theseus::BC::ComputeBdrFaceGradFlux(dc, bc, state1, fluxN);
1060
1061 for(int idim = 0;idim < dim;idim++){
1062 for(int q = 0;q < neq;q++){
1063 flux_dir[q] = fluxN[q]*nor_point[idim];
1064 }
1065 Theseus::Kernels::el_scatter_add(flux_dir, nfp, neq, fp, scale, rhs_face[idim]);
1066 }
1067 }
1068
1069 };
1070}
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
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