Lines Matching +full:- +full:j
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()
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()
186 for (CeedInt j = 0; j < 3; j++) { // Rows of each Jacobian matrix in ConvectiveFluxJacobian_Euler() local
187 dF[i][j + 1][0] = ((i == j) ? ((gamma - 1.) * (u_sq / 2.)) : 0.) - u[i] * u[j]; in ConvectiveFluxJacobian_Euler()
190 …dF[i][j + 1][k + 1] = ((j == k) ? u[i] : 0.) + ((i == k) ? u[j] : 0.) - ((i == j) ? u[k] : 0.) * (… 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()
240 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; in ICsEuler() local
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()
300 for (CeedInt j = 0; j < 3; j++) { in Euler() local
302 drhodx[j] += drho[k] * dXdx[k][j]; in Euler()
303 dEdx[j] += dE[k] * dXdx[k][j]; in Euler()
305 dUdx[j][k] += dU[j][l] * dXdx[l][k]; in Euler()
306 dXdxdXdxT[j][k] += dXdx[j][l] * dXdx[k][l]; // dXdx_j,k * dXdx_k,j 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()
316 for (CeedInt j = 0; j < 5; j++) { in Euler() local
317 v[j][i] = 0.; in Euler()
318 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] = 0.; in Euler()
321 // -- Density in Euler()
322 // ---- u rho in Euler()
323 …for (CeedInt j = 0; j < 3; j++) dv[j][0][i] += wdetJ * (rho * u[0] * dXdx[j][0] + rho * u[1] * dXd… in Euler() local
324 // -- Momentum in Euler()
325 // ---- rho (u x u) + P I3 in Euler()
326 for (CeedInt j = 0; j < 3; j++) { in Euler() local
328 …dv[k][j + 1][i] += wdetJ * ((rho * u[j] * u[0] + (j == 0 ? P : 0.)) * dXdx[k][0] + (rho * u[j] * u… in Euler()
329 (rho * u[j] * u[2] + (j == 2 ? P : 0.)) * dXdx[k][2]); in Euler()
332 // -- Total Energy Density in Euler()
333 // ---- (E + P) u in Euler()
334 …for (CeedInt j = 0; j < 3; j++) dv[j][4][i] += wdetJ * (E + P) * (u[0] * dXdx[j][0] + u[1] * dXdx[… in Euler() local
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()
343 for (CeedInt j = 0; j < 3; j++) { in Euler() local
344 dqdx[0][j] = drhodx[j]; in Euler()
345 dqdx[4][j] = dEdx[j]; in Euler()
346 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j]; in Euler()
349 // ---- strong_conv = dF/dq * dq/dx (Strong convection) in Euler()
351 for (CeedInt j = 0; j < 3; j++) { in Euler() local
353 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; in Euler()
358 // -- Tau elements in Euler()
363 // -- Stabilization method: none or SU in Euler()
365 switch (context->stabilization) { in Euler()
369 for (CeedInt j = 0; j < 3; j++) { in Euler() local
371 … for (CeedInt l = 0; l < 5; l++) stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; in Euler()
375 for (CeedInt j = 0; j < 5; j++) { in Euler() local
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()
423 for (CeedInt j = 0; j < 3; j++) { in IFunction_Euler() local
425 drhodx[j] += drho[k] * dXdx[k][j]; in IFunction_Euler()
426 dEdx[j] += dE[k] * dXdx[k][j]; in IFunction_Euler()
428 dUdx[j][k] += dU[j][l] * dXdx[l][k]; in IFunction_Euler()
429 dXdxdXdxT[j][k] += dXdx[j][l] * dXdx[k][l]; // dXdx_j,k * dXdx_k,j 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()
438 for (CeedInt j = 0; j < 5; j++) { in IFunction_Euler() local
439 v[j][i] = 0.; in IFunction_Euler()
440 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] = 0.; in IFunction_Euler()
442 //-----mass matrix in IFunction_Euler()
443 for (CeedInt j = 0; j < 5; j++) v[j][i] += wdetJ * q_dot[j][i]; in IFunction_Euler() local
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() local
448 // -- Momentum in IFunction_Euler()
449 // ---- rho (u x u) + P I3 in IFunction_Euler()
450 for (CeedInt j = 0; j < 3; j++) { in IFunction_Euler() local
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()
453 (rho * u[j] * u[2] + (j == 2 ? P : 0.)) * dXdx[k][2]); 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() local
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()
467 for (CeedInt j = 0; j < 3; j++) { in IFunction_Euler() local
468 dqdx[0][j] = drhodx[j]; in IFunction_Euler()
469 dqdx[4][j] = dEdx[j]; in IFunction_Euler()
470 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j]; in IFunction_Euler()
473 // ---- strong_conv = dF/dq * dq/dx (Strong convection) in IFunction_Euler()
475 for (CeedInt j = 0; j < 3; j++) { in IFunction_Euler() local
477 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; in IFunction_Euler()
481 // ---- Strong residual in IFunction_Euler()
483 for (CeedInt j = 0; j < 5; j++) strong_res[j] = q_dot[j][i] + strong_conv[j]; in IFunction_Euler() local
486 // -- Tau elements in IFunction_Euler()
491 // -- Stabilization method: none, SU, or SUPG in IFunction_Euler()
493 switch (context->stabilization) { in IFunction_Euler()
497 for (CeedInt j = 0; j < 3; j++) { in IFunction_Euler() local
499 … for (CeedInt l = 0; l < 5; l++) stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; in IFunction_Euler()
503 for (CeedInt j = 0; j < 5; j++) { in IFunction_Euler() local
504 …for (CeedInt k = 0; k < 3; k++) dv[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXd… in IFunction_Euler()
508 for (CeedInt j = 0; j < 3; j++) { in IFunction_Euler() local
510 … for (CeedInt l = 0; l < 5; l++) stab[k][j] = jacob_F_conv[j][k][l] * Tau_x[j] * strong_res[l]; in IFunction_Euler()
514 for (CeedInt j = 0; j < 5; j++) { in IFunction_Euler() local
515 …for (CeedInt k = 0; k < 3; k++) dv[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXd… 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()
559 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.; in TravelingVortex_Inflow() local
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() local
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()
615 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0; in Euler_Outflow() local
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() local
629 // -- Total Energy Density in Euler_Outflow()
630 v[4][i] -= wdetJb * u_normal * (E + P); in Euler_Outflow()