blasius.c (cce3ee4efcaa7914b8a0b517c757fe3477daefe6) blasius.c (ff9b3c0e2ebd8a09dedd1c00be3d2c5b29de65cc)
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

--- 7 unchanged lines hidden (view full) ---

16
17#include "../navierstokes.h"
18#include "stg_shur14.h"
19
20PetscErrorCode CompressibleBlasiusResidual(SNES snes, Vec X, Vec R, void *ctx) {
21 const BlasiusContext blasius = (BlasiusContext)ctx;
22 const PetscScalar *Tf, *Th; // Chebyshev coefficients
23 PetscScalar *r, f[4], h[4];
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

--- 7 unchanged lines hidden (view full) ---

16
17#include "../navierstokes.h"
18#include "stg_shur14.h"
19
20PetscErrorCode CompressibleBlasiusResidual(SNES snes, Vec X, Vec R, void *ctx) {
21 const BlasiusContext blasius = (BlasiusContext)ctx;
22 const PetscScalar *Tf, *Th; // Chebyshev coefficients
23 PetscScalar *r, f[4], h[4];
24 PetscInt N = blasius->n_cheb;
24 PetscInt N = blasius->n_cheb;
25 State S_infty = blasius->S_infty;
26 CeedScalar U_infty = sqrt(Dot3(S_infty.Y.velocity, S_infty.Y.velocity));
25
26 PetscFunctionBeginUser;
27
28 PetscFunctionBeginUser;
27 PetscScalar Ma = Mach(&blasius->newtonian_ctx, blasius->T_inf, blasius->U_inf), Pr = Prandtl(&blasius->newtonian_ctx),
29 PetscScalar Ma = Mach(&blasius->newtonian_ctx, S_infty.Y.temperature, U_infty), Pr = Prandtl(&blasius->newtonian_ctx),
28 gamma = HeatCapacityRatio(&blasius->newtonian_ctx);
29
30 PetscCall(VecGetArrayRead(X, &Tf));
31 Th = Tf + N;
32 PetscCall(VecGetArray(R, &r));
33
34 // Left boundary conditions f = f' = 0
35 ChebyshevEval(N, Tf, -1., blasius->eta_max, f);

--- 18 unchanged lines hidden (view full) ---

54 mu_tilde[1] * rho_tilde[0] + mu_tilde[0] * rho_tilde[1],
55 };
56 r[3 + i] = 2 * (mu_rho_tilde[0] * f[3] + mu_rho_tilde[1] * f[2]) + f[2] * f[0];
57 r[N + 2 + i] = (mu_rho_tilde[0] * h[2] + mu_rho_tilde[1] * h[1]) + Pr * f[0] * h[1] + Pr * (gamma - 1) * mu_rho_tilde[0] * PetscSqr(Ma * f[2]);
58 }
59
60 // h - left end boundary condition
61 ChebyshevEval(N - 1, Th, -1., blasius->eta_max, h);
30 gamma = HeatCapacityRatio(&blasius->newtonian_ctx);
31
32 PetscCall(VecGetArrayRead(X, &Tf));
33 Th = Tf + N;
34 PetscCall(VecGetArray(R, &r));
35
36 // Left boundary conditions f = f' = 0
37 ChebyshevEval(N, Tf, -1., blasius->eta_max, f);

--- 18 unchanged lines hidden (view full) ---

56 mu_tilde[1] * rho_tilde[0] + mu_tilde[0] * rho_tilde[1],
57 };
58 r[3 + i] = 2 * (mu_rho_tilde[0] * f[3] + mu_rho_tilde[1] * f[2]) + f[2] * f[0];
59 r[N + 2 + i] = (mu_rho_tilde[0] * h[2] + mu_rho_tilde[1] * h[1]) + Pr * f[0] * h[1] + Pr * (gamma - 1) * mu_rho_tilde[0] * PetscSqr(Ma * f[2]);
60 }
61
62 // h - left end boundary condition
63 ChebyshevEval(N - 1, Th, -1., blasius->eta_max, h);
62 r[N] = h[0] - blasius->T_wall / blasius->T_inf;
64 r[N] = h[0] - blasius->T_wall / S_infty.Y.temperature;
63
64 // h - right end boundary condition
65 ChebyshevEval(N - 1, Th, 1., blasius->eta_max, h);
66 r[N + 1] = h[0] - 1.;
67
68 // Restore vectors
69 PetscCall(VecRestoreArrayRead(X, &Tf));
70 PetscCall(VecRestoreArray(R, &r));

--- 176 unchanged lines hidden (view full) ---

247 // ------------------------------------------------------
248 problem->ics.qfunction = ICsBlasius;
249 problem->ics.qfunction_loc = ICsBlasius_loc;
250
251 CeedScalar U_inf = 40; // m/s
252 CeedScalar T_inf = 288.; // K
253 CeedScalar T_wall = 288.; // K
254 CeedScalar delta0 = 4.2e-3; // m
65
66 // h - right end boundary condition
67 ChebyshevEval(N - 1, Th, 1., blasius->eta_max, h);
68 r[N + 1] = h[0] - 1.;
69
70 // Restore vectors
71 PetscCall(VecRestoreArrayRead(X, &Tf));
72 PetscCall(VecRestoreArray(R, &r));

--- 176 unchanged lines hidden (view full) ---

249 // ------------------------------------------------------
250 problem->ics.qfunction = ICsBlasius;
251 problem->ics.qfunction_loc = ICsBlasius_loc;
252
253 CeedScalar U_inf = 40; // m/s
254 CeedScalar T_inf = 288.; // K
255 CeedScalar T_wall = 288.; // K
256 CeedScalar delta0 = 4.2e-3; // m
255 CeedScalar P0 = 1.01e5; // Pa
257 CeedScalar P_inf = 1.01e5; // Pa
256 PetscInt N = 20; // Number of Chebyshev terms
257 PetscBool weakT = PETSC_FALSE; // weak density or temperature
258 PetscReal mesh_refine_height = 5.9e-4; // m
259 PetscReal mesh_growth = 1.08; // [-]
260 PetscInt mesh_Ndelta = 45; // [-]
261 PetscReal mesh_top_angle = 5; // degrees
262 char mesh_ynodes_path[PETSC_MAX_PATH_LEN] = "";
258 PetscInt N = 20; // Number of Chebyshev terms
259 PetscBool weakT = PETSC_FALSE; // weak density or temperature
260 PetscReal mesh_refine_height = 5.9e-4; // m
261 PetscReal mesh_growth = 1.08; // [-]
262 PetscInt mesh_Ndelta = 45; // [-]
263 PetscReal mesh_top_angle = 5; // degrees
264 char mesh_ynodes_path[PETSC_MAX_PATH_LEN] = "";
265 PetscBool flg;
263
264 PetscOptionsBegin(comm, NULL, "Options for BLASIUS problem", NULL);
265 PetscCall(PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", NULL, weakT, &weakT, NULL));
266 PetscCall(PetscOptionsScalar("-velocity_infinity", "Velocity at boundary layer edge", NULL, U_inf, &U_inf, NULL));
267 PetscCall(PetscOptionsScalar("-temperature_infinity", "Temperature at boundary layer edge", NULL, T_inf, &T_inf, NULL));
266
267 PetscOptionsBegin(comm, NULL, "Options for BLASIUS problem", NULL);
268 PetscCall(PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", NULL, weakT, &weakT, NULL));
269 PetscCall(PetscOptionsScalar("-velocity_infinity", "Velocity at boundary layer edge", NULL, U_inf, &U_inf, NULL));
270 PetscCall(PetscOptionsScalar("-temperature_infinity", "Temperature at boundary layer edge", NULL, T_inf, &T_inf, NULL));
271 PetscCall(PetscOptionsScalar("-pressure_infinity", "Pressure at boundary layer edge", NULL, P_inf, &P_inf, &flg));
272 PetscCall(PetscOptionsDeprecated("-P0", "-pressure_infinity", "libCEED 0.12.0", "Use -pressure_infinity to set pressure at boundary layer edge"));
273 if (!flg) PetscCall(PetscOptionsScalar("-P0", "Pressure at boundary layer edge", NULL, P_inf, &P_inf, &flg));
268 PetscCall(PetscOptionsScalar("-temperature_wall", "Temperature at wall", NULL, T_wall, &T_wall, NULL));
269 PetscCall(PetscOptionsScalar("-delta0", "Boundary layer height at inflow", NULL, delta0, &delta0, NULL));
274 PetscCall(PetscOptionsScalar("-temperature_wall", "Temperature at wall", NULL, T_wall, &T_wall, NULL));
275 PetscCall(PetscOptionsScalar("-delta0", "Boundary layer height at inflow", NULL, delta0, &delta0, NULL));
270 PetscCall(PetscOptionsScalar("-P0", "Pressure at outflow", NULL, P0, &P0, NULL));
271 PetscCall(PetscOptionsInt("-n_chebyshev", "Number of Chebyshev terms", NULL, N, &N, NULL));
272 PetscCheck(3 <= N && N <= BLASIUS_MAX_N_CHEBYSHEV, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n_chebyshev %" PetscInt_FMT " must be in range [3, %d]", N,
273 BLASIUS_MAX_N_CHEBYSHEV);
274 if (user->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) {
275 PetscCall(PetscOptionsBoundedInt("-platemesh_Ndelta", "Velocity at boundary layer edge", NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1));
276 PetscCall(PetscOptionsScalar("-platemesh_refine_height", "Height of boundary layer mesh refinement", NULL, mesh_refine_height,
277 &mesh_refine_height, NULL));
278 PetscCall(PetscOptionsScalar("-platemesh_growth", "Geometric growth rate of boundary layer mesh", NULL, mesh_growth, &mesh_growth, NULL));

--- 9 unchanged lines hidden (view full) ---

288
289 PetscScalar meter = user->units->meter;
290 PetscScalar second = user->units->second;
291 PetscScalar Kelvin = user->units->Kelvin;
292 PetscScalar Pascal = user->units->Pascal;
293
294 T_inf *= Kelvin;
295 T_wall *= Kelvin;
276 PetscCall(PetscOptionsInt("-n_chebyshev", "Number of Chebyshev terms", NULL, N, &N, NULL));
277 PetscCheck(3 <= N && N <= BLASIUS_MAX_N_CHEBYSHEV, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n_chebyshev %" PetscInt_FMT " must be in range [3, %d]", N,
278 BLASIUS_MAX_N_CHEBYSHEV);
279 if (user->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) {
280 PetscCall(PetscOptionsBoundedInt("-platemesh_Ndelta", "Velocity at boundary layer edge", NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1));
281 PetscCall(PetscOptionsScalar("-platemesh_refine_height", "Height of boundary layer mesh refinement", NULL, mesh_refine_height,
282 &mesh_refine_height, NULL));
283 PetscCall(PetscOptionsScalar("-platemesh_growth", "Geometric growth rate of boundary layer mesh", NULL, mesh_growth, &mesh_growth, NULL));

--- 9 unchanged lines hidden (view full) ---

293
294 PetscScalar meter = user->units->meter;
295 PetscScalar second = user->units->second;
296 PetscScalar Kelvin = user->units->Kelvin;
297 PetscScalar Pascal = user->units->Pascal;
298
299 T_inf *= Kelvin;
300 T_wall *= Kelvin;
296 P0 *= Pascal;
301 P_inf *= Pascal;
297 U_inf *= meter / second;
298 delta0 *= meter;
299
300 if (user->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) {
301 PetscReal *mesh_ynodes = NULL;
302 PetscInt mesh_nynodes = 0;
303 if (strcmp(mesh_ynodes_path, "")) PetscCall(GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes));
304 PetscCall(ModifyMesh(comm, dm, problem->dim, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle, &mesh_ynodes, &mesh_nynodes));
305 PetscCall(PetscFree(mesh_ynodes));
306 }
307
308 // Some properties depend on parameters from NewtonianIdealGas
309 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));
310
302 U_inf *= meter / second;
303 delta0 *= meter;
304
305 if (user->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) {
306 PetscReal *mesh_ynodes = NULL;
307 PetscInt mesh_nynodes = 0;
308 if (strcmp(mesh_ynodes_path, "")) PetscCall(GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes));
309 PetscCall(ModifyMesh(comm, dm, problem->dim, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle, &mesh_ynodes, &mesh_nynodes));
310 PetscCall(PetscFree(mesh_ynodes));
311 }
312
313 // Some properties depend on parameters from NewtonianIdealGas
314 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));
315
311 blasius_ctx->weakT = weakT;
312 blasius_ctx->U_inf = U_inf;
313 blasius_ctx->T_inf = T_inf;
314 blasius_ctx->T_wall = T_wall;
315 blasius_ctx->delta0 = delta0;
316 blasius_ctx->P0 = P0;
317 blasius_ctx->n_cheb = N;
318 newtonian_ig_ctx->P0 = P0;
319 blasius_ctx->implicit = user->phys->implicit;
320 blasius_ctx->newtonian_ctx = *newtonian_ig_ctx;
316 StatePrimitive Y_inf = {
317 .pressure = P_inf, .velocity = {U_inf, 0, 0},
318 .temperature = T_inf
319 };
320 State S_infty = StateFromPrimitive(newtonian_ig_ctx, Y_inf);
321
321
322 blasius_ctx->weakT = weakT;
323 blasius_ctx->T_wall = T_wall;
324 blasius_ctx->delta0 = delta0;
325 blasius_ctx->S_infty = S_infty;
326 blasius_ctx->n_cheb = N;
327 newtonian_ig_ctx->idl_pressure = P_inf;
328 blasius_ctx->implicit = user->phys->implicit;
329 blasius_ctx->newtonian_ctx = *newtonian_ig_ctx;
330
322 {
323 PetscReal domain_min[3], domain_max[3];
324 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
325 blasius_ctx->x_inflow = domain_min[0];
326 blasius_ctx->eta_max = 5 * domain_max[1] / blasius_ctx->delta0;
327 }
328 PetscBool diff_filter_mms = PETSC_FALSE;
329 PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter_mms, NULL));
330 if (!use_stg && !diff_filter_mms) PetscCall(ComputeChebyshevCoefficients(blasius_ctx));
331
332 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx));
333
334 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &blasius_context));
335 PetscCallCeed(ceed, CeedQFunctionContextSetData(blasius_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*blasius_ctx), blasius_ctx));
336 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(blasius_context, CEED_MEM_HOST, FreeContextPetsc));
337
338 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context));
339 problem->ics.qfunction_context = blasius_context;
340 if (use_stg) {
331 {
332 PetscReal domain_min[3], domain_max[3];
333 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
334 blasius_ctx->x_inflow = domain_min[0];
335 blasius_ctx->eta_max = 5 * domain_max[1] / blasius_ctx->delta0;
336 }
337 PetscBool diff_filter_mms = PETSC_FALSE;
338 PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter_mms, NULL));
339 if (!use_stg && !diff_filter_mms) PetscCall(ComputeChebyshevCoefficients(blasius_ctx));
340
341 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx));
342
343 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &blasius_context));
344 PetscCallCeed(ceed, CeedQFunctionContextSetData(blasius_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*blasius_ctx), blasius_ctx));
345 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(blasius_context, CEED_MEM_HOST, FreeContextPetsc));
346
347 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context));
348 problem->ics.qfunction_context = blasius_context;
349 if (use_stg) {
341 PetscCall(SetupStg(comm, dm, problem, user, weakT, T_inf, P0));
350 PetscCall(SetupStg(comm, dm, problem, user, weakT, S_infty.Y.temperature, S_infty.Y.pressure));
342 } else if (diff_filter_mms) {
343 PetscCall(DifferentialFilterMmsICSetup(problem));
344 } else {
345 problem->apply_inflow.qfunction = Blasius_Inflow;
346 problem->apply_inflow.qfunction_loc = Blasius_Inflow_loc;
347 problem->apply_inflow_jacobian.qfunction = Blasius_Inflow_Jacobian;
348 problem->apply_inflow_jacobian.qfunction_loc = Blasius_Inflow_Jacobian_loc;
349 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_context, &problem->apply_inflow.qfunction_context));
350 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_context, &problem->apply_inflow_jacobian.qfunction_context));
351 }
352 PetscFunctionReturn(PETSC_SUCCESS);
353}
351 } else if (diff_filter_mms) {
352 PetscCall(DifferentialFilterMmsICSetup(problem));
353 } else {
354 problem->apply_inflow.qfunction = Blasius_Inflow;
355 problem->apply_inflow.qfunction_loc = Blasius_Inflow_loc;
356 problem->apply_inflow_jacobian.qfunction = Blasius_Inflow_Jacobian;
357 problem->apply_inflow_jacobian.qfunction_loc = Blasius_Inflow_Jacobian_loc;
358 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_context, &problem->apply_inflow.qfunction_context));
359 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_context, &problem->apply_inflow_jacobian.qfunction_context));
360 }
361 PetscFunctionReturn(PETSC_SUCCESS);
362}