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 /// Density current initial condition and operator for Navier-Stokes example using PETSc
12 //   Semi-Implicit Formulations of the Navier-Stokes Equations: Application to
45 //       delta_theta = r <= rc : thetaC(1 + cos(pi r/rc)) / 2
46 //                     r > rc : 0
47 //         r        = sqrt( (x - xc)**2 + (y - yc)**2 + (z - zc)**2 )
51 //       Pibar      = 1. + g**2 (exp( - N**2 z / g ) - 1) / (cp theta0 N**2)
72 //   N               ,  Brunt-Vaisala frequency
75 //   Rd     = cp - cv,  Specific heat difference
83 // This helper function provides support for the exact, time-dependent solution
89   const CeedScalar            theta0  = context->theta0;  in Exact_DC()
90   const CeedScalar            thetaC  = context->thetaC;  in Exact_DC()
91   const CeedScalar            P0      = context->P0;  in Exact_DC()
92   const CeedScalar            N       = context->N;  in Exact_DC()
93   const CeedScalar            rc      = context->rc;  in Exact_DC()
94   const CeedScalar           *center  = context->center;  in Exact_DC()
95   const CeedScalar           *dc_axis = context->dc_axis;  in Exact_DC()
96   NewtonianIdealGasContext    gas     = &context->newtonian_ctx;  in Exact_DC()
97   const CeedScalar            cp      = gas->cp;  in Exact_DC()
98   const CeedScalar            cv      = gas->cv;  in Exact_DC()
99   const CeedScalar            Rd      = cp - cv;  in Exact_DC()
100   const CeedScalar           *g_vec   = gas->g;  in Exact_DC()
101   const CeedScalar            g       = -g_vec[2];  in Exact_DC()
104   // -- Coordinates  in Exact_DC()
109   // -- Potential temperature, density current  in Exact_DC()
110   CeedScalar rr[3] = {x - center[0], y - center[1], z - center[2]};  in Exact_DC()
111   // (I - q q^T) r: distance from dc_axis (or from center if dc_axis is the zero vector)  in Exact_DC()
112   for (CeedInt i = 0; i < 3; i++) rr[i] -= dc_axis[i] * Dot3(dc_axis, rr);  in Exact_DC()
113   const CeedScalar r           = sqrt(Dot3(rr, rr));  in Exact_DC()  local
114   const CeedScalar delta_theta = r <= rc ? thetaC * (1. + cos(M_PI * r / rc)) / 2. : 0.;  in Exact_DC()
117   // -- Exner pressure, hydrostatic balance  in Exact_DC()
118   const CeedScalar Pi = 1. + Square(g) * (exp(-Square(N) * z / g) - 1.) / (cp * theta0 * Square(N));  in Exact_DC()
139   const NewtonianIdealGasContext gas     = &context->newtonian_ctx;  in ICsDC()
145     StateToQ(gas, s, q, gas->state_var);  in ICsDC()