1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Utility functions for setting up Blasius Boundary Layer 10 11 #include "../navierstokes.h" 12 #include "../qfunctions/blasius.h" 13 #include "stg_shur14.h" 14 15 static PetscErrorCode GetYNodeLocs(const MPI_Comm comm, 16 const char path[PETSC_MAX_PATH_LEN], PetscReal **pynodes, 17 PetscInt *nynodes) { 18 PetscErrorCode ierr; 19 PetscInt ndims, dims[2]; 20 FILE *fp; 21 const PetscInt char_array_len = 512; 22 char line[char_array_len]; 23 char **array; 24 PetscReal *node_locs; 25 PetscFunctionBeginUser; 26 27 ierr = PetscFOpen(comm, path, "r", &fp); CHKERRQ(ierr); 28 ierr = PetscSynchronizedFGets(comm, fp, char_array_len, line); CHKERRQ(ierr); 29 ierr = PetscStrToArray(line, ' ', &ndims, &array); CHKERRQ(ierr); 30 31 for (PetscInt i=0; i<ndims; i++) dims[i] = atoi(array[i]); 32 if (ndims<2) dims[1] = 1; // Assume 1 column of data is not otherwise specified 33 *nynodes = dims[0]; 34 ierr = PetscMalloc1(*nynodes, &node_locs); CHKERRQ(ierr); 35 36 for (PetscInt i=0; i<dims[0]; i++) { 37 ierr = PetscSynchronizedFGets(comm, fp, char_array_len, line); CHKERRQ(ierr); 38 ierr = PetscStrToArray(line, ' ', &ndims, &array); CHKERRQ(ierr); 39 if (ndims < dims[1]) SETERRQ(comm, -1, 40 "Line %" PetscInt_FMT" of %s does not contain enough columns (%" 41 PetscInt_FMT" instead of %" PetscInt_FMT ")", i, 42 path, ndims, dims[1]); 43 44 node_locs[i] = (PetscReal) atof(array[0]); 45 } 46 ierr = PetscFClose(comm, fp); CHKERRQ(ierr); 47 *pynodes = node_locs; 48 PetscFunctionReturn(0); 49 } 50 51 /* \brief Modify the domain and mesh for blasius 52 * 53 * Modifies mesh such that `N` elements are within `refine_height` with a 54 * geometric growth ratio of `growth`. Excess elements are then distributed 55 * linearly in logspace to the top surface. 56 * 57 * The top surface is also angled downwards, so that it may be used as an 58 * outflow. It's angle is controlled by `top_angle` (in units of degrees). 59 * 60 * If `node_locs` is not NULL, then the nodes will be placed at `node_locs` 61 * locations. If it is NULL, then the modified coordinate values will be set in 62 * the array, along with `num_node_locs`. 63 */ 64 static PetscErrorCode ModifyMesh(MPI_Comm comm, DM dm, PetscInt dim, 65 PetscReal growth, PetscInt N, 66 PetscReal refine_height, PetscReal top_angle, 67 PetscReal *node_locs[], PetscInt *num_node_locs) { 68 PetscInt ierr, narr, ncoords; 69 PetscReal domain_min[3], domain_max[3], domain_size[3]; 70 PetscScalar *arr_coords; 71 Vec vec_coords; 72 PetscFunctionBeginUser; 73 74 PetscReal angle_coeff = tan(top_angle*(M_PI/180)); 75 76 // Get domain boundary information 77 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 78 for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 79 80 // Get coords array from DM 81 ierr = DMGetCoordinatesLocal(dm, &vec_coords); CHKERRQ(ierr); 82 ierr = VecGetLocalSize(vec_coords, &narr); CHKERRQ(ierr); 83 ierr = VecGetArray(vec_coords, &arr_coords); CHKERRQ(ierr); 84 85 PetscScalar (*coords)[dim] = (PetscScalar(*)[dim]) arr_coords; 86 ncoords = narr/dim; 87 88 // Get mesh information 89 PetscInt nmax = 3, faces[3]; 90 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, 91 NULL); CHKERRQ(ierr); 92 // Get element size of the box mesh, for indexing each node 93 const PetscReal dybox = domain_size[1]/faces[1]; 94 95 if (!*node_locs) { 96 // Calculate the first element height 97 PetscReal dy1 = refine_height*(growth-1)/(pow(growth, N)-1); 98 99 // Calculate log of sizing outside BL 100 PetscReal logdy = (log(domain_max[1]) - log(refine_height)) / (faces[1] - N); 101 102 *num_node_locs = faces[1] + 1; 103 PetscReal *temp_node_locs; 104 ierr = PetscMalloc1(*num_node_locs, &temp_node_locs); CHKERRQ(ierr); 105 106 for (PetscInt i=0; i<ncoords; i++) { 107 PetscInt y_box_index = round(coords[i][1]/dybox); 108 if (y_box_index <= N) { 109 coords[i][1] = (1 - (coords[i][0] - domain_min[0])*angle_coeff/domain_max[1]) 110 * dy1 * (pow(growth, coords[i][1]/dybox)-1)/(growth-1); 111 } else { 112 PetscInt j = y_box_index - N; 113 coords[i][1] = (1 - (coords[i][0] - domain_min[0])*angle_coeff/domain_max[1]) 114 * exp(log(refine_height) + logdy*j); 115 } 116 if (coords[i][0] == domain_min[0] && coords[i][2] == domain_min[2]) 117 temp_node_locs[y_box_index] = coords[i][1]; 118 } 119 120 *node_locs = temp_node_locs; 121 } else { 122 // Error checking 123 if (*num_node_locs < faces[1] +1) 124 SETERRQ(comm, -1, "The y_node_locs_path has too few locations; " 125 "There are %d + 1 nodes, but only %d locations given", 126 faces[1]+1, *num_node_locs); 127 if (*num_node_locs > faces[1] +1) { 128 ierr = PetscPrintf(comm, "WARNING: y_node_locs_path has more locations (%d) " 129 "than the mesh has nodes (%d). This maybe unintended.\n", 130 *num_node_locs, faces[1]+1); CHKERRQ(ierr); 131 } 132 PetscScalar max_y = (*node_locs)[faces[1]]; 133 134 for (PetscInt i=0; i<ncoords; i++) { 135 // Determine which y-node we're at 136 PetscInt y_box_index = round(coords[i][1]/dybox); 137 coords[i][1] = (1 - (coords[i][0] - domain_min[0])*angle_coeff/max_y) 138 * (*node_locs)[y_box_index]; 139 } 140 } 141 142 ierr = VecRestoreArray(vec_coords, &arr_coords); CHKERRQ(ierr); 143 ierr = DMSetCoordinatesLocal(dm, vec_coords); CHKERRQ(ierr); 144 145 PetscFunctionReturn(0); 146 } 147 148 PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx) { 149 150 PetscInt ierr; 151 User user = *(User *)ctx; 152 MPI_Comm comm = PETSC_COMM_WORLD; 153 PetscBool use_stg = PETSC_FALSE; 154 BlasiusContext blasius_ctx; 155 NewtonianIdealGasContext newtonian_ig_ctx; 156 CeedQFunctionContext blasius_context; 157 158 PetscFunctionBeginUser; 159 ierr = NS_NEWTONIAN_IG(problem, dm, ctx); CHKERRQ(ierr); 160 ierr = PetscCalloc1(1, &blasius_ctx); CHKERRQ(ierr); 161 162 // ------------------------------------------------------ 163 // SET UP Blasius 164 // ------------------------------------------------------ 165 problem->ics.qfunction = ICsBlasius; 166 problem->ics.qfunction_loc = ICsBlasius_loc; 167 168 CeedScalar Uinf = 40; // m/s 169 CeedScalar delta0 = 4.2e-4; // m 170 CeedScalar theta0 = 288.; // K 171 CeedScalar P0 = 1.01e5; // Pa 172 PetscBool weakT = PETSC_FALSE; // weak density or temperature 173 PetscReal mesh_refine_height = 5.9e-4; // m 174 PetscReal mesh_growth = 1.08; // [-] 175 PetscInt mesh_Ndelta = 45; // [-] 176 PetscReal mesh_top_angle = 5; // degrees 177 char mesh_ynodes_path[PETSC_MAX_PATH_LEN] = ""; 178 179 PetscOptionsBegin(comm, NULL, "Options for BLASIUS problem", NULL); 180 ierr = PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", 181 NULL, weakT, &weakT, NULL); CHKERRQ(ierr); 182 ierr = PetscOptionsScalar("-Uinf", "Velocity at boundary layer edge", 183 NULL, Uinf, &Uinf, NULL); CHKERRQ(ierr); 184 ierr = PetscOptionsScalar("-delta0", "Boundary layer height at inflow", 185 NULL, delta0, &delta0, NULL); CHKERRQ(ierr); 186 ierr = PetscOptionsScalar("-theta0", "Wall temperature", 187 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 188 ierr = PetscOptionsScalar("-P0", "Pressure at outflow", 189 NULL, P0, &P0, NULL); CHKERRQ(ierr); 190 ierr = PetscOptionsBoundedInt("-platemesh_Ndelta", 191 "Velocity at boundary layer edge", 192 NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1); CHKERRQ(ierr); 193 ierr = PetscOptionsScalar("-platemesh_refine_height", 194 "Height of boundary layer mesh refinement", 195 NULL, mesh_refine_height, &mesh_refine_height, NULL); CHKERRQ(ierr); 196 ierr = PetscOptionsScalar("-platemesh_growth", 197 "Geometric growth rate of boundary layer mesh", 198 NULL, mesh_growth, &mesh_growth, NULL); CHKERRQ(ierr); 199 ierr = PetscOptionsScalar("-platemesh_top_angle", 200 "Geometric top_angle rate of boundary layer mesh", 201 NULL, mesh_top_angle, &mesh_top_angle, NULL); CHKERRQ(ierr); 202 ierr = PetscOptionsString("-platemesh_y_node_locs_path", 203 "Path to file with y node locations. " 204 "If empty, will use the algorithmic mesh warping.", NULL, 205 mesh_ynodes_path, mesh_ynodes_path, 206 sizeof(mesh_ynodes_path), NULL); CHKERRQ(ierr); 207 ierr = PetscOptionsBool("-stg_use", "Use STG inflow boundary condition", 208 NULL, use_stg, &use_stg, NULL); CHKERRQ(ierr); 209 PetscOptionsEnd(); 210 211 PetscScalar meter = user->units->meter; 212 PetscScalar second = user->units->second; 213 PetscScalar Kelvin = user->units->Kelvin; 214 PetscScalar Pascal = user->units->Pascal; 215 216 theta0 *= Kelvin; 217 P0 *= Pascal; 218 Uinf *= meter / second; 219 delta0 *= meter; 220 221 PetscReal *mesh_ynodes = NULL; 222 PetscInt mesh_nynodes = 0; 223 if (strcmp(mesh_ynodes_path, "")) { 224 ierr = GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes); 225 CHKERRQ(ierr); 226 } 227 ierr = ModifyMesh(comm, dm, problem->dim, mesh_growth, mesh_Ndelta, 228 mesh_refine_height, mesh_top_angle, &mesh_ynodes, 229 &mesh_nynodes); CHKERRQ(ierr); 230 231 // Some properties depend on parameters from NewtonianIdealGas 232 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, 233 CEED_MEM_HOST, &newtonian_ig_ctx); 234 235 blasius_ctx->weakT = weakT; 236 blasius_ctx->Uinf = Uinf; 237 blasius_ctx->delta0 = delta0; 238 blasius_ctx->theta0 = theta0; 239 blasius_ctx->P0 = P0; 240 newtonian_ig_ctx->P0 = P0; 241 blasius_ctx->implicit = user->phys->implicit; 242 blasius_ctx->newtonian_ctx = *newtonian_ig_ctx; 243 244 { 245 PetscReal domain_min[3]; 246 ierr = DMGetBoundingBox(dm, domain_min, NULL); CHKERRQ(ierr); 247 blasius_ctx->x_inflow = domain_min[0]; 248 } 249 250 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, 251 &newtonian_ig_ctx); 252 253 CeedQFunctionContextCreate(user->ceed, &blasius_context); 254 CeedQFunctionContextSetData(blasius_context, CEED_MEM_HOST, 255 CEED_USE_POINTER, 256 sizeof(*blasius_ctx), blasius_ctx); 257 CeedQFunctionContextSetDataDestroy(blasius_context, CEED_MEM_HOST, 258 FreeContextPetsc); 259 260 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 261 problem->ics.qfunction_context = blasius_context; 262 if (use_stg) { 263 ierr = SetupSTG(comm, dm, problem, user, weakT, theta0, P0, mesh_ynodes, 264 mesh_nynodes); CHKERRQ(ierr); 265 } else { 266 problem->apply_inflow.qfunction = Blasius_Inflow; 267 problem->apply_inflow.qfunction_loc = Blasius_Inflow_loc; 268 problem->apply_inflow_jacobian.qfunction = Blasius_Inflow_Jacobian; 269 problem->apply_inflow_jacobian.qfunction_loc = Blasius_Inflow_Jacobian_loc; 270 CeedQFunctionContextReferenceCopy(blasius_context, 271 &problem->apply_inflow.qfunction_context); 272 CeedQFunctionContextReferenceCopy(blasius_context, 273 &problem->apply_inflow_jacobian.qfunction_context); 274 } 275 ierr = PetscFree(mesh_ynodes); CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278