1 // SPDX-FileCopyrightText: Copyright (c) 2017-2025, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 #include <honee.h> 5 #include <petscdmplex.h> 6 7 #ifndef M_PI 8 #define M_PI 3.14159265358979323846 9 #endif 10 11 static PetscErrorCode GetYNodeLocs(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], PetscReal **pynodes, PetscInt *nynodes) { 12 int ndims, dims[2] = {0.}; 13 FILE *fp; 14 const PetscInt char_array_len = 512; 15 char line[char_array_len]; 16 PetscReal *node_locs; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscFOpen(comm, path, "r", &fp)); 20 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 21 { 22 char **array; 23 24 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 25 for (PetscInt i = 0; i < ndims; i++) dims[i] = atoi(array[i]); 26 PetscCall(PetscStrToArrayDestroy(ndims, array)); 27 } 28 if (ndims < 2) dims[1] = 1; // Assume 1 column of data is not otherwise specified 29 *nynodes = dims[0]; 30 PetscCall(PetscMalloc1(*nynodes, &node_locs)); 31 32 for (PetscInt i = 0; i < dims[0]; i++) { 33 char **array; 34 35 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 36 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 37 PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED, 38 "Line %" PetscInt_FMT " of %s does not contain correct number of columns (%d instead of %d)", i, path, ndims, dims[1]); 39 40 node_locs[i] = (PetscReal)atof(array[0]); 41 PetscCall(PetscStrToArrayDestroy(ndims, array)); 42 } 43 PetscCall(PetscFClose(comm, fp)); 44 *pynodes = node_locs; 45 PetscFunctionReturn(PETSC_SUCCESS); 46 } 47 48 /* @brief Modify the domain and mesh for blasius 49 * 50 * Modifies mesh such that `N` elements are within `refine_height` with a geometric growth ratio of `growth`. Excess elements are then distributed 51 * linearly in logspace to the top surface. 52 * 53 * The top surface is also angled downwards, so that it may be used as an outflow. 54 * It's angle is controlled by `top_angle` (in units of degrees). 55 * 56 * If `node_locs` is not NULL, then the nodes will be placed at `node_locs` locations. 57 * If it is NULL, then the modified coordinate values will be set in the array, along with `num_node_locs`. 58 */ 59 static PetscErrorCode HoneeMeshTransform_PlateMesh(MPI_Comm comm, DM dm, PetscReal growth, PetscInt N, PetscReal refine_height, PetscReal top_angle, 60 PetscReal *node_locs[], PetscInt *num_node_locs) { 61 PetscInt narr, ncoords, dim; 62 PetscReal domain_min[3], domain_max[3], domain_size[3]; 63 PetscScalar *arr_coords; 64 Vec vec_coords; 65 66 PetscFunctionBeginUser; 67 PetscCall(DMGetDimension(dm, &dim)); 68 PetscReal angle_coeff = tan(top_angle * (M_PI / 180)); 69 // Get domain boundary information 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 // Get coords array from DM 74 PetscCall(DMGetCoordinatesLocal(dm, &vec_coords)); 75 PetscCall(VecGetLocalSize(vec_coords, &narr)); 76 PetscCall(VecGetArray(vec_coords, &arr_coords)); 77 78 PetscScalar(*coords)[dim] = (PetscScalar(*)[dim])arr_coords; 79 ncoords = narr / dim; 80 81 // Get mesh information 82 PetscInt nmax = 3, faces[3]; 83 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL)); 84 // Get element size of the box mesh, for indexing each node 85 const PetscReal dybox = domain_size[1] / faces[1]; 86 87 if (!*node_locs) { 88 // Calculate the first element height 89 PetscReal dy1 = refine_height * (growth - 1) / (pow(growth, N) - 1); 90 91 // Calculate log of sizing outside BL 92 PetscReal logdy = (log(domain_max[1]) - log(refine_height)) / (faces[1] - N); 93 94 *num_node_locs = faces[1] + 1; 95 PetscReal *temp_node_locs; 96 PetscCall(PetscMalloc1(*num_node_locs, &temp_node_locs)); 97 98 for (PetscInt i = 0; i < ncoords; i++) { 99 PetscInt y_box_index = round(coords[i][1] / dybox); 100 if (y_box_index <= N) { 101 coords[i][1] = 102 (1 - (coords[i][0] - domain_min[0]) * angle_coeff / domain_max[1]) * dy1 * (pow(growth, coords[i][1] / dybox) - 1) / (growth - 1); 103 } else { 104 PetscInt j = y_box_index - N; 105 coords[i][1] = (1 - (coords[i][0] - domain_min[0]) * angle_coeff / domain_max[1]) * exp(log(refine_height) + logdy * j); 106 } 107 if (coords[i][0] == domain_min[0] && coords[i][2] == domain_min[2]) temp_node_locs[y_box_index] = coords[i][1]; 108 } 109 110 *node_locs = temp_node_locs; 111 } else { 112 PetscCheck(*num_node_locs >= faces[1] + 1, comm, PETSC_ERR_FILE_UNEXPECTED, 113 "The y_node_locs_path has too few locations; There are %" PetscInt_FMT " + 1 nodes, but only %" PetscInt_FMT " locations given", 114 faces[1] + 1, *num_node_locs); 115 if (*num_node_locs > faces[1] + 1) { 116 PetscCall(PetscPrintf(comm, 117 "WARNING: y_node_locs_path has more locations (%" PetscInt_FMT ") " 118 "than the mesh has nodes (%" PetscInt_FMT "). This maybe unintended.\n", 119 *num_node_locs, faces[1] + 1)); 120 } 121 PetscScalar max_y = (*node_locs)[faces[1]]; 122 123 for (PetscInt i = 0; i < ncoords; i++) { 124 // Determine which y-node we're at 125 PetscInt y_box_index = round(coords[i][1] / dybox); 126 coords[i][1] = (1 - (coords[i][0] - domain_min[0]) * angle_coeff / max_y) * (*node_locs)[y_box_index]; 127 } 128 } 129 130 PetscCall(VecRestoreArray(vec_coords, &arr_coords)); 131 PetscCall(DMSetCoordinatesLocal(dm, vec_coords)); 132 PetscFunctionReturn(PETSC_SUCCESS); 133 } 134 135 /** 136 @brief Transform the mesh based on pass options 137 138 Use the option `-mesh_transform` to determine the type of mesh transformation to apply 139 140 @param dm `DM` of mesh to be transformed 141 **/ 142 PetscErrorCode HoneeMeshTransformFromOptions(DM dm) { 143 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 144 145 PetscFunctionBeginUser; 146 PetscOptionsBegin(comm, NULL, "Options for HONEE Mesh Transformation", NULL); 147 MeshTransform mesh_transform_type = MESH_TRANSFORM_NONE; 148 PetscCall(PetscOptionsDeprecated("-mesh_transform", "-meshtransform_type", "HONEE 0.0", NULL)); 149 PetscCall(PetscOptionsEnum("-meshtransform_type", "Mesh transform to perform", NULL, MeshTransforms, (PetscEnum)mesh_transform_type, 150 (PetscEnum *)&mesh_transform_type, NULL)); 151 152 switch (mesh_transform_type) { 153 case MESH_TRANSFORM_PLATEMESH: { 154 PetscReal mesh_refine_height = 5.9e-4; // m 155 PetscReal mesh_growth = 1.08; // [-] 156 PetscInt mesh_Ndelta = 45; // [-] 157 PetscReal mesh_top_angle = 5; // degrees 158 char mesh_ynodes_path[PETSC_MAX_PATH_LEN] = ""; 159 160 PetscCall(PetscOptionsDeprecated("-platemesh_Ndelta", "-meshtransform_platemesh_Ndelta", "HONEE 0.0", NULL)); 161 PetscCall(PetscOptionsDeprecated("-platemesh_refine_height", "-meshtransform_platemesh_refine_height", "HONEE 0.0", NULL)); 162 PetscCall(PetscOptionsDeprecated("-platemesh_growth", "-meshtransform_platemesh_growth", "HONEE 0.0", NULL)); 163 PetscCall(PetscOptionsDeprecated("-platemesh_top_angle", "-meshtransform_platemesh_top_angle", "HONEE 0.0", NULL)); 164 PetscCall(PetscOptionsDeprecated("-platemesh_y_node_locs_path", "-meshtransform_platemesh_y_node_locs_path", "HONEE 0.0", NULL)); 165 166 PetscCall(PetscOptionsBoundedInt("-meshtransform_platemesh_Ndelta", "Velocity at boundary layer edge", NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 167 1)); 168 PetscCall(PetscOptionsScalar("-meshtransform_platemesh_refine_height", "Height of boundary layer mesh refinement", NULL, mesh_refine_height, 169 &mesh_refine_height, NULL)); 170 PetscCall(PetscOptionsScalar("-meshtransform_platemesh_growth", "Geometric growth rate of boundary layer mesh", NULL, mesh_growth, &mesh_growth, 171 NULL)); 172 PetscCall(PetscOptionsScalar("-meshtransform_platemesh_top_angle", "Geometric top_angle rate of boundary layer mesh", NULL, mesh_top_angle, 173 &mesh_top_angle, NULL)); 174 PetscCall(PetscOptionsString("-meshtransform_platemesh_y_node_locs_path", 175 "Path to file with y node locations. " 176 "If empty, will use the algorithmic mesh warping.", 177 NULL, mesh_ynodes_path, mesh_ynodes_path, sizeof(mesh_ynodes_path), NULL)); 178 PetscReal *mesh_ynodes = NULL; 179 PetscInt mesh_nynodes = 0; 180 if (strcmp(mesh_ynodes_path, "")) PetscCall(GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes)); 181 PetscCall(HoneeMeshTransform_PlateMesh(comm, dm, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle, &mesh_ynodes, &mesh_nynodes)); 182 PetscCall(PetscFree(mesh_ynodes)); 183 } break; 184 case MESH_TRANSFORM_NONE: 185 break; 186 } 187 PetscOptionsEnd(); 188 189 { // Scale coordinates based on DM length scale 190 // We want to make sure to do the scaling *after* the mesh transformation, so units in the transform are physical 191 Vec coords; 192 PetscReal meter; 193 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &meter)); 194 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 195 PetscCall(VecScale(coords, meter)); 196 PetscCall(DMSetCoordinatesLocal(dm, coords)); 197 } 198 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201