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