1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Density current initial condition and operator for Navier-Stokes example using PETSc 19 20 // Model from: 21 // Semi-Implicit Formulations of the Navier-Stokes Equations: Application to 22 // Nonhydrostatic Atmospheric Modeling, Giraldo, Restelli, and Lauter (2010). 23 24 #ifndef densitycurrent_h 25 #define densitycurrent_h 26 27 #include <math.h> 28 #include <ceed.h> 29 30 #ifndef M_PI 31 #define M_PI 3.14159265358979323846 32 #endif 33 34 #ifndef setup_context_struct 35 #define setup_context_struct 36 typedef struct SetupContext_ *SetupContext; 37 struct SetupContext_ { 38 CeedScalar theta0; 39 CeedScalar thetaC; 40 CeedScalar P0; 41 CeedScalar N; 42 CeedScalar cv; 43 CeedScalar cp; 44 CeedScalar g; 45 CeedScalar rc; 46 CeedScalar lx; 47 CeedScalar ly; 48 CeedScalar lz; 49 CeedScalar center[3]; 50 CeedScalar dc_axis[3]; 51 CeedScalar wind[3]; 52 CeedScalar time; 53 int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION 54 int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER 55 int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK 56 }; 57 #endif 58 59 // ***************************************************************************** 60 // This function sets the initial conditions and the boundary conditions 61 // 62 // These initial conditions are given in terms of potential temperature and 63 // Exner pressure and then converted to density and total energy. 64 // Initial momentum density is zero. 65 // 66 // Initial Conditions: 67 // Potential Temperature: 68 // theta = thetabar + delta_theta 69 // thetabar = theta0 exp( N**2 z / g ) 70 // delta_theta = r <= rc : thetaC(1 + cos(pi r/rc)) / 2 71 // r > rc : 0 72 // r = sqrt( (x - xc)**2 + (y - yc)**2 + (z - zc)**2 ) 73 // with (xc,yc,zc) center of domain, rc characteristic radius of thermal bubble 74 // Exner Pressure: 75 // Pi = Pibar + deltaPi 76 // Pibar = 1. + g**2 (exp( - N**2 z / g ) - 1) / (cp theta0 N**2) 77 // deltaPi = 0 (hydrostatic balance) 78 // Velocity/Momentum Density: 79 // Ui = ui = 0 80 // 81 // Conversion to Conserved Variables: 82 // rho = P0 Pi**(cv/Rd) / (Rd theta) 83 // E = rho (cv T + (u u)/2 + g z) 84 // 85 // Boundary Conditions: 86 // Mass Density: 87 // 0.0 flux 88 // Momentum Density: 89 // 0.0 90 // Energy Density: 91 // 0.0 flux 92 // 93 // Constants: 94 // theta0 , Potential temperature constant 95 // thetaC , Potential temperature perturbation 96 // P0 , Pressure at the surface 97 // N , Brunt-Vaisala frequency 98 // cv , Specific heat, constant volume 99 // cp , Specific heat, constant pressure 100 // Rd = cp - cv, Specific heat difference 101 // g , Gravity 102 // rc , Characteristic radius of thermal bubble 103 // center , Location of bubble center 104 // dc_axis , Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric 105 // ***************************************************************************** 106 107 // ***************************************************************************** 108 // This helper function provides support for the exact, time-dependent solution 109 // (currently not implemented) and IC formulation for density current 110 // ***************************************************************************** 111 CEED_QFUNCTION_HELPER int Exact_DC(CeedInt dim, CeedScalar time, 112 const CeedScalar X[], CeedInt Nf, CeedScalar q[], 113 void *ctx) { 114 // Context 115 const SetupContext context = (SetupContext)ctx; 116 const CeedScalar theta0 = context->theta0; 117 const CeedScalar thetaC = context->thetaC; 118 const CeedScalar P0 = context->P0; 119 const CeedScalar N = context->N; 120 const CeedScalar cv = context->cv; 121 const CeedScalar cp = context->cp; 122 const CeedScalar g = context->g; 123 const CeedScalar rc = context->rc; 124 const CeedScalar *center = context->center; 125 const CeedScalar *dc_axis = context->dc_axis; 126 const CeedScalar Rd = cp - cv; 127 128 // Setup 129 // -- Coordinates 130 const CeedScalar x = X[0]; 131 const CeedScalar y = X[1]; 132 const CeedScalar z = X[2]; 133 134 // -- Potential temperature, density current 135 CeedScalar rr[3] = {x - center[0], y - center[1], z - center[2]}; 136 // (I - q q^T) r: distance from dc_axis (or from center if dc_axis is the zero vector) 137 for (CeedInt i=0; i<3; i++) 138 rr[i] -= dc_axis[i] * 139 (dc_axis[0]*rr[0] + dc_axis[1]*rr[1] + dc_axis[2]*rr[2]); 140 const CeedScalar r = sqrt(rr[0]*rr[0] + rr[1]*rr[1] + rr[2]*rr[2]); 141 const CeedScalar delta_theta = r <= rc ? thetaC*(1. + cos(M_PI*r/rc))/2. : 0.; 142 const CeedScalar theta = theta0*exp(N*N*z/g) + delta_theta; 143 144 // -- Exner pressure, hydrostatic balance 145 const CeedScalar Pi = 1. + g*g*(exp(-N*N*z/g) - 1.) / (cp*theta0*N*N); 146 // -- Density 147 148 const CeedScalar rho = P0 * pow(Pi, cv/Rd) / (Rd*theta); 149 150 // Initial Conditions 151 q[0] = rho; 152 q[1] = 0.0; 153 q[2] = 0.0; 154 q[3] = 0.0; 155 q[4] = rho * (cv*theta*Pi + g*z); 156 157 return 0; 158 } 159 160 // ***************************************************************************** 161 // This QFunction sets the initial conditions for density current 162 // ***************************************************************************** 163 CEED_QFUNCTION(ICsDC)(void *ctx, CeedInt Q, 164 const CeedScalar *const *in, CeedScalar *const *out) { 165 // Inputs 166 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 167 168 // Outputs 169 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 170 171 CeedPragmaSIMD 172 // Quadrature Point Loop 173 for (CeedInt i=0; i<Q; i++) { 174 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 175 CeedScalar q[5] = {0.}; 176 177 Exact_DC(3, 0., x, 5, q, ctx); 178 179 for (CeedInt j=0; j<5; j++) 180 q0[j][i] = q[j]; 181 } // End of Quadrature Point Loop 182 183 return 0; 184 } 185 186 #endif // densitycurrent_h 187