PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
Boundaries.h File Reference
#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscsys.h>
#include <petscdmcomposite.h>
#include <petscsystypes.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "interpolation.h"
#include "AnalyticalSolution.h"
#include "ParticleMotion.h"
#include "BC_Handlers.h"
#include "wallfunction.h"
Include dependency graph for Boundaries.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

PetscErrorCode BoundarySystem_Create (UserCtx *user, const char *bcs_filename)
 Initializes the entire boundary system.
 
PetscErrorCode BoundarySystem_ExecuteStep (UserCtx *user)
 Executes one full boundary condition update cycle for a time step.
 
PetscErrorCode BoundarySystem_Destroy (UserCtx *user)
 Cleans up and destroys all boundary system resources.
 
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, making it responsible for placing particles on that portion of the surface.
 
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.
 
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 rank on that face) and random logical coordinates within that cell, suitable for placing a particle on the inlet surface.
 
PetscErrorCode TranslateModernBCsToLegacy (UserCtx *user)
 
PetscErrorCode InflowFlux (UserCtx *user)
 
PetscErrorCode OutflowFlux (UserCtx *user)
 
PetscErrorCode FormBCS (UserCtx *user)
 
PetscErrorCode BoundarySystem_ExecuteStep_Legacy (UserCtx *user)
 Acts as a temporary bridge to the legacy boundary condition implementation.
 

Function Documentation

◆ BoundarySystem_Create()

PetscErrorCode BoundarySystem_Create ( UserCtx *  user,
const char *  bcs_filename 
)

Initializes the entire boundary system.

Parameters
userThe main UserCtx struct.
bcs_filenameThe path to the boundary conditions configuration file.

◆ BoundarySystem_ExecuteStep()

PetscErrorCode BoundarySystem_ExecuteStep ( UserCtx *  user)

Executes one full boundary condition update cycle for a time step.

Parameters
userThe main UserCtx struct.

◆ BoundarySystem_Destroy()

PetscErrorCode BoundarySystem_Destroy ( UserCtx *  user)

Cleans up and destroys all boundary system resources.

Parameters
userThe main UserCtx struct.

◆ 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, making it responsible for placing particles on that portion of the surface.

The determination is based on the rank's owned nodes (from DMDALocalInfo) and the global node counts, in conjunction with the user->identifiedInletBCFace. A rank can service an inlet face if it owns the cells adjacent to that global boundary and has a non-zero extent (owns cells) in the tangential dimensions of that face.

Parameters
userPointer to the UserCtx structure, containing identifiedInletBCFace.
infoPointer to the DMDALocalInfo for the current rank's DA (node-based).
IM_nodes_globalGlobal number of nodes in the I-direction (e.g., user->IM + 1 if user->IM is cell count).
JM_nodes_globalGlobal number of nodes in the J-direction.
KM_nodes_globalGlobal number of nodes in the K-direction.
[out]can_service_inlet_outPointer to a PetscBool; set to PETSC_TRUE if the rank services (part of) the inlet, PETSC_FALSE otherwise.
Returns
PetscErrorCode 0 on success, non-zero on failure.

Definition at line 24 of file Boundaries.c.

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}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
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
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
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
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:649
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:650
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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.

This function is a general utility for parallel boundary operations. It checks if the local domain of the current MPI rank is adjacent to a specified global boundary face. A rank "services" a face if it owns the cells adjacent to that global boundary and has a non-zero extent (i.e., owns at least one cell) in the tangential dimensions of that face.

Parameters
infoPointer to the DMDALocalInfo for the current rank's DA.
face_idThe specific global face (e.g., BC_FACE_NEG_Z) to check.
[out]can_service_outPointer to a PetscBool; set to PETSC_TRUE if the rank services the face, PETSC_FALSE otherwise.
Returns
PetscErrorCode 0 on success.

Definition at line 130 of file Boundaries.c.

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}
Here is the call graph for this function:

◆ 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 rank on that face) and random logical coordinates within that cell, suitable for placing a particle on the inlet surface.

It is the caller's responsibility to ensure CanRankServiceInletFace returned true.

Parameters
userPointer to UserCtx.
infoPointer to DMDALocalInfo for the current rank (node-based).
xs_gnode,ys_gnode,zs_gnodeLocal starting node indices (incl. ghosts) for the rank's DA.
IM_nodes_global,JM_nodes_global,KM_nodes_globalGlobal node counts.
rand_logic_i_ptr,rand_logic_j_ptr,rand_logic_k_ptrPointers to RNGs for logical coords.
[out]ci_metric_lnode_out,cj_metric_lnode_out,ck_metric_lnode_outLocal node indices of the selected cell's origin (these are local to the rank's DA including ghosts).
[out]xi_metric_logic_out,eta_metric_logic_out,zta_metric_logic_outLogical coords [0,1] within the cell.
Returns
PetscErrorCode

Definition at line 225 of file Boundaries.c.

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}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ TranslateModernBCsToLegacy()

PetscErrorCode TranslateModernBCsToLegacy ( UserCtx *  user)

Definition at line 391 of file Boundaries.c.

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}
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
BCTypeToString
const char * BCTypeToString(BCType type)
Helper function to convert BCType enum to a string representation.
Definition logging.c:662
BCType
BCType
Defines the general mathematical/physical category of a boundary.
Definition variables.h:206
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:648
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:653
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:264
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:199
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InflowFlux()

PetscErrorCode InflowFlux ( UserCtx *  user)

Definition at line 819 of file Boundaries.c.

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}
InletInterpolation
static PetscReal InletInterpolation(PetscReal r2, UserCtx *user)
Definition Boundaries.c:456
M_PI
#define M_PI
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
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::lNvert
Vec lNvert
Definition variables.h:657
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:633
SimCtx::CMy_c
PetscReal CMy_c
Definition variables.h:557
UserCtx::lZet
Vec lZet
Definition variables.h:673
SimCtx::inletprofile
PetscInt inletprofile
Definition variables.h:562
UserCtx::Ucont
Vec Ucont
Definition variables.h:657
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::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::lUcont
Vec lUcont
Definition variables.h:657
UserCtx::info
DMDALocalInfo info
Definition variables.h:637
Cmpnts::y
PetscScalar y
Definition variables.h:100
UserCtx::lEta
Vec lEta
Definition variables.h:673
UserCtx::Cent
Vec Cent
Definition variables.h:673
SimCtx::CMx_c
PetscReal CMx_c
Definition variables.h:557
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ OutflowFlux()

PetscErrorCode OutflowFlux ( UserCtx *  user)

Definition at line 1165 of file Boundaries.c.

1165 {
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}
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:571
UserCtx::_this
PetscInt _this
Definition variables.h:643
Here is the caller graph for this function:

◆ FormBCS()

PetscErrorCode FormBCS ( UserCtx *  user)

Definition at line 1261 of file Boundaries.c.

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 }
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:1785
SimCtx::fish
PetscInt fish
Definition variables.h:540
SimCtx::block_number
PetscInt block_number
Definition variables.h:562
SimCtx::channelz
PetscInt channelz
Definition variables.h:562
UserCtx::Centz
Vec Centz
Definition variables.h:674
SimCtx::rstart_flg
PetscBool rstart_flg
Definition variables.h:528
SimCtx::ren
PetscReal ren
Definition variables.h:549
SimCtx::dt
PetscReal dt
Definition variables.h:527
UserCtx::lUstar
Vec lUstar
Definition variables.h:652
SimCtx::ccc
PetscInt ccc
Definition variables.h:574
SimCtx::ratio
PetscReal ratio
Definition variables.h:575
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:523
UserCtx::Centx
Vec Centx
Definition variables.h:674
UserCtx::FluxIntfcSum
PetscReal FluxIntfcSum
Definition variables.h:654
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:579
SimCtx::Fluxsum
PetscReal Fluxsum
Definition variables.h:571
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:573
SimCtx::step
PetscInt step
Definition variables.h:521
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:572
UserCtx::lAj
Vec lAj
Definition variables.h:673
UserCtx::lUcat
Vec lUcat
Definition variables.h:657
UserCtx::FluxIntpSum
PetscReal FluxIntpSum
Definition variables.h:654
UserCtx::Nvert
Vec Nvert
Definition variables.h:657
UserCtx::Centy
Vec Centy
Definition variables.h:674
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:121
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundarySystem_ExecuteStep_Legacy()

PetscErrorCode BoundarySystem_ExecuteStep_Legacy ( UserCtx *  user)

Acts as a temporary bridge to the legacy boundary condition implementation.

This function is the key to our integration strategy. It matches the signature of the modern BoundarySystem_ExecuteStep function that SetEulerianFields expects to call.

However, instead of containing new handler-based logic, it simply calls the monolithic legacy FormBCS function. This allows the modern orchestrator to drive the old solver logic without modification.

Parameters
userThe UserCtx for a single block.
Returns
PetscErrorCode

Definition at line 440 of file Boundaries.c.

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}
FormBCS
PetscErrorCode FormBCS(UserCtx *user)
Definition Boundaries.c:1261
Here is the call graph for this function: