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