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, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, BCFace face_id, PetscBool *can_service_out)
 Determines if the current MPI rank owns any part of a specified global face.
 
PetscErrorCode GetDeterministicFaceGridLocation (UserCtx *user, const DMDALocalInfo *info, PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, PetscInt IM_cells_global, PetscInt JM_cells_global, PetscInt KM_cells_global, PetscInt64 particle_global_id, PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out, PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out, PetscBool *placement_successful_out)
 Places particles in a deterministic grid/raster pattern on a specified domain face.
 
PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace (UserCtx *user, const DMDALocalInfo *info, PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, PetscRandom *rand_logic_i_ptr, PetscRandom *rand_logic_j_ptr, PetscRandom *rand_logic_k_ptr, PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out, PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out)
 Assuming the current rank services the inlet face, this function selects a random cell (owned by this 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)
 Applies inlet boundary conditions based on the modern BC handling system.
 
PetscErrorCode OutflowFlux (UserCtx *user)
 Calculates the total outgoing flux through all OUTLET faces for reporting.
 
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 25 of file Boundaries.c.

28{
29 PetscErrorCode ierr;
30 PetscMPIInt rank_for_logging; // For detailed debugging logs
31 PetscFunctionBeginUser;
32 PROFILE_FUNCTION_BEGIN;
33
34 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
35
36 *can_service_inlet_out = PETSC_FALSE; // Default to no service
37
38 if (!user->inletFaceDefined) {
39 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Inlet face not defined in user context. Cannot service.\n", rank_for_logging);
40 PetscFunctionReturn(0);
41 }
42
43 // Get the range of cells owned by this rank in each dimension
44 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
45 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
46 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
47
48 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
49 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
50 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
51
52 // Determine the global index of the last cell (0-indexed) in each direction.
53 // Example: If IM_nodes_global = 11 (nodes 0-10), there are 10 cells (0-9). Last cell index is 9.
54 // Formula: global_nodes - 1 (num cells) - 1 (0-indexed) = global_nodes - 2.
55 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1; // -1 if 0 or 1 node (i.e., 0 cells)
56 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
57 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
58
59 switch (user->identifiedInletBCFace) {
60 case BC_FACE_NEG_X: // Inlet on the global I-minimum face (face of cell C_i=0)
61 // Rank services if its first owned node is global node 0 (info->xs == 0),
62 // and it owns cells in I, J, and K directions.
63 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
64 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
65 *can_service_inlet_out = PETSC_TRUE;
66 }
67 break;
68 case BC_FACE_POS_X: // Inlet on the global I-maximum face (face of cell C_i=last_global_cell_idx_i)
69 // Rank services if it owns the last cell in I-direction,
70 // and has extent in J and K.
71 if (last_global_cell_idx_i >= 0 && /* Check for valid global domain */
72 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i && /* Rank's last cell is the global last cell */
73 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
74 *can_service_inlet_out = PETSC_TRUE;
75 }
76 break;
77 case BC_FACE_NEG_Y:
78 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
79 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
80 *can_service_inlet_out = PETSC_TRUE;
81 }
82 break;
83 case BC_FACE_POS_Y:
84 if (last_global_cell_idx_j >= 0 &&
85 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
86 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
87 *can_service_inlet_out = PETSC_TRUE;
88 }
89 break;
90 case BC_FACE_NEG_Z:
91 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
92 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
93 *can_service_inlet_out = PETSC_TRUE;
94 }
95 break;
96 case BC_FACE_POS_Z:
97 if (last_global_cell_idx_k >= 0 &&
98 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
99 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
100 *can_service_inlet_out = PETSC_TRUE;
101 }
102 break;
103 default:
104 LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d]: Unknown inlet face %s.\n", rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace));
105 break;
106 }
107
108
109 LOG_ALLOW(LOCAL, LOG_TRACE,
110 "[Rank %d] Check Service for Inlet %s:\n"
111 " - Local Domain: starts at cell (%d,%d,%d), has (%d,%d,%d) cells.\n"
112 " - Global Domain: has (%d,%d,%d) nodes, so last cell is (%d,%d,%d).\n",
113 rank_for_logging,
114 BCFaceToString((BCFace)user->identifiedInletBCFace),
115 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
116 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
117 IM_nodes_global, JM_nodes_global, KM_nodes_global,
118 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
119
120 LOG_ALLOW(LOCAL, LOG_DEBUG,"[Rank %d] Inlet Face %s Service Check Result: %s\n",
121 rank_for_logging,
122 BCFaceToString((BCFace)user->identifiedInletBCFace),
123 (*can_service_inlet_out) ? "CAN SERVICE" : "CANNOT SERVICE");
124
125 PROFILE_FUNCTION_END;
126
127 PetscFunctionReturn(0);
128}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
BCFaceToString
const char * BCFaceToString(BCFace face)
Helper function to convert BCFace enum to a string representation.
Definition logging.c:643
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:201
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:34
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global, PetscInt *xm_cell_local)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:1663
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:680
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:681
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:200
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:201
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:203
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:202
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:203
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:201
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:202
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CanRankServiceFace()

PetscErrorCode CanRankServiceFace ( const DMDALocalInfo *  info,
PetscInt  IM_nodes_global,
PetscInt  JM_nodes_global,
PetscInt  KM_nodes_global,
BCFace  face_id,
PetscBool *  can_service_out 
)

Determines if the current MPI rank owns any part of a specified global face.

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.
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.
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 150 of file Boundaries.c.

152{
153 PetscErrorCode ierr;
154 PetscMPIInt rank_for_logging;
155 PetscFunctionBeginUser;
156
157 PROFILE_FUNCTION_BEGIN;
158
159 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
160
161 *can_service_out = PETSC_FALSE; // Default to no service
162
163 // Get the range of cells owned by this rank
164 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
165 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
166 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
167 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
168 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
169 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
170
171 // Determine the global index of the last cell (0-indexed) in each direction.
172 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
173 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
174 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
175
176 switch (face_id) {
177 case BC_FACE_NEG_X:
178 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
179 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
180 *can_service_out = PETSC_TRUE;
181 }
182 break;
183 case BC_FACE_POS_X:
184 if (last_global_cell_idx_i >= 0 &&
185 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i &&
186 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
187 *can_service_out = PETSC_TRUE;
188 }
189 break;
190 case BC_FACE_NEG_Y:
191 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
192 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
193 *can_service_out = PETSC_TRUE;
194 }
195 break;
196 case BC_FACE_POS_Y:
197 if (last_global_cell_idx_j >= 0 &&
198 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
199 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
200 *can_service_out = PETSC_TRUE;
201 }
202 break;
203 case BC_FACE_NEG_Z:
204 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
205 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
206 *can_service_out = PETSC_TRUE;
207 }
208 break;
209 case BC_FACE_POS_Z:
210 if (last_global_cell_idx_k >= 0 &&
211 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
212 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
213 *can_service_out = PETSC_TRUE;
214 }
215 break;
216 default:
217 LOG_ALLOW(LOCAL, LOG_WARNING, "Rank %d: Unknown face enum %d. \n", rank_for_logging, face_id);
218 break;
219 }
220
221 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d check for face %s: Result=%s. \n",
222 rank_for_logging, BCFaceToString((BCFace)face_id), (*can_service_out ? "TRUE" : "FALSE"));
223
224 PROFILE_FUNCTION_END;
225
226 PetscFunctionReturn(0);
227}
Here is the call graph for this function:

◆ GetDeterministicFaceGridLocation()

PetscErrorCode GetDeterministicFaceGridLocation ( UserCtx *  user,
const DMDALocalInfo *  info,
PetscInt  xs_gnode_rank,
PetscInt  ys_gnode_rank,
PetscInt  zs_gnode_rank,
PetscInt  IM_cells_global,
PetscInt  JM_cells_global,
PetscInt  KM_cells_global,
PetscInt64  particle_global_id,
PetscInt *  ci_metric_lnode_out,
PetscInt *  cj_metric_lnode_out,
PetscInt *  ck_metric_lnode_out,
PetscReal *  xi_metric_logic_out,
PetscReal *  eta_metric_logic_out,
PetscReal *  zta_metric_logic_out,
PetscBool *  placement_successful_out 
)

Places particles in a deterministic grid/raster pattern on a specified domain face.

This function creates a set of equidistant, parallel lines of particles near the four edges of the face specified by user->identifiedInletBCFace. The number of lines drawn from each edge is hardcoded within this function (default is 2).

For example, if grid_layers=2 on face BC_FACE_NEG_X, the function will create particle lines at:

  • y ~ 0*dy, y ~ 1*dy (parallel to the Z-axis, starting from the J=0 edge)
  • y ~ y_max, y ~ y_max-dy (parallel to the Z-axis, starting from the J=max edge)
  • z ~ 0*dz, z ~ 1*dz (parallel to the Y-axis, starting from the K=0 edge)
  • z ~ z_max, z ~ z_max-dz (parallel to the Y-axis, starting from the K=max edge)

The particle's final position is set just inside the target cell face to ensure it is correctly located. The total number of particles (simCtx->np) is distributed as evenly as possible among all generated lines.

The function includes extensive validation to stop with an error if the requested grid placement is geometrically impossible (e.g., in a 2D domain or if layers would overlap). It also issues warnings for non-fatal but potentially unintended configurations.

Parameters
userPointer to UserCtx, which must contain a valid identifiedInletBCFace.
infoPointer to DMDALocalInfo for the current rank's grid layout.
xs_gnode_rank,ys_gnode_rank,zs_gnode_rankLocal starting node indices (incl. ghosts) for the rank's DA.
IM_cells_global,JM_cells_global,KM_cells_globalGlobal cell counts.
particle_global_idThe unique global ID of the particle being placed (from 0 to np-1).
[out]ci_metric_lnode_outLocal I-node index of the selected cell's origin.
[out]cj_metric_lnode_outLocal J-node index of the selected cell's origin.
[out]ck_metric_lnode_outLocal K-node index of the selected cell's origin.
[out]xi_metric_logic_outLogical xi-coordinate [0,1] within the cell.
[out]eta_metric_logic_outLogical eta-coordinate [0,1] within the cell.
[out]zta_metric_logic_outLogical zta-coordinate [0,1] within the cell.
[out]placement_successful_outPETSC_TRUE if the point belongs to this rank, PETSC_FALSE otherwise.
Returns
PetscErrorCode

Definition at line 267 of file Boundaries.c.

275{
276 SimCtx *simCtx = user->simCtx;
277 PetscReal global_logic_i, global_logic_j, global_logic_k;
278 PetscErrorCode ierr;
279 PetscMPIInt rank_for_logging;
280
281 PetscFunctionBeginUser;
282 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
283
284 *placement_successful_out = PETSC_FALSE; // Default to failure
285
286 // --- Step 1: Configuration and Input Validation ---
287
288 // *** Hardcoded number of grid layers. Change this value to alter the pattern. ***
289 const PetscInt grid_layers = 2;
290
291 LOG_ALLOW(LOCAL, LOG_DEBUG,
292 "[Rank %d] Placing particle %lld on face %s with grid_layers=%d in global domain (%d,%d,%d) cells.\n",
293 rank_for_logging, (long long)particle_global_id, BCFaceToString(user->identifiedInletBCFace), grid_layers,
294 IM_cells_global, JM_cells_global, KM_cells_global);
295
296 const char *face_name = BCFaceToString(user->identifiedInletBCFace);
297
298 // Fatal Error Checks: Ensure the requested grid is geometrically possible.
299 // The total layers from opposite faces (2 * grid_layers) must be less than the domain size.
300 switch (user->identifiedInletBCFace) {
301 case BC_FACE_NEG_X: case BC_FACE_POS_X:
302 if (JM_cells_global <= 1 || KM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (J-cells=%d, K-cells=%d).", face_name, JM_cells_global, KM_cells_global);
303 if (2 * grid_layers >= JM_cells_global || 2 * grid_layers >= KM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing J/K faces would overlap in this domain (J-cells=%d, K-cells=%d).", grid_layers, JM_cells_global, KM_cells_global);
304 break;
305 case BC_FACE_NEG_Y: case BC_FACE_POS_Y:
306 if (IM_cells_global <= 1 || KM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (I-cells=%d, K-cells=%d).", face_name, IM_cells_global, KM_cells_global);
307 if (2 * grid_layers >= IM_cells_global || 2 * grid_layers >= KM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing I/K faces would overlap in this domain (I-cells=%d, K-cells=%d).", grid_layers, IM_cells_global, KM_cells_global);
308 break;
309 case BC_FACE_NEG_Z: case BC_FACE_POS_Z:
310 if (IM_cells_global <= 1 || JM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (I-cells=%d, J-cells=%d).", face_name, IM_cells_global, JM_cells_global);
311 if (2 * grid_layers >= IM_cells_global || 2 * grid_layers >= JM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing I/J faces would overlap in this domain (I-cells=%d, J-cells=%d).", grid_layers, IM_cells_global, JM_cells_global);
312 break;
313 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid identifiedInletBCFace specified: %d", user->identifiedInletBCFace);
314 }
315
316 const PetscInt num_lines_total = 4 * grid_layers;
317 if (simCtx->np < num_lines_total) {
318 LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Total particle count (%lld) is less than the number of grid lines requested (%d). Some lines may be empty.\n", (long long)simCtx->np, num_lines_total);
319 }
320 if (simCtx->np > 0 && simCtx->np % num_lines_total != 0) {
321 LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Total particle count (%lld) is not evenly divisible by the number of grid lines (%d). Distribution will be uneven.\n", (long long)simCtx->np, num_lines_total);
322 }
323
324 // --- Step 2: Map global particle ID to a line and a point on that line ---
325 if (simCtx->np == 0) PetscFunctionReturn(0); // Nothing to do
326
327 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Distributing %lld particles over %d lines on face %s.\n",
328 rank_for_logging, (long long)simCtx->np, num_lines_total, face_name);
329
330 const PetscInt points_per_line = PetscMax(1, simCtx->np / num_lines_total);
331 PetscInt line_index = particle_global_id / points_per_line;
332 PetscInt point_index_on_line = particle_global_id % points_per_line;
333 line_index = PetscMin(line_index, num_lines_total - 1); // Clamp to handle uneven division
334
335 // Decode the line_index into an edge group (0-3) and a layer within that group (0 to grid_layers-1)
336 const PetscInt edge_group = line_index / grid_layers;
337 const PetscInt layer_index = line_index % grid_layers;
338
339 // --- Step 3: Calculate placement coordinates based on the decoded indices ---
340 const PetscReal epsilon = 1.0e-6; // Small offset to keep particles off exact cell boundaries
341 const PetscReal layer_spacing_norm_i = (IM_cells_global > 0) ? 1.0 / (PetscReal)IM_cells_global : 0.0;
342 const PetscReal layer_spacing_norm_j = (JM_cells_global > 0) ? 1.0 / (PetscReal)JM_cells_global : 0.0;
343 const PetscReal layer_spacing_norm_k = (KM_cells_global > 0) ? 1.0 / (PetscReal)KM_cells_global : 0.0;
344
345 PetscReal variable_coord; // The coordinate that varies along a line
346 if (points_per_line <= 1) {
347 variable_coord = 0.5; // Place single point in the middle
348 } else {
349 variable_coord = ((PetscReal)point_index_on_line + 0.5)/ (PetscReal)(points_per_line);
350 }
351 variable_coord = PetscMin(1.0 - epsilon, PetscMax(epsilon, variable_coord)); // Clamp within [eps, 1-eps]
352
353 // Main logic switch to determine the three global logical coordinates
354 switch (user->identifiedInletBCFace) {
355 case BC_FACE_NEG_X:
356 global_logic_i = 0.5 * layer_spacing_norm_i; // Place near the face, in the middle of the first cell
357 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
358 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
359 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
360 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
361 break;
362 case BC_FACE_POS_X:
363 global_logic_i = 1.0 - (0.5 * layer_spacing_norm_i); // Place near the face, in the middle of the last cell
364 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
365 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
366 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
367 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
368 break;
369 case BC_FACE_NEG_Y:
370 global_logic_j = 0.5 * layer_spacing_norm_j;
371 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
372 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
373 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
374 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
375 break;
376 case BC_FACE_POS_Y:
377 global_logic_j = 1.0 - (0.5 * layer_spacing_norm_j);
378 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
379 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
380 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
381 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
382 break;
383 case BC_FACE_NEG_Z:
384 global_logic_k = 0.5 * layer_spacing_norm_k;
385 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
386 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
387 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
388 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
389 break;
390 case BC_FACE_POS_Z:
391 global_logic_k = 1.0 - (0.5 * layer_spacing_norm_k);
392 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
393 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
394 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
395 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
396 break;
397 }
398
399 LOG_ALLOW(LOCAL, LOG_TRACE,
400 "[Rank %d] Particle %lld assigned to line %d (edge group %d, layer %d) with variable_coord=%.4f.\n"
401 " -> Global logical coords: (i,j,k) = (%.6f, %.6f, %.6f)\n",
402 rank_for_logging, (long long)particle_global_id, line_index, edge_group, layer_index, variable_coord,
403 global_logic_i, global_logic_j, global_logic_k);
404
405 // --- Step 4: Convert global logical coordinate to global cell index and intra-cell logicals ---
406 PetscReal global_cell_coord_i = global_logic_i * IM_cells_global;
407 PetscInt I_g = (PetscInt)global_cell_coord_i;
408 *xi_metric_logic_out = global_cell_coord_i - I_g;
409
410 PetscReal global_cell_coord_j = global_logic_j * JM_cells_global;
411 PetscInt J_g = (PetscInt)global_cell_coord_j;
412 *eta_metric_logic_out = global_cell_coord_j - J_g;
413
414 PetscReal global_cell_coord_k = global_logic_k * KM_cells_global;
415 PetscInt K_g = (PetscInt)global_cell_coord_k;
416 *zta_metric_logic_out = global_cell_coord_k - K_g;
417
418 // --- Step 5: Check if this rank owns the target cell and finalize outputs ---
419 if ((I_g >= info->xs && I_g < info->xs + info->xm) &&
420 (J_g >= info->ys && J_g < info->ys + info->ym) &&
421 (K_g >= info->zs && K_g < info->zs + info->zm))
422 {
423 // Convert global cell index to the local node index for this rank's DA patch
424 *ci_metric_lnode_out = (I_g - info->xs) + xs_gnode_rank;
425 *cj_metric_lnode_out = (J_g - info->ys) + ys_gnode_rank;
426 *ck_metric_lnode_out = (K_g - info->zs) + zs_gnode_rank;
427 *placement_successful_out = PETSC_TRUE;
428 }
429
430 LOG_ALLOW(LOCAL, LOG_VERBOSE,
431 "[Rank %d] Particle %lld placement %s. Local cell origin node: (I,J,K) = (%d,%d,%d), intra-cell logicals: (xi,eta,zta)=(%.6f,%.6f,%.6f)\n",
432 rank_for_logging, (long long)particle_global_id,
433 (*placement_successful_out ? "SUCCESSFUL" : "NOT ON THIS RANK"),
434 *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
435 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
436
437 PetscFunctionReturn(0);
438}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
SimCtx::np
PetscInt np
Definition variables.h:616
SimCtx
The master context for the entire simulation.
Definition variables.h:538
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetRandomCellAndLogicalCoordsOnInletFace()

PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace ( UserCtx *  user,
const DMDALocalInfo *  info,
PetscInt  xs_gnode_rank,
PetscInt  ys_gnode_rank,
PetscInt  zs_gnode_rank,
PetscInt  IM_nodes_global,
PetscInt  JM_nodes_global,
PetscInt  KM_nodes_global,
PetscRandom *  rand_logic_i_ptr,
PetscRandom *  rand_logic_j_ptr,
PetscRandom *  rand_logic_k_ptr,
PetscInt *  ci_metric_lnode_out,
PetscInt *  cj_metric_lnode_out,
PetscInt *  ck_metric_lnode_out,
PetscReal *  xi_metric_logic_out,
PetscReal *  eta_metric_logic_out,
PetscReal *  zta_metric_logic_out 
)

Assuming the current rank services the inlet face, this function selects a random cell (owned by this 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 459 of file Boundaries.c.

466{
467 PetscErrorCode ierr = 0;
468 PetscReal r_val_i_sel, r_val_j_sel, r_val_k_sel;
469 PetscInt local_cell_idx_on_face_dim1 = 0; // 0-indexed relative to owned cells on face
470 PetscInt local_cell_idx_on_face_dim2 = 0;
471 PetscMPIInt rank_for_logging;
472
473 PetscFunctionBeginUser;
474
475 PROFILE_FUNCTION_BEGIN;
476
477 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
478
479 // Defaults for cell origin node (local index for the rank's DA patch, including ghosts)
480 *ci_metric_lnode_out = xs_gnode_rank; *cj_metric_lnode_out = ys_gnode_rank; *ck_metric_lnode_out = zs_gnode_rank;
481 // Defaults for logical coordinates
482 *xi_metric_logic_out = 0.5; *eta_metric_logic_out = 0.5; *zta_metric_logic_out = 0.5;
483
484 // Get number of cells this rank owns in each dimension (tangential to the face mainly)
485 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
486 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
487 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
488
489 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
490 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
491 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
492
493 // Index of the last cell (0-indexed) in each global direction
494 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
495 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
496 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
497
498 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d: Inlet face %s. Owned cells(i,j,k):(%d,%d,%d). GlobNodes(I,J,K):(%d,%d,%d). Rank's DA node starts at (%d,%d,%d).\n",
499 rank_for_logging, user->identifiedInletBCFace, num_owned_cells_on_rank_i,num_owned_cells_on_rank_j,num_owned_cells_on_rank_k,
500 IM_nodes_global,JM_nodes_global,KM_nodes_global, xs_gnode_rank,ys_gnode_rank,zs_gnode_rank);
501
502
503 switch (user->identifiedInletBCFace) {
504 case BC_FACE_NEG_X: // Particle on -X face of cell C_0 (origin node N_0)
505 // Cell origin node is the first owned node in I by this rank (global index info->xs).
506 // Its local index within the rank's DA (incl ghosts) is xs_gnode_rank.
507 *ci_metric_lnode_out = xs_gnode_rank;
508 *xi_metric_logic_out = 1.0e-6;
509
510 // Tangential dimensions are J and K. Select an owned cell randomly on this face.
511 // num_owned_cells_on_rank_j/k must be > 0 (checked by CanRankServiceInletFace)
512 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
513 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j); // Index among owned J-cells
514 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
515 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1; // Offset from start of rank's J-nodes
516
517 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
518 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
519 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
520 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
521
522 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
523 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
524 break;
525
526 case BC_FACE_POS_X: // Particle on +X face of cell C_last_I (origin node N_last_I_origin)
527 // Origin node of the last I-cell is global_node_idx = last_global_cell_idx_i.
528 // Its local index in rank's DA: (last_global_cell_idx_i - info->xs) + xs_gnode_rank
529 *ci_metric_lnode_out = xs_gnode_rank + (last_global_cell_idx_i - info->xs);
530 *xi_metric_logic_out = 1.0 - 1.0e-6;
531
532 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
533 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
534 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
535 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1;
536
537 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
538 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
539 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
540 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
541
542 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
543 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
544 break;
545 // ... (Cases for Y and Z faces, following the same pattern) ...
546 case BC_FACE_NEG_Y:
547 *cj_metric_lnode_out = ys_gnode_rank;
548 *eta_metric_logic_out = 1.0e-6;
549 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
550 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
551 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
552 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
553 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
554 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
555 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
556 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
557 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
558 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
559 break;
560 case BC_FACE_POS_Y:
561 *cj_metric_lnode_out = ys_gnode_rank + (last_global_cell_idx_j - info->ys);
562 *eta_metric_logic_out = 1.0 - 1.0e-6;
563 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
564 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
565 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
566 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
567 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
568 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
569 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
570 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
571 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
572 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
573 break;
574 case BC_FACE_NEG_Z: // Your example case
575 *ck_metric_lnode_out = zs_gnode_rank; // Cell origin is the first owned node in K by this rank
576 *zta_metric_logic_out = 1.0e-6; // Place particle slightly inside this cell from its -Z face
577 // Tangential dimensions are I and J
578 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
579 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
580 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
581 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
582
583 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
584 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
585 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
586 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
587
588 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for I
589 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for J
590 break;
591 case BC_FACE_POS_Z:
592 *ck_metric_lnode_out = zs_gnode_rank + (last_global_cell_idx_k - info->zs);
593 *zta_metric_logic_out = 1.0 - 1.0e-6;
594 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
595 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
596 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
597 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
598 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
599 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
600 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
601 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
602 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
603 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
604 break;
605 default:
606 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "GetRandomCellAndLogicOnInletFace: Invalid user->identifiedInletBCFace %d. \n", user->identifiedInletBCFace);
607 }
608
609 PetscReal eps = 1.0e-7;
610 if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
611 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
612 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
613 } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
614 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
615 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
616 } else {
617 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
618 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
619 }
620
621 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Rank %d: Target Cell Node =(%d,%d,%d). (xi,et,zt)=(%.2e,%.2f,%.2f). \n",
622 rank_for_logging, *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
623 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
624
625 PROFILE_FUNCTION_END;
626
627 PetscFunctionReturn(0);
628}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ TranslateModernBCsToLegacy()

PetscErrorCode TranslateModernBCsToLegacy ( UserCtx *  user)

Definition at line 635 of file Boundaries.c.

636{
637 PetscFunctionBeginUser;
638 PROFILE_FUNCTION_BEGIN;
639 LOG_ALLOW(GLOBAL,LOG_DEBUG," Translating modern BC config to legacy integer codes...\n");
640
641 for (int i = 0; i < 6; i++) {
642 BCType modern_type = user->boundary_faces[i].mathematical_type;
643 user->bctype[i] = (int)modern_type;
644
645
646 BCFace current_face = (BCFace)i;
647 const char* face_str = BCFaceToString(current_face);
648 const char* bc_type_str = BCTypeToString(modern_type);
649 LOG_ALLOW(GLOBAL,LOG_TRACE," for face %s(%d), legacy type = %d & modern type = %s .\n",face_str,i,user->bctype[i],bc_type_str);
650 }
651 LOG_ALLOW(GLOBAL,LOG_DEBUG," -> Translation complete.\n");
652 PROFILE_FUNCTION_END;
653 PetscFunctionReturn(0);
654}
BCTypeToString
const char * BCTypeToString(BCType type)
Helper function to convert BCType enum to a string representation.
Definition logging.c:705
BCType
BCType
Defines the general mathematical/physical category of a boundary.
Definition variables.h:207
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:679
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:684
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:273
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InflowFlux()

PetscErrorCode InflowFlux ( UserCtx *  user)

Applies inlet boundary conditions based on the modern BC handling system.

This function iterates through all 6 domain faces. For each face identified as an INLET, it applies the velocity profile specified by its assigned handler and parameters (e.g., 'constant_velocity' with vx,vy,vz or 'parabolic' with u_max).

It calculates the contravariant flux (Ucont), Cartesian velocity on the face (Ubcs), and the staggered Cartesian velocity (Ucat). It also computes the total incoming flux and area across all MPI ranks.

Parameters
userThe main UserCtx struct containing the BC configuration and PETSc objects.
Returns
PetscErrorCode 0 on success.

Definition at line 706 of file Boundaries.c.

707{
708 PetscErrorCode ierr;
709 PetscReal lFluxIn = 0.0, lAreaIn = 0.0, AreaSumIn;
710 Vec lCoor;
711 Cmpnts ***ucont, ***ubcs, ***ucat, ***coor, ***csi, ***eta, ***zet, ***cent;
712 PetscReal ***nvert;
713
714 DM da = user->da, fda = user->fda;
715 DMDALocalInfo info = user->info;
716 PetscInt xs = info.xs, xe = info.xs + info.xm;
717 PetscInt ys = info.ys, ye = info.ys + info.ym;
718 PetscInt zs = info.zs, ze = info.zs + info.zm;
719 PetscInt mx = info.mx, my = info.my, mz = info.mz;
720 PetscInt lxs, lxe, lys, lye, lzs, lze;
721
722 // Get context for coordinate transformation if needed by a handler
723 SimCtx *simCtx = user->simCtx;
724 PetscReal CMx_c = simCtx->CMx_c;
725 PetscReal CMy_c = simCtx->CMy_c;
726 PetscReal CMz_c = simCtx->CMz_c;
727
728 PetscFunctionBeginUser;
729
730 PROFILE_FUNCTION_BEGIN;
731
732 lxs = xs; lxe = xe; lys = ys; lye = ye; lzs = zs; lze = ze;
733 if (xs==0) lxs = xs+1; if (ys==0) lys = ys+1; if (zs==0) lzs = zs+1;
734 if (xe==mx) lxe = xe-1; if (ye==my) lye = ye-1; if (ze==mz) lze = ze-1;
735
736 // --- Get PETSc arrays ---
737 DMGetCoordinatesLocal(da,&lCoor);
738 DMDAVecGetArray(fda, lCoor, &coor); // Use local coordinates
739 DMDAVecGetArray(fda, user->Ucont, &ucont);
740 DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs);
741 DMDAVecGetArray(fda, user->Ucat, &ucat);
742 DMDAVecGetArray(da, user->lNvert, &nvert);
743 DMDAVecGetArray(fda, user->Cent, &cent);
744 DMDAVecGetArray(fda, user->lCsi, &csi);
745 DMDAVecGetArray(fda, user->lEta, &eta);
746 DMDAVecGetArray(fda, user->lZet, &zet);
747
748 // --- Main Loop over all 6 faces ---
749 for (PetscInt fn=0; fn<6; fn++) {
750 BoundaryFaceConfig *face_config = &user->boundary_faces[fn];
751
752 if (face_config->mathematical_type != INLET) {
753 continue; // Skip non-inlet faces
754 }
755
756 // This processor only acts if it is on the boundary of the global domain
757 PetscBool is_on_boundary = ( (fn==0 && xs==0) || (fn==1 && xe==mx) ||
758 (fn==2 && ys==0) || (fn==3 && ye==my) ||
759 (fn==4 && zs==0) || (fn==5 && ze==mz) );
760 if (!is_on_boundary) continue;
761
762 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Applying INLET handler for face: %s \n", BCFaceToString(face_config->face_id));
763
764 // --- Loop over the specific face geometry ---
765 switch(face_config->face_id) {
766 case BC_FACE_NEG_X: // -Xi
767 case BC_FACE_POS_X: // +Xi
768 {
769 PetscReal sign = (face_config->face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
770 PetscInt i = (face_config->face_id == BC_FACE_NEG_X) ? xs : mx - 2;
771 for (PetscInt k=lzs; k<lze; k++) {
772 for (PetscInt j=lys; j<lye; j++) {
773 if ( (sign > 0 && nvert[k][j][i+1] > 0.1) || (sign < 0 && nvert[k][j][i] > 0.1) ) continue;
774
775 PetscReal uin_this_point = 0.0;
776 // --- Determine velocity based on the handler for this point ---
777 if (face_config->handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY) {
778 PetscBool found;
779 ierr = GetBCParamReal(face_config->params, "vx", &uin_this_point, &found); CHKERRQ(ierr);
780 } else if(face_config->handler_type== BC_HANDLER_INLET_PARABOLIC){
781 PetscBool found;
782 PetscReal umax,diameter=1.0;
783 ierr = GetBCParamReal(face_config->params,"u_max",&umax,&found); CHKERRQ(ierr);
784 ierr = GetBCParamReal(face_config->params,"diameter",&diameter,&found); CHKERRQ(ierr);
785
786 // Radius is in the YZ-plane for an X-face inlet
787 PetscReal yc = cent[k][j][i + (sign>0)].y - CMy_c;
788 PetscReal zc = cent[k][j][i + (sign>0)].z - CMz_c;
789 PetscReal r = sqrt(yc*yc + zc*zc);
790 PetscReal r_norm = 2.0 * r / diameter;
791 uin_this_point = umax * (1.0 - r_norm * r_norm);
792 if (r_norm > 1.0) uin_this_point = 0.0;
793 }
794 // Add other X-face handlers like 'else if (handler == ...)' here
795
796 // --- Apply the calculated velocity ---
797 PetscReal CellArea = sqrt(csi[k][j][i].z*csi[k][j][i].z + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].x*csi[k][j][i].x);
798 lAreaIn += CellArea;
799 ucont[k][j][i].x = sign * uin_this_point * CellArea;
800 lFluxIn += ucont[k][j][i].x;
801 ubcs[k][j][i + (sign < 0)].x = sign * uin_this_point * csi[k][j][i].x / CellArea;
802 ubcs[k][j][i + (sign < 0)].y = sign * uin_this_point * csi[k][j][i].y / CellArea;
803 ubcs[k][j][i + (sign < 0)].z = sign * uin_this_point * csi[k][j][i].z / CellArea;
804 ucat[k][j][i + (sign > 0)] = ubcs[k][j][i + (sign < 0)];
805 }
806 }
807 } break;
808
809 case BC_FACE_NEG_Y: // -Eta
810 case BC_FACE_POS_Y: // +Eta
811 {
812 PetscReal sign = (face_config->face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
813 PetscInt j = (face_config->face_id == BC_FACE_NEG_Y) ? ys : my - 2;
814 for (PetscInt k=lzs; k<lze; k++) {
815 for (PetscInt i=lxs; i<lxe; i++) {
816 if ( (sign > 0 && nvert[k][j+1][i] > 0.1) || (sign < 0 && nvert[k][j][i] > 0.1) ) continue;
817
818 PetscReal uin_this_point = 0.0;
819 if (face_config->handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY) {
820 PetscBool found;
821 ierr = GetBCParamReal(face_config->params, "vy", &uin_this_point, &found); CHKERRQ(ierr);
822 }else if(face_config->handler_type == BC_HANDLER_INLET_PARABOLIC){
823 PetscBool found;
824 PetscReal umax,diameter=1.0;
825 ierr = GetBCParamReal(face_config->params,"umax",&umax,&found); CHKERRQ(ierr);
826 ierr = GetBCParamReal(face_config->params,"diameter",&diameter,&found); CHKERRQ(ierr);
827
828 // Radius is in the XZ-plane for a Y-face inlet
829 PetscReal xc = cent[k][j + (sign>0)][i].x - CMx_c;
830 PetscReal zc = cent[k][j + (sign>0)][i].z - CMz_c;
831 PetscReal r = sqrt(xc*xc + zc*zc);
832 PetscReal r_norm = 2.0 * r / diameter;
833 uin_this_point = umax * (1.0 - r_norm * r_norm);
834 if (r_norm > 1.0) uin_this_point = 0.0;
835
836 }
837
838 PetscReal CellArea = sqrt(eta[k][j][i].z*eta[k][j][i].z + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].x*eta[k][j][i].x);
839 lAreaIn += CellArea;
840 ucont[k][j][i].y = sign * uin_this_point * CellArea;
841 lFluxIn += ucont[k][j][i].y;
842 ubcs[k][j + (sign < 0)][i].x = sign * uin_this_point * eta[k][j][i].x / CellArea;
843 ubcs[k][j + (sign < 0)][i].y = sign * uin_this_point * eta[k][j][i].y / CellArea;
844 ubcs[k][j + (sign < 0)][i].z = sign * uin_this_point * eta[k][j][i].z / CellArea;
845 ucat[k][j + (sign > 0)][i] = ubcs[k][j + (sign < 0)][i];
846 }
847 }
848 } break;
849
850 case BC_FACE_NEG_Z: // -Zeta
851 case BC_FACE_POS_Z: // +Zeta
852 {
853 PetscReal sign = (face_config->face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
854 PetscInt k = (face_config->face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
855 for (PetscInt j=lys; j<lye; j++) {
856 for (PetscInt i=lxs; i<lxe; i++) {
857 if ( (sign > 0 && nvert[k+1][j][i] > 0.1) || (sign < 0 && nvert[k][j][i] > 0.1) ) continue;
858
859 PetscReal uin_this_point = 0.0;
860 if (face_config->handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY) {
861 PetscBool found;
862 ierr = GetBCParamReal(face_config->params, "vz", &uin_this_point, &found); CHKERRQ(ierr);
863 } else if (face_config->handler_type == BC_HANDLER_INLET_PARABOLIC) {
864 PetscBool found;
865 PetscReal umax, diameter=1.0; // Default diameter for r_norm = 2*r
866 ierr = GetBCParamReal(face_config->params, "u_max", &umax, &found); CHKERRQ(ierr);
867 ierr = GetBCParamReal(face_config->params, "diameter", &diameter, &found); CHKERRQ(ierr); // Optional
868
869 // Radius in the XY-plane for a Z-face inlet.
870 PetscReal xc = cent[k + (sign>0)][j][i].x - CMx_c;
871 PetscReal yc = cent[k + (sign>0)][j][i].y - CMy_c;
872 PetscReal r = sqrt(xc*xc + yc*yc);
873 PetscReal r_norm = 2.0 * r / diameter; // Normalized radius
874 uin_this_point = umax * (1.0 - r_norm * r_norm);
875 if (r_norm > 1.0) uin_this_point = 0.0;
876 }
877
878 PetscReal CellArea = sqrt(zet[k][j][i].z*zet[k][j][i].z + zet[k][j][i].y*zet[k][j][i].y + zet[k][j][i].x*zet[k][j][i].x);
879 lAreaIn += CellArea;
880 ucont[k][j][i].z = sign * uin_this_point * CellArea;
881 lFluxIn += ucont[k][j][i].z;
882 ubcs[k + (sign < 0)][j][i].x = sign * uin_this_point * zet[k][j][i].x / CellArea;
883 ubcs[k + (sign < 0)][j][i].y = sign * uin_this_point * zet[k][j][i].y / CellArea;
884 ubcs[k + (sign < 0)][j][i].z = sign * uin_this_point * zet[k][j][i].z / CellArea;
885 ucat[k + (sign > 0)][j][i] = ubcs[k + (sign < 0)][j][i];
886 }
887 }
888 } break;
889 } // end switch(face_id)
890 } // end for(fn)
891
892 // --- Finalize: Sum flux and area from all processes ---
893 ierr = MPI_Allreduce(&lFluxIn, &simCtx->FluxInSum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
894 ierr = MPI_Allreduce(&lAreaIn, &AreaSumIn, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
895
896 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Inflow: Flux - Area: %le - %le \n", simCtx->FluxInSum, AreaSumIn);
897
898 // --- Restore PETSc arrays ---
899 DMDAVecRestoreArray(fda, lCoor, &coor);
900 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
901 DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs);
902 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
903 DMDAVecRestoreArray(da, user->lNvert, &nvert);
904 DMDAVecRestoreArray(fda, user->Cent, &cent);
905 DMDAVecRestoreArray(fda, user->lCsi, &csi);
906 DMDAVecRestoreArray(fda, user->lEta, &eta);
907 DMDAVecRestoreArray(fda, user->lZet, &zet);
908
909 // --- Update local vectors for subsequent computations ---
910 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
911 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
912
913 PROFILE_FUNCTION_END;
914
915 PetscFunctionReturn(0);
916}
GetBCParamReal
PetscErrorCode GetBCParamReal(BC_Param *params, const char *key, PetscReal *value_out, PetscBool *found)
Searches a BC_Param linked list for a key and returns its value as a double.
Definition io.c:356
INLET
@ INLET
Definition variables.h:214
UserCtx::lNvert
Vec lNvert
Definition variables.h:688
SimCtx::CMy_c
PetscReal CMy_c
Definition variables.h:589
UserCtx::lZet
Vec lZet
Definition variables.h:705
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:231
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:232
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:274
UserCtx::Ucont
Vec Ucont
Definition variables.h:688
BCS::Ubcs
Vec Ubcs
Boundary condition velocity values. (Comment: "An ugly hack, waste of memory")
Definition variables.h:121
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::Bcs
BCS Bcs
Definition variables.h:682
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:602
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
SimCtx::CMz_c
PetscReal CMz_c
Definition variables.h:589
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:275
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Ucat
Vec Ucat
Definition variables.h:688
UserCtx::lUcont
Vec lUcont
Definition variables.h:688
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:705
UserCtx::Cent
Vec Cent
Definition variables.h:705
SimCtx::CMx_c
PetscReal CMx_c
Definition variables.h:589
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:271
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
sign
double sign(double a)
Definition wallfunction.c:227
Here is the call graph for this function:
Here is the caller graph for this function:

◆ OutflowFlux()

PetscErrorCode OutflowFlux ( UserCtx *  user)

Calculates the total outgoing flux through all OUTLET faces for reporting.

NOTE: In a mixed modern/legacy environment, this function is for DIAGNOSTICS ONLY. It reads the contravariant velocities and calculates the total flux passing through faces marked as OUTLET. It does NOT apply any boundary conditions itself, as that is still the responsibility of the legacy FormBCS function.

Parameters
userThe main UserCtx struct containing BC config and PETSc vectors.
Returns
PetscErrorCode 0 on success.

Definition at line 932 of file Boundaries.c.

932 {
933
934 PetscErrorCode ierr;
935 PetscReal lFluxOut = 0.0;
936 Cmpnts ***ucont;
937
938 DM fda = user->fda;
939 DMDALocalInfo info = user->info;
940 PetscInt xs = info.xs, xe = info.xs + info.xm;
941 PetscInt ys = info.ys, ye = info.ys + info.ym;
942 PetscInt zs = info.zs, ze = info.zs + info.zm;
943 PetscInt mx = info.mx, my = info.my, mz = info.mz;
944
945 PetscFunctionBeginUser;
946
947 PROFILE_FUNCTION_BEGIN;
948
949 ierr = DMDAVecGetArrayRead(fda, user->Ucont, &ucont); CHKERRQ(ierr);
950
951 // --- Loop over all 6 faces to find OUTLETS ---
952 for (PetscInt fn = 0; fn < 6; fn++) {
953 if (user->boundary_faces[fn].mathematical_type != OUTLET) {
954 continue;
955 }
956
957 PetscBool is_on_boundary = ( (fn==0 && xs==0) || (fn==1 && xe==mx) ||
958 (fn==2 && ys==0) || (fn==3 && ye==my) ||
959 (fn==4 && zs==0) || (fn==5 && ze==mz) );
960 if (!is_on_boundary) continue;
961
962 // --- Sum the flux for the appropriate face and component ---
963 switch ((BCFace)fn) {
964 case BC_FACE_NEG_X: case BC_FACE_POS_X: {
965 PetscInt i = (fn == 0) ? xs : mx - 2;
966 for (PetscInt k=info.zs; k<info.zs+info.zm; k++) for (PetscInt j=info.ys; j<info.ys+info.ym; j++) {
967 lFluxOut += ucont[k][j][i].x;
968 }
969 } break;
970
971 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: {
972 PetscInt j = (fn == 2) ? ys : my - 2;
973 for (PetscInt k=info.zs; k<info.zs+info.zm; k++) for (PetscInt i=info.xs; i<info.xs+info.xm; i++) {
974 lFluxOut += ucont[k][j][i].y;
975 }
976 } break;
977
978 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: {
979 PetscInt k = (fn == 4) ? zs : mz - 2;
980 for (PetscInt j=info.ys; j<info.ys+info.ym; j++) for (PetscInt i=info.xs; i<info.xs+info.xm; i++) {
981 lFluxOut += ucont[k][j][i].z;
982 }
983 } break;
984 } // end switch
985 } // end for loop
986
987 ierr = DMDAVecRestoreArrayRead(fda, user->Ucont, &ucont); CHKERRQ(ierr);
988
989 // --- Finalize: Sum and store the global total flux ---
990 ierr = MPI_Allreduce(&lFluxOut, &user->simCtx->FluxOutSum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
991
992 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Reported Global FluxOutSum = %.6f\n", user->simCtx->FluxOutSum);
993
994 PROFILE_FUNCTION_END;
995
996 PetscFunctionReturn(0);
997}
OUTLET
@ OUTLET
Definition variables.h:213
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:602
Here is the caller graph for this function:

◆ FormBCS()

PetscErrorCode FormBCS ( UserCtx *  user)

Definition at line 1022 of file Boundaries.c.

1023{
1024 DM da = user->da, fda = user->fda;
1025 DMDALocalInfo info = user->info;
1026 PetscInt xs = info.xs, xe = info.xs + info.xm;
1027 PetscInt ys = info.ys, ye = info.ys + info.ym;
1028 PetscInt zs = info.zs, ze = info.zs + info.zm;
1029 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1030 PetscInt lxs, lxe, lys, lye, lzs, lze;
1031 PetscInt i, j, k;
1032
1033 PetscReal ***nvert,***lnvert; //local working array
1034
1035 PetscReal ***p,***lp;
1036 Cmpnts ***ucont, ***ubcs, ***ucat,***lucat, ***csi, ***eta, ***zet;
1037 Cmpnts ***cent,***centx,***centy,***centz,***coor;
1038 PetscScalar FluxIn, FluxOut,ratio;
1039 PetscScalar lArea, AreaSum;
1040
1041 PetscScalar FarFluxIn=0., FarFluxOut=0., FarFluxInSum, FarFluxOutSum;
1042 PetscScalar FarAreaIn=0., FarAreaOut=0., FarAreaInSum, FarAreaOutSum;
1043 PetscScalar FluxDiff, VelDiffIn, VelDiffOut;
1044 Cmpnts V_frame;
1045
1046 PetscReal Un, nx,ny,nz,A;
1047
1048 SimCtx *simCtx = user->simCtx;
1049
1050 lxs = xs; lxe = xe;
1051 lys = ys; lye = ye;
1052 lzs = zs; lze = ze;
1053
1054 PetscInt gxs, gxe, gys, gye, gzs, gze;
1055
1056 gxs = info.gxs; gxe = gxs + info.gxm;
1057 gys = info.gys; gye = gys + info.gym;
1058 gzs = info.gzs; gze = gzs + info.gzm;
1059
1060 if (xs==0) lxs = xs+1;
1061 if (ys==0) lys = ys+1;
1062 if (zs==0) lzs = zs+1;
1063
1064 if (xe==mx ) lxe = xe-1;
1065 if (ye==my ) lye = ye-1;
1066 if (ze==mz ) lze = ze-1;
1067
1068 PetscFunctionBeginUser;
1069
1070 PROFILE_FUNCTION_BEGIN;
1071
1072 DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs);
1073
1074 DMDAVecGetArray(fda, user->lCsi, &csi);
1075 DMDAVecGetArray(fda, user->lEta, &eta);
1076 DMDAVecGetArray(fda, user->lZet, &zet);
1077
1078 PetscInt ttemp;
1079 for (ttemp=0; ttemp<3; ttemp++) {
1080 DMDAVecGetArray(da, user->Nvert, &nvert);
1081 DMDAVecGetArray(fda, user->lUcat, &ucat);
1082 DMDAVecGetArray(fda, user->Ucont, &ucont);
1083/* ================================================================================== */
1084/* FAR-FIELD BC */
1085/* ================================================================================== */
1086
1087
1088 // reset FAR FLUXES
1089 FarFluxIn = 0.; FarFluxOut=0.;
1090 FarAreaIn = 0.; FarAreaOut=0.;
1091
1092 PetscReal lFlux_abs=0.0,FluxSum_abs=0.0,ratio=0.0;
1093
1094 V_frame.x=0.;
1095 V_frame.y=0.;
1096 V_frame.z=0.;
1097
1098
1099
1100 if (user->bctype[0]==6) {
1101 if (xs == 0) {
1102 i= xs;
1103 for (k=lzs; k<lze; k++) {
1104 for (j=lys; j<lye; j++) {
1105 ubcs[k][j][i].x = ucat[k][j][i+1].x;
1106 ubcs[k][j][i].y = ucat[k][j][i+1].y;
1107 ubcs[k][j][i].z = ucat[k][j][i+1].z;
1108 ucont[k][j][i].x = ubcs[k][j][i].x * csi[k][j][i].x;
1109 FarFluxIn += ucont[k][j][i].x;
1110 lFlux_abs += fabs(ucont[k][j][i].x);
1111 FarAreaIn += csi[k][j][i].x;
1112 }
1113 }
1114 }
1115 }
1116
1117 if (user->bctype[1]==6) {
1118 if (xe==mx) {
1119 i= xe-1;
1120 for (k=lzs; k<lze; k++) {
1121 for (j=lys; j<lye; j++) {
1122 ubcs[k][j][i].x = ucat[k][j][i-1].x;
1123 ubcs[k][j][i].y = ucat[k][j][i-1].y;
1124 ubcs[k][j][i].z = ucat[k][j][i-1].z;
1125 ucont[k][j][i-1].x = ubcs[k][j][i].x * csi[k][j][i-1].x;
1126 FarFluxOut += ucont[k][j][i-1].x;
1127 lFlux_abs += fabs(ucont[k][j][i-1].x);
1128 FarAreaOut += csi[k][j][i-1].x;
1129 }
1130 }
1131 }
1132 }
1133
1134 if (user->bctype[2]==6) {
1135 if (ys==0) {
1136 j= ys;
1137 for (k=lzs; k<lze; k++) {
1138 for (i=lxs; i<lxe; i++) {
1139 ubcs[k][j][i].x = ucat[k][j+1][i].x;
1140 ubcs[k][j][i].y = ucat[k][j+1][i].y;
1141 ubcs[k][j][i].z = ucat[k][j+1][i].z;
1142 ucont[k][j][i].y = ubcs[k][j][i].y * eta[k][j][i].y;
1143 FarFluxIn += ucont[k][j][i].y;
1144 lFlux_abs += fabs(ucont[k][j][i].y);
1145 FarAreaIn += eta[k][j][i].y;
1146 }
1147 }
1148 }
1149 }
1150
1151 if (user->bctype[3]==6) {
1152 if (ye==my) {
1153 j=ye-1;
1154 for (k=lzs; k<lze; k++) {
1155 for (i=lxs; i<lxe; i++) {
1156 ubcs[k][j][i].x = ucat[k][j-1][i].x;
1157 ubcs[k][j][i].y = ucat[k][j-1][i].y;
1158 ubcs[k][j][i].z = ucat[k][j-1][i].z;
1159 ucont[k][j-1][i].y = ubcs[k][j][i].y * eta[k][j-1][i].y;
1160 FarFluxOut += ucont[k][j-1][i].y;
1161 lFlux_abs += fabs(ucont[k][j-1][i].y);
1162 FarAreaOut += eta[k][j-1][i].y;
1163 }
1164 }
1165 }
1166 }
1167
1168 if (user->bctype[4]==6 || user->bctype[4]==10) {
1169 if (zs==0) {
1170 k = 0;
1171 for (j=lys; j<lye; j++) {
1172 for (i=lxs; i<lxe; i++) {
1173 ubcs[k][j][i].x = ucat[k+1][j][i].x;
1174 ubcs[k][j][i].y = ucat[k+1][j][i].y;
1175 ubcs[k][j][i].z = ucat[k+1][j][i].z;
1176 ucont[k][j][i].z = ubcs[k][j][i].z * zet[k][j][i].z;
1177 FarFluxIn += ucont[k][j][i].z;
1178 lFlux_abs += fabs(ucont[k][j][i].z);
1179 FarAreaIn += zet[k][j][i].z;
1180 }
1181 }
1182 }
1183 }
1184
1185 if (user->bctype[5]==6 || user->bctype[5]==10) {
1186 if (ze==mz) {
1187 k = ze-1;
1188 for (j=lys; j<lye; j++) {
1189 for (i=lxs; i<lxe; i++) {
1190 ubcs[k][j][i].x = ucat[k-1][j][i].x;
1191 ubcs[k][j][i].y = ucat[k-1][j][i].y;
1192 ubcs[k][j][i].z = ucat[k-1][j][i].z;
1193 ucont[k-1][j][i].z = ubcs[k][j][i].z * zet[k-1][j][i].z;
1194 FarFluxOut += ucont[k-1][j][i].z;
1195 lFlux_abs += fabs(ucont[k-1][j][i].z);
1196 FarAreaOut += zet[k-1][j][i].z;
1197 }
1198 }
1199 }
1200 }
1201
1202 MPI_Allreduce(&FarFluxIn,&FarFluxInSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1203 MPI_Allreduce(&FarFluxOut,&FarFluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1204 MPI_Allreduce(&lFlux_abs,&FluxSum_abs,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1205 MPI_Allreduce(&FarAreaIn,&FarAreaInSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1206 MPI_Allreduce(&FarAreaOut,&FarAreaOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1207
1208 if (user->bctype[5]==6 || user->bctype[3]==6 || user->bctype[1]==6) {
1209
1210 ratio=(FarFluxInSum - FarFluxOutSum)/FluxSum_abs;
1211 if (fabs(FluxSum_abs) <1.e-10) ratio = 0.;
1212
1213 LOG_ALLOW(GLOBAL,LOG_DEBUG, "/FluxSum_abs %le ratio %le \n", FluxSum_abs,ratio);
1214
1215 FluxDiff = 0.5*(FarFluxInSum - FarFluxOutSum) ;
1216 VelDiffIn = FluxDiff / FarAreaInSum ;
1217
1218 if (fabs(FarAreaInSum) <1.e-6) VelDiffIn = 0.;
1219
1220 VelDiffOut = FluxDiff / FarAreaOutSum ;
1221
1222 if (fabs(FarAreaOutSum) <1.e-6) VelDiffOut = 0.;
1223
1224 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Far Flux Diff %d %le %le %le %le %le %le %le\n", simCtx->step, FarFluxInSum, FarFluxOutSum, FluxDiff, FarAreaInSum, FarAreaOutSum, VelDiffIn, VelDiffOut);
1225
1226 }
1227
1228 if (user->bctype[5]==10) {
1229 FluxDiff = simCtx->FluxInSum -( FarFluxOutSum -FarFluxInSum) ;
1230 VelDiffIn = 1/3.*FluxDiff / (FarAreaInSum);// + FarsimCtx->AreaOutSum);
1231 if (fabs(FarAreaInSum) <1.e-6) VelDiffIn = 0.;
1232
1233 VelDiffOut = 2./3.* FluxDiff / (FarAreaOutSum) ;//(FarAreaInSum + FarsimCtx->AreaOutSum) ;
1234
1235 if (fabs(FarAreaOutSum) <1.e-6) VelDiffOut = 0.;
1236
1237 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Far Flux Diff %d %le %le %le %le %le %le %le\n", simCtx->step, FarFluxInSum, FarFluxOutSum, FluxDiff, FarAreaInSum, FarAreaOutSum, VelDiffIn, VelDiffOut);
1238
1239 }
1240
1241
1242 // scale global mass conservation
1243
1244 if (user->bctype[5]==6 || user->bctype[5]==10) {
1245 if (ze==mz) {
1246 k = ze-1;
1247 for (j=lys; j<lye; j++) {
1248 for (i=lxs; i<lxe; i++) {
1249 ucont[k-1][j][i].z += ratio*fabs(ucont[k-1][j][i].z);
1250 ubcs[k][j][i].z = ucont[k-1][j][i].z/zet[k-1][j][i].z;
1251 // ubcs[k][j][i].z = ucat[k-1][j][i].z + VelDiffOut ;//+ V_frame.z;
1252 // ucont[k-1][j][i].z = ubcs[k][j][i].z * zet[k-1][j][i].z;
1253 }
1254 }
1255 }
1256 }
1257
1258 if (user->bctype[3]==6) {
1259 if (ye==my) {
1260 j=ye-1;
1261 for (k=lzs; k<lze; k++) {
1262 for (i=lxs; i<lxe; i++) {
1263
1264 ucont[k][j-1][i].y +=ratio*fabs(ucont[k][j-1][i].y);
1265 ubcs[k][j][i].y = ucont[k][j-1][i].y /eta[k][j-1][i].y;
1266 // ubcs[k][j][i].y = ucat[k][j-1][i].y + VelDiffOut;// + V_frame.y;
1267 // ucont[k][j-1][i].y = ubcs[k][j][i].y * eta[k][j-1][i].y;
1268
1269 }
1270 }
1271 }
1272 }
1273
1274 if (user->bctype[1]==6) {
1275 if (xe==mx) {
1276 i= xe-1;
1277 for (k=lzs; k<lze; k++) {
1278 for (j=lys; j<lye; j++) {
1279 ucont[k][j][i-1].x +=ratio*fabs(ucont[k][j][i-1].x);
1280 ubcs[k][j][i].x = ucont[k][j][i-1].x / csi[k][j][i-1].x ;
1281 // ubcs[k][j][i].x = ucat[k][j][i-1].x + VelDiffOut;// + V_frame.x;
1282 // ucont[k][j][i-1].x = ubcs[k][j][i].x * csi[k][j][i-1].x;
1283 }
1284 }
1285 }
1286 }
1287
1288
1289 if (user->bctype[0]==6) {
1290 if (xs == 0) {
1291 i= xs;
1292 for (k=lzs; k<lze; k++) {
1293 for (j=lys; j<lye; j++) {
1294 ucont[k][j][i].x -=ratio*fabs(ucont[k][j][i].x);
1295 ubcs[k][j][i].x = ucont[k][j][i].x / csi[k][j][i].x;
1296 // ubcs[k][j][i].x = ucat[k][j][i+1].x - VelDiffIn;// + V_frame.x;
1297 // ucont[k][j][i].x = ubcs[k][j][i].x * csi[k][j][i].x;
1298 }
1299 }
1300 }
1301 }
1302
1303
1304 if (user->bctype[2]==6) {
1305 if (ys==0) {
1306 j= ys;
1307 for (k=lzs; k<lze; k++) {
1308 for (i=lxs; i<lxe; i++) {
1309 ucont[k][j][i].y -=ratio*fabs(ucont[k][j][i].y);
1310 ubcs[k][j][i].y = ucont[k][j][i].y / eta[k][j][i].y;
1311 // ubcs[k][j][i].y = ucat[k][j+1][i].y - VelDiffIn;// + V_frame.y;
1312 // ucont[k][j][i].y = ubcs[k][j][i].y * eta[k][j][i].y;
1313 }
1314 }
1315 }
1316 }
1317
1318
1319 if (user->bctype[4]==6 || user->bctype[5]==10) {
1320 if (zs==0) {
1321 k = 0;
1322 for (j=lys; j<lye; j++) {
1323 for (i=lxs; i<lxe; i++) {
1324 ucont[k][j][i].z -=ratio*fabs(ucont[k][j][i].z);
1325 ubcs[k][j][i].z =ucont[k][j][i].z / zet[k][j][i].z;
1326 // ubcs[k][j][i].z = ucat[k+1][j][i].z - VelDiffIn;// + V_frame.z;
1327 // ucont[k][j][i].z = ubcs[k][j][i].z * zet[k][j][i].z;
1328
1329 }
1330 }
1331 }
1332 }
1333
1334//// Amir wall Ogrid
1335
1336if (user->bctype[2]==1 || user->bctype[2]==-1) {
1337 if (ys==0) {
1338 j= ys;
1339 for (k=lzs; k<lze; k++) {
1340 for (i=lxs; i<lxe; i++) {
1341 A=sqrt(eta[k][j][i].z*eta[k][j][i].z +
1342 eta[k][j][i].y*eta[k][j][i].y +
1343 eta[k][j][i].x*eta[k][j][i].x);
1344 nx=eta[k][j][i].x/A;
1345 ny=eta[k][j][i].y/A;
1346 nz=eta[k][j][i].z/A;
1347 Un=ucat[k][j+1][i].x*nx+ucat[k][j+1][i].y*ny+ucat[k][j+1][i].z*nz;
1348 ubcs[k][j][i].x = 0.0;
1349 ubcs[k][j][i].y = 0.0;
1350 ubcs[k][j][i].z = 0.0;
1351 ucont[k][j][i].y = 0.;
1352 }
1353 }
1354 }
1355 }
1356
1357/* ================================================================================== */
1358/* SOLID WALL BC (NO-SLIP / NO-PENETRATION) */
1359/* ================================================================================== */
1360
1361// NOTE: This block is added to explicitly handle bctype=1 (solid wall) for all faces.
1362// It ensures both no-slip (ubcs=0) and no-penetration (ucont_normal=0).
1363// ubcs is handled by the implicit-zero assumption, but ucont must be set explicitly.
1364
1365// -X Face (i=0)
1366if (user->bctype[0]==1 || user->bctype[0]==-1) {
1367 if (xs==0) {
1368 i= xs;
1369 for (k=lzs; k<lze; k++) {
1370 for (j=lys; j<lye; j++) {
1371 ubcs[k][j][i].x = 0.0;
1372 ubcs[k][j][i].y = 0.0;
1373 ubcs[k][j][i].z = 0.0;
1374 ucont[k][j][i].x = 0.0; // Enforce no-penetration
1375 }
1376 }
1377 }
1378}
1379
1380// +X Face (i=mx-1)
1381if (user->bctype[1]==1 || user->bctype[1]==-1) {
1382 if (xe==mx) {
1383 i= xe-1;
1384 for (k=lzs; k<lze; k++) {
1385 for (j=lys; j<lye; j++) {
1386 ubcs[k][j][i].x = 0.0;
1387 ubcs[k][j][i].y = 0.0;
1388 ubcs[k][j][i].z = 0.0;
1389 // The relevant ucont is at the face, index i-1
1390 ucont[k][j][i-1].x = 0.0; // Enforce no-penetration
1391 }
1392 }
1393 }
1394}
1395
1396// -Y Face (j=0)
1397if (user->bctype[2]==1 || user->bctype[2]==-1) {
1398 if (ys==0) {
1399 j= ys;
1400 for (k=lzs; k<lze; k++) {
1401 for (i=lxs; i<lxe; i++) {
1402 ubcs[k][j][i].x = 0.0;
1403 ubcs[k][j][i].y = 0.0;
1404 ubcs[k][j][i].z = 0.0;
1405 ucont[k][j][i].y = 0.0; // Enforce no-penetration
1406 }
1407 }
1408 }
1409}
1410
1411// +Y Face (j=my-1)
1412if (user->bctype[3]==1 || user->bctype[3]==-1) {
1413 if (ye==my) {
1414 j= ye-1;
1415 for (k=lzs; k<lze; k++) {
1416 for (i=lxs; i<lxe; i++) {
1417 ubcs[k][j][i].x = 0.0;
1418 ubcs[k][j][i].y = 0.0;
1419 ubcs[k][j][i].z = 0.0;
1420 // The relevant ucont is at the face, index j-1
1421 ucont[k][j-1][i].y = 0.0; // Enforce no-penetration
1422 }
1423 }
1424 }
1425}
1426
1427/* Original "Amir wall Ogrid" block can now be removed or commented out
1428 as its functionality is included above.
1429if (user->bctype[2]==1 || user->bctype[2]==-1) { ... }
1430*/
1431
1432/* ================================================================================== */
1433/* SYMMETRY BC */
1434/* ================================================================================== */
1435
1436 if (user->bctype[0]==3) {
1437 if (xs==0) {
1438 i= xs;
1439
1440 for (k=lzs; k<lze; k++) {
1441 for (j=lys; j<lye; j++) {
1442 A=sqrt(csi[k][j][i].z*csi[k][j][i].z +
1443 csi[k][j][i].y*csi[k][j][i].y +
1444 csi[k][j][i].x*csi[k][j][i].x);
1445 nx=csi[k][j][i].x/A;
1446 ny=csi[k][j][i].y/A;
1447 nz=csi[k][j][i].z/A;
1448 Un=ucat[k][j][i+1].x*nx+ucat[k][j][i+1].y*ny+ucat[k][j][i+1].z*nz;
1449 ubcs[k][j][i].x = ucat[k][j][i+1].x-Un*nx;//-V_frame.x;
1450 ubcs[k][j][i].y = ucat[k][j][i+1].y-Un*ny;
1451 ubcs[k][j][i].z = ucat[k][j][i+1].z-Un*nz;
1452 ucont[k][j][i].x = 0.;
1453 }
1454 }
1455 }
1456 }
1457
1458 if (user->bctype[1]==3) {
1459 if (xe==mx) {
1460 i= xe-1;
1461
1462 for (k=lzs; k<lze; k++) {
1463 for (j=lys; j<lye; j++) {
1464 A=sqrt(csi[k][j][i-1].z*csi[k][j][i-1].z +
1465 csi[k][j][i-1].y*csi[k][j][i-1].y +
1466 csi[k][j][i-1].x*csi[k][j][i-1].x);
1467 nx=csi[k][j][i-1].x/A;
1468 ny=csi[k][j][i-1].y/A;
1469 nz=csi[k][j][i-1].z/A;
1470 Un=ucat[k][j][i-1].x*nx+ucat[k][j][i-1].y*ny+ucat[k][j][i-1].z*nz;
1471 ubcs[k][j][i].x = ucat[k][j][i-1].x-Un*nx;
1472 ubcs[k][j][i].y = ucat[k][j][i-1].y-Un*ny;
1473 ubcs[k][j][i].z = ucat[k][j][i-1].z-Un*nz;
1474 ucont[k][j][i-1].x = 0.;
1475 }
1476 }
1477 }
1478 }
1479
1480 if (user->bctype[2]==3) {
1481 if (ys==0) {
1482 j= ys;
1483
1484 for (k=lzs; k<lze; k++) {
1485 for (i=lxs; i<lxe; i++) {
1486 A=sqrt(eta[k][j][i].z*eta[k][j][i].z +
1487 eta[k][j][i].y*eta[k][j][i].y +
1488 eta[k][j][i].x*eta[k][j][i].x);
1489 nx=eta[k][j][i].x/A;
1490 ny=eta[k][j][i].y/A;
1491 nz=eta[k][j][i].z/A;
1492 Un=ucat[k][j+1][i].x*nx+ucat[k][j+1][i].y*ny+ucat[k][j+1][i].z*nz;
1493 ubcs[k][j][i].x = ucat[k][j+1][i].x-Un*nx;
1494 ubcs[k][j][i].y = ucat[k][j+1][i].y-Un*ny;
1495 ubcs[k][j][i].z = ucat[k][j+1][i].z-Un*nz;
1496 ucont[k][j][i].y = 0.;
1497 }
1498 }
1499 }
1500 }
1501
1502 if (user->bctype[3]==3) {
1503 if (ye==my) {
1504 j=ye-1;
1505
1506 for (k=lzs; k<lze; k++) {
1507 for (i=lxs; i<lxe; i++) {
1508 A=sqrt(eta[k][j-1][i].z*eta[k][j-1][i].z +
1509 eta[k][j-1][i].y*eta[k][j-1][i].y +
1510 eta[k][j-1][i].x*eta[k][j-1][i].x);
1511 nx=eta[k][j-1][i].x/A;
1512 ny=eta[k][j-1][i].y/A;
1513 nz=eta[k][j-1][i].z/A;
1514 Un=ucat[k][j-1][i].x*nx+ucat[k][j-1][i].y*ny+ucat[k][j-1][i].z*nz;
1515 ubcs[k][j][i].x = ucat[k][j-1][i].x-Un*nx;
1516 ubcs[k][j][i].y = ucat[k][j-1][i].y-Un*ny;
1517 ubcs[k][j][i].z = ucat[k][j-1][i].z-Un*nz;
1518 ucont[k][j-1][i].y = 0.;
1519 }
1520 }
1521 }
1522 }
1523
1524
1525 if (user->bctype[4]==3) {
1526 if (zs==0) {
1527 k = 0;
1528 for (j=lys; j<lye; j++) {
1529 for (i=lxs; i<lxe; i++) {
1530 A=sqrt(zet[k][j][i].z*zet[k][j][i].z +
1531 zet[k][j][i].y*zet[k][j][i].y +
1532 zet[k][j][i].x*zet[k][j][i].x);
1533 nx=zet[k][j][i].x/A;
1534 ny=zet[k][j][i].y/A;
1535 nz=zet[k][j][i].z/A;
1536 Un=ucat[k+1][j][i].x*nx+ucat[k+1][j][i].y*ny+ucat[k+1][j][i].z*nz;
1537 ubcs[k][j][i].x = ucat[k+1][j][i].x-Un*nx;
1538 ubcs[k][j][i].y = ucat[k+1][j][i].y-Un*ny;
1539 ubcs[k][j][i].z = ucat[k+1][j][i].z-Un*nz;
1540 ucont[k][j][i].z = 0.;
1541 }
1542 }
1543 }
1544 }
1545
1546 if (user->bctype[5]==3) {
1547 if (ze==mz) {
1548 k = ze-1;
1549 for (j=lys; j<lye; j++) {
1550 for (i=lxs; i<lxe; i++) {
1551 A=sqrt(zet[k-1][j][i].z*zet[k-1][j][i].z +
1552 zet[k-1][j][i].y*zet[k-1][j][i].y +
1553 zet[k-1][j][i].x*zet[k-1][j][i].x);
1554 nx=zet[k-1][j][i].x/A;
1555 ny=zet[k-1][j][i].y/A;
1556 nz=zet[k-1][j][i].z/A;
1557 Un=ucat[k-1][j][i].x*nx+ucat[k-1][j][i].y*ny+ucat[k-1][j][i].z*nz;
1558 ubcs[k][j][i].x = ucat[k-1][j][i].x-Un*nx;
1559 ubcs[k][j][i].y = ucat[k-1][j][i].y-Un*ny;
1560 ubcs[k][j][i].z = ucat[k-1][j][i].z-Un*nz;
1561 ucont[k-1][j][i].z = 0.;
1562 }
1563 }
1564 }
1565 }
1566
1567/* ================================================================================== */
1568/* CHARACTERISTIC OUTLET BC :8 */
1569/* ================================================================================== */
1570
1571 if (user->bctype[5]==8) {
1572 if (ze == mz) {
1573 k = ze-2;
1574 FluxOut = 0;
1575 for (j=lys; j<lye; j++) {
1576 for (i=lxs; i<lxe; i++) {
1577 FluxOut += ucont[k][j][i].z;
1578 }
1579 }
1580 }
1581 else {
1582 FluxOut = 0.;
1583 }
1584
1585 FluxIn = simCtx->FluxInSum + FarFluxInSum;
1586
1587 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1588 // PetscGlobalSum(PETSC_COMM_WORLD,&FluxOut, &FluxOutSum);
1589
1590 //ratio = FluxInSum / FluxOutSum;
1591 ratio = FluxIn / simCtx->FluxOutSum;
1592 if (fabs(simCtx->FluxOutSum) < 1.e-6) ratio = 1.;
1593 //if (fabs(FluxInSum) <1.e-6) ratio = 0.;
1594 if (fabs(FluxIn) <1.e-6) ratio = 0.;
1595 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Char Ratio %d %le %le %le %le %d %d\n", simCtx->step, ratio, FluxIn, simCtx->FluxOutSum, FarFluxInSum,zs, ze);
1596
1597 if (ze==mz) {
1598 k = ze-1;
1599 for (j=lys; j<lye; j++) {
1600 for (i=lxs; i<lxe; i++) {
1601 ubcs[k][j][i].x = ucat[k-1][j][i].x;
1602 ubcs[k][j][i].y = ucat[k-1][j][i].y;
1603 if (simCtx->step==0 || simCtx->step==1)
1604 if (simCtx->inletprofile<0)
1605 ubcs[k][j][i].z = -1.;
1606 else if (user->bctype[4]==6)
1607 ubcs[k][j][i].z = 0.;
1608 else
1609 ubcs[k][j][i].z = 1.;//ubcs[0][j][i].z;//-1.;//1.;
1610
1611 else
1612 ucont[k-1][j][i].z = ucont[k-1][j][i].z*ratio;
1613
1614 ubcs[k][j][i].z = ucont[k-1][j][i].z / zet[k-1][j][i].z;
1615 }
1616 }
1617 }
1618 }
1619
1620
1621/* ================================================================================== */
1622/* OUTLETBC :4 */
1623/* ================================================================================== */
1624
1625
1626 if (user->bctype[5]==OUTLET || user->bctype[5]==14 || user->bctype[5]==20) {
1627 lArea=0.;
1628 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"+Zeta Outlet \n");
1629 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"FluxOutSum before FormBCS applied = %.6f \n",simCtx->FluxOutSum);
1630 if (ze == mz) {
1631 // k = ze-3;
1632 k=ze-1;
1633 FluxOut = 0;
1634 for (j=lys; j<lye; j++) {
1635 for (i=lxs; i<lxe; i++) {
1636
1637 if ((nvert[k-1][j][i])<0.1) {
1638 FluxOut += (ucat[k-1][j][i].x * (zet[k-1][j][i].x) +
1639 ucat[k-1][j][i].y * (zet[k-1][j][i].y) +
1640 ucat[k-1][j][i].z * (zet[k-1][j][i].z));
1641
1642 lArea += sqrt( (zet[k-1][j][i].x) * (zet[k-1][j][i].x) +
1643 (zet[k-1][j][i].y) * (zet[k-1][j][i].y) +
1644 (zet[k-1][j][i].z) * (zet[k-1][j][i].z));
1645
1646 }
1647 }
1648 }
1649 }
1650 else {
1651 FluxOut = 0.;
1652 }
1653
1654 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1655 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1656
1657 user->simCtx->AreaOutSum = AreaSum;
1658
1659 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"AreaOutSum = %.6f | FluxOutSum = %.6f \n",AreaSum,simCtx->FluxOutSum);
1660
1661 if (simCtx->block_number>1 && user->bctype[5]==14) {
1662 simCtx->FluxOutSum += user->FluxIntfcSum;
1663 // AreaSum += user->AreaIntfcSum;
1664 }
1665
1666 FluxIn = simCtx->FluxInSum + FarFluxInSum + user->FluxIntpSum;
1667 if (user->bctype[5]==20)
1668 ratio = (FluxIn / simCtx->FluxOutSum);
1669 else
1670 ratio = (FluxIn - simCtx->FluxOutSum) / AreaSum;
1671
1672 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Ratio for momentum correction = %.6f \n",ratio);
1673
1674 /* user->FluxOutSum += ratio*user->simCtx->AreaOutSum; */
1675 simCtx->FluxOutSum =0.0;
1676 FluxOut=0.0;
1677 if (ze==mz) {
1678 k = ze-1;
1679 for (j=lys; j<lye; j++) {
1680 for (i=lxs; i<lxe; i++) {
1681 if ((nvert[k-1][j][i])<0.1) {
1682
1683 ubcs[k][j][i].x = ucat[k-1][j][i].x;//+ratio;
1684 ubcs[k][j][i].y = ucat[k-1][j][i].y;
1685 ubcs[k][j][i].z = ucat[k-1][j][i].z;// + ratio;//*n_z;
1686
1687 // ucont[k-1][j][i].z = ubcs[k][j][i].z * zet[k-1][j][i].z;
1688 if (user->bctype[5]==20)
1689 ucont[k-1][j][i].z = (ubcs[k][j][i].x * (zet[k-1][j][i].x ) +
1690 ubcs[k][j][i].y * (zet[k-1][j][i].y ) +
1691 ubcs[k][j][i].z * (zet[k-1][j][i].z ))*ratio;
1692
1693 else{
1694 ucont[k-1][j][i].z = (ubcs[k][j][i].x * (zet[k-1][j][i].x ) +
1695 ubcs[k][j][i].y * (zet[k-1][j][i].y ) +
1696 ubcs[k][j][i].z * (zet[k-1][j][i].z ))
1697 + ratio * sqrt( (zet[k-1][j][i].x) * (zet[k-1][j][i].x) +
1698 (zet[k-1][j][i].y) * (zet[k-1][j][i].y) +
1699 (zet[k-1][j][i].z) * (zet[k-1][j][i].z));
1700
1701 FluxOut += ucont[k-1][j][i].z;
1702 }
1703 }//if
1704 }
1705 }
1706 }
1707
1708 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1709 LOG_ALLOW(GLOBAL,LOG_TRACE, "Timestep = %d | FluxInSum = %.6f | FlucOutSum = %.6f | FluxIntfcSum = %.6f | FluxIntpSum = %.6f \n", simCtx->step, simCtx->FluxInSum, simCtx->FluxOutSum, user->FluxIntfcSum,user->FluxIntpSum);
1710
1711 } else if (user->bctype[5]==2) {
1712 /* Designed for driven cavity problem (top(k=kmax) wall moving)
1713 u_x = 1 at k==kmax */
1714 if (ze==mz) {
1715 k = ze-1;
1716 for (j=lys; j<lye; j++) {
1717 for (i=lxs; i<lxe; i++) {
1718 ubcs[k][j][i].x = 0.;// - ucat[k-1][j][i].x;
1719 ubcs[k][j][i].y = 1;//sin(2*3.14*simCtx->step*simCtx->dt);//1.;//- ucat[k-1][j][i].y;
1720 ubcs[k][j][i].z = 0.;//- ucat[k-1][j][i].z;
1721 }
1722 }
1723 }
1724 }
1725
1726
1727 /* OUTLET at k==0 */
1728 if (user->bctype[4]==OUTLET) {
1729 lArea=0.;
1730 if (zs == 0) {
1731 k = zs;
1732 // k= zs + 1;
1733 FluxOut = 0;
1734 for (j=lys; j<lye; j++) {
1735 for (i=lxs; i<lxe; i++) {
1736
1737
1738 FluxOut += ucat[k+1][j][i].z * zet[k][j][i].z ;
1739
1740 lArea += zet[k][j][i].z;
1741
1742
1743
1744 }
1745 }
1746 }
1747 else {
1748 FluxOut = 0.;
1749 }
1750
1751 FluxIn = simCtx->FluxInSum + FarFluxInSum;
1752
1753 MPI_Allreduce(&FluxOut,&simCtx->FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1754 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1755
1756
1757 ratio = (simCtx->FluxInSum - simCtx->FluxOutSum) / AreaSum;
1758 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Ratio b %d %le %le %le %le %d %d\n", simCtx->step,ratio, simCtx->FluxInSum, simCtx->FluxOutSum, AreaSum,zs, ze);
1759
1760 if (zs==0) {
1761 k = 0;
1762 for (j=lys; j<lye; j++) {
1763 for (i=lxs; i<lxe; i++) {
1764 ubcs[k][j][i].x = ucat[k+1][j][i].x;
1765 ubcs[k][j][i].y = ucat[k+1][j][i].y;
1766 ubcs[k][j][i].z = ucat[k+1][j][i].z;
1767 ucont[k][j][i].z = (ubcs[k][j][i].z+ratio) * zet[k][j][i].z;
1768 }
1769 }
1770 }
1771 }
1772
1773
1774
1775/* ================================================================================== */
1776/* Ogrid :77 */
1777/* ================================================================================== */
1778 /*
1779 if (user->bctype[3]=77 && Ogrid)
1780 {Cmpnts ***coor;
1781 lArea=0.;
1782 FluxOut=0.0;
1783 // k = ze-3;
1784
1785 Vec Coor; DMGetCoordinatesLocal(da, &Coor);
1786 DMDAVecGetArray(fda,Coor,&coor);
1787 if (ye==my) {
1788 j=my-2;
1789 for (k=lzs; k<lze; k++){
1790 for (i=lxs; i<lxe; i++) {
1791
1792 FluxOut += ucont[k][j][i].y;
1793 lArea += sqrt( (eta[k-1][j][i].x) * (eta[k-1][j][i].x) +
1794 (eta[k-1][j][i].y) * (eta[k-1][j][i].y) +
1795 (eta[k-1][j][i].z) * (eta[k-1][j][i].z));
1796
1797 }
1798 }
1799 }
1800
1801 else {
1802 FluxOut = 0.;
1803 }
1804
1805 MPI_Allreduce(&FluxOut,&FluxOutSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1806 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1807
1808 user->FluxOutSum = FluxOutSum;
1809 simCtx->AreaOutSum = AreaSum;
1810 export PL="$HOME/CSBL/Codes/fdf_project/pic_les"
1811 ratio=2*(FluxIn-FluxOutSum)/AreaSum;
1812 ratio=0.0;
1813 if (ye==my){
1814 j=my-2;
1815 for (k=lzs; k<lze; k++){
1816 for (i=lxs; i<lxe; i++) {
1817 if ((nvert[k-1][j][i])<0.1 && coor[k][j][i].z >= 800.) {
1818 ucont[k][j][i].y *= (1+ratio);
1819
1820 }
1821 }
1822 }
1823 }
1824
1825 if (ye==my){
1826 j=my-2;
1827 for (k=lzs; k<lze; k++){
1828 for (i=lxs; i<lxe; i++) {
1829
1830 ubcs[k][j+1][i].z=ucat[k][j][i].z;
1831 ubcs[k][j+1][i].y=ucat[k][j][i].y;
1832 ubcs[k][j+1][i].x=ucat[k][j][i].x;
1833 }
1834 }
1835 }
1836
1837 // LOG_ALLOW(GLOBAL,LOG_DEBUG, " ratio %le ",ratio);
1838
1839 }
1840
1841 */
1842
1843
1844/* ================================================================================== */
1845/* Channelz */
1846/* ================================================================================== */
1847 // Amir channel flow correction
1848 if (user->bctype[4]==7 && simCtx->channelz==1) {
1849 Vec Coor; DMGetCoordinatesLocal(da, &Coor);
1850 DMDAVecGetArray(fda,Coor,&coor);
1851 Cmpnts ***uch;
1852 DMDAVecGetArray(fda, user->Bcs.Uch, &uch);
1853
1854 lArea=0.0;
1855 // if (zs==0) {
1856 // k=0;
1857 FluxIn=0.0;
1858
1859 double Fluxbcs=0.0, Fluxbcssum,ratiobcs;
1860
1861 if (zs==0) {
1862 k=0;
1863 Fluxbcs=0.0;
1864 for (j=lys; j<lye; j++) {
1865 for (i=lxs; i<lxe; i++) {
1866 if (nvert[k][j][i]<0.1){
1867 Fluxbcs += ucont[k][j][i].z;
1868
1869 lArea += sqrt((zet[k][j][i].x) * (zet[k][j][i].x) +
1870 (zet[k][j][i].y) * (zet[k][j][i].y) +
1871 (zet[k][j][i].z) * (zet[k][j][i].z));
1872 }
1873 }
1874 }
1875 }
1876
1877
1878 //int kk=(simCtx->step % (mz-2))+2;
1879
1880 for (k=zs;k<lze;k++){
1881 for (j=lys; j<lye; j++) {
1882 for (i=lxs; i<lxe; i++) {
1883 if (nvert[k][j][i]<0.1){
1884 FluxIn += ucont[k][j][i].z /((mz)-1);
1885 }
1886 }
1887 }
1888 }
1889 // else {
1890 // FluxIn=0.0;
1891// }
1892
1893 MPI_Allreduce(&FluxIn,&simCtx->Fluxsum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1894 MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1895 MPI_Allreduce(&Fluxbcs,&Fluxbcssum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1896
1897 if (simCtx->step==simCtx->StartStep && simCtx->StartStep > 0 && simCtx->ccc==0) {
1898 simCtx->ccc=1;
1899 simCtx->FluxInSum=Fluxbcssum;
1900// simCtx->FluxInSum=6.3908;
1901 LOG_ALLOW(LOCAL,LOG_DEBUG, " FluxInSum %le .\n ",simCtx->FluxInSum);
1902 }
1903
1904 simCtx->FluxInSum=6.35066;
1905 ratio=(simCtx->FluxInSum-simCtx->Fluxsum)/AreaSum;
1906 simCtx->ratio=ratio;
1907 ratiobcs=(simCtx->FluxInSum-Fluxbcssum)/AreaSum;
1908
1909 if (zs==0) {
1910 k=0;
1911 for (j=lys; j<lye; j++) {
1912 for (i=lxs; i<lxe; i++) {
1913 if (nvert[k+1][j][i]<0.1){
1914 if (simCtx->fish) {
1915 ucont[k][j][i].z+=ratiobcs * /* (1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt( (zet[k+1][j][i].x) * (zet[k+1][j][i].x) +
1916 (zet[k+1][j][i].y) * (zet[k+1][j][i].y) +
1917 (zet[k+1][j][i].z) * (zet[k+1][j][i].z));
1918 }
1919
1920 uch[k][j][i].z=ratiobcs * /* (1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt( (zet[k+1][j][i].x) * (zet[k+1][j][i].x) +
1921 (zet[k+1][j][i].y) * (zet[k+1][j][i].y) +
1922 (zet[k+1][j][i].z) * (zet[k+1][j][i].z));
1923 }
1924 }
1925 }
1926 }
1927
1928 if (ze==mz) {
1929 k=mz-1;
1930 for (j=lys; j<lye; j++) {
1931 for (i=lxs; i<lxe; i++) {
1932 if (nvert[k-1][j][i]<0.1){
1933 if (simCtx->fish){
1934 ucont[k-1][j][i].z+=ratiobcs * /*(1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt((zet[k-1][j][i].x) * (zet[k-1][j][i].x) +
1935 (zet[k-1][j][i].y) * (zet[k-1][j][i].y) +
1936 (zet[k-1][j][i].z) * (zet[k-1][j][i].z));
1937 }
1938 uch[k][j][i].z=ratiobcs * /*(1.-exp(-500. * (1.-fabs(coor[k][j][i].y)))) */ sqrt( (zet[k+1][j][i].x) * (zet[k+1][j][i].x) +
1939 (zet[k+1][j][i].y) * (zet[k+1][j][i].y) +
1940 (zet[k+1][j][i].z) * (zet[k+1][j][i].z));
1941 }
1942 }
1943 }
1944 }
1945 DMDAVecRestoreArray(fda, user->Bcs.Uch, &uch);
1946 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Ratio %le %.15le %.15le %.15le \n", ratio, ratiobcs, simCtx->FluxInSum,AreaSum);
1947 DMDAVecRestoreArray(fda,Coor,&coor);
1948
1949 ////.................////
1950
1951
1952
1953
1954 //just for check
1955/* if (zs==0) { */
1956/* k=10; */
1957/* FluxIn=0.0; */
1958/* for (j=lys; j<lye; j++) { */
1959/* for (i=lxs; i<lxe; i++) { */
1960/* if (nvert[k+1][j][i]<0.1){ */
1961/* FluxIn += ucont[k][j][i].z; */
1962/* lArea += sqrt((zet[k+1][j][i].x) * (zet[k+1][j][i].x) + */
1963/* (zet[k+1][j][i].y) * (zet[k+1][j][i].y) + */
1964/* (zet[k+1][j][i].z) * (zet[k+1][j][i].z)); */
1965/* } */
1966/* } */
1967/* } */
1968/* } */
1969/* else { */
1970/* FluxIn=0.0; */
1971/* } */
1972/* // LOG_ALLOW(PETSC_COMM_SELF, " Fluxsum %le ",FluxIn); */
1973
1974/* MPI_Allreduce(&FluxIn,&Fluxsum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); */
1975/* MPI_Allreduce(&lArea,&AreaSum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); */
1976
1977
1978
1979 }//channel
1980
1981
1982
1983 /* ==================================================================================== */
1984 /* Cylinder O-grid */
1985 /* ==================================================================================== */
1986 if (user->bctype[3]==12) {
1987 /* Designed to test O-grid for flow over a cylinder at jmax velocity is 1 (horizontal)
1988 u_x = 1 at k==kmax */
1989 if (ye==my) {
1990 j = ye-1;
1991 for (k=lzs; k<lze; k++) {
1992 for (i=lxs; i<lxe; i++) {
1993 ubcs[k][j][i].x = 0.;
1994 ubcs[k][j][i].y = 0.;
1995 ubcs[k][j][i].z = 1.;
1996 }
1997 }
1998 }
1999 }
2000 /* ==================================================================================== */
2001 /* Annulus */
2002 /* ==================================================================================== */
2003 /* designed to test periodic boundary condition for O-grid j=0 rotates */
2004 DMDAVecGetArray(fda, user->Cent, &cent);
2005 if (user->bctype[2]==11) {
2006 if (ys==0){
2007 j=0;
2008 for (k=lzs; k<lze; k++) {
2009 for (i=lxs; i<lxe; i++) {
2010
2011 /* ubcs[k][j][i].x=0.0; */
2012
2013/* ubcs[k][j][i].y = -cent[k][j+1][i].z/sqrt(cent[k][j+1][i].z*cent[k][j+1][i].z+cent[k][j+1][i].y*cent[k][j+1][i].y); */
2014
2015/* ubcs[k][j][i].z =cent[k][j+1][i].y/sqrt(cent[k][j+1][i].z*cent[k][j+1][i].z+cent[k][j+1][i].y*cent[k][j+1][i].y); */
2016 ubcs[k][j][i].x = cent[k][j+1][i].y/sqrt(cent[k][j+1][i].x*cent[k][j+1][i].x+cent[k][j+1][i].y*cent[k][j+1][i].y);;
2017 ubcs[k][j][i].y =-cent[k][j+1][i].x/sqrt(cent[k][j+1][i].x*cent[k][j+1][i].x+cent[k][j+1][i].y*cent[k][j+1][i].y);
2018 ubcs[k][j][i].z =0.0;
2019 // if(k==1) LOG_ALLOW(PETSC_COMM_SELF, "@ i= %d j=%d k=%d ubcs.y is %le\n",i,j,k,ubcs[k][j][i].y);
2020 }
2021 }
2022 }
2023 }
2024 /* ==================================================================================== */
2025 /* Rheology */
2026 /* ==================================================================================== */
2027
2028 if(simCtx->rheology && (user->bctype[2]==13 || user->bctype[3]==13 || user->bctype[4]==13 || user->bctype[5]==13)){
2029 LOG_ALLOW(GLOBAL,LOG_DEBUG, "moving plate velocity for rheology setup is %le \n",simCtx->U_bc);
2030 }
2031
2032 if (user->bctype[2]==13){
2033 if (ys==0){
2034 j=0;
2035 for (k=lzs; k<lze; k++) {
2036 for (i=lxs; i<lxe; i++) {
2037 ubcs[k][j][i].x = 0.;
2038 // ubcs[k][j][i].x = -simCtx->U_bc;
2039 ubcs[k][j][i].y = 0.;
2040 ubcs[k][j][i].z = -simCtx->U_bc;
2041 //ubcs[k][j][i].z =0.0;
2042 }
2043 }
2044 }
2045 }
2046 if (user->bctype[3]==13){
2047 if (ye==my){
2048 j=ye-1;
2049 for (k=lzs; k<lze; k++) {
2050 for (i=lxs; i<lxe; i++) {
2051 ubcs[k][j][i].x = 0.;
2052 // ubcs[k][j][i].x = simCtx->U_bc;
2053 ubcs[k][j][i].y = 0.;
2054 ubcs[k][j][i].z = simCtx->U_bc;
2055 //ubcs[k][j][i].z =0.0;
2056 }
2057 }
2058 }
2059 }
2060
2061 if (user->bctype[4]==13){
2062 if (zs==0){
2063 k=0;
2064 for (j=lys; j<lye; j++) {
2065 for (i=lxs; i<lxe; i++) {
2066 ubcs[k][j][i].x =-simCtx->U_bc;
2067 ubcs[k][j][i].y = 0.;
2068 ubcs[k][j][i].z = 0.;
2069 }
2070 }
2071 }
2072 }
2073 if (user->bctype[5]==13){
2074 if (ze==mz){
2075 k=ze-1;
2076 for (j=lys; j<lye; j++) {
2077 for (i=lxs; i<lxe; i++) {
2078 ubcs[k][j][i].x = simCtx->U_bc;
2079 ubcs[k][j][i].y = 0.;
2080 ubcs[k][j][i].z = 0.;
2081 }
2082 }
2083 }
2084 }
2085 DMDAVecRestoreArray(fda, user->Cent, &cent);
2086 DMDAVecRestoreArray(fda, user->lUcat, &ucat);
2087 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
2088 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2089 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2090 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2091
2092
2093
2094 Contra2Cart(user); // it also does global to local for Ucat
2095
2096/* ================================================================================== */
2097/* WALL FUNCTION */
2098/* ================================================================================== */
2099
2100 if (simCtx->wallfunction && user->bctype[2]==-1) {
2101 PetscReal ***ustar, ***aj;
2102 //Mohsen Dec 2015
2103 Vec Aj = user->lAj;
2104 DMDAVecGetArray(fda, user->Ucat, &ucat);
2105 DMDAVecGetArray(fda, user->Ucont, &ucont);
2106 DMDAVecGetArray(da, Aj, &aj);
2107// DMDAVecGetArray(da, user->Nvert, &nvert);
2108 DMDAVecGetArray(da, user->lUstar, &ustar);
2109
2110 // wall function for boundary
2111 for (k=lzs; k<lze; k++)
2112 for (j=lys; j<lye; j++)
2113 for (i=lxs; i<lxe; i++) {
2114
2115 if( nvert[k][j][i]<1.1 && user->bctype[2]==-1 && j==1 )
2116 {
2117 double area = sqrt( eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z );
2118 double sb, sc;
2119 double ni[3], nj[3], nk[3];
2120 Cmpnts Uc, Ua;
2121
2122 Ua.x = Ua.y = Ua.z = 0;
2123 sb = 0.5/aj[k][j][i]/area;
2124
2125
2126 sc = 2*sb + 0.5/aj[k][j+1][i]/area;
2127 Uc = ucat[k][j+1][i];
2128
2129
2130 //Calculate_normal(csi[k][j][i], eta[k][j][i], zet[k][j][i], ni, nj, nk);
2131 //if(j==my-2) nj[0]*=-1, nj[1]*=-1, nj[2]*=-1;
2132
2133 double AA=sqrt(eta[k][j][i].z*eta[k][j][i].z +
2134 eta[k][j][i].y*eta[k][j][i].y +
2135 eta[k][j][i].x*eta[k][j][i].x);
2136 nj[0]=eta[k][j][i].x/AA;
2137 nj[1]=eta[k][j][i].y/AA;
2138 nj[2]=eta[k][j][i].z/AA;
2139 noslip (user, sc, sb, Ua, Uc, &ucat[k][j][i], nj[0], nj[1], nj[2]);
2140 wall_function_loglaw(user, 1.e-16, sc, sb, Ua, Uc, &ucat[k][j][i], &ustar[k][j][i], nj[0], nj[1], nj[2]);
2141
2142 // nvert[k][j][i]=1.; /* set nvert to 1 to exclude from rhs */
2143 // if (k==1)
2144 // LOG_ALLOW(GLOBAL,LOG_DEBUG, " %d %le %le %le %le %le %le %le %le %le\n",i, sb,aj[k][j][i],AA, nj[0], nj[1], nj[2],ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z);
2145
2146 }
2147 }
2148 if (ys==0) {
2149 j= ys;
2150
2151 for (k=lzs; k<lze; k++) {
2152 for (i=lxs; i<lxe; i++) {
2153 ubcs[k][j][i].x = 0.0;
2154 ubcs[k][j][i].y = 0.0;
2155 ubcs[k][j][i].z = 0.0;
2156 ucont[k][j][i].y = 0.;
2157 }
2158 }
2159 }
2160
2161 DMDAVecRestoreArray(da, Aj, &aj);
2162 DMDAVecRestoreArray(da, user->lUstar, &ustar);
2163 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2164 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
2165 // DMDAVecRestoreArray(da, user->Nvert, &nvert);
2166 // DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2167 // DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2168
2169 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2170 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2171// DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat);
2172// DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat);
2173
2174
2175 }
2176
2177/* ================================================================================== */
2178
2179 DMDAVecGetArray(fda, user->Ucat, &ucat);
2180
2181/* ================================================================================== */
2182
2183 // boundary conditions on ghost nodes
2184 if (xs==0 && user->bctype[0]!=7) {
2185 i = xs;
2186 for (k=lzs; k<lze; k++) {
2187 for (j=lys; j<lye; j++) {
2188 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i+1].x;
2189 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i+1].y;
2190 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i+1].z;
2191 }
2192 }
2193 }
2194
2195 if (xe==mx && user->bctype[0]!=7) {
2196 i = xe-1;
2197 for (k=lzs; k<lze; k++) {
2198 for (j=lys; j<lye; j++) {
2199 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i-1].x;
2200 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i-1].y;
2201 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i-1].z;
2202 }
2203 }
2204 }
2205
2206
2207 if (ys==0 && user->bctype[2]!=7) {
2208 j = ys;
2209 for (k=lzs; k<lze; k++) {
2210 for (i=lxs; i<lxe; i++) {
2211 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j+1][i].x;
2212 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j+1][i].y;
2213 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j+1][i].z;
2214 }
2215 }
2216 }
2217
2218 if (ye==my && user->bctype[2]!=7) {
2219 j = ye-1;
2220 for (k=lzs; k<lze; k++) {
2221 for (i=lxs; i<lxe; i++) {
2222 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j-1][i].x;
2223 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j-1][i].y;
2224 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j-1][i].z;
2225 }
2226 }
2227 }
2228
2229 if (zs==0 && user->bctype[4]!=7) {
2230 k = zs;
2231 for (j=lys; j<lye; j++) {
2232 for (i=lxs; i<lxe; i++) {
2233 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k+1][j][i].x;
2234 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k+1][j][i].y;
2235 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k+1][j][i].z;
2236 }
2237 }
2238 }
2239
2240 if (ze==mz && user->bctype[4]!=7) {
2241 k = ze-1;
2242 for (j=lys; j<lye; j++) {
2243 for (i=lxs; i<lxe; i++) {
2244 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k-1][j][i].x;
2245 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k-1][j][i].y;
2246 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k-1][j][i].z;
2247 }
2248 }
2249 }
2250
2251 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2252 /* ================================================================================== */
2253 /* Periodic BC *///Mohsen
2254/* ================================================================================== */
2255 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7){
2256 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2257 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2258 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2259 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2260 /* /\* DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert); *\/ */
2261 /* /\* DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert); *\/ */
2262 //Mohsen Dec 2015
2263 DMDAVecGetArray(da, user->lP, &lp);
2264 DMDAVecGetArray(da, user->lNvert, &lnvert);
2265 DMDAVecGetArray(fda, user->lUcat, &lucat);
2266 DMDAVecGetArray(da, user->P, &p);
2267 DMDAVecGetArray(da, user->Nvert, &nvert);
2268 DMDAVecGetArray(fda, user->Ucat, &ucat);
2269
2270 if (user->bctype[0]==7 || user->bctype[1]==7){
2271 if (xs==0){
2272 i=xs;
2273 for (k=zs; k<ze; k++) {
2274 for (j=ys; j<ye; j++) {
2275 if(j>0 && k>0 && j<user->JM && k<user->KM){
2276 ucat[k][j][i]=lucat[k][j][i-2];
2277 p[k][j][i]=lp[k][j][i-2];
2278 nvert[k][j][i]=lnvert[k][j][i-2];
2279 }
2280 }
2281 }
2282 }
2283 }
2284 if (user->bctype[2]==7 || user->bctype[3]==7){
2285 if (ys==0){
2286 j=ys;
2287 for (k=zs; k<ze; k++) {
2288 for (i=xs; i<xe; i++) {
2289 if(i>0 && k>0 && i<user->IM && k<user->KM){
2290 ucat[k][j][i]=lucat[k][j-2][i];
2291 p[k][j][i]=lp[k][j-2][i];
2292 nvert[k][j][i]=lnvert[k][j-2][i];
2293 }
2294 }
2295 }
2296 }
2297 }
2298 if (user->bctype[4]==7 || user->bctype[5]==7){
2299 if (zs==0){
2300 k=zs;
2301 for (j=ys; j<ye; j++) {
2302 for (i=xs; i<xe; i++) {
2303 if(i>0 && j>0 && i<user->IM && j<user->JM){
2304 ucat[k][j][i]=lucat[k-2][j][i];
2305 nvert[k][j][i]=lnvert[k-2][j][i];
2306 //amir
2307 p[k][j][i]=lp[k-2][j][i];
2308 }
2309 }
2310 }
2311 }
2312 }
2313 if (user->bctype[0]==7 || user->bctype[1]==7){
2314 if (xe==mx){
2315 i=mx-1;
2316 for (k=zs; k<ze; k++) {
2317 for (j=ys; j<ye; j++) {
2318 if(j>0 && k>0 && j<user->JM && k<user->KM){
2319 ucat[k][j][i]=lucat[k][j][i+2];
2320 p[k][j][i]=lp[k][j][i+2];
2321 nvert[k][j][i]=lnvert[k][j][i+2];
2322 }
2323 }
2324 }
2325 }
2326 }
2327 if (user->bctype[2]==7 || user->bctype[3]==7){
2328 if (ye==my){
2329 j=my-1;
2330 for (k=zs; k<ze; k++) {
2331 for (i=xs; i<xe; i++) {
2332 if(i>0 && k>0 && i<user->IM && k<user->KM){
2333 ucat[k][j][i]=lucat[k][j+2][i];
2334 p[k][j][i]=lp[k][j+2][i];
2335 nvert[k][j][i]=lnvert[k][j+2][i];
2336 }
2337 }
2338 }
2339 }
2340 }
2341
2342 if (user->bctype[4]==7 || user->bctype[5]==7){
2343 if (ze==mz){
2344 k=mz-1;
2345 for (j=ys; j<ye; j++) {
2346 for (i=xs; i<xe; i++) {
2347 if(i>0 && j>0 && i<user->IM && j<user->JM){
2348 ucat[k][j][i]=lucat[k+2][j][i];
2349 nvert[k][j][i]=lnvert[k+2][j][i];
2350 //amir
2351 p[k][j][i]=lp[k+2][j][i];
2352 }
2353 }
2354 }
2355 }
2356 }
2357
2358
2359
2360
2361
2362 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2363 DMDAVecRestoreArray(da, user->lP, &lp);
2364 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2365 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2366 DMDAVecRestoreArray(da, user->P, &p);
2367 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2368
2369 /* /\* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); *\/ */
2370 /* /\* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); *\/ */
2371
2372 /* /\* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); *\/ */
2373 /* /\* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); *\/ */
2374
2375 /* /\* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); *\/ */
2376 /* /\* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); *\/ */
2377
2378 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2379 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2380
2381 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2382 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2383
2384 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2385 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2386}
2387 // 0 velocity on the corner point
2388
2389 DMDAVecGetArray(fda, user->Ucat, &ucat);
2390 DMDAVecGetArray(da, user->P, &p);
2391
2392 if (zs==0) {
2393 k=0;
2394 if (xs==0) {
2395 i=0;
2396 for (j=ys; j<ye; j++) {
2397 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i+1].x);
2398 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i+1].y);
2399 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i+1].z);
2400 p[k][j][i]= 0.5*(p[k+1][j][i]+p[k][j][i+1]);
2401 }
2402 }
2403 if (xe == mx) {
2404 i=mx-1;
2405 for (j=ys; j<ye; j++) {
2406 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i-1].x);
2407 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i-1].y);
2408 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i-1].z);
2409 p[k][j][i] = 0.5*(p[k+1][j][i]+p[k][j][i-1]);
2410 }
2411 }
2412
2413 if (ys==0) {
2414 j=0;
2415 for (i=xs; i<xe; i++) {
2416 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j+1][i].x);
2417 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j+1][i].y);
2418 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j+1][i].z);
2419 p[k][j][i] = 0.5*(p[k+1][j][i]+p[k][j+1][i]);
2420 }
2421 }
2422
2423 if (ye==my) {
2424 j=my-1;
2425 for (i=xs; i<xe; i++) {
2426 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j-1][i].x);
2427 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j-1][i].y);
2428 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j-1][i].z);
2429 p[k][j][i] = 0.5*(p[k+1][j][i]+p[k][j-1][i]);
2430 }
2431 }
2432 }
2433
2434 if (ze==mz) {
2435 k=mz-1;
2436 if (xs==0) {
2437 i=0;
2438 for (j=ys; j<ye; j++) {
2439 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i+1].x);
2440 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i+1].y);
2441 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i+1].z);
2442 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j][i+1]);
2443 }
2444 }
2445 if (xe == mx) {
2446 i=mx-1;
2447 for (j=ys; j<ye; j++) {
2448 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i-1].x);
2449 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i-1].y);
2450 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i-1].z);
2451 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j][i-1]);
2452 }
2453 }
2454
2455 if (ys==0) {
2456 j=0;
2457 for (i=xs; i<xe; i++) {
2458 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j+1][i].x);
2459 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j+1][i].y);
2460 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j+1][i].z);
2461 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j+1][i]);
2462 }
2463 }
2464
2465 if (ye==my) {
2466 j=my-1;
2467 for (i=xs; i<xe; i++) {
2468 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j-1][i].x);
2469 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j-1][i].y);
2470 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j-1][i].z);
2471 p[k][j][i] = 0.5*(p[k-1][j][i]+p[k][j-1][i]);
2472 }
2473 }
2474 }
2475
2476 if (ys==0) {
2477 j=0;
2478 if (xs==0) {
2479 i=0;
2480 for (k=zs; k<ze; k++) {
2481 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i+1].x);
2482 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i+1].y);
2483 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i+1].z);
2484 p[k][j][i]= 0.5*(p[k][j+1][i]+p[k][j][i+1]);
2485 }
2486 }
2487
2488 if (xe==mx) {
2489 i=mx-1;
2490 for (k=zs; k<ze; k++) {
2491 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i-1].x);
2492 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i-1].y);
2493 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i-1].z);
2494 p[k][j][i] = 0.5*(p[k][j+1][i]+p[k][j][i-1]);
2495 }
2496 }
2497 }
2498
2499 if (ye==my) {
2500 j=my-1;
2501 if (xs==0) {
2502 i=0;
2503 for (k=zs; k<ze; k++) {
2504 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i+1].x);
2505 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i+1].y);
2506 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i+1].z);
2507 p[k][j][i] = 0.5*(p[k][j-1][i]+p[k][j][i+1]);
2508 }
2509 }
2510
2511 if (xe==mx) {
2512 i=mx-1;
2513 for (k=zs; k<ze; k++) {
2514 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i-1].x);
2515 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i-1].y);
2516 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i-1].z);
2517 p[k][j][i] = 0.5*(p[k][j-1][i]+p[k][j][i-1]);
2518 }
2519 }
2520 }
2521
2522 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2523 DMDAVecRestoreArray(da, user->P, &p);
2524
2525 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2526 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2527
2528 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2529 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2530
2531 //Mohsen Nov 2012
2532 //Velocity and Presurre at corners for Periodic BC's
2533
2534 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7){
2535 //i-direction
2536
2537 DMDAVecGetArray(fda, user->lUcat, &lucat);
2538 DMDAVecGetArray(da, user->lP, &lp);
2539 DMDAVecGetArray(da, user->lNvert, &lnvert);
2540 DMDAVecGetArray(fda, user->Ucat, &ucat);
2541 DMDAVecGetArray(da, user->P, &p);
2542 DMDAVecGetArray(da, user->Nvert, &nvert);
2543
2544 if (user->bctype[0]==7){
2545 if (xs==0){
2546 i=xs;
2547 for (k=zs; k<ze; k++) {
2548 for (j=ys; j<ye; j++) {
2549 ucat[k][j][i]=lucat[k][j][i-2];
2550 p[k][j][i]=lp[k][j][i-2];
2551 nvert[k][j][i]=lnvert[k][j][i-2];
2552 }
2553 }
2554 }
2555 }
2556 if (user->bctype[1]==7){
2557 if (xe==mx){
2558 i=xe-1;
2559 for (k=zs; k<ze; k++) {
2560 for (j=ys; j<ye; j++) {
2561 ucat[k][j][i].x=lucat[k][j][i+2].x;
2562 p[k][j][i]=lp[k][j][i+2];
2563 nvert[k][j][i]=lnvert[k][j][i+2];
2564 }
2565 }
2566 }
2567 }
2568 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2569 DMDAVecRestoreArray(da, user->lP, &lp);
2570 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2571 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2572 DMDAVecRestoreArray(da, user->P, &p);
2573 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2574
2575 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2576/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2577
2578/* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); */
2579/* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); */
2580
2581/* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2582/* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2583
2584 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2585 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2586
2587 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2588 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2589
2590 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2591 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2592
2593 //j-direction
2594 DMDAVecGetArray(fda, user->lUcat, &lucat);
2595 DMDAVecGetArray(da, user->lP, &lp);
2596 DMDAVecGetArray(da, user->lNvert, &lnvert);
2597 DMDAVecGetArray(fda, user->Ucat, &ucat);
2598 DMDAVecGetArray(da, user->P, &p);
2599 DMDAVecGetArray(da, user->Nvert, &nvert);
2600
2601 if (user->bctype[2]==7){
2602 if (ys==0){
2603 j=ys;
2604 for (k=zs; k<ze; k++) {
2605 for (i=xs; i<xe; i++) {
2606 ucat[k][j][i]=lucat[k][j-2][i];
2607 p[k][j][i]=lp[k][j-2][i];
2608 nvert[k][j][i]=lnvert[k][j-2][i];
2609 }
2610 }
2611 }
2612 }
2613
2614 if (user->bctype[3]==7){
2615 if (ye==my){
2616 j=my-1;
2617 for (k=zs; k<ze; k++) {
2618 for (i=xs; i<xe; i++) {
2619 ucat[k][j][i]=lucat[k][j+2][i];
2620 p[k][j][i]=lp[k][j+2][i];
2621 nvert[k][j][i]=lnvert[k][j+2][i];
2622 }
2623 }
2624 }
2625 }
2626
2627 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2628 DMDAVecRestoreArray(da, user->lP, &lp);
2629 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2630 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2631 DMDAVecRestoreArray(da, user->P, &p);
2632 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2633
2634/* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2635/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2636
2637/* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); */
2638/* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); */
2639
2640/* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2641/* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2642
2643 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2644 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2645
2646 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2647 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2648
2649 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2650 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2651
2652 //k-direction
2653 DMDAVecGetArray(fda, user->lUcat, &lucat);
2654 DMDAVecGetArray(da, user->lP, &lp);
2655 DMDAVecGetArray(da, user->lNvert, &lnvert);
2656 DMDAVecGetArray(fda, user->Ucat, &ucat);
2657 DMDAVecGetArray(da, user->P, &p);
2658 DMDAVecGetArray(da, user->Nvert, &nvert);
2659
2660 if (user->bctype[4]==7){
2661 if (zs==0){
2662 k=zs;
2663 for (j=ys; j<ye; j++) {
2664 for (i=xs; i<xe; i++) {
2665 ucat[k][j][i]=lucat[k-2][j][i];
2666 nvert[k][j][i]=lnvert[k-2][j][i];
2667 //amir
2668 p[k][j][i]=lp[k-2][j][i];
2669 }
2670 }
2671 }
2672 }
2673 if (user->bctype[5]==7){
2674 if (ze==mz){
2675 k=mz-1;
2676 for (j=ys; j<ye; j++) {
2677 for (i=xs; i<xe; i++) {
2678 ucat[k][j][i]=lucat[k+2][j][i];
2679 nvert[k][j][i]=lnvert[k+2][j][i];
2680 p[k][j][i]=lp[k+2][j][i];
2681 }
2682 }
2683 }
2684 }
2685
2686 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
2687 DMDAVecRestoreArray(da, user->lP, &lp);
2688 DMDAVecRestoreArray(da, user->lNvert, &lnvert);
2689 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2690 DMDAVecRestoreArray(da, user->P, &p);
2691 DMDAVecRestoreArray(da, user->Nvert, &nvert);
2692
2693/* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2694/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
2695
2696/* DMLocalToGlobalBegin(da, user->lP, INSERT_VALUES, user->P); */
2697/* DMLocalToGlobalEnd(da, user->lP, INSERT_VALUES, user->P); */
2698
2699/* DMLocalToGlobalBegin(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2700/* DMLocalToGlobalEnd(da, user->lNvert, INSERT_VALUES, user->Nvert); */
2701
2702 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2703 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2704
2705 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2706 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2707
2708 DMGlobalToLocalBegin(da, user->Nvert, INSERT_VALUES, user->lNvert);
2709 DMGlobalToLocalEnd(da, user->Nvert, INSERT_VALUES, user->lNvert);
2710 }
2711 /* ================================================================================== */
2712/* Analytical Vortex BC */
2713/* ================================================================================== */
2714
2715 DMDAVecGetArray(fda, user->Ucat, &ucat);
2716 DMDAVecGetArray(fda, user->Ucont, &ucont);
2717 DMDAVecGetArray(fda, user->Cent, &cent);
2718 DMDAVecGetArray(fda, user->Centx, &centx);
2719 DMDAVecGetArray(fda, user->Centy, &centy);
2720 DMDAVecGetArray(fda, user->Centz, &centz);
2721
2722 if (user->bctype[0]==9) {
2723 if (xs==0) {
2724 i= xs;
2725 for (k=lzs; k<lze; k++) {
2726 for (j=lys; j<lye; j++) {
2727 ucat[k][j][i].x=-cos(cent[k][j][i+1].x)*sin(cent[k][j][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2728 ucat[k][j][i].y=-sin(cent[k][j][i+1].x)*cos(cent[k][j][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2729 ucat[k][j][i].z =0.0;
2730
2731 ucont[k][j][i].x =-(cos(centx[k][j][i].x)*sin(centx[k][j][i].y)*csi[k][j][i].x)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2732 }
2733 }
2734 if (ys==0) {
2735 j=ys;
2736 for (k=lzs; k<lze; k++) {
2737 ucat[k][j][i].x=cos(cent[k][j+1][i+1].x)*sin(cent[k][j+1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2738 ucat[k][j][i].y=-sin(cent[k][j+1][i+1].x)*cos(cent[k][j+1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2739 ucat[k][j][i].z =0.0;
2740 }
2741 }
2742 if (ye==my) {
2743 j=ye-1;
2744 for (k=lzs; k<lze; k++) {
2745 ucat[k][j][i].x=cos(cent[k][j-1][i+1].x)*sin(cent[k][j-1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2746 ucat[k][j][i].y=-sin(cent[k][j-1][i+1].x)*cos(cent[k][j-1][i+1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2747 ucat[k][j][i].z =0.0;
2748 }
2749 }
2750 }
2751 }
2752 if (user->bctype[1]==9) {
2753 if (xe==mx) {
2754 i= xe-1;
2755 for (k=lzs; k<lze; k++) {
2756 for (j=lys; j<lye; j++) {
2757 ucat[k][j][i].x=-cos(cent[k][j][i-1].x)*sin(cent[k][j][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2758 ucat[k][j][i].y=-sin(cent[k][j][i-1].x)*cos(cent[k][j][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2759 ucat[k][j][i].z =0.0;
2760
2761 ucont[k][j][i-1].x =(-cos(centx[k][j][i-1].x)*sin(centx[k][j][i-1].y)*csi[k][j][i-1].x)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2762 }
2763 }
2764 if (ys==0) {
2765 j=ys;
2766 for (k=lzs; k<lze; k++) {
2767 ucat[k][j][i].x=cos(cent[k][j+1][i-1].x)*sin(cent[k][j+1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2768 ucat[k][j][i].y=-sin(cent[k][j+1][i-1].x)*cos(cent[k][j+1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2769 ucat[k][j][i].z =0.0;
2770 }
2771 }
2772 if (ye==my) {
2773 j=ye-1;
2774 for (k=lzs; k<lze; k++) {
2775 ucat[k][j][i].x=cos(cent[k][j-1][i-1].x)*sin(cent[k][j-1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2776 ucat[k][j][i].y=-sin(cent[k][j-1][i-1].x)*cos(cent[k][j-1][i-1].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2777 ucat[k][j][i].z =0.0;
2778 }
2779 }
2780 }
2781 }
2782
2783 if (user->bctype[2]==9) {
2784 if (ys==0) {
2785 j= ys;
2786 for (k=lzs; k<lze; k++) {
2787 for (i=lxs; i<lxe; i++) {
2788 ucat[k][j][i].x=cos(cent[k][j+1][i].x)*sin(cent[k][j+1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2789 ucat[k][j][i].y=sin(cent[k][j+1][i].x)*cos(cent[k][j+1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2790 ucat[k][j][i].z =0.0;
2791
2792 ucont[k][j][i].y=(sin(centy[k][j][i].x)*cos(centy[k][j][i].y)*eta[k][j][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2793 }
2794 }
2795 }
2796 }
2797
2798
2799 if (user->bctype[3]==9) {
2800 if (ye==my) {
2801 j= ye-1;
2802 for (k=lzs; k<lze; k++) {
2803 for (i=lxs; i<lxe; i++) {
2804 ucat[k][j][i].x=cos(cent[k][j-1][i].x)*sin(cent[k][j-1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2805 ucat[k][j][i].y=sin(cent[k][j-1][i].x)*cos(cent[k][j-1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2806 ucat[k][j][i].z =0.0;
2807
2808 ucont[k][j-1][i].y=(sin(centy[k][j-1][i].x)*cos(centy[k][j-1][i].y)*eta[k][j-1][i].y)*exp(-2.0*simCtx->dt*(simCtx->step+1)/simCtx->ren);
2809 }
2810 }
2811 }
2812 }
2813 if (user->bctype[4]==9) {
2814 if (zs==0) {
2815 k= zs;
2816 for (j=ys; j<ye; j++) {
2817 for (i=xs; i<xe; i++) {
2818 ucat[k][j][i].x=ucat[k+1][j][i].x;
2819 ucat[k][j][i].y=ucat[k+1][j][i].y;
2820 ucat[k][j][i].z=ucat[k+1][j][i].z;
2821
2822 ucont[k][j][i].z=0.0;
2823 }
2824 }
2825 }
2826 }
2827 if (user->bctype[5]==9) {
2828 if (ze==mz) {
2829 k= ze-1;
2830 for (j=ys; j<ye; j++) {
2831 for (i=xs; i<xe; i++) {
2832 ucat[k][j][i].x=ucat[k-1][j][i].x;
2833 ucat[k][j][i].y=ucat[k-1][j][i].y;
2834 ucat[k][j][i].z=ucat[k-1][j][i].z;
2835
2836 ucont[k-1][j][i].z=0.0;
2837 }
2838 }
2839 }
2840 }
2841
2842 DMDAVecRestoreArray(fda, user->Cent, &cent);
2843 DMDAVecRestoreArray(fda, user->Centx, &centx);
2844 DMDAVecRestoreArray(fda, user->Centy, &centy);
2845 DMDAVecRestoreArray(fda, user->Centz, &centz);
2846
2847 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
2848 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
2849
2850 } // ttemp
2851
2852
2853 DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs);
2854 DMDAVecRestoreArray(fda, user->lCsi, &csi);
2855 DMDAVecRestoreArray(fda, user->lEta, &eta);
2856 DMDAVecRestoreArray(fda, user->lZet, &zet);
2857
2858 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
2859 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
2860 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2861 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
2862 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2863 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2864
2865 PROFILE_FUNCTION_END;
2866
2867 PetscFunctionReturn(0);
2868 }
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2084
SimCtx::fish
PetscInt fish
Definition variables.h:572
SimCtx::block_number
PetscInt block_number
Definition variables.h:593
SimCtx::channelz
PetscInt channelz
Definition variables.h:593
UserCtx::Centz
Vec Centz
Definition variables.h:706
SimCtx::ren
PetscReal ren
Definition variables.h:581
SimCtx::dt
PetscReal dt
Definition variables.h:552
SimCtx::inletprofile
PetscInt inletprofile
Definition variables.h:593
UserCtx::lUstar
Vec lUstar
Definition variables.h:683
SimCtx::ccc
PetscInt ccc
Definition variables.h:605
SimCtx::ratio
PetscReal ratio
Definition variables.h:606
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:548
UserCtx::Centx
Vec Centx
Definition variables.h:706
UserCtx::FluxIntfcSum
PetscReal FluxIntfcSum
Definition variables.h:685
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:610
SimCtx::rheology
PetscInt rheology
Definition variables.h:567
SimCtx::Fluxsum
PetscReal Fluxsum
Definition variables.h:602
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:604
SimCtx::step
PetscInt step
Definition variables.h:546
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:603
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
UserCtx::FluxIntpSum
PetscReal FluxIntpSum
Definition variables.h:685
UserCtx::Nvert
Vec Nvert
Definition variables.h:688
UserCtx::Centy
Vec Centy
Definition variables.h:706
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:122
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 673 of file Boundaries.c.

674{
675 PetscErrorCode ierr;
676 PetscFunctionBeginUser;
677 PROFILE_FUNCTION_BEGIN;
678 // The sole purpose of this function is to call the old logic for the
679 // specific block context that was passed in.
680 ierr = FormBCS(user); CHKERRQ(ierr);
681
682 // NOTE: The legacy `main` called Block_Interface_U after the FormBCS loop.
683 // We will handle that at a higher level (in our AdvanceSimulation loop)
684 // after all blocks have had their BCs applied for the step.
685 PROFILE_FUNCTION_END;
686 PetscFunctionReturn(0);
687}
FormBCS
PetscErrorCode FormBCS(UserCtx *user)
Definition Boundaries.c:1022
Here is the call graph for this function: