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 /// Advection initial condition and operator for Navier-Stokes example using PETSc
26 // -- ROTATION (default)
34 //            increases to (1.-r/rc), then 0. everywhere else
44 // -- TRANSLATION
52 //            increases to (1.-r/rc), then 0. everywhere else
70 // This helper function provides the exact, time-dependent solution and IC formulation for 2D advec…
74   const CeedScalar      rc      = context->rc;  in Exact_AdvectionGeneric()
75   const CeedScalar      lx      = context->lx;  in Exact_AdvectionGeneric()
76   const CeedScalar      ly      = context->ly;  in Exact_AdvectionGeneric()
77   const CeedScalar      lz      = dim == 2 ? 0. : context->lz;  in Exact_AdvectionGeneric()
78   const CeedScalar     *wind    = context->wind;  in Exact_AdvectionGeneric()
86   CeedScalar r = 0.;  in Exact_AdvectionGeneric()  local
87   switch (context->initial_condition_type) {  in Exact_AdvectionGeneric()
90       r = sqrt(Square(x - x0[0]) + Square(y - x0[1]) + Square(z - x0[2]));  in Exact_AdvectionGeneric()
93       r = sqrt(Square(x - center[0]) + Square(y - center[1]));  in Exact_AdvectionGeneric()
99   switch (context->wind_type) {  in Exact_AdvectionGeneric()
102       q[1] = -(y - center[1]);  in Exact_AdvectionGeneric()
103       q[2] = (x - center[0]);  in Exact_AdvectionGeneric()
116   switch (context->initial_condition_type) {  in Exact_AdvectionGeneric()
119       switch (context->bubble_continuity_type) {  in Exact_AdvectionGeneric()
122           q[4] = r <= rc ? (1. - r / rc) : 0.;  in Exact_AdvectionGeneric()
126           q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) : 0.;  in Exact_AdvectionGeneric()
130 …    q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) * fmin(1.0, (center[1] - y) / 1.25) : 0.;  in Exact_AdvectionGeneric()
133           q[4] = r <= rc ? .5 + .5 * cos(r * M_PI / rc) : 0;  in Exact_AdvectionGeneric()
138       CeedScalar half_width = context->lx / 2;  in Exact_AdvectionGeneric()
139       q[4]                  = r > half_width ? 0. : cos(2 * M_PI * r / half_width + M_PI) + 1.;  in Exact_AdvectionGeneric()
143       CeedScalar       inflow_to_point[3] = {x - context->lx / 2, y, 0};  in Exact_AdvectionGeneric()
149 …x < boundary_threshold && wind[0] < boundary_threshold) ||                // outflow at -x boundary  in Exact_AdvectionGeneric()
150 …y < boundary_threshold && wind[1] < boundary_threshold) ||                // outflow at -y boundary  in Exact_AdvectionGeneric()
151 …(x > context->lx - boundary_threshold && wind[0] > boundary_threshold) ||  // outflow at +x bounda…  in Exact_AdvectionGeneric()
152 …(y > context->ly - boundary_threshold && wind[1] > boundary_threshold)     // outflow at +y bounda…  in Exact_AdvectionGeneric()
190     Exact_AdvectionGeneric(2, context->time, x, 5, q, ctx);  in ICsAdvection2d()
258   switch (context->stabilization_tau) {  in Tau()
263       return context->CtauS / sqrt(DotN(uX, uX, dim));  in Tau()
270 …return 1 / sqrt(Square(2 * context->Ctau_t / context->dt) + DotN(s.Y.velocity, gij_uj, dim) * cont…  in Tau()
326     if (context->strong_form) {  // Strong Galerkin convection term: v div(E u)  in IFunction_AdvectionGeneric()
328     } else {  // Weak Galerkin convection term: -dv \cdot (E u)  in IFunction_AdvectionGeneric()
329       for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = -wdetJ * s.U.E_total * uX[j];  in IFunction_AdvectionGeneric()
335       for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total;  in IFunction_AdvectionGeneric()
337       for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] -= wdetJ * Fe_dXdx[k];  in IFunction_AdvectionGeneric()
341     for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) {  in IFunction_AdvectionGeneric()
395     for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) {  in MassFunction_AdvectionGeneric()
459     if (context->strong_form) {  // Strong Galerkin convection term: v div(E u)  in RHSFunction_AdvectionGeneric()
460       v[4][i] = -wdetJ * strong_conv;  in RHSFunction_AdvectionGeneric()
462     } else {  // Weak Galerkin convection term: -dv \cdot (E u)  in RHSFunction_AdvectionGeneric()
470       for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total;  in RHSFunction_AdvectionGeneric()
476     for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) {  in RHSFunction_AdvectionGeneric()
481           grad_v[j][4][i] -= wdetJ * TauS * strong_conv * uX[j];  in RHSFunction_AdvectionGeneric()
517   const CeedScalar E_wind      = context->E_wind;  in Advection_InOutFlowGeneric()
518   const CeedScalar strong_form = context->strong_form;  in Advection_InOutFlowGeneric()
519   const bool       is_implicit = context->implicit;  in Advection_InOutFlowGeneric()
528     wdetJb *= is_implicit ? -1. : 1.;  in Advection_InOutFlowGeneric()
538       v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal;  in Advection_InOutFlowGeneric()
540       v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal;  in Advection_InOutFlowGeneric()