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 case STATEVAR_ENTROPY: 82 switch (riemann) { 83 case RIEMANN_HLL: 84 problem->apply_freestream.qfunction = Freestream_Entropy_HLL; 85 problem->apply_freestream.qfunction_loc = Freestream_Entropy_HLL_loc; 86 problem->apply_freestream_jacobian.qfunction = Freestream_Jacobian_Entropy_HLL; 87 problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Entropy_HLL_loc; 88 break; 89 case RIEMANN_HLLC: 90 problem->apply_freestream.qfunction = Freestream_Entropy_HLLC; 91 problem->apply_freestream.qfunction_loc = Freestream_Entropy_HLLC_loc; 92 problem->apply_freestream_jacobian.qfunction = Freestream_Jacobian_Entropy_HLLC; 93 problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Entropy_HLLC_loc; 94 break; 95 } 96 break; 97 } 98 99 Y_inf.pressure *= Pascal; 100 for (int i = 0; i < 3; i++) Y_inf.velocity[i] *= meter / second; 101 Y_inf.temperature *= Kelvin; 102 103 State S_infty = StateFromPrimitive(newtonian_ig_ctx, Y_inf); 104 105 // -- Set freestream_ctx struct values 106 PetscCall(PetscCalloc1(1, &freestream_ctx)); 107 freestream_ctx->newtonian_ctx = *newtonian_ig_ctx; 108 freestream_ctx->S_infty = S_infty; 109 110 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &freestream_context)); 111 PetscCallCeed(ceed, CeedQFunctionContextSetData(freestream_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*freestream_ctx), freestream_ctx)); 112 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(freestream_context, CEED_MEM_HOST, FreeContextPetsc)); 113 problem->apply_freestream.qfunction_context = freestream_context; 114 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(freestream_context, &problem->apply_freestream_jacobian.qfunction_context)); 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 static const char *const OutflowTypes[] = {"RIEMANN", "PRESSURE", "OutflowType", "OUTFLOW_", NULL}; 119 typedef enum { 120 OUTFLOW_RIEMANN, 121 OUTFLOW_PRESSURE, 122 } OutflowType; 123 124 PetscErrorCode OutflowBCSetup(ProblemData problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference) { 125 User user = *(User *)ctx; 126 Ceed ceed = user->ceed; 127 OutflowContext outflow_ctx; 128 OutflowType outflow_type = OUTFLOW_RIEMANN; 129 CeedQFunctionContext outflow_context; 130 const PetscScalar Kelvin = user->units->Kelvin; 131 const PetscScalar Pascal = user->units->Pascal; 132 133 PetscFunctionBeginUser; 134 CeedScalar pressure = reference->pressure / Pascal; 135 CeedScalar temperature = reference->temperature / Kelvin; 136 CeedScalar recirc = 1, softplus_velocity = 1e-2; 137 PetscOptionsBegin(user->comm, NULL, "Options for Outflow boundary condition", NULL); 138 PetscCall( 139 PetscOptionsEnum("-outflow_type", "Type of outflow condition", NULL, OutflowTypes, (PetscEnum)outflow_type, (PetscEnum *)&outflow_type, NULL)); 140 PetscCall(PetscOptionsScalar("-outflow_pressure", "Pressure at outflow condition", NULL, pressure, &pressure, NULL)); 141 if (outflow_type == OUTFLOW_RIEMANN) { 142 PetscCall(PetscOptionsScalar("-outflow_temperature", "Temperature at outflow condition", NULL, temperature, &temperature, NULL)); 143 PetscCall( 144 PetscOptionsReal("-outflow_recirc", "Fraction of recirculation to allow in exterior velocity state [0,1]", NULL, recirc, &recirc, NULL)); 145 PetscCall(PetscOptionsReal("-outflow_softplus_velocity", "Characteristic velocity of softplus regularization", NULL, softplus_velocity, 146 &softplus_velocity, NULL)); 147 } 148 PetscOptionsEnd(); 149 pressure *= Pascal; 150 temperature *= Kelvin; 151 152 switch (outflow_type) { 153 case OUTFLOW_RIEMANN: 154 switch (user->phys->state_var) { 155 case STATEVAR_CONSERVATIVE: 156 problem->apply_outflow.qfunction = RiemannOutflow_Conserv; 157 problem->apply_outflow.qfunction_loc = RiemannOutflow_Conserv_loc; 158 problem->apply_outflow_jacobian.qfunction = RiemannOutflow_Jacobian_Conserv; 159 problem->apply_outflow_jacobian.qfunction_loc = RiemannOutflow_Jacobian_Conserv_loc; 160 break; 161 case STATEVAR_PRIMITIVE: 162 problem->apply_outflow.qfunction = RiemannOutflow_Prim; 163 problem->apply_outflow.qfunction_loc = RiemannOutflow_Prim_loc; 164 problem->apply_outflow_jacobian.qfunction = RiemannOutflow_Jacobian_Prim; 165 problem->apply_outflow_jacobian.qfunction_loc = RiemannOutflow_Jacobian_Prim_loc; 166 break; 167 case STATEVAR_ENTROPY: 168 problem->apply_outflow.qfunction = RiemannOutflow_Entropy; 169 problem->apply_outflow.qfunction_loc = RiemannOutflow_Entropy_loc; 170 problem->apply_outflow_jacobian.qfunction = RiemannOutflow_Jacobian_Entropy; 171 problem->apply_outflow_jacobian.qfunction_loc = RiemannOutflow_Jacobian_Entropy_loc; 172 break; 173 } 174 break; 175 case OUTFLOW_PRESSURE: 176 switch (user->phys->state_var) { 177 case STATEVAR_CONSERVATIVE: 178 problem->apply_outflow.qfunction = PressureOutflow_Conserv; 179 problem->apply_outflow.qfunction_loc = PressureOutflow_Conserv_loc; 180 problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian_Conserv; 181 problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_Conserv_loc; 182 break; 183 case STATEVAR_PRIMITIVE: 184 problem->apply_outflow.qfunction = PressureOutflow_Prim; 185 problem->apply_outflow.qfunction_loc = PressureOutflow_Prim_loc; 186 problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian_Prim; 187 problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_Prim_loc; 188 break; 189 case STATEVAR_ENTROPY: 190 problem->apply_outflow.qfunction = PressureOutflow_Entropy; 191 problem->apply_outflow.qfunction_loc = PressureOutflow_Entropy_loc; 192 problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian_Entropy; 193 problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_Entropy_loc; 194 break; 195 } 196 break; 197 } 198 PetscCall(PetscCalloc1(1, &outflow_ctx)); 199 outflow_ctx->gas = *newtonian_ig_ctx; 200 outflow_ctx->recirc = recirc; 201 outflow_ctx->softplus_velocity = softplus_velocity; 202 outflow_ctx->pressure = pressure; 203 outflow_ctx->temperature = temperature; 204 205 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &outflow_context)); 206 PetscCallCeed(ceed, CeedQFunctionContextSetData(outflow_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*outflow_ctx), outflow_ctx)); 207 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(outflow_context, CEED_MEM_HOST, FreeContextPetsc)); 208 problem->apply_outflow.qfunction_context = outflow_context; 209 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(outflow_context, &problem->apply_outflow_jacobian.qfunction_context)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212