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.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Density current initial condition and operator for Navier-Stokes example using PETSc
10
11 // Model from:
12 // Semi-Implicit Formulations of the Navier-Stokes Equations: Application to
13 // Nonhydrostatic Atmospheric Modeling, Giraldo, Restelli, and Lauter (2010).
14 #include <ceed/types.h>
15 #ifndef CEED_RUNNING_JIT_PASS
16 #include <math.h>
17 #endif
18
19 #include "newtonian_state.h"
20 #include "newtonian_types.h"
21 #include "utils.h"
22
23 typedef struct DensityCurrentContext_ *DensityCurrentContext;
24 struct DensityCurrentContext_ {
25 CeedScalar theta0;
26 CeedScalar thetaC;
27 CeedScalar P0;
28 CeedScalar N;
29 CeedScalar rc;
30 CeedScalar center[3];
31 CeedScalar dc_axis[3];
32 struct NewtonianIdealGasContext_ newtonian_ctx;
33 };
34
35 // *****************************************************************************
36 // This function sets the initial conditions and the boundary conditions
37 //
38 // These initial conditions are given in terms of potential temperature and Exner pressure and then converted to density and total energy.
39 // Initial momentum density is zero.
40 //
41 // Initial Conditions:
42 // Potential Temperature:
43 // theta = thetabar + delta_theta
44 // thetabar = theta0 exp( N**2 z / g )
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 )
48 // with (xc,yc,zc) center of domain, rc characteristic radius of thermal bubble
49 // Exner Pressure:
50 // Pi = Pibar + deltaPi
51 // Pibar = 1. + g**2 (exp( - N**2 z / g ) - 1) / (cp theta0 N**2)
52 // deltaPi = 0 (hydrostatic balance)
53 // Velocity/Momentum Density:
54 // Ui = ui = 0
55 //
56 // Conversion to Conserved Variables:
57 // rho = P0 Pi**(cv/Rd) / (Rd theta)
58 // E = rho (cv T + (u u)/2 + g z)
59 //
60 // Boundary Conditions:
61 // Mass Density:
62 // 0.0 flux
63 // Momentum Density:
64 // 0.0
65 // Energy Density:
66 // 0.0 flux
67 //
68 // Constants:
69 // theta0 , Potential temperature constant
70 // thetaC , Potential temperature perturbation
71 // P0 , Pressure at the surface
72 // N , Brunt-Vaisala frequency
73 // cv , Specific heat, constant volume
74 // cp , Specific heat, constant pressure
75 // Rd = cp - cv, Specific heat difference
76 // g , Gravity
77 // rc , Characteristic radius of thermal bubble
78 // center , Location of bubble center
79 // dc_axis , Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric
80 // *****************************************************************************
81
82 // *****************************************************************************
83 // This helper function provides support for the exact, time-dependent solution
84 // (currently not implemented) and IC formulation for density current
85 // *****************************************************************************
Exact_DC(CeedInt dim,CeedScalar time,const CeedScalar X[],CeedInt Nf,void * ctx)86 CEED_QFUNCTION_HELPER State Exact_DC(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, void *ctx) {
87 // Context
88 const DensityCurrentContext context = (DensityCurrentContext)ctx;
89 const CeedScalar theta0 = context->theta0;
90 const CeedScalar thetaC = context->thetaC;
91 const CeedScalar P0 = context->P0;
92 const CeedScalar N = context->N;
93 const CeedScalar rc = context->rc;
94 const CeedScalar *center = context->center;
95 const CeedScalar *dc_axis = context->dc_axis;
96 NewtonianIdealGasContext gas = &context->newtonian_ctx;
97 const CeedScalar cp = gas->cp;
98 const CeedScalar cv = gas->cv;
99 const CeedScalar Rd = cp - cv;
100 const CeedScalar *g_vec = gas->g;
101 const CeedScalar g = -g_vec[2];
102
103 // Setup
104 // -- Coordinates
105 const CeedScalar x = X[0];
106 const CeedScalar y = X[1];
107 const CeedScalar z = X[2];
108
109 // -- Potential temperature, density current
110 CeedScalar rr[3] = {x - center[0], y - center[1], z - center[2]};
111 // (I - q q^T) r: distance from dc_axis (or from center if dc_axis is the zero vector)
112 for (CeedInt i = 0; i < 3; i++) rr[i] -= dc_axis[i] * Dot3(dc_axis, rr);
113 const CeedScalar r = sqrt(Dot3(rr, rr));
114 const CeedScalar delta_theta = r <= rc ? thetaC * (1. + cos(M_PI * r / rc)) / 2. : 0.;
115 const CeedScalar theta = theta0 * exp(Square(N) * z / g) + delta_theta;
116
117 // -- Exner pressure, hydrostatic balance
118 const CeedScalar Pi = 1. + Square(g) * (exp(-Square(N) * z / g) - 1.) / (cp * theta0 * Square(N));
119
120 // Initial Conditions
121 CeedScalar Y[5] = {0.};
122 Y[0] = P0 * pow(Pi, cp / Rd);
123 Y[1] = 0.0;
124 Y[2] = 0.0;
125 Y[3] = 0.0;
126 Y[4] = Pi * theta;
127
128 return StateFromY(gas, Y);
129 }
130
131 // *****************************************************************************
132 // This QFunction sets the initial conditions for density current
133 // *****************************************************************************
ICsDC(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)134 CEED_QFUNCTION(ICsDC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
135 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
136 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
137
138 const DensityCurrentContext context = (DensityCurrentContext)ctx;
139 const NewtonianIdealGasContext gas = &context->newtonian_ctx;
140
141 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
142 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
143 State s = Exact_DC(3, 0., x, 5, ctx);
144 CeedScalar q[5] = {0};
145 StateToQ(gas, s, q, gas->state_var);
146 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
147 }
148 return 0;
149 }
150