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 PetscErrorCode NS_CHANNEL(ProblemData problem, DM dm, void *ctx, SimpleBC bc) { 28 Honee honee = *(Honee *)ctx; 29 MPI_Comm comm = honee->comm; 30 Ceed ceed = honee->ceed; 31 ChannelContext channel_ctx; 32 NewtonianIdealGasContext newtonian_ig_ctx; 33 CeedQFunctionContext channel_qfctx; 34 PetscBool use_divdiff_verify_mesh = PETSC_FALSE; 35 36 PetscFunctionBeginUser; 37 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx, bc)); 38 PetscCall(PetscCalloc1(1, &channel_ctx)); 39 40 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mesh_transform_channel_div_diff_projection_verify", &use_divdiff_verify_mesh, NULL)); 41 if (use_divdiff_verify_mesh) PetscCall(DivDiffFluxVerifyMesh(dm)); 42 43 // ------------------------------------------------------ 44 // SET UP Channel 45 // ------------------------------------------------------ 46 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx)); 47 problem->ics.qf_func_ptr = ICsChannel; 48 problem->ics.qf_loc = ICsChannel_loc; 49 if (honee->phys->state_var == STATEVAR_CONSERVATIVE) { 50 problem->apply_inflow.qf_func_ptr = Channel_Inflow; 51 problem->apply_inflow.qf_loc = Channel_Inflow_loc; 52 } 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 PetscScalar meter = honee->units->meter; 67 PetscScalar second = honee->units->second; 68 PetscScalar Kelvin = honee->units->Kelvin; 69 PetscScalar Pascal = honee->units->Pascal; 70 71 theta0 *= Kelvin; 72 P0 *= Pascal; 73 umax *= meter / second; 74 75 //-- Setup Problem information 76 CeedScalar H, center; 77 { 78 PetscReal domain_min[3], domain_max[3], domain_size[3]; 79 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 80 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 81 82 H = 0.5 * domain_size[1] * meter; 83 center = H + domain_min[1] * meter; 84 } 85 86 // Some properties depend on parameters from NewtonianIdealGas 87 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 88 89 channel_ctx->center = center; 90 channel_ctx->H = H; 91 channel_ctx->theta0 = theta0; 92 channel_ctx->P0 = P0; 93 channel_ctx->umax = umax; 94 channel_ctx->implicit = honee->phys->implicit; 95 channel_ctx->B = body_force_scale * 2 * umax * newtonian_ig_ctx->mu / (H * H); 96 97 { 98 // Calculate Body force 99 CeedScalar cv = newtonian_ig_ctx->cv, cp = newtonian_ig_ctx->cp; 100 CeedScalar Rd = cp - cv; 101 CeedScalar rho = P0 / (Rd * theta0); 102 CeedScalar g[] = {channel_ctx->B / rho, 0., 0.}; 103 PetscCall(PetscArraycpy(newtonian_ig_ctx->g, g, 3)); 104 } 105 channel_ctx->newtonian_ctx = *newtonian_ig_ctx; 106 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 107 108 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &channel_qfctx)); 109 PetscCallCeed(ceed, CeedQFunctionContextSetData(channel_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*channel_ctx), channel_ctx)); 110 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(channel_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 111 112 problem->ics.qfctx = channel_qfctx; 113 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &problem->apply_inflow.qfctx)); 114 115 for (PetscCount b = 0; b < problem->num_bc_defs; b++) { 116 BCDefinition bc_def = problem->bc_defs[b]; 117 const char *name; 118 119 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL)); 120 if (honee->phys->state_var == STATEVAR_CONSERVATIVE && !strcmp(name, "outflow")) { 121 HoneeBCStruct honee_bc; 122 123 PetscCall(PetscPrintf(comm, "WARNING! Channel flow with Inflow and Outflow is currently broken.\n")); 124 PetscCall(PetscNew(&honee_bc)); 125 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &honee_bc->qfctx)); 126 honee_bc->honee = honee; 127 honee_bc->jac_data_size_sur = honee->phys->implicit ? problem->jac_data_size_sur : 0; 128 PetscCall(BCDefinitionSetContext(bc_def, HoneeBCDestroy, honee_bc)); 129 130 PetscCall(BCDefinitionSetIFunction(bc_def, ChannelOutflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp)); 131 PetscCall(BCDefinitionSetIJacobian(bc_def, NULL, NULL)); 132 } 133 } 134 PetscFunctionReturn(PETSC_SUCCESS); 135 } 136 137 // This function transforms the mesh coordinates to mimic the mesh used in 138 // *A better consistency for low-order stabilized finite element methods* Jansen et. al. 1999 139 // which is used to verify the projection of divergence of diffusive flux. See !27 and !94 for more details. 140 static PetscErrorCode DivDiffFluxVerifyMesh(DM dm) { 141 PetscInt narr, ncoords, dim; 142 PetscReal domain_min[3], domain_max[3], domain_size[3]; 143 PetscScalar *arr_coords; 144 Vec vec_coords; 145 146 PetscFunctionBeginUser; 147 PetscCall(DMGetDimension(dm, &dim)); 148 // Get domain boundary information 149 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 150 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 151 152 // Get coords array from DM 153 PetscCall(DMGetCoordinatesLocal(dm, &vec_coords)); 154 PetscCall(VecGetLocalSize(vec_coords, &narr)); 155 PetscCall(VecGetArray(vec_coords, &arr_coords)); 156 157 PetscScalar(*coords)[dim] = (PetscScalar(*)[dim])arr_coords; 158 ncoords = narr / dim; 159 160 // Get mesh information 161 PetscInt nmax = 3, faces[3]; 162 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL)); 163 // Get element size of the box mesh, for indexing each node 164 const PetscReal dxbox = domain_size[0] / (faces[0]); 165 166 for (PetscInt i = 0; i < ncoords; i++) { 167 PetscInt x_box_index = round(coords[i][0] / dxbox); 168 if (x_box_index % 2) { 169 coords[i][0] = (x_box_index - 1) * dxbox + 0.5 * dxbox; 170 } 171 } 172 173 PetscCall(VecRestoreArray(vec_coords, &arr_coords)); 174 PetscCall(DMSetCoordinatesLocal(dm, vec_coords)); 175 PetscFunctionReturn(0); 176 } 177