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 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 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 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, 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, 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. 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