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