1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Utility functions for setting up Channel flow
6
7 #include "../qfunctions/channel.h"
8
9 #include <ceed.h>
10 #include <petscdm.h>
11
12 #include <navierstokes.h>
13
14 static PetscErrorCode DivDiffFluxVerifyMesh(DM dm);
15
ChannelOutflowBCSetup_CreateIFunctionQF(BCDefinition bc_def,CeedQFunction * qf)16 static PetscErrorCode ChannelOutflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) {
17 HoneeBCStruct honee_bc;
18
19 PetscFunctionBeginUser;
20 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
21 PetscCheck(honee_bc->honee->phys->state_var == STATEVAR_CONSERVATIVE, PETSC_COMM_WORLD, PETSC_ERR_SUP,
22 "Channel outflow only valid for Conservative variables, recieved %s", StateVariables[honee_bc->honee->phys->state_var]);
23 PetscCall(HoneeBCCreateIFunctionQF(bc_def, Channel_Outflow, Channel_Outflow_loc, honee_bc->qfctx, qf));
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
ChannelInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def,CeedQFunction * qf)27 static PetscErrorCode ChannelInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) {
28 HoneeBCStruct honee_bc;
29
30 PetscFunctionBeginUser;
31 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
32 PetscCheck(honee_bc->honee->phys->state_var == STATEVAR_CONSERVATIVE, PETSC_COMM_WORLD, PETSC_ERR_SUP,
33 "Channel inflow only valid for Conservative variables, recieved %s", StateVariables[honee_bc->honee->phys->state_var]);
34 PetscCall(HoneeBCCreateIFunctionQF(bc_def, Channel_Inflow, Channel_Inflow_loc, honee_bc->qfctx, qf));
35 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
NS_CHANNEL(ProblemData problem,DM dm,void * ctx)38 PetscErrorCode NS_CHANNEL(ProblemData problem, DM dm, void *ctx) {
39 Honee honee = *(Honee *)ctx;
40 MPI_Comm comm = honee->comm;
41 Ceed ceed = honee->ceed;
42 ChannelContext channel_ctx;
43 NewtonianIdealGasContext newtonian_ig_ctx;
44 CeedQFunctionContext channel_qfctx;
45 PetscBool use_divdiff_verify_mesh = PETSC_FALSE;
46
47 PetscFunctionBeginUser;
48 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx));
49 PetscCall(PetscNew(&channel_ctx));
50
51 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mesh_transform_channel_div_diff_projection_verify", &use_divdiff_verify_mesh, NULL));
52 if (use_divdiff_verify_mesh) PetscCall(DivDiffFluxVerifyMesh(dm));
53
54 // -- Command Line Options
55 CeedScalar umax = 10.; // m/s
56 CeedScalar theta0 = 300.; // K
57 CeedScalar P0 = 1.e5; // Pa
58 PetscReal body_force_scale = 1.;
59 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL);
60 PetscCall(PetscOptionsScalar("-umax", "Centerline velocity of the Channel", NULL, umax, &umax, NULL));
61 PetscCall(PetscOptionsScalar("-theta0", "Wall temperature", NULL, theta0, &theta0, NULL));
62 PetscCall(PetscOptionsScalar("-P0", "Pressure at outflow", NULL, P0, &P0, NULL));
63 PetscCall(PetscOptionsReal("-body_force_scale", "Multiplier for body force", NULL, body_force_scale = 1, &body_force_scale, NULL));
64 PetscOptionsEnd();
65
66 Units units = honee->units;
67
68 theta0 *= units->Kelvin;
69 P0 *= units->Pascal;
70 umax *= units->meter / units->second;
71
72 //-- Setup Problem information
73 CeedScalar H, center;
74 {
75 PetscReal domain_min[3], domain_max[3], domain_size[3];
76 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
77 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
78
79 H = 0.5 * domain_size[1];
80 center = H + domain_min[1];
81 }
82
83 // Some properties depend on parameters from NewtonianIdealGas
84 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx));
85
86 channel_ctx->center = center;
87 channel_ctx->H = H;
88 channel_ctx->theta0 = theta0;
89 channel_ctx->P0 = P0;
90 channel_ctx->umax = umax;
91 channel_ctx->implicit = honee->phys->implicit;
92 channel_ctx->B = body_force_scale * 2 * umax * newtonian_ig_ctx->gas.mu / (H * H);
93
94 {
95 // Calculate Body force
96 CeedScalar cv = newtonian_ig_ctx->gas.cv, cp = newtonian_ig_ctx->gas.cp;
97 CeedScalar Rd = cp - cv;
98 CeedScalar rho = P0 / (Rd * theta0);
99 CeedScalar g[] = {channel_ctx->B / rho, 0., 0.};
100 PetscCall(PetscArraycpy(newtonian_ig_ctx->g, g, 3));
101 }
102 channel_ctx->newt_ctx = *newtonian_ig_ctx;
103 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx));
104
105 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &channel_qfctx));
106 PetscCallCeed(ceed, CeedQFunctionContextSetData(channel_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*channel_ctx), channel_ctx));
107 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(channel_qfctx, CEED_MEM_HOST, FreeContextPetsc));
108
109 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx));
110 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsChannel, .qf_loc = ICsChannel_loc, .qfctx = channel_qfctx};
111
112 for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
113 BCDefinition bc_def = problem->bc_defs[b];
114 const char *name;
115
116 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL));
117 if (honee->phys->state_var == STATEVAR_CONSERVATIVE && !strcmp(name, "outflow")) {
118 HoneeBCStruct honee_bc;
119
120 PetscCall(PetscPrintf(comm, "WARNING! Channel flow with Inflow and Outflow is currently broken.\n"));
121 PetscCall(PetscNew(&honee_bc));
122 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &honee_bc->qfctx));
123 honee_bc->honee = honee;
124 honee_bc->num_comps_jac_data = honee->phys->implicit ? 11 : 0;
125 PetscCall(BCDefinitionSetContext(bc_def, (PetscCtxDestroyFn *)HoneeBCDestroy, honee_bc));
126
127 PetscCall(BCDefinitionSetIFunction(bc_def, ChannelOutflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp));
128 PetscCall(BCDefinitionSetIJacobian(bc_def, NULL, NULL));
129 } else if (honee->phys->state_var == STATEVAR_CONSERVATIVE && !strcmp(name, "inflow")) {
130 HoneeBCStruct honee_bc;
131
132 PetscCall(PetscPrintf(comm, "WARNING! Channel flow with Inflow and Outflow is currently broken.\n"));
133 PetscCall(PetscNew(&honee_bc));
134 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &honee_bc->qfctx));
135 honee_bc->honee = honee;
136 honee_bc->num_comps_jac_data = honee->phys->implicit ? 11 : 0;
137 PetscCall(BCDefinitionSetContext(bc_def, (PetscCtxDestroyFn *)HoneeBCDestroy, honee_bc));
138
139 PetscCall(BCDefinitionSetIFunction(bc_def, ChannelInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp));
140 PetscCall(BCDefinitionSetIJacobian(bc_def, NULL, NULL));
141 }
142 }
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
146 // This function transforms the mesh coordinates to mimic the mesh used in
147 // *A better consistency for low-order stabilized finite element methods* Jansen et. al. 1999
148 // which is used to verify the projection of divergence of diffusive flux. See !27 and !94 for more details.
DivDiffFluxVerifyMesh(DM dm)149 static PetscErrorCode DivDiffFluxVerifyMesh(DM dm) {
150 PetscInt narr, ncoords, dim;
151 PetscReal domain_min[3], domain_max[3], domain_size[3];
152 PetscScalar *arr_coords;
153 Vec vec_coords;
154
155 PetscFunctionBeginUser;
156 PetscCall(DMGetDimension(dm, &dim));
157 // Get domain boundary information
158 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
159 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
160
161 // Get coords array from DM
162 PetscCall(DMGetCoordinatesLocal(dm, &vec_coords));
163 PetscCall(VecGetLocalSize(vec_coords, &narr));
164 PetscCall(VecGetArray(vec_coords, &arr_coords));
165
166 PetscScalar(*coords)[dim] = (PetscScalar(*)[dim])arr_coords;
167 ncoords = narr / dim;
168
169 // Get mesh information
170 PetscInt nmax = 3, faces[3];
171 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL));
172 // Get element size of the box mesh, for indexing each node
173 const PetscReal dxbox = domain_size[0] / (faces[0]);
174
175 for (PetscInt i = 0; i < ncoords; i++) {
176 PetscInt x_box_index = round(coords[i][0] / dxbox);
177 if (x_box_index % 2) {
178 coords[i][0] = (x_box_index - 1) * dxbox + 0.5 * dxbox;
179 }
180 }
181
182 PetscCall(VecRestoreArray(vec_coords, &arr_coords));
183 PetscCall(DMSetCoordinatesLocal(dm, vec_coords));
184 PetscFunctionReturn(0);
185 }
186