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 // ------------------------------------------------------ 55 // SET UP Channel 56 // ------------------------------------------------------ 57 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx)); 58 problem->ics.qf_func_ptr = ICsChannel; 59 problem->ics.qf_loc = ICsChannel_loc; 60 61 // -- Command Line Options 62 CeedScalar umax = 10.; // m/s 63 CeedScalar theta0 = 300.; // K 64 CeedScalar P0 = 1.e5; // Pa 65 PetscReal body_force_scale = 1.; 66 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL); 67 PetscCall(PetscOptionsScalar("-umax", "Centerline velocity of the Channel", NULL, umax, &umax, NULL)); 68 PetscCall(PetscOptionsScalar("-theta0", "Wall temperature", NULL, theta0, &theta0, NULL)); 69 PetscCall(PetscOptionsScalar("-P0", "Pressure at outflow", NULL, P0, &P0, NULL)); 70 PetscCall(PetscOptionsReal("-body_force_scale", "Multiplier for body force", NULL, body_force_scale = 1, &body_force_scale, NULL)); 71 PetscOptionsEnd(); 72 73 Units units = honee->units; 74 75 theta0 *= units->Kelvin; 76 P0 *= units->Pascal; 77 umax *= units->meter / units->second; 78 79 //-- Setup Problem information 80 CeedScalar H, center; 81 { 82 PetscReal domain_min[3], domain_max[3], domain_size[3]; 83 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 84 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 85 86 H = 0.5 * domain_size[1] * units->meter; 87 center = H + domain_min[1] * units->meter; 88 } 89 90 // Some properties depend on parameters from NewtonianIdealGas 91 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 92 93 channel_ctx->center = center; 94 channel_ctx->H = H; 95 channel_ctx->theta0 = theta0; 96 channel_ctx->P0 = P0; 97 channel_ctx->umax = umax; 98 channel_ctx->implicit = honee->phys->implicit; 99 channel_ctx->B = body_force_scale * 2 * umax * newtonian_ig_ctx->mu / (H * H); 100 101 { 102 // Calculate Body force 103 CeedScalar cv = newtonian_ig_ctx->cv, cp = newtonian_ig_ctx->cp; 104 CeedScalar Rd = cp - cv; 105 CeedScalar rho = P0 / (Rd * theta0); 106 CeedScalar g[] = {channel_ctx->B / rho, 0., 0.}; 107 PetscCall(PetscArraycpy(newtonian_ig_ctx->g, g, 3)); 108 } 109 channel_ctx->newtonian_ctx = *newtonian_ig_ctx; 110 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 111 112 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &channel_qfctx)); 113 PetscCallCeed(ceed, CeedQFunctionContextSetData(channel_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*channel_ctx), channel_ctx)); 114 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(channel_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 115 problem->ics.qfctx = channel_qfctx; 116 117 for (PetscCount b = 0; b < problem->num_bc_defs; b++) { 118 BCDefinition bc_def = problem->bc_defs[b]; 119 const char *name; 120 121 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL)); 122 if (honee->phys->state_var == STATEVAR_CONSERVATIVE && !strcmp(name, "outflow")) { 123 HoneeBCStruct honee_bc; 124 125 PetscCall(PetscPrintf(comm, "WARNING! Channel flow with Inflow and Outflow is currently broken.\n")); 126 PetscCall(PetscNew(&honee_bc)); 127 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &honee_bc->qfctx)); 128 honee_bc->honee = honee; 129 honee_bc->num_comps_jac_data = honee->phys->implicit ? 11 : 0; 130 PetscCall(BCDefinitionSetContext(bc_def, HoneeBCDestroy, honee_bc)); 131 132 PetscCall(BCDefinitionSetIFunction(bc_def, ChannelOutflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp)); 133 PetscCall(BCDefinitionSetIJacobian(bc_def, NULL, NULL)); 134 } else if (honee->phys->state_var == STATEVAR_CONSERVATIVE && !strcmp(name, "inflow")) { 135 HoneeBCStruct honee_bc; 136 137 PetscCall(PetscPrintf(comm, "WARNING! Channel flow with Inflow and Outflow is currently broken.\n")); 138 PetscCall(PetscNew(&honee_bc)); 139 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &honee_bc->qfctx)); 140 honee_bc->honee = honee; 141 honee_bc->num_comps_jac_data = honee->phys->implicit ? 11 : 0; 142 PetscCall(BCDefinitionSetContext(bc_def, HoneeBCDestroy, honee_bc)); 143 144 PetscCall(BCDefinitionSetIFunction(bc_def, ChannelInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp)); 145 PetscCall(BCDefinitionSetIJacobian(bc_def, NULL, NULL)); 146 } 147 } 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 // This function transforms the mesh coordinates to mimic the mesh used in 152 // *A better consistency for low-order stabilized finite element methods* Jansen et. al. 1999 153 // which is used to verify the projection of divergence of diffusive flux. See !27 and !94 for more details. 154 static PetscErrorCode DivDiffFluxVerifyMesh(DM dm) { 155 PetscInt narr, ncoords, dim; 156 PetscReal domain_min[3], domain_max[3], domain_size[3]; 157 PetscScalar *arr_coords; 158 Vec vec_coords; 159 160 PetscFunctionBeginUser; 161 PetscCall(DMGetDimension(dm, &dim)); 162 // Get domain boundary information 163 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 164 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 165 166 // Get coords array from DM 167 PetscCall(DMGetCoordinatesLocal(dm, &vec_coords)); 168 PetscCall(VecGetLocalSize(vec_coords, &narr)); 169 PetscCall(VecGetArray(vec_coords, &arr_coords)); 170 171 PetscScalar(*coords)[dim] = (PetscScalar(*)[dim])arr_coords; 172 ncoords = narr / dim; 173 174 // Get mesh information 175 PetscInt nmax = 3, faces[3]; 176 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL)); 177 // Get element size of the box mesh, for indexing each node 178 const PetscReal dxbox = domain_size[0] / (faces[0]); 179 180 for (PetscInt i = 0; i < ncoords; i++) { 181 PetscInt x_box_index = round(coords[i][0] / dxbox); 182 if (x_box_index % 2) { 183 coords[i][0] = (x_box_index - 1) * dxbox + 0.5 * dxbox; 184 } 185 } 186 187 PetscCall(VecRestoreArray(vec_coords, &arr_coords)); 188 PetscCall(DMSetCoordinatesLocal(dm, vec_coords)); 189 PetscFunctionReturn(0); 190 } 191