Lines Matching +full:- +full:- +full:substitute
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 /// Euler traveling vortex initial condition and operator for Navier-Stokes
38 // T = 1 - (gamma - 1) vortex_strength**2 exp(1 - r**2) / (8 gamma pi**2)
40 // rho = (T/S_vortex)^(1 / (gamma - 1))
44 // ui = 1 + vortex_strength exp((1 - r**2)/2.) [yc - y, x - xc] / (2 pi)
45 // r = sqrt( (x - xc)**2 + (y - yc)**2 )
49 // E = P / (gamma - 1) + rho (u u)/2
61 // This helper function provides support for the exact, time-dependent solution (currently not impl…
67 const CeedScalar vortex_strength = context->vortex_strength; in Exact_Euler()
68 const CeedScalar *center = context->center; // Center of the domain in Exact_Euler()
69 const CeedScalar *mean_velocity = context->mean_velocity; in Exact_Euler()
80 const CeedScalar x0 = x - xc; in Exact_Euler()
81 const CeedScalar y0 = y - yc; in Exact_Euler()
83 const CeedScalar C = vortex_strength * exp((1. - r * r) / 2.) / (2. * M_PI); in Exact_Euler()
84 …const CeedScalar delta_T = -(gamma - 1.) * vortex_strength * vortex_strength * exp(1 - r * r) / (… in Exact_Euler()
86 …const CeedScalar S_bubble = (gamma - 1.) * vortex_strength * vortex_strength / (8. * gamma * M_PI … in Exact_Euler()
90 switch (context->euler_test) { in Exact_Euler()
95 // Solve for rho, then substitute for P in Exact_Euler()
96 rho = pow(T / S_vortex, 1 / (gamma - 1.)); in Exact_Euler()
98 u[0] = mean_velocity[0] - C * y0; in Exact_Euler()
106 q[4] = P / (gamma - 1.) + rho * (u[0] * u[0] + u[1] * u[1]) / 2.; in Exact_Euler()
132 …case 3: // Velocity zero, pressure constant (so density and internal energy will be non-constant)… in Exact_Euler()
136 T = 1. - S_bubble * exp(1. - r * r); in Exact_Euler()
146 …Constant nonzero velocity, pressure constant (so density and internal energy will be non-constant), in Exact_Euler()
149 T = 1. - S_bubble * exp(1. - r * r); in Exact_Euler()
161 case 5: // non-smooth thermal bubble - cylinder in Exact_Euler()
163 T = 1. - (r < 1. ? S_bubble : 0.); in Exact_Euler()
187 dF[i][j + 1][0] = ((i == j) ? ((gamma - 1.) * (u_sq / 2.)) : 0.) - u[i] * u[j]; in ConvectiveFluxJacobian_Euler()
190 …][k + 1] = ((j == k) ? u[i] : 0.) + ((i == k) ? u[j] : 0.) - ((i == j) ? u[k] : 0.) * (gamma - 1.); in ConvectiveFluxJacobian_Euler()
191 …dF[i][4][k + 1] = ((i == k) ? (E * gamma / rho - (gamma - 1.) * u_sq / 2.) : 0.) - (gamma - 1.… in ConvectiveFluxJacobian_Euler()
193 dF[i][j + 1][4] = ((i == j) ? (gamma - 1.) : 0.); in ConvectiveFluxJacobian_Euler()
195 dF[i][4][0] = u[i] * ((gamma - 1.) * u_sq - E * gamma / rho); in ConvectiveFluxJacobian_Euler()
205 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix
212 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number )
239 Exact_Euler(3, context->curr_time, x, 5, q, ctx); in ICsEuler()
251 // rho - Mass Density
252 // Ui - Momentum Density, Ui = rho ui
253 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2
261 // P = (gamma - 1) (E - rho (u u) / 2)
277 const CeedScalar c_tau = context->c_tau; in Euler()
282 // -- Interp in in Euler()
311 …alar E_kinetic = 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]), E_internal = E - E_kinetic, in Euler()
312 P = E_internal * (gamma - 1.); // P = pressure in Euler()
321 // -- Density in Euler()
322 // ---- u rho in Euler()
324 // -- Momentum in Euler()
325 // ---- rho (u x u) + P I3 in Euler()
332 // -- Total Energy Density in Euler()
333 // ---- (E + P) u in Euler()
336 // --Stabilization terms in Euler()
337 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction in Euler()
341 // ---- dqdx collects drhodx, dUdx and dEdx in one vector in Euler()
349 // ---- strong_conv = dF/dq * dq/dx (Strong convection) in Euler()
358 // -- Tau elements in Euler()
363 // -- Stabilization method: none or SU in Euler()
365 switch (context->stabilization) { in Euler()
376 …for (CeedInt k = 0; k < 3; k++) dv[k][j][i] -= wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXd… in Euler()
399 const CeedScalar c_tau = context->c_tau; in IFunction_Euler()
405 // -- Interp in in IFunction_Euler()
433 …alar E_kinetic = 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]), E_internal = E - E_kinetic, in IFunction_Euler()
434 P = E_internal * (gamma - 1.); // P = pressure in IFunction_Euler()
442 //-----mass matrix in IFunction_Euler()
445 // -- Density in IFunction_Euler()
446 // ---- u rho in IFunction_Euler()
447 …for (CeedInt j = 0; j < 3; j++) dv[j][0][i] -= wdetJ * (rho * u[0] * dXdx[j][0] + rho * u[1] * dXd… in IFunction_Euler()
448 // -- Momentum in IFunction_Euler()
449 // ---- rho (u x u) + P I3 in IFunction_Euler()
452 …dv[k][j + 1][i] -= wdetJ * ((rho * u[j] * u[0] + (j == 0 ? P : 0.)) * dXdx[k][0] + (rho * u[j] * u… in IFunction_Euler()
456 // -- Total Energy Density in IFunction_Euler()
457 // ---- (E + P) u in IFunction_Euler()
458 …for (CeedInt j = 0; j < 3; j++) dv[j][4][i] -= wdetJ * (E + P) * (u[0] * dXdx[j][0] + u[1] * dXdx[… in IFunction_Euler()
460 // -- Stabilization terms in IFunction_Euler()
461 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction in IFunction_Euler()
465 // ---- dqdx collects drhodx, dUdx and dEdx in one vector in IFunction_Euler()
473 // ---- strong_conv = dF/dq * dq/dx (Strong convection) in IFunction_Euler()
481 // ---- Strong residual in IFunction_Euler()
486 // -- Tau elements in IFunction_Euler()
491 // -- Stabilization method: none, SU, or SUPG in IFunction_Euler()
493 switch (context->stabilization) { in IFunction_Euler()
533 const int euler_test = context->euler_test; in TravelingVortex_Inflow()
534 const bool is_implicit = context->implicit; in TravelingVortex_Inflow()
535 CeedScalar *mean_velocity = context->mean_velocity; in TravelingVortex_Inflow()
553 wdetJb *= is_implicit ? -1. : 1.; in TravelingVortex_Inflow()
570 // -- Density in TravelingVortex_Inflow()
571 v[0][i] -= wdetJb * rho_inlet * face_normal; in TravelingVortex_Inflow()
573 // -- Momentum in TravelingVortex_Inflow()
574 …for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho_inlet * face_normal * mean_velocity[j… in TravelingVortex_Inflow()
576 // -- Total Energy Density in TravelingVortex_Inflow()
577 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet); in TravelingVortex_Inflow()
595 const bool is_implicit = context->implicit; in Euler_Outflow()
596 CeedScalar *mean_velocity = context->mean_velocity; in Euler_Outflow()
602 // -- Interp in in Euler_Outflow()
609 wdetJb *= is_implicit ? -1. : 1.; in Euler_Outflow()
620 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure in Euler_Outflow()
623 // -- Density in Euler_Outflow()
624 v[0][i] -= wdetJb * rho * u_normal; in Euler_Outflow()
626 // -- Momentum in Euler_Outflow()
627 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho * u_normal * u[j] + norm[j] * P); in Euler_Outflow()
629 // -- Total Energy Density in Euler_Outflow()
630 v[4][i] -= wdetJb * u_normal * (E + P); in Euler_Outflow()