Skip to content

Commit

Permalink
Merge pull request #1566 from CEED/jrwrigh/blasius_state
Browse files Browse the repository at this point in the history 
fluids: Use State in blasius context
  • Loading branch information
@jrwrigh
jrwrigh committed Apr 25, 2024
2 parents 5dd8ea6 + ff9b3c0 commit 8e16760
Show file tree
Hide file tree
Showing 6 changed files with 66 additions and 54 deletions.
10 changes: 5 additions & 5 deletions examples/fluids/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -948,6 +948,11 @@ The Blasius problem has the following command-line options in addition to the Ne
- `288`
- `K`

* - `-pressure_infinity`
- Atmospheric pressure, also sets IDL reference pressure
- `1.01E5`
- `Pa`

* - `-temperature_wall`
- Wall temperature
- `288`
Expand All @@ -958,11 +963,6 @@ The Blasius problem has the following command-line options in addition to the Ne
- `4.2e-3`
- `m`

* - `-P0`
- Atmospheric pressure
- `1.01E5`
- `Pa`

* - `-platemesh_modify_mesh`
- Whether to modify the mesh using the given options below.
- `false`
Expand Down
43 changes: 26 additions & 17 deletions examples/fluids/problems/blasius.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,12 @@ PetscErrorCode CompressibleBlasiusResidual(SNES snes, Vec X, Vec R, void *ctx) {
const BlasiusContext blasius = (BlasiusContext)ctx;
const PetscScalar *Tf, *Th; // Chebyshev coefficients
PetscScalar *r, f[4], h[4];
PetscInt N = blasius->n_cheb;
PetscInt N = blasius->n_cheb;
State S_infty = blasius->S_infty;
CeedScalar U_infty = sqrt(Dot3(S_infty.Y.velocity, S_infty.Y.velocity));

PetscFunctionBeginUser;
PetscScalar Ma = Mach(&blasius->newtonian_ctx, blasius->T_inf, blasius->U_inf), Pr = Prandtl(&blasius->newtonian_ctx),
PetscScalar Ma = Mach(&blasius->newtonian_ctx, S_infty.Y.temperature, U_infty), Pr = Prandtl(&blasius->newtonian_ctx),
gamma = HeatCapacityRatio(&blasius->newtonian_ctx);

PetscCall(VecGetArrayRead(X, &Tf));
Expand Down Expand Up @@ -59,7 +61,7 @@ PetscErrorCode CompressibleBlasiusResidual(SNES snes, Vec X, Vec R, void *ctx) {

// h - left end boundary condition
ChebyshevEval(N - 1, Th, -1., blasius->eta_max, h);
r[N] = h[0] - blasius->T_wall / blasius->T_inf;
r[N] = h[0] - blasius->T_wall / S_infty.Y.temperature;

// h - right end boundary condition
ChebyshevEval(N - 1, Th, 1., blasius->eta_max, h);
Expand Down Expand Up @@ -252,22 +254,25 @@ PetscErrorCode NS_BLASIUS(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
CeedScalar T_inf = 288.; // K
CeedScalar T_wall = 288.; // K
CeedScalar delta0 = 4.2e-3; // m
CeedScalar P0 = 1.01e5; // Pa
CeedScalar P_inf = 1.01e5; // Pa
PetscInt N = 20; // Number of Chebyshev terms
PetscBool weakT = PETSC_FALSE; // weak density or temperature
PetscReal mesh_refine_height = 5.9e-4; // m
PetscReal mesh_growth = 1.08; // [-]
PetscInt mesh_Ndelta = 45; // [-]
PetscReal mesh_top_angle = 5; // degrees
char mesh_ynodes_path[PETSC_MAX_PATH_LEN] = "";
PetscBool flg;

PetscOptionsBegin(comm, NULL, "Options for BLASIUS problem", NULL);
PetscCall(PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", NULL, weakT, &weakT, NULL));
PetscCall(PetscOptionsScalar("-velocity_infinity", "Velocity at boundary layer edge", NULL, U_inf, &U_inf, NULL));
PetscCall(PetscOptionsScalar("-temperature_infinity", "Temperature at boundary layer edge", NULL, T_inf, &T_inf, NULL));
PetscCall(PetscOptionsScalar("-pressure_infinity", "Pressure at boundary layer edge", NULL, P_inf, &P_inf, &flg));
PetscCall(PetscOptionsDeprecated("-P0", "-pressure_infinity", "libCEED 0.12.0", "Use -pressure_infinity to set pressure at boundary layer edge"));
if (!flg) PetscCall(PetscOptionsScalar("-P0", "Pressure at boundary layer edge", NULL, P_inf, &P_inf, &flg));
PetscCall(PetscOptionsScalar("-temperature_wall", "Temperature at wall", NULL, T_wall, &T_wall, NULL));
PetscCall(PetscOptionsScalar("-delta0", "Boundary layer height at inflow", NULL, delta0, &delta0, NULL));
PetscCall(PetscOptionsScalar("-P0", "Pressure at outflow", NULL, P0, &P0, NULL));
PetscCall(PetscOptionsInt("-n_chebyshev", "Number of Chebyshev terms", NULL, N, &N, NULL));
PetscCheck(3 <= N && N <= BLASIUS_MAX_N_CHEBYSHEV, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n_chebyshev %" PetscInt_FMT " must be in range [3, %d]", N,
BLASIUS_MAX_N_CHEBYSHEV);
Expand All @@ -293,7 +298,7 @@ PetscErrorCode NS_BLASIUS(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {

T_inf *= Kelvin;
T_wall *= Kelvin;
P0 *= Pascal;
P_inf *= Pascal;
U_inf *= meter / second;
delta0 *= meter;

Expand All @@ -308,16 +313,20 @@ PetscErrorCode NS_BLASIUS(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
// Some properties depend on parameters from NewtonianIdealGas
PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));

blasius_ctx->weakT = weakT;
blasius_ctx->U_inf = U_inf;
blasius_ctx->T_inf = T_inf;
blasius_ctx->T_wall = T_wall;
blasius_ctx->delta0 = delta0;
blasius_ctx->P0 = P0;
blasius_ctx->n_cheb = N;
newtonian_ig_ctx->P0 = P0;
blasius_ctx->implicit = user->phys->implicit;
blasius_ctx->newtonian_ctx = *newtonian_ig_ctx;
StatePrimitive Y_inf = {
.pressure = P_inf, .velocity = {U_inf, 0, 0},
.temperature = T_inf
};
State S_infty = StateFromPrimitive(newtonian_ig_ctx, Y_inf);

blasius_ctx->weakT = weakT;
blasius_ctx->T_wall = T_wall;
blasius_ctx->delta0 = delta0;
blasius_ctx->S_infty = S_infty;
blasius_ctx->n_cheb = N;
newtonian_ig_ctx->idl_pressure = P_inf;
blasius_ctx->implicit = user->phys->implicit;
blasius_ctx->newtonian_ctx = *newtonian_ig_ctx;

{
PetscReal domain_min[3], domain_max[3];
Expand All @@ -338,7 +347,7 @@ PetscErrorCode NS_BLASIUS(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context));
problem->ics.qfunction_context = blasius_context;
if (use_stg) {
PetscCall(SetupStg(comm, dm, problem, user, weakT, T_inf, P0));
PetscCall(SetupStg(comm, dm, problem, user, weakT, S_infty.Y.temperature, S_infty.Y.pressure));
} else if (diff_filter_mms) {
PetscCall(DifferentialFilterMmsICSetup(problem));
} else {
Expand Down
3 changes: 1 addition & 2 deletions examples/fluids/problems/newtonian.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,9 +322,8 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData problem, DM dm, void *ctx, SimpleBC b
newtonian_ig_ctx->Ctau_C = Ctau_C;
newtonian_ig_ctx->Ctau_M = Ctau_M;
newtonian_ig_ctx->Ctau_E = Ctau_E;
newtonian_ig_ctx->P0 = reference.pressure;
newtonian_ig_ctx->idl_pressure = reference.pressure;
newtonian_ig_ctx->stabilization = stab;
newtonian_ig_ctx->P0 = reference.pressure;
newtonian_ig_ctx->is_implicit = implicit;
newtonian_ig_ctx->state_var = state_var;
newtonian_ig_ctx->idl_enable = idl_enable;
Expand Down
60 changes: 32 additions & 28 deletions examples/fluids/qfunctions/blasius.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,11 @@

typedef struct BlasiusContext_ *BlasiusContext;
struct BlasiusContext_ {
bool implicit; // !< Using implicit timesteping or not
bool weakT; // !< flag to set Temperature weakly at inflow
CeedScalar delta0; // !< Boundary layer height at inflow
CeedScalar U_inf; // !< Velocity at boundary layer edge
CeedScalar T_inf; // !< Temperature at boundary layer edge
bool implicit; // !< Using implicit timesteping or not
bool weakT; // !< flag to set Temperature weakly at inflow
CeedScalar delta0; // !< Boundary layer height at inflow
State S_infty;
CeedScalar T_wall; // !< Temperature at the wall
CeedScalar P0; // !< Pressure at outflow
CeedScalar x_inflow; // !< Location of inflow in x
CeedScalar n_cheb; // !< Number of Chebyshev terms
CeedScalar *X; // !< Chebyshev polynomial coordinate vector (CPU only)
Expand Down Expand Up @@ -72,25 +70,26 @@ CEED_QFUNCTION_HELPER void ChebyshevEval(int N, const double *Tf, double x, doub
// *****************************************************************************
State CEED_QFUNCTION_HELPER(BlasiusSolution)(const BlasiusContext blasius, const CeedScalar x[3], const CeedScalar x0, const CeedScalar x_inflow,
const CeedScalar rho_infty, CeedScalar *t12) {
CeedInt N = blasius->n_cheb;
CeedScalar mu = blasius->newtonian_ctx.mu;
CeedScalar nu = mu / rho_infty;
CeedScalar eta = x[1] * sqrt(blasius->U_inf / (nu * (x0 + x[0] - x_inflow)));
CeedScalar X = 2 * (eta / blasius->eta_max) - 1.;
CeedScalar U_inf = blasius->U_inf;
CeedScalar Rd = GasConstant(&blasius->newtonian_ctx);
CeedInt N = blasius->n_cheb;
CeedScalar mu = blasius->newtonian_ctx.mu;
State S_infty = blasius->S_infty;
CeedScalar nu = mu / rho_infty;
CeedScalar U_infty = sqrt(Dot3(S_infty.Y.velocity, S_infty.Y.velocity));
CeedScalar eta = x[1] * sqrt(U_infty / (nu * (x0 + x[0] - x_inflow)));
CeedScalar X = 2 * (eta / blasius->eta_max) - 1.;
CeedScalar Rd = GasConstant(&blasius->newtonian_ctx);

CeedScalar f[4], h[4];
ChebyshevEval(N, blasius->Tf_cheb, X, blasius->eta_max, f);
ChebyshevEval(N - 1, blasius->Th_cheb, X, blasius->eta_max, h);

*t12 = mu * U_inf * f[2] * sqrt(U_inf / (nu * (x0 + x[0] - x_inflow)));
*t12 = mu * U_infty * f[2] * sqrt(U_infty / (nu * (x0 + x[0] - x_inflow)));

CeedScalar Y[5];
Y[1] = U_inf * f[1];
Y[2] = 0.5 * sqrt(nu * U_inf / (x0 + x[0] - x_inflow)) * (eta * f[1] - f[0]);
Y[1] = U_infty * f[1];
Y[2] = 0.5 * sqrt(nu * U_infty / (x0 + x[0] - x_inflow)) * (eta * f[1] - f[0]);
Y[3] = 0.;
Y[4] = blasius->T_inf * h[0];
Y[4] = S_infty.Y.temperature * h[0];
Y[0] = rho_infty / h[0] * Rd * Y[4];
return StateFromY(&blasius->newtonian_ctx, Y);
}
Expand All @@ -109,14 +108,14 @@ CEED_QFUNCTION(ICsBlasius)(void *ctx, CeedInt Q, const CeedScalar *const *in, Ce
const CeedScalar x_inflow = context->x_inflow;
CeedScalar t12;

const CeedScalar Y_inf[5] = {context->P0, context->U_inf, 0, 0, context->T_inf};
const State s_inf = StateFromY(gas, Y_inf);
const State S_infty = context->S_infty;
const CeedScalar U_infty = sqrt(Dot3(S_infty.Y.velocity, S_infty.Y.velocity));

const CeedScalar x0 = context->U_inf * s_inf.U.density / (mu * 25 / Square(delta0));
const CeedScalar x0 = U_infty * S_infty.U.density / (mu * 25 / Square(delta0));

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]};
State s = BlasiusSolution(context, x, x0, x_inflow, s_inf.U.density, &t12);
State s = BlasiusSolution(context, x, x0, x_inflow, S_infty.U.density, &t12);
CeedScalar q[5] = {0};

switch (gas->state_var) {
Expand All @@ -143,8 +142,10 @@ CEED_QFUNCTION(Blasius_Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in

const bool is_implicit = context->implicit;
const NewtonianIdealGasContext gas = &context->newtonian_ctx;
const CeedScalar rho_0 = context->P0 / (GasConstant(gas) * context->T_inf);
const CeedScalar x0 = context->U_inf * rho_0 / (gas->mu * 25 / Square(context->delta0));
State S_infty = context->S_infty;
const CeedScalar rho_0 = S_infty.U.density;
const CeedScalar U_infty = sqrt(Dot3(S_infty.Y.velocity, S_infty.Y.velocity));
const CeedScalar x0 = U_infty * rho_0 / (gas->mu * 25 / Square(context->delta0));
const CeedScalar zeros[11] = {0.};

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
Expand Down Expand Up @@ -198,8 +199,10 @@ CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, const CeedScalar *
const bool is_implicit = context->implicit;
const CeedScalar Rd = GasConstant(gas);
const CeedScalar gamma = HeatCapacityRatio(gas);
const CeedScalar rho_0 = context->P0 / (Rd * context->T_inf);
const CeedScalar x0 = context->U_inf * rho_0 / (gas->mu * 25 / (Square(context->delta0)));
const State S_infty = context->S_infty;
const CeedScalar rho_0 = S_infty.U.density;
const CeedScalar U_infty = sqrt(Dot3(S_infty.Y.velocity, S_infty.Y.velocity));
const CeedScalar x0 = U_infty * rho_0 / (gas->mu * 25 / Square(context->delta0));

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
CeedScalar wdetJb, norm[3];
Expand All @@ -216,11 +219,12 @@ CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, const CeedScalar *
if (context->weakT) {
// rho should be from the current solution
drho = dq[0][i];
CeedScalar dE_internal = drho * gas->cv * context->T_inf;
CeedScalar dE_internal = drho * gas->cv * S_infty.Y.temperature;
CeedScalar dE_kinetic = .5 * drho * Dot3(s.Y.velocity, s.Y.velocity);
dE = dE_internal + dE_kinetic;
dP = drho * Rd * context->T_inf; // interior rho with exterior T
} else { // rho specified, E_internal from solution
dP = drho * Rd * S_infty.Y.temperature; // interior rho with exterior T
} else {
// rho specified, E_internal from solution
drho = 0;
dE = dq[4][i];
dP = dE * (gamma - 1.);
Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/qfunctions/newtonian.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,7 @@ CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedSc
NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
const CeedScalar *g = context->g;
const CeedScalar dt = context->dt;
const CeedScalar P0 = context->P0;
const CeedScalar P0 = context->idl_pressure;

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/qfunctions/newtonian_types.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,11 @@ struct NewtonianIdealGasContext_ {
CeedScalar dt;
CeedScalar time;
CeedScalar ijacobian_time_shift;
CeedScalar P0;
bool is_implicit;
StateVariable state_var;
StabilizationType stabilization;
bool idl_enable;
CeedScalar idl_pressure;
CeedScalar idl_amplitude;
CeedScalar idl_start;
CeedScalar idl_length;
Expand Down

0 comments on commit 8e16760

Please sign in to comment.