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