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