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 PetscErrorCode FreestreamBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { 22 User user = *(User *)ctx; 23 MPI_Comm comm = user->comm; 24 Ceed ceed = user->ceed; 25 FreestreamContext freestream_ctx; 26 CeedQFunctionContext freestream_context; 27 RiemannFluxType riemann = RIEMANN_HLLC; 28 PetscScalar meter = user->units->meter; 29 PetscScalar second = user->units->second; 30 PetscScalar Kelvin = user->units->Kelvin; 31 PetscScalar Pascal = user->units->Pascal; 32 33 PetscFunctionBeginUser; 34 // Freestream inherits reference state. We re-dimensionalize so the defaults 35 // in -help will be visible in SI units. 36 StatePrimitive Y_inf = {.pressure = reference->pressure / Pascal, .velocity = {0}, .temperature = reference->temperature / Kelvin}; 37 for (int i = 0; i < 3; i++) Y_inf.velocity[i] = reference->velocity[i] * second / meter; 38 39 PetscOptionsBegin(comm, NULL, "Options for Freestream boundary condition", NULL); 40 PetscCall(PetscOptionsEnum("-freestream_riemann", "Riemann solver to use in freestream boundary condition", NULL, RiemannSolverTypes, 41 (PetscEnum)riemann, (PetscEnum *)&riemann, NULL)); 42 PetscCall(PetscOptionsScalar("-freestream_pressure", "Pressure at freestream condition", NULL, Y_inf.pressure, &Y_inf.pressure, NULL)); 43 PetscInt narray = 3; 44 PetscCall(PetscOptionsScalarArray("-freestream_velocity", "Velocity at freestream condition", NULL, Y_inf.velocity, &narray, NULL)); 45 PetscCall(PetscOptionsScalar("-freestream_temperature", "Temperature at freestream condition", NULL, Y_inf.temperature, &Y_inf.temperature, NULL)); 46 PetscOptionsEnd(); 47 48 switch (user->phys->state_var) { 49 case STATEVAR_CONSERVATIVE: 50 switch (riemann) { 51 case RIEMANN_HLL: 52 problem->apply_freestream.qfunction = Freestream_Conserv_HLL; 53 problem->apply_freestream.qfunction_loc = Freestream_Conserv_HLL_loc; 54 problem->apply_freestream_jacobian.qfunction = Freestream_Jacobian_Conserv_HLL; 55 problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Conserv_HLL_loc; 56 break; 57 case RIEMANN_HLLC: 58 problem->apply_freestream.qfunction = Freestream_Conserv_HLLC; 59 problem->apply_freestream.qfunction_loc = Freestream_Conserv_HLLC_loc; 60 problem->apply_freestream_jacobian.qfunction = Freestream_Jacobian_Conserv_HLLC; 61 problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Conserv_HLLC_loc; 62 break; 63 } 64 break; 65 case STATEVAR_PRIMITIVE: 66 switch (riemann) { 67 case RIEMANN_HLL: 68 problem->apply_freestream.qfunction = Freestream_Prim_HLL; 69 problem->apply_freestream.qfunction_loc = Freestream_Prim_HLL_loc; 70 problem->apply_freestream_jacobian.qfunction = Freestream_Jacobian_Prim_HLL; 71 problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Prim_HLL_loc; 72 break; 73 case RIEMANN_HLLC: 74 problem->apply_freestream.qfunction = Freestream_Prim_HLLC; 75 problem->apply_freestream.qfunction_loc = Freestream_Prim_HLLC_loc; 76 problem->apply_freestream_jacobian.qfunction = Freestream_Jacobian_Prim_HLLC; 77 problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Prim_HLLC_loc; 78 break; 79 } 80 break; 81 } 82 83 Y_inf.pressure *= Pascal; 84 for (int i = 0; i < 3; i++) Y_inf.velocity[i] *= meter / second; 85 Y_inf.temperature *= Kelvin; 86 87 State S_infty = StateFromPrimitive(newtonian_ig_ctx, Y_inf); 88 89 // -- Set freestream_ctx struct values 90 PetscCall(PetscCalloc1(1, &freestream_ctx)); 91 freestream_ctx->newtonian_ctx = *newtonian_ig_ctx; 92 freestream_ctx->S_infty = S_infty; 93 94 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &freestream_context)); 95 PetscCallCeed(ceed, CeedQFunctionContextSetData(freestream_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*freestream_ctx), freestream_ctx)); 96 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(freestream_context, CEED_MEM_HOST, FreeContextPetsc)); 97 problem->apply_freestream.qfunction_context = freestream_context; 98 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(freestream_context, &problem->apply_freestream_jacobian.qfunction_context)); 99 PetscFunctionReturn(PETSC_SUCCESS); 100 } 101 102 static const char *const OutflowTypes[] = {"RIEMANN", "PRESSURE", "OutflowType", "OUTFLOW_", NULL}; 103 typedef enum { 104 OUTFLOW_RIEMANN, 105 OUTFLOW_PRESSURE, 106 } OutflowType; 107 108 PetscErrorCode OutflowBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { 109 User user = *(User *)ctx; 110 Ceed ceed = user->ceed; 111 OutflowContext outflow_ctx; 112 OutflowType outflow_type = OUTFLOW_RIEMANN; 113 CeedQFunctionContext outflow_context; 114 const PetscScalar Kelvin = user->units->Kelvin; 115 const PetscScalar Pascal = user->units->Pascal; 116 117 PetscFunctionBeginUser; 118 CeedScalar pressure = reference->pressure / Pascal; 119 CeedScalar temperature = reference->temperature / Kelvin; 120 CeedScalar recirc = 1, softplus_velocity = 1e-2; 121 PetscOptionsBegin(user->comm, NULL, "Options for Outflow boundary condition", NULL); 122 PetscCall( 123 PetscOptionsEnum("-outflow_type", "Type of outflow condition", NULL, OutflowTypes, (PetscEnum)outflow_type, (PetscEnum *)&outflow_type, NULL)); 124 PetscCall(PetscOptionsScalar("-outflow_pressure", "Pressure at outflow condition", NULL, pressure, &pressure, NULL)); 125 if (outflow_type == OUTFLOW_RIEMANN) { 126 PetscCall(PetscOptionsScalar("-outflow_temperature", "Temperature at outflow condition", NULL, temperature, &temperature, NULL)); 127 PetscCall( 128 PetscOptionsReal("-outflow_recirc", "Fraction of recirculation to allow in exterior velocity state [0,1]", NULL, recirc, &recirc, NULL)); 129 PetscCall(PetscOptionsReal("-outflow_softplus_velocity", "Characteristic velocity of softplus regularization", NULL, softplus_velocity, 130 &softplus_velocity, NULL)); 131 } 132 PetscOptionsEnd(); 133 pressure *= Pascal; 134 temperature *= Kelvin; 135 136 switch (outflow_type) { 137 case OUTFLOW_RIEMANN: 138 switch (user->phys->state_var) { 139 case STATEVAR_CONSERVATIVE: 140 problem->apply_outflow.qfunction = RiemannOutflow_Conserv; 141 problem->apply_outflow.qfunction_loc = RiemannOutflow_Conserv_loc; 142 problem->apply_outflow_jacobian.qfunction = RiemannOutflow_Jacobian_Conserv; 143 problem->apply_outflow_jacobian.qfunction_loc = RiemannOutflow_Jacobian_Conserv_loc; 144 break; 145 case STATEVAR_PRIMITIVE: 146 problem->apply_outflow.qfunction = RiemannOutflow_Prim; 147 problem->apply_outflow.qfunction_loc = RiemannOutflow_Prim_loc; 148 problem->apply_outflow_jacobian.qfunction = RiemannOutflow_Jacobian_Prim; 149 problem->apply_outflow_jacobian.qfunction_loc = RiemannOutflow_Jacobian_Prim_loc; 150 break; 151 } 152 break; 153 case OUTFLOW_PRESSURE: 154 switch (user->phys->state_var) { 155 case STATEVAR_CONSERVATIVE: 156 problem->apply_outflow.qfunction = PressureOutflow_Conserv; 157 problem->apply_outflow.qfunction_loc = PressureOutflow_Conserv_loc; 158 problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian_Conserv; 159 problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_Conserv_loc; 160 break; 161 case STATEVAR_PRIMITIVE: 162 problem->apply_outflow.qfunction = PressureOutflow_Prim; 163 problem->apply_outflow.qfunction_loc = PressureOutflow_Prim_loc; 164 problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian_Prim; 165 problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_Prim_loc; 166 break; 167 } 168 break; 169 } 170 PetscCall(PetscCalloc1(1, &outflow_ctx)); 171 outflow_ctx->gas = *newtonian_ig_ctx; 172 outflow_ctx->recirc = recirc; 173 outflow_ctx->softplus_velocity = softplus_velocity; 174 outflow_ctx->pressure = pressure; 175 outflow_ctx->temperature = temperature; 176 177 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &outflow_context)); 178 PetscCallCeed(ceed, CeedQFunctionContextSetData(outflow_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*outflow_ctx), outflow_ctx)); 179 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(outflow_context, CEED_MEM_HOST, FreeContextPetsc)); 180 problem->apply_outflow.qfunction_context = outflow_context; 181 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(outflow_context, &problem->apply_outflow_jacobian.qfunction_context)); 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184