Theseus
Compressible flow solver
Loading...
Searching...
No Matches
NSOperator_impl.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
7namespace Theseus
8{
9
10 // Device version of ComputeGlobalEntropyVector
11 template<typename PhysicsT>
12 void NSOperator<PhysicsT>::ComputeEntropyState(const mfem::Vector &u, mfem::Vector &e) const
13 {
14 Theseus::ScopedTimer timer("ComputeEntropyState");
15
16 // This block is executed by the host
17 const int nval_restr = operator_cache.restr_v->Height();
18
19 // Copy the device cache so that it is not member data
20 auto dc = device_cache;
21
22 // Device cache parameters
23 const int ne = dc.num_elements;
24 const int ndof = dc.ndof_scalar_el;
25 const int neq = dc.num_equations;
26 const int npts = ndof * ne;
27
28 MFEM_ASSERT(nval_restr == npts*neq, "Unexpected size in ComputeEntropyState");
29
30 auto gas = dc.gas;
31
32 if (operator_cache.uVol.Size() != nval_restr){
33 operator_cache.uVol.SetSize(nval_restr);
34 operator_cache.uVol.UseDevice();
35 }
36 mfem::Vector &restrU(operator_cache.uVol);
37 if (operator_cache.uVol.Size() != nval_restr){
38 operator_cache.uVol.SetSize(nval_restr);
39 operator_cache.uVol.UseDevice();
40 }
41 mfem::Vector &restrE(operator_cache.sVol);
42 if(restrE.Size() != nval_restr){
43 restrE.SetSize(nval_restr);
44 restrE.UseDevice();
45 }
46 mfem::real_t *eState_d = restrE.Write();
47 operator_cache.restr_v->Mult(u, restrU);
48
49 if(e.Size() != u.Size()){
50 e.SetSize(u.Size());
51 e.UseDevice();
52 }
53
54 const mfem::real_t *restrU_d = restrU.Read();
55 const int estride = ndof*neq;
56
57 // Inside the FORALL below, executed on device
58 mfem::forall(npts, [=] MFEM_HOST_DEVICE (int pt)
59 {
60 const int elno = pt / ndof;
61 const int ept = pt % ndof;
62 const int eoff = elno * estride;
63 const mfem::real_t *u_el = restrU_d + eoff;
64
65 mfem::real_t elUstate[Theseus::MAXEQ];
66 Theseus::Kernels::el_gather_state(u_el, ndof, neq, ept, elUstate);
67 Theseus::PointStateView S{elUstate};
68
69 mfem::real_t elEState[Theseus::MAXEQ];
70 Theseus::PointStateViewRW E{elEState};
71 gas.entropy_state(S, E);
72 mfem::real_t *e_el = eState_d + eoff;
73 Theseus::Kernels::el_scatter_assign(elEState, ndof, neq, ept, 1.0, e_el);
74 });
75
76 operator_cache.restr_v->MultTranspose(restrE, e);
77
78 }
79
80 // This routine replaces the gradient of the entropy stored in gradState with the gradient
81 // of the primitive variables so that all the data in gradState is replaced.
82 template<typename PhysicsT>
83 void NSOperator<PhysicsT>::ComputeGradPrimFromGradEntropy(const mfem::Vector &u, std::vector<mfem::Vector *> &gradEntropy) const
84 {
85 Theseus::ScopedTimer timer("GradEntropyToGradPrim");
86 // This block is executed by the host
87 const int nval_restr = operator_cache.restr_v->Height();
88
89 // Copy the device cache so that it is not member data
90 auto dc = device_cache;
91
92 // Device cache parameters
93 const int ne = dc.num_elements;
94 const int ndof = dc.ndof_scalar_el;
95 const int neq = dc.num_equations;
96 const int npts = ndof * ne;
97 const int dim = dc.dim;
98
99 MFEM_ASSERT(nval_restr == npts*neq, "Unexpected size in ComputeEntropyState");
100
101 auto gas = dc.gas;
102
103 if (operator_cache.uVol.Size() != nval_restr){
104 operator_cache.uVol.SetSize(nval_restr);
105 operator_cache.uVol.UseDevice();
106 }
107 mfem::Vector &restr_state(operator_cache.uVol);
108 if(!operator_cache.u_vol_restr_ready){
109 operator_cache.restr_v->Mult(u, restr_state);
110 operator_cache.u_vol_restr_ready = true;
111 }
112 const mfem::real_t *restr_state_d = restr_state.Read();
113 const int estride = ndof*neq;
114
115 // Leave this temporary for now
116 if(operator_cache.volAux.Size() != nval_restr){
117 operator_cache.volAux.SetSize(nval_restr);
118 operator_cache.volAux.UseDevice();
119 }
120 mfem::Vector &restr_grad_prim_dir(operator_cache.volAux);
121
122 for(int idim = 0;idim < dim;idim++){
123
124 mfem::Vector &grad_state_dir(*gradEntropy[idim]);
125 operator_cache.restr_v->Mult(grad_state_dir, restr_grad_prim_dir);
126 mfem::real_t *grad_prim_dir_d = restr_grad_prim_dir.Write();
127
128 // Inside the FORALL below, executed on device
129 mfem::forall(npts, [=] MFEM_HOST_DEVICE (int pt)
130 {
131 const int e = pt / ndof;
132 const int ept = pt % ndof;
133 const int eoff = e * estride;
134 const mfem::real_t *u_el = restr_state_d + eoff;
135 mfem::real_t *grad_prim_el = grad_prim_dir_d + eoff;
136
137 mfem::real_t el_U[Theseus::MAXEQ];
138 Theseus::Kernels::el_gather_state(u_el, ndof, neq, ept, el_U);
139 Theseus::PointStateView CV{el_U};
140
141 mfem::real_t el_gradS[Theseus::MAXEQ];
142 Theseus::Kernels::el_gather_state(grad_prim_el, ndof, neq, ept, el_gradS);
143 Theseus::PointStateView dS{el_gradS};
144
145 mfem::real_t el_gradP[Theseus::MAXEQ];
146 Theseus::PointStateViewRW dP{el_gradP};
147 gas.grad_entropy_to_grad_prim(CV, dS, dP);
148 if (dc.axisymmetric)
149 {
150 ProjectAxisPrimitiveGradientDirection(gas, dc.elRadius_d[pt], idim, el_gradP);
151 }
152 Kernels::el_scatter_assign(el_gradP, ndof, neq, ept, 1.0, grad_prim_el);
153
154 });
155
156 operator_cache.restr_v->MultTranspose(restr_grad_prim_dir, grad_state_dir);
157
158 }
159 }
160
161 template<typename PhysicsT>
162 mfem::real_t NSOperator<PhysicsT>::FlowMult(const mfem::Vector &u, mfem::Vector &pdudt) const
163 {
164 Theseus::ScopedTimer timer("NSRHS");
165
166 const int dim = operator_cache.dim;
167 {
168 mfem::Vector &entropyState(operator_cache.entropyState);
169 {
170 if (entropyState.Size() != u.Size()){
171 entropyState.SetSize(u.Size());
172 entropyState.UseDevice();
173 }
174 ComputeEntropyState(u, entropyState);
175 }
176 std::vector<mfem::Vector *> gradPrim(dim);
177 {
178 // grad_u is a vector of parallel grid functions
179 // this bit grabs an mfem::Vector ref.
180 // Note that incoming grad_u is really grad of entropy,
181 // which we pack into gradPrim, and then call a
182 // function which overwrites the entropy gradient
183 // with the primitive gradient.
184 for(int idim = 0;idim < dim;idim++){
185 gradPrim[idim] = &(*grad_u[idim]);
186 }
187 GradOperator(entropyState, gradPrim);
188 ComputeGradPrimFromGradEntropy(u, gradPrim);
189 }
190 max_char_speed = MultCNS(u, gradPrim, pdudt);
191 }
192 return max_char_speed;
193 }
194
195 template<typename PhysicsT>
196 void NSOperator<PhysicsT>::GradOperator_Volume(const mfem::Vector &pu,
197 std::vector<mfem::Vector *> &p_grad_u) const
198 {
199 Theseus::ScopedTimer timer("GradOperator_Volume");
200
201 const int dim = operator_cache.dim;
202 const int restr_size = operator_cache.restr_v->Height();
203
204 if(operator_cache.sVol.Size() != restr_size){
205 operator_cache.sVol.SetSize(restr_size);
206 operator_cache.sVol.UseDevice();
207 }
208 mfem::Vector &Ue(operator_cache.sVol);
209 mfem::real_t *dU_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
210 mfem::real_t *pgrad_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
211 if (operator_cache.gradVol.size() != dim){
212 operator_cache.gradVol.resize(dim);
213 for(int idim = 0;idim < dim;idim++){
214 operator_cache.gradVol[idim].SetSize(restr_size);
215 operator_cache.gradVol[idim].UseDevice();
216 }
217 }
218 std::vector<mfem::Vector> &dUe(operator_cache.gradVol);
219 for(int idim = 0;idim < dim;idim++){
220 dU_d[idim] = dUe[idim].Write();
221 pgrad_d[idim] = p_grad_u[idim]->Write();
222 }
223
224 mfem::forall(restr_size, [=] MFEM_HOST_DEVICE (int i)
225 {
226 for(int idim = 0;idim < dim;idim++){
227 dU_d[idim][i] = mfem::real_t(0);
228 pgrad_d[idim][i] = mfem::real_t(0);
229 }
230 });
231
232 operator_cache.restr_v->Mult(pu, Ue);
233 const mfem::real_t *Ue_d = Ue.Read();
234
235 auto dc = device_cache;
236
237 const int ne = dc.num_elements;
238 const int ndof = dc.ndof_scalar_el;
239 const int neq = dc.num_equations;
240 const int estride = ndof * neq;
241 const int jac_stride = ndof;
242 const int metric_stride = ndof * dc.dim * dc.dim;
243
244 const mfem::real_t *elJac_d = dc.elJac_d;
245 const mfem::real_t *elMetric_d = dc.elMetric_d;
246
247 mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
248 {
249 const mfem::real_t *u_el = Ue_d + e * estride;
250 mfem::real_t *du_el_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
251 for(int idim = 0;idim < dim;idim++){
252 du_el_d[idim] = dU_d[idim] + e*estride;
253 }
254
255 const mfem::real_t *jac_el = elJac_d + e * jac_stride;
256 const mfem::real_t *metric_el = elMetric_d + e * metric_stride;
257
258 Theseus::DGSEMIntegrator::AssembleGradElementVolumeKernel(dc, u_el, jac_el, metric_el,
259 du_el_d);
260 });
261
262 for(int idim = 0;idim < dim;idim++){
263 operator_cache.restr_v->AddMultTranspose(dUe[idim], *p_grad_u[idim]);
264 }
265
266 }
267
268 template<typename PhysicsT>
269 void NSOperator<PhysicsT>::GradOperator_BoundaryFaces(const mfem::Vector &pu,
270 std::vector<mfem::Vector *> &p_grad_u) const
271 {
272 Theseus::ScopedTimer timer("GradOperator_BoundaryFaces");
273
274 auto dc = device_cache;
275 const int dim = dc.dim;
276 const int neq = dc.num_equations;
277 const int nfp = dc.num_face_points;
278 const int face_size = nfp * neq;
279 const int restr_size = operator_cache.restr_b->Height();
280 const int nfaces_restr = restr_size / face_size;
281 const int norm_size = nfp * dim;
282 const int npoints_bnd = nfaces_restr * nfp;
283 const int psize = pu.Size();
284
285 if(restr_size == 0){
286 return;
287 }
288
289 mfem::Vector &u_faces(operator_cache.sBnd);
290 if(u_faces.Size() != restr_size){
291 u_faces.SetSize(restr_size);
292 u_faces.UseDevice();
293 }
294
295 std::vector<mfem::Vector> &rhs_faces(operator_cache.gradBnd);
296 if(rhs_faces.size() != dim){
297 rhs_faces.resize(dim);
298 for(int idim = 0;idim < dim;idim++){
299 rhs_faces[idim].SetSize(restr_size);
300 rhs_faces[idim].UseDevice();
301 }
302 }
303
304 mfem::Vector &duBnd(operator_cache.duBnd);
305 if(duBnd.Size() != psize){
306 duBnd.SetSize(psize);
307 duBnd.UseDevice();
308 }
309
310 operator_cache.restr_b->Mult(pu, u_faces);
311
312 const mfem::real_t *u_d = u_faces.Read();
313 const mfem::real_t *nor_d = dc.bnd_nor_d;
314 const mfem::real_t *wt_d = dc.bnd_wt_d;
315 const int *bnd_marker_index_d = dc.bnd_marker_index_d;
316
317 mfem::real_t *rhs_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
318 mfem::real_t *du_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
319 for(int idim = 0;idim < dim;idim++){
320 rhs_d[idim] = rhs_faces[idim].Write();
321 du_d[idim] = duBnd.Write();
322 }
323
324 for (int idim = 0; idim < dim; ++idim) {
325 mfem::real_t *rd = rhs_d[idim];
326 mfem::forall(restr_size, [=] MFEM_HOST_DEVICE (int i) { rd[i] = mfem::real_t(0); });
327 }
328 for (int idim = 0; idim < 1; ++idim) {
329 mfem::real_t *dud = du_d[idim];
330 mfem::forall(psize, [=] MFEM_HOST_DEVICE (int i) { dud[i] = mfem::real_t(0); });
331 }
332
333 mfem::forall(npoints_bnd, [=] MFEM_HOST_DEVICE (int p)
334 {
335 const int f = p / nfp;
336 const int fp = p % nfp;
337
338 const int bnd_face_marker_index = bnd_marker_index_d[f];
339 if (bnd_face_marker_index < 0)
340 {
341 return;
342 }
343
344 const int bc_index = bnd_face_marker_index; // same convention as inviscid device path for now
345 if (bc_index < 0)
346 {
347 return;
348 }
349
350 const Theseus::BCDescriptor &bc = dc.bc_descr_d[bc_index];
351 if (bc.type == int(Theseus::BCType::Invalid))
352 {
353 return;
354 }
355
356 const int face_offset = f * face_size;
357 const int norm_offset = f * norm_size;
358 const int w_offset = f * nfp;
359
360 const mfem::real_t *u_face_d = u_d + face_offset;
361
362 const mfem::real_t *nor_face_d = nor_d + norm_offset;
363 const mfem::real_t *nor_point = nor_face_d + fp * dim;
364
365 // Legacy one-sided boundary lifting uses +1/(w0*J1)
366 const mfem::real_t scale = wt_d[w_offset + fp];
367
368 mfem::real_t *rhs_face[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
369 for(int idim = 0;idim < dim;idim++){
370 rhs_face[idim] = rhs_d[idim] + face_offset;
371 }
372
373 Theseus::DGSEMIntegrator::AssembleGradBoundaryPointKernel(dc, bc,
374 u_face_d,
375 nor_point,
376 scale,
377 fp,
378 rhs_face);
379 });
380
381 for(int idim = 0;idim < dim;idim++){
382 operator_cache.restr_b->MultTranspose(rhs_faces[idim], duBnd);
383 *p_grad_u[idim] += duBnd;
384 }
385
386 }
387
388 template<typename PhysicsT>
389 void NSOperator<PhysicsT>::GradOperator_InteriorFaces(const mfem::Vector &pu,
390 std::vector<mfem::Vector *> &p_grad_u) const
391 {
392 Theseus::ScopedTimer timer("GradOperator_InteriorFaces");
393
394 auto dc = device_cache;
395 const int dim = dc.dim;
396 const int psize = pu.Size();
397 const int restr_size = operator_cache.restr_f->Height();
398 const int neq = dc.num_equations;
399 const int nfp = dc.num_face_points;
400 const int nfaces = restr_size / (2 * nfp * neq);
401 const int face_size = 2 * nfp * neq;
402 const int norm_size = nfp * dim;
403
404 mfem::Vector &u_faces(operator_cache.sInt);
405 if(u_faces.Size() != restr_size){
406 u_faces.SetSize(restr_size);
407 u_faces.UseDevice();
408 }
409
410 std::vector<mfem::Vector> &rhs_faces(operator_cache.gradInt);
411 if(rhs_faces.size() != dim){
412 rhs_faces.resize(dim);
413 for(int idim = 0;idim < dim;idim++){
414 rhs_faces[idim].SetSize(restr_size);
415 rhs_faces[idim].UseDevice();
416 }
417 }
418
419 mfem::Vector &duInt(operator_cache.duInt);
420 if(duInt.Size() != psize){
421 duInt.SetSize(psize);
422 duInt.UseDevice();
423 }
424
425 operator_cache.restr_f->Mult(pu, u_faces);
426 const mfem::real_t *u_d = u_faces.Read();
427
428 mfem::real_t *rhs_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
429 mfem::real_t *du_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
430 for(int idim = 0;idim < dim;idim++){
431 rhs_d[idim] = rhs_faces[idim].Write();
432 // We only need 1dim at a time, lets only use 1 temp
433 du_d[idim] = duInt.Write();
434 }
435
436 for (int idim = 0; idim < dim; ++idim) {
437 mfem::real_t *rd = rhs_d[idim];
438 mfem::forall(restr_size, [=] MFEM_HOST_DEVICE (int i) { rd[i] = mfem::real_t(0); });
439 }
440 // HardCode to dim=1 for now - temporary is only 1d
441 for (int idim = 0; idim < 1; ++idim) {
442 mfem::real_t *dud = du_d[idim];
443 mfem::forall(psize, [=] MFEM_HOST_DEVICE (int i) { dud[i] = mfem::real_t(0); });
444 }
445
446 const mfem::real_t *nor_d = dc.nor_d;
447 const mfem::real_t *wm_d = dc.fw_minus_d;
448 const mfem::real_t *wp_d = dc.fw_plus_d;
449
450 mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int f)
451 {
452 const int face_offset = f * face_size;
453 const int norm_offset = f * norm_size;
454 const int w_offset = f * nfp;
455
456 const mfem::real_t *u_face_d = u_d + face_offset;
457 const mfem::real_t *nor_face_d = nor_d + norm_offset;
458 const mfem::real_t *w_minus_d = wm_d + w_offset;
459 const mfem::real_t *w_plus_d = wp_d + w_offset;
460
461 mfem::real_t *rhs_face[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
462 for(int idim = 0;idim < dim;idim++){
463 rhs_face[idim] = rhs_d[idim] + face_offset;
464 }
465
466 Theseus::DGSEMIntegrator::AssembleGradInteriorFaceKernel(dc,
467 u_face_d,
468 nor_face_d,
469 w_minus_d,
470 w_plus_d,
471 rhs_face);
472 });
473
474 for(int idim = 0;idim < dim;idim++){
475 operator_cache.restr_f->MultTranspose(rhs_faces[idim], duInt);
476 *p_grad_u[idim] += duInt;
477 }
478
479 }
480
481 template<typename PhysicsT>
482 void NSOperator<PhysicsT>::GradOperator(const mfem::Vector &u,
483 std::vector<mfem::Vector *> &grad_u) const
484 {
485 Theseus::ScopedTimer timer("GradOperator");
486 const int dim = operator_cache.dim;
487 const mfem::Vector &pu = this->Prolongate(u);
488 std::vector<mfem::Vector *> p_grad_(dim);
489 if (this->P)
490 {
491 const int psize = this->P->Height();
492 if(operator_cache.pGrad.size() != dim){
493 operator_cache.pGrad.resize(dim);
494 for(int idim = 0;idim < dim;idim++){
495 operator_cache.pGrad[idim].SetSize(psize);
496 operator_cache.pGrad[idim].UseDevice();
497 }
498 }
499 for(int idim = 0;idim < dim;idim++){
500 p_grad_[idim] = &(operator_cache.pGrad[idim]);
501 }
502 }
503 std::vector<mfem::Vector *> &p_grad_u = this->P ? p_grad_ : grad_u;
504
505 MFEM_ASSERT(p_grad_u.size() == dim, "Size mismatch for gradient storage");
506 MFEM_ASSERT(grad_u.size() == dim, "Size mismatch for gradient storage");
507
508 GradOperator_Volume(pu, p_grad_u);
509
510 GradOperator_InteriorFaces(pu, p_grad_u);
511
512 GradOperator_BoundaryFaces(pu, p_grad_u);
513
514 if (this->Serial())
515 {
516 if (this->cP)
517 {
518 for(int idim = 0;idim < dim;idim++){
519 this->cP->MultTranspose(*p_grad_u[idim], *grad_u[idim]);
520 }
521 }
522 }
523 else
524 {
525 if(this->P){
526 for(int idim = 0;idim < dim;idim++){
527 this->P->MultTranspose(*p_grad_u[idim], *grad_u[idim]);
528 }
529 }
530 }
531
532 const int N = this->ess_tdof_list.Size();
533 const auto idx = this->ess_tdof_list.Read();
534 // std::cout << "N ZERO = " << N << std::endl;
535 for(int idim = 0;idim < dim;idim++){
536 auto gradu_dim_d = grad_u[idim]->ReadWrite();
537 mfem::forall(N, [=] MFEM_HOST_DEVICE (int i) { gradu_dim_d[idx[i]] = 0.0; });
538 }
539 }
540
541 template<typename PhysicsT>
542 mfem::real_t NSOperator<PhysicsT>::MultCNS_InteriorFaces(const mfem::Vector &pu,
543 const std::vector<mfem::Vector *> &p_grad_prim,
544 mfem::Vector &pdudt) const
545 {
546 Theseus::ScopedTimer timer("MultCNS_InteriorFaces");
547
548 auto dc = device_cache;
549 const int dim = dc.dim;
550 const int neq = dc.num_equations;
551 const int nfp = dc.num_face_points;
552 const int nfaces = operator_cache.restr_f->Height() / (nfp * neq * 2); // (+/-)
553 const int face_size = 2*nfp*neq;
554 const int norm_size = nfp*dim;
555
556 const int restr_size = operator_cache.restr_f->Height();
557 mfem::Vector &int_u(operator_cache.uInt);
558 if(int_u.Size() != restr_size){
559 int_u.SetSize(restr_size);
560 int_u.UseDevice();
561 }
562
563 mfem::Vector &rhs_faces(operator_cache.rhsInt);
564 if(rhs_faces.Size() != restr_size){
565 rhs_faces.SetSize(restr_size);
566 rhs_faces.UseDevice();
567 }
568
569 mfem::Vector &faces_dudt(operator_cache.dudtInt);
570 if(faces_dudt.Size() != pdudt.Size()){
571 faces_dudt.SetSize(pdudt.Size());
572 faces_dudt.UseDevice();
573 }
574 // faces_dudt = pdudt;
575
576 const mfem::real_t *grad_prim_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
577 std::vector<mfem::Vector> &int_grad_prim(operator_cache.gradInt);
578
579 if(int_grad_prim.size() != dim){
580 int_grad_prim.resize(dim);
581 for(int idim = 0;idim < dim;idim++){
582 int_grad_prim[idim].SetSize(restr_size);
583 int_grad_prim[idim].UseDevice();
584 }
585 }
586
587 for(int idim = 0;idim < dim;idim++){
588 operator_cache.restr_f->Mult(*p_grad_prim[idim], int_grad_prim[idim]);
589 grad_prim_d[idim] = int_grad_prim[idim].Read();
590 }
591
592 // If zeroed before accumulation, do it explicitly on device:
593 // Potentially, this is not needed at all since I think we overwrite everything
594 {
595 mfem::real_t *d = rhs_faces.Write();
596 mfem::forall(rhs_faces.Size(), [=] MFEM_HOST_DEVICE (int i) { d[i] = mfem::real_t(0); });
597 }
598
599 if(!operator_cache.u_int_restr_ready){
600 operator_cache.restr_f->Mult(pu, int_u);
601 operator_cache.u_int_restr_ready = true;
602 }
603
604 const mfem::real_t *u_d = int_u.Read();
605 mfem::real_t *rhs_d = rhs_faces.Write();
606 const mfem::real_t *nor_d = dc.nor_d; // size nfaces*nfp*dim
607 const mfem::real_t *inv1_d = dc.fw_minus_d; // size nfaces*nfp
608 const mfem::real_t *inv2_d = dc.fw_plus_d; // size nfaces*nfp
609 const mfem::real_t *face_radius_d = dc.face_radius_d;
610
611 mfem::real_t *ws_d = dc.ifWaveSpeed_d;
612
613 mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int i)
614 {
615 const int face_offset = i*face_size;
616 const int n_offset = i*norm_size;
617 const int w_offset = i*nfp;
618
619 const mfem::real_t *u_face_d = u_d + face_offset;
620 mfem::real_t *rhs_face_d = rhs_d + face_offset;
621 const mfem::real_t *nor_face_d = nor_d + n_offset;
622 const mfem::real_t *w_minus_d = inv1_d + w_offset;
623 const mfem::real_t *w_plus_d = inv2_d + w_offset;
624 const mfem::real_t *radius_face_d = dc.axisymmetric ?
625 face_radius_d + w_offset : nullptr;
626 const mfem::real_t *dprim_face_x = (dim > 0) ? grad_prim_d[0] + face_offset : nullptr;
627 const mfem::real_t *dprim_face_y = (dim > 1) ? grad_prim_d[1] + face_offset : nullptr;
628 const mfem::real_t *dprim_face_z = (dim > 2) ? grad_prim_d[2] + face_offset : nullptr;
629
630 // Call one fused kernel for inviscid and viscous facial terms
631 mfem::real_t ws = Theseus::DGSEMIntegrator::AssembleViscousElementFaceKernel(dc, u_face_d, nor_face_d,
632 w_minus_d, w_plus_d,
633 dprim_face_x,
634 dprim_face_y,
635 dprim_face_z,
636 radius_face_d,
637 rhs_face_d);
638 ws_d[i] = ws;
639 });
640
641 operator_cache.restr_f->MultTranspose(rhs_faces, faces_dudt);
642 pdudt += faces_dudt; // on device?
643
644 // Finish up on the host:
645 // - Reduce for rank-local max_char_speed
646 const mfem::real_t *ws = operator_cache.ifWaveSpeed.HostRead();
647 mfem::real_t max_char_speed_facial = 0.0;
648 for(int f = 0;f < operator_cache.num_interior_faces;f++)
649 {
650 max_char_speed_facial = std::max(max_char_speed_facial, ws[f]);
651 }
652
653 return max_char_speed_facial;
654 }
655
656
657 template<typename PhysicsT>
658 mfem::real_t NSOperator<PhysicsT>::MultCNS_BoundaryFaces(const mfem::Vector &pu,
659 const std::vector<mfem::Vector *> &p_grad_prim,
660 mfem::Vector &pdudt) const
661 {
662 Theseus::ScopedTimer timer("MultCNS_BoundaryFaces");
663
664 auto dc = device_cache;
665 const int dim = dc.dim;
666 const int neq = dc.num_equations;
667 const int nfp = dc.num_face_points;
668 const int face_size = nfp * neq;
669 const int restr_size = operator_cache.restr_b->Height();
670 const int nfaces_restr = restr_size / face_size;
671 const int norm_size = nfp * dc.dim;
672 const int npoints_bnd = nfaces_restr * nfp;
673 const int psize = pdudt.Size();
674
675 if(restr_size == 0){
676 return 0.0;
677 }
678
679 mfem::Vector &rhs_faces(operator_cache.rhsBnd);
680 mfem::Vector &faces_dudt(operator_cache.dudtBnd);
681 if(rhs_faces.Size() != restr_size){
682 rhs_faces.SetSize(restr_size);
683 rhs_faces.UseDevice();
684 }
685 if(faces_dudt.Size() != psize){
686 faces_dudt.SetSize(psize);
687 faces_dudt.UseDevice();
688 }
689 mfem::Vector &bnd_u(operator_cache.uBnd);
690 if(bnd_u.Size() != restr_size){
691 bnd_u.SetSize(restr_size);
692 bnd_u.UseDevice();
693 }
694 const mfem::real_t *grad_prim_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
695 std::vector<mfem::Vector> &bnd_grad_prim(operator_cache.gradBnd);
696 if(bnd_grad_prim.size() != dim){
697 bnd_grad_prim.resize(dim);
698 for(int idim = 0;idim < dim;idim++){
699 bnd_grad_prim[idim].SetSize(restr_size);
700 bnd_grad_prim[idim].UseDevice();
701 }
702 }
703 for(int idim = 0;idim < dim;idim++){
704 operator_cache.restr_b->Mult(*p_grad_prim[idim], bnd_grad_prim[idim]);
705 grad_prim_d[idim] = bnd_grad_prim[idim].Read();
706 }
707
708 // If zeroed before accumulation, do it explicitly on device:
709 // Potentially, this is not needed at all since I think we overwrite everything
710 {
711 mfem::real_t *rd = rhs_faces.Write();
712 mfem::forall(rhs_faces.Size(), [=] MFEM_HOST_DEVICE (int i)
713 { rd[i] = mfem::real_t(0);});
714 mfem::real_t *fd = faces_dudt.Write();
715 mfem::forall(faces_dudt.Size(), [=] MFEM_HOST_DEVICE (int i)
716 { fd[i] = mfem::real_t(0);});
717 }
718
719 if(!operator_cache.u_bnd_restr_ready){
720 operator_cache.restr_b->Mult(pu, bnd_u);
721 operator_cache.u_bnd_restr_ready = true;
722 }
723
724 const mfem::real_t *u_d = bnd_u.Read();
725 mfem::real_t *rhs_d = rhs_faces.Write();
726
727 const mfem::real_t *nor_d = dc.bnd_nor_d; // size nfaces*nfp*dim
728 const mfem::real_t *radius_d = dc.bnd_radius_d;
729 const mfem::real_t *inv1_d = dc.bnd_wt_d; // size nfaces*nfp
730 const int *bnd_marker_index_d = dc.bnd_marker_index_d;
731 mfem::real_t *ws_d = dc.bndWaveSpeed_d;
732
733 mfem::forall(npoints_bnd, [=] MFEM_HOST_DEVICE (int p)
734 {
735 const int f = p / nfp;
736 const int fp = p % nfp;
737
738 int bnd_face_marker_index = bnd_marker_index_d[f];
739 if(bnd_face_marker_index < 0){
740 ws_d[p] = 0.0;
741 return;
742 }
743 int bc_index = bnd_face_marker_index; // no mapping atm
744 if(bc_index < 0){
745 ws_d[p] = 0.0;
746 return;
747 }
748 const Theseus::BCDescriptor &bc = dc.bc_descr_d[bc_index];
749 if (bc.type == int(Theseus::BCType::Invalid))
750 {
751 ws_d[p] = 0.0;
752 return;
753 }
754
755 const int face_offset = f * face_size;
756 const int n_offset = f * norm_size;
757 const int w_offset = f * nfp;
758
759 const mfem::real_t *u_face_d = u_d + face_offset;
760 mfem::real_t *rhs_face_d = rhs_d + face_offset;
761 const mfem::real_t *nor_face_d = nor_d + n_offset;
762 const mfem::real_t *w_minus_d = inv1_d + w_offset;
763 const mfem::real_t *nor_point = nor_face_d + fp*dim;
764 const mfem::real_t radius = dc.axisymmetric ? radius_d[w_offset + fp] : 0.0;
765 mfem::real_t scale = -w_minus_d[fp];
766 mfem::real_t state1[Theseus::MAXEQ];
767 mfem::real_t fluxN[Theseus::MAXEQ];
768 mfem::real_t gradPrim_x[Theseus::MAXEQ];
769 mfem::real_t gradPrim_y[Theseus::MAXEQ];
770 mfem::real_t gradPrim_z[Theseus::MAXEQ];
771 const mfem::real_t *dprim_face_x = (dim > 0) ? grad_prim_d[0] + face_offset : nullptr;
772 const mfem::real_t *dprim_face_y = (dim > 1) ? grad_prim_d[1] + face_offset : nullptr;
773 const mfem::real_t *dprim_face_z = (dim > 2) ? grad_prim_d[2] + face_offset : nullptr;
774 Theseus::Kernels::el_gather_grad_state(dprim_face_x, dprim_face_y, dprim_face_z,
775 dim, nfp, neq, fp, gradPrim_x, gradPrim_y,
776 gradPrim_z);
777 Theseus::Kernels::el_gather_state(u_face_d, nfp, neq, fp, state1);
778
779 const mfem::real_t ws = \
780 Theseus::BC::ApplyViscousBoundaryCondition(dc, bc, state1, gradPrim_x, gradPrim_y,
781 gradPrim_z, nor_point, radius, fluxN);
782 Theseus::Kernels::el_scatter_add(fluxN, nfp, neq, fp, scale, rhs_face_d);
783 ws_d[p] = ws;
784 });
785
786 operator_cache.restr_b->MultTranspose(rhs_faces, faces_dudt);
787 pdudt += faces_dudt; // on device? (likely yes)
788
789 // Finish up on the host:
790 // - Reduce for rank-local max_char_speed
791 const mfem::real_t *ws = operator_cache.bndWaveSpeed.HostRead();
792 mfem::real_t max_char_speed_facial = 0.0;
793 for(int p = 0;p < npoints_bnd;p++)
794 {
795 max_char_speed_facial = std::max(max_char_speed_facial, ws[p]);
796 }
797
798 return max_char_speed_facial;
799 }
800
801 template<typename PhysicsT>
802 mfem::real_t NSOperator<PhysicsT>::MultCNS_Volume(const mfem::Vector &pu, const std::vector<mfem::Vector *> &p_grad_prim,
803 mfem::Vector &pdudt) const
804 {
805 Theseus::ScopedTimer timer("MultCNS_Volume");
806 // Copy the device cache so that it is not member data
807 auto dc = device_cache;
808 const int dim = dc.dim;
809 const int restr_size = operator_cache.restr_v->Height();
810
811 mfem::Vector &vol_u(operator_cache.uVol);
812 if(vol_u.Size() != restr_size){
813 vol_u.SetSize(restr_size);
814 vol_u.UseDevice();
815 }
816
817 mfem::Vector &dUe(operator_cache.rhsVol);
818 if(dUe.Size() != restr_size){
819 dUe.SetSize(restr_size);
820 dUe.UseDevice();
821 }
822
823 std::vector<mfem::Vector> &vol_grad_prim(operator_cache.gradVol);
824 if(vol_grad_prim.size() != dim){
825 vol_grad_prim.resize(dim);
826 for(int idim = 0;idim < dim;idim++){
827 vol_grad_prim[idim].SetSize(restr_size);
828 vol_grad_prim[idim].UseDevice();
829 }
830 }
831
832 if(!operator_cache.u_vol_restr_ready){
833 operator_cache.restr_v->Mult(pu, vol_u);
834 operator_cache.u_vol_restr_ready = true;
835 }
836 for(int idim = 0;idim < dim;idim++){
837 operator_cache.restr_v->Mult(*p_grad_prim[idim], vol_grad_prim[idim]);
838 }
839
840 // Zero the RHS array on-device
841 {
842 mfem::real_t *d = dUe.Write();
843 mfem::forall(dUe.Size(), [=] MFEM_HOST_DEVICE (int i) { d[i] = mfem::real_t(0); });
844 }
845
846 // Set up the read-only pointers for restr inputs
847 const mfem::real_t *Ue_d = vol_u.Read();
848 const mfem::real_t *gradPrim_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
849 for(int idim = 0;idim < dim;idim++){
850 gradPrim_d[idim] = vol_grad_prim[idim].Read();
851 }
852
853 // Write-only for RHS
854 mfem::real_t *dUe_d = dUe.Write();
855
856 // Device cache parameters
857 const int ne = dc.num_elements;
858 const int ndof = dc.ndof_scalar_el;
859 const int neq = dc.num_equations;
860
861#ifdef SUBCELL_FV_BLENDING
862 const int Np_x = dc.Np_x;
863 const int Np_y = dc.Np_y;
864 const int Np_z = dc.Np_z;
865 const int npe_metric_xi = (Np_x + 1)*Np_y*Np_z;
866 const int npe_metric_eta = Np_x*(Np_y + 1)*Np_z;
867 const int npe_metric_zeta = Np_x * Np_y * (Np_z + 1);
868 const mfem::real_t *metric_xi_d = dc.subcell_metric_xi_d;
869 const mfem::real_t *metric_eta_d = (dim > 1 ? dc.subcell_metric_eta_d : nullptr);
870 const mfem::real_t *metric_zeta_d = (dim > 2 ? dc.subcell_metric_zeta_d : nullptr);
871
872 mfem::Vector dUfv(operator_cache.restr_v->Height());
873 dUfv.UseDevice();
874 mfem::real_t *dUfv_d = dUfv.Write();
875 // zero the array on-device
876 {
877 mfem::real_t *d = dUfv_d;
878 mfem::forall(dUfv.Size(), [=] MFEM_HOST_DEVICE (int i) { d[i] = mfem::real_t(0); });
879 }
880
881 const mfem::real_t *alpha_d = operator_cache.alpha->Read();
882#endif
883
884 // Derived parameters
885 const int metric_stride = ndof * dim * dim;
886 const int jac_stride = ndof;
887 const int estride = ndof*neq;
888
889 // Device cache data/arrays
890 const int *elem_attr_d = dc.elem_attr_d;
891 const int *attr_marker_d = dc.attr_marker_d;
892 const mfem::real_t *elJac_d = dc.elJac_d;
893 const mfem::real_t *elMetric_d = dc.elMetric_d;
894 const mfem::real_t *elRadius_d = dc.elRadius_d;
895
896 mfem::real_t *ws_d = dc.elWaveSpeed_d;
897
898 // Inside the FORALL below, executed on device
899 mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
900 {
901
902 const mfem::real_t *jac_el = elJac_d + e * jac_stride;
903 const mfem::real_t *metric_el = elMetric_d + e * metric_stride;
904 const mfem::real_t *radius_el = dc.axisymmetric ?
905 elRadius_d + e * jac_stride : nullptr;
906
907 const int attr = elem_attr_d[e];
908 if (attr_marker_d[attr-1] == 0) {
909 ws_d[e] = 0.0;
910 return;
911 }
912
913 // Element-specific inputs and outputs
914 const int eoff = e * estride;
915 const mfem::real_t *u_el = Ue_d + eoff;
916 mfem::real_t *du_el = dUe_d + eoff;
917
918 mfem::real_t cs_el = \
919 Theseus::DGSEMIntegrator::AssembleElementVolumeKernel(dc, u_el,
920 jac_el, metric_el, du_el);
921#ifdef SUBCELL_FV_BLENDING
922 mfem::real_t alpha_fv = alpha_d[e];
923 if(alpha_fv > 1e-16){
924 mfem::real_t alpha_dg = (1.0 - alpha_fv);
925 mfem::real_t *du_fv = dUfv_d + eoff;
926 const mfem::real_t *el_metric_xi = metric_xi_d + e * npe_metric_xi * dim;
927 const mfem::real_t *el_metric_eta = (dim > 1 ? metric_eta_d + e * npe_metric_eta * dim :
928 nullptr);
929 const mfem::real_t *el_metric_zeta = (dim > 2 ? metric_zeta_d + e * npe_metric_zeta * dim :
930 nullptr);
931 const mfem::real_t cs_fv = \
932 Theseus::DGSEMIntegrator::ComputeFVFluxesKernel(dc, u_el, jac_el, el_metric_xi, el_metric_eta, el_metric_zeta, du_fv);
933
934 for(int ipt = 0;ipt < estride;ipt++){
935 du_el[ipt] = alpha_dg * du_el[ipt] + alpha_fv * du_fv[ipt];
936 }
937
938 cs_el = Kernels::rmax(cs_el, cs_fv);
939 }
940#endif
941 AddAxisymmetricEulerElementSource(
942 dc, u_el, radius_el, jac_el, metric_el, du_el);
943 ws_d[e] = cs_el;
944
945 // Inviscid part is done: dUe currrently holds the inviscid RHS
946 // Host code mixes inviscid and viscous assembly, we need separate.
947 // Call the Viscous Assembly routine
948 const mfem::real_t *grad_prim_el[Theseus::MAXDIM] = {nullptr, nullptr, nullptr};
949 for(int idim = 0;idim < dim;idim++){
950 grad_prim_el[idim] = gradPrim_d[idim] + eoff;
951 }
952
953 Theseus::DGSEMIntegrator::AssembleViscousElementVolumeKernel(dc, u_el, jac_el, metric_el,
954 radius_el,
955 grad_prim_el[0], grad_prim_el[1],
956 grad_prim_el[2], du_el);
957
958 });
959
960 // The rest is identical to Euler operator
961 // Scatter RHS back to storage
962 operator_cache.restr_v->AddMultTranspose(dUe, pdudt);
963
964 // Finish up on the host:
965 // - Reduce for rank-local max_char_speed
966 const mfem::real_t *ws = operator_cache.elWaveSpeed.HostRead();
967 mfem::real_t max_char_speed = 0.0;
968 for(int e = 0;e < operator_cache.num_elements;e++)
969 {
970 max_char_speed = std::max(max_char_speed, ws[e]);
971 }
972
973 return max_char_speed;
974 }
975
976 template<typename PhysicsT>
977 mfem::real_t NSOperator<PhysicsT>::MultCNS(const mfem::Vector &u, const std::vector<mfem::Vector *> &grad_prim,
978 mfem::Vector &pdudt) const
979 {
980 const int dim = operator_cache.dim;
981 std::vector<mfem::Vector *> p_grad_(dim);
982 if (this->P)
983 {
984 const int psize = this->P->Height();
985 if(operator_cache.pGrad.size() != dim){
986 operator_cache.pGrad.resize(dim);
987 for(int idim = 0;idim < dim;idim++){
988 operator_cache.pGrad[idim].SetSize(psize);
989 operator_cache.pGrad[idim].UseDevice();
990 }
991 }
992 for(int idim = 0;idim < dim;idim++){
993 p_grad_[idim] = &(operator_cache.pGrad[idim]);
994 this->P->Mult(*grad_prim[idim], *p_grad_[idim]);
995 }
996 if(operator_cache.pdudt.Size() != psize){
997 operator_cache.pdudt.SetSize(psize);
998 operator_cache.pdudt.UseDevice();
999 }
1000 }
1001 const std::vector<mfem::Vector *> &pGradPrim = this->P ? p_grad_ : grad_prim;
1002
1003 mfem::real_t max_char_speed = MultCNS_Volume(u, pGradPrim, pdudt);
1004
1005 mfem::real_t max_char_speed_faces = MultCNS_InteriorFaces(u, pGradPrim, pdudt);
1006 max_char_speed = std::max(max_char_speed, max_char_speed_faces);
1007
1008 mfem::real_t max_char_speed_bnd = MultCNS_BoundaryFaces(u, pGradPrim, pdudt);
1009 max_char_speed = std::max(max_char_speed, max_char_speed_bnd);
1010
1011 return max_char_speed;
1012 }
1013
1014}
Theseus::NSOperator::ComputeEntropyState
void ComputeEntropyState(const mfem::Vector &u, mfem::Vector &e) const
Definition NSOperator_impl.hpp:12
Theseus::NSOperator::MultCNS
mfem::real_t MultCNS(const mfem::Vector &u, const std::vector< mfem::Vector * > &grad_prim, mfem::Vector &pdudt) const
Definition NSOperator_impl.hpp:977
Theseus::NSOperator::GradOperator
void GradOperator(const mfem::Vector &u, std::vector< mfem::Vector * > &grad_u) const
Definition NSOperator_impl.hpp:482
Theseus::NSOperator::MultCNS_InteriorFaces
mfem::real_t MultCNS_InteriorFaces(const mfem::Vector &pu, const std::vector< mfem::Vector * > &p_grad_prim, mfem::Vector &pdudt) const
Definition NSOperator_impl.hpp:542
Theseus::NSOperator::MultCNS_Volume
mfem::real_t MultCNS_Volume(const mfem::Vector &pu, const std::vector< mfem::Vector * > &p_grad_prim, mfem::Vector &pdudt) const
Definition NSOperator_impl.hpp:802
Theseus::NSOperator::GradOperator_Volume
void GradOperator_Volume(const mfem::Vector &pu, std::vector< mfem::Vector * > &p_grad_u) const
Definition NSOperator_impl.hpp:196
Theseus::NSOperator::GradOperator_InteriorFaces
void GradOperator_InteriorFaces(const mfem::Vector &pu, std::vector< mfem::Vector * > &p_grad_u) const
Definition NSOperator_impl.hpp:389
Theseus::NSOperator::FlowMult
mfem::real_t FlowMult(const mfem::Vector &u, mfem::Vector &dudt) const override
Definition NSOperator_impl.hpp:162
Theseus::NSOperator::ComputeGradPrimFromGradEntropy
void ComputeGradPrimFromGradEntropy(const mfem::Vector &u, std::vector< mfem::Vector * > &gradEntropy) const
Definition NSOperator_impl.hpp:83
Theseus::NSOperator::GradOperator_BoundaryFaces
void GradOperator_BoundaryFaces(const mfem::Vector &pu, std::vector< mfem::Vector * > &p_grad_u) const
Definition NSOperator_impl.hpp:269
Theseus::NSOperator::MultCNS_BoundaryFaces
mfem::real_t MultCNS_BoundaryFaces(const mfem::Vector &pu, const std::vector< mfem::Vector * > &p_grad_prim, mfem::Vector &pdudt) const
Definition NSOperator_impl.hpp:658
Theseus::ScopedTimer
Definition timer.hpp:17
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::AddAxisymmetricEulerElementSource
MFEM_HOST_DEVICE void AddAxisymmetricEulerElementSource(const ContextT &ctx, const mfem::real_t *element_state, const mfem::real_t *element_radius, const mfem::real_t *element_jacobian, const mfem::real_t *element_metric, mfem::real_t *element_rate)
Definition AxisymmetricSource.hpp:221
Theseus::MAXDIM
constexpr const int MAXDIM
Definition theseus_kernels.hpp:14
Theseus::MAXEQ
constexpr const int MAXEQ
Definition theseus_kernels.hpp:13
Theseus::ProjectAxisPrimitiveGradientDirection
MFEM_HOST_DEVICE bool ProjectAxisPrimitiveGradientDirection(const GasT &gas, mfem::real_t radius, int derivative_direction, mfem::real_t *dprim)
Definition AxisymmetricSource.hpp:261
Theseus::BCDescriptor
Definition bc_cache_utilities.hpp:35
Theseus::BCDescriptor::type
int type
Definition bc_cache_utilities.hpp:36
Theseus::PointStateViewRW
Definition GasState.hpp:218
Theseus::PointStateView
Definition GasState.hpp:127