1 // Copyright (c) 2017-2026, 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 Channel flow
10
11 #include "../qfunctions/channel.h"
12
13 #include <ceed.h>
14 #include <petscdm.h>
15
16 #include "../navierstokes.h"
17
NS_CHANNEL(ProblemData problem,DM dm,void * ctx,SimpleBC bc)18 PetscErrorCode NS_CHANNEL(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
19 User user = *(User *)ctx;
20 MPI_Comm comm = user->comm;
21 Ceed ceed = user->ceed;
22 ChannelContext channel_ctx;
23 NewtonianIdealGasContext newtonian_ig_ctx;
24 CeedQFunctionContext channel_context;
25
26 PetscFunctionBeginUser;
27 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx, bc));
28 PetscCall(PetscCalloc1(1, &channel_ctx));
29
30 // ------------------------------------------------------
31 // SET UP Channel
32 // ------------------------------------------------------
33 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context));
34 problem->ics.qfunction = ICsChannel;
35 problem->ics.qfunction_loc = ICsChannel_loc;
36 if (user->phys->state_var == STATEVAR_CONSERVATIVE) {
37 problem->apply_inflow.qfunction = Channel_Inflow;
38 problem->apply_inflow.qfunction_loc = Channel_Inflow_loc;
39 problem->apply_outflow.qfunction = Channel_Outflow;
40 problem->apply_outflow.qfunction_loc = Channel_Outflow_loc;
41 }
42
43 // -- Command Line Options
44 CeedScalar umax = 10.; // m/s
45 CeedScalar theta0 = 300.; // K
46 CeedScalar P0 = 1.e5; // Pa
47 PetscReal body_force_scale = 1.;
48 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL);
49 PetscCall(PetscOptionsScalar("-umax", "Centerline velocity of the Channel", NULL, umax, &umax, NULL));
50 PetscCall(PetscOptionsScalar("-theta0", "Wall temperature", NULL, theta0, &theta0, NULL));
51 PetscCall(PetscOptionsScalar("-P0", "Pressure at outflow", NULL, P0, &P0, NULL));
52 PetscCall(PetscOptionsReal("-body_force_scale", "Multiplier for body force", NULL, body_force_scale = 1, &body_force_scale, NULL));
53 PetscOptionsEnd();
54
55 PetscScalar meter = user->units->meter;
56 PetscScalar second = user->units->second;
57 PetscScalar Kelvin = user->units->Kelvin;
58 PetscScalar Pascal = user->units->Pascal;
59
60 theta0 *= Kelvin;
61 P0 *= Pascal;
62 umax *= meter / second;
63
64 //-- Setup Problem information
65 CeedScalar H, center;
66 {
67 PetscReal domain_min[3], domain_max[3], domain_size[3];
68 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
69 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
70
71 H = 0.5 * domain_size[1] * meter;
72 center = H + domain_min[1] * meter;
73 }
74
75 // Some properties depend on parameters from NewtonianIdealGas
76 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));
77
78 channel_ctx->center = center;
79 channel_ctx->H = H;
80 channel_ctx->theta0 = theta0;
81 channel_ctx->P0 = P0;
82 channel_ctx->umax = umax;
83 channel_ctx->implicit = user->phys->implicit;
84 channel_ctx->B = body_force_scale * 2 * umax * newtonian_ig_ctx->mu / (H * H);
85
86 {
87 // Calculate Body force
88 CeedScalar cv = newtonian_ig_ctx->cv, cp = newtonian_ig_ctx->cp;
89 CeedScalar Rd = cp - cv;
90 CeedScalar rho = P0 / (Rd * theta0);
91 CeedScalar g[] = {channel_ctx->B / rho, 0., 0.};
92 PetscCall(PetscArraycpy(newtonian_ig_ctx->g, g, 3));
93 }
94 channel_ctx->newtonian_ctx = *newtonian_ig_ctx;
95 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx));
96
97 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &channel_context));
98 PetscCallCeed(ceed, CeedQFunctionContextSetData(channel_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*channel_ctx), channel_ctx));
99 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(channel_context, CEED_MEM_HOST, FreeContextPetsc));
100
101 problem->ics.qfunction_context = channel_context;
102 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_context, &problem->apply_inflow.qfunction_context));
103 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_context, &problem->apply_outflow.qfunction_context));
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106