PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Boundaries.c
Go to the documentation of this file.
1#include "Boundaries.h" // The main header for our project
2#include <string.h> // For strcasecmp
3#include <ctype.h> // For isspace
4
5#undef __FUNCT__
6#define __FUNCT__ "CanRankServiceInletFace"
7/**
8 * @brief Determines if the current MPI rank owns any part of the globally defined inlet face,
9 * making it responsible for placing particles on that portion of the surface.
10 *
11 * The determination is based on the rank's owned nodes (from `DMDALocalInfo`) and
12 * the global node counts, in conjunction with the `user->identifiedInletBCFace`.
13 * A rank can service an inlet face if it owns the cells adjacent to that global boundary
14 * and has a non-zero extent (owns cells) in the tangential dimensions of that face.
15 *
16 * @param user Pointer to the UserCtx structure, containing `identifiedInletBCFace`.
17 * @param info Pointer to the DMDALocalInfo for the current rank's DA (node-based).
18 * @param IM_nodes_global Global number of nodes in the I-direction (e.g., user->IM + 1 if user->IM is cell count).
19 * @param JM_nodes_global Global number of nodes in the J-direction.
20 * @param KM_nodes_global Global number of nodes in the K-direction.
21 * @param[out] can_service_inlet_out Pointer to a PetscBool; set to PETSC_TRUE if the rank
22 * services (part of) the inlet, PETSC_FALSE otherwise.
23 * @return PetscErrorCode 0 on success, non-zero on failure.
24 */
25PetscErrorCode CanRankServiceInletFace(UserCtx *user, const DMDALocalInfo *info,
26 PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
27 PetscBool *can_service_inlet_out)
28{
29 PetscErrorCode ierr;
30 PetscMPIInt rank_for_logging; // For detailed debugging logs
31 PetscFunctionBeginUser;
32 PROFILE_FUNCTION_BEGIN;
33
34 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
35
36 *can_service_inlet_out = PETSC_FALSE; // Default to no service
37
38 if (!user->inletFaceDefined) {
39 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Inlet face not defined in user context. Cannot service.\n", rank_for_logging);
40 PetscFunctionReturn(0);
41 }
42
43 // Get the range of cells owned by this rank in each dimension
44 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
45 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
46 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
47
48 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
49 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
50 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
51
52 // Determine the global index of the last cell (0-indexed) in each direction.
53 // Example: If IM_nodes_global = 11 (nodes 0-10), there are 10 cells (0-9). Last cell index is 9.
54 // Formula: global_nodes - 1 (num cells) - 1 (0-indexed) = global_nodes - 2.
55 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1; // -1 if 0 or 1 node (i.e., 0 cells)
56 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
57 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
58
59 switch (user->identifiedInletBCFace) {
60 case BC_FACE_NEG_X: // Inlet on the global I-minimum face (face of cell C_i=0)
61 // Rank services if its first owned node is global node 0 (info->xs == 0),
62 // and it owns cells in I, J, and K directions.
63 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
64 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
65 *can_service_inlet_out = PETSC_TRUE;
66 }
67 break;
68 case BC_FACE_POS_X: // Inlet on the global I-maximum face (face of cell C_i=last_global_cell_idx_i)
69 // Rank services if it owns the last cell in I-direction,
70 // and has extent in J and K.
71 if (last_global_cell_idx_i >= 0 && /* Check for valid global domain */
72 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i && /* Rank's last cell is the global last cell */
73 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
74 *can_service_inlet_out = PETSC_TRUE;
75 }
76 break;
77 case BC_FACE_NEG_Y:
78 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
79 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
80 *can_service_inlet_out = PETSC_TRUE;
81 }
82 break;
83 case BC_FACE_POS_Y:
84 if (last_global_cell_idx_j >= 0 &&
85 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
86 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
87 *can_service_inlet_out = PETSC_TRUE;
88 }
89 break;
90 case BC_FACE_NEG_Z:
91 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
92 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
93 *can_service_inlet_out = PETSC_TRUE;
94 }
95 break;
96 case BC_FACE_POS_Z:
97 if (last_global_cell_idx_k >= 0 &&
98 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
99 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
100 *can_service_inlet_out = PETSC_TRUE;
101 }
102 break;
103 default:
104 LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d]: Unknown inlet face %s.\n", rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace));
105 break;
106 }
107
108
109 LOG_ALLOW(LOCAL, LOG_TRACE,
110 "[Rank %d] Check Service for Inlet %s:\n"
111 " - Local Domain: starts at cell (%d,%d,%d), has (%d,%d,%d) cells.\n"
112 " - Global Domain: has (%d,%d,%d) nodes, so last cell is (%d,%d,%d).\n",
113 rank_for_logging,
114 BCFaceToString((BCFace)user->identifiedInletBCFace),
115 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
116 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
117 IM_nodes_global, JM_nodes_global, KM_nodes_global,
118 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
119
120 LOG_ALLOW(LOCAL, LOG_DEBUG,"[Rank %d] Inlet Face %s Service Check Result: %s\n",
121 rank_for_logging,
122 BCFaceToString((BCFace)user->identifiedInletBCFace),
123 (*can_service_inlet_out) ? "CAN SERVICE" : "CANNOT SERVICE");
124
125 PROFILE_FUNCTION_END;
126
127 PetscFunctionReturn(0);
128}
129
130#undef __FUNCT__
131#define __FUNCT__ "CanRankServiceFace"
132
133/**
134 * @brief Determines if the current MPI rank owns any part of a specified global face.
135 *
136 * This function is a general utility for parallel boundary operations. It checks if the
137 * local domain of the current MPI rank is adjacent to a specified global boundary face.
138 * A rank "services" a face if it owns the cells adjacent to that global boundary and has
139 * a non-zero extent (i.e., owns at least one cell) in the tangential dimensions of that face.
140 *
141 * @param info Pointer to the DMDALocalInfo for the current rank's DA.
142 * @param IM_nodes_global Global number of nodes in the I-direction (e.g., user->IM + 1 if user->IM is cell count).
143 * @param JM_nodes_global Global number of nodes in the J-direction.
144 * @param KM_nodes_global Global number of nodes in the K-direction.
145 * @param face_id The specific global face (e.g., BC_FACE_NEG_Z) to check.
146 * @param[out] can_service_out Pointer to a PetscBool; set to PETSC_TRUE if the rank
147 * services the face, PETSC_FALSE otherwise.
148 * @return PetscErrorCode 0 on success.
149 */
150PetscErrorCode CanRankServiceFace(const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
151 BCFace face_id, PetscBool *can_service_out)
152{
153 PetscErrorCode ierr;
154 PetscMPIInt rank_for_logging;
155 PetscFunctionBeginUser;
156
157 PROFILE_FUNCTION_BEGIN;
158
159 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
160
161 *can_service_out = PETSC_FALSE; // Default to no service
162
163 // Get the range of cells owned by this rank
164 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
165 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
166 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
167 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
168 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
169 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
170
171 // Determine the global index of the last cell (0-indexed) in each direction.
172 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
173 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
174 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
175
176 switch (face_id) {
177 case BC_FACE_NEG_X:
178 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
179 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
180 *can_service_out = PETSC_TRUE;
181 }
182 break;
183 case BC_FACE_POS_X:
184 if (last_global_cell_idx_i >= 0 &&
185 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i &&
186 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
187 *can_service_out = PETSC_TRUE;
188 }
189 break;
190 case BC_FACE_NEG_Y:
191 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
192 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
193 *can_service_out = PETSC_TRUE;
194 }
195 break;
196 case BC_FACE_POS_Y:
197 if (last_global_cell_idx_j >= 0 &&
198 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
199 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
200 *can_service_out = PETSC_TRUE;
201 }
202 break;
203 case BC_FACE_NEG_Z:
204 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
205 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
206 *can_service_out = PETSC_TRUE;
207 }
208 break;
209 case BC_FACE_POS_Z:
210 if (last_global_cell_idx_k >= 0 &&
211 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
212 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
213 *can_service_out = PETSC_TRUE;
214 }
215 break;
216 default:
217 LOG_ALLOW(LOCAL, LOG_WARNING, "Rank %d: Unknown face enum %d. \n", rank_for_logging, face_id);
218 break;
219 }
220
221 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d check for face %s: Result=%s. \n",
222 rank_for_logging, BCFaceToString((BCFace)face_id), (*can_service_out ? "TRUE" : "FALSE"));
223
224 PROFILE_FUNCTION_END;
225
226 PetscFunctionReturn(0);
227}
228
229#undef __FUNCT__
230#define __FUNCT__ "GetDeterministicFaceGridLocation"
231
232/**
233 * @brief Places particles in a deterministic grid/raster pattern on a specified domain face.
234 *
235 * This function creates a set of equidistant, parallel lines of particles near the four
236 * edges of the face specified by user->identifiedInletBCFace. The number of lines drawn
237 * from each edge is hardcoded within this function (default is 2).
238 *
239 * For example, if grid_layers=2 on face BC_FACE_NEG_X, the function will create particle lines at:
240 * - y ~ 0*dy, y ~ 1*dy (parallel to the Z-axis, starting from the J=0 edge)
241 * - y ~ y_max, y ~ y_max-dy (parallel to the Z-axis, starting from the J=max edge)
242 * - z ~ 0*dz, z ~ 1*dz (parallel to the Y-axis, starting from the K=0 edge)
243 * - z ~ z_max, z ~ z_max-dz (parallel to the Y-axis, starting from the K=max edge)
244 *
245 * The particle's final position is set just inside the target cell face to ensure it is
246 * correctly located. The total number of particles (simCtx->np) is distributed as evenly
247 * as possible among all generated lines.
248 *
249 * The function includes extensive validation to stop with an error if the requested grid
250 * placement is geometrically impossible (e.g., in a 2D domain or if layers would overlap).
251 * It also issues warnings for non-fatal but potentially unintended configurations.
252 *
253 * @param user Pointer to UserCtx, which must contain a valid identifiedInletBCFace.
254 * @param info Pointer to DMDALocalInfo for the current rank's grid layout.
255 * @param xs_gnode_rank, ys_gnode_rank, zs_gnode_rank Local starting node indices (incl. ghosts) for the rank's DA.
256 * @param IM_cells_global, JM_cells_global, KM_cells_global Global cell counts.
257 * @param particle_global_id The unique global ID of the particle being placed (from 0 to np-1).
258 * @param[out] ci_metric_lnode_out Local I-node index of the selected cell's origin.
259 * @param[out] cj_metric_lnode_out Local J-node index of the selected cell's origin.
260 * @param[out] ck_metric_lnode_out Local K-node index of the selected cell's origin.
261 * @param[out] xi_metric_logic_out Logical xi-coordinate [0,1] within the cell.
262 * @param[out] eta_metric_logic_out Logical eta-coordinate [0,1] within the cell.
263 * @param[out] zta_metric_logic_out Logical zta-coordinate [0,1] within the cell.
264 * @param[out] placement_successful_out PETSC_TRUE if the point belongs to this rank, PETSC_FALSE otherwise.
265 * @return PetscErrorCode
266 */
267PetscErrorCode GetDeterministicFaceGridLocation(
268 UserCtx *user, const DMDALocalInfo *info,
269 PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank,
270 PetscInt IM_cells_global, PetscInt JM_cells_global, PetscInt KM_cells_global,
271 PetscInt64 particle_global_id,
272 PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out,
273 PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out,
274 PetscBool *placement_successful_out)
275{
276 SimCtx *simCtx = user->simCtx;
277 PetscReal global_logic_i, global_logic_j, global_logic_k;
278 PetscErrorCode ierr;
279 PetscMPIInt rank_for_logging;
280
281 PetscFunctionBeginUser;
282 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
283
284 *placement_successful_out = PETSC_FALSE; // Default to failure
285
286 // --- Step 1: Configuration and Input Validation ---
287
288 // *** Hardcoded number of grid layers. Change this value to alter the pattern. ***
289 const PetscInt grid_layers = 2;
290
291 LOG_ALLOW(LOCAL, LOG_DEBUG,
292 "[Rank %d] Placing particle %lld on face %s with grid_layers=%d in global domain (%d,%d,%d) cells.\n",
293 rank_for_logging, (long long)particle_global_id, BCFaceToString(user->identifiedInletBCFace), grid_layers,
294 IM_cells_global, JM_cells_global, KM_cells_global);
295
296 const char *face_name = BCFaceToString(user->identifiedInletBCFace);
297
298 // Fatal Error Checks: Ensure the requested grid is geometrically possible.
299 // The total layers from opposite faces (2 * grid_layers) must be less than the domain size.
300 switch (user->identifiedInletBCFace) {
301 case BC_FACE_NEG_X: case BC_FACE_POS_X:
302 if (JM_cells_global <= 1 || KM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (J-cells=%d, K-cells=%d).", face_name, JM_cells_global, KM_cells_global);
303 if (2 * grid_layers >= JM_cells_global || 2 * grid_layers >= KM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing J/K faces would overlap in this domain (J-cells=%d, K-cells=%d).", grid_layers, JM_cells_global, KM_cells_global);
304 break;
305 case BC_FACE_NEG_Y: case BC_FACE_POS_Y:
306 if (IM_cells_global <= 1 || KM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (I-cells=%d, K-cells=%d).", face_name, IM_cells_global, KM_cells_global);
307 if (2 * grid_layers >= IM_cells_global || 2 * grid_layers >= KM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing I/K faces would overlap in this domain (I-cells=%d, K-cells=%d).", grid_layers, IM_cells_global, KM_cells_global);
308 break;
309 case BC_FACE_NEG_Z: case BC_FACE_POS_Z:
310 if (IM_cells_global <= 1 || JM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (I-cells=%d, J-cells=%d).", face_name, IM_cells_global, JM_cells_global);
311 if (2 * grid_layers >= IM_cells_global || 2 * grid_layers >= JM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing I/J faces would overlap in this domain (I-cells=%d, J-cells=%d).", grid_layers, IM_cells_global, JM_cells_global);
312 break;
313 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid identifiedInletBCFace specified: %d", user->identifiedInletBCFace);
314 }
315
316 const PetscInt num_lines_total = 4 * grid_layers;
317 if (simCtx->np < num_lines_total) {
318 LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Total particle count (%lld) is less than the number of grid lines requested (%d). Some lines may be empty.\n", (long long)simCtx->np, num_lines_total);
319 }
320 if (simCtx->np > 0 && simCtx->np % num_lines_total != 0) {
321 LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Total particle count (%lld) is not evenly divisible by the number of grid lines (%d). Distribution will be uneven.\n", (long long)simCtx->np, num_lines_total);
322 }
323
324 // --- Step 2: Map global particle ID to a line and a point on that line ---
325 if (simCtx->np == 0) PetscFunctionReturn(0); // Nothing to do
326
327 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Distributing %lld particles over %d lines on face %s.\n",
328 rank_for_logging, (long long)simCtx->np, num_lines_total, face_name);
329
330 const PetscInt points_per_line = PetscMax(1, simCtx->np / num_lines_total);
331 PetscInt line_index = particle_global_id / points_per_line;
332 PetscInt point_index_on_line = particle_global_id % points_per_line;
333 line_index = PetscMin(line_index, num_lines_total - 1); // Clamp to handle uneven division
334
335 // Decode the line_index into an edge group (0-3) and a layer within that group (0 to grid_layers-1)
336 const PetscInt edge_group = line_index / grid_layers;
337 const PetscInt layer_index = line_index % grid_layers;
338
339 // --- Step 3: Calculate placement coordinates based on the decoded indices ---
340 const PetscReal epsilon = 1.0e-6; // Small offset to keep particles off exact cell boundaries
341 const PetscReal layer_spacing_norm_i = (IM_cells_global > 0) ? 1.0 / (PetscReal)IM_cells_global : 0.0;
342 const PetscReal layer_spacing_norm_j = (JM_cells_global > 0) ? 1.0 / (PetscReal)JM_cells_global : 0.0;
343 const PetscReal layer_spacing_norm_k = (KM_cells_global > 0) ? 1.0 / (PetscReal)KM_cells_global : 0.0;
344
345 PetscReal variable_coord; // The coordinate that varies along a line
346 if (points_per_line <= 1) {
347 variable_coord = 0.5; // Place single point in the middle
348 } else {
349 variable_coord = ((PetscReal)point_index_on_line + 0.5)/ (PetscReal)(points_per_line);
350 }
351 variable_coord = PetscMin(1.0 - epsilon, PetscMax(epsilon, variable_coord)); // Clamp within [eps, 1-eps]
352
353 // Main logic switch to determine the three global logical coordinates
354 switch (user->identifiedInletBCFace) {
355 case BC_FACE_NEG_X:
356 global_logic_i = 0.5 * layer_spacing_norm_i; // Place near the face, in the middle of the first cell
357 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
358 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
359 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
360 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
361 break;
362 case BC_FACE_POS_X:
363 global_logic_i = 1.0 - (0.5 * layer_spacing_norm_i); // Place near the face, in the middle of the last cell
364 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
365 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
366 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
367 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
368 break;
369 case BC_FACE_NEG_Y:
370 global_logic_j = 0.5 * layer_spacing_norm_j;
371 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
372 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
373 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
374 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
375 break;
376 case BC_FACE_POS_Y:
377 global_logic_j = 1.0 - (0.5 * layer_spacing_norm_j);
378 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
379 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
380 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
381 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
382 break;
383 case BC_FACE_NEG_Z:
384 global_logic_k = 0.5 * layer_spacing_norm_k;
385 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
386 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
387 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
388 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
389 break;
390 case BC_FACE_POS_Z:
391 global_logic_k = 1.0 - (0.5 * layer_spacing_norm_k);
392 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
393 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
394 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
395 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
396 break;
397 }
398
399 LOG_ALLOW(LOCAL, LOG_TRACE,
400 "[Rank %d] Particle %lld assigned to line %d (edge group %d, layer %d) with variable_coord=%.4f.\n"
401 " -> Global logical coords: (i,j,k) = (%.6f, %.6f, %.6f)\n",
402 rank_for_logging, (long long)particle_global_id, line_index, edge_group, layer_index, variable_coord,
403 global_logic_i, global_logic_j, global_logic_k);
404
405 // --- Step 4: Convert global logical coordinate to global cell index and intra-cell logicals ---
406 PetscReal global_cell_coord_i = global_logic_i * IM_cells_global;
407 PetscInt I_g = (PetscInt)global_cell_coord_i;
408 *xi_metric_logic_out = global_cell_coord_i - I_g;
409
410 PetscReal global_cell_coord_j = global_logic_j * JM_cells_global;
411 PetscInt J_g = (PetscInt)global_cell_coord_j;
412 *eta_metric_logic_out = global_cell_coord_j - J_g;
413
414 PetscReal global_cell_coord_k = global_logic_k * KM_cells_global;
415 PetscInt K_g = (PetscInt)global_cell_coord_k;
416 *zta_metric_logic_out = global_cell_coord_k - K_g;
417
418 // --- Step 5: Check if this rank owns the target cell and finalize outputs ---
419 if ((I_g >= info->xs && I_g < info->xs + info->xm) &&
420 (J_g >= info->ys && J_g < info->ys + info->ym) &&
421 (K_g >= info->zs && K_g < info->zs + info->zm))
422 {
423 // Convert global cell index to the local node index for this rank's DA patch
424 *ci_metric_lnode_out = (I_g - info->xs) + xs_gnode_rank;
425 *cj_metric_lnode_out = (J_g - info->ys) + ys_gnode_rank;
426 *ck_metric_lnode_out = (K_g - info->zs) + zs_gnode_rank;
427 *placement_successful_out = PETSC_TRUE;
428 }
429
430 LOG_ALLOW(LOCAL, LOG_VERBOSE,
431 "[Rank %d] Particle %lld placement %s. Local cell origin node: (I,J,K) = (%d,%d,%d), intra-cell logicals: (xi,eta,zta)=(%.6f,%.6f,%.6f)\n",
432 rank_for_logging, (long long)particle_global_id,
433 (*placement_successful_out ? "SUCCESSFUL" : "NOT ON THIS RANK"),
434 *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
435 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
436
437 PetscFunctionReturn(0);
438}
439
440#undef __FUNCT__
441#define __FUNCT__ "GetRandomFCellAndLogicOnInletFace"
442
443/**
444 * @brief Assuming the current rank services the inlet face, this function selects a random
445 * cell (owned by this rank on that face) and random logical coordinates within that cell,
446 * suitable for placing a particle on the inlet surface.
447 *
448 * It is the caller's responsibility to ensure CanRankServiceInletFace returned true.
449 *
450 * @param user Pointer to UserCtx.
451 * @param info Pointer to DMDALocalInfo for the current rank (node-based).
452 * @param xs_gnode, ys_gnode, zs_gnode Local starting node indices (incl. ghosts) for the rank's DA.
453 * @param IM_nodes_global, JM_nodes_global, KM_nodes_global Global node counts.
454 * @param rand_logic_i_ptr, rand_logic_j_ptr, rand_logic_k_ptr Pointers to RNGs for logical coords.
455 * @param[out] ci_metric_lnode_out, cj_metric_lnode_out, ck_metric_lnode_out Local node indices of the selected cell's origin (these are local to the rank's DA including ghosts).
456 * @param[out] xi_metric_logic_out, eta_metric_logic_out, zta_metric_logic_out Logical coords [0,1] within the cell.
457 * @return PetscErrorCode
458 */
459PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace(
460 UserCtx *user, const DMDALocalInfo *info,
461 PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, // Local starting node index (with ghosts) of the rank's DA patch
462 PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
463 PetscRandom *rand_logic_i_ptr, PetscRandom *rand_logic_j_ptr, PetscRandom *rand_logic_k_ptr,
464 PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out,
465 PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out)
466{
467 PetscErrorCode ierr = 0;
468 PetscReal r_val_i_sel, r_val_j_sel, r_val_k_sel;
469 PetscInt local_cell_idx_on_face_dim1 = 0; // 0-indexed relative to owned cells on face
470 PetscInt local_cell_idx_on_face_dim2 = 0;
471 PetscMPIInt rank_for_logging;
472
473 PetscFunctionBeginUser;
474
475 PROFILE_FUNCTION_BEGIN;
476
477 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
478
479 // Defaults for cell origin node (local index for the rank's DA patch, including ghosts)
480 *ci_metric_lnode_out = xs_gnode_rank; *cj_metric_lnode_out = ys_gnode_rank; *ck_metric_lnode_out = zs_gnode_rank;
481 // Defaults for logical coordinates
482 *xi_metric_logic_out = 0.5; *eta_metric_logic_out = 0.5; *zta_metric_logic_out = 0.5;
483
484 // Get number of cells this rank owns in each dimension (tangential to the face mainly)
485 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
486 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
487 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
488
489 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
490 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
491 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
492
493 // Index of the last cell (0-indexed) in each global direction
494 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
495 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
496 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
497
498 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d: Inlet face %s. Owned cells(i,j,k):(%d,%d,%d). GlobNodes(I,J,K):(%d,%d,%d). Rank's DA node starts at (%d,%d,%d).\n",
499 rank_for_logging, user->identifiedInletBCFace, num_owned_cells_on_rank_i,num_owned_cells_on_rank_j,num_owned_cells_on_rank_k,
500 IM_nodes_global,JM_nodes_global,KM_nodes_global, xs_gnode_rank,ys_gnode_rank,zs_gnode_rank);
501
502
503 switch (user->identifiedInletBCFace) {
504 case BC_FACE_NEG_X: // Particle on -X face of cell C_0 (origin node N_0)
505 // Cell origin node is the first owned node in I by this rank (global index info->xs).
506 // Its local index within the rank's DA (incl ghosts) is xs_gnode_rank.
507 *ci_metric_lnode_out = xs_gnode_rank;
508 *xi_metric_logic_out = 1.0e-6;
509
510 // Tangential dimensions are J and K. Select an owned cell randomly on this face.
511 // num_owned_cells_on_rank_j/k must be > 0 (checked by CanRankServiceInletFace)
512 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
513 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j); // Index among owned J-cells
514 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
515 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1; // Offset from start of rank's J-nodes
516
517 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
518 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
519 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
520 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
521
522 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
523 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
524 break;
525
526 case BC_FACE_POS_X: // Particle on +X face of cell C_last_I (origin node N_last_I_origin)
527 // Origin node of the last I-cell is global_node_idx = last_global_cell_idx_i.
528 // Its local index in rank's DA: (last_global_cell_idx_i - info->xs) + xs_gnode_rank
529 *ci_metric_lnode_out = xs_gnode_rank + (last_global_cell_idx_i - info->xs);
530 *xi_metric_logic_out = 1.0 - 1.0e-6;
531
532 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
533 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
534 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
535 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1;
536
537 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
538 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
539 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
540 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
541
542 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
543 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
544 break;
545 // ... (Cases for Y and Z faces, following the same pattern) ...
546 case BC_FACE_NEG_Y:
547 *cj_metric_lnode_out = ys_gnode_rank;
548 *eta_metric_logic_out = 1.0e-6;
549 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
550 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
551 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
552 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
553 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
554 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
555 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
556 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
557 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
558 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
559 break;
560 case BC_FACE_POS_Y:
561 *cj_metric_lnode_out = ys_gnode_rank + (last_global_cell_idx_j - info->ys);
562 *eta_metric_logic_out = 1.0 - 1.0e-6;
563 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
564 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
565 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
566 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
567 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
568 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
569 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
570 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
571 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
572 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
573 break;
574 case BC_FACE_NEG_Z: // Your example case
575 *ck_metric_lnode_out = zs_gnode_rank; // Cell origin is the first owned node in K by this rank
576 *zta_metric_logic_out = 1.0e-6; // Place particle slightly inside this cell from its -Z face
577 // Tangential dimensions are I and J
578 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
579 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
580 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
581 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
582
583 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
584 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
585 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
586 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
587
588 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for I
589 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for J
590 break;
591 case BC_FACE_POS_Z:
592 *ck_metric_lnode_out = zs_gnode_rank + (last_global_cell_idx_k - info->zs);
593 *zta_metric_logic_out = 1.0 - 1.0e-6;
594 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
595 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
596 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
597 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
598 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
599 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
600 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
601 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
602 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
603 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
604 break;
605 default:
606 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "GetRandomCellAndLogicOnInletFace: Invalid user->identifiedInletBCFace %d. \n", user->identifiedInletBCFace);
607 }
608
609 PetscReal eps = 1.0e-7;
610 if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
611 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
612 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
613 } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
614 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
615 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
616 } else {
617 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
618 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
619 }
620
621 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Rank %d: Target Cell Node =(%d,%d,%d). (xi,et,zt)=(%.2e,%.2f,%.2f). \n",
622 rank_for_logging, *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
623 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
624
625 PROFILE_FUNCTION_END;
626
627 PetscFunctionReturn(0);
628}
629
630
631
632#undef __FUNCT__
633#define __FUNCT__ "TranslateModernBCsToLegacy"
634
635PetscErrorCode TranslateModernBCsToLegacy(UserCtx *user)
636{
637 PetscFunctionBeginUser;
638 PROFILE_FUNCTION_BEGIN;
639 LOG_ALLOW(GLOBAL,LOG_DEBUG," Translating modern BC config to legacy integer codes...\n");
640
641 for (int i = 0; i < 6; i++) {
642 BCType modern_type = user->boundary_faces[i].mathematical_type;
643 user->bctype[i] = (int)modern_type;
644
645
646 BCFace current_face = (BCFace)i;
647 const char* face_str = BCFaceToString(current_face);
648 const char* bc_type_str = BCTypeToString(modern_type);
649 LOG_ALLOW(GLOBAL,LOG_TRACE," for face %s(%d), legacy type = %d & modern type = %s .\n",face_str,i,user->bctype[i],bc_type_str);
650 }
651 LOG_ALLOW(GLOBAL,LOG_DEBUG," -> Translation complete.\n");
652 PROFILE_FUNCTION_END;
653 PetscFunctionReturn(0);
654}
655
656#undef __FUNCT__
657#define __FUNCT__ "BoundarySystem_ExecuteStep_Legacy"
658
659/**
660 * @brief Acts as a temporary bridge to the legacy boundary condition implementation.
661 *
662 * This function is the key to our integration strategy. It matches the signature
663 * of the modern `BoundarySystem_ExecuteStep` function that SetEulerianFields
664 * expects to call.
665 *
666 * However, instead of containing new handler-based logic, it simply calls the
667 * monolithic legacy `FormBCS` function. This allows the modern orchestrator to
668 * drive the old solver logic without modification.
669 *
670 * @param user The UserCtx for a single block.
671 * @return PetscErrorCode
672 */
673PetscErrorCode BoundarySystem_ExecuteStep_Legacy(UserCtx *user)
674{
675 PetscErrorCode ierr;
676 PetscFunctionBeginUser;
677 PROFILE_FUNCTION_BEGIN;
678 // The sole purpose of this function is to call the old logic for the
679 // specific block context that was passed in.
680 ierr = FormBCS(user); CHKERRQ(ierr);
681
682 // NOTE: The legacy `main` called Block_Interface_U after the FormBCS loop.
683 // We will handle that at a higher level (in our AdvanceSimulation loop)
684 // after all blocks have had their BCs applied for the step.
685 PROFILE_FUNCTION_END;
686 PetscFunctionReturn(0);
687}
688
689#undef __FUNCT__
690#define __FUNCT__ "InflowFlux"
691
692/**
693 * @brief Applies inlet boundary conditions based on the modern BC handling system.
694 *
695 * This function iterates through all 6 domain faces. For each face identified as an
696 * INLET, it applies the velocity profile specified by its assigned handler and
697 * parameters (e.g., 'constant_velocity' with vx,vy,vz or 'parabolic' with u_max).
698 *
699 * It calculates the contravariant flux (Ucont), Cartesian velocity on the face (Ubcs),
700 * and the staggered Cartesian velocity (Ucat). It also computes the total incoming
701 * flux and area across all MPI ranks.
702 *
703 * @param user The main UserCtx struct containing the BC configuration and PETSc objects.
704 * @return PetscErrorCode 0 on success.
705 */
706PetscErrorCode InflowFlux(UserCtx *user)
707{
708 PetscErrorCode ierr;
709 PetscReal lFluxIn = 0.0, lAreaIn = 0.0, AreaSumIn;
710 Vec lCoor;
711 Cmpnts ***ucont, ***ubcs, ***ucat, ***coor, ***csi, ***eta, ***zet, ***cent;
712 PetscReal ***nvert;
713
714 DM da = user->da, fda = user->fda;
715 DMDALocalInfo info = user->info;
716 PetscInt xs = info.xs, xe = info.xs + info.xm;
717 PetscInt ys = info.ys, ye = info.ys + info.ym;
718 PetscInt zs = info.zs, ze = info.zs + info.zm;
719 PetscInt mx = info.mx, my = info.my, mz = info.mz;
720 PetscInt lxs, lxe, lys, lye, lzs, lze;
721
722 // Get context for coordinate transformation if needed by a handler
723 SimCtx *simCtx = user->simCtx;
724 PetscReal CMx_c = simCtx->CMx_c;
725 PetscReal CMy_c = simCtx->CMy_c;
726 PetscReal CMz_c = simCtx->CMz_c;
727
728 PetscFunctionBeginUser;
729
730 PROFILE_FUNCTION_BEGIN;
731
732 lxs = xs; lxe = xe; lys = ys; lye = ye; lzs = zs; lze = ze;
733 if (xs==0) lxs = xs+1; if (ys==0) lys = ys+1; if (zs==0) lzs = zs+1;
734 if (xe==mx) lxe = xe-1; if (ye==my) lye = ye-1; if (ze==mz) lze = ze-1;
735
736 // --- Get PETSc arrays ---
737 DMGetCoordinatesLocal(da,&lCoor);
738 DMDAVecGetArray(fda, lCoor, &coor); // Use local coordinates
739 DMDAVecGetArray(fda, user->Ucont, &ucont);
740 DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs);
741 DMDAVecGetArray(fda, user->Ucat, &ucat);
742 DMDAVecGetArray(da, user->lNvert, &nvert);
743 DMDAVecGetArray(fda, user->Cent, &cent);
744 DMDAVecGetArray(fda, user->lCsi, &csi);
745 DMDAVecGetArray(fda, user->lEta, &eta);
746 DMDAVecGetArray(fda, user->lZet, &zet);
747
748 // --- Main Loop over all 6 faces ---
749 for (PetscInt fn=0; fn<6; fn++) {
750 BoundaryFaceConfig *face_config = &user->boundary_faces[fn];
751
752 if (face_config->mathematical_type != INLET) {
753 continue; // Skip non-inlet faces
754 }
755
756 // This processor only acts if it is on the boundary of the global domain
757 PetscBool is_on_boundary = ( (fn==0 && xs==0) || (fn==1 && xe==mx) ||
758 (fn==2 && ys==0) || (fn==3 && ye==my) ||
759 (fn==4 && zs==0) || (fn==5 && ze==mz) );
760 if (!is_on_boundary) continue;
761
762 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Applying INLET handler for face: %s \n", BCFaceToString(face_config->face_id));
763
764 // --- Loop over the specific face geometry ---
765 switch(face_config->face_id) {
766 case BC_FACE_NEG_X: // -Xi
767 case BC_FACE_POS_X: // +Xi
768 {
769 PetscReal sign = (face_config->face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
770 PetscInt i = (face_config->face_id == BC_FACE_NEG_X) ? xs : mx - 2;
771 for (PetscInt k=lzs; k<lze; k++) {
772 for (PetscInt j=lys; j<lye; j++) {
773 if ( (sign > 0 && nvert[k][j][i+1] > 0.1) || (sign < 0 && nvert[k][j][i] > 0.1) ) continue;
774
775 PetscReal uin_this_point = 0.0;
776 // --- Determine velocity based on the handler for this point ---
777 if (face_config->handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY) {
778 PetscBool found;
779 ierr = GetBCParamReal(face_config->params, "vx", &uin_this_point, &found); CHKERRQ(ierr);
780 } else if(face_config->handler_type== BC_HANDLER_INLET_PARABOLIC){
781 PetscBool found;
782 PetscReal umax,diameter=1.0;
783 ierr = GetBCParamReal(face_config->params,"u_max",&umax,&found); CHKERRQ(ierr);
784 ierr = GetBCParamReal(face_config->params,"diameter",&diameter,&found); CHKERRQ(ierr);
785
786 // Radius is in the YZ-plane for an X-face inlet
787 PetscReal yc = cent[k][j][i + (sign>0)].y - CMy_c;
788 PetscReal zc = cent[k][j][i + (sign>0)].z - CMz_c;
789 PetscReal r = sqrt(yc*yc + zc*zc);
790 PetscReal r_norm = 2.0 * r / diameter;
791 uin_this_point = umax * (1.0 - r_norm * r_norm);
792 if (r_norm > 1.0) uin_this_point = 0.0;
793 }
794 // Add other X-face handlers like 'else if (handler == ...)' here
795
796 // --- Apply the calculated velocity ---
797 PetscReal CellArea = sqrt(csi[k][j][i].z*csi[k][j][i].z + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].x*csi[k][j][i].x);
798 lAreaIn += CellArea;
799 ucont[k][j][i].x = sign * uin_this_point * CellArea;
800 lFluxIn += ucont[k][j][i].x;
801 ubcs[k][j][i + (sign < 0)].x = sign * uin_this_point * csi[k][j][i].x / CellArea;
802 ubcs[k][j][i + (sign < 0)].y = sign * uin_this_point * csi[k][j][i].y / CellArea;
803 ubcs[k][j][i + (sign < 0)].z = sign * uin_this_point * csi[k][j][i].z / CellArea;
804 ucat[k][j][i + (sign > 0)] = ubcs[k][j][i + (sign < 0)];
805 }
806 }
807 } break;
808
809 case BC_FACE_NEG_Y: // -Eta
810 case BC_FACE_POS_Y: // +Eta
811 {
812 PetscReal sign = (face_config->face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
813 PetscInt j = (face_config->face_id == BC_FACE_NEG_Y) ? ys : my - 2;
814 for (PetscInt k=lzs; k<lze; k++) {
815 for (PetscInt i=lxs; i<lxe; i++) {
816 if ( (sign > 0 && nvert[k][j+1][i] > 0.1) || (sign < 0 && nvert[k][j][i] > 0.1) ) continue;
817
818 PetscReal uin_this_point = 0.0;
819 if (face_config->handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY) {
820 PetscBool found;
821 ierr = GetBCParamReal(face_config->params, "vy", &uin_this_point, &found); CHKERRQ(ierr);
822 }else if(face_config->handler_type == BC_HANDLER_INLET_PARABOLIC){
823 PetscBool found;
824 PetscReal umax,diameter=1.0;
825 ierr = GetBCParamReal(face_config->params,"umax",&umax,&found); CHKERRQ(ierr);
826 ierr = GetBCParamReal(face_config->params,"diameter",&diameter,&found); CHKERRQ(ierr);
827
828 // Radius is in the XZ-plane for a Y-face inlet
829 PetscReal xc = cent[k][j + (sign>0)][i].x - CMx_c;
830 PetscReal zc = cent[k][j + (sign>0)][i].z - CMz_c;
831 PetscReal r = sqrt(xc*xc + zc*zc);
832 PetscReal r_norm = 2.0 * r / diameter;
833 uin_this_point = umax * (1.0 - r_norm * r_norm);
834 if (r_norm > 1.0) uin_this_point = 0.0;
835
836 }
837
838 PetscReal CellArea = sqrt(eta[k][j][i].z*eta[k][j][i].z + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].x*eta[k][j][i].x);
839 lAreaIn += CellArea;
840 ucont[k][j][i].y = sign * uin_this_point * CellArea;
841 lFluxIn += ucont[k][j][i].y;
842 ubcs[k][j + (sign < 0)][i].x = sign * uin_this_point * eta[k][j][i].x / CellArea;
843 ubcs[k][j + (sign < 0)][i].y = sign * uin_this_point * eta[k][j][i].y / CellArea;
844 ubcs[k][j + (sign < 0)][i].z = sign * uin_this_point * eta[k][j][i].z / CellArea;
845 ucat[k][j + (sign > 0)][i] = ubcs[k][j + (sign < 0)][i];
846 }
847 }
848 } break;
849
850 case BC_FACE_NEG_Z: // -Zeta
851 case BC_FACE_POS_Z: // +Zeta
852 {
853 PetscReal sign = (face_config->face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
854 PetscInt k = (face_config->face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
855 for (PetscInt j=lys; j<lye; j++) {
856 for (PetscInt i=lxs; i<lxe; i++) {
857 if ( (sign > 0 && nvert[k+1][j][i] > 0.1) || (sign < 0 && nvert[k][j][i] > 0.1) ) continue;
858
859 PetscReal uin_this_point = 0.0;
860 if (face_config->handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY) {
861 PetscBool found;
862 ierr = GetBCParamReal(face_config->params, "vz", &uin_this_point, &found); CHKERRQ(ierr);
863 } else if (face_config->handler_type == BC_HANDLER_INLET_PARABOLIC) {
864 PetscBool found;
865 PetscReal umax, diameter=1.0; // Default diameter for r_norm = 2*r
866 ierr = GetBCParamReal(face_config->params, "u_max", &umax, &found); CHKERRQ(ierr);
867 ierr = GetBCParamReal(face_config->params, "diameter", &diameter, &found); CHKERRQ(ierr); // Optional
868
869 // Radius in the XY-plane for a Z-face inlet.
870 PetscReal xc = cent[k + (sign>0)][j][i].x - CMx_c;
871 PetscReal yc = cent[k + (sign>0)][j][i].y - CMy_c;
872 PetscReal r = sqrt(xc*xc + yc*yc);
873 PetscReal r_norm = 2.0 * r / diameter; // Normalized radius
874 uin_this_point = umax * (1.0 - r_norm * r_norm);
875 if (r_norm > 1.0) uin_this_point = 0.0;
876 }
877
878 PetscReal CellArea = sqrt(zet[k][j][i].z*zet[k][j][i].z + zet[k][j][i].y*zet[k][j][i].y + zet[k][j][i].x*zet[k][j][i].x);
879 lAreaIn += CellArea;
880 ucont[k][j][i].z = sign * uin_this_point * CellArea;
881 lFluxIn += ucont[k][j][i].z;
882 ubcs[k + (sign < 0)][j][i].x = sign * uin_this_point * zet[k][j][i].x / CellArea;
883 ubcs[k + (sign < 0)][j][i].y = sign * uin_this_point * zet[k][j][i].y / CellArea;
884 ubcs[k + (sign < 0)][j][i].z = sign * uin_this_point * zet[k][j][i].z / CellArea;
885 ucat[k + (sign > 0)][j][i] = ubcs[k + (sign < 0)][j][i];
886 }
887 }
888 } break;
889 } // end switch(face_id)
890 } // end for(fn)
891
892 // --- Finalize: Sum flux and area from all processes ---
893 ierr = MPI_Allreduce(&lFluxIn, &simCtx->FluxInSum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
894 ierr = MPI_Allreduce(&lAreaIn, &AreaSumIn, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
895
896 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Inflow: Flux - Area: %le - %le \n", simCtx->FluxInSum, AreaSumIn);
897
898 // --- Restore PETSc arrays ---
899 DMDAVecRestoreArray(fda, lCoor, &coor);
900 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
901 DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs);
902 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
903 DMDAVecRestoreArray(da, user->lNvert, &nvert);
904 DMDAVecRestoreArray(fda, user->Cent, &cent);
905 DMDAVecRestoreArray(fda, user->lCsi, &csi);
906 DMDAVecRestoreArray(fda, user->lEta, &eta);
907 DMDAVecRestoreArray(fda, user->lZet, &zet);
908
909 // --- Update local vectors for subsequent computations ---
910 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
911 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
912
913 PROFILE_FUNCTION_END;
914
915 PetscFunctionReturn(0);
916}
917
918#undef __FUNCT__
919#define __FUNCT__ "OutflowFlux"
920
921/**
922 * @brief Calculates the total outgoing flux through all OUTLET faces for reporting.
923 *
924 * NOTE: In a mixed modern/legacy environment, this function is for DIAGNOSTICS ONLY.
925 * It reads the contravariant velocities and calculates the total flux passing through
926 * faces marked as OUTLET. It does NOT apply any boundary conditions itself, as that
927 * is still the responsibility of the legacy FormBCS function.
928 *
929 * @param user The main UserCtx struct containing BC config and PETSc vectors.
930 * @return PetscErrorCode 0 on success.
931 */
932PetscErrorCode OutflowFlux(UserCtx *user) {
933
934 PetscErrorCode ierr;
935 PetscReal lFluxOut = 0.0;
936 Cmpnts ***ucont;
937
938 DM fda = user->fda;
939 DMDALocalInfo info = user->info;
940 PetscInt xs = info.xs, xe = info.xs + info.xm;
941 PetscInt ys = info.ys, ye = info.ys + info.ym;
942 PetscInt zs = info.zs, ze = info.zs + info.zm;
943 PetscInt mx = info.mx, my = info.my, mz = info.mz;
944
945 PetscFunctionBeginUser;
946
947 PROFILE_FUNCTION_BEGIN;
948
949 ierr = DMDAVecGetArrayRead(fda, user->Ucont, &ucont); CHKERRQ(ierr);
950
951 // --- Loop over all 6 faces to find OUTLETS ---
952 for (PetscInt fn = 0; fn < 6; fn++) {
953 if (user->boundary_faces[fn].mathematical_type != OUTLET) {
954 continue;
955 }
956
957 PetscBool is_on_boundary = ( (fn==0 && xs==0) || (fn==1 && xe==mx) ||
958 (fn==2 && ys==0) || (fn==3 && ye==my) ||
959 (fn==4 && zs==0) || (fn==5 && ze==mz) );
960 if (!is_on_boundary) continue;
961
962 // --- Sum the flux for the appropriate face and component ---
963 switch ((BCFace)fn) {
964 case BC_FACE_NEG_X: case BC_FACE_POS_X: {
965 PetscInt i = (fn == 0) ? xs : mx - 2;
966 for (PetscInt k=info.zs; k<info.zs+info.zm; k++) for (PetscInt j=info.ys; j<info.ys+info.ym; j++) {
967 lFluxOut += ucont[k][j][i].x;
968 }
969 } break;
970
971 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: {
972 PetscInt j = (fn == 2) ? ys : my - 2;
973 for (PetscInt k=info.zs; k<info.zs+info.zm; k++) for (PetscInt i=info.xs; i<info.xs+info.xm; i++) {
974 lFluxOut += ucont[k][j][i].y;
975 }
976 } break;
977
978 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: {
979 PetscInt k = (fn == 4) ? zs : mz - 2;
980 for (PetscInt j=info.ys; j<info.ys+info.ym; j++) for (PetscInt i=info.xs; i<info.xs+info.xm; i++) {
981 lFluxOut += ucont[k][j][i].z;
982 }
983 } break;
984 } // end switch
985 } // end for loop
986
987 ierr = DMDAVecRestoreArrayRead(fda, user->Ucont, &ucont); CHKERRQ(ierr);
988
989 // --- Finalize: Sum and store the global total flux ---
990 ierr = MPI_Allreduce(&lFluxOut, &user->simCtx->FluxOutSum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
991
992 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Reported Global FluxOutSum = %.6f\n", user->simCtx->FluxOutSum);
993
994 PROFILE_FUNCTION_END;
995
996 PetscFunctionReturn(0);
997}
998
999#undef __FUNCT__
1000#define __FUNCT__ "FormBCS"
1001
1002/* Boundary condition defination (array user->bctype[0-5]):
1003 0: interpolation/interface
1004 -1: wallfunction
1005 1: solid wall (not moving)
1006 2: moving solid wall (U=1)
1007 3: slip wall/symmetry
1008 5: Inlet
1009 4: Outlet
1010 6: farfield
1011 7: periodic
1012 8: Characteristic BC
1013 9: Analytical Vortex
1014 10: Oulet Junction
1015 11: Annulus
1016 12: Ogrid
1017 13: Rheology
1018 14: Outlet with Interface
1019 15: No Gradient (Similar to Farfield)
1020*/
1021
1022PetscErrorCode FormBCS(UserCtx *user)
1023{
1024 DM da = user->da, fda = user->fda;
1025 DMDALocalInfo info = user->info;
1026 PetscInt xs = info.xs, xe = info.xs + info.xm;
1027 PetscInt ys = info.ys, ye = info.ys + info.ym;
1028 PetscInt zs = info.zs, ze = info.zs + info.zm;
1029 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1030 PetscInt lxs, lxe, lys, lye, lzs, lze;
1031 PetscInt i, j, k;
1032
1033 PetscReal ***nvert,***lnvert; //local working array
1034
1035 PetscReal ***p,***lp;
1036 Cmpnts ***ucont, ***ubcs, ***ucat,***lucat, ***csi, ***eta, ***zet;
1037 Cmpnts ***cent,***centx,***centy,***centz,***coor;
1038 PetscScalar FluxIn, FluxOut,ratio;
1039 PetscScalar lArea, AreaSum;
1040
1041 PetscScalar FarFluxIn=0., FarFluxOut=0., FarFluxInSum, FarFluxOutSum;
1042 PetscScalar FarAreaIn=0., FarAreaOut=0., FarAreaInSum, FarAreaOutSum;
1043 PetscScalar FluxDiff, VelDiffIn, VelDiffOut;
1044 Cmpnts V_frame;
1045
1046 PetscReal Un, nx,ny,nz,A;
1047
1048 SimCtx *simCtx = user->simCtx;
1049
1050 lxs = xs; lxe = xe;
1051 lys = ys; lye = ye;
1052 lzs = zs; lze = ze;
1053
1054 PetscInt gxs, gxe, gys, gye, gzs, gze;
1055
1056 gxs = info.gxs; gxe = gxs + info.gxm;
1057 gys = info.gys; gye = gys + info.gym;
1058 gzs = info.gzs; gze = gzs + info.gzm;
1059
1060 if (xs==0) lxs = xs+1;
1061 if (ys==0) lys = ys+1;
1062 if (zs==0) lzs = zs+1;
1063
1064 if (xe==mx ) lxe = xe-1;
1065 if (ye==my ) lye = ye-1;
1066 if (ze==mz ) lze = ze-1;
1067
1068 PetscFunctionBeginUser;
1069
1070 PROFILE_FUNCTION_BEGIN;
1071
1072 DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs);
1073
1074 DMDAVecGetArray(fda, user->lCsi, &csi);
1075 DMDAVecGetArray(fda, user->lEta, &eta);
1076 DMDAVecGetArray(fda, user->lZet, &zet);
1077
1078 PetscInt ttemp;
1079 for (ttemp=0; ttemp<3; ttemp++) {
1080 DMDAVecGetArray(da, user->Nvert, &nvert);
1081 DMDAVecGetArray(fda, user->lUcat, &ucat);
1082 DMDAVecGetArray(fda, user->Ucont, &ucont);
1083/* ================================================================================== */
1084/* FAR-FIELD BC */
1085/* ================================================================================== */
1086
1087
1088 // reset FAR FLUXES
1089 FarFluxIn = 0.; FarFluxOut=0.;
1090 FarAreaIn = 0.; FarAreaOut=0.;
1091
1092 PetscReal lFlux_abs=0.0,FluxSum_abs=0.0,ratio=0.0;
1093
1094 V_frame.x=0.;
1095 V_frame.y=0.;
1096 V_frame.z=0.;
1097
1098
1099
1100 if (user->bctype[0]==6) {
1101 if (xs == 0) {
1102 i= xs;
1103 for (k=lzs; k<lze; k++) {
1104 for (j=lys; j<lye; j++) {
1105 ubcs[k][j][i].x = ucat[k][j][i+1].x;
1106 ubcs[k][j][i].y = ucat[k][j][i+1].y;
1107 ubcs[k][j][i].z = ucat[k][j][i+1].z;
1108 ucont[k][j][i].x = ubcs[k][j][i].x * csi[k][j][i].x;
1109 FarFluxIn += ucont[k][j][i].x;
1110 lFlux_abs += fabs(ucont[k][j][i].x);
1111 FarAreaIn += csi[k][j][i].x;
1112 }
1113 }
1114 }
1115 }
1116
1117 if (user->bctype[1]==6) {
1118 if (xe==mx) {
1119 i= xe-1;
1120 for (k=lzs; k<lze; k++) {
1121 for (j=lys; j<lye; j++) {
1122 ubcs[k][j][i].x = ucat[k][j][i-1].x;
1123 ubcs[k][j][i].y = ucat[k][j][i-1].y;
1124 ubcs[k][j][i].z = ucat[k][j][i-1].z;
1125 ucont[k][j][i-1].x = ubcs[k][j][i].x * csi[k][j][i-1].x;
1126 FarFluxOut += ucont[k][j][i-1].x;
1127 lFlux_abs += fabs(ucont[k][j][i-1].x);
1128 FarAreaOut += csi[k][j][i-1].x;
1129 }
1130 }
1131 }
1132 }
1133
1134 if (user->bctype[2]==6) {
1135 if (ys==0) {
1136 j= ys;
1137 for (k=lzs; k<lze; k++) {
1138 for (i=lxs; i<lxe; i++) {
1139 ubcs[k][j][i].x = ucat[k][j+1][i].x;
1140 ubcs[k][j][i].y = ucat[k][j+1][i].y;
1141 ubcs[k][j][i].z = ucat[k][j+1][i].z;
1142 ucont[k][j][i].y = ubcs[k][j][i].y * eta[k][j][i].y;
1143 FarFluxIn += ucont[k][j][i].y;
1144 lFlux_abs += fabs(ucont[k][j][i].y);
1145 FarAreaIn += eta[k][j][i].y;
1146 }
1147 }
1148 }
1149 }
1150
1151 if (user->bctype[3]==6) {
1152 if (ye==my) {
1153 j=ye-1;
1154 for (k=lzs; k<lze; k++) {
1155 for (i=lxs; i<lxe; i++) {
1156 ubcs[k][j][i].x = ucat[k][j-1][i].x;
1157 ubcs[k][j][i].y = ucat[k][j-1][i].y;
1158 ubcs[k][j][i].z = ucat[k][j-1][i].z;
1159 ucont[k][j-1][i].y = ubcs[k][j][i].y * eta[k][j-1][i].y;
1160 FarFluxOut += ucont[k][j-1][i].y;
1161 lFlux_abs += fabs(ucont[k][j-1][i].y);
1162 FarAreaOut += eta[k][j-1][i].y;
1163 }
1164 }
1165 }
1166 }
1167
1168 if (user->bctype[4]==6 || user->bctype[4]==10) {
1169 if (zs==0) {
1170 k = 0;
1171 for (j=lys; j<lye; j++) {
1172 for (i=lxs; i<lxe; i++) {
1173 ubcs[k][j][i].x = ucat[k+1][j][i].x;
1174 ubcs[k][j][i].y = ucat[k+1][j][i].y;
1175 ubcs[k][j][i].z = ucat[k+1][j][i].z;
1176 ucont[k][j][i].z = ubcs[k][j][i].z * zet[k][j][i].z;
1177 FarFluxIn += ucont[k][j][i].z;
1178 lFlux_abs += fabs(ucont[k][j][i].z);
1179 FarAreaIn += zet[k][j][i].z;
1180 }
1181 }
1182 }
1183 }
1184
1185 if (user->bctype[5]==6 || user->bctype[5]==10) {
1186 if (ze==mz) {
1187 k = ze-1;
1188 for (j=lys; j<lye; j++) {
1189 for (i=lxs; i<lxe; i++) {
1190 ubcs[k][j][i].x = ucat[k-1][j][i].x;
1191 ubcs[k][j][i].y = ucat[k-1][j][i].y;
1192 ubcs[k][j][i].z = ucat[k-1][j][i].z;
1193 ucont[k-1][j][i].z = ubcs[k][j][i].z * zet[k-1][j][i].z;
1194 FarFluxOut += ucont[k-1][j][i].z;
1195 lFlux_abs += fabs(ucont[k-1][j][i].z);
1196 FarAreaOut += zet[k-1][j][i].z;
1197 }
1198 }
1199 }
1200 }
1201
1202 MPI_Allreduce(&FarFluxIn,&FarFluxInSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1203 MPI_Allreduce(&FarFluxOut,&FarFluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1204 MPI_Allreduce(&lFlux_abs,&FluxSum_abs,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1205 MPI_Allreduce(&FarAreaIn,&FarAreaInSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1206 MPI_Allreduce(&FarAreaOut,&FarAreaOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1207
1208 if (user->bctype[5]==6 || user->bctype[3]==6 || user->bctype[1]==6) {
1209
1210 ratio=(FarFluxInSum - FarFluxOutSum)/FluxSum_abs;
1211 if (fabs(FluxSum_abs) <1.e-10) ratio = 0.;
1212
1213 LOG_ALLOW(GLOBAL,LOG_DEBUG, "/FluxSum_abs %le ratio %le \n", FluxSum_abs,ratio);
1214
1215 FluxDiff = 0.5*(FarFluxInSum - FarFluxOutSum) ;
1216 VelDiffIn = FluxDiff / FarAreaInSum ;
1217
1218 if (fabs(FarAreaInSum) <1.e-6) VelDiffIn = 0.;
1219
1220 VelDiffOut = FluxDiff / FarAreaOutSum ;
1221
1222 if (fabs(FarAreaOutSum) <1.e-6) VelDiffOut = 0.;
1223
1224 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Far Flux Diff %d %le %le %le %le %le %le %le\n", simCtx->step, FarFluxInSum, FarFluxOutSum, FluxDiff, FarAreaInSum, FarAreaOutSum, VelDiffIn, VelDiffOut);
1225
1226 }
1227
1228 if (user->bctype[5]==10) {
1229 FluxDiff = simCtx->FluxInSum -( FarFluxOutSum -FarFluxInSum) ;
1230 VelDiffIn = 1/3.*FluxDiff / (FarAreaInSum);// + FarsimCtx->AreaOutSum);
1231 if (fabs(FarAreaInSum) <1.e-6) VelDiffIn = 0.;
1232
1233 VelDiffOut = 2./3.* FluxDiff / (FarAreaOutSum) ;//(FarAreaInSum + FarsimCtx->AreaOutSum) ;
1234
1235 if (fabs(FarAreaOutSum) <1.e-6) VelDiffOut = 0.;
1236
1237 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Far Flux Diff %d %le %le %le %le %le %le %le\n", simCtx->step, FarFluxInSum, FarFluxOutSum, FluxDiff, FarAreaInSum, FarAreaOutSum, VelDiffIn, VelDiffOut);
1238
1239 }
1240
1241
1242 // scale global mass conservation
1243
1244 if (user->bctype[5]==6 || user->bctype[5]==10) {
1245 if (ze==mz) {
1246 k = ze-1;
1247 for (j=lys; j<lye; j++) {
1248 for (i=lxs; i<lxe; i++) {
1249 ucont[k-1][j][i].z += ratio*fabs(ucont[k-1][j][i].z);
1250 ubcs[k][j][i].z = ucont[k-1][j][i].z/zet[k-1][j][i].z;
1251 // ubcs[k][j][i].z = ucat[k-1][j][i].z + VelDiffOut ;//+ V_frame.z;
1252 // ucont[k-1][j][i].z = ubcs[k][j][i].z * zet[k-1][j][i].z;
1253 }
1254 }
1255 }
1256 }
1257
1258 if (user->bctype[3]==6) {
1259 if (ye==my) {
1260 j=ye-1;
1261 for (k=lzs; k<lze; k++) {
1262 for (i=lxs; i<lxe; i++) {
1263
1264 ucont[k][j-1][i].y +=ratio*fabs(ucont[k][j-1][i].y);
1265 ubcs[k][j][i].y = ucont[k][j-1][i].y /eta[k][j-1][i].y;
1266 // ubcs[k][j][i].y = ucat[k][j-1][i].y + VelDiffOut;// + V_frame.y;
1267 // ucont[k][j-1][i].y = ubcs[k][j][i].y * eta[k][j-1][i].y;
1268
1269 }
1270 }
1271 }
1272 }
1273
1274 if (user->bctype[1]==6) {
1275 if (xe==mx) {
1276 i= xe-1;
1277 for (k=lzs; k<lze; k++) {
1278 for (j=lys; j<lye; j++) {
1279 ucont[k][j][i-1].x +=ratio*fabs(ucont[k][j][i-1].x);
1280 ubcs[k][j][i].x = ucont[k][j][i-1].x / csi[k][j][i-1].x ;
1281 // ubcs[k][j][i].x = ucat[k][j][i-1].x + VelDiffOut;// + V_frame.x;
1282 // ucont[k][j][i-1].x = ubcs[k][j][i].x * csi[k][j][i-1].x;
1283 }
1284 }
1285 }
1286 }
1287
1288
1289 if (user->bctype[0]==6) {
1290 if (xs == 0) {
1291 i= xs;
1292 for (k=lzs; k<lze; k++) {
1293 for (j=lys; j<lye; j++) {
1294 ucont[k][j][i].x -=ratio*fabs(ucont[k][j][i].x);
1295 ubcs[k][j][i].x = ucont[k][j][i].x / csi[k][j][i].x;
1296 // ubcs[k][j][i].x = ucat[k][j][i+1].x - VelDiffIn;// + V_frame.x;
1297 // ucont[k][j][i].x = ubcs[k][j][i].x * csi[k][j][i].x;
1298 }
1299 }
1300 }
1301 }
1302
1303
1304 if (user->bctype[2]==6) {
1305 if (ys==0) {
1306 j= ys;
1307 for (k=lzs; k<lze; k++) {
1308 for (i=lxs; i<lxe; i++) {
1309 ucont[k][j][i].y -=ratio*fabs(ucont[k][j][i].y);
1310 ubcs[k][j][i].y = ucont[k][j][i].y / eta[k][j][i].y;
1311 // ubcs[k][j][i].y = ucat[k][j+1][i].y - VelDiffIn;// + V_frame.y;
1312 // ucont[k][j][i].y = ubcs[k][j][i].y * eta[k][j][i].y;
1313 }
1314 }
1315 }
1316 }
1317
1318
1319 if (user->bctype[4]==6 || user->bctype[5]==10) {
1320 if (zs==0) {
1321 k = 0;
1322 for (j=lys; j<lye; j++) {
1323 for (i=lxs; i<lxe; i++) {
1324 ucont[k][j][i].z -=ratio*fabs(ucont[k][j][i].z);
1325 ubcs[k][j][i].z =ucont[k][j][i].z / zet[k][j][i].z;
1326 // ubcs[k][j][i].z = ucat[k+1][j][i].z - VelDiffIn;// + V_frame.z;
1327 // ucont[k][j][i].z = ubcs[k][j][i].z * zet[k][j][i].z;
1328
1329 }
1330 }
1331 }
1332 }
1333
1334//// Amir wall Ogrid
1335
1336if (user->bctype[2]==1 || user->bctype[2]==-1) {
1337 if (ys==0) {
1338 j= ys;
1339 for (k=lzs; k<lze; k++) {
1340 for (i=lxs; i<lxe; i++) {
1341 A=sqrt(eta[k][j][i].z*eta[k][j][i].z +
1342 eta[k][j][i].y*eta[k][j][i].y +
1343 eta[k][j][i].x*eta[k][j][i].x);
1344 nx=eta[k][j][i].x/A;
1345 ny=eta[k][j][i].y/A;
1346 nz=eta[k][j][i].z/A;
1347 Un=ucat[k][j+1][i].x*nx+ucat[k][j+1][i].y*ny+ucat[k][j+1][i].z*nz;
1348 ubcs[k][j][i].x = 0.0;
1349 ubcs[k][j][i].y = 0.0;
1350 ubcs[k][j][i].z = 0.0;
1351 ucont[k][j][i].y = 0.;
1352 }
1353 }
1354 }
1355 }
1356
1357/* ================================================================================== */
1358/* SOLID WALL BC (NO-SLIP / NO-PENETRATION) */
1359/* ================================================================================== */
1360
1361// NOTE: This block is added to explicitly handle bctype=1 (solid wall) for all faces.
1362// It ensures both no-slip (ubcs=0) and no-penetration (ucont_normal=0).
1363// ubcs is handled by the implicit-zero assumption, but ucont must be set explicitly.
1364
1365// -X Face (i=0)
1366if (user->bctype[0]==1 || user->bctype[0]==-1) {
1367 if (xs==0) {
1368 i= xs;
1369 for (k=lzs; k<lze; k++) {
1370 for (j=lys; j<lye; j++) {
1371 ubcs[k][j][i].x = 0.0;
1372 ubcs[k][j][i].y = 0.0;
1373 ubcs[k][j][i].z = 0.0;
1374 ucont[k][j][i].x = 0.0; // Enforce no-penetration
1375 }
1376 }
1377 }
1378}
1379
1380// +X Face (i=mx-1)
1381if (user->bctype[1]==1 || user->bctype[1]==-1) {
1382 if (xe==mx) {
1383 i= xe-1;
1384 for (k=lzs; k<lze; k++) {
1385 for (j=lys; j<lye; j++) {
1386 ubcs[k][j][i].x = 0.0;
1387 ubcs[k][j][i].y = 0.0;
1388 ubcs[k][j][i].z = 0.0;
1389 // The relevant ucont is at the face, index i-1
1390 ucont[k][j][i-1].x = 0.0; // Enforce no-penetration
1391 }
1392 }
1393 }
1394}
1395
1396// -Y Face (j=0)
1397if (user->bctype[2]==1 || user->bctype[2]==-1) {
1398 if (ys==0) {
1399 j= ys;
1400 for (k=lzs; k<lze; k++) {
1401 for (i=lxs; i<lxe; i++) {
1402 ubcs[k][j][i].x = 0.0;
1403 ubcs[k][j][i].y = 0.0;
1404 ubcs[k][j][i].z = 0.0;
1405 ucont[k][j][i].y = 0.0; // Enforce no-penetration
1406 }
1407 }
1408 }
1409}
1410
1411// +Y Face (j=my-1)
1412if (user->bctype[3]==1 || user->bctype[3]==-1) {
1413 if (ye==my) {
1414 j= ye-1;
1415 for (k=lzs; k<lze; k++) {
1416 for (i=lxs; i<lxe; i++) {
1417 ubcs[k][j][i].x = 0.0;
1418 ubcs[k][j][i].y = 0.0;
1419 ubcs[k][j][i].z = 0.0;
1420 // The relevant ucont is at the face, index j-1
1421 ucont[k][j-1][i].y = 0.0; // Enforce no-penetration
1422 }
1423 }
1424 }
1425}
1426
1427/* Original "Amir wall Ogrid" block can now be removed or commented out
1428 as its functionality is included above.
1429if (user->bctype[2]==1 || user->bctype[2]==-1) { ... }
1430*/
1431
1432/* ================================================================================== */
1433/* SYMMETRY BC */
1434/* ================================================================================== */
1435
1436 if (user->bctype[0]==3) {
1437 if (xs==0) {
1438 i= xs;
1439
1440 for (k=lzs; k<lze; k++) {
1441 for (j=lys; j<lye; j++) {
1442 A=sqrt(csi[k][j][i].z*csi[k][j][i].z +
1443 csi[k][j][i].y*csi[k][j][i].y +
1444 csi[k][j][i].x*csi[k][j][i].x);
1445 nx=csi[k][j][i].x/A;
1446 ny=csi[k][j][i].y/A;
1447 nz=csi[k][j][i].z/A;
1448 Un=ucat[k][j][i+1].x*nx+ucat[k][j][i+1].y*ny+ucat[k][j][i+1].z*nz;
1449 ubcs[k][j][i].x = ucat[k][j][i+1].x-Un*nx;//-V_frame.x;
1450 ubcs[k][j][i].y = ucat[k][j][i+1].y-Un*ny;
1451 ubcs[k][j][i].z = ucat[k][j][i+1].z-Un*nz;
1452 ucont[k][j][i].x = 0.;
1453 }
1454 }
1455 }
1456 }
1457
1458 if (user->bctype[1]==3) {
1459 if (xe==mx) {
1460 i= xe-1;
1461
1462 for (k=lzs; k<lze; k++) {
1463 for (j=lys; j<lye; j++) {
1464 A=sqrt(csi[k][j][i-1].z*csi[k][j][i-1].z +
1465 csi[k][j][i-1].y*csi[k][j][i-1].y +
1466 csi[k][j][i-1].x*csi[k][j][i-1].x);
1467 nx=csi[k][j][i-1].x/A;
1468 ny=csi[k][j][i-1].y/A;
1469 nz=csi[k][j][i-1].z/A;
1470 Un=ucat[k][j][i-1].x*nx+ucat[k][j][i-1].y*ny+ucat[k][j][i-1].z*nz;
1471 ubcs[k][j][i].x = ucat[k][j][i-1].x-Un*nx;
1472 ubcs[k][j][i].y = ucat[k][j][i-1].y-Un*ny;
1473 ubcs[k][j][i].z = ucat[k][j][i-1].z-Un*nz;
1474 ucont[k][j][i-1].x = 0.;
1475 }
1476 }
1477 }
1478 }
1479
1480 if (user->bctype[2]==3) {
1481 if (ys==0) {
1482 j= ys;
1483
1484 for (k=lzs; k<lze; k++) {
1485 for (i=lxs; i<lxe; i++) {
1486 A=sqrt(eta[k][j][i].z*eta[k][j][i].z +
1487 eta[k][j][i].y*eta[k][j][i].y +
1488 eta[k][j][i].x*eta[k][j][i].x);
1489 nx=eta[k][j][i].x/A;
1490 ny=eta[k][j][i].y/A;
1491 nz=eta[k][j][i].z/A;
1492 Un=ucat[k][j+1][i].x*nx+ucat[k][j+1][i].y*ny+ucat[k][j+1][i].z*nz;
1493 ubcs[k][j][i].x = ucat[k][j+1][i].x-Un*nx;
1494 ubcs[k][j][i].y = ucat[k][j+1][i].y-Un*ny;
1495 ubcs[k][j][i].z = ucat[k][j+1][i].z-Un*nz;
1496 ucont[k][j][i].y = 0.;
1497 }
1498 }
1499 }
1500 }
1501
1502 if (user->bctype[3]==3) {
1503 if (ye==my) {
1504 j=ye-1;
1505
1506 for (k=lzs; k<lze; k++) {
1507 for (i=lxs; i<lxe; i++) {
1508 A=sqrt(eta[k][j-1][i].z*eta[k][j-1][i].z +
1509 eta[k][j-1][i].y*eta[k][j-1][i].y +
1510 eta[k][j-1][i].x*eta[k][j-1][i].x);
1511 nx=eta[k][j-1][i].x/A;
1512 ny=eta[k][j-1][i].y/A;
1513 nz=eta[k][j-1][i].z/A;
1514 Un=ucat[k][j-1][i].x*nx+ucat[k][j-1][i].y*ny+ucat[k][j-1][i].z*nz;
1515 ubcs[k][j][i].x = ucat[k][j-1][i].x-Un*nx;
1516 ubcs[k][j][i].y = ucat[k][j-1][i].y-Un*ny;
1517 ubcs[k][j][i].z = ucat[k][j-1][i].z-Un*nz;
1518 ucont[k][j-1][i].y = 0.;
1519 }
1520 }
1521 }
1522 }
1523
1524
1525 if (user->bctype[4]==3) {
1526 if (zs==0) {
1527 k = 0;
1528 for (j=lys; j<lye; j++) {
1529 for (i=lxs; i<lxe; i++) {
1530 A=sqrt(zet[k][j][i].z*zet[k][j][i].z +
1531 zet[k][j][i].y*zet[k][j][i].y +
1532 zet[k][j][i].x*zet[k][j][i].x);
1533 nx=zet[k][j][i].x/A;
1534 ny=zet[k][j][i].y/A;
1535 nz=zet[k][j][i].z/A;
1536 Un=ucat[k+1][j][i].x*nx+ucat[k+1][j][i].y*ny+ucat[k+1][j][i].z*nz;
1537 ubcs[k][j][i].x = ucat[k+1][j][i].x-Un*nx;
1538 ubcs[k][j][i].y = ucat[k+1][j][i].y-Un*ny;
1539 ubcs[k][j][i].z = ucat[k+1][j][i].z-Un*nz;
1540 ucont[k][j][i].z = 0.;
1541 }
1542 }
1543 }
1544 }
1545
1546 if (user->bctype[5]==3) {
1547 if (ze==mz) {
1548 k = ze-1;
1549 for (j=lys; j<lye; j++) {
1550 for (i=lxs; i<lxe; i++) {
1551 A=sqrt(zet[k-1][j][i].z*zet[k-1][j][i].z +
1552 zet[k-1][j][i].y*zet[k-1][j][i].y +
1553 zet[k-1][j][i].x*zet[k-1][j][i].x);
1554 nx=zet[k-1][j][i].x/A;
1555 ny=zet[k-1][j][i].y/A;
1556 nz=zet[k-1][j][i].z/A;
1557 Un=ucat[k-1][j][i].x*nx+ucat[k-1][j][i].y*ny+ucat[k-1][j][i].z*nz;
1558 ubcs[k][j][i].x = ucat[k-1][j][i].x-Un*nx;
1559 ubcs[k][j][i].y = ucat[k-1][j][i].y-Un*ny;
1560 ubcs[k][j][i].z = ucat[k-1][j][i].z-Un*nz;
1561 ucont[k-1][j][i].z = 0.;
1562 }
1563 }
1564 }
1565 }
1566
1567/* ================================================================================== */
1568/* CHARACTERISTIC OUTLET BC :8 */
1569/* ================================================================================== */
1570
1571 if (user->bctype[5]==8) {
1572 if (ze == mz) {
1573 k = ze-2;
1574 FluxOut = 0;
1575 for (j=lys; j<lye; j++) {
1576 for (i=lxs; i<lxe; i++) {
1577 FluxOut += ucont[k][j][i].z;
1578 }
1579 }
1580 }
1581 else {
1582 FluxOut = 0.;
1583 }
1584
1585 FluxIn = simCtx->FluxInSum + FarFluxInSum;
1586
1587 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1588 // PetscGlobalSum(PETSC_COMM_WORLD,&FluxOut, &FluxOutSum);
1589
1590 //ratio = FluxInSum / FluxOutSum;
1591 ratio = FluxIn / simCtx->FluxOutSum;
1592 if (fabs(simCtx->FluxOutSum) < 1.e-6) ratio = 1.;
1593 //if (fabs(FluxInSum) <1.e-6) ratio = 0.;
1594 if (fabs(FluxIn) <1.e-6) ratio = 0.;
1595 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Char Ratio %d %le %le %le %le %d %d\n", simCtx->step, ratio, FluxIn, simCtx->FluxOutSum, FarFluxInSum,zs, ze);
1596
1597 if (ze==mz) {
1598 k = ze-1;
1599 for (j=lys; j<lye; j++) {
1600 for (i=lxs; i<lxe; i++) {
1601 ubcs[k][j][i].x = ucat[k-1][j][i].x;
1602 ubcs[k][j][i].y = ucat[k-1][j][i].y;
1603 if (simCtx->step==0 || simCtx->step==1)
1604 if (simCtx->inletprofile<0)
1605 ubcs[k][j][i].z = -1.;
1606 else if (user->bctype[4]==6)
1607 ubcs[k][j][i].z = 0.;
1608 else
1609 ubcs[k][j][i].z = 1.;//ubcs[0][j][i].z;//-1.;//1.;
1610
1611 else
1612 ucont[k-1][j][i].z = ucont[k-1][j][i].z*ratio;
1613
1614 ubcs[k][j][i].z = ucont[k-1][j][i].z / zet[k-1][j][i].z;
1615 }
1616 }
1617 }
1618 }
1619
1620
1621/* ================================================================================== */
1622/* OUTLETBC :4 */
1623/* ================================================================================== */
1624
1625
1626 if (user->bctype[5]==OUTLET || user->bctype[5]==14 || user->bctype[5]==20) {
1627 lArea=0.;
1628 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"+Zeta Outlet \n");
1629 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"FluxOutSum before FormBCS applied = %.6f \n",simCtx->FluxOutSum);
1630 if (ze == mz) {
1631 // k = ze-3;
1632 k=ze-1;
1633 FluxOut = 0;
1634 for (j=lys; j<lye; j++) {
1635 for (i=lxs; i<lxe; i++) {
1636
1637 if ((nvert[k-1][j][i])<0.1) {
1638 FluxOut += (ucat[k-1][j][i].x * (zet[k-1][j][i].x) +
1639 ucat[k-1][j][i].y * (zet[k-1][j][i].y) +
1640 ucat[k-1][j][i].z * (zet[k-1][j][i].z));
1641
1642 lArea += sqrt( (zet[k-1][j][i].x) * (zet[k-1][j][i].x) +
1643 (zet[k-1][j][i].y) * (zet[k-1][j][i].y) +
1644 (zet[k-1][j][i].z) * (zet[k-1][j][i].z));
1645
1646 }
1647 }
1648 }
1649 }
1650 else {
1651 FluxOut = 0.;
1652 }
1653
1654 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1655 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1656
1657 user->simCtx->AreaOutSum = AreaSum;
1658
1659 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"AreaOutSum = %.6f | FluxOutSum = %.6f \n",AreaSum,simCtx->FluxOutSum);
1660
1661 if (simCtx->block_number>1 && user->bctype[5]==14) {
1662 simCtx->FluxOutSum += user->FluxIntfcSum;
1663 // AreaSum += user->AreaIntfcSum;
1664 }
1665
1666 FluxIn = simCtx->FluxInSum + FarFluxInSum + user->FluxIntpSum;
1667 if (user->bctype[5]==20)
1668 ratio = (FluxIn / simCtx->FluxOutSum);
1669 else
1670 ratio = (FluxIn - simCtx->FluxOutSum) / AreaSum;
1671
1672 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Ratio for momentum correction = %.6f \n",ratio);
1673
1674 /* user->FluxOutSum += ratio*user->simCtx->AreaOutSum; */
1675 simCtx->FluxOutSum =0.0;
1676 FluxOut=0.0;
1677 if (ze==mz) {
1678 k = ze-1;
1679 for (j=lys; j<lye; j++) {
1680 for (i=lxs; i<lxe; i++) {
1681 if ((nvert[k-1][j][i])<0.1) {
1682
1683 ubcs[k][j][i].x = ucat[k-1][j][i].x;//+ratio;
1684 ubcs[k][j][i].y = ucat[k-1][j][i].y;
1685 ubcs[k][j][i].z = ucat[k-1][j][i].z;// + ratio;//*n_z;
1686
1687 // ucont[k-1][j][i].z = ubcs[k][j][i].z * zet[k-1][j][i].z;
1688 if (user->bctype[5]==20)
1689 ucont[k-1][j][i].z = (ubcs[k][j][i].x * (zet[k-1][j][i].x ) +
1690 ubcs[k][j][i].y * (zet[k-1][j][i].y ) +
1691 ubcs[k][j][i].z * (zet[k-1][j][i].z ))*ratio;
1692
1693 else{
1694 ucont[k-1][j][i].z = (ubcs[k][j][i].x * (zet[k-1][j][i].x ) +
1695 ubcs[k][j][i].y * (zet[k-1][j][i].y ) +
1696 ubcs[k][j][i].z * (zet[k-1][j][i].z ))
1697 + ratio * sqrt( (zet[k-1][j][i].x) * (zet[k-1][j][i].x) +
1698 (zet[k-1][j][i].y) * (zet[k-1][j][i].y) +
1699 (zet[k-1][j][i].z) * (zet[k-1][j][i].z));
1700
1701 FluxOut += ucont[k-1][j][i].z;
1702 }
1703 }//if
1704 }
1705 }
1706 }
1707
1708 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1709 LOG_ALLOW(GLOBAL,LOG_TRACE, "Timestep = %d | FluxInSum = %.6f | FlucOutSum = %.6f | FluxIntfcSum = %.6f | FluxIntpSum = %.6f \n", simCtx->step, simCtx->FluxInSum, simCtx->FluxOutSum, user->FluxIntfcSum,user->FluxIntpSum);
1710
1711 } else if (user->bctype[5]==2) {
1712 /* Designed for driven cavity problem (top(k=kmax) wall moving)
1713 u_x = 1 at k==kmax */
1714 if (ze==mz) {
1715 k = ze-1;
1716 for (j=lys; j<lye; j++) {
1717 for (i=lxs; i<lxe; i++) {
1718 ubcs[k][j][i].x = 0.;// - ucat[k-1][j][i].x;
1719 ubcs[k][j][i].y = 1;//sin(2*3.14*simCtx->step*simCtx->dt);//1.;//- ucat[k-1][j][i].y;
1720 ubcs[k][j][i].z = 0.;//- ucat[k-1][j][i].z;
1721 }
1722 }
1723 }
1724 }
1725
1726
1727 /* OUTLET at k==0 */
1728 if (user->bctype[4]==OUTLET) {
1729 lArea=0.;
1730 if (zs == 0) {
1731 k = zs;
1732 // k= zs + 1;
1733 FluxOut = 0;
1734 for (j=lys; j<lye; j++) {
1735 for (i=lxs; i<lxe; i++) {
1736
1737
1738 FluxOut += ucat[k+1][j][i].z * zet[k][j][i].z ;
1739
1740 lArea += zet[k][j][i].z;
1741
1742
1743
1744 }
1745 }
1746 }
1747 else {
1748 FluxOut = 0.;
1749 }
1750
1751 FluxIn = simCtx->FluxInSum + FarFluxInSum;
1752
1753 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1754 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1755
1756
1757 ratio = (simCtx->FluxInSum - simCtx->FluxOutSum) / AreaSum;
1758 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Ratio b %d %le %le %le %le %d %d\n", simCtx->step,ratio, simCtx->FluxInSum, simCtx->FluxOutSum, AreaSum,zs, ze);
1759
1760 if (zs==0) {
1761 k = 0;
1762 for (j=lys; j<lye; j++) {
1763 for (i=lxs; i<lxe; i++) {
1764 ubcs[k][j][i].x = ucat[k+1][j][i].x;
1765 ubcs[k][j][i].y = ucat[k+1][j][i].y;
1766 ubcs[k][j][i].z = ucat[k+1][j][i].z;
1767 ucont[k][j][i].z = (ubcs[k][j][i].z+ratio) * zet[k][j][i].z;
1768 }
1769 }
1770 }
1771 }
1772
1773
1774
1775/* ================================================================================== */
1776/* Ogrid :77 */
1777/* ================================================================================== */
1778 /*
1779 if (user->bctype[3]=77 && Ogrid)
1780 {Cmpnts ***coor;
1781 lArea=0.;
1782 FluxOut=0.0;
1783 // k = ze-3;
1784
1785 Vec Coor; DMGetCoordinatesLocal(da, &Coor);
1786 DMDAVecGetArray(fda,Coor,&coor);
1787 if (ye==my) {
1788 j=my-2;
1789 for (k=lzs; k<lze; k++){
1790 for (i=lxs; i<lxe; i++) {
1791
1792 FluxOut += ucont[k][j][i].y;
1793 lArea += sqrt( (eta[k-1][j][i].x) * (eta[k-1][j][i].x) +
1794 (eta[k-1][j][i].y) * (eta[k-1][j][i].y) +
1795 (eta[k-1][j][i].z) * (eta[k-1][j][i].z));
1796
1797 }
1798 }
1799 }
1800
1801 else {
1802 FluxOut = 0.;
1803 }
1804
1805 MPI_Allreduce(&FluxOut,&FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1806 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1807
1808 user->FluxOutSum = FluxOutSum;
1809 simCtx->AreaOutSum = AreaSum;
1810 export PL="$HOME/CSBL/Codes/fdf_project/pic_les"
1811 ratio=2*(FluxIn-FluxOutSum)/AreaSum;
1812 ratio=0.0;
1813 if (ye==my){
1814 j=my-2;
1815 for (k=lzs; k<lze; k++){
1816 for (i=lxs; i<lxe; i++) {
1817 if ((nvert[k-1][j][i])<0.1 && coor[k][j][i].z >= 800.) {
1818 ucont[k][j][i].y *= (1+ratio);
1819
1820 }
1821 }
1822 }
1823 }
1824
1825 if (ye==my){
1826 j=my-2;
1827 for (k=lzs; k<lze; k++){
1828 for (i=lxs; i<lxe; i++) {
1829
1830 ubcs[k][j+1][i].z=ucat[k][j][i].z;
1831 ubcs[k][j+1][i].y=ucat[k][j][i].y;
1832 ubcs[k][j+1][i].x=ucat[k][j][i].x;
1833 }
1834 }
1835 }
1836
1837 // LOG_ALLOW(GLOBAL,LOG_DEBUG, " ratio %le ",ratio);
1838
1839 }
1840
1841 */
1842
1843
1844/* ================================================================================== */
1845/* Channelz */
1846/* ================================================================================== */
1847 // Amir channel flow correction
1848 if (user->bctype[4]==7 && simCtx->channelz==1) {
1849 Vec Coor; DMGetCoordinatesLocal(da, &Coor);
1850 DMDAVecGetArray(fda,Coor,&coor);
1851 Cmpnts ***uch;
1852 DMDAVecGetArray(fda, user->Bcs.Uch, &uch);
1853
1854 lArea=0.0;
1855 // if (zs==0) {
1856 // k=0;
1857 FluxIn=0.0;
1858
1859 double Fluxbcs=0.0, Fluxbcssum,ratiobcs;
1860
1861 if (zs==0) {
1862 k=0;
1863 Fluxbcs=0.0;
1864 for (j=lys; j<lye; j++) {
1865 for (i=lxs; i<lxe; i++) {
1866 if (nvert[k][j][i]<0.1){
1867 Fluxbcs += ucont[k][j][i].z;
1868
1869 lArea += sqrt((zet[k][j][i].x) * (zet[k][j][i].x) +
1870 (zet[k][j][i].y) * (zet[k][j][i].y) +
1871 (zet[k][j][i].z) * (zet[k][j][i].z));
1872 }
1873 }
1874 }
1875 }
1876
1877
1878 //int kk=(simCtx->step % (mz-2))+2;
1879
1880 for (k=zs;k<lze;k++){
1881 for (j=lys; j<lye; j++) {
1882 for (i=lxs; i<lxe; i++) {
1883 if (nvert[k][j][i]<0.1){
1884 FluxIn += ucont[k][j][i].z /((mz)-1);
1885 }
1886 }
1887 }
1888 }
1889 // else {
1890 // FluxIn=0.0;
1891// }
1892
1893 MPI_Allreduce(&FluxIn,&simCtx->Fluxsum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1894 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1895 MPI_Allreduce(&Fluxbcs,&Fluxbcssum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1896
1897 if (simCtx->step==simCtx->StartStep && simCtx->StartStep > 0 && simCtx->ccc==0) {
1898 simCtx->ccc=1;
1899 simCtx->FluxInSum=Fluxbcssum;
1900// simCtx->FluxInSum=6.3908;
1901 LOG_ALLOW(LOCAL,LOG_DEBUG, " FluxInSum %le .\n ",simCtx->FluxInSum);
1902 }
1903
1904 simCtx->FluxInSum=6.35066;
1905 ratio=(simCtx->FluxInSum-simCtx->Fluxsum)/AreaSum;
1906 simCtx->ratio=ratio;
1907 ratiobcs=(simCtx->FluxInSum-Fluxbcssum)/AreaSum;
1908
1909 if (zs==0) {
1910 k=0;
1911 for (j=lys; j<lye; j++) {
1912 for (i=lxs; i<lxe; i++) {
1913 if (nvert[k+1][j][i]<0.1){
1914 if (simCtx->fish) {
1915 ucont[k][j][i].z+=ratiobcs * /* (1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt( (zet[k+1][j][i].x) * (zet[k+1][j][i].x) +
1916 (zet[k+1][j][i].y) * (zet[k+1][j][i].y) +
1917 (zet[k+1][j][i].z) * (zet[k+1][j][i].z));
1918 }
1919
1920 uch[k][j][i].z=ratiobcs * /* (1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt( (zet[k+1][j][i].x) * (zet[k+1][j][i].x) +
1921 (zet[k+1][j][i].y) * (zet[k+1][j][i].y) +
1922 (zet[k+1][j][i].z) * (zet[k+1][j][i].z));
1923 }
1924 }
1925 }
1926 }
1927
1928 if (ze==mz) {
1929 k=mz-1;
1930 for (j=lys; j<lye; j++) {
1931 for (i=lxs; i<lxe; i++) {
1932 if (nvert[k-1][j][i]<0.1){
1933 if (simCtx->fish){
1934 ucont[k-1][j][i].z+=ratiobcs * /*(1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt((zet[k-1][j][i].x) * (zet[k-1][j][i].x) +
1935 (zet[k-1][j][i].y) * (zet[k-1][j][i].y) +
1936 (zet[k-1][j][i].z) * (zet[k-1][j][i].z));
1937 }
1938 uch[k][j][i].z=ratiobcs * /*(1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt( (zet[k+1][j][i].x) * (zet[k+1][j][i].x) +
1939 (zet[k+1][j][i].y) * (zet[k+1][j][i].y) +
1940 (zet[k+1][j][i].z) * (zet[k+1][j][i].z));
1941 }
1942 }
1943 }
1944 }
1945 DMDAVecRestoreArray(fda, user->Bcs.Uch, &uch);
1946 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Ratio %le %.15le %.15le %.15le \n", ratio, ratiobcs, simCtx->FluxInSum,AreaSum);
1947 DMDAVecRestoreArray(fda,Coor,&coor);
1948
1949 ////.................////
1950
1951
1952
1953
1954 //just for check
1955/* if (zs==0) { */
1956/* k=10; */
1957/* FluxIn=0.0; */
1958/* for (j=lys; j<lye; j++) { */
1959/* for (i=lxs; i<lxe; i++) { */
1960/* if (nvert[k+1][j][i]<0.1){ */
1961/* FluxIn += ucont[k][j][i].z; */
1962/* lArea += sqrt((zet[k+1][j][i].x) * (zet[k+1][j][i].x) + */
1963/* (zet[k+1][j][i].y) * (zet[k+1][j][i].y) + */
1964/* (zet[k+1][j][i].z) * (zet[k+1][j][i].z)); */
1965/* } */
1966/* } */
1967/* } */
1968/* } */
1969/* else { */
1970/* FluxIn=0.0; */
1971/* } */
1972/* // LOG_ALLOW(PETSC_COMM_SELF, " Fluxsum %le ",FluxIn); */
1973
1974/* MPI_Allreduce(&FluxIn,&Fluxsum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); */
1975/* MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); */
1976
1977
1978
1979 }//channel
1980
1981
1982
1983 /* ==================================================================================== */
1984 /* Cylinder O-grid */
1985 /* ==================================================================================== */
1986 if (user->bctype[3]==12) {
1987 /* Designed to test O-grid for flow over a cylinder at jmax velocity is 1 (horizontal)
1988 u_x = 1 at k==kmax */
1989 if (ye==my) {
1990 j = ye-1;
1991 for (k=lzs; k<lze; k++) {
1992 for (i=lxs; i<lxe; i++) {
1993 ubcs[k][j][i].x = 0.;
1994 ubcs[k][j][i].y = 0.;
1995 ubcs[k][j][i].z = 1.;
1996 }
1997 }
1998 }
1999 }
2000 /* ==================================================================================== */
2001 /* Annulus */
2002 /* ==================================================================================== */
2003 /* designed to test periodic boundary condition for O-grid j=0 rotates */
2004 DMDAVecGetArray(fda, user->Cent, &cent);
2005 if (user->bctype[2]==11) {
2006 if (ys==0){
2007 j=0;
2008 for (k=lzs; k<lze; k++) {
2009 for (i=lxs; i<lxe; i++) {
2010
2011 /* ubcs[k][j][i].x=0.0; */
2012
2013/* ubcs[k][j][i].y = -cent[k][j+1][i].z/sqrt(cent[k][j+1][i].z*cent[k][j+1][i].z+cent[k][j+1][i].y*cent[k][j+1][i].y); */
2014
2015/* ubcs[k][j][i].z =cent[k][j+1][i].y/sqrt(cent[k][j+1][i].z*cent[k][j+1][i].z+cent[k][j+1][i].y*cent[k][j+1][i].y); */
2016 ubcs[k][j][i].x = cent[k][j+1][i].y/sqrt(cent[k][j+1][i].x*cent[k][j+1][i].x+cent[k][j+1][i].y*cent[k][j+1][i].y);;
2017 ubcs[k][j][i].y =-cent[k][j+1][i].x/sqrt(cent[k][j+1][i].x*cent[k][j+1][i].x+cent[k][j+1][i].y*cent[k][j+1][i].y);
2018 ubcs[k][j][i].z =0.0;
2019 // if(k==1) LOG_ALLOW(PETSC_COMM_SELF, "@ i= %d j=%d k=%d ubcs.y is %le\n",i,j,k,ubcs[k][j][i].y);
2020 }
2021 }
2022 }
2023 }
2024 /* ==================================================================================== */
2025 /* Rheology */
2026 /* ==================================================================================== */
2027
2028 if(simCtx->rheology && (user->bctype[2]==13 || user->bctype[3]==13 || user->bctype[4]==13 || user->bctype[5]==13)){
2029 LOG_ALLOW(GLOBAL,LOG_DEBUG, "moving plate velocity for rheology setup is %le \n",simCtx->U_bc);
2030 }
2031
2032 if (user->bctype[2]==13){
2033 if (ys==0){
2034 j=0;
2035 for (k=lzs; k<lze; k++) {
2036 for (i=lxs; i<lxe; i++) {
2037 ubcs[k][j][i].x = 0.;
2038 // ubcs[k][j][i].x = -simCtx->U_bc;
2039 ubcs[k][j][i].y = 0.;
2040 ubcs[k][j][i].z = -simCtx->U_bc;
2041 //ubcs[k][j][i].z =0.0;
2042 }
2043 }
2044 }
2045 }
2046 if (user->bctype[3]==13){
2047 if (ye==my){
2048 j=ye-1;
2049 for (k=lzs; k<lze; k++) {
2050 for (i=lxs; i<lxe; i++) {
2051 ubcs[k][j][i].x = 0.;
2052 // ubcs[k][j][i].x = simCtx->U_bc;
2053 ubcs[k][j][i].y = 0.;
2054 ubcs[k][j][i].z = simCtx->U_bc;
2055 //ubcs[k][j][i].z =0.0;
2056 }
2057 }
2058 }
2059 }
2060
2061 if (user->bctype[4]==13){
2062 if (zs==0){
2063 k=0;
2064 for (j=lys; j<lye; j++) {
2065 for (i=lxs; i<lxe; i++) {
2066 ubcs[k][j][i].x =-simCtx->U_bc;
2067 ubcs[k][j][i].y = 0.;
2068 ubcs[k][j][i].z = 0.;
2069 }
2070 }
2071 }
2072 }
2073 if (user->bctype[5]==13){
2074 if (ze==mz){
2075 k=ze-1;
2076 for (j=lys; j<lye; j++) {
2077 for (i=lxs; i<lxe; i++) {
2078 ubcs[k][j][i].x = simCtx->U_bc;
2079 ubcs[k][j][i].y = 0.;
2080 ubcs[k][j][i].z = 0.;
2081 }
2082 }
2083 }
2084 }
2085 DMDAVecRestoreArray(fda, user->Cent, &cent);
2086 DMDAVecRestoreArray(fda, user->lUcat, &ucat);
2087 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
2088 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2089 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2090 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2091
2092
2093
2094 Contra2Cart(user); // it also does global to local for Ucat
2095
2096/* ================================================================================== */
2097/* WALL FUNCTION */
2098/* ================================================================================== */
2099
2100 if (simCtx->wallfunction && user->bctype[2]==-1) {
2101 PetscReal ***ustar, ***aj;
2102 //Mohsen Dec 2015
2103 Vec Aj = user->lAj;
2104 DMDAVecGetArray(fda, user->Ucat, &ucat);
2105 DMDAVecGetArray(fda, user->Ucont, &ucont);
2106 DMDAVecGetArray(da, Aj, &aj);
2107// DMDAVecGetArray(da, user->Nvert, &nvert);
2108 DMDAVecGetArray(da, user->lUstar, &ustar);
2109
2110 // wall function for boundary
2111 for (k=lzs; k<lze; k++)
2112 for (j=lys; j<lye; j++)
2113 for (i=lxs; i<lxe; i++) {
2114
2115 if( nvert[k][j][i]<1.1 && user->bctype[2]==-1 && j==1 )
2116 {
2117 double area = sqrt( eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z );
2118 double sb, sc;
2119 double ni[3], nj[3], nk[3];
2120 Cmpnts Uc, Ua;
2121
2122 Ua.x = Ua.y = Ua.z = 0;
2123 sb = 0.5/aj[k][j][i]/area;
2124
2125
2126 sc = 2*sb + 0.5/aj[k][j+1][i]/area;
2127 Uc = ucat[k][j+1][i];
2128
2129
2130 //Calculate_normal(csi[k][j][i], eta[k][j][i], zet[k][j][i], ni, nj, nk);
2131 //if(j==my-2) nj[0]*=-1, nj[1]*=-1, nj[2]*=-1;
2132
2133 double AA=sqrt(eta[k][j][i].z*eta[k][j][i].z +
2134 eta[k][j][i].y*eta[k][j][i].y +
2135 eta[k][j][i].x*eta[k][j][i].x);
2136 nj[0]=eta[k][j][i].x/AA;
2137 nj[1]=eta[k][j][i].y/AA;
2138 nj[2]=eta[k][j][i].z/AA;
2139 noslip (user, sc, sb, Ua, Uc, &ucat[k][j][i], nj[0], nj[1], nj[2]);
2140 wall_function_loglaw(user, 1.e-16, sc, sb, Ua, Uc, &ucat[k][j][i], &ustar[k][j][i], nj[0], nj[1], nj[2]);
2141
2142 // nvert[k][j][i]=1.; /* set nvert to 1 to exclude from rhs */
2143 // if (k==1)
2144 // LOG_ALLOW(GLOBAL,LOG_DEBUG, " %d %le %le %le %le %le %le %le %le %le\n",i, sb,aj[k][j][i],AA, nj[0], nj[1], nj[2],ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z);
2145
2146 }
2147 }
2148 if (ys==0) {
2149 j= ys;
2150
2151 for (k=lzs; k<lze; k++) {
2152 for (i=lxs; i<lxe; i++) {
2153 ubcs[k][j][i].x = 0.0;
2154 ubcs[k][j][i].y = 0.0;
2155 ubcs[k][j][i].z = 0.0;
2156 ucont[k][j][i].y = 0.;
2157 }
2158 }
2159 }
2160
2161 DMDAVecRestoreArray(da, Aj, &aj);
2162 DMDAVecRestoreArray(da, user->lUstar, &ustar);
2163 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2164 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
2165 // DMDAVecRestoreArray(da, user->Nvert, &nvert);
2166 // DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2167 // DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2168
2169 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2170 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2171// DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat);
2172// DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat);
2173
2174
2175 }
2176
2177/* ================================================================================== */
2178
2179 DMDAVecGetArray(fda, user->Ucat, &ucat);
2180
2181/* ================================================================================== */
2182
2183 // boundary conditions on ghost nodes
2184 if (xs==0 && user->bctype[0]!=7) {
2185 i = xs;
2186 for (k=lzs; k<lze; k++) {
2187 for (j=lys; j<lye; j++) {
2188 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i+1].x;
2189 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i+1].y;
2190 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i+1].z;
2191 }
2192 }
2193 }
2194
2195 if (xe==mx && user->bctype[0]!=7) {
2196 i = xe-1;
2197 for (k=lzs; k<lze; k++) {
2198 for (j=lys; j<lye; j++) {
2199 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i-1].x;
2200 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i-1].y;
2201 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i-1].z;
2202 }
2203 }
2204 }
2205
2206
2207 if (ys==0 && user->bctype[2]!=7) {
2208 j = ys;
2209 for (k=lzs; k<lze; k++) {
2210 for (i=lxs; i<lxe; i++) {
2211 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j+1][i].x;
2212 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j+1][i].y;
2213 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j+1][i].z;
2214 }
2215 }
2216 }
2217
2218 if (ye==my && user->bctype[2]!=7) {
2219 j = ye-1;
2220 for (k=lzs; k<lze; k++) {
2221 for (i=lxs; i<lxe; i++) {
2222 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j-1][i].x;
2223 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j-1][i].y;
2224 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j-1][i].z;
2225 }
2226 }
2227 }
2228
2229 if (zs==0 && user->bctype[4]!=7) {
2230 k = zs;
2231 for (j=lys; j<lye; j++) {
2232 for (i=lxs; i<lxe; i++) {
2233 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k+1][j][i].x;
2234 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k+1][j][i].y;
2235 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k+1][j][i].z;
2236 }
2237 }
2238 }
2239
2240 if (ze==mz && user->bctype[4]!=7) {
2241 k = ze-1;
2242 for (j=lys; j<lye; j++) {
2243 for (i=lxs; i<lxe; i++) {
2244 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k-1][j][i].x;
2245 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k-1][j][i].y;
2246 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k-1][j][i].z;
2247 }
2248 }
2249 }
2250
2251 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2252 /* ================================================================================== */
2253 /* Periodic BC *///Mohsen
2254/* ================================================================================== */
2255 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7){
2256 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2257 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2258 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2259 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2260 /* /\* DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert); *\/ */
2261 /* /\* DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert); *\/ */
2262 //Mohsen Dec 2015
2263 DMDAVecGetArray(da, user->lP, &lp);
2264 DMDAVecGetArray(da, user->lNvert, &lnvert);
2265 DMDAVecGetArray(fda, user->lUcat, &lucat);
2266 DMDAVecGetArray(da, user->P, &p);
2267 DMDAVecGetArray(da, user->Nvert, &nvert);
2268 DMDAVecGetArray(fda, user->Ucat, &ucat);
2269
2270 if (user->bctype[0]==7 || user->bctype[1]==7){
2271 if (xs==0){
2272 i=xs;
2273 for (k=zs; k<ze; k++) {
2274 for (j=ys; j<ye; j++) {
2275 if(j>0 && k>0 && j<user->JM && k<user->KM){
2276 ucat[k][j][i]=lucat[k][j][i-2];
2277 p[k][j][i]=lp[k][j][i-2];
2278 nvert[k][j][i]=lnvert[k][j][i-2];
2279 }
2280 }
2281 }
2282 }
2283 }
2284 if (user->bctype[2]==7 || user->bctype[3]==7){
2285 if (ys==0){
2286 j=ys;
2287 for (k=zs; k<ze; k++) {
2288 for (i=xs; i<xe; i++) {
2289 if(i>0 && k>0 && i<user->IM && k<user->KM){
2290 ucat[k][j][i]=lucat[k][j-2][i];
2291 p[k][j][i]=lp[k][j-2][i];
2292 nvert[k][j][i]=lnvert[k][j-2][i];
2293 }
2294 }
2295 }
2296 }
2297 }
2298 if (user->bctype[4]==7 || user->bctype[5]==7){
2299 if (zs==0){
2300 k=zs;
2301 for (j=ys; j<ye; j++) {
2302 for (i=xs; i<xe; i++) {
2303 if(i>0 && j>0 && i<user->IM && j<user->JM){
2304 ucat[k][j][i]=lucat[k-2][j][i];
2305 nvert[k][j][i]=lnvert[k-2][j][i];
2306 //amir
2307 p[k][j][i]=lp[k-2][j][i];
2308 }
2309 }
2310 }
2311 }
2312 }
2313 if (user->bctype[0]==7 || user->bctype[1]==7){
2314 if (xe==mx){
2315 i=mx-1;
2316 for (k=zs; k<ze; k++) {
2317 for (j=ys; j<ye; j++) {
2318 if(j>0 && k>0 && j<user->JM && k<user->KM){
2319 ucat[k][j][i]=lucat[k][j][i+2];
2320 p[k][j][i]=lp[k][j][i+2];
2321 nvert[k][j][i]=lnvert[k][j][i+2];
2322 }
2323 }
2324 }
2325 }
2326 }
2327 if (user->bctype[2]==7 || user->bctype[3]==7){
2328 if (ye==my){
2329 j=my-1;
2330 for (k=zs; k<ze; k++) {
2331 for (i=xs; i<xe; i++) {
2332 if(i>0 && k>0 && i<user->IM && k<user->KM){
2333 ucat[k][j][i]=lucat[k][j+2][i];
2334 p[k][j][i]=lp[k][j+2][i];
2335 nvert[k][j][i]=lnvert[k][j+2][i];
2336 }
2337 }
2338 }
2339 }
2340 }
2341
2342 if (user->bctype[4]==7 || user->bctype[5]==7){
2343 if (ze==mz){
2344 k=mz-1;
2345 for (j=ys; j<ye; j++) {
2346 for (i=xs; i<xe; i++) {
2347 if(i>0 && j>0 && i<user->IM && j<user->JM){
2348 ucat[k][j][i]=lucat[k+2][j][i];
2349 nvert[k][j][i]=lnvert[k+2][j][i];
2350 //amir
2351 p[k][j][i]=lp[k+2][j][i];
2352 }
2353 }
2354 }
2355 }
2356 }
2357
2358
2359
2360
2361
2362 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2363 DMDAVecRestoreArray(da, user->lP, &lp);
2364 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2365 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2366 DMDAVecRestoreArray(da, user->P, &p);
2367 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2368
2369 /* /\* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); *\/ */
2370 /* /\* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); *\/ */
2371
2372 /* /\* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); *\/ */
2373 /* /\* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); *\/ */
2374
2375 /* /\* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); *\/ */
2376 /* /\* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); *\/ */
2377
2378 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2379 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2380
2381 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2382 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2383
2384 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2385 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2386}
2387 // 0 velocity on the corner point
2388
2389 DMDAVecGetArray(fda, user->Ucat, &ucat);
2390 DMDAVecGetArray(da, user->P, &p);
2391
2392 if (zs==0) {
2393 k=0;
2394 if (xs==0) {
2395 i=0;
2396 for (j=ys; j<ye; j++) {
2397 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i+1].x);
2398 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i+1].y);
2399 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i+1].z);
2400 p[k][j][i]= 0.5*(p[k+1][j][i]+p[k][j][i+1]);
2401 }
2402 }
2403 if (xe == mx) {
2404 i=mx-1;
2405 for (j=ys; j<ye; j++) {
2406 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i-1].x);
2407 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i-1].y);
2408 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i-1].z);
2409 p[k][j][i] = 0.5*(p[k+1][j][i]+p[k][j][i-1]);
2410 }
2411 }
2412
2413 if (ys==0) {
2414 j=0;
2415 for (i=xs; i<xe; i++) {
2416 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j+1][i].x);
2417 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j+1][i].y);
2418 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j+1][i].z);
2419 p[k][j][i] = 0.5*(p[k+1][j][i]+p[k][j+1][i]);
2420 }
2421 }
2422
2423 if (ye==my) {
2424 j=my-1;
2425 for (i=xs; i<xe; i++) {
2426 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j-1][i].x);
2427 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j-1][i].y);
2428 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j-1][i].z);
2429 p[k][j][i] = 0.5*(p[k+1][j][i]+p[k][j-1][i]);
2430 }
2431 }
2432 }
2433
2434 if (ze==mz) {
2435 k=mz-1;
2436 if (xs==0) {
2437 i=0;
2438 for (j=ys; j<ye; j++) {
2439 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i+1].x);
2440 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i+1].y);
2441 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i+1].z);
2442 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j][i+1]);
2443 }
2444 }
2445 if (xe == mx) {
2446 i=mx-1;
2447 for (j=ys; j<ye; j++) {
2448 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i-1].x);
2449 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i-1].y);
2450 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i-1].z);
2451 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j][i-1]);
2452 }
2453 }
2454
2455 if (ys==0) {
2456 j=0;
2457 for (i=xs; i<xe; i++) {
2458 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j+1][i].x);
2459 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j+1][i].y);
2460 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j+1][i].z);
2461 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j+1][i]);
2462 }
2463 }
2464
2465 if (ye==my) {
2466 j=my-1;
2467 for (i=xs; i<xe; i++) {
2468 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j-1][i].x);
2469 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j-1][i].y);
2470 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j-1][i].z);
2471 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j-1][i]);
2472 }
2473 }
2474 }
2475
2476 if (ys==0) {
2477 j=0;
2478 if (xs==0) {
2479 i=0;
2480 for (k=zs; k<ze; k++) {
2481 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i+1].x);
2482 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i+1].y);
2483 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i+1].z);
2484 p[k][j][i]= 0.5*(p[k][j+1][i]+p[k][j][i+1]);
2485 }
2486 }
2487
2488 if (xe==mx) {
2489 i=mx-1;
2490 for (k=zs; k<ze; k++) {
2491 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i-1].x);
2492 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i-1].y);
2493 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i-1].z);
2494 p[k][j][i] = 0.5*(p[k][j+1][i]+p[k][j][i-1]);
2495 }
2496 }
2497 }
2498
2499 if (ye==my) {
2500 j=my-1;
2501 if (xs==0) {
2502 i=0;
2503 for (k=zs; k<ze; k++) {
2504 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i+1].x);
2505 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i+1].y);
2506 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i+1].z);
2507 p[k][j][i] = 0.5*(p[k][j-1][i]+p[k][j][i+1]);
2508 }
2509 }
2510
2511 if (xe==mx) {
2512 i=mx-1;
2513 for (k=zs; k<ze; k++) {
2514 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i-1].x);
2515 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i-1].y);
2516 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i-1].z);
2517 p[k][j][i] = 0.5*(p[k][j-1][i]+p[k][j][i-1]);
2518 }
2519 }
2520 }
2521
2522 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2523 DMDAVecRestoreArray(da, user->P, &p);
2524
2525 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2526 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2527
2528 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2529 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2530
2531 //Mohsen Nov 2012
2532 //Velocity and Presurre at corners for Periodic BC's
2533
2534 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7){
2535 //i-direction
2536
2537 DMDAVecGetArray(fda, user->lUcat, &lucat);
2538 DMDAVecGetArray(da, user->lP, &lp);
2539 DMDAVecGetArray(da, user->lNvert, &lnvert);
2540 DMDAVecGetArray(fda, user->Ucat, &ucat);
2541 DMDAVecGetArray(da, user->P, &p);
2542 DMDAVecGetArray(da, user->Nvert, &nvert);
2543
2544 if (user->bctype[0]==7){
2545 if (xs==0){
2546 i=xs;
2547 for (k=zs; k<ze; k++) {
2548 for (j=ys; j<ye; j++) {
2549 ucat[k][j][i]=lucat[k][j][i-2];
2550 p[k][j][i]=lp[k][j][i-2];
2551 nvert[k][j][i]=lnvert[k][j][i-2];
2552 }
2553 }
2554 }
2555 }
2556 if (user->bctype[1]==7){
2557 if (xe==mx){
2558 i=xe-1;
2559 for (k=zs; k<ze; k++) {
2560 for (j=ys; j<ye; j++) {
2561 ucat[k][j][i].x=lucat[k][j][i+2].x;
2562 p[k][j][i]=lp[k][j][i+2];
2563 nvert[k][j][i]=lnvert[k][j][i+2];
2564 }
2565 }
2566 }
2567 }
2568 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2569 DMDAVecRestoreArray(da, user->lP, &lp);
2570 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2571 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2572 DMDAVecRestoreArray(da, user->P, &p);
2573 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2574
2575 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2576/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2577
2578/* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); */
2579/* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); */
2580
2581/* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2582/* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2583
2584 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2585 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2586
2587 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2588 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2589
2590 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2591 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2592
2593 //j-direction
2594 DMDAVecGetArray(fda, user->lUcat, &lucat);
2595 DMDAVecGetArray(da, user->lP, &lp);
2596 DMDAVecGetArray(da, user->lNvert, &lnvert);
2597 DMDAVecGetArray(fda, user->Ucat, &ucat);
2598 DMDAVecGetArray(da, user->P, &p);
2599 DMDAVecGetArray(da, user->Nvert, &nvert);
2600
2601 if (user->bctype[2]==7){
2602 if (ys==0){
2603 j=ys;
2604 for (k=zs; k<ze; k++) {
2605 for (i=xs; i<xe; i++) {
2606 ucat[k][j][i]=lucat[k][j-2][i];
2607 p[k][j][i]=lp[k][j-2][i];
2608 nvert[k][j][i]=lnvert[k][j-2][i];
2609 }
2610 }
2611 }
2612 }
2613
2614 if (user->bctype[3]==7){
2615 if (ye==my){
2616 j=my-1;
2617 for (k=zs; k<ze; k++) {
2618 for (i=xs; i<xe; i++) {
2619 ucat[k][j][i]=lucat[k][j+2][i];
2620 p[k][j][i]=lp[k][j+2][i];
2621 nvert[k][j][i]=lnvert[k][j+2][i];
2622 }
2623 }
2624 }
2625 }
2626
2627 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2628 DMDAVecRestoreArray(da, user->lP, &lp);
2629 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2630 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2631 DMDAVecRestoreArray(da, user->P, &p);
2632 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2633
2634/* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2635/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2636
2637/* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); */
2638/* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); */
2639
2640/* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2641/* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2642
2643 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2644 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2645
2646 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2647 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2648
2649 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2650 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2651
2652 //k-direction
2653 DMDAVecGetArray(fda, user->lUcat, &lucat);
2654 DMDAVecGetArray(da, user->lP, &lp);
2655 DMDAVecGetArray(da, user->lNvert, &lnvert);
2656 DMDAVecGetArray(fda, user->Ucat, &ucat);
2657 DMDAVecGetArray(da, user->P, &p);
2658 DMDAVecGetArray(da, user->Nvert, &nvert);
2659
2660 if (user->bctype[4]==7){
2661 if (zs==0){
2662 k=zs;
2663 for (j=ys; j<ye; j++) {
2664 for (i=xs; i<xe; i++) {
2665 ucat[k][j][i]=lucat[k-2][j][i];
2666 nvert[k][j][i]=lnvert[k-2][j][i];
2667 //amir
2668 p[k][j][i]=lp[k-2][j][i];
2669 }
2670 }
2671 }
2672 }
2673 if (user->bctype[5]==7){
2674 if (ze==mz){
2675 k=mz-1;
2676 for (j=ys; j<ye; j++) {
2677 for (i=xs; i<xe; i++) {
2678 ucat[k][j][i]=lucat[k+2][j][i];
2679 nvert[k][j][i]=lnvert[k+2][j][i];
2680 p[k][j][i]=lp[k+2][j][i];
2681 }
2682 }
2683 }
2684 }
2685
2686 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2687 DMDAVecRestoreArray(da, user->lP, &lp);
2688 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2689 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2690 DMDAVecRestoreArray(da, user->P, &p);
2691 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2692
2693/* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2694/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2695
2696/* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); */
2697/* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); */
2698
2699/* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2700/* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2701
2702 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2703 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2704
2705 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2706 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2707
2708 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2709 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2710 }
2711 /* ================================================================================== */
2712/* Analytical Vortex BC */
2713/* ================================================================================== */
2714
2715 DMDAVecGetArray(fda, user->Ucat, &ucat);
2716 DMDAVecGetArray(fda, user->Ucont, &ucont);
2717 DMDAVecGetArray(fda, user->Cent, &cent);
2718 DMDAVecGetArray(fda, user->Centx, &centx);
2719 DMDAVecGetArray(fda, user->Centy, &centy);
2720 DMDAVecGetArray(fda, user->Centz, &centz);
2721
2722 if (user->bctype[0]==9) {
2723 if (xs==0) {
2724 i= xs;
2725 for (k=lzs; k<lze; k++) {
2726 for (j=lys; j<lye; j++) {
2727 ucat[k][j][i].x=-cos(cent[k][j][i+1].x)*sin(cent[k][j][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2728 ucat[k][j][i].y=-sin(cent[k][j][i+1].x)*cos(cent[k][j][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2729 ucat[k][j][i].z =0.0;
2730
2731 ucont[k][j][i].x =-(cos(centx[k][j][i].x)*sin(centx[k][j][i].y)*csi[k][j][i].x)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2732 }
2733 }
2734 if (ys==0) {
2735 j=ys;
2736 for (k=lzs; k<lze; k++) {
2737 ucat[k][j][i].x=cos(cent[k][j+1][i+1].x)*sin(cent[k][j+1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2738 ucat[k][j][i].y=-sin(cent[k][j+1][i+1].x)*cos(cent[k][j+1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2739 ucat[k][j][i].z =0.0;
2740 }
2741 }
2742 if (ye==my) {
2743 j=ye-1;
2744 for (k=lzs; k<lze; k++) {
2745 ucat[k][j][i].x=cos(cent[k][j-1][i+1].x)*sin(cent[k][j-1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2746 ucat[k][j][i].y=-sin(cent[k][j-1][i+1].x)*cos(cent[k][j-1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2747 ucat[k][j][i].z =0.0;
2748 }
2749 }
2750 }
2751 }
2752 if (user->bctype[1]==9) {
2753 if (xe==mx) {
2754 i= xe-1;
2755 for (k=lzs; k<lze; k++) {
2756 for (j=lys; j<lye; j++) {
2757 ucat[k][j][i].x=-cos(cent[k][j][i-1].x)*sin(cent[k][j][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2758 ucat[k][j][i].y=-sin(cent[k][j][i-1].x)*cos(cent[k][j][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2759 ucat[k][j][i].z =0.0;
2760
2761 ucont[k][j][i-1].x =(-cos(centx[k][j][i-1].x)*sin(centx[k][j][i-1].y)*csi[k][j][i-1].x)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2762 }
2763 }
2764 if (ys==0) {
2765 j=ys;
2766 for (k=lzs; k<lze; k++) {
2767 ucat[k][j][i].x=cos(cent[k][j+1][i-1].x)*sin(cent[k][j+1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2768 ucat[k][j][i].y=-sin(cent[k][j+1][i-1].x)*cos(cent[k][j+1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2769 ucat[k][j][i].z =0.0;
2770 }
2771 }
2772 if (ye==my) {
2773 j=ye-1;
2774 for (k=lzs; k<lze; k++) {
2775 ucat[k][j][i].x=cos(cent[k][j-1][i-1].x)*sin(cent[k][j-1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2776 ucat[k][j][i].y=-sin(cent[k][j-1][i-1].x)*cos(cent[k][j-1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2777 ucat[k][j][i].z =0.0;
2778 }
2779 }
2780 }
2781 }
2782
2783 if (user->bctype[2]==9) {
2784 if (ys==0) {
2785 j= ys;
2786 for (k=lzs; k<lze; k++) {
2787 for (i=lxs; i<lxe; i++) {
2788 ucat[k][j][i].x=cos(cent[k][j+1][i].x)*sin(cent[k][j+1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2789 ucat[k][j][i].y=sin(cent[k][j+1][i].x)*cos(cent[k][j+1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2790 ucat[k][j][i].z =0.0;
2791
2792 ucont[k][j][i].y=(sin(centy[k][j][i].x)*cos(centy[k][j][i].y)*eta[k][j][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2793 }
2794 }
2795 }
2796 }
2797
2798
2799 if (user->bctype[3]==9) {
2800 if (ye==my) {
2801 j= ye-1;
2802 for (k=lzs; k<lze; k++) {
2803 for (i=lxs; i<lxe; i++) {
2804 ucat[k][j][i].x=cos(cent[k][j-1][i].x)*sin(cent[k][j-1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2805 ucat[k][j][i].y=sin(cent[k][j-1][i].x)*cos(cent[k][j-1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2806 ucat[k][j][i].z =0.0;
2807
2808 ucont[k][j-1][i].y=(sin(centy[k][j-1][i].x)*cos(centy[k][j-1][i].y)*eta[k][j-1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2809 }
2810 }
2811 }
2812 }
2813 if (user->bctype[4]==9) {
2814 if (zs==0) {
2815 k= zs;
2816 for (j=ys; j<ye; j++) {
2817 for (i=xs; i<xe; i++) {
2818 ucat[k][j][i].x=ucat[k+1][j][i].x;
2819 ucat[k][j][i].y=ucat[k+1][j][i].y;
2820 ucat[k][j][i].z=ucat[k+1][j][i].z;
2821
2822 ucont[k][j][i].z=0.0;
2823 }
2824 }
2825 }
2826 }
2827 if (user->bctype[5]==9) {
2828 if (ze==mz) {
2829 k= ze-1;
2830 for (j=ys; j<ye; j++) {
2831 for (i=xs; i<xe; i++) {
2832 ucat[k][j][i].x=ucat[k-1][j][i].x;
2833 ucat[k][j][i].y=ucat[k-1][j][i].y;
2834 ucat[k][j][i].z=ucat[k-1][j][i].z;
2835
2836 ucont[k-1][j][i].z=0.0;
2837 }
2838 }
2839 }
2840 }
2841
2842 DMDAVecRestoreArray(fda, user->Cent, &cent);
2843 DMDAVecRestoreArray(fda, user->Centx, &centx);
2844 DMDAVecRestoreArray(fda, user->Centy, &centy);
2845 DMDAVecRestoreArray(fda, user->Centz, &centz);
2846
2847 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2848 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
2849
2850 } // ttemp
2851
2852
2853 DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs);
2854 DMDAVecRestoreArray(fda, user->lCsi, &csi);
2855 DMDAVecRestoreArray(fda, user->lEta, &eta);
2856 DMDAVecRestoreArray(fda, user->lZet, &zet);
2857
2858 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2859 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2860 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2861 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2862 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2863 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2864
2865 PROFILE_FUNCTION_END;
2866
2867 PetscFunctionReturn(0);
2868 }
GetRandomCellAndLogicalCoordsOnInletFace
PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace(UserCtx *user, const DMDALocalInfo *info, PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, PetscRandom *rand_logic_i_ptr, PetscRandom *rand_logic_j_ptr, PetscRandom *rand_logic_k_ptr, PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out, PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out)
Assuming the current rank services the inlet face, this function selects a random cell (owned by this...
Definition Boundaries.c:459
BoundarySystem_ExecuteStep_Legacy
PetscErrorCode BoundarySystem_ExecuteStep_Legacy(UserCtx *user)
Acts as a temporary bridge to the legacy boundary condition implementation.
Definition Boundaries.c:673
OutflowFlux
PetscErrorCode OutflowFlux(UserCtx *user)
Calculates the total outgoing flux through all OUTLET faces for reporting.
Definition Boundaries.c:932
CanRankServiceInletFace
PetscErrorCode CanRankServiceInletFace(UserCtx *user, const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, PetscBool *can_service_inlet_out)
Determines if the current MPI rank owns any part of the globally defined inlet face,...
Definition Boundaries.c:25
InflowFlux
PetscErrorCode InflowFlux(UserCtx *user)
Applies inlet boundary conditions based on the modern BC handling system.
Definition Boundaries.c:706
CanRankServiceFace
PetscErrorCode CanRankServiceFace(const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, BCFace face_id, PetscBool *can_service_out)
Determines if the current MPI rank owns any part of a specified global face.
Definition Boundaries.c:150
GetDeterministicFaceGridLocation
PetscErrorCode GetDeterministicFaceGridLocation(UserCtx *user, const DMDALocalInfo *info, PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, PetscInt IM_cells_global, PetscInt JM_cells_global, PetscInt KM_cells_global, PetscInt64 particle_global_id, PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out, PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out, PetscBool *placement_successful_out)
Places particles in a deterministic grid/raster pattern on a specified domain face.
Definition Boundaries.c:267
FormBCS
PetscErrorCode FormBCS(UserCtx *user)
Definition Boundaries.c:1022
TranslateModernBCsToLegacy
PetscErrorCode TranslateModernBCsToLegacy(UserCtx *user)
Definition Boundaries.c:635
GetBCParamReal
PetscErrorCode GetBCParamReal(BC_Param *params, const char *key, PetscReal *value_out, PetscBool *found)
Searches a BC_Param linked list for a key and returns its value as a double.
Definition io.c:356
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
BCFaceToString
const char * BCFaceToString(BCFace face)
Helper function to convert BCFace enum to a string representation.
Definition logging.c:643
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:201
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
BCTypeToString
const char * BCTypeToString(BCType type)
Helper function to convert BCType enum to a string representation.
Definition logging.c:705
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:34
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2084
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global, PetscInt *xm_cell_local)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:1663
BCType
BCType
Defines the general mathematical/physical category of a boundary.
Definition variables.h:207
INLET
@ INLET
Definition variables.h:214
OUTLET
@ OUTLET
Definition variables.h:213
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:680
SimCtx::fish
PetscInt fish
Definition variables.h:572
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:679
SimCtx::block_number
PetscInt block_number
Definition variables.h:593
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:681
SimCtx::channelz
PetscInt channelz
Definition variables.h:593
UserCtx::lNvert
Vec lNvert
Definition variables.h:688
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:602
SimCtx::CMy_c
PetscReal CMy_c
Definition variables.h:589
UserCtx::Centz
Vec Centz
Definition variables.h:706
UserCtx::lZet
Vec lZet
Definition variables.h:705
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:231
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:232
SimCtx::ren
PetscReal ren
Definition variables.h:581
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:274
SimCtx::dt
PetscReal dt
Definition variables.h:552
SimCtx::inletprofile
PetscInt inletprofile
Definition variables.h:593
UserCtx::lUstar
Vec lUstar
Definition variables.h:683
SimCtx::np
PetscInt np
Definition variables.h:616
SimCtx::ccc
PetscInt ccc
Definition variables.h:605
SimCtx::ratio
PetscReal ratio
Definition variables.h:606
UserCtx::Ucont
Vec Ucont
Definition variables.h:688
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:548
BCS::Ubcs
Vec Ubcs
Boundary condition velocity values. (Comment: "An ugly hack, waste of memory")
Definition variables.h:121
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::Centx
Vec Centx
Definition variables.h:706
UserCtx::Bcs
BCS Bcs
Definition variables.h:682
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:602
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
SimCtx::CMz_c
PetscReal CMz_c
Definition variables.h:589
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:275
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Ucat
Vec Ucat
Definition variables.h:688
UserCtx::FluxIntfcSum
PetscReal FluxIntfcSum
Definition variables.h:685
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:610
SimCtx::rheology
PetscInt rheology
Definition variables.h:567
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:684
SimCtx::Fluxsum
PetscReal Fluxsum
Definition variables.h:602
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:604
UserCtx::lUcont
Vec lUcont
Definition variables.h:688
SimCtx::step
PetscInt step
Definition variables.h:546
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:603
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::FluxIntpSum
PetscReal FluxIntpSum
Definition variables.h:685
UserCtx::lEta
Vec lEta
Definition variables.h:705
UserCtx::Cent
Vec Cent
Definition variables.h:705
UserCtx::Nvert
Vec Nvert
Definition variables.h:688
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:273
UserCtx::Centy
Vec Centy
Definition variables.h:706
SimCtx::CMx_c
PetscReal CMx_c
Definition variables.h:589
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:122
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:200
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:201
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:203
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:202
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:203
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:201
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:202
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:271
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:538
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:661
sign
double sign(double a)
Definition wallfunction.c:227
noslip
void noslip(UserCtx *user, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, double nx, double ny, double nz)
Definition wallfunction.c:13
wall_function_loglaw
void wall_function_loglaw(UserCtx *user, double ks, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, PetscReal *ustar, double nx, double ny, double nz)
Definition wallfunction.c:276