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