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