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