1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
4 // SPDX-License-Identifier: BSD-2-Clause
9 /// Operator for Navier-Stokes example using PETSc
41     State      s    = StateFromPrimitive(&context->gas, context->reference);  in ICsNewtonianIG()
42     StateToQ(&context->gas, s, q, state_var);  in ICsNewtonianIG()
43     for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];  in ICsNewtonianIG()  local
87     Tau_diagPrim(context, s, dXdx, context->dt, Tau_d);  in MassFunction_Newtonian()
91     for (CeedInt j = 0; j < 5; j++) {  in MassFunction_Newtonian()  local
93 …Grad_v[k][j][i] = wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k…  in MassFunction_Newtonian()
105 // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping…
107 // This is 3D compressible Navier-Stokes in conservation form with state variables of density, mome…
110 //   rho - Mass Density
111 //   Ui  - Momentum Density,      Ui = rho ui
112 //   E   - Total Energy Density,  E  = rho (cv T + (u u)/2 + g z)
114 // Navier-Stokes Equations:
125 //   P = (gamma - 1) (E - rho (u u) / 2 - rho g z)
135 //  SU   = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) )
138 //   lambda = - 2 / 3,  From Stokes hypothesis
146 // We require the product of the inverse of the Jacobian (dXdx_j,k) and its transpose (dXdx_k,j) to…
158   const CeedScalar        *g       = context->g;  in RHSFunction_Newtonian()
159   const CeedScalar         dt      = context->dt;  in RHSFunction_Newtonian()
160   const CeedScalar         P0      = context->idl_pressure;  in RHSFunction_Newtonian()
165     for (int j = 0; j < 5; j++) U[j] = q[j][i];  in RHSFunction_Newtonian()  local
185     for (CeedInt j = 0; j < 5; j++) {  in RHSFunction_Newtonian()  local
186 …nt k = 0; k < 3; k++) Grad_v[k][j][i] = wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1]…  in RHSFunction_Newtonian()
190     for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j];  in RHSFunction_Newtonian()  local
192     if (context->idl_enable) {  in RHSFunction_Newtonian()
193 …dScalar sigma         = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context…  in RHSFunction_Newtonian()
194       CeedScalar       damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.};  in RHSFunction_Newtonian()
196       for (int j = 0; j < 5; j++) v[j][i] -= wdetJ * idl_residual[j];  in RHSFunction_Newtonian()  local
199     // -- Stabilization method: none (Galerkin), SU, or SUPG  in RHSFunction_Newtonian()
204     for (CeedInt j = 0; j < 5; j++) {  in RHSFunction_Newtonian()  local
205 …(CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] -= wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXd…  in RHSFunction_Newtonian()
212 // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepp…
214 //  SU   = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) )
215 //  SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) )
229   const CeedScalar        *g       = context->g;  in IFunction_Newtonian()
230   const CeedScalar         dt      = context->dt;  in IFunction_Newtonian()
231   const CeedScalar         P0      = context->idl_pressure;  in IFunction_Newtonian()
256     for (CeedInt j = 0; j < 5; j++) {  in IFunction_Newtonian()  local
258 …Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[…  in IFunction_Newtonian()
264     // -- Stabilization method: none (Galerkin), SU, or SUPG  in IFunction_Newtonian()
266     for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i];  in IFunction_Newtonian()  local
270     for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]);  in IFunction_Newtonian()  local
271     if (context->idl_enable) {  in IFunction_Newtonian()
272 …onst CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context…  in IFunction_Newtonian()
274       CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.};  in IFunction_Newtonian()
276       for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j];  in IFunction_Newtonian()  local
282     for (CeedInt j = 0; j < 5; j++) {  in IFunction_Newtonian()  local
284 …Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[…  in IFunction_Newtonian()
307 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping…
318   const CeedScalar        *g       = context->g;  in IJacobian_Newtonian()
331     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];  in IJacobian_Newtonian()  local
351     for (int j = 0; j < 5; j++) {  in IJacobian_Newtonian()  local
352 … k = 0; k < 3; k++) Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * dFlux[j][0] + dXdx[k][1] * dFlux[j][1…  in IJacobian_Newtonian()
358 …for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force…  in IJacobian_Newtonian()  local
360     if (context->idl_enable) {  in IJacobian_Newtonian()
363 …// This is a Picard-type linearization of the damping and could be replaced by an InternalDampingL…  in IJacobian_Newtonian()
365       for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j];  in IJacobian_Newtonian()  local
368     // -- Stabilization method: none (Galerkin), SU, or SUPG  in IJacobian_Newtonian()
370     for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j];  in IJacobian_Newtonian()  local
373     for (int j = 0; j < 5; j++) {  in IJacobian_Newtonian()  local
374 …or (int k = 0; k < 3; k++) Grad_v[k][j][i] += wdetJ * (dstab[j][0] * dXdx[k][0] + dstab[j][1] * dX…  in IJacobian_Newtonian()
401   CeedScalar(*jac_data_sur)              = context->is_implicit ? out[1] : NULL;  in BoundaryIntegral()
403   const bool is_implicit = context->is_implicit;  in BoundaryIntegral()
411     wdetJb *= is_implicit ? -1. : 1.;  in BoundaryIntegral()
428     for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];  in BoundaryIntegral()  local
462   const bool                     is_implicit = context->is_implicit;  in BoundaryIntegral_Jacobian()
467     wdetJb *= is_implicit ? -1. : 1.;  in BoundaryIntegral_Jacobian()
472     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];  in BoundaryIntegral_Jacobian()  local
493     for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j];  in BoundaryIntegral_Jacobian()  local