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 PetscErrorCode NS_CHANNEL(ProblemData problem, DM dm, void *ctx, SimpleBC bc) { 17 Honee honee = *(Honee *)ctx; 18 MPI_Comm comm = honee->comm; 19 Ceed ceed = honee->ceed; 20 ChannelContext channel_ctx; 21 NewtonianIdealGasContext newtonian_ig_ctx; 22 CeedQFunctionContext channel_qfctx; 23 PetscBool use_divdiff_verify_mesh = PETSC_FALSE; 24 25 PetscFunctionBeginUser; 26 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx, bc)); 27 PetscCall(PetscCalloc1(1, &channel_ctx)); 28 29 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mesh_transform_channel_div_diff_projection_verify", &use_divdiff_verify_mesh, NULL)); 30 if (use_divdiff_verify_mesh) PetscCall(DivDiffFluxVerifyMesh(dm)); 31 32 // ------------------------------------------------------ 33 // SET UP Channel 34 // ------------------------------------------------------ 35 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx)); 36 problem->ics.qf_func_ptr = ICsChannel; 37 problem->ics.qf_loc = ICsChannel_loc; 38 if (honee->phys->state_var == STATEVAR_CONSERVATIVE) { 39 problem->apply_inflow.qf_func_ptr = Channel_Inflow; 40 problem->apply_inflow.qf_loc = Channel_Inflow_loc; 41 problem->apply_outflow.qf_func_ptr = Channel_Outflow; 42 problem->apply_outflow.qf_loc = Channel_Outflow_loc; 43 } 44 45 // -- Command Line Options 46 CeedScalar umax = 10.; // m/s 47 CeedScalar theta0 = 300.; // K 48 CeedScalar P0 = 1.e5; // Pa 49 PetscReal body_force_scale = 1.; 50 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL); 51 PetscCall(PetscOptionsScalar("-umax", "Centerline velocity of the Channel", NULL, umax, &umax, NULL)); 52 PetscCall(PetscOptionsScalar("-theta0", "Wall temperature", NULL, theta0, &theta0, NULL)); 53 PetscCall(PetscOptionsScalar("-P0", "Pressure at outflow", NULL, P0, &P0, NULL)); 54 PetscCall(PetscOptionsReal("-body_force_scale", "Multiplier for body force", NULL, body_force_scale = 1, &body_force_scale, NULL)); 55 PetscOptionsEnd(); 56 57 PetscScalar meter = honee->units->meter; 58 PetscScalar second = honee->units->second; 59 PetscScalar Kelvin = honee->units->Kelvin; 60 PetscScalar Pascal = honee->units->Pascal; 61 62 theta0 *= Kelvin; 63 P0 *= Pascal; 64 umax *= meter / second; 65 66 //-- Setup Problem information 67 CeedScalar H, center; 68 { 69 PetscReal domain_min[3], domain_max[3], domain_size[3]; 70 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 71 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 72 73 H = 0.5 * domain_size[1] * meter; 74 center = H + domain_min[1] * meter; 75 } 76 77 // Some properties depend on parameters from NewtonianIdealGas 78 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 79 80 channel_ctx->center = center; 81 channel_ctx->H = H; 82 channel_ctx->theta0 = theta0; 83 channel_ctx->P0 = P0; 84 channel_ctx->umax = umax; 85 channel_ctx->implicit = honee->phys->implicit; 86 channel_ctx->B = body_force_scale * 2 * umax * newtonian_ig_ctx->mu / (H * H); 87 88 { 89 // Calculate Body force 90 CeedScalar cv = newtonian_ig_ctx->cv, cp = newtonian_ig_ctx->cp; 91 CeedScalar Rd = cp - cv; 92 CeedScalar rho = P0 / (Rd * theta0); 93 CeedScalar g[] = {channel_ctx->B / rho, 0., 0.}; 94 PetscCall(PetscArraycpy(newtonian_ig_ctx->g, g, 3)); 95 } 96 channel_ctx->newtonian_ctx = *newtonian_ig_ctx; 97 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 98 99 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &channel_qfctx)); 100 PetscCallCeed(ceed, CeedQFunctionContextSetData(channel_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*channel_ctx), channel_ctx)); 101 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(channel_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 102 103 problem->ics.qfctx = channel_qfctx; 104 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &problem->apply_inflow.qfctx)); 105 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(channel_qfctx, &problem->apply_outflow.qfctx)); 106 PetscFunctionReturn(PETSC_SUCCESS); 107 } 108 109 // This function transforms the mesh coordinates to mimic the mesh used in 110 // *A better consistency for low-order stabilized finite element methods* Jansen et. al. 1999 111 // which is used to verify the projection of divergence of diffusive flux. See !27 and !94 for more details. 112 static PetscErrorCode DivDiffFluxVerifyMesh(DM dm) { 113 PetscInt narr, ncoords, dim; 114 PetscReal domain_min[3], domain_max[3], domain_size[3]; 115 PetscScalar *arr_coords; 116 Vec vec_coords; 117 118 PetscFunctionBeginUser; 119 PetscCall(DMGetDimension(dm, &dim)); 120 // Get domain boundary information 121 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 122 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 123 124 // Get coords array from DM 125 PetscCall(DMGetCoordinatesLocal(dm, &vec_coords)); 126 PetscCall(VecGetLocalSize(vec_coords, &narr)); 127 PetscCall(VecGetArray(vec_coords, &arr_coords)); 128 129 PetscScalar(*coords)[dim] = (PetscScalar(*)[dim])arr_coords; 130 ncoords = narr / dim; 131 132 // Get mesh information 133 PetscInt nmax = 3, faces[3]; 134 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL)); 135 // Get element size of the box mesh, for indexing each node 136 const PetscReal dxbox = domain_size[0] / (faces[0]); 137 138 for (PetscInt i = 0; i < ncoords; i++) { 139 PetscInt x_box_index = round(coords[i][0] / dxbox); 140 if (x_box_index % 2) { 141 coords[i][0] = (x_box_index - 1) * dxbox + 0.5 * dxbox; 142 } 143 } 144 145 PetscCall(VecRestoreArray(vec_coords, &arr_coords)); 146 PetscCall(DMSetCoordinatesLocal(dm, vec_coords)); 147 PetscFunctionReturn(0); 148 } 149