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
GetYNodeLocs(const MPI_Comm comm,const char path[PETSC_MAX_PATH_LEN],PetscReal ** pynodes,PetscInt * nynodes)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 */
HoneeMeshTransform_PlateMesh(MPI_Comm comm,DM dm,PetscReal growth,PetscInt N,PetscReal refine_height,PetscReal top_angle,PetscReal * node_locs[],PetscInt * num_node_locs)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 **/
HoneeMeshTransformFromOptions(DM dm)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