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