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