// SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause /// @file /// Utility functions for setting up Freestream boundary condition #include "../qfunctions/bc_freestream.h" #include #include #include #include "../qfunctions/newtonian_types.h" static const char *const RiemannSolverTypes[] = {"HLL", "HLLC", "RiemannSolverTypes", "RIEMANN_", NULL}; static PetscErrorCode RiemannSolverUnitTests(NewtonianIdealGasContext gas, CeedScalar rtol); PetscErrorCode FreestreamBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { Honee honee = *(Honee *)ctx; MPI_Comm comm = honee->comm; Ceed ceed = honee->ceed; FreestreamContext freestream_ctx; CeedQFunctionContext freestream_qfctx; RiemannFluxType riemann = RIEMANN_HLLC; PetscScalar meter = honee->units->meter; PetscScalar second = honee->units->second; PetscScalar Kelvin = honee->units->Kelvin; PetscScalar Pascal = honee->units->Pascal; PetscFunctionBeginUser; // Freestream inherits reference state. We re-dimensionalize so the defaults // in -help will be visible in SI units. StatePrimitive Y_inf = {.pressure = reference->pressure / Pascal, .velocity = {0}, .temperature = reference->temperature / Kelvin}; for (int i = 0; i < 3; i++) Y_inf.velocity[i] = reference->velocity[i] * second / meter; PetscOptionsBegin(comm, NULL, "Options for Freestream boundary condition", NULL); PetscCall(PetscOptionsEnum("-freestream_riemann", "Riemann solver to use in freestream boundary condition", NULL, RiemannSolverTypes, (PetscEnum)riemann, (PetscEnum *)&riemann, NULL)); PetscCall(PetscOptionsScalar("-freestream_pressure", "Pressure at freestream condition", NULL, Y_inf.pressure, &Y_inf.pressure, NULL)); PetscInt narray = 3; PetscCall(PetscOptionsScalarArray("-freestream_velocity", "Velocity at freestream condition", NULL, Y_inf.velocity, &narray, NULL)); PetscCall(PetscOptionsScalar("-freestream_temperature", "Temperature at freestream condition", NULL, Y_inf.temperature, &Y_inf.temperature, NULL)); PetscOptionsEnd(); switch (honee->phys->state_var) { case STATEVAR_CONSERVATIVE: switch (riemann) { case RIEMANN_HLL: problem->apply_freestream.qf_func_ptr = Freestream_Conserv_HLL; problem->apply_freestream.qf_loc = Freestream_Conserv_HLL_loc; problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Conserv_HLL; problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Conserv_HLL_loc; break; case RIEMANN_HLLC: problem->apply_freestream.qf_func_ptr = Freestream_Conserv_HLLC; problem->apply_freestream.qf_loc = Freestream_Conserv_HLLC_loc; problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Conserv_HLLC; problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Conserv_HLLC_loc; break; } break; case STATEVAR_PRIMITIVE: switch (riemann) { case RIEMANN_HLL: problem->apply_freestream.qf_func_ptr = Freestream_Prim_HLL; problem->apply_freestream.qf_loc = Freestream_Prim_HLL_loc; problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Prim_HLL; problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Prim_HLL_loc; break; case RIEMANN_HLLC: problem->apply_freestream.qf_func_ptr = Freestream_Prim_HLLC; problem->apply_freestream.qf_loc = Freestream_Prim_HLLC_loc; problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Prim_HLLC; problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Prim_HLLC_loc; break; } break; case STATEVAR_ENTROPY: switch (riemann) { case RIEMANN_HLL: problem->apply_freestream.qf_func_ptr = Freestream_Entropy_HLL; problem->apply_freestream.qf_loc = Freestream_Entropy_HLL_loc; problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Entropy_HLL; problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Entropy_HLL_loc; break; case RIEMANN_HLLC: problem->apply_freestream.qf_func_ptr = Freestream_Entropy_HLLC; problem->apply_freestream.qf_loc = Freestream_Entropy_HLLC_loc; problem->apply_freestream_jacobian.qf_func_ptr = Freestream_Jacobian_Entropy_HLLC; problem->apply_freestream_jacobian.qf_loc = Freestream_Jacobian_Entropy_HLLC_loc; break; } break; } Y_inf.pressure *= Pascal; for (int i = 0; i < 3; i++) Y_inf.velocity[i] *= meter / second; Y_inf.temperature *= Kelvin; State S_infty = StateFromPrimitive(newtonian_ig_ctx, Y_inf); // -- Set freestream_ctx struct values PetscCall(PetscCalloc1(1, &freestream_ctx)); freestream_ctx->newtonian_ctx = *newtonian_ig_ctx; freestream_ctx->S_infty = S_infty; PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &freestream_qfctx)); PetscCallCeed(ceed, CeedQFunctionContextSetData(freestream_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*freestream_ctx), freestream_ctx)); PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(freestream_qfctx, CEED_MEM_HOST, FreeContextPetsc)); problem->apply_freestream.qfctx = freestream_qfctx; PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(freestream_qfctx, &problem->apply_freestream_jacobian.qfctx)); { PetscBool run_unit_tests = PETSC_FALSE; PetscCall(PetscOptionsGetBool(NULL, NULL, "-riemann_solver_unit_tests", &run_unit_tests, NULL)); if (run_unit_tests) PetscCall(RiemannSolverUnitTests(newtonian_ig_ctx, 5e-7)); } PetscFunctionReturn(PETSC_SUCCESS); } static const char *const OutflowTypes[] = {"RIEMANN", "PRESSURE", "OutflowType", "OUTFLOW_", NULL}; typedef enum { OUTFLOW_RIEMANN, OUTFLOW_PRESSURE, } OutflowType; PetscErrorCode OutflowBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { Honee honee = *(Honee *)ctx; Ceed ceed = honee->ceed; OutflowContext outflow_ctx; OutflowType outflow_type = OUTFLOW_RIEMANN; CeedQFunctionContext outflow_qfctx; const PetscScalar Kelvin = honee->units->Kelvin; const PetscScalar Pascal = honee->units->Pascal; PetscFunctionBeginUser; CeedScalar pressure = reference->pressure / Pascal; CeedScalar temperature = reference->temperature / Kelvin; CeedScalar recirc = 1, softplus_velocity = 1e-2; PetscOptionsBegin(honee->comm, NULL, "Options for Outflow boundary condition", NULL); PetscCall( PetscOptionsEnum("-outflow_type", "Type of outflow condition", NULL, OutflowTypes, (PetscEnum)outflow_type, (PetscEnum *)&outflow_type, NULL)); PetscCall(PetscOptionsScalar("-outflow_pressure", "Pressure at outflow condition", NULL, pressure, &pressure, NULL)); if (outflow_type == OUTFLOW_RIEMANN) { PetscCall(PetscOptionsScalar("-outflow_temperature", "Temperature at outflow condition", NULL, temperature, &temperature, NULL)); PetscCall( PetscOptionsReal("-outflow_recirc", "Fraction of recirculation to allow in exterior velocity state [0,1]", NULL, recirc, &recirc, NULL)); PetscCall(PetscOptionsReal("-outflow_softplus_velocity", "Characteristic velocity of softplus regularization", NULL, softplus_velocity, &softplus_velocity, NULL)); } PetscOptionsEnd(); pressure *= Pascal; temperature *= Kelvin; switch (outflow_type) { case OUTFLOW_RIEMANN: switch (honee->phys->state_var) { case STATEVAR_CONSERVATIVE: problem->apply_outflow.qf_func_ptr = RiemannOutflow_Conserv; problem->apply_outflow.qf_loc = RiemannOutflow_Conserv_loc; problem->apply_outflow_jacobian.qf_func_ptr = RiemannOutflow_Jacobian_Conserv; problem->apply_outflow_jacobian.qf_loc = RiemannOutflow_Jacobian_Conserv_loc; break; case STATEVAR_PRIMITIVE: problem->apply_outflow.qf_func_ptr = RiemannOutflow_Prim; problem->apply_outflow.qf_loc = RiemannOutflow_Prim_loc; problem->apply_outflow_jacobian.qf_func_ptr = RiemannOutflow_Jacobian_Prim; problem->apply_outflow_jacobian.qf_loc = RiemannOutflow_Jacobian_Prim_loc; break; case STATEVAR_ENTROPY: problem->apply_outflow.qf_func_ptr = RiemannOutflow_Entropy; problem->apply_outflow.qf_loc = RiemannOutflow_Entropy_loc; problem->apply_outflow_jacobian.qf_func_ptr = RiemannOutflow_Jacobian_Entropy; problem->apply_outflow_jacobian.qf_loc = RiemannOutflow_Jacobian_Entropy_loc; break; } break; case OUTFLOW_PRESSURE: switch (honee->phys->state_var) { case STATEVAR_CONSERVATIVE: problem->apply_outflow.qf_func_ptr = PressureOutflow_Conserv; problem->apply_outflow.qf_loc = PressureOutflow_Conserv_loc; problem->apply_outflow_jacobian.qf_func_ptr = PressureOutflow_Jacobian_Conserv; problem->apply_outflow_jacobian.qf_loc = PressureOutflow_Jacobian_Conserv_loc; break; case STATEVAR_PRIMITIVE: problem->apply_outflow.qf_func_ptr = PressureOutflow_Prim; problem->apply_outflow.qf_loc = PressureOutflow_Prim_loc; problem->apply_outflow_jacobian.qf_func_ptr = PressureOutflow_Jacobian_Prim; problem->apply_outflow_jacobian.qf_loc = PressureOutflow_Jacobian_Prim_loc; break; case STATEVAR_ENTROPY: problem->apply_outflow.qf_func_ptr = PressureOutflow_Entropy; problem->apply_outflow.qf_loc = PressureOutflow_Entropy_loc; problem->apply_outflow_jacobian.qf_func_ptr = PressureOutflow_Jacobian_Entropy; problem->apply_outflow_jacobian.qf_loc = PressureOutflow_Jacobian_Entropy_loc; break; } break; } PetscCall(PetscCalloc1(1, &outflow_ctx)); outflow_ctx->gas = *newtonian_ig_ctx; outflow_ctx->recirc = recirc; outflow_ctx->softplus_velocity = softplus_velocity; outflow_ctx->pressure = pressure; outflow_ctx->temperature = temperature; PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &outflow_qfctx)); PetscCallCeed(ceed, CeedQFunctionContextSetData(outflow_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*outflow_ctx), outflow_ctx)); PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(outflow_qfctx, CEED_MEM_HOST, FreeContextPetsc)); problem->apply_outflow.qfctx = outflow_qfctx; PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(outflow_qfctx, &problem->apply_outflow_jacobian.qfctx)); PetscFunctionReturn(PETSC_SUCCESS); } // @brief Calculate relative error, (A - B) / S // If S < threshold, then set S=1 static inline CeedScalar RelativeError(CeedScalar S, CeedScalar A, CeedScalar B, CeedScalar threshold) { return (A - B) / (fabs(S) > threshold ? S : 1); } // @brief Check errors of a State vector and print if above tolerance static PetscErrorCode CheckQWithTolerance(const CeedScalar Q_s[5], const CeedScalar Q_a[5], const CeedScalar Q_b[5], const char *name, PetscReal rtol_0, PetscReal rtol_u, PetscReal rtol_4) { CeedScalar relative_error[5]; // relative error CeedScalar divisor_threshold = 10 * CEED_EPSILON; PetscFunctionBeginUser; relative_error[0] = RelativeError(Q_s[0], Q_a[0], Q_b[0], divisor_threshold); relative_error[4] = RelativeError(Q_s[4], Q_a[4], Q_b[4], divisor_threshold); CeedScalar u_magnitude = sqrt(Square(Q_s[1]) + Square(Q_s[2]) + Square(Q_s[3])); for (int i = 1; i < 4; i++) { relative_error[i] = RelativeError(u_magnitude, Q_a[i], Q_b[i], divisor_threshold); } if (fabs(relative_error[0]) >= rtol_0) { printf("%s[0] error %g (expected %.10e, got %.10e)\n", name, relative_error[0], Q_s[0], Q_a[0]); } for (int i = 1; i < 4; i++) { if (fabs(relative_error[i]) >= rtol_u) { printf("%s[%d] error %g (expected %.10e, got %.10e)\n", name, i, relative_error[i], Q_s[i], Q_a[i]); } } if (fabs(relative_error[4]) >= rtol_4) { printf("%s[4] error %g (expected %.10e, got %.10e)\n", name, relative_error[4], Q_s[4], Q_a[4]); } PetscFunctionReturn(PETSC_SUCCESS); } // @brief Verify RiemannFlux_HLL_fwd function against finite-difference approximation static PetscErrorCode TestRiemannHLL_fwd(NewtonianIdealGasContext gas, CeedScalar rtol_0, CeedScalar rtol_u, CeedScalar rtol_4) { CeedScalar eps = 4e-7; // Finite difference step char buf[128]; const CeedScalar T = 200; const CeedScalar rho = 1.2; const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; const CeedScalar u_base = 40; const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; const CeedScalar Y0_left[5] = {p, u[0], u[1], u[2], T}; const CeedScalar Y0_right[5] = {1.2 * p, 1.2 * u[0], 1.2 * u[1], 1.2 * u[2], 1.2 * T}; CeedScalar normal[3] = {1, 2, 3}; PetscFunctionBeginUser; State left0 = StateFromY(gas, Y0_left); State right0 = StateFromY(gas, Y0_right); ScaleN(normal, 1 / Norm3(normal), 3); for (int i = 0; i < 10; i++) { CeedScalar dFlux[5] = {0.}, dFlux_fd[5] = {0.}; { // Calculate dFlux using *_fwd function CeedScalar dY_right[5] = {0}; CeedScalar dY_left[5] = {0}; if (i < 5) { dY_left[i] = Y0_left[i]; } else { dY_right[i % 5] = Y0_right[i % 5]; } State dleft0 = StateFromY_fwd(gas, left0, dY_left); State dright0 = StateFromY_fwd(gas, right0, dY_right); StateConservative dFlux_state = RiemannFlux_HLL_fwd(gas, left0, dleft0, right0, dright0, normal); UnpackState_U(dFlux_state, dFlux); } { // Calculate dFlux_fd via finite difference approximation CeedScalar Y1_left[5] = {Y0_left[0], Y0_left[1], Y0_left[2], Y0_left[3], Y0_left[4]}; CeedScalar Y1_right[5] = {Y0_right[0], Y0_right[1], Y0_right[2], Y0_right[3], Y0_right[4]}; CeedScalar Flux0[5], Flux1[5]; if (i < 5) { Y1_left[i] *= 1 + eps; } else { Y1_right[i % 5] *= 1 + eps; } State left1 = StateFromY(gas, Y1_left); State right1 = StateFromY(gas, Y1_right); StateConservative Flux0_state = RiemannFlux_HLL(gas, left0, right0, normal); StateConservative Flux1_state = RiemannFlux_HLL(gas, left1, right1, normal); UnpackState_U(Flux0_state, Flux0); UnpackState_U(Flux1_state, Flux1); for (int j = 0; j < 5; j++) dFlux_fd[j] = (Flux1[j] - Flux0[j]) / eps; } snprintf(buf, sizeof buf, "RiemannFlux_HLL i=%d: Flux", i); PetscCall(CheckQWithTolerance(dFlux_fd, dFlux, dFlux_fd, buf, rtol_0, rtol_u, rtol_4)); } PetscFunctionReturn(PETSC_SUCCESS); } // @brief Verify RiemannFlux_HLLC_fwd function against finite-difference approximation static PetscErrorCode TestRiemannHLLC_fwd(NewtonianIdealGasContext gas, CeedScalar rtol_0, CeedScalar rtol_u, CeedScalar rtol_4) { CeedScalar eps = 4e-7; // Finite difference step char buf[128]; const CeedScalar T = 200; const CeedScalar rho = 1.2; const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; const CeedScalar u_base = 40; const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; const CeedScalar Y0_left[5] = {p, u[0], u[1], u[2], T}; const CeedScalar Y0_right[5] = {1.2 * p, 1.2 * u[0], 1.2 * u[1], 1.2 * u[2], 1.2 * T}; CeedScalar normal[3] = {1, 2, 3}; PetscFunctionBeginUser; State left0 = StateFromY(gas, Y0_left); State right0 = StateFromY(gas, Y0_right); ScaleN(normal, 1 / Norm3(normal), 3); for (int i = 0; i < 10; i++) { CeedScalar dFlux[5] = {0.}, dFlux_fd[5] = {0.}; { // Calculate dFlux using *_fwd function CeedScalar dY_right[5] = {0}; CeedScalar dY_left[5] = {0}; if (i < 5) { dY_left[i] = Y0_left[i]; } else { dY_right[i % 5] = Y0_right[i % 5]; } State dleft0 = StateFromY_fwd(gas, left0, dY_left); State dright0 = StateFromY_fwd(gas, right0, dY_right); StateConservative dFlux_state = RiemannFlux_HLLC_fwd(gas, left0, dleft0, right0, dright0, normal); UnpackState_U(dFlux_state, dFlux); } { // Calculate dFlux_fd via finite difference approximation CeedScalar Y1_left[5] = {Y0_left[0], Y0_left[1], Y0_left[2], Y0_left[3], Y0_left[4]}; CeedScalar Y1_right[5] = {Y0_right[0], Y0_right[1], Y0_right[2], Y0_right[3], Y0_right[4]}; CeedScalar Flux0[5], Flux1[5]; if (i < 5) { Y1_left[i] *= 1 + eps; } else { Y1_right[i % 5] *= 1 + eps; } State left1 = StateFromY(gas, Y1_left); State right1 = StateFromY(gas, Y1_right); StateConservative Flux0_state = RiemannFlux_HLLC(gas, left0, right0, normal); StateConservative Flux1_state = RiemannFlux_HLLC(gas, left1, right1, normal); UnpackState_U(Flux0_state, Flux0); UnpackState_U(Flux1_state, Flux1); for (int j = 0; j < 5; j++) dFlux_fd[j] = (Flux1[j] - Flux0[j]) / eps; } snprintf(buf, sizeof buf, "RiemannFlux_HLLC i=%d: Flux", i); PetscCall(CheckQWithTolerance(dFlux_fd, dFlux, dFlux_fd, buf, rtol_0, rtol_u, rtol_4)); } PetscFunctionReturn(PETSC_SUCCESS); } // @brief Verify ComputeHLLSpeeds_Roe_fwd function against finite-difference approximation static PetscErrorCode TestComputeHLLSpeeds_Roe_fwd(NewtonianIdealGasContext gas, CeedScalar rtol) { CeedScalar eps = 4e-7; // Finite difference step char buf[128]; const CeedScalar T = 200; const CeedScalar rho = 1.2; const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; const CeedScalar u_base = 40; const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; const CeedScalar Y0_left[5] = {p, u[0], u[1], u[2], T}; const CeedScalar Y0_right[5] = {1.2 * p, 1.2 * u[0], 1.2 * u[1], 1.2 * u[2], 1.2 * T}; CeedScalar normal[3] = {1, 2, 3}; PetscFunctionBeginUser; State left0 = StateFromY(gas, Y0_left); State right0 = StateFromY(gas, Y0_right); ScaleN(normal, 1 / Norm3(normal), 3); CeedScalar u_left0 = Dot3(left0.Y.velocity, normal); CeedScalar u_right0 = Dot3(right0.Y.velocity, normal); for (int i = 0; i < 10; i++) { CeedScalar ds_left, ds_right, ds_left_fd, ds_right_fd; { // Calculate ds_{left,right} using *_fwd function CeedScalar dY_right[5] = {0}; CeedScalar dY_left[5] = {0}; if (i < 5) { dY_left[i] = Y0_left[i]; } else { dY_right[i % 5] = Y0_right[i % 5]; } State dleft0 = StateFromY_fwd(gas, left0, dY_left); State dright0 = StateFromY_fwd(gas, right0, dY_right); CeedScalar du_left = Dot3(dleft0.Y.velocity, normal); CeedScalar du_right = Dot3(dright0.Y.velocity, normal); CeedScalar s_left, s_right; // Throw away ComputeHLLSpeeds_Roe_fwd(gas, left0, dleft0, u_left0, du_left, right0, dright0, u_right0, du_right, &s_left, &ds_left, &s_right, &ds_right); } { // Calculate ds_{left,right}_fd via finite difference approximation CeedScalar Y1_left[5] = {Y0_left[0], Y0_left[1], Y0_left[2], Y0_left[3], Y0_left[4]}; CeedScalar Y1_right[5] = {Y0_right[0], Y0_right[1], Y0_right[2], Y0_right[3], Y0_right[4]}; if (i < 5) { Y1_left[i] *= 1 + eps; } else { Y1_right[i % 5] *= 1 + eps; } State left1 = StateFromY(gas, Y1_left); State right1 = StateFromY(gas, Y1_right); CeedScalar u_left1 = Dot3(left1.Y.velocity, normal); CeedScalar u_right1 = Dot3(right1.Y.velocity, normal); CeedScalar s_left0, s_right0, s_left1, s_right1; ComputeHLLSpeeds_Roe(gas, left0, u_left0, right0, u_right0, &s_left0, &s_right0); ComputeHLLSpeeds_Roe(gas, left1, u_left1, right1, u_right1, &s_left1, &s_right1); ds_left_fd = (s_left1 - s_left0) / eps; ds_right_fd = (s_right1 - s_right0) / eps; } snprintf(buf, sizeof buf, "ComputeHLLSpeeds_Roe i=%d:", i); { CeedScalar divisor_threshold = 10 * CEED_EPSILON; CeedScalar ds_left_err, ds_right_err; ds_left_err = RelativeError(ds_left_fd, ds_left, ds_left_fd, divisor_threshold); ds_right_err = RelativeError(ds_right_fd, ds_right, ds_right_fd, divisor_threshold); 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); 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); } } PetscFunctionReturn(PETSC_SUCCESS); } // @brief Verify TotalSpecificEnthalpy_fwd function against finite-difference approximation static PetscErrorCode TestTotalSpecificEnthalpy_fwd(NewtonianIdealGasContext gas, CeedScalar rtol) { CeedScalar eps = 4e-7; // Finite difference step char buf[128]; const CeedScalar T = 200; const CeedScalar rho = 1.2; const CeedScalar p = (HeatCapacityRatio(gas) - 1) * rho * gas->cv * T; const CeedScalar u_base = 40; const CeedScalar u[3] = {u_base, u_base * 1.1, u_base * 1.2}; const CeedScalar Y0[5] = {p, u[0], u[1], u[2], T}; PetscFunctionBeginUser; State state0 = StateFromY(gas, Y0); for (int i = 0; i < 5; i++) { CeedScalar dH, dH_fd; { // Calculate dH using *_fwd function CeedScalar dY[5] = {0}; dY[i] = Y0[i]; State dstate0 = StateFromY_fwd(gas, state0, dY); dH = TotalSpecificEnthalpy_fwd(gas, state0, dstate0); } { // Calculate dH_fd via finite difference approximation CeedScalar H0, H1; CeedScalar Y1[5] = {Y0[0], Y0[1], Y0[2], Y0[3], Y0[4]}; Y1[i] *= 1 + eps; State state1 = StateFromY(gas, Y1); H0 = TotalSpecificEnthalpy(gas, state0); H1 = TotalSpecificEnthalpy(gas, state1); dH_fd = (H1 - H0) / eps; } snprintf(buf, sizeof buf, "TotalSpecificEnthalpy i=%d:", i); { CeedScalar divisor_threshold = 10 * CEED_EPSILON; CeedScalar dH_err; dH_err = RelativeError(dH_fd, dH, dH_fd, divisor_threshold); if (fabs(dH_err) >= rtol) printf("%s dH error %g (expected %.10e, got %.10e)\n", buf, dH_err, dH_fd, dH); } } PetscFunctionReturn(PETSC_SUCCESS); } // @brief Verify RoeSetup_fwd function against finite-difference approximation static PetscErrorCode TestRowSetup_fwd(NewtonianIdealGasContext gas, CeedScalar rtol) { CeedScalar eps = 4e-7; // Finite difference step char buf[128]; const CeedScalar rho0[2] = {1.2, 1.4}; PetscFunctionBeginUser; for (int i = 0; i < 2; i++) { RoeWeights dR, dR_fd; { // Calculate using *_fwd function CeedScalar drho[5] = {0}; drho[i] = rho0[i]; dR = RoeSetup_fwd(rho0[0], rho0[1], drho[0], drho[1]); } { // Calculate via finite difference approximation RoeWeights R0, R1; CeedScalar rho1[5] = {rho0[0], rho0[1]}; rho1[i] *= 1 + eps; R0 = RoeSetup(rho0[0], rho0[1]); R1 = RoeSetup(rho1[0], rho1[1]); dR_fd.left = (R1.left - R0.left) / eps; dR_fd.right = (R1.right - R0.right) / eps; } snprintf(buf, sizeof buf, "RoeSetup i=%d:", i); { CeedScalar divisor_threshold = 10 * CEED_EPSILON; RoeWeights dR_err; dR_err.left = RelativeError(dR_fd.left, dR.left, dR_fd.left, divisor_threshold); dR_err.right = RelativeError(dR_fd.right, dR.right, dR_fd.right, divisor_threshold); 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); 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); } } PetscFunctionReturn(PETSC_SUCCESS); } // @brief Test Riemann solver related `*_fwd` functions via finite-difference approximation static PetscErrorCode RiemannSolverUnitTests(NewtonianIdealGasContext gas, CeedScalar rtol) { PetscFunctionBeginUser; PetscCall(TestRiemannHLL_fwd(gas, rtol, rtol, rtol)); PetscCall(TestRiemannHLLC_fwd(gas, rtol, rtol, rtol)); PetscCall(TestComputeHLLSpeeds_Roe_fwd(gas, rtol)); PetscCall(TestTotalSpecificEnthalpy_fwd(gas, rtol)); PetscCall(TestRowSetup_fwd(gas, rtol)); PetscFunctionReturn(PETSC_SUCCESS); }