1dc936754SJeremy L Thompson // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2727da7e7SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 33a8779fbSJames Wright // 4727da7e7SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 53a8779fbSJames Wright // 6727da7e7SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 73a8779fbSJames Wright 83a8779fbSJames Wright /// @file 93a8779fbSJames Wright /// Operator for Navier-Stokes example using PETSc 103a8779fbSJames Wright #include <ceed.h> 11d0cce58aSJeremy L Thompson #include <math.h> 127b530f2aSAdelekeBankole #include <stdlib.h> 132b916ea7SJeremy L Thompson 14475b2820SJames Wright #include "newtonian_state.h" 15d0cce58aSJeremy L Thompson #include "newtonian_types.h" 16d1b9ef12SLeila Ghaffari #include "stabilization.h" 17d0cce58aSJeremy L Thompson #include "utils.h" 18bb8a0c61SJames Wright 1994a7b3d2SKenneth E. Jansen CEED_QFUNCTION_HELPER void InternalDampingLayer(const NewtonianIdealGasContext context, const State s, const CeedScalar sigma, CeedScalar damp_Y[5], 20e7754af5SKenneth E. Jansen CeedScalar damp_residual[5]) { 21e7754af5SKenneth E. Jansen ScaleN(damp_Y, sigma, 5); 22edcfef1bSKenneth E. Jansen State damp_s = StateFromY_fwd(context, s, damp_Y); 23e7754af5SKenneth E. Jansen 24e7754af5SKenneth E. Jansen CeedScalar U[5]; 25e7754af5SKenneth E. Jansen UnpackState_U(damp_s.U, U); 26e7754af5SKenneth E. Jansen for (int i = 0; i < 5; i++) damp_residual[i] += U[i]; 27e7754af5SKenneth E. Jansen } 28e7754af5SKenneth E. Jansen 29bb8a0c61SJames Wright // ***************************************************************************** 303a8779fbSJames Wright // This QFunction sets a "still" initial condition for generic Newtonian IG problems 313a8779fbSJames Wright // ***************************************************************************** 328fff8293SJames Wright CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 333a8779fbSJames Wright // Inputs 343a8779fbSJames Wright 353a8779fbSJames Wright // Outputs 363a8779fbSJames Wright CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 373a8779fbSJames Wright 38bb8a0c61SJames Wright // Context 39bb8a0c61SJames Wright const SetupContext context = (SetupContext)ctx; 40bb8a0c61SJames Wright 413a8779fbSJames Wright // Quadrature Point Loop 422b916ea7SJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 433a8779fbSJames Wright CeedScalar q[5] = {0.}; 44edcfef1bSKenneth E. Jansen State s = StateFromPrimitive(&context->gas, context->reference); 458fff8293SJames Wright StateToQ(&context->gas, s, q, state_var); 462b916ea7SJeremy L Thompson for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 473a8779fbSJames Wright } // End of Quadrature Point Loop 483a8779fbSJames Wright return 0; 493a8779fbSJames Wright } 503a8779fbSJames Wright 512b916ea7SJeremy L Thompson CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 528fff8293SJames Wright return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_PRIMITIVE); 53b8fb7609SAdeleke O. Bankole } 54b8fb7609SAdeleke O. Bankole CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 558fff8293SJames Wright return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 56cbe60e31SLeila Ghaffari } 57cbe60e31SLeila Ghaffari 58cbe60e31SLeila Ghaffari // ***************************************************************************** 5904e40bb6SJeremy L Thompson // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method 603a8779fbSJames Wright // 6104e40bb6SJeremy L Thompson // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density. 623a8779fbSJames Wright // 633a8779fbSJames Wright // State Variables: q = ( rho, U1, U2, U3, E ) 643a8779fbSJames Wright // rho - Mass Density 653a8779fbSJames Wright // Ui - Momentum Density, Ui = rho ui 663a8779fbSJames Wright // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 673a8779fbSJames Wright // 683a8779fbSJames Wright // Navier-Stokes Equations: 693a8779fbSJames Wright // drho/dt + div( U ) = 0 703a8779fbSJames Wright // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 713a8779fbSJames Wright // dE/dt + div( (E + P) u ) = div( Fe ) 723a8779fbSJames Wright // 733a8779fbSJames Wright // Viscous Stress: 743a8779fbSJames Wright // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 753a8779fbSJames Wright // 763a8779fbSJames Wright // Thermal Stress: 773a8779fbSJames Wright // Fe = u Fu + k grad( T ) 78bb8a0c61SJames Wright // Equation of State 793a8779fbSJames Wright // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 803a8779fbSJames Wright // 813a8779fbSJames Wright // Stabilization: 823a8779fbSJames Wright // Tau = diag(TauC, TauM, TauM, TauM, TauE) 833a8779fbSJames Wright // f1 = rho sqrt(ui uj gij) 843a8779fbSJames Wright // gij = dXi/dX * dXi/dX 853a8779fbSJames Wright // TauC = Cc f1 / (8 gii) 863a8779fbSJames Wright // TauM = min( 1 , 1 / f1 ) 873a8779fbSJames Wright // TauE = TauM / (Ce cv) 883a8779fbSJames Wright // 893a8779fbSJames Wright // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 903a8779fbSJames Wright // 913a8779fbSJames Wright // Constants: 923a8779fbSJames Wright // lambda = - 2 / 3, From Stokes hypothesis 933a8779fbSJames Wright // mu , Dynamic viscosity 943a8779fbSJames Wright // k , Thermal conductivity 953a8779fbSJames Wright // cv , Specific heat, constant volume 963a8779fbSJames Wright // cp , Specific heat, constant pressure 973a8779fbSJames Wright // g , Gravity 983a8779fbSJames Wright // gamma = cp / cv, Specific heat ratio 993a8779fbSJames Wright // 10004e40bb6SJeremy L Thompson // We require the product of the inverse of the Jacobian (dXdx_j,k) and its transpose (dXdx_k,j) to properly compute integrals of the form: int( gradv 10104e40bb6SJeremy L Thompson // gradu ) 1023a8779fbSJames Wright // ***************************************************************************** 1032b916ea7SJeremy L Thompson CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 1043a8779fbSJames Wright // Inputs 1053d65b166SJames Wright const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 10687bd45e7SJames Wright const CeedScalar(*Grad_q) = in[1]; 107ade49511SJames Wright const CeedScalar(*q_data) = in[2]; 1083d65b166SJames Wright 1093a8779fbSJames Wright // Outputs 1103d65b166SJames Wright CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 1113d65b166SJames Wright CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 1123a8779fbSJames Wright 1133a8779fbSJames Wright // Context 1143a8779fbSJames Wright NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 115bb8a0c61SJames Wright const CeedScalar *g = context->g; 116bb8a0c61SJames Wright const CeedScalar dt = context->dt; 1173a8779fbSJames Wright 1183a8779fbSJames Wright // Quadrature Point Loop 1193d65b166SJames Wright CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 120ade49511SJames Wright CeedScalar U[5], wdetJ, dXdx[3][3]; 121c1a52365SJed Brown for (int j = 0; j < 5; j++) U[j] = q[j][i]; 122ade49511SJames Wright StoredValuesUnpack(Q, i, 0, 1, q_data, &wdetJ); 123ade49511SJames Wright StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx); 124edcfef1bSKenneth E. Jansen State s = StateFromU(context, U); 125c1a52365SJed Brown 126c1a52365SJed Brown State grad_s[3]; 127edcfef1bSKenneth E. Jansen StatePhysicalGradientFromReference(Q, i, context, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s); 128c1a52365SJed Brown 129c1a52365SJed Brown CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 13040a33f2dSJames Wright KMStrainRate_State(grad_s, strain_rate); 131c1a52365SJed Brown NewtonianStress(context, strain_rate, kmstress); 132c1a52365SJed Brown KMUnpack(kmstress, stress); 133c1a52365SJed Brown ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 134c1a52365SJed Brown 135c1a52365SJed Brown StateConservative F_inviscid[3]; 136c1a52365SJed Brown FluxInviscid(context, s, F_inviscid); 137c1a52365SJed Brown 138c1a52365SJed Brown // Total flux 139c1a52365SJed Brown CeedScalar Flux[5][3]; 140d1b9ef12SLeila Ghaffari FluxTotal(F_inviscid, stress, Fe, Flux); 141c1a52365SJed Brown 1427523f6aaSJames Wright for (CeedInt j = 0; j < 5; j++) { 1437523f6aaSJames Wright for (CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] = wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 1442b916ea7SJeremy L Thompson } 145c1a52365SJed Brown 14660dbb574SKenneth E. Jansen const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], Dot3(s.U.momentum, g)}; 1472b916ea7SJeremy L Thompson for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 1483a8779fbSJames Wright 149d1b9ef12SLeila Ghaffari // -- Stabilization method: none (Galerkin), SU, or SUPG 150d1b9ef12SLeila Ghaffari CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 151d1b9ef12SLeila Ghaffari Tau_diagPrim(context, s, dXdx, dt, Tau_d); 152edcfef1bSKenneth E. Jansen Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 1533a8779fbSJames Wright 1542b916ea7SJeremy L Thompson for (CeedInt j = 0; j < 5; j++) { 1552b916ea7SJeremy L Thompson for (CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] -= wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 1562b916ea7SJeremy L Thompson } 1573a8779fbSJames Wright } // End Quadrature Point Loop 1583a8779fbSJames Wright 1593a8779fbSJames Wright // Return 1603a8779fbSJames Wright return 0; 1613a8779fbSJames Wright } 1623a8779fbSJames Wright 1633a8779fbSJames Wright // ***************************************************************************** 16404e40bb6SJeremy L Thompson // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method 1653a8779fbSJames Wright // 1663a8779fbSJames Wright // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 1673a8779fbSJames Wright // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 16804e40bb6SJeremy L Thompson // (diffusive terms will be added later) 1693a8779fbSJames Wright // ***************************************************************************** 1708fff8293SJames Wright CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 1713a8779fbSJames Wright // Inputs 1723d65b166SJames Wright const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 17387bd45e7SJames Wright const CeedScalar(*Grad_q) = in[1]; 1743d65b166SJames Wright const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 175ade49511SJames Wright const CeedScalar(*q_data) = in[3]; 1763d65b166SJames Wright const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1773d65b166SJames Wright 1783a8779fbSJames Wright // Outputs 1793d65b166SJames Wright CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 1803d65b166SJames Wright CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 181ade49511SJames Wright CeedScalar(*jac_data) = out[2]; 1823d65b166SJames Wright 1833a8779fbSJames Wright // Context 1843a8779fbSJames Wright NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 185bb8a0c61SJames Wright const CeedScalar *g = context->g; 186bb8a0c61SJames Wright const CeedScalar dt = context->dt; 187e7754af5SKenneth E. Jansen const CeedScalar P0 = context->P0; 1883a8779fbSJames Wright 1893a8779fbSJames Wright // Quadrature Point Loop 1903d65b166SJames Wright CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 1913d65b166SJames Wright const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 192c1a52365SJed Brown const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 193edcfef1bSKenneth E. Jansen const State s = StateFromQ(context, qi, state_var); 194c1a52365SJed Brown 195ade49511SJames Wright CeedScalar wdetJ, dXdx[3][3]; 196ade49511SJames Wright QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 197c1a52365SJed Brown State grad_s[3]; 198edcfef1bSKenneth E. Jansen StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 199c1a52365SJed Brown 200c1a52365SJed Brown CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 20140a33f2dSJames Wright KMStrainRate_State(grad_s, strain_rate); 202c1a52365SJed Brown NewtonianStress(context, strain_rate, kmstress); 203c1a52365SJed Brown KMUnpack(kmstress, stress); 204c1a52365SJed Brown ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 205c1a52365SJed Brown 206c1a52365SJed Brown StateConservative F_inviscid[3]; 207c1a52365SJed Brown FluxInviscid(context, s, F_inviscid); 208c1a52365SJed Brown 209c1a52365SJed Brown // Total flux 210c1a52365SJed Brown CeedScalar Flux[5][3]; 211d1b9ef12SLeila Ghaffari FluxTotal(F_inviscid, stress, Fe, Flux); 212c1a52365SJed Brown 2137523f6aaSJames Wright for (CeedInt j = 0; j < 5; j++) { 2147523f6aaSJames Wright for (CeedInt k = 0; k < 3; k++) { 2157523f6aaSJames Wright Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 2163d65b166SJames Wright } 2172b916ea7SJeremy L Thompson } 218c1a52365SJed Brown 21960dbb574SKenneth E. Jansen const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], Dot3(s.U.momentum, g)}; 2203a8779fbSJames Wright 221d1b9ef12SLeila Ghaffari // -- Stabilization method: none (Galerkin), SU, or SUPG 222edcfef1bSKenneth E. Jansen CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5]; 22376555becSJames Wright for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 224edcfef1bSKenneth E. Jansen State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var); 22576555becSJames Wright UnpackState_U(s_dot.U, U_dot); 22676555becSJames Wright 2272b916ea7SJeremy L Thompson for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 228e7754af5SKenneth E. Jansen if (context->idl_enable) { 22994a7b3d2SKenneth E. Jansen const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 23094a7b3d2SKenneth E. Jansen StoredValuesPack(Q, i, 14, 1, &sigma, jac_data); 231e7754af5SKenneth E. Jansen CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 23294a7b3d2SKenneth E. Jansen InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 233e7754af5SKenneth E. Jansen for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 234e7754af5SKenneth E. Jansen } 235e7754af5SKenneth E. Jansen 236d1b9ef12SLeila Ghaffari Tau_diagPrim(context, s, dXdx, dt, Tau_d); 237edcfef1bSKenneth E. Jansen Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 2383a8779fbSJames Wright 2392b916ea7SJeremy L Thompson for (CeedInt j = 0; j < 5; j++) { 2403d65b166SJames Wright for (CeedInt k = 0; k < 3; k++) { 2413d65b166SJames Wright Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 2423d65b166SJames Wright } 2432b916ea7SJeremy L Thompson } 244ade49511SJames Wright StoredValuesPack(Q, i, 0, 5, qi, jac_data); 245ade49511SJames Wright StoredValuesPack(Q, i, 5, 6, kmstress, jac_data); 246ade49511SJames Wright StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data); 2473a8779fbSJames Wright 2483a8779fbSJames Wright } // End Quadrature Point Loop 2493a8779fbSJames Wright 2503a8779fbSJames Wright // Return 2513a8779fbSJames Wright return 0; 2523a8779fbSJames Wright } 253f0b65372SJed Brown 2542b916ea7SJeremy L Thompson CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 2558fff8293SJames Wright return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 25676555becSJames Wright } 25776555becSJames Wright 2582b916ea7SJeremy L Thompson CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 2598fff8293SJames Wright return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 26076555becSJames Wright } 26176555becSJames Wright 262cbe60e31SLeila Ghaffari // ***************************************************************************** 26304e40bb6SJeremy L Thompson // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 264cbe60e31SLeila Ghaffari // ***************************************************************************** 2658fff8293SJames Wright CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 266f0b65372SJed Brown // Inputs 2673d65b166SJames Wright const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 26887bd45e7SJames Wright const CeedScalar(*Grad_dq) = in[1]; 269ade49511SJames Wright const CeedScalar(*q_data) = in[2]; 27094a7b3d2SKenneth E. Jansen const CeedScalar(*jac_data) = in[3]; 2713d65b166SJames Wright 272f0b65372SJed Brown // Outputs 2733d65b166SJames Wright CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 2743d65b166SJames Wright CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 2753d65b166SJames Wright 276f0b65372SJed Brown // Context 277f0b65372SJed Brown NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 278f0b65372SJed Brown const CeedScalar *g = context->g; 279f0b65372SJed Brown 280f0b65372SJed Brown // Quadrature Point Loop 2813d65b166SJames Wright CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 282ade49511SJames Wright CeedScalar wdetJ, dXdx[3][3]; 283ade49511SJames Wright QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 284f0b65372SJed Brown 2858789e95fSJames Wright CeedScalar qi[5], kmstress[6], Tau_d[3]; 286ade49511SJames Wright StoredValuesUnpack(Q, i, 0, 5, jac_data, qi); 287ade49511SJames Wright StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress); 288ade49511SJames Wright StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d); 289edcfef1bSKenneth E. Jansen State s = StateFromQ(context, qi, state_var); 290f0b65372SJed Brown 291edcfef1bSKenneth E. Jansen CeedScalar dqi[5]; 29276555becSJames Wright for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 293edcfef1bSKenneth E. Jansen State ds = StateFromQ_fwd(context, s, dqi, state_var); 294f0b65372SJed Brown 295f0b65372SJed Brown State grad_ds[3]; 296edcfef1bSKenneth E. Jansen StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 297f0b65372SJed Brown 298f0b65372SJed Brown CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 29940a33f2dSJames Wright KMStrainRate_State(grad_ds, dstrain_rate); 300f0b65372SJed Brown NewtonianStress(context, dstrain_rate, dkmstress); 301f0b65372SJed Brown KMUnpack(dkmstress, dstress); 302f0b65372SJed Brown KMUnpack(kmstress, stress); 303f0b65372SJed Brown ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 304f0b65372SJed Brown 305f0b65372SJed Brown StateConservative dF_inviscid[3]; 306f0b65372SJed Brown FluxInviscid_fwd(context, s, ds, dF_inviscid); 307f0b65372SJed Brown 308f0b65372SJed Brown // Total flux 309f0b65372SJed Brown CeedScalar dFlux[5][3]; 310d1b9ef12SLeila Ghaffari FluxTotal(dF_inviscid, dstress, dFe, dFlux); 311f0b65372SJed Brown 31222387d3aSJames Wright for (int j = 0; j < 5; j++) { 31322387d3aSJames Wright for (int k = 0; k < 3; k++) Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * dFlux[j][0] + dXdx[k][1] * dFlux[j][1] + dXdx[k][2] * dFlux[j][2]); 3142b916ea7SJeremy L Thompson } 315f0b65372SJed Brown 31660dbb574SKenneth E. Jansen const CeedScalar dbody_force[5] = {0, ds.U.density * g[0], ds.U.density * g[1], ds.U.density * g[2], Dot3(ds.U.momentum, g)}; 31776555becSJames Wright CeedScalar dU[5] = {0.}; 31876555becSJames Wright UnpackState_U(ds.U, dU); 3192b916ea7SJeremy L Thompson for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 320f0b65372SJed Brown 321e7754af5SKenneth E. Jansen if (context->idl_enable) { 32294a7b3d2SKenneth E. Jansen const CeedScalar sigma = jac_data[14 * Q + i]; 323e7754af5SKenneth E. Jansen CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 324e7754af5SKenneth E. Jansen // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 32594a7b3d2SKenneth E. Jansen InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 326e7754af5SKenneth E. Jansen for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 327e7754af5SKenneth E. Jansen } 328e7754af5SKenneth E. Jansen 329d1b9ef12SLeila Ghaffari // -- Stabilization method: none (Galerkin), SU, or SUPG 330d1b9ef12SLeila Ghaffari CeedScalar dstab[5][3], U_dot[5] = {0}; 331d1b9ef12SLeila Ghaffari for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 332edcfef1bSKenneth E. Jansen Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, dstab); 333d1b9ef12SLeila Ghaffari 3342b916ea7SJeremy L Thompson for (int j = 0; j < 5; j++) { 3352b916ea7SJeremy L Thompson for (int k = 0; k < 3; k++) Grad_v[k][j][i] += wdetJ * (dstab[j][0] * dXdx[k][0] + dstab[j][1] * dXdx[k][1] + dstab[j][2] * dXdx[k][2]); 3362b916ea7SJeremy L Thompson } 337f0b65372SJed Brown } // End Quadrature Point Loop 338f0b65372SJed Brown return 0; 339f0b65372SJed Brown } 3408085925cSJames Wright 3412b916ea7SJeremy L Thompson CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 3428fff8293SJames Wright return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 34376555becSJames Wright } 34476555becSJames Wright 3452b916ea7SJeremy L Thompson CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 3468fff8293SJames Wright return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 34776555becSJames Wright } 34876555becSJames Wright 349d1b9ef12SLeila Ghaffari // ***************************************************************************** 3508085925cSJames Wright // Compute boundary integral (ie. for strongly set inflows) 351d1b9ef12SLeila Ghaffari // ***************************************************************************** 3528fff8293SJames Wright CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 353*4b96a86bSJames Wright const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 3543d65b166SJames Wright const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 35587bd45e7SJames Wright const CeedScalar(*Grad_q) = in[1]; 356ade49511SJames Wright const CeedScalar(*q_data_sur) = in[2]; 3573d65b166SJames Wright CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 358*4b96a86bSJames Wright CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL; 3598085925cSJames Wright 360d3b25f3aSJames Wright const bool is_implicit = context->is_implicit; 3618085925cSJames Wright 3622b916ea7SJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 36341e73928SJames Wright const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 364edcfef1bSKenneth E. Jansen State s = StateFromQ(context, qi, state_var); 3658085925cSJames Wright 366ade49511SJames Wright CeedScalar wdetJb, dXdx[2][3], norm[3]; 367ade49511SJames Wright QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 368ade49511SJames Wright wdetJb *= is_implicit ? -1. : 1.; 3698085925cSJames Wright 370d3b25f3aSJames Wright State grad_s[3]; 371edcfef1bSKenneth E. Jansen StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 3728085925cSJames Wright 373d3b25f3aSJames Wright CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 37440a33f2dSJames Wright KMStrainRate_State(grad_s, strain_rate); 375d3b25f3aSJames Wright NewtonianStress(context, strain_rate, kmstress); 376d3b25f3aSJames Wright KMUnpack(kmstress, stress); 377d3b25f3aSJames Wright ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 378d3b25f3aSJames Wright 379d3b25f3aSJames Wright StateConservative F_inviscid[3]; 380d3b25f3aSJames Wright FluxInviscid(context, s, F_inviscid); 381d3b25f3aSJames Wright 382c5740391SJames Wright CeedScalar Flux[5]; 383c5740391SJames Wright FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 384d3b25f3aSJames Wright 385c5740391SJames Wright for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 3868085925cSJames Wright 387*4b96a86bSJames Wright if (is_implicit) { 388ade49511SJames Wright StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 389ade49511SJames Wright StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 3908085925cSJames Wright } 391*4b96a86bSJames Wright } 3928085925cSJames Wright return 0; 3938085925cSJames Wright } 3948085925cSJames Wright 3952b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 3968fff8293SJames Wright return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 397d4559bbeSJames Wright } 398d4559bbeSJames Wright 3992b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 4008fff8293SJames Wright return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 401d4559bbeSJames Wright } 402d4559bbeSJames Wright 403d1b9ef12SLeila Ghaffari // ***************************************************************************** 40468ae065aSJames Wright // Jacobian for "set nothing" boundary integral 405d1b9ef12SLeila Ghaffari // ***************************************************************************** 4062b916ea7SJeremy L Thompson CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 4078fff8293SJames Wright StateVariable state_var) { 40868ae065aSJames Wright // Inputs 4093d65b166SJames Wright const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 41087bd45e7SJames Wright const CeedScalar(*Grad_dq) = in[1]; 411ade49511SJames Wright const CeedScalar(*q_data_sur) = in[2]; 412c1484fadSKenneth E. Jansen const CeedScalar(*jac_data_sur) = in[4]; 4133d65b166SJames Wright 41468ae065aSJames Wright // Outputs 41568ae065aSJames Wright CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 41668ae065aSJames Wright 41768ae065aSJames Wright const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 418ade49511SJames Wright const bool is_implicit = context->is_implicit; 41968ae065aSJames Wright 42068ae065aSJames Wright // Quadrature Point Loop 4213d65b166SJames Wright CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 422ade49511SJames Wright CeedScalar wdetJb, dXdx[2][3], norm[3]; 423ade49511SJames Wright QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 424ade49511SJames Wright wdetJb *= is_implicit ? -1. : 1.; 42568ae065aSJames Wright 426edcfef1bSKenneth E. Jansen CeedScalar qi[5], kmstress[6], dqi[5]; 427ade49511SJames Wright StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 428ade49511SJames Wright StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 42941e73928SJames Wright for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 4303934e2b1SJames Wright 431edcfef1bSKenneth E. Jansen State s = StateFromQ(context, qi, state_var); 432edcfef1bSKenneth E. Jansen State ds = StateFromQ_fwd(context, s, dqi, state_var); 43368ae065aSJames Wright 43468ae065aSJames Wright State grad_ds[3]; 435edcfef1bSKenneth E. Jansen StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 43668ae065aSJames Wright 43768ae065aSJames Wright CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 43840a33f2dSJames Wright KMStrainRate_State(grad_ds, dstrain_rate); 43968ae065aSJames Wright NewtonianStress(context, dstrain_rate, dkmstress); 44068ae065aSJames Wright KMUnpack(dkmstress, dstress); 44168ae065aSJames Wright KMUnpack(kmstress, stress); 44268ae065aSJames Wright ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 44368ae065aSJames Wright 44468ae065aSJames Wright StateConservative dF_inviscid[3]; 44568ae065aSJames Wright FluxInviscid_fwd(context, s, ds, dF_inviscid); 44668ae065aSJames Wright 447c5740391SJames Wright CeedScalar dFlux[5]; 448c5740391SJames Wright FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 44968ae065aSJames Wright 450c5740391SJames Wright for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 45168ae065aSJames Wright } // End Quadrature Point Loop 45268ae065aSJames Wright return 0; 45368ae065aSJames Wright } 45468ae065aSJames Wright 4552b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 4568fff8293SJames Wright return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 457d4559bbeSJames Wright } 458d4559bbeSJames Wright 4592b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 4608fff8293SJames Wright return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 461d4559bbeSJames Wright } 462