1 // Copyright (c) 2017-2022, 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 /// QFunctions for the `bc_freestream` and `bc_outflow` boundary conditions 10 11 #include "bc_freestream_type.h" 12 #include "newtonian_state.h" 13 #include "newtonian_types.h" 14 #include "riemann_solver.h" 15 16 // ***************************************************************************** 17 // Freestream Boundary Condition 18 // ***************************************************************************** 19 CEED_QFUNCTION_HELPER int Freestream(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var, 20 RiemannFluxType flux_type) { 21 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 22 const CeedScalar(*q_data_sur) = in[2]; 23 24 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 25 CeedScalar(*jac_data_sur) = out[1]; 26 27 const FreestreamContext context = (FreestreamContext)ctx; 28 const NewtonianIdealGasContext newt_ctx = &context->newtonian_ctx; 29 const bool is_implicit = newt_ctx->is_implicit; 30 const CeedScalar zeros[6] = {0.}; 31 32 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 33 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 34 State s = StateFromQ(newt_ctx, qi, state_var); 35 36 CeedScalar wdetJb, norm[3]; 37 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 38 wdetJb *= is_implicit ? -1. : 1.; 39 40 StateConservative flux; 41 switch (flux_type) { 42 case RIEMANN_HLL: 43 flux = RiemannFlux_HLL(newt_ctx, s, context->S_infty, norm); 44 break; 45 case RIEMANN_HLLC: 46 flux = RiemannFlux_HLLC(newt_ctx, s, context->S_infty, norm); 47 break; 48 } 49 CeedScalar Flux[5]; 50 UnpackState_U(flux, Flux); 51 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 52 53 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 54 StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set 55 } 56 return 0; 57 } 58 59 CEED_QFUNCTION(Freestream_Conserv_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 60 return Freestream(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLL); 61 } 62 63 CEED_QFUNCTION(Freestream_Prim_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 64 return Freestream(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLL); 65 } 66 67 CEED_QFUNCTION(Freestream_Conserv_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 68 return Freestream(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLLC); 69 } 70 71 CEED_QFUNCTION(Freestream_Prim_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 72 return Freestream(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLLC); 73 } 74 75 CEED_QFUNCTION_HELPER int Freestream_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var, 76 RiemannFluxType flux_type) { 77 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 78 const CeedScalar(*q_data_sur) = in[2]; 79 const CeedScalar(*jac_data_sur) = in[4]; 80 81 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 82 83 const FreestreamContext context = (FreestreamContext)ctx; 84 const NewtonianIdealGasContext newt_ctx = &context->newtonian_ctx; 85 const bool is_implicit = newt_ctx->is_implicit; 86 const State dS_infty = {0}; 87 88 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 89 CeedScalar wdetJb, norm[3]; 90 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 91 wdetJb *= is_implicit ? -1. : 1.; 92 93 CeedScalar qi[5], dqi[5]; 94 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 95 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 96 State s = StateFromQ(newt_ctx, qi, state_var); 97 State ds = StateFromQ_fwd(newt_ctx, s, dqi, state_var); 98 99 StateConservative dflux; 100 switch (flux_type) { 101 case RIEMANN_HLL: 102 dflux = RiemannFlux_HLL_fwd(newt_ctx, s, ds, context->S_infty, dS_infty, norm); 103 break; 104 case RIEMANN_HLLC: 105 dflux = RiemannFlux_HLLC_fwd(newt_ctx, s, ds, context->S_infty, dS_infty, norm); 106 break; 107 } 108 CeedScalar dFlux[5]; 109 UnpackState_U(dflux, dFlux); 110 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 111 } 112 return 0; 113 } 114 115 CEED_QFUNCTION(Freestream_Jacobian_Conserv_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 116 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLL); 117 } 118 119 CEED_QFUNCTION(Freestream_Jacobian_Prim_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 120 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLL); 121 } 122 123 CEED_QFUNCTION(Freestream_Jacobian_Conserv_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 124 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLLC); 125 } 126 127 CEED_QFUNCTION(Freestream_Jacobian_Prim_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 128 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLLC); 129 } 130 131 // Note the identity 132 // 133 // softplus(x) - x = log(1 + exp(x)) - x 134 // = log(1 + exp(x)) + log(exp(-x)) 135 // = log((1 + exp(x)) * exp(-x)) 136 // = log(exp(-x) + 1) 137 // = softplus(-x) 138 CEED_QFUNCTION_HELPER CeedScalar Softplus(CeedScalar x, CeedScalar width) { 139 if (x > 40 * width) return x; 140 return width * log1p(exp(x / width)); 141 } 142 143 CEED_QFUNCTION_HELPER CeedScalar Softplus_fwd(CeedScalar x, CeedScalar dx, CeedScalar width) { 144 if (x > 40 * width) return 1; 145 const CeedScalar t = exp(x / width); 146 return t / (1 + t); 147 } 148 149 // Viscous Outflow boundary condition, setting a constant exterior pressure and 150 // temperature as input for a Riemann solve. This condition is stable even in 151 // recirculating flow so long as the exterior temperature is sensible. 152 // 153 // The velocity in the exterior state has optional softplus regularization to 154 // keep it outflow. These parameters have been finnicky in practice and provide 155 // little or no benefit in the tests we've run thus far, thus we recommend 156 // skipping this feature and just allowing recirculation. 157 CEED_QFUNCTION_HELPER int RiemannOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 158 // Inputs 159 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 160 const CeedScalar(*Grad_q) = in[1]; 161 const CeedScalar(*q_data_sur) = in[2]; 162 163 // Outputs 164 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 165 CeedScalar(*jac_data_sur) = out[1]; 166 167 const OutflowContext outflow = (OutflowContext)ctx; 168 const NewtonianIdealGasContext gas = &outflow->gas; 169 const bool is_implicit = gas->is_implicit; 170 171 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 172 CeedScalar wdetJb, dXdx[2][3], norm[3]; 173 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 174 wdetJb *= is_implicit ? -1. : 1.; 175 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 176 State s_int = StateFromQ(gas, qi, state_var); 177 178 StatePrimitive y_ext = s_int.Y; 179 y_ext.pressure = outflow->pressure; 180 y_ext.temperature = outflow->temperature; 181 const CeedScalar u_normal = Dot3(y_ext.velocity, norm); 182 const CeedScalar proj = (1 - outflow->recirc) * Softplus(-u_normal, outflow->softplus_velocity); 183 for (CeedInt j = 0; j < 3; j++) { 184 y_ext.velocity[j] += norm[j] * proj; // (I - n n^T) projects into the plane tangent to the normal 185 } 186 State s_ext = StateFromPrimitive(gas, y_ext); 187 188 State grad_s[3]; 189 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s_int, state_var, Grad_q, dXdx, grad_s); 190 191 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 192 KMStrainRate_State(grad_s, strain_rate); 193 NewtonianStress(gas, strain_rate, kmstress); 194 KMUnpack(kmstress, stress); 195 ViscousEnergyFlux(gas, s_int.Y, grad_s, stress, Fe); 196 197 StateConservative F_inviscid_normal = RiemannFlux_HLLC(gas, s_int, s_ext, norm); 198 199 CeedScalar Flux[5]; 200 FluxTotal_RiemannBoundary(F_inviscid_normal, stress, Fe, norm, Flux); 201 202 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 203 204 // Save values for Jacobian 205 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 206 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 207 } 208 return 0; 209 } 210 211 CEED_QFUNCTION(RiemannOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 212 return RiemannOutflow(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 213 } 214 215 CEED_QFUNCTION(RiemannOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 216 return RiemannOutflow(ctx, Q, in, out, STATEVAR_PRIMITIVE); 217 } 218 219 // ***************************************************************************** 220 // Jacobian for Riemann pressure/temperature outflow boundary condition 221 // ***************************************************************************** 222 CEED_QFUNCTION_HELPER int RiemannOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 223 StateVariable state_var) { 224 // Inputs 225 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 226 const CeedScalar(*Grad_dq) = in[1]; 227 const CeedScalar(*q_data_sur) = in[2]; 228 const CeedScalar(*jac_data_sur) = in[4]; 229 230 // Outputs 231 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 232 233 const OutflowContext outflow = (OutflowContext)ctx; 234 const NewtonianIdealGasContext gas = &outflow->gas; 235 const bool is_implicit = gas->is_implicit; 236 237 // Quadrature Point Loop 238 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 239 CeedScalar wdetJb, dXdx[2][3], norm[3]; 240 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 241 wdetJb *= is_implicit ? -1. : 1.; 242 243 CeedScalar qi[5], kmstress[6], dqi[5]; 244 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 245 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 246 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 247 248 State s_int = StateFromQ(gas, qi, state_var); 249 const State ds_int = StateFromQ_fwd(gas, s_int, dqi, state_var); 250 StatePrimitive y_ext = s_int.Y, dy_ext = ds_int.Y; 251 y_ext.pressure = outflow->pressure; 252 y_ext.temperature = outflow->temperature; 253 dy_ext.pressure = 0; 254 dy_ext.temperature = 0; 255 const CeedScalar u_normal = Dot3(s_int.Y.velocity, norm); 256 const CeedScalar du_normal = Dot3(ds_int.Y.velocity, norm); 257 const CeedScalar proj = (1 - outflow->recirc) * Softplus(-u_normal, outflow->softplus_velocity); 258 const CeedScalar dproj = (1 - outflow->recirc) * Softplus_fwd(-u_normal, -du_normal, outflow->softplus_velocity); 259 for (CeedInt j = 0; j < 3; j++) { 260 y_ext.velocity[j] += norm[j] * proj; 261 dy_ext.velocity[j] += norm[j] * dproj; 262 } 263 264 State s_ext = StateFromPrimitive(gas, y_ext); 265 State ds_ext = StateFromPrimitive_fwd(gas, s_ext, dy_ext); 266 267 State grad_ds[3]; 268 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s_int, state_var, Grad_dq, dXdx, grad_ds); 269 270 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 271 KMStrainRate_State(grad_ds, dstrain_rate); 272 NewtonianStress(gas, dstrain_rate, dkmstress); 273 KMUnpack(dkmstress, dstress); 274 KMUnpack(kmstress, stress); 275 ViscousEnergyFlux_fwd(gas, s_int.Y, ds_int.Y, grad_ds, stress, dstress, dFe); 276 277 StateConservative dF_inviscid_normal = RiemannFlux_HLLC_fwd(gas, s_int, ds_int, s_ext, ds_ext, norm); 278 279 CeedScalar dFlux[5]; 280 FluxTotal_RiemannBoundary(dF_inviscid_normal, dstress, dFe, norm, dFlux); 281 282 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 283 } // End Quadrature Point Loop 284 return 0; 285 } 286 287 CEED_QFUNCTION(RiemannOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 288 return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 289 } 290 291 CEED_QFUNCTION(RiemannOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 292 return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 293 } 294 295 // ***************************************************************************** 296 // Outflow boundary condition, weakly setting a constant pressure. This is the 297 // classic outflow condition used by PHASTA-C and retained largely for 298 // comparison. In our experiments, it is never better than RiemannOutflow, and 299 // will crash if outflow ever becomes an inflow, as occurs with strong 300 // acoustics, vortices, etc. 301 // ***************************************************************************** 302 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 303 // Inputs 304 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 305 const CeedScalar(*Grad_q) = in[1]; 306 const CeedScalar(*q_data_sur) = in[2]; 307 308 // Outputs 309 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 310 CeedScalar(*jac_data_sur) = out[1]; 311 312 const OutflowContext outflow = (OutflowContext)ctx; 313 const NewtonianIdealGasContext gas = &outflow->gas; 314 const bool is_implicit = gas->is_implicit; 315 316 // Quadrature Point Loop 317 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 318 // Setup 319 // -- Interp in 320 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 321 State s = StateFromQ(gas, qi, state_var); 322 s.Y.pressure = outflow->pressure; 323 324 CeedScalar wdetJb, dXdx[2][3], norm[3]; 325 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 326 wdetJb *= is_implicit ? -1. : 1.; 327 328 State grad_s[3]; 329 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s); 330 331 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 332 KMStrainRate_State(grad_s, strain_rate); 333 NewtonianStress(gas, strain_rate, kmstress); 334 KMUnpack(kmstress, stress); 335 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe); 336 337 StateConservative F_inviscid[3]; 338 FluxInviscid(gas, s, F_inviscid); 339 340 CeedScalar Flux[5]; 341 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 342 343 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 344 345 // Save values for Jacobian 346 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 347 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 348 } // End Quadrature Point Loop 349 return 0; 350 } 351 352 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 353 return PressureOutflow(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 354 } 355 356 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 357 return PressureOutflow(ctx, Q, in, out, STATEVAR_PRIMITIVE); 358 } 359 360 // ***************************************************************************** 361 // Jacobian for weak-pressure outflow boundary condition 362 // ***************************************************************************** 363 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 364 StateVariable state_var) { 365 // Inputs 366 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 367 const CeedScalar(*Grad_dq) = in[1]; 368 const CeedScalar(*q_data_sur) = in[2]; 369 const CeedScalar(*jac_data_sur) = in[4]; 370 371 // Outputs 372 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 373 374 const OutflowContext outflow = (OutflowContext)ctx; 375 const NewtonianIdealGasContext gas = &outflow->gas; 376 const bool is_implicit = gas->is_implicit; 377 378 // Quadrature Point Loop 379 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 380 CeedScalar wdetJb, dXdx[2][3], norm[3]; 381 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 382 wdetJb *= is_implicit ? -1. : 1.; 383 384 CeedScalar qi[5], kmstress[6], dqi[5]; 385 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 386 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 387 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 388 389 State s = StateFromQ(gas, qi, state_var); 390 State ds = StateFromQ_fwd(gas, s, dqi, state_var); 391 s.Y.pressure = outflow->pressure; 392 ds.Y.pressure = 0.; 393 394 State grad_ds[3]; 395 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_dq, dXdx, grad_ds); 396 397 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 398 KMStrainRate_State(grad_ds, dstrain_rate); 399 NewtonianStress(gas, dstrain_rate, dkmstress); 400 KMUnpack(dkmstress, dstress); 401 KMUnpack(kmstress, stress); 402 ViscousEnergyFlux_fwd(gas, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 403 404 StateConservative dF_inviscid[3]; 405 FluxInviscid_fwd(gas, s, ds, dF_inviscid); 406 407 CeedScalar dFlux[5]; 408 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 409 410 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 411 } // End Quadrature Point Loop 412 return 0; 413 } 414 415 CEED_QFUNCTION(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 416 return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 417 } 418 419 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 420 return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 421 } 422