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 /// Utility functions for setting up problems using the Newtonian Qfunction 10 11 #include "../navierstokes.h" 12 #include "../qfunctions/setupgeo.h" 13 #include "../qfunctions/newtonian.h" 14 15 PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *setup_ctx, 16 void *ctx) { 17 SetupContext setup_context = *(SetupContext *)setup_ctx; 18 User user = *(User *)ctx; 19 StabilizationType stab; 20 MPI_Comm comm = PETSC_COMM_WORLD; 21 PetscBool implicit; 22 PetscBool has_curr_time = PETSC_FALSE; 23 PetscInt ierr; 24 PetscFunctionBeginUser; 25 26 ierr = PetscCalloc1(1, &user->phys->newtonian_ig_ctx); CHKERRQ(ierr); 27 28 // ------------------------------------------------------ 29 // Setup Generic Newtonian IG Problem 30 // ------------------------------------------------------ 31 problem->dim = 3; 32 problem->q_data_size_vol = 10; 33 problem->q_data_size_sur = 4; 34 problem->setup_vol = Setup; 35 problem->setup_vol_loc = Setup_loc; 36 problem->ics = ICsNewtonianIG; 37 problem->ics_loc = ICsNewtonianIG_loc; 38 problem->setup_sur = SetupBoundary; 39 problem->setup_sur_loc = SetupBoundary_loc; 40 problem->apply_vol_rhs = Newtonian; 41 problem->apply_vol_rhs_loc = Newtonian_loc; 42 problem->apply_vol_ifunction = IFunction_Newtonian; 43 problem->apply_vol_ifunction_loc = IFunction_Newtonian_loc; 44 problem->setup_ctx = SetupContext_DENSITY_CURRENT; 45 problem->non_zero_time = PETSC_FALSE; 46 problem->print_info = PRINT_DENSITY_CURRENT; 47 48 // ------------------------------------------------------ 49 // Create the libCEED context 50 // ------------------------------------------------------ 51 CeedScalar theta0 = 300.; // K 52 CeedScalar thetaC = -15.; // K 53 CeedScalar P0 = 1.e5; // Pa 54 CeedScalar N = 0.01; // 1/s 55 CeedScalar cv = 717.; // J/(kg K) 56 CeedScalar cp = 1004.; // J/(kg K) 57 CeedScalar g = 9.81; // m/s^2 58 CeedScalar lambda = -2./3.; // - 59 CeedScalar mu = 75.; // Pa s, dynamic viscosity 60 // mu = 75 is not physical for air, but is good for numerical stability 61 CeedScalar k = 0.02638; // W/(m K) 62 CeedScalar c_tau = 0.5; // - 63 // c_tau = 0.5 is reported as "optimal" in Hughes et al 2010 64 PetscReal domain_min[3], domain_max[3], domain_size[3]; 65 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 66 for (int i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 67 68 // ------------------------------------------------------ 69 // Create the PETSc context 70 // ------------------------------------------------------ 71 PetscScalar meter = 1e-2; // 1 meter in scaled length units 72 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 73 PetscScalar second = 1e-2; // 1 second in scaled time units 74 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 75 PetscScalar W_per_m_K, Pascal, J_per_kg_K, m_per_squared_s; 76 77 // ------------------------------------------------------ 78 // Command line Options 79 // ------------------------------------------------------ 80 ierr = PetscOptionsBegin(comm, NULL, 81 "Options for Newtonian Ideal Gas based problem", 82 NULL); CHKERRQ(ierr); 83 // -- Physics 84 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 85 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 86 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 87 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 88 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 89 NULL, P0, &P0, NULL); CHKERRQ(ierr); 90 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 91 NULL, N, &N, NULL); CHKERRQ(ierr); 92 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 93 NULL, cv, &cv, NULL); CHKERRQ(ierr); 94 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 95 NULL, cp, &cp, NULL); CHKERRQ(ierr); 96 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 97 NULL, g, &g, NULL); CHKERRQ(ierr); 98 ierr = PetscOptionsScalar("-lambda", 99 "Stokes hypothesis second viscosity coefficient", 100 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 101 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 102 NULL, mu, &mu, NULL); CHKERRQ(ierr); 103 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 104 NULL, k, &k, NULL); CHKERRQ(ierr); 105 106 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 107 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 108 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 109 ierr = PetscOptionsScalar("-c_tau", "Stabilization constant", 110 NULL, c_tau, &c_tau, NULL); CHKERRQ(ierr); 111 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 112 NULL, implicit=PETSC_FALSE, &implicit, NULL); 113 CHKERRQ(ierr); 114 115 // -- Units 116 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 117 NULL, meter, &meter, NULL); CHKERRQ(ierr); 118 meter = fabs(meter); 119 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 120 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 121 kilogram = fabs(kilogram); 122 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 123 NULL, second, &second, NULL); CHKERRQ(ierr); 124 second = fabs(second); 125 ierr = PetscOptionsScalar("-units_Kelvin", 126 "1 Kelvin in scaled temperature units", 127 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 128 Kelvin = fabs(Kelvin); 129 130 // -- Warnings 131 if (stab == STAB_SUPG && !implicit) { 132 ierr = PetscPrintf(comm, 133 "Warning! Use -stab supg only with -implicit\n"); 134 CHKERRQ(ierr); 135 } 136 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 137 138 // ------------------------------------------------------ 139 // Set up the PETSc context 140 // ------------------------------------------------------ 141 // -- Define derived units 142 Pascal = kilogram / (meter * PetscSqr(second)); 143 J_per_kg_K = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 144 m_per_squared_s = meter / PetscSqr(second); 145 W_per_m_K = kilogram * meter / (pow(second,3) * Kelvin); 146 147 user->units->meter = meter; 148 user->units->kilogram = kilogram; 149 user->units->second = second; 150 user->units->Kelvin = Kelvin; 151 user->units->Pascal = Pascal; 152 user->units->J_per_kg_K = J_per_kg_K; 153 user->units->m_per_squared_s = m_per_squared_s; 154 user->units->W_per_m_K = W_per_m_K; 155 156 // ------------------------------------------------------ 157 // Set up the libCEED context 158 // ------------------------------------------------------ 159 // -- Scale variables to desired units 160 theta0 *= Kelvin; 161 thetaC *= Kelvin; 162 P0 *= Pascal; 163 N *= (1./second); 164 cv *= J_per_kg_K; 165 cp *= J_per_kg_K; 166 g *= m_per_squared_s; 167 mu *= Pascal * second; 168 k *= W_per_m_K; 169 for (int i=0; i<3; i++) domain_size[i] *= meter; 170 problem->dm_scale = meter; 171 172 // -- Setup Context 173 setup_context->theta0 = theta0; 174 setup_context->thetaC = thetaC; 175 setup_context->P0 = P0; 176 setup_context->N = N; 177 setup_context->cv = cv; 178 setup_context->cp = cp; 179 setup_context->g = g; 180 setup_context->lx = domain_size[0]; 181 setup_context->ly = domain_size[1]; 182 setup_context->lz = domain_size[2]; 183 setup_context->time = 0; 184 185 // -- Solver Settings 186 user->phys->stab = stab; 187 user->phys->implicit = implicit; 188 user->phys->has_curr_time = has_curr_time; 189 190 // -- QFunction Context 191 user->phys->newtonian_ig_ctx->lambda = lambda; 192 user->phys->newtonian_ig_ctx->mu = mu; 193 user->phys->newtonian_ig_ctx->k = k; 194 user->phys->newtonian_ig_ctx->cv = cv; 195 user->phys->newtonian_ig_ctx->cp = cp; 196 user->phys->newtonian_ig_ctx->g = g; 197 user->phys->newtonian_ig_ctx->c_tau = c_tau; 198 user->phys->newtonian_ig_ctx->stabilization = stab; 199 200 PetscFunctionReturn(0); 201 } 202 203 PetscErrorCode SetupContext_NEWTONIAN_IG(Ceed ceed, CeedData ceed_data, 204 AppCtx app_ctx, SetupContext setup_ctx, Physics phys) { 205 PetscFunctionBeginUser; 206 CeedQFunctionContextCreate(ceed, &ceed_data->setup_context); 207 CeedQFunctionContextSetData(ceed_data->setup_context, CEED_MEM_HOST, 208 CEED_USE_POINTER, sizeof(*setup_ctx), setup_ctx); 209 CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->setup_context); 210 CeedQFunctionContextCreate(ceed, &ceed_data->newt_ig_context); 211 CeedQFunctionContextSetData(ceed_data->newt_ig_context, CEED_MEM_HOST, 212 CEED_USE_POINTER, 213 sizeof(*phys->newtonian_ig_ctx), phys->newtonian_ig_ctx); 214 if (ceed_data->qf_rhs_vol) 215 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->newt_ig_context); 216 if (ceed_data->qf_ifunction_vol) 217 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, 218 ceed_data->newt_ig_context); 219 PetscFunctionReturn(0); 220 } 221