1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Utility functions for setting up Freestream boundary condition 6 7 #include "../qfunctions/bc_freestream.h" 8 9 #include <ceed.h> 10 #include <petscdm.h> 11 12 #include <navierstokes.h> 13 #include "../qfunctions/newtonian_types.h" 14 15 static const char *const RiemannSolverTypes[] = {"HLL", "HLLC", "RiemannSolverTypes", "RIEMANN_", NULL}; 16 17 static PetscErrorCode RiemannSolverUnitTests(NewtonianIdealGasContext gas, CeedScalar rtol); 18 19 PetscErrorCode FreestreamBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { 20 User user = *(User *)ctx; 21 MPI_Comm comm = user->comm; 22 Ceed ceed = user->ceed; 23 FreestreamContext freestream_ctx; 24 CeedQFunctionContext freestream_qfctx; 25 RiemannFluxType riemann = RIEMANN_HLLC; 26 PetscScalar meter = user->units->meter; 27 PetscScalar second = user->units->second; 28 PetscScalar Kelvin = user->units->Kelvin; 29 PetscScalar Pascal = user->units->Pascal; 30 31 PetscFunctionBeginUser; 32 // Freestream inherits reference state. We re-dimensionalize so the defaults 33 // in -help will be visible in SI units. 34 StatePrimitive Y_inf = {.pressure = reference->pressure / Pascal, .velocity = {0}, .temperature = reference->temperature / Kelvin}; 35 for (int i = 0; i < 3; i++) Y_inf.velocity[i] = reference->velocity[i] * second / meter; 36 37 PetscOptionsBegin(comm, NULL, "Options for Freestream boundary condition", NULL); 38 PetscCall(PetscOptionsEnum("-freestream_riemann", "Riemann solver to use in freestream boundary condition", NULL, RiemannSolverTypes, 39 (PetscEnum)riemann, (PetscEnum *)&riemann, NULL)); 40 PetscCall(PetscOptionsScalar("-freestream_pressure", "Pressure at freestream condition", NULL, Y_inf.pressure, &Y_inf.pressure, NULL)); 41 PetscInt narray = 3; 42 PetscCall(PetscOptionsScalarArray("-freestream_velocity", "Velocity at freestream condition", NULL, Y_inf.velocity, &narray, NULL)); 43 PetscCall(PetscOptionsScalar("-freestream_temperature", "Temperature at freestream condition", NULL, Y_inf.temperature, &Y_inf.temperature, NULL)); 44 PetscOptionsEnd(); 45 46 switch (user->phys->state_var) { 47 case STATEVAR_CONSERVATIVE: 48 switch (riemann) { 49 case RIEMANN_HLL: 50 problem->apply_freestream.qf_func_ptr = Freestream_Conserv_HLL; 51 problem->apply_freestream.qf_loc = Freestream_Conserv_HLL_loc; 52 problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Conserv_HLL; 53 problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Conserv_HLL_loc; 54 break; 55 case RIEMANN_HLLC: 56 problem->apply_freestream.qf_func_ptr = Freestream_Conserv_HLLC; 57 problem->apply_freestream.qf_loc = Freestream_Conserv_HLLC_loc; 58 problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Conserv_HLLC; 59 problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Conserv_HLLC_loc; 60 break; 61 } 62 break; 63 case STATEVAR_PRIMITIVE: 64 switch (riemann) { 65 case RIEMANN_HLL: 66 problem->apply_freestream.qf_func_ptr = Freestream_Prim_HLL; 67 problem->apply_freestream.qf_loc = Freestream_Prim_HLL_loc; 68 problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Prim_HLL; 69 problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Prim_HLL_loc; 70 break; 71 case RIEMANN_HLLC: 72 problem->apply_freestream.qf_func_ptr = Freestream_Prim_HLLC; 73 problem->apply_freestream.qf_loc = Freestream_Prim_HLLC_loc; 74 problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Prim_HLLC; 75 problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Prim_HLLC_loc; 76 break; 77 } 78 break; 79 case STATEVAR_ENTROPY: 80 switch (riemann) { 81 case RIEMANN_HLL: 82 problem->apply_freestream.qf_func_ptr = Freestream_Entropy_HLL; 83 problem->apply_freestream.qf_loc = Freestream_Entropy_HLL_loc; 84 problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Entropy_HLL; 85 problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Entropy_HLL_loc; 86 break; 87 case RIEMANN_HLLC: 88 problem->apply_freestream.qf_func_ptr = Freestream_Entropy_HLLC; 89 problem->apply_freestream.qf_loc = Freestream_Entropy_HLLC_loc; 90 problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Entropy_HLLC; 91 problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Entropy_HLLC_loc; 92 break; 93 } 94 break; 95 } 96 97 Y_inf.pressure *= Pascal; 98 for (int i = 0; i < 3; i++) Y_inf.velocity[i] *= meter / second; 99 Y_inf.temperature *= Kelvin; 100 101 State S_infty = StateFromPrimitive(newtonian_ig_ctx, Y_inf); 102 103 // -- Set freestream_ctx struct values 104 PetscCall(PetscCalloc1(1, &freestream_ctx)); 105 freestream_ctx->newtonian_ctx = *newtonian_ig_ctx; 106 freestream_ctx->S_infty = S_infty; 107 108 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &freestream_qfctx)); 109 PetscCallCeed(ceed, CeedQFunctionContextSetData(freestream_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*freestream_ctx), freestream_ctx)); 110 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(freestream_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 111 problem->apply_freestream.qfctx = freestream_qfctx; 112 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(freestream_qfctx, &problem->apply_freestream_jacobian.qfctx)); 113 114 { 115 PetscBool run_unit_tests = PETSC_FALSE; 116 117 PetscCall(PetscOptionsGetBool(NULL, NULL, "-riemann_solver_unit_tests", &run_unit_tests, NULL)); 118 if (run_unit_tests) PetscCall(RiemannSolverUnitTests(newtonian_ig_ctx, 5e-7)); 119 } 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 static const char *const OutflowTypes[] = {"RIEMANN", "PRESSURE", "OutflowType", "OUTFLOW_", NULL}; 124 typedef enum { 125 OUTFLOW_RIEMANN, 126 OUTFLOW_PRESSURE, 127 } OutflowType; 128 129 PetscErrorCode OutflowBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { 130 User user = *(User *)ctx; 131 Ceed ceed = user->ceed; 132 OutflowContext outflow_ctx; 133 OutflowType outflow_type = OUTFLOW_RIEMANN; 134 CeedQFunctionContext outflow_qfctx; 135 const PetscScalar Kelvin = user->units->Kelvin; 136 const PetscScalar Pascal = user->units->Pascal; 137 138 PetscFunctionBeginUser; 139 CeedScalar pressure = reference->pressure / Pascal; 140 CeedScalar temperature = reference->temperature / Kelvin; 141 CeedScalar recirc = 1, softplus_velocity = 1e-2; 142 PetscOptionsBegin(user->comm, NULL, "Options for Outflow boundary condition", NULL); 143 PetscCall( 144 PetscOptionsEnum("-outflow_type", "Type of outflow condition", NULL, OutflowTypes, (PetscEnum)outflow_type, (PetscEnum *)&outflow_type, NULL)); 145 PetscCall(PetscOptionsScalar("-outflow_pressure", "Pressure at outflow condition", NULL, pressure, &pressure, NULL)); 146 if (outflow_type == OUTFLOW_RIEMANN) { 147 PetscCall(PetscOptionsScalar("-outflow_temperature", "Temperature at outflow condition", NULL, temperature, &temperature, NULL)); 148 PetscCall( 149 PetscOptionsReal("-outflow_recirc", "Fraction of recirculation to allow in exterior velocity state [0,1]", NULL, recirc, &recirc, NULL)); 150 PetscCall(PetscOptionsReal("-outflow_softplus_velocity", "Characteristic velocity of softplus regularization", NULL, softplus_velocity, 151 &softplus_velocity, NULL)); 152 } 153 PetscOptionsEnd(); 154 pressure *= Pascal; 155 temperature *= Kelvin; 156 157 switch (outflow_type) { 158 case OUTFLOW_RIEMANN: 159 switch (user->phys->state_var) { 160 case STATEVAR_CONSERVATIVE: 161 problem->apply_outflow.qf_func_ptr = RiemannOutflow_Conserv; 162 problem->apply_outflow.qf_loc = RiemannOutflow_Conserv_loc; 163 problem->apply_outflow_jacobian.qf_func_ptr = RiemannOutflow_Jacobian_Conserv; 164 problem->apply_outflow_jacobian.qf_loc = RiemannOutflow_Jacobian_Conserv_loc; 165 break; 166 case STATEVAR_PRIMITIVE: 167 problem->apply_outflow.qf_func_ptr = RiemannOutflow_Prim; 168 problem->apply_outflow.qf_loc = RiemannOutflow_Prim_loc; 169 problem->apply_outflow_jacobian.qf_func_ptr = RiemannOutflow_Jacobian_Prim; 170 problem->apply_outflow_jacobian.qf_loc = RiemannOutflow_Jacobian_Prim_loc; 171 break; 172 case STATEVAR_ENTROPY: 173 problem->apply_outflow.qf_func_ptr = RiemannOutflow_Entropy; 174 problem->apply_outflow.qf_loc = RiemannOutflow_Entropy_loc; 175 problem->apply_outflow_jacobian.qf_func_ptr = RiemannOutflow_Jacobian_Entropy; 176 problem->apply_outflow_jacobian.qf_loc = RiemannOutflow_Jacobian_Entropy_loc; 177 break; 178 } 179 break; 180 case OUTFLOW_PRESSURE: 181 switch (user->phys->state_var) { 182 case STATEVAR_CONSERVATIVE: 183 problem->apply_outflow.qf_func_ptr = PressureOutflow_Conserv; 184 problem->apply_outflow.qf_loc = PressureOutflow_Conserv_loc; 185 problem->apply_outflow_jacobian.qf_func_ptr = PressureOutflow_Jacobian_Conserv; 186 problem->apply_outflow_jacobian.qf_loc = PressureOutflow_Jacobian_Conserv_loc; 187 break; 188 case STATEVAR_PRIMITIVE: 189 problem->apply_outflow.qf_func_ptr = PressureOutflow_Prim; 190 problem->apply_outflow.qf_loc = PressureOutflow_Prim_loc; 191 problem->apply_outflow_jacobian.qf_func_ptr = PressureOutflow_Jacobian_Prim; 192 problem->apply_outflow_jacobian.qf_loc = PressureOutflow_Jacobian_Prim_loc; 193 break; 194 case STATEVAR_ENTROPY: 195 problem->apply_outflow.qf_func_ptr = PressureOutflow_Entropy; 196 problem->apply_outflow.qf_loc = PressureOutflow_Entropy_loc; 197 problem->apply_outflow_jacobian.qf_func_ptr = PressureOutflow_Jacobian_Entropy; 198 problem->apply_outflow_jacobian.qf_loc = PressureOutflow_Jacobian_Entropy_loc; 199 break; 200 } 201 break; 202 } 203 PetscCall(PetscCalloc1(1, &outflow_ctx)); 204 outflow_ctx->gas = *newtonian_ig_ctx; 205 outflow_ctx->recirc = recirc; 206 outflow_ctx->softplus_velocity = softplus_velocity; 207 outflow_ctx->pressure = pressure; 208 outflow_ctx->temperature = temperature; 209 210 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &outflow_qfctx)); 211 PetscCallCeed(ceed, CeedQFunctionContextSetData(outflow_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*outflow_ctx), outflow_ctx)); 212 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(outflow_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 213 problem->apply_outflow.qfctx = outflow_qfctx; 214 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(outflow_qfctx, &problem->apply_outflow_jacobian.qfctx)); 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217 218 // @brief Calculate relative error, (A - B) / S 219 // If S < threshold, then set S=1 220 static inline CeedScalar RelativeError(CeedScalar S, CeedScalar A, CeedScalar B, CeedScalar threshold) { 221 return (A - B) / (fabs(S) > threshold ? S : 1); 222 } 223 224 // @brief Check errors of a State vector and print if above tolerance 225 static PetscErrorCode CheckQWithTolerance(const CeedScalar Q_s[5], const CeedScalar Q_a[5], const CeedScalar Q_b[5], const char *name, 226 PetscReal rtol_0, PetscReal rtol_u, PetscReal rtol_4) { 227 CeedScalar relative_error[5]; // relative error 228 CeedScalar divisor_threshold = 10 * CEED_EPSILON; 229 230 PetscFunctionBeginUser; 231 relative_error[0] = RelativeError(Q_s[0], Q_a[0], Q_b[0], divisor_threshold); 232 relative_error[4] = RelativeError(Q_s[4], Q_a[4], Q_b[4], divisor_threshold); 233 234 CeedScalar u_magnitude = sqrt(Square(Q_s[1]) + Square(Q_s[2]) + Square(Q_s[3])); 235 for (int i = 1; i < 4; i++) { 236 relative_error[i] = RelativeError(u_magnitude, Q_a[i], Q_b[i], divisor_threshold); 237 } 238 239 if (fabs(relative_error[0]) >= rtol_0) { 240 printf("%s[0] error %g (expected %.10e, got %.10e)\n", name, relative_error[0], Q_s[0], Q_a[0]); 241 } 242 for (int i = 1; i < 4; i++) { 243 if (fabs(relative_error[i]) >= rtol_u) { 244 printf("%s[%d] error %g (expected %.10e, got %.10e)\n", name, i, relative_error[i], Q_s[i], Q_a[i]); 245 } 246 } 247 if (fabs(relative_error[4]) >= rtol_4) { 248 printf("%s[4] error %g (expected %.10e, got %.10e)\n", name, relative_error[4], Q_s[4], Q_a[4]); 249 } 250 PetscFunctionReturn(PETSC_SUCCESS); 251 } 252 253 // @brief Verify RiemannFlux_HLL_fwd function against finite-difference approximation 254 static PetscErrorCode TestRiemannHLL_fwd(NewtonianIdealGasContext gas, CeedScalar rtol_0, CeedScalar rtol_u, CeedScalar rtol_4) { 255 CeedScalar eps = 4e-7; // Finite difference step 256 char buf[128]; 257 const CeedScalar T = 200; 258 const CeedScalar rho = 1.2; 259 const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; 260 const CeedScalar u_base = 40; 261 const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; 262 const CeedScalar Y0_left[5] = {p, u[0], u[1], u[2], T}; 263 const CeedScalar Y0_right[5] = {1.2 * p, 1.2 * u[0], 1.2 * u[1], 1.2 * u[2], 1.2 * T}; 264 CeedScalar normal[3] = {1, 2, 3}; 265 266 PetscFunctionBeginUser; 267 State left0 = StateFromY(gas, Y0_left); 268 State right0 = StateFromY(gas, Y0_right); 269 ScaleN(normal, 1 / Norm3(normal), 3); 270 271 for (int i = 0; i < 10; i++) { 272 CeedScalar dFlux[5] = {0.}, dFlux_fd[5] = {0.}; 273 { // Calculate dFlux using *_fwd function 274 CeedScalar dY_right[5] = {0}; 275 CeedScalar dY_left[5] = {0}; 276 277 if (i < 5) { 278 dY_left[i] = Y0_left[i]; 279 } else { 280 dY_right[i % 5] = Y0_right[i % 5]; 281 } 282 State dleft0 = StateFromY_fwd(gas, left0, dY_left); 283 State dright0 = StateFromY_fwd(gas, right0, dY_right); 284 285 StateConservative dFlux_state = RiemannFlux_HLL_fwd(gas, left0, dleft0, right0, dright0, normal); 286 UnpackState_U(dFlux_state, dFlux); 287 } 288 289 { // Calculate dFlux_fd via finite difference approximation 290 CeedScalar Y1_left[5] = {Y0_left[0], Y0_left[1], Y0_left[2], Y0_left[3], Y0_left[4]}; 291 CeedScalar Y1_right[5] = {Y0_right[0], Y0_right[1], Y0_right[2], Y0_right[3], Y0_right[4]}; 292 CeedScalar Flux0[5], Flux1[5]; 293 294 if (i < 5) { 295 Y1_left[i] *= 1 + eps; 296 } else { 297 Y1_right[i % 5] *= 1 + eps; 298 } 299 State left1 = StateFromY(gas, Y1_left); 300 State right1 = StateFromY(gas, Y1_right); 301 302 StateConservative Flux0_state = RiemannFlux_HLL(gas, left0, right0, normal); 303 StateConservative Flux1_state = RiemannFlux_HLL(gas, left1, right1, normal); 304 UnpackState_U(Flux0_state, Flux0); 305 UnpackState_U(Flux1_state, Flux1); 306 for (int j = 0; j < 5; j++) dFlux_fd[j] = (Flux1[j] - Flux0[j]) / eps; 307 } 308 309 snprintf(buf, sizeof buf, "RiemannFlux_HLL i=%d: Flux", i); 310 PetscCall(CheckQWithTolerance(dFlux_fd, dFlux, dFlux_fd, buf, rtol_0, rtol_u, rtol_4)); 311 } 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 315 // @brief Verify RiemannFlux_HLLC_fwd function against finite-difference approximation 316 static PetscErrorCode TestRiemannHLLC_fwd(NewtonianIdealGasContext gas, CeedScalar rtol_0, CeedScalar rtol_u, CeedScalar rtol_4) { 317 CeedScalar eps = 4e-7; // Finite difference step 318 char buf[128]; 319 const CeedScalar T = 200; 320 const CeedScalar rho = 1.2; 321 const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; 322 const CeedScalar u_base = 40; 323 const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; 324 const CeedScalar Y0_left[5] = {p, u[0], u[1], u[2], T}; 325 const CeedScalar Y0_right[5] = {1.2 * p, 1.2 * u[0], 1.2 * u[1], 1.2 * u[2], 1.2 * T}; 326 CeedScalar normal[3] = {1, 2, 3}; 327 328 PetscFunctionBeginUser; 329 State left0 = StateFromY(gas, Y0_left); 330 State right0 = StateFromY(gas, Y0_right); 331 ScaleN(normal, 1 / Norm3(normal), 3); 332 333 for (int i = 0; i < 10; i++) { 334 CeedScalar dFlux[5] = {0.}, dFlux_fd[5] = {0.}; 335 { // Calculate dFlux using *_fwd function 336 CeedScalar dY_right[5] = {0}; 337 CeedScalar dY_left[5] = {0}; 338 339 if (i < 5) { 340 dY_left[i] = Y0_left[i]; 341 } else { 342 dY_right[i % 5] = Y0_right[i % 5]; 343 } 344 State dleft0 = StateFromY_fwd(gas, left0, dY_left); 345 State dright0 = StateFromY_fwd(gas, right0, dY_right); 346 347 StateConservative dFlux_state = RiemannFlux_HLLC_fwd(gas, left0, dleft0, right0, dright0, normal); 348 UnpackState_U(dFlux_state, dFlux); 349 } 350 351 { // Calculate dFlux_fd via finite difference approximation 352 CeedScalar Y1_left[5] = {Y0_left[0], Y0_left[1], Y0_left[2], Y0_left[3], Y0_left[4]}; 353 CeedScalar Y1_right[5] = {Y0_right[0], Y0_right[1], Y0_right[2], Y0_right[3], Y0_right[4]}; 354 CeedScalar Flux0[5], Flux1[5]; 355 356 if (i < 5) { 357 Y1_left[i] *= 1 + eps; 358 } else { 359 Y1_right[i % 5] *= 1 + eps; 360 } 361 State left1 = StateFromY(gas, Y1_left); 362 State right1 = StateFromY(gas, Y1_right); 363 364 StateConservative Flux0_state = RiemannFlux_HLLC(gas, left0, right0, normal); 365 StateConservative Flux1_state = RiemannFlux_HLLC(gas, left1, right1, normal); 366 UnpackState_U(Flux0_state, Flux0); 367 UnpackState_U(Flux1_state, Flux1); 368 for (int j = 0; j < 5; j++) dFlux_fd[j] = (Flux1[j] - Flux0[j]) / eps; 369 } 370 371 snprintf(buf, sizeof buf, "RiemannFlux_HLLC i=%d: Flux", i); 372 PetscCall(CheckQWithTolerance(dFlux_fd, dFlux, dFlux_fd, buf, rtol_0, rtol_u, rtol_4)); 373 } 374 PetscFunctionReturn(PETSC_SUCCESS); 375 } 376 377 // @brief Verify ComputeHLLSpeeds_Roe_fwd function against finite-difference approximation 378 static PetscErrorCode TestComputeHLLSpeeds_Roe_fwd(NewtonianIdealGasContext gas, CeedScalar rtol) { 379 CeedScalar eps = 4e-7; // Finite difference step 380 char buf[128]; 381 const CeedScalar T = 200; 382 const CeedScalar rho = 1.2; 383 const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; 384 const CeedScalar u_base = 40; 385 const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; 386 const CeedScalar Y0_left[5] = {p, u[0], u[1], u[2], T}; 387 const CeedScalar Y0_right[5] = {1.2 * p, 1.2 * u[0], 1.2 * u[1], 1.2 * u[2], 1.2 * T}; 388 CeedScalar normal[3] = {1, 2, 3}; 389 390 PetscFunctionBeginUser; 391 State left0 = StateFromY(gas, Y0_left); 392 State right0 = StateFromY(gas, Y0_right); 393 ScaleN(normal, 1 / Norm3(normal), 3); 394 CeedScalar u_left0 = Dot3(left0.Y.velocity, normal); 395 CeedScalar u_right0 = Dot3(right0.Y.velocity, normal); 396 397 for (int i = 0; i < 10; i++) { 398 CeedScalar ds_left, ds_right, ds_left_fd, ds_right_fd; 399 { // Calculate ds_{left,right} using *_fwd function 400 CeedScalar dY_right[5] = {0}; 401 CeedScalar dY_left[5] = {0}; 402 403 if (i < 5) { 404 dY_left[i] = Y0_left[i]; 405 } else { 406 dY_right[i % 5] = Y0_right[i % 5]; 407 } 408 State dleft0 = StateFromY_fwd(gas, left0, dY_left); 409 State dright0 = StateFromY_fwd(gas, right0, dY_right); 410 CeedScalar du_left = Dot3(dleft0.Y.velocity, normal); 411 CeedScalar du_right = Dot3(dright0.Y.velocity, normal); 412 413 CeedScalar s_left, s_right; // Throw away 414 ComputeHLLSpeeds_Roe_fwd(gas, left0, dleft0, u_left0, du_left, right0, dright0, u_right0, du_right, &s_left, &ds_left, &s_right, &ds_right); 415 } 416 417 { // Calculate ds_{left,right}_fd via finite difference approximation 418 CeedScalar Y1_left[5] = {Y0_left[0], Y0_left[1], Y0_left[2], Y0_left[3], Y0_left[4]}; 419 CeedScalar Y1_right[5] = {Y0_right[0], Y0_right[1], Y0_right[2], Y0_right[3], Y0_right[4]}; 420 421 if (i < 5) { 422 Y1_left[i] *= 1 + eps; 423 } else { 424 Y1_right[i % 5] *= 1 + eps; 425 } 426 State left1 = StateFromY(gas, Y1_left); 427 State right1 = StateFromY(gas, Y1_right); 428 CeedScalar u_left1 = Dot3(left1.Y.velocity, normal); 429 CeedScalar u_right1 = Dot3(right1.Y.velocity, normal); 430 431 CeedScalar s_left0, s_right0, s_left1, s_right1; 432 ComputeHLLSpeeds_Roe(gas, left0, u_left0, right0, u_right0, &s_left0, &s_right0); 433 ComputeHLLSpeeds_Roe(gas, left1, u_left1, right1, u_right1, &s_left1, &s_right1); 434 ds_left_fd = (s_left1 - s_left0) / eps; 435 ds_right_fd = (s_right1 - s_right0) / eps; 436 } 437 438 snprintf(buf, sizeof buf, "ComputeHLLSpeeds_Roe i=%d:", i); 439 { 440 CeedScalar divisor_threshold = 10 * CEED_EPSILON; 441 CeedScalar ds_left_err, ds_right_err; 442 443 ds_left_err = RelativeError(ds_left_fd, ds_left, ds_left_fd, divisor_threshold); 444 ds_right_err = RelativeError(ds_right_fd, ds_right, ds_right_fd, divisor_threshold); 445 if (fabs(ds_left_err) >= rtol) printf("%s ds_left error %g (expected %.10e, got %.10e)\n", buf, ds_left_err, ds_left_fd, ds_left); 446 if (fabs(ds_right_err) >= rtol) printf("%s ds_right error %g (expected %.10e, got %.10e)\n", buf, ds_right_err, ds_right_fd, ds_right); 447 } 448 } 449 PetscFunctionReturn(PETSC_SUCCESS); 450 } 451 452 // @brief Verify TotalSpecificEnthalpy_fwd function against finite-difference approximation 453 static PetscErrorCode TestTotalSpecificEnthalpy_fwd(NewtonianIdealGasContext gas, CeedScalar rtol) { 454 CeedScalar eps = 4e-7; // Finite difference step 455 char buf[128]; 456 const CeedScalar T = 200; 457 const CeedScalar rho = 1.2; 458 const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; 459 const CeedScalar u_base = 40; 460 const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; 461 const CeedScalar Y0[5] = {p, u[0], u[1], u[2], T}; 462 463 PetscFunctionBeginUser; 464 State state0 = StateFromY(gas, Y0); 465 466 for (int i = 0; i < 5; i++) { 467 CeedScalar dH, dH_fd; 468 { // Calculate dH using *_fwd function 469 CeedScalar dY[5] = {0}; 470 471 dY[i] = Y0[i]; 472 State dstate0 = StateFromY_fwd(gas, state0, dY); 473 dH = TotalSpecificEnthalpy_fwd(gas, state0, dstate0); 474 } 475 476 { // Calculate dH_fd via finite difference approximation 477 CeedScalar H0, H1; 478 CeedScalar Y1[5] = {Y0[0], Y0[1], Y0[2], Y0[3], Y0[4]}; 479 Y1[i] *= 1 + eps; 480 State state1 = StateFromY(gas, Y1); 481 482 H0 = TotalSpecificEnthalpy(gas, state0); 483 H1 = TotalSpecificEnthalpy(gas, state1); 484 dH_fd = (H1 - H0) / eps; 485 } 486 487 snprintf(buf, sizeof buf, "TotalSpecificEnthalpy i=%d:", i); 488 { 489 CeedScalar divisor_threshold = 10 * CEED_EPSILON; 490 CeedScalar dH_err; 491 492 dH_err = RelativeError(dH_fd, dH, dH_fd, divisor_threshold); 493 if (fabs(dH_err) >= rtol) printf("%s dH error %g (expected %.10e, got %.10e)\n", buf, dH_err, dH_fd, dH); 494 } 495 } 496 PetscFunctionReturn(PETSC_SUCCESS); 497 } 498 499 // @brief Verify RoeSetup_fwd function against finite-difference approximation 500 static PetscErrorCode TestRowSetup_fwd(NewtonianIdealGasContext gas, CeedScalar rtol) { 501 CeedScalar eps = 4e-7; // Finite difference step 502 char buf[128]; 503 const CeedScalar rho0[2] = {1.2, 1.4}; 504 505 PetscFunctionBeginUser; 506 for (int i = 0; i < 2; i++) { 507 RoeWeights dR, dR_fd; 508 { // Calculate using *_fwd function 509 CeedScalar drho[5] = {0}; 510 511 drho[i] = rho0[i]; 512 dR = RoeSetup_fwd(rho0[0], rho0[1], drho[0], drho[1]); 513 } 514 515 { // Calculate via finite difference approximation 516 RoeWeights R0, R1; 517 CeedScalar rho1[5] = {rho0[0], rho0[1]}; 518 rho1[i] *= 1 + eps; 519 520 R0 = RoeSetup(rho0[0], rho0[1]); 521 R1 = RoeSetup(rho1[0], rho1[1]); 522 dR_fd.left = (R1.left - R0.left) / eps; 523 dR_fd.right = (R1.right - R0.right) / eps; 524 } 525 526 snprintf(buf, sizeof buf, "RoeSetup i=%d:", i); 527 { 528 CeedScalar divisor_threshold = 10 * CEED_EPSILON; 529 RoeWeights dR_err; 530 531 dR_err.left = RelativeError(dR_fd.left, dR.left, dR_fd.left, divisor_threshold); 532 dR_err.right = RelativeError(dR_fd.right, dR.right, dR_fd.right, divisor_threshold); 533 if (fabs(dR_err.left) >= rtol) printf("%s dR.left error %g (expected %.10e, got %.10e)\n", buf, dR_err.left, dR_fd.left, dR.left); 534 if (fabs(dR_err.right) >= rtol) printf("%s dR.right error %g (expected %.10e, got %.10e)\n", buf, dR_err.right, dR_fd.right, dR.right); 535 } 536 } 537 PetscFunctionReturn(PETSC_SUCCESS); 538 } 539 540 // @brief Test Riemann solver related `*_fwd` functions via finite-difference approximation 541 static PetscErrorCode RiemannSolverUnitTests(NewtonianIdealGasContext gas, CeedScalar rtol) { 542 PetscFunctionBeginUser; 543 PetscCall(TestRiemannHLL_fwd(gas, rtol, rtol, rtol)); 544 PetscCall(TestRiemannHLLC_fwd(gas, rtol, rtol, rtol)); 545 PetscCall(TestComputeHLLSpeeds_Roe_fwd(gas, rtol)); 546 PetscCall(TestTotalSpecificEnthalpy_fwd(gas, rtol)); 547 PetscCall(TestRowSetup_fwd(gas, rtol)); 548 PetscFunctionReturn(PETSC_SUCCESS); 549 } 550