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 16 #ifndef newtonian_context_struct 17 #define newtonian_context_struct 18 typedef struct NewtonianIdealGasContext_ *NewtonianIdealGasContext; 19 struct NewtonianIdealGasContext_ { 20 CeedScalar lambda; 21 CeedScalar mu; 22 CeedScalar k; 23 CeedScalar cv; 24 CeedScalar cp; 25 CeedScalar g[3]; 26 CeedScalar c_tau; 27 CeedScalar Ctau_t; 28 CeedScalar Ctau_v; 29 CeedScalar Ctau_C; 30 CeedScalar Ctau_M; 31 CeedScalar Ctau_E; 32 CeedScalar dt; 33 StabilizationType stabilization; 34 }; 35 #endif 36 37 PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *setup_ctx, 38 void *ctx) { 39 SetupContext setup_context = *(SetupContext *)setup_ctx; 40 User user = *(User *)ctx; 41 StabilizationType stab; 42 MPI_Comm comm = PETSC_COMM_WORLD; 43 PetscBool implicit; 44 PetscBool has_curr_time = PETSC_FALSE; 45 PetscInt ierr; 46 PetscFunctionBeginUser; 47 48 ierr = PetscCalloc1(1, &user->phys->newtonian_ig_ctx); CHKERRQ(ierr); 49 50 // ------------------------------------------------------ 51 // Setup Generic Newtonian IG Problem 52 // ------------------------------------------------------ 53 problem->dim = 3; 54 problem->q_data_size_vol = 10; 55 problem->q_data_size_sur = 4; 56 problem->setup_vol.qfunction = Setup; 57 problem->setup_vol.qfunction_loc = Setup_loc; 58 problem->ics.qfunction = ICsNewtonianIG; 59 problem->ics.qfunction_loc = ICsNewtonianIG_loc; 60 problem->setup_sur.qfunction = SetupBoundary; 61 problem->setup_sur.qfunction_loc = SetupBoundary_loc; 62 problem->apply_vol_rhs.qfunction = Newtonian; 63 problem->apply_vol_rhs.qfunction_loc = Newtonian_loc; 64 problem->apply_vol_ifunction.qfunction = IFunction_Newtonian; 65 problem->apply_vol_ifunction.qfunction_loc = IFunction_Newtonian_loc; 66 problem->setup_ctx = SetupContext_DENSITY_CURRENT; 67 problem->non_zero_time = PETSC_FALSE; 68 problem->print_info = PRINT_DENSITY_CURRENT; 69 70 // ------------------------------------------------------ 71 // Create the libCEED context 72 // ------------------------------------------------------ 73 CeedScalar cv = 717.; // J/(kg K) 74 CeedScalar cp = 1004.; // J/(kg K) 75 CeedScalar g[3] = {0, 0, -9.81}; // m/s^2 76 CeedScalar lambda = -2./3.; // - 77 CeedScalar mu = 1.8e-5; // Pa s, dynamic viscosity 78 CeedScalar k = 0.02638; // W/(m K) 79 CeedScalar c_tau = 0.5; // - 80 CeedScalar Ctau_t = 1.0; // - 81 CeedScalar Ctau_v = 36.0; // TODO make function of degree 82 CeedScalar Ctau_C = 1.0; // TODO make function of degree 83 CeedScalar Ctau_M = 1.0; // TODO make function of degree 84 CeedScalar Ctau_E = 1.0; // TODO make function of degree 85 PetscReal domain_min[3], domain_max[3], domain_size[3]; 86 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 87 for (int i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 88 89 // ------------------------------------------------------ 90 // Create the PETSc context 91 // ------------------------------------------------------ 92 PetscScalar meter = 1; // 1 meter in scaled length units 93 PetscScalar kilogram = 1; // 1 kilogram in scaled mass units 94 PetscScalar second = 1; // 1 second in scaled time units 95 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 96 PetscScalar W_per_m_K, Pascal, J_per_kg_K, m_per_squared_s; 97 98 // ------------------------------------------------------ 99 // Command line Options 100 // ------------------------------------------------------ 101 PetscOptionsBegin(comm, NULL, "Options for Newtonian Ideal Gas based problem", 102 NULL); 103 104 // -- Physics 105 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 106 NULL, cv, &cv, NULL); CHKERRQ(ierr); 107 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 108 NULL, cp, &cp, NULL); CHKERRQ(ierr); 109 ierr = PetscOptionsScalar("-lambda", 110 "Stokes hypothesis second viscosity coefficient", 111 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 112 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 113 NULL, mu, &mu, NULL); CHKERRQ(ierr); 114 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 115 NULL, k, &k, NULL); CHKERRQ(ierr); 116 117 PetscInt dim = problem->dim; 118 ierr = PetscOptionsRealArray("-g", "Gravitational acceleration", 119 NULL, g, &dim, NULL); CHKERRQ(ierr); 120 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 121 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 122 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 123 ierr = PetscOptionsScalar("-c_tau", "Stabilization constant", 124 NULL, c_tau, &c_tau, NULL); CHKERRQ(ierr); 125 ierr = PetscOptionsScalar("-Ctau_t", "Stabilization time constant", 126 NULL, Ctau_t, &Ctau_t, NULL); CHKERRQ(ierr); 127 ierr = PetscOptionsScalar("-Ctau_v", "Stabilization viscous constant", 128 NULL, Ctau_v, &Ctau_v, NULL); CHKERRQ(ierr); 129 ierr = PetscOptionsScalar("-Ctau_C", "Stabilization continuity constant", 130 NULL, Ctau_C, &Ctau_C, NULL); CHKERRQ(ierr); 131 ierr = PetscOptionsScalar("-Ctau_M", "Stabilization momentum constant", 132 NULL, Ctau_M, &Ctau_M, NULL); CHKERRQ(ierr); 133 ierr = PetscOptionsScalar("-Ctau_E", "Stabilization energy constant", 134 NULL, Ctau_E, &Ctau_E, NULL); CHKERRQ(ierr); 135 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 136 NULL, implicit=PETSC_FALSE, &implicit, NULL); 137 CHKERRQ(ierr); 138 139 // -- Units 140 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 141 NULL, meter, &meter, NULL); CHKERRQ(ierr); 142 meter = fabs(meter); 143 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 144 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 145 kilogram = fabs(kilogram); 146 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 147 NULL, second, &second, NULL); CHKERRQ(ierr); 148 second = fabs(second); 149 ierr = PetscOptionsScalar("-units_Kelvin", 150 "1 Kelvin in scaled temperature units", 151 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 152 Kelvin = fabs(Kelvin); 153 154 // -- Warnings 155 if (stab == STAB_SUPG && !implicit) { 156 ierr = PetscPrintf(comm, 157 "Warning! Use -stab supg only with -implicit\n"); 158 CHKERRQ(ierr); 159 } 160 PetscOptionsEnd(); 161 162 // ------------------------------------------------------ 163 // Set up the PETSc context 164 // ------------------------------------------------------ 165 // -- Define derived units 166 Pascal = kilogram / (meter * PetscSqr(second)); 167 J_per_kg_K = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 168 m_per_squared_s = meter / PetscSqr(second); 169 W_per_m_K = kilogram * meter / (pow(second,3) * Kelvin); 170 171 user->units->meter = meter; 172 user->units->kilogram = kilogram; 173 user->units->second = second; 174 user->units->Kelvin = Kelvin; 175 user->units->Pascal = Pascal; 176 user->units->J_per_kg_K = J_per_kg_K; 177 user->units->m_per_squared_s = m_per_squared_s; 178 user->units->W_per_m_K = W_per_m_K; 179 180 // ------------------------------------------------------ 181 // Set up the libCEED context 182 // ------------------------------------------------------ 183 // -- Scale variables to desired units 184 cv *= J_per_kg_K; 185 cp *= J_per_kg_K; 186 mu *= Pascal * second; 187 k *= W_per_m_K; 188 for (int i=0; i<3; i++) domain_size[i] *= meter; 189 for (int i=0; i<3; i++) g[i] *= m_per_squared_s; 190 problem->dm_scale = meter; 191 192 // -- Setup Context 193 setup_context->cv = cv; 194 setup_context->cp = cp; 195 setup_context->lx = domain_size[0]; 196 setup_context->ly = domain_size[1]; 197 setup_context->lz = domain_size[2]; 198 setup_context->time = 0; 199 ierr = PetscArraycpy(setup_context->g, g, 3); CHKERRQ(ierr); 200 201 // -- Solver Settings 202 user->phys->stab = stab; 203 user->phys->implicit = implicit; 204 user->phys->has_curr_time = has_curr_time; 205 206 // -- QFunction Context 207 user->phys->newtonian_ig_ctx->lambda = lambda; 208 user->phys->newtonian_ig_ctx->mu = mu; 209 user->phys->newtonian_ig_ctx->k = k; 210 user->phys->newtonian_ig_ctx->cv = cv; 211 user->phys->newtonian_ig_ctx->cp = cp; 212 user->phys->newtonian_ig_ctx->c_tau = c_tau; 213 user->phys->newtonian_ig_ctx->Ctau_t = Ctau_t; 214 user->phys->newtonian_ig_ctx->Ctau_v = Ctau_v; 215 user->phys->newtonian_ig_ctx->Ctau_C = Ctau_C; 216 user->phys->newtonian_ig_ctx->Ctau_M = Ctau_M; 217 user->phys->newtonian_ig_ctx->Ctau_E = Ctau_E; 218 user->phys->newtonian_ig_ctx->stabilization = stab; 219 ierr = PetscArraycpy(user->phys->newtonian_ig_ctx->g, g, 3); CHKERRQ(ierr); 220 221 PetscFunctionReturn(0); 222 } 223 224 PetscErrorCode SetupContext_NEWTONIAN_IG(Ceed ceed, CeedData ceed_data, 225 AppCtx app_ctx, SetupContext setup_ctx, Physics phys) { 226 PetscFunctionBeginUser; 227 CeedQFunctionContextCreate(ceed, &ceed_data->setup_context); 228 CeedQFunctionContextSetData(ceed_data->setup_context, CEED_MEM_HOST, 229 CEED_USE_POINTER, sizeof(*setup_ctx), setup_ctx); 230 CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->setup_context); 231 232 CeedQFunctionContextCreate(ceed, &ceed_data->newt_ig_context); 233 CeedQFunctionContextSetData(ceed_data->newt_ig_context, CEED_MEM_HOST, 234 CEED_USE_POINTER, 235 sizeof(*phys->newtonian_ig_ctx), phys->newtonian_ig_ctx); 236 CeedQFunctionContextRegisterDouble(ceed_data->newt_ig_context, "timestep size", 237 offsetof(struct NewtonianIdealGasContext_, dt), 1, "Size of timestep, delta t"); 238 239 if (ceed_data->qf_rhs_vol) 240 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->newt_ig_context); 241 if (ceed_data->qf_ifunction_vol) 242 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, 243 ceed_data->newt_ig_context); 244 PetscFunctionReturn(0); 245 } 246