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