PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
Boundaries.c File Reference
#include "Boundaries.h"
#include <string.h>
#include <ctype.h>
Include dependency graph for Boundaries.c:

Go to the source code of this file.

Macros

#define __FUNCT__   "CanRankServiceInletFace"
 
#define __FUNCT__   "CanRankServiceFace"
 
#define __FUNCT__   "GetDeterministicFaceGridLocation"
 
#define __FUNCT__   "GetRandomFCellAndLogicOnInletFace"
 
#define __FUNCT__   "ClassifyMomentumRow"
 
#define __FUNCT__   "EnforceRHSBoundaryConditions"
 
#define __FUNCT__   "BoundaryCondition_Create"
 
#define __FUNCT__   "BoundarySystem_Validate"
 
#define __FUNCT__   "BoundarySystem_Initialize"
 
#define __FUNCT__   "PropagateBoundaryConfigToCoarserLevels"
 
#define __FUNCT__   "BoundarySystem_ExecuteStep"
 
#define __FUNCT__   "BoundarySystem_RefreshUbcs"
 
#define __FUNCT__   "BoundarySystem_Destroy"
 
#define __FUNCT__   "TransferPeriodicFieldByDirection"
 
#define __FUNCT__   "SynchronizePeriodicCellFields"
 
#define __FUNCT__   "GetPersistentFaceField"
 
#define __FUNCT__   "TransferPeriodicFaceFieldByDirection"
 
#define __FUNCT__   "SynchronizePeriodicFaceFields"
 
#define __FUNCT__   "GetPersistentStaggeredField"
 
#define __FUNCT__   "TransferPeriodicStaggeredFieldByDirection"
 
#define __FUNCT__   "SynchronizePeriodicStaggeredFields"
 
#define __FUNCT__   "PreparePeriodicQuickStencilFields"
 
#define __FUNCT__   "SynchronizePeriodicLocalStaggeredField"
 
#define __FUNCT__   "ApplyMetricsPeriodicBCs"
 
#define __FUNCT__   "ApplyPeriodicBCs"
 
#define __FUNCT__   "UpdateDummyCells"
 
#define __FUNCT__   "UpdateCornerNodes"
 
#define __FUNCT__   "ApplyWallFunction"
 
#define __FUNCT__   "FinalizePostProjectionCellFields"
 
#define __FUNCT__   "ApplyBoundaryConditions"
 

Functions

PetscErrorCode CanRankServiceInletFace (UserCtx *user, const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, PetscBool *can_service_inlet_out)
 Internal helper implementation: CanRankServiceInletFace().
 
PetscErrorCode CanRankServiceFace (const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, BCFace face_id, PetscBool *can_service_out)
 Implementation of CanRankServiceFace().
 
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)
 Internal helper implementation: GetDeterministicFaceGridLocation().
 
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)
 Internal helper implementation: GetRandomCellAndLogicalCoordsOnInletFace().
 
MomentumRowType ClassifyMomentumRow (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscInt component, PetscInt *ri, PetscInt *rj, PetscInt *rk)
 Implementation of ClassifyMomentumRow().
 
PetscErrorCode EnforceRHSBoundaryConditions (UserCtx *user)
 Implementation of EnforceRHSBoundaryConditions().
 
PetscErrorCode BoundaryCondition_Create (BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
 Internal helper implementation: BoundaryCondition_Create().
 
PetscErrorCode BoundarySystem_Validate (UserCtx *user)
 Internal helper implementation: BoundarySystem_Validate().
 
PetscErrorCode BoundarySystem_Initialize (UserCtx *user, const char *bcs_filename)
 Implementation of BoundarySystem_Initialize().
 
PetscErrorCode PropagateBoundaryConfigToCoarserLevels (SimCtx *simCtx)
 Internal helper implementation: PropagateBoundaryConfigToCoarserLevels().
 
PetscErrorCode BoundarySystem_ExecuteStep (UserCtx *user)
 Implementation of BoundarySystem_ExecuteStep().
 
PetscErrorCode BoundarySystem_RefreshUbcs (UserCtx *user)
 Internal helper implementation: BoundarySystem_RefreshUbcs().
 
PetscErrorCode BoundarySystem_Destroy (UserCtx *user)
 Implementation of BoundarySystem_Destroy().
 
static PetscErrorCode TransferPeriodicFieldByDirection (UserCtx *user, FieldId field_id, char direction)
 Copies one cell field's wrapped local values onto the owned periodic duplicate plane.
 
PetscErrorCode SynchronizePeriodicCellFields (UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
 Implementation of SynchronizePeriodicCellFields().
 
static PetscErrorCode GetPersistentFaceField (UserCtx *user, FieldId field_id, char face_direction, DM *dm, Vec *global_vec, Vec *local_vec, PetscInt *dof)
 Resolves one registered persistent single-face-family field.
 
static PetscErrorCode IsFaceCenterCoordinateField (FieldId field_id, PetscBool *is_coordinate)
 Returns whether a registered face field stores physical coordinates.
 
static PetscErrorCode TranslatePeriodicFaceCenterGhosts (UserCtx *user, Vec local_vec)
 Applies geometric translations to wrapped face-center ghost coordinates.
 
static PetscErrorCode TransferPeriodicFaceFieldByDirection (UserCtx *user, FieldId field_id, char face_direction, char periodic_direction)
 Transfers one registered face-family field along one periodic axis.
 
PetscErrorCode SynchronizePeriodicFaceFields (UserCtx *user, char face_direction, PetscInt num_fields, const FieldId field_ids[])
 Synchronizes persistent fields belonging to one face family.
 
static PetscErrorCode GetPersistentStaggeredField (UserCtx *user, FieldId field_id, DM *dm, Vec *global_vec, Vec *local_vec)
 Resolves one registered persistent component-staggered field.
 
static PetscErrorCode TransferPeriodicStaggeredFieldByDirection (UserCtx *user, FieldId field_id, char periodic_direction)
 Transfers one component-staggered field along one periodic axis.
 
PetscErrorCode SynchronizePeriodicStaggeredFields (UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
 Implementation of SynchronizePeriodicStaggeredFields().
 
PetscErrorCode PreparePeriodicQuickStencilFields (UserCtx *user, Vec local_vector_field, Vec local_scalar_field)
 Implementation of PreparePeriodicQuickStencilFields().
 
PetscErrorCode SynchronizePeriodicLocalStaggeredField (UserCtx *user, Vec local_field)
 Implementation of SynchronizePeriodicLocalStaggeredField().
 
PetscErrorCode ApplyMetricsPeriodicBCs (UserCtx *user)
 Internal helper implementation: ApplyMetricsPeriodicBCs().
 
PetscErrorCode ApplyPeriodicBCs (UserCtx *user)
 Internal helper implementation: ApplyPeriodicBCs().
 
PetscErrorCode UpdateDummyCells (UserCtx *user)
 Internal helper implementation: UpdateDummyCells().
 
PetscErrorCode UpdateCornerNodes (UserCtx *user)
 Internal helper implementation: UpdateCornerNodes().
 
PetscErrorCode ApplyWallFunction (UserCtx *user)
 Internal helper implementation: ApplyWallFunction().
 
PetscErrorCode FinalizePostProjectionCellFields (UserCtx *user)
 Implementation of FinalizePostProjectionCellFields().
 
PetscErrorCode ApplyBoundaryConditions (UserCtx *user)
 Implementation of ApplyBoundaryConditions().
 

Macro Definition Documentation

◆ __FUNCT__ [1/30]

#define __FUNCT__   "CanRankServiceInletFace"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [2/30]

#define __FUNCT__   "CanRankServiceFace"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [3/30]

#define __FUNCT__   "GetDeterministicFaceGridLocation"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [4/30]

#define __FUNCT__   "GetRandomFCellAndLogicOnInletFace"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [5/30]

#define __FUNCT__   "ClassifyMomentumRow"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [6/30]

#define __FUNCT__   "EnforceRHSBoundaryConditions"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [7/30]

#define __FUNCT__   "BoundaryCondition_Create"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [8/30]

#define __FUNCT__   "BoundarySystem_Validate"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [9/30]

#define __FUNCT__   "BoundarySystem_Initialize"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [10/30]

#define __FUNCT__   "PropagateBoundaryConfigToCoarserLevels"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [11/30]

#define __FUNCT__   "BoundarySystem_ExecuteStep"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [12/30]

#define __FUNCT__   "BoundarySystem_RefreshUbcs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [13/30]

#define __FUNCT__   "BoundarySystem_Destroy"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [14/30]

#define __FUNCT__   "TransferPeriodicFieldByDirection"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [15/30]

#define __FUNCT__   "SynchronizePeriodicCellFields"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [16/30]

#define __FUNCT__   "GetPersistentFaceField"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [17/30]

#define __FUNCT__   "TransferPeriodicFaceFieldByDirection"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [18/30]

#define __FUNCT__   "SynchronizePeriodicFaceFields"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [19/30]

#define __FUNCT__   "GetPersistentStaggeredField"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [20/30]

#define __FUNCT__   "TransferPeriodicStaggeredFieldByDirection"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [21/30]

#define __FUNCT__   "SynchronizePeriodicStaggeredFields"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [22/30]

#define __FUNCT__   "PreparePeriodicQuickStencilFields"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [23/30]

#define __FUNCT__   "SynchronizePeriodicLocalStaggeredField"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [24/30]

#define __FUNCT__   "ApplyMetricsPeriodicBCs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [25/30]

#define __FUNCT__   "ApplyPeriodicBCs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [26/30]

#define __FUNCT__   "UpdateDummyCells"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [27/30]

#define __FUNCT__   "UpdateCornerNodes"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [28/30]

#define __FUNCT__   "ApplyWallFunction"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [29/30]

#define __FUNCT__   "FinalizePostProjectionCellFields"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [30/30]

#define __FUNCT__   "ApplyBoundaryConditions"

Definition at line 6 of file Boundaries.c.

Function Documentation

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

Internal helper implementation: CanRankServiceInletFace().

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.

Local to this translation unit.

Definition at line 11 of file Boundaries.c.

14{
15 PetscErrorCode ierr;
16 PetscMPIInt rank_for_logging; // For detailed debugging logs
17 PetscFunctionBeginUser;
18 PROFILE_FUNCTION_BEGIN;
19
20 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
21
22 *can_service_inlet_out = PETSC_FALSE; // Default to no service
23
24 if (!user->inletFaceDefined) {
25 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Inlet face not defined in user context. Cannot service.\n", rank_for_logging);
26 PROFILE_FUNCTION_END;
27 PetscFunctionReturn(0);
28 }
29
30 // Get the range of cells owned by this rank in each dimension
31 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
32 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
33 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
34
35 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
36 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
37 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
38
39 // Determine the global index of the last cell (0-indexed) in each direction.
40 // Example: If IM_nodes_global = 11 (nodes 0-10), there are 10 cells (0-9). Last cell index is 9.
41 // Formula: global_nodes - 1 (num cells) - 1 (0-indexed) = global_nodes - 2.
42 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)
43 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
44 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
45
46 switch (user->identifiedInletBCFace) {
47 case BC_FACE_NEG_X: // Inlet on the global I-minimum face (face of cell C_i=0)
48 // Rank services if its first owned node is global node 0 (info->xs == 0),
49 // and it owns cells in I, J, and K directions.
50 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
51 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
52 *can_service_inlet_out = PETSC_TRUE;
53 }
54 break;
55 case BC_FACE_POS_X: // Inlet on the global I-maximum face (face of cell C_i=last_global_cell_idx_i)
56 // Rank services if it owns the last cell in I-direction,
57 // and has extent in J and K.
58 if (last_global_cell_idx_i >= 0 && /* Check for valid global domain */
59 (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 */
60 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
61 *can_service_inlet_out = PETSC_TRUE;
62 }
63 break;
64 case BC_FACE_NEG_Y:
65 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
66 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
67 *can_service_inlet_out = PETSC_TRUE;
68 }
69 break;
70 case BC_FACE_POS_Y:
71 if (last_global_cell_idx_j >= 0 &&
72 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
73 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
74 *can_service_inlet_out = PETSC_TRUE;
75 }
76 break;
77 case BC_FACE_NEG_Z:
78 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
79 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
80 *can_service_inlet_out = PETSC_TRUE;
81 }
82 break;
83 case BC_FACE_POS_Z:
84 if (last_global_cell_idx_k >= 0 &&
85 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
86 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
87 *can_service_inlet_out = PETSC_TRUE;
88 }
89 break;
90 default:
91 LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d]: Unknown inlet face %s.\n", rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace));
92 break;
93 }
94
95 LOG_ALLOW(LOCAL, LOG_TRACE,
96 "[Rank %d] Check Service for Inlet %s:\n"
97 " - Local Domain: starts at cell (%d,%d,%d), has (%d,%d,%d) cells.\n"
98 " - Global Domain: has (%d,%d,%d) nodes, so last cell is (%d,%d,%d).\n",
99 rank_for_logging,
100 BCFaceToString((BCFace)user->identifiedInletBCFace),
101 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
102 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
103 IM_nodes_global, JM_nodes_global, KM_nodes_global,
104 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
105
106 LOG_ALLOW(LOCAL, LOG_INFO,"[Rank %d] Inlet Face %s Service Check Result: %s | Owned Cells (I,J,K): (%d,%d,%d) | Starts at Cell (%d,%d,%d)\n",
107 rank_for_logging,
108 BCFaceToString((BCFace)user->identifiedInletBCFace),
109 (*can_service_inlet_out) ? "CAN SERVICE" : "CANNOT SERVICE",
110 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
111 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k);
112
113 PROFILE_FUNCTION_END;
114
115 PetscFunctionReturn(0);
116}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
BCFaceToString
const char * BCFaceToString(BCFace face)
Returns the canonical log token for a boundary-face enum value.
Definition logging.c:671
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:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:859
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:33
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:850
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:2285
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:932
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:933
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:261
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:262
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:264
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:263
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:264
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:262
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:263
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 
)

Implementation of CanRankServiceFace().

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

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Boundaries.h.

See also
CanRankServiceFace()

Definition at line 127 of file Boundaries.c.

129{
130 PetscErrorCode ierr;
131 PetscMPIInt rank_for_logging;
132 PetscFunctionBeginUser;
133
134 PROFILE_FUNCTION_BEGIN;
135
136 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
137
138 *can_service_out = PETSC_FALSE; // Default to no service
139
140 // Get the range of cells owned by this rank
141 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
142 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
143 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
144 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
145 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
146 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
147
148 // Determine the global index of the last cell (0-indexed) in each direction.
149 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
150 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
151 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
152
153 switch (face_id) {
154 case BC_FACE_NEG_X:
155 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
156 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
157 *can_service_out = PETSC_TRUE;
158 }
159 break;
160 case BC_FACE_POS_X:
161 if (last_global_cell_idx_i >= 0 &&
162 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i &&
163 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
164 *can_service_out = PETSC_TRUE;
165 }
166 break;
167 case BC_FACE_NEG_Y:
168 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
169 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
170 *can_service_out = PETSC_TRUE;
171 }
172 break;
173 case BC_FACE_POS_Y:
174 if (last_global_cell_idx_j >= 0 &&
175 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
176 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
177 *can_service_out = PETSC_TRUE;
178 }
179 break;
180 case BC_FACE_NEG_Z:
181 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
182 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
183 *can_service_out = PETSC_TRUE;
184 }
185 break;
186 case BC_FACE_POS_Z:
187 if (last_global_cell_idx_k >= 0 &&
188 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
189 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
190 *can_service_out = PETSC_TRUE;
191 }
192 break;
193 default:
194 LOG_ALLOW(LOCAL, LOG_WARNING, "Rank %d: Unknown face enum %d. \n", rank_for_logging, face_id);
195 break;
196 }
197
198 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d check for face %s: Result=%s. \n",
199 rank_for_logging, BCFaceToString((BCFace)face_id), (*can_service_out ? "TRUE" : "FALSE"));
200
201 PROFILE_FUNCTION_END;
202
203 PetscFunctionReturn(0);
204}
Here is the call graph for this function:
Here is the caller 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 
)

Internal helper implementation: GetDeterministicFaceGridLocation().

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

Local to this translation unit.

Definition at line 213 of file Boundaries.c.

221{
222 SimCtx *simCtx = user->simCtx;
223 PetscReal global_logic_i = 0.0, global_logic_j = 0.0, global_logic_k = 0.0;
224 PetscErrorCode ierr;
225 PetscMPIInt rank_for_logging;
226
227 PetscFunctionBeginUser;
228 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
229
230 *placement_successful_out = PETSC_FALSE; // Default to failure
231
232 // --- Step 1: Configuration and Input Validation ---
233
234 // *** Hardcoded number of grid layers. Change this value to alter the pattern. ***
235 const PetscInt grid_layers = 2;
236
237 LOG_ALLOW(LOCAL, LOG_DEBUG,
238 "[Rank %d] Placing particle %lld on face %s with grid_layers=%d in global domain (%d,%d,%d) cells.\n",
239 rank_for_logging, (long long)particle_global_id, BCFaceToString(user->identifiedInletBCFace), grid_layers,
240 IM_cells_global, JM_cells_global, KM_cells_global);
241
242 const char *face_name = BCFaceToString(user->identifiedInletBCFace);
243
244 // Fatal Error Checks: Ensure the requested grid is geometrically possible.
245 // The total layers from opposite faces (2 * grid_layers) must be less than the domain size.
246 switch (user->identifiedInletBCFace) {
247 case BC_FACE_NEG_X: case BC_FACE_POS_X:
248 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);
249 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);
250 break;
251 case BC_FACE_NEG_Y: case BC_FACE_POS_Y:
252 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);
253 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);
254 break;
255 case BC_FACE_NEG_Z: case BC_FACE_POS_Z:
256 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);
257 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);
258 break;
259 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid identifiedInletBCFace specified: %d", user->identifiedInletBCFace);
260 }
261
262 const PetscInt num_lines_total = 4 * grid_layers;
263 if (simCtx->np < num_lines_total) {
264 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);
265 }
266 if (simCtx->np > 0 && simCtx->np % num_lines_total != 0) {
267 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);
268 }
269
270 // --- Step 2: Map global particle ID to a line and a point on that line ---
271 if (simCtx->np == 0) PetscFunctionReturn(0); // Nothing to do
272
273 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Distributing %lld particles over %d lines on face %s.\n",
274 rank_for_logging, (long long)simCtx->np, num_lines_total, face_name);
275
276 const PetscInt points_per_line = PetscMax(1, simCtx->np / num_lines_total);
277 PetscInt line_index = particle_global_id / points_per_line;
278 PetscInt point_index_on_line = particle_global_id % points_per_line;
279 line_index = PetscMin(line_index, num_lines_total - 1); // Clamp to handle uneven division
280
281 // Decode the line_index into an edge group (0-3) and a layer within that group (0 to grid_layers-1)
282 const PetscInt edge_group = line_index / grid_layers;
283 const PetscInt layer_index = line_index % grid_layers;
284
285 // --- Step 3: Calculate placement coordinates based on the decoded indices ---
286 const PetscReal layer_spacing_norm_i = (IM_cells_global > 0) ? 1.0 / (PetscReal)IM_cells_global : 0.0;
287 const PetscReal layer_spacing_norm_j = (JM_cells_global > 0) ? 1.0 / (PetscReal)JM_cells_global : 0.0;
288 const PetscReal layer_spacing_norm_k = (KM_cells_global > 0) ? 1.0 / (PetscReal)KM_cells_global : 0.0;
289
290 // Grid-aware epsilon: scale with minimum cell size to keep particles away from rank boundaries
291 const PetscReal min_layer_spacing = PetscMin(layer_spacing_norm_i, PetscMin(layer_spacing_norm_j, layer_spacing_norm_k));
292 const PetscReal epsilon = 0.5 * min_layer_spacing; // Keep particles 10% of cell width from boundaries
293
294 PetscReal variable_coord; // The coordinate that varies along a line
295 if (points_per_line <= 1) {
296 variable_coord = 0.5; // Place single point in the middle
297 } else {
298 variable_coord = ((PetscReal)point_index_on_line + 0.5)/ (PetscReal)(points_per_line);
299 }
300 variable_coord = PetscMin(1.0 - epsilon, PetscMax(epsilon, variable_coord)); // Clamp within [eps, 1-eps]
301
302 // Main logic switch to determine the three global logical coordinates
303 switch (user->identifiedInletBCFace) {
304 case BC_FACE_NEG_X:
305 global_logic_i = 0.5 * layer_spacing_norm_i; // Place near the face, in the middle of the first cell
306 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
307 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
308 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
309 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
310 break;
311 case BC_FACE_POS_X:
312 global_logic_i = 1.0 - (0.5 * layer_spacing_norm_i); // Place near the face, in the middle of the last cell
313 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
314 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
315 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
316 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
317 break;
318 case BC_FACE_NEG_Y:
319 global_logic_j = 0.5 * layer_spacing_norm_j;
320 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
321 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
322 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
323 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
324 break;
325 case BC_FACE_POS_Y:
326 global_logic_j = 1.0 - (0.5 * layer_spacing_norm_j);
327 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
328 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
329 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
330 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
331 break;
332 case BC_FACE_NEG_Z:
333 global_logic_k = 0.5 * layer_spacing_norm_k;
334 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
335 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
336 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
337 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
338 break;
339 case BC_FACE_POS_Z:
340 global_logic_k = 1.0 - (0.5 * layer_spacing_norm_k);
341 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
342 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
343 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
344 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
345 break;
346 }
347
348 LOG_ALLOW(LOCAL, LOG_TRACE,
349 "[Rank %d] Particle %lld assigned to line %d (edge group %d, layer %d) with variable_coord=%.4f.\n"
350 " -> Global logical coords: (i,j,k) = (%.6f, %.6f, %.6f)\n",
351 rank_for_logging, (long long)particle_global_id, line_index, edge_group, layer_index, variable_coord,
352 global_logic_i, global_logic_j, global_logic_k);
353
354 // --- Step 4: Convert global logical coordinate to global cell index and intra-cell logicals ---
355 PetscReal global_cell_coord_i = global_logic_i * IM_cells_global;
356 PetscInt I_g = (PetscInt)global_cell_coord_i;
357 *xi_metric_logic_out = global_cell_coord_i - I_g;
358
359 PetscReal global_cell_coord_j = global_logic_j * JM_cells_global;
360 PetscInt J_g = (PetscInt)global_cell_coord_j;
361 *eta_metric_logic_out = global_cell_coord_j - J_g;
362
363 PetscReal global_cell_coord_k = global_logic_k * KM_cells_global;
364 PetscInt K_g = (PetscInt)global_cell_coord_k;
365 *zta_metric_logic_out = global_cell_coord_k - K_g;
366
367 // --- Step 5: Check if this rank owns the target cell and finalize outputs ---
368 if ((I_g >= info->xs && I_g < info->xs + info->xm) &&
369 (J_g >= info->ys && J_g < info->ys + info->ym) &&
370 (K_g >= info->zs && K_g < info->zs + info->zm))
371 {
372 // Convert global cell index to the local node index for this rank's DA patch
373 *ci_metric_lnode_out = (I_g - info->xs) + xs_gnode_rank;
374 *cj_metric_lnode_out = (J_g - info->ys) + ys_gnode_rank;
375 *ck_metric_lnode_out = (K_g - info->zs) + zs_gnode_rank;
376 *placement_successful_out = PETSC_TRUE;
377 }
378
379 LOG_ALLOW(LOCAL, LOG_VERBOSE,
380 "[Rank %d] Particle %lld placement %s.\n",
381 rank_for_logging, (long long)particle_global_id,
382 (*placement_successful_out ? "SUCCESSFUL" : "NOT ON THIS RANK"));
383
384 if(*placement_successful_out){
385 LOG_ALLOW(LOCAL,LOG_TRACE,"Local cell origin node: (I,J,K) = (%d,%d,%d), intra-cell logicals: (xi,eta,zta)=(%.6f,%.6f,%.6f)\n",
386 *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
387 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
388 }
389
390 PetscFunctionReturn(0);
391}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:34
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
SimCtx::np
PetscInt np
Definition variables.h:827
SimCtx
The master context for the entire simulation.
Definition variables.h:695
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 
)

Internal helper implementation: GetRandomCellAndLogicalCoordsOnInletFace().

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.

Local to this translation unit.

Definition at line 400 of file Boundaries.c.

407{
408 PetscErrorCode ierr = 0;
409 PetscReal r_val_i_sel, r_val_j_sel, r_val_k_sel;
410 PetscInt local_cell_idx_on_face_dim1 = 0; // 0-indexed relative to owned cells on face
411 PetscInt local_cell_idx_on_face_dim2 = 0;
412 PetscMPIInt rank_for_logging;
413
414 PetscFunctionBeginUser;
415
416 PROFILE_FUNCTION_BEGIN;
417
418 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
419
420 // Get number of cells this rank owns in each dimension (tangential to the face mainly)
421 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
422 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
423 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
424
425 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
426 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
427 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
428
429 // Defaults for cell origin node (local index for the rank's DA patch, including ghosts)
430 *ci_metric_lnode_out = xs_gnode_rank; *cj_metric_lnode_out = ys_gnode_rank; *ck_metric_lnode_out = zs_gnode_rank;
431 // Defaults for logical coordinates
432 *xi_metric_logic_out = 0.5; *eta_metric_logic_out = 0.5; *zta_metric_logic_out = 0.5;
433
434 // Index of the last cell (0-indexed) in each global direction
435 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
436 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
437 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
438
439 LOG_ALLOW(LOCAL, LOG_INFO, "PARTICLE_INIT_DEBUG Rank %d: Inlet face %s.\n"
440 " Owned cells (i,j,k): (%d,%d,%d)\n"
441 " Global nodes (I,J,K): (%d,%d,%d)\n"
442 " info->xs,ys,zs (first owned node GLOBAL): (%d,%d,%d)\n"
443 " info->xm,ym,zm (num owned nodes GLOBAL): (%d,%d,%d)\n"
444 " xs_gnode_rank,ys_gnode_rank,zs_gnode_rank (DMDAGetCorners): (%d,%d,%d)\n"
445 " owned_start_cell (i,j,k) GLOBAL: (%d,%d,%d)\n"
446 " last_global_cell_idx (i,j,k): (%d,%d,%d)\n",
447 rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace),
448 num_owned_cells_on_rank_i,num_owned_cells_on_rank_j,num_owned_cells_on_rank_k,
449 IM_nodes_global,JM_nodes_global,KM_nodes_global,
450 info->xs, info->ys, info->zs,
451 info->xm, info->ym, info->zm,
452 xs_gnode_rank,ys_gnode_rank,zs_gnode_rank,
453 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
454 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
455
456
457 switch (user->identifiedInletBCFace) {
458 case BC_FACE_NEG_X: // Particle on -X face of cell C_0 (origin node N_0)
459 // Cell origin node is the first owned node in I by this rank (global index info->xs).
460 // Its local index within the rank's DA (incl ghosts) is xs_gnode_rank.
461 *ci_metric_lnode_out = xs_gnode_rank;
462 *xi_metric_logic_out = 1.0e-6;
463
464 // Tangential dimensions are J and K. Select an owned cell randomly on this face.
465 // num_owned_cells_on_rank_j/k must be > 0 (checked by CanRankServiceInletFace)
466 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
467 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j); // Index among owned J-cells
468 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
469 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1; // Offset from start of rank's J-nodes
470
471 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
472 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
473 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
474 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
475
476 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
477 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
478 break;
479
480 case BC_FACE_POS_X: // Particle on +X face of cell C_last_I (origin node N_last_I_origin)
481 // Origin node of the last I-cell is global_node_idx = last_global_cell_idx_i.
482 // Its local index in rank's DA: (last_global_cell_idx_i - info->xs) + xs_gnode_rank
483 *ci_metric_lnode_out = xs_gnode_rank + (last_global_cell_idx_i - info->xs);
484 *xi_metric_logic_out = 1.0 - 1.0e-6;
485
486 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
487 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
488 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
489 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1;
490
491 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
492 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
493 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
494 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
495
496 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
497 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
498 break;
499 // ... (Cases for Y and Z faces, following the same pattern) ...
500 case BC_FACE_NEG_Y:
501 *cj_metric_lnode_out = ys_gnode_rank;
502 *eta_metric_logic_out = 1.0e-6;
503 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
504 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
505 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
506 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
507 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
508 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
509 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
510 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
511 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
512 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
513 break;
514 case BC_FACE_POS_Y:
515 *cj_metric_lnode_out = ys_gnode_rank + (last_global_cell_idx_j - info->ys);
516 *eta_metric_logic_out = 1.0 - 1.0e-6;
517 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
518 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
519 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
520 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
521 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
522 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
523 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
524 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
525 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
526 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
527 break;
528 case BC_FACE_NEG_Z: // Your example case
529 *ck_metric_lnode_out = zs_gnode_rank; // Cell origin is the first owned node in K by this rank
530 *zta_metric_logic_out = 1.0e-6; // Place particle slightly inside this cell from its -Z face
531 // Tangential dimensions are I and J
532 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
533 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
534 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
535 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
536
537 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
538 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
539 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
540 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
541
542 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for I
543 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for J
544 break;
545 case BC_FACE_POS_Z:
546 *ck_metric_lnode_out = zs_gnode_rank + (last_global_cell_idx_k - info->zs);
547 *zta_metric_logic_out = 1.0 - 1.0e-6;
548 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
549 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
550 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
551 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
552 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
553 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
554 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
555 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
556 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
557 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
558 break;
559 default:
560 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "GetRandomCellAndLogicOnInletFace: Invalid user->identifiedInletBCFace %d. \n", user->identifiedInletBCFace);
561 }
562
563 PetscReal eps = 1.0e-7;
564 if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
565 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
566 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
567 } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
568 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
569 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
570 } else {
571 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
572 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
573 }
574
575 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Rank %d: Target Cell Node =(%d,%d,%d). (xi,et,zt)=(%.2e,%.2f,%.2f). \n",
576 rank_for_logging, *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
577 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
578
579 PROFILE_FUNCTION_END;
580
581 PetscFunctionReturn(0);
582}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ClassifyMomentumRow()

MomentumRowType ClassifyMomentumRow ( UserCtx *  user,
PetscInt  i,
PetscInt  j,
PetscInt  k,
PetscInt  component,
PetscInt *  ri,
PetscInt *  rj,
PetscInt *  rk 
)

Implementation of ClassifyMomentumRow().

Single source of truth for "which staggered momentum rows are unknowns".

Pure index/boundary-type arithmetic: no field reads, no communication. Precedence matters. A conditioned row is reported first because its explicit Dirichlet value is more specific than the homogeneous fallback, and a periodic duplicate is reported before the homogeneous case because the Newton path needs its representative index to build F = X_dup - X_rep.

Definition at line 597 of file Boundaries.c.

599{
600 const PetscInt mx = user->info.mx, my = user->info.my, mz = user->info.mz;
601 const PetscInt coord[3] = {i, j, k};
602 const PetscInt size[3] = {mx, my, mz};
603 const BCFace neg_face[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
604 PetscBool periodic[3], periodic_duplicate = PETSC_FALSE;
605 PetscBool residual_zeroed = PETSC_FALSE, conditioned = PETSC_FALSE;
606
607 *ri = i; *rj = j; *rk = k;
608 for (PetscInt axis = 0; axis < 3; ++axis) {
609 periodic[axis] = (PetscBool)(
610 user->boundary_faces[neg_face[axis]].mathematical_type == PERIODIC);
611 if (periodic[axis] && coord[axis] == 0) {
612 periodic_duplicate = PETSC_TRUE;
613 if (axis == 0) *ri = -2;
614 else if (axis == 1) *rj = -2;
615 else *rk = -2;
616 }
617 if (periodic[axis] && coord[axis] == size[axis] - 1) {
618 periodic_duplicate = PETSC_TRUE;
619 if (axis == 0) *ri = mx + 1;
620 else if (axis == 1) *rj = my + 1;
621 else *rk = mz + 1;
622 }
623
624 if (!periodic[axis] && coord[axis] == 0) residual_zeroed = PETSC_TRUE;
625 if (coord[axis] == size[axis] - 1) residual_zeroed = PETSC_TRUE;
626 if (!periodic[axis] && coord[axis] == size[axis] - 2 && component == axis)
627 residual_zeroed = PETSC_TRUE;
628 }
629
630 if (!periodic[component] &&
631 (coord[component] == 0 || coord[component] == size[component] - 2)) {
632 PetscBool tangential_interior = PETSC_TRUE;
633 for (PetscInt axis = 0; axis < 3; ++axis) {
634 if (axis == component) continue;
635 if (coord[axis] < 1 || coord[axis] > size[axis] - 2)
636 tangential_interior = PETSC_FALSE;
637 }
638 conditioned = tangential_interior;
639 }
640
641 if (conditioned) return MOM_ROW_FIXED_CONDITIONED;
642 if (periodic_duplicate) return MOM_ROW_PERIODIC_DUPLICATE;
643 if (residual_zeroed) return MOM_ROW_FIXED_HOMOGENEOUS;
644 return MOM_ROW_PHYSICAL;
645}
MOM_ROW_FIXED_HOMOGENEOUS
@ MOM_ROW_FIXED_HOMOGENEOUS
Dummy/tangential row carrying no unknown at all.
Definition Boundaries.h:255
MOM_ROW_PHYSICAL
@ MOM_ROW_PHYSICAL
Independent unknown governed by the momentum equation.
Definition Boundaries.h:253
MOM_ROW_PERIODIC_DUPLICATE
@ MOM_ROW_PERIODIC_DUPLICATE
Duplicate of a wrapped representative row (see ri, rj, rk).
Definition Boundaries.h:256
MOM_ROW_FIXED_CONDITIONED
@ MOM_ROW_FIXED_CONDITIONED
Strong Dirichlet row; the value comes from ApplyBoundaryConditions().
Definition Boundaries.h:254
PERIODIC
@ PERIODIC
Definition variables.h:292
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:931
UserCtx::info
DMDALocalInfo info
Definition variables.h:918
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:368
Here is the caller graph for this function:

◆ EnforceRHSBoundaryConditions()

PetscErrorCode EnforceRHSBoundaryConditions ( UserCtx *  user)

Implementation of EnforceRHSBoundaryConditions().

Zeroes every momentum RHS row that does not carry an independent unknown.

The sweep is deliberately expressed over every owned location rather than over the six boundary slabs: restating "which indices can be non-physical" here is exactly the duplication that let the periodic duplicate column go unzeroed. ClassifyMomentumRow() is a handful of integer comparisons and the walk is a single pass with no stencil access, which is negligible next to the several ghosted stencil passes ComputeRHS() has already made over the same range. MomentumNewtonKrylov_ApplyConstraints() walks the same range the same way.

Definition at line 660 of file Boundaries.c.

661{
662 PetscErrorCode ierr;
663 DMDALocalInfo info = user->info;
664 Cmpnts ***rhs;
665
666 PetscFunctionBeginUser;
667 PROFILE_FUNCTION_BEGIN;
668
669 // Get a writable pointer to the local data of the global RHS vector.
670 ierr = DMDAVecGetArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
671
672 for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
673 for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
674 for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
675 PetscScalar *row = &rhs[k][j][i].x; /* .x/.y/.z are contiguous */
676 for (PetscInt component = 0; component < 3; component++) {
677 PetscInt ri, rj, rk;
678 if (ClassifyMomentumRow(user, i, j, k, component, &ri, &rj, &rk) != MOM_ROW_PHYSICAL)
679 row[component] = 0.0;
680 }
681 }
682 }
683 }
684
685 // --- Release the pointer to the local data ---
686 ierr = DMDAVecRestoreArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
687
688 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d, Block %d: Finished enforcing RHS boundary conditions.\n",
689 user->simCtx->rank, user->_this);
690
691 PROFILE_FUNCTION_END;
692
693 PetscFunctionReturn(0);
694}
ClassifyMomentumRow
MomentumRowType ClassifyMomentumRow(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscInt component, PetscInt *ri, PetscInt *rj, PetscInt *rk)
Implementation of ClassifyMomentumRow().
Definition Boundaries.c:597
UserCtx::Rhs
Vec Rhs
Definition variables.h:947
SimCtx::rank
PetscMPIInt rank
Definition variables.h:698
UserCtx::_this
PetscInt _this
Definition variables.h:924
Cmpnts::x
PetscScalar x
Definition variables.h:103
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundaryCondition_Create()

PetscErrorCode BoundaryCondition_Create ( BCHandlerType  handler_type,
BoundaryCondition **  new_bc_ptr 
)

Internal helper implementation: BoundaryCondition_Create().

(Private) Creates and configures a specific BoundaryCondition handler object.

Local to this translation unit.

Definition at line 703 of file Boundaries.c.

704{
705 PetscErrorCode ierr;
706 PetscFunctionBeginUser;
707
708 const char* handler_name = BCHandlerTypeToString(handler_type);
709 LOG_ALLOW(LOCAL, LOG_DEBUG, "Factory called for handler type %s. \n", handler_name);
710
711 ierr = PetscMalloc1(1, new_bc_ptr); CHKERRQ(ierr);
712 BoundaryCondition *bc = *new_bc_ptr;
713
714 bc->type = handler_type;
715 bc->priority = -1; // Default priority; can be overridden in specific handlers
716 bc->data = NULL;
717 bc->Initialize = NULL;
718 bc->PreStep = NULL;
719 bc->Apply = NULL;
720 bc->PostStep = NULL;
721 bc->UpdateUbcs = NULL;
722 bc->Destroy = NULL;
723
724 LOG_ALLOW(LOCAL, LOG_DEBUG, "Allocated generic handler object at address %p.\n", (void*)bc);
725
726 switch (handler_type) {
727
728 case BC_HANDLER_OUTLET_CONSERVATION:
729 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_OutletConservation().\n");
730 ierr = Create_OutletConservation(bc); CHKERRQ(ierr);
731 break;
732
733 case BC_HANDLER_WALL_NOSLIP:
734 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_WallNoSlip().\n");
735 ierr = Create_WallNoSlip(bc); CHKERRQ(ierr);
736 break;
737
738 case BC_HANDLER_INLET_CONSTANT_VELOCITY:
739 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletConstantVelocity().\n");
740 ierr = Create_InletConstantVelocity(bc); CHKERRQ(ierr);
741 break;
742
743 case BC_HANDLER_PERIODIC_GEOMETRIC:
744 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicGeometric().\n");
745 ierr = Create_PeriodicGeometric(bc);
746 break;
747
748 case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX:
749 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenConstant().\n");
750 ierr = Create_PeriodicDrivenConstant(bc);
751 break;
752
753 case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX:
754 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenInitial().\n");
755 ierr = Create_PeriodicDrivenInitial(bc);
756 break;
757
758 case BC_HANDLER_INLET_PARABOLIC:
759 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletParabolicProfile().\n");
760 ierr = Create_InletParabolicProfile(bc); CHKERRQ(ierr);
761 break;
762
763 case BC_HANDLER_INLET_PROFILE_FROM_FILE:
764 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletProfileFromFile().\n");
765 ierr = Create_InletProfileFromFile(bc); CHKERRQ(ierr);
766 break;
767 //Add cases for other handlers here in future phases
768
769 default:
770 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type (%s) is not recognized or implemented in the factory.\n", handler_name);
771 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) not recognized in factory.\n", handler_type, handler_name);
772 }
773
774 if(bc->priority < 0) {
775 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
776 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
777 }
778
779 LOG_ALLOW(LOCAL, LOG_DEBUG, "Successfully created and configured handler for %s.\n", handler_name);
780 PetscFunctionReturn(0);
781}
Create_InletConstantVelocity
PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
Configures a BoundaryCondition object to behave as a constant velocity inlet.
Definition BC_Handlers.c:322
Create_InletProfileFromFile
PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
Configures a BoundaryCondition object for a file-prescribed inlet profile.
Definition BC_Handlers.c:1292
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Configures a BoundaryCondition object for geometric periodic coupling.
Definition BC_Handlers.c:2159
Create_InletParabolicProfile
PetscErrorCode Create_InletParabolicProfile(BoundaryCondition *bc)
Configures a BoundaryCondition object for a parabolic inlet profile.
Definition BC_Handlers.c:715
Create_PeriodicDrivenInitial
PetscErrorCode Create_PeriodicDrivenInitial(BoundaryCondition *bc)
Configures a BoundaryCondition object for initial-flux periodic driving.
Definition BC_Handlers.c:2669
Create_PeriodicDrivenConstant
PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
Configures a BoundaryCondition object for periodic driven-flow forcing.
Definition BC_Handlers.c:2325
Create_WallNoSlip
PetscErrorCode Create_WallNoSlip(BoundaryCondition *bc)
Configures a BoundaryCondition object to behave as a no-slip, stationary wall.
Definition BC_Handlers.c:114
Create_OutletConservation
PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
Configures a BoundaryCondition object for conservative outlet treatment.
Definition BC_Handlers.c:1644
BCHandlerTypeToString
const char * BCHandlerTypeToString(BCHandlerType handler_type)
Converts a BCHandlerType enum to its string representation.
Definition logging.c:793
LOG_ERROR
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:29
BoundaryCondition
The "virtual table" struct for a boundary condition handler object.
Definition variables.h:353
BoundaryCondition::PostStep
PetscErrorCode(* PostStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:360
BoundaryCondition::PreStep
PetscErrorCode(* PreStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:358
BoundaryCondition::type
BCHandlerType type
Definition variables.h:354
BoundaryCondition::Destroy
PetscErrorCode(* Destroy)(BoundaryCondition *self)
Definition variables.h:362
BoundaryCondition::Initialize
PetscErrorCode(* Initialize)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:357
BoundaryCondition::UpdateUbcs
PetscErrorCode(* UpdateUbcs)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:361
BoundaryCondition::Apply
PetscErrorCode(* Apply)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:359
BoundaryCondition::priority
BCPriorityType priority
Definition variables.h:355
BC_HANDLER_PERIODIC_GEOMETRIC
@ BC_HANDLER_PERIODIC_GEOMETRIC
Definition variables.h:316
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:309
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:308
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:319
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:318
BC_HANDLER_INLET_PROFILE_FROM_FILE
@ BC_HANDLER_INLET_PROFILE_FROM_FILE
Definition variables.h:310
BC_HANDLER_WALL_NOSLIP
@ BC_HANDLER_WALL_NOSLIP
Definition variables.h:305
BC_HANDLER_OUTLET_CONSERVATION
@ BC_HANDLER_OUTLET_CONSERVATION
Definition variables.h:314
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundarySystem_Validate()

PetscErrorCode BoundarySystem_Validate ( UserCtx *  user)

Internal helper implementation: BoundarySystem_Validate().

(Public) Validates the consistency and compatibility of the parsed boundary condition system.

Local to this translation unit.

Definition at line 789 of file Boundaries.c.

790{
791 PetscErrorCode ierr;
792 const BCFace neg_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
793 const BCFace pos_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
794 const char axis_names[3] = {'X', 'Y', 'Z'};
795 DMBoundaryType bx, by, bz;
796 PetscBool dm_periodic[3];
797 PetscFunctionBeginUser;
798
799 LOG_ALLOW(GLOBAL, LOG_INFO, "Validating parsed boundary condition configuration...\n");
800 ierr = DMDAGetInfo(user->da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
801 &bx, &by, &bz, NULL); CHKERRQ(ierr);
802 dm_periodic[0] = (PetscBool)(bx == DM_BOUNDARY_PERIODIC);
803 dm_periodic[1] = (PetscBool)(by == DM_BOUNDARY_PERIODIC);
804 dm_periodic[2] = (PetscBool)(bz == DM_BOUNDARY_PERIODIC);
805
806 // --- Rule Set 1: Geometric periodic faces must be paired and match the DM topology. ---
807 for (PetscInt axis = 0; axis < 3; axis++) {
808 const PetscBool neg_periodic =
809 user->boundary_faces[neg_faces[axis]].mathematical_type == PERIODIC;
810 const PetscBool pos_periodic =
811 user->boundary_faces[pos_faces[axis]].mathematical_type == PERIODIC;
812
813 PetscCheck(neg_periodic == pos_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
814 "Configuration Error: Periodic boundaries in the %c direction must be paired; "
815 "%s is %s while %s is %s.",
816 axis_names[axis],
817 BCFaceToString(neg_faces[axis]), neg_periodic ? "PERIODIC" : "not periodic",
818 BCFaceToString(pos_faces[axis]), pos_periodic ? "PERIODIC" : "not periodic");
819 PetscCheck(dm_periodic[axis] == neg_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
820 "Configuration Error: The %c-direction DM periodic flag (%d) does not match "
821 "the paired boundary configuration (%s).",
822 axis_names[axis], (int)dm_periodic[axis], neg_periodic ? "PERIODIC" : "not periodic");
823 }
824
825 // --- Rule Set 2: Driven Flow Handler Consistency ---
826 // This specialized validator will check all rules related to driven flow handlers.
827 ierr = Validate_DrivenFlowConfiguration(user); CHKERRQ(ierr);
828
829 // --- Rule Set 3: (Future Extension) Overset Interface Consistency ---
830 // ierr = Validate_OversetConfiguration(user); CHKERRQ(ierr);
831
832 LOG_ALLOW(GLOBAL, LOG_INFO, "Boundary configuration is valid.\n");
833
834 PetscFunctionReturn(0);
835}
Validate_DrivenFlowConfiguration
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
(Private) Validates all consistency rules for a driven flow (channel/pipe) setup.
Definition BC_Handlers.c:15
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundarySystem_Initialize()

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

Implementation of BoundarySystem_Initialize().

Initializes the entire boundary system.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Boundaries.h.

See also
BoundarySystem_Initialize()

Definition at line 850 of file Boundaries.c.

851{
852 PetscErrorCode ierr;
853 PetscFunctionBeginUser;
854
855 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting creation and initialization of all boundary handlers.\n");
856
857 // =========================================================================
858 // Step 0: Clear any existing boundary handlers (if re-initializing).
859 // This ensures no memory leaks if this function is called multiple times.
860 // =========================================================================
861 for (int i = 0; i < 6; i++) {
862 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
863 if (face_cfg->handler) {
864 LOG_ALLOW(LOCAL, LOG_DEBUG, "Destroying existing handler on Face %s before re-initialization.\n", BCFaceToString((BCFace)i));
865 if (face_cfg->handler->Destroy) {
866 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
867 }
868 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
869 face_cfg->handler = NULL;
870 }
871 }
872 // =========================================================================
873
874 // Step 0.1: Initiate flux sums to zero
875 user->simCtx->FluxInSum = 0.0;
876 user->simCtx->FluxOutSum = 0.0;
877 user->simCtx->FarFluxInSum = 0.0;
878 user->simCtx->FarFluxOutSum = 0.0;
879 // =========================================================================
880
881 // Step 1: Parse the configuration file to determine user intent.
882 // This function, defined in io.c, populates the configuration enums and parameter
883 // lists within the user->boundary_faces array on all MPI ranks.
884 ierr = ParseAllBoundaryConditions(user, bcs_filename); CHKERRQ(ierr);
885 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration file '%s' parsed successfully.\n", bcs_filename);
886
887 // Step 1.1: Validate the parsed configuration to ensure there are no Boundary Condition conflicts
888 ierr = BoundarySystem_Validate(user); CHKERRQ(ierr);
889
890 // Step 2: Create and Initialize the handler object for each of the 6 faces.
891 for (int i = 0; i < 6; i++) {
892 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
893
894 const char *face_name = BCFaceToString(face_cfg->face_id);
895 const char *type_name = BCTypeToString(face_cfg->mathematical_type);
896 const char *handler_name = BCHandlerTypeToString(face_cfg->handler_type);
897
898 LOG_ALLOW(LOCAL, LOG_DEBUG, "Creating handler for Face %s with Type %s and handler '%s'.\n", face_name, type_name,handler_name);
899
900 // Use the private factory to construct the correct handler object based on the parsed type.
901 // The factory returns a pointer to the new handler object, which we store in the config struct.
902 ierr = BoundaryCondition_Create(face_cfg->handler_type, &face_cfg->handler); CHKERRQ(ierr);
903
904 // Step 3: Call the specific Initialize() method for the newly created handler.
905 // This allows the handler to perform its own setup, like reading parameters from the
906 // face_cfg->params list and setting the initial field values on its face.
907 if (face_cfg->handler && face_cfg->handler->Initialize) {
908 LOG_ALLOW(LOCAL, LOG_DEBUG, "Calling Initialize() method for handler %s(%s) on Face %s.\n",type_name,handler_name,face_name);
909
910 // Prepare the context needed by the Initialize() function.
911 BCContext ctx = {
912 .user = user,
913 .face_id = face_cfg->face_id,
914 .global_inflow_sum = &user->simCtx->FluxInSum, // Global flux sums are not relevant during initialization.
915 .global_outflow_sum = &user->simCtx->FluxOutSum,
916 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
917 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
918 };
919
920 ierr = face_cfg->handler->Initialize(face_cfg->handler, &ctx); CHKERRQ(ierr);
921 } else {
922 LOG_ALLOW(LOCAL, LOG_DEBUG, "Handler %s(%s) for Face %s has no Initialize() method, skipping.\n", type_name,handler_name,face_name);
923 }
924 }
925 // =========================================================================
926 // NO SYNCHRONIZATION NEEDED HERE
927 // =========================================================================
928 // Initialize() only reads parameters and allocates memory.
929 // It does NOT modify field values (Ucat, Ucont, Ubcs).
930 // Field values are set by:
931 // 1. Initial conditions (before this function)
932 // 2. Apply() during timestepping (after this function)
933 // The first call to ApplyBoundaryConditions() will handle synchronization.
934 // =========================================================================
935
936 LOG_ALLOW(GLOBAL, LOG_INFO, "All boundary handlers created and initialized successfully.\n");
937 PetscFunctionReturn(0);
938}
BoundarySystem_Validate
PetscErrorCode BoundarySystem_Validate(UserCtx *user)
Internal helper implementation: BoundarySystem_Validate().
Definition Boundaries.c:789
BoundaryCondition_Create
PetscErrorCode BoundaryCondition_Create(BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
Internal helper implementation: BoundaryCondition_Create().
Definition Boundaries.c:703
ParseAllBoundaryConditions
PetscErrorCode ParseAllBoundaryConditions(UserCtx *user, const char *bcs_input_filename)
Parses the boundary conditions file to configure the type, handler, and any associated parameters for...
Definition io.c:837
BCTypeToString
const char * BCTypeToString(BCType type)
Returns the canonical log token for a boundary mathematical type.
Definition logging.c:773
SimCtx::FarFluxInSum
PetscReal FarFluxInSum
Definition variables.h:799
SimCtx::FarFluxOutSum
PetscReal FarFluxOutSum
Definition variables.h:799
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:799
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:369
BCContext::user
UserCtx * user
Definition variables.h:344
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:799
BoundaryFaceConfig::handler
BoundaryCondition * handler
Definition variables.h:371
BCContext
Provides execution context for a boundary condition handler.
Definition variables.h:343
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:366
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PropagateBoundaryConfigToCoarserLevels()

PetscErrorCode PropagateBoundaryConfigToCoarserLevels ( SimCtx *  simCtx)

Internal helper implementation: PropagateBoundaryConfigToCoarserLevels().

Propagates boundary condition configuration from finest to all coarser multigrid levels.

Local to this translation unit.

Definition at line 947 of file Boundaries.c.

948{
949 PetscErrorCode ierr;
950 UserMG *usermg = &simCtx->usermg;
951
952 PetscFunctionBeginUser;
953 PROFILE_FUNCTION_BEGIN;
954
955 LOG_ALLOW(GLOBAL, LOG_INFO, "Propagating BC configuration from finest to coarser multigrid levels...\n");
956
957 // Loop from second-finest down to coarsest
958 for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
959 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
960 UserCtx *user_coarse = &usermg->mgctx[level].user[bi];
961 UserCtx *user_fine = &usermg->mgctx[level + 1].user[bi];
962
963 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Copying BC config from level %d to level %d, block %d\n",
964 simCtx->rank, level + 1, level, bi);
965
966 // Copy the 6 boundary face configurations
967 for (int face_i = 0; face_i < 6; face_i++) {
968 user_coarse->boundary_faces[face_i].face_id = user_fine->boundary_faces[face_i].face_id;
969 user_coarse->boundary_faces[face_i].mathematical_type = user_fine->boundary_faces[face_i].mathematical_type;
970 user_coarse->boundary_faces[face_i].handler_type = user_fine->boundary_faces[face_i].handler_type;
971
972 // Copy parameter list (deep copy)
973 FreeBC_ParamList(user_coarse->boundary_faces[face_i].params); // Clear any existing
974 user_coarse->boundary_faces[face_i].params = NULL;
975
976 BC_Param **dst_next = &user_coarse->boundary_faces[face_i].params;
977 for (BC_Param *src = user_fine->boundary_faces[face_i].params; src; src = src->next) {
978 BC_Param *new_param;
979 ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
980 ierr = PetscStrallocpy(src->key, &new_param->key); CHKERRQ(ierr);
981 ierr = PetscStrallocpy(src->value, &new_param->value); CHKERRQ(ierr);
982 new_param->next = NULL;
983 *dst_next = new_param;
984 dst_next = &new_param->next;
985 }
986
987 // IMPORTANT: Do NOT create handler objects for coarser levels
988 // Handlers are only needed at finest level for timestepping Apply() calls
989 user_coarse->boundary_faces[face_i].handler = NULL;
990 }
991
992 // Propagate the particle inlet lookup fields to coarse levels as well.
993 user_coarse->inletFaceDefined = user_fine->inletFaceDefined;
994 user_coarse->identifiedInletBCFace = user_fine->identifiedInletBCFace;
995 }
996 }
997
998 LOG_ALLOW(GLOBAL, LOG_INFO, "BC configuration propagation complete.\n");
999
1000 PROFILE_FUNCTION_END;
1001 PetscFunctionReturn(0);
1002}
FreeBC_ParamList
void FreeBC_ParamList(BC_Param *head)
Frees an entire linked list of boundary-condition parameters.
Definition io.c:664
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:253
MGCtx::user
UserCtx * user
Definition variables.h:571
SimCtx::block_number
PetscInt block_number
Definition variables.h:790
BC_Param_s::next
struct BC_Param_s * next
Definition variables.h:339
BC_Param_s::key
char * key
Definition variables.h:337
SimCtx::usermg
UserMG usermg
Definition variables.h:852
BC_Param_s::value
char * value
Definition variables.h:338
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:370
UserMG::mglevels
PetscInt mglevels
Definition variables.h:578
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:581
BC_Param_s
A node in a linked list for storing key-value parameters from the bcs.dat file.
Definition variables.h:336
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:906
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:577
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundarySystem_ExecuteStep()

PetscErrorCode BoundarySystem_ExecuteStep ( UserCtx *  user)

Implementation of BoundarySystem_ExecuteStep().

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

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Boundaries.h.

See also
BoundarySystem_ExecuteStep()

Definition at line 1018 of file Boundaries.c.

1019{
1020 PetscErrorCode ierr;
1021 PetscFunctionBeginUser;
1022 PROFILE_FUNCTION_BEGIN;
1023
1024 LOG_ALLOW(LOCAL, LOG_DEBUG, "Starting.\n");
1025
1026 // =========================================================================
1027 // PRIORITY 0: INLETS
1028 // =========================================================================
1029
1030 PetscReal local_inflow_pre = 0.0;
1031 PetscReal local_inflow_post = 0.0;
1032 PetscReal global_inflow_pre = 0.0;
1033 PetscReal global_inflow_post = 0.0;
1034 PetscInt num_handlers[3] = {0,0,0};
1035
1036 LOG_ALLOW(LOCAL, LOG_TRACE, " (INLETS): Begin.\n");
1037
1038 // Phase 1: PreStep - Preparation (e.g., calculate profiles, read files)
1039 for (int i = 0; i < 6; i++) {
1040 BoundaryCondition *handler = user->boundary_faces[i].handler;
1041 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1042 if (!handler->PreStep) continue;
1043
1044 num_handlers[0]++;
1045 BCContext ctx = {
1046 .user = user,
1047 .face_id = (BCFace)i,
1048 .global_inflow_sum = NULL,
1049 .global_outflow_sum = NULL,
1050 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1051 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1052 };
1053
1054 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1055 ierr = handler->PreStep(handler, &ctx, &local_inflow_pre, NULL); CHKERRQ(ierr);
1056 }
1057
1058 // Optional: Global communication for PreStep (for debugging)
1059 if (local_inflow_pre != 0.0) {
1060 ierr = MPI_Allreduce(&local_inflow_pre, &global_inflow_pre, 1, MPIU_REAL,
1061 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1062 LOG_ALLOW(GLOBAL, LOG_TRACE, " PreStep predicted flux: %.6e\n", global_inflow_pre);
1063 }
1064
1065 // Phase 2: Apply - Set boundary conditions
1066 for (int i = 0; i < 6; i++) {
1067 BoundaryCondition *handler = user->boundary_faces[i].handler;
1068 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1069 if(!handler->Apply) continue; // For example Periodic BCs
1070
1071 num_handlers[1]++;
1072
1073 BCContext ctx = {
1074 .user = user,
1075 .face_id = (BCFace)i,
1076 .global_inflow_sum = NULL,
1077 .global_outflow_sum = NULL,
1078 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1079 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1080 };
1081
1082 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1083 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1084 }
1085
1086 // Phase 3: PostStep - Measure actual flux
1087 for (int i = 0; i < 6; i++) {
1088 BoundaryCondition *handler = user->boundary_faces[i].handler;
1089 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1090 if (!handler->PostStep) continue;
1091
1092 num_handlers[2]++;
1093
1094 BCContext ctx = {
1095 .user = user,
1096 .face_id = (BCFace)i,
1097 .global_inflow_sum = NULL,
1098 .global_outflow_sum = NULL,
1099 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1100 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1101 };
1102
1103 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1104 ierr = handler->PostStep(handler, &ctx, &local_inflow_post, NULL); CHKERRQ(ierr);
1105 }
1106
1107 // Phase 4: Global communication - Sum flux for other priorities to use
1108 ierr = MPI_Allreduce(&local_inflow_post, &global_inflow_post, 1, MPIU_REAL,
1109 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1110
1111 // Store for next priority levels
1112 user->simCtx->FluxInSum = global_inflow_post;
1113
1114 LOG_ALLOW(GLOBAL, LOG_INFO,
1115 " (INLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxInSum = %.6e\n",
1116 num_handlers[0],num_handlers[1],num_handlers[2], global_inflow_post);
1117
1118 // =========================================================================
1119 // PRIORITY 1: FARFIELD
1120 // =========================================================================
1121
1122 PetscReal local_farfield_in_pre = 0.0;
1123 PetscReal local_farfield_out_pre = 0.0;
1124 PetscReal local_farfield_in_post = 0.0;
1125 PetscReal local_farfield_out_post = 0.0;
1126 PetscReal global_farfield_in_pre = 0.0;
1127 PetscReal global_farfield_out_pre = 0.0;
1128 PetscReal global_farfield_in_post = 0.0;
1129 PetscReal global_farfield_out_post = 0.0;
1130 memset(num_handlers,0,sizeof(num_handlers));
1131
1132 LOG_ALLOW(LOCAL, LOG_TRACE, " (FARFIELD): Begin.\n");
1133
1134 // Phase 1: PreStep - Analyze flow direction, measure initial flux
1135 for (int i = 0; i < 6; i++) {
1136 BoundaryCondition *handler = user->boundary_faces[i].handler;
1137 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1138 if (!handler->PreStep) continue;
1139
1140 num_handlers[0]++;
1141 BCContext ctx = {
1142 .user = user,
1143 .face_id = (BCFace)i,
1144 .global_inflow_sum = &user->simCtx->FluxInSum, // Available from Priority 0
1145 .global_outflow_sum = NULL,
1146 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1147 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1148 };
1149
1150 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1151 ierr = handler->PreStep(handler, &ctx, &local_farfield_in_pre, &local_farfield_out_pre);
1152 CHKERRQ(ierr);
1153 }
1154
1155 // Phase 2: Global communication (optional, for debugging)
1156 if (local_farfield_in_pre != 0.0 || local_farfield_out_pre != 0.0) {
1157 ierr = MPI_Allreduce(&local_farfield_in_pre, &global_farfield_in_pre, 1, MPIU_REAL,
1158 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1159 ierr = MPI_Allreduce(&local_farfield_out_pre, &global_farfield_out_pre, 1, MPIU_REAL,
1160 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1161
1162 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1163 " Farfield pre-analysis: In=%.6e, Out=%.6e\n",
1164 global_farfield_in_pre, global_farfield_out_pre);
1165 }
1166
1167 // Phase 3: Apply - Set farfield boundary conditions
1168 for (int i = 0; i < 6; i++) {
1169 BoundaryCondition *handler = user->boundary_faces[i].handler;
1170 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1171 if(!handler->Apply) continue; // For example Periodic BCs
1172
1173 num_handlers[1]++;
1174
1175 BCContext ctx = {
1176 .user = user,
1177 .face_id = (BCFace)i,
1178 .global_inflow_sum = &user->simCtx->FluxInSum,
1179 .global_outflow_sum = NULL,
1180 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1181 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1182 };
1183
1184 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1185 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1186 }
1187
1188 // Phase 4: PostStep - Measure actual farfield fluxes
1189 for (int i = 0; i < 6; i++) {
1190 BoundaryCondition *handler = user->boundary_faces[i].handler;
1191 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1192 if (!handler->PostStep) continue;
1193
1194 num_handlers[2]++;
1195
1196 BCContext ctx = {
1197 .user = user,
1198 .face_id = (BCFace)i,
1199 .global_inflow_sum = &user->simCtx->FluxInSum,
1200 .global_outflow_sum = NULL,
1201 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1202 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1203 };
1204
1205 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1206 ierr = handler->PostStep(handler, &ctx, &local_farfield_in_post, &local_farfield_out_post);
1207 CHKERRQ(ierr);
1208 }
1209
1210 // Phase 5: Global communication - Store for outlet priority
1211 if (num_handlers > 0) {
1212 ierr = MPI_Allreduce(&local_farfield_in_post, &global_farfield_in_post, 1, MPIU_REAL,
1213 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1214 ierr = MPI_Allreduce(&local_farfield_out_post, &global_farfield_out_post, 1, MPIU_REAL,
1215 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1216
1217 // Store for outlet handlers to use
1218 user->simCtx->FarFluxInSum = global_farfield_in_post;
1219 user->simCtx->FarFluxOutSum = global_farfield_out_post;
1220
1221 LOG_ALLOW(GLOBAL, LOG_INFO,
1222 " (FARFIELD): %d Prestep(s), %d Application(s), %d Poststep(s) , InFlux=%.6e, OutFlux=%.6e\n",
1223 num_handlers[0],num_handlers[1],num_handlers[2], global_farfield_in_post, global_farfield_out_post);
1224 } else {
1225 // No farfield handlers - zero out the fluxes
1226 user->simCtx->FarFluxInSum = 0.0;
1227 user->simCtx->FarFluxOutSum = 0.0;
1228 }
1229
1230
1231 // =========================================================================
1232 // PRIORITY 2: WALLS
1233 // =========================================================================
1234
1235 memset(num_handlers,0,sizeof(num_handlers));
1236
1237 LOG_ALLOW(LOCAL, LOG_TRACE, " (WALLS): Begin.\n");
1238
1239 // Phase 1: PreStep - Preparation (usually no-op for walls)
1240 for (int i = 0; i < 6; i++) {
1241 BoundaryCondition *handler = user->boundary_faces[i].handler;
1242 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1243 if (!handler->PreStep) continue;
1244
1245 num_handlers[0]++;
1246 BCContext ctx = {
1247 .user = user,
1248 .face_id = (BCFace)i,
1249 .global_inflow_sum = &user->simCtx->FluxInSum,
1250 .global_outflow_sum = NULL,
1251 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1252 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1253 };
1254
1255 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1256 ierr = handler->PreStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1257 }
1258
1259 // No global communication needed for walls
1260
1261 // Phase 2: Apply - Set boundary conditions
1262 for (int i = 0; i < 6; i++) {
1263 BoundaryCondition *handler = user->boundary_faces[i].handler;
1264 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1265 if(!handler->Apply) continue; // For example Periodic BCs
1266
1267 num_handlers[1]++;
1268
1269 BCContext ctx = {
1270 .user = user,
1271 .face_id = (BCFace)i,
1272 .global_inflow_sum = &user->simCtx->FluxInSum,
1273 .global_outflow_sum = NULL,
1274 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1275 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1276 };
1277
1278 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1279 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1280 }
1281
1282 // Phase 3: PostStep - Post-application processing (usually no-op for walls)
1283 for (int i = 0; i < 6; i++) {
1284 BoundaryCondition *handler = user->boundary_faces[i].handler;
1285 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1286 if (!handler->PostStep) continue;
1287
1288 num_handlers[2]++;
1289
1290 BCContext ctx = {
1291 .user = user,
1292 .face_id = (BCFace)i,
1293 .global_inflow_sum = &user->simCtx->FluxInSum,
1294 .global_outflow_sum = NULL,
1295 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1296 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1297 };
1298
1299 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1300 ierr = handler->PostStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1301 }
1302
1303 // No global communication needed for walls
1304
1305 LOG_ALLOW(GLOBAL, LOG_INFO, " (WALLS): %d Prestep(s), %d Application(s), %d Poststep(s) applied.\n",
1306 num_handlers[0],num_handlers[1],num_handlers[2]);
1307
1308
1309 // =========================================================================
1310 // PRIORITY 3: OUTLETS
1311 // =========================================================================
1312
1313 PetscReal local_outflow_pre = 0.0;
1314 PetscReal local_outflow_post = 0.0;
1315 PetscReal global_outflow_pre = 0.0;
1316 PetscReal global_outflow_post = 0.0;
1317 memset(num_handlers,0,sizeof(num_handlers));
1318
1319 LOG_ALLOW(LOCAL, LOG_TRACE, " (OUTLETS): Begin.\n");
1320
1321 // Phase 1: PreStep - Measure uncorrected outflow (from ucat)
1322 for (int i = 0; i < 6; i++) {
1323 BoundaryCondition *handler = user->boundary_faces[i].handler;
1324 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1325 if (!handler->PreStep) continue;
1326
1327 num_handlers[0]++;
1328 BCContext ctx = {
1329 .user = user,
1330 .face_id = (BCFace)i,
1331 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1332 .global_outflow_sum = NULL,
1333 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1334 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1335 };
1336
1337 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1338 ierr = handler->PreStep(handler, &ctx, NULL, &local_outflow_pre); CHKERRQ(ierr);
1339 }
1340
1341 // Phase 2: Global communication - Get uncorrected outflow sum
1342 ierr = MPI_Allreduce(&local_outflow_pre, &global_outflow_pre, 1, MPIU_REAL,
1343 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1344
1345 // Calculate total inflow (inlet + farfield inflow)
1346 PetscReal total_inflow = user->simCtx->FluxInSum + user->simCtx->FarFluxInSum;
1347
1348 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1349 " Uncorrected outflow: %.6e, Total inflow: %.6e (Inlet: %.6e + Farfield: %.6e)\n",
1350 global_outflow_pre, total_inflow, user->simCtx->FluxInSum,
1351 user->simCtx->FarFluxInSum);
1352
1353 // Phase 3: Apply - Set corrected boundary conditions
1354 for (int i = 0; i < 6; i++) {
1355 BoundaryCondition *handler = user->boundary_faces[i].handler;
1356 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1357 if(!handler->Apply) continue; // For example Periodic BCs
1358
1359 num_handlers[1]++;
1360
1361 BCContext ctx = {
1362 .user = user,
1363 .face_id = (BCFace)i,
1364 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1365 .global_outflow_sum = &global_outflow_pre, // From PreStep above
1366 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1367 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1368 };
1369
1370 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1371 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1372 }
1373
1374 // Phase 4: PostStep - Measure corrected outflow (verification)
1375 for (int i = 0; i < 6; i++) {
1376 BoundaryCondition *handler = user->boundary_faces[i].handler;
1377 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1378 if (!handler->PostStep) continue;
1379
1380 num_handlers[2]++;
1381
1382 BCContext ctx = {
1383 .user = user,
1384 .face_id = (BCFace)i,
1385 .global_inflow_sum = &user->simCtx->FluxInSum,
1386 .global_outflow_sum = &global_outflow_pre,
1387 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1388 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1389 };
1390
1391 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1392 ierr = handler->PostStep(handler, &ctx, NULL, &local_outflow_post); CHKERRQ(ierr);
1393 }
1394
1395 // Phase 5: Global communication - Verify conservation
1396 ierr = MPI_Allreduce(&local_outflow_post, &global_outflow_post, 1, MPIU_REAL,
1397 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1398
1399 // Store for global reporting.
1400 user->simCtx->FluxOutSum = global_outflow_post;
1401
1402 // Conservation check (compare total outflow vs total inflow)
1403 PetscReal total_outflow = global_outflow_post + user->simCtx->FarFluxOutSum;
1404 PetscReal flux_error = PetscAbsReal(total_outflow - total_inflow);
1405 PetscReal relative_error = (total_inflow > 1e-16) ?
1406 flux_error / total_inflow : flux_error;
1407
1408 LOG_ALLOW(GLOBAL, LOG_INFO,
1409 " (OUTLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxOutSum = %.6e\n",
1410 num_handlers[0],num_handlers[1],num_handlers[2], global_outflow_post);
1411 LOG_ALLOW(GLOBAL, LOG_INFO,
1412 " Conservation: Total In=%.6e, Total Out=%.6e, Error=%.3e (%.2e)%%)\n",
1413 total_inflow, total_outflow, flux_error, relative_error * 100.0);
1414
1415 if (relative_error > 1e-6) {
1416 LOG_ALLOW(GLOBAL, LOG_WARNING,
1417 " WARNING: Large mass conservation error (%.2e%%)!\n",
1418 relative_error * 100.0);
1419 }
1420
1421
1422 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Complete.\n");
1423
1424 PROFILE_FUNCTION_END;
1425 PetscFunctionReturn(0);
1426}
BC_PRIORITY_OUTLET
@ BC_PRIORITY_OUTLET
Definition variables.h:328
BC_PRIORITY_FARFIELD
@ BC_PRIORITY_FARFIELD
Definition variables.h:326
BC_PRIORITY_WALL
@ BC_PRIORITY_WALL
Definition variables.h:327
BC_PRIORITY_INLET
@ BC_PRIORITY_INLET
Definition variables.h:325
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundarySystem_RefreshUbcs()

PetscErrorCode BoundarySystem_RefreshUbcs ( UserCtx *  user)

Internal helper implementation: BoundarySystem_RefreshUbcs().

(Private) A lightweight execution engine that calls the UpdateUbcs() method on all relevant handlers.

Local to this translation unit.

Definition at line 1440 of file Boundaries.c.

1441{
1442 PetscErrorCode ierr;
1443 PetscFunctionBeginUser;
1444
1445 LOG_ALLOW(GLOBAL, LOG_TRACE, "Refreshing `ubcs` targets for flow-dependent boundaries...\n");
1446
1447 // Loop through all 6 faces of the domain
1448 for (int i = 0; i < 6; i++) {
1449 BoundaryCondition *handler = user->boundary_faces[i].handler;
1450
1451 // THE FILTER:
1452 // This is the core logic. We only act if a handler exists for the face
1453 // AND that handler has explicitly implemented the `UpdateUbcs` method.
1454 if (handler && handler->UpdateUbcs) {
1455
1456 const char *face_name = BCFaceToString((BCFace)i);
1457 LOG_ALLOW(LOCAL, LOG_TRACE, " Calling UpdateUbcs() for handler on Face %s.\n", face_name);
1458
1459 // Prepare the context. For this refresh step, we don't need to pass flux sums.
1460 BCContext ctx = {
1461 .user = user,
1462 .face_id = (BCFace)i,
1463 .global_inflow_sum = NULL,
1464 .global_outflow_sum = NULL,
1465 .global_farfield_inflow_sum = NULL,
1466 .global_farfield_outflow_sum = NULL
1467 };
1468
1469 // Call the handler's specific UpdateUbcs function pointer.
1470 ierr = handler->UpdateUbcs(handler, &ctx); CHKERRQ(ierr);
1471 }
1472 }
1473
1474 PetscFunctionReturn(0);
1475}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BoundarySystem_Destroy()

PetscErrorCode BoundarySystem_Destroy ( UserCtx *  user)

Implementation of BoundarySystem_Destroy().

Cleans up and destroys all boundary system resources.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Boundaries.h.

See also
BoundarySystem_Destroy()

Definition at line 1490 of file Boundaries.c.

1491{
1492 PetscErrorCode ierr;
1493 PetscFunctionBeginUser;
1494
1495
1496
1497 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting destruction of all boundary handlers. \n");
1498
1499 for (int i = 0; i < 6; i++) {
1500 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
1501 const char *face_name = BCFaceToString(face_cfg->face_id);
1502
1503 // --- Step 1: Free the parameter linked list associated with this face ---
1504 if (face_cfg->params) {
1505 LOG_ALLOW(LOCAL, LOG_DEBUG, " Freeing parameter list for Face %d (%s). \n", i, face_name);
1506 FreeBC_ParamList(face_cfg->params);
1507 face_cfg->params = NULL; // Good practice to nullify dangling pointers
1508 }
1509
1510 // --- Step 2: Destroy the handler object itself ---
1511 if (face_cfg->handler) {
1512 const char *handler_name = BCHandlerTypeToString(face_cfg->handler->type);
1513 LOG_ALLOW(LOCAL, LOG_DEBUG, " Destroying handler '%s' on Face %d (%s).\n", handler_name, i, face_name);
1514
1515 // Call the handler's specific cleanup function first, if it exists.
1516 // This will free any memory stored in the handler's private `data` pointer.
1517 if (face_cfg->handler->Destroy) {
1518 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
1519 }
1520
1521 // Finally, free the generic BoundaryCondition object itself.
1522 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
1523 face_cfg->handler = NULL;
1524 }
1525 }
1526
1527 LOG_ALLOW(GLOBAL, LOG_INFO, "Destruction complete.\n");
1528 PetscFunctionReturn(0);
1529}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ TransferPeriodicFieldByDirection()

static PetscErrorCode TransferPeriodicFieldByDirection ( UserCtx *  user,
FieldId  field_id,
char  direction 
)
static

Copies one cell field's wrapped local values onto the owned periodic duplicate plane.

Handles scalar and three-component storage for one selected logical axis.

Definition at line 1537 of file Boundaries.c.

1538{
1539 PetscErrorCode ierr;
1540 DMDALocalInfo info = user->info;
1541 PetscInt xs = info.xs, xe = info.xs + info.xm;
1542 PetscInt ys = info.ys, ye = info.ys + info.ym;
1543 PetscInt zs = info.zs, ze = info.zs + info.zm;
1544 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1545
1546 FieldView field_view;
1547 DM dm;
1548 Vec global_vec;
1549 Vec local_vec;
1550 PetscInt dof;
1551
1552 PetscFunctionBeginUser;
1553 PetscCall(FieldGetView(user, field_id, &field_view));
1554 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_CELL_SYNC) != 0u,
1555 PETSC_COMM_SELF, PETSC_ERR_SUP,
1556 "Field '%s' is not registered for periodic cell synchronization.",
1557 field_view.descriptor->canonical_name);
1558 dm = field_view.dm;
1559 global_vec = field_view.global_vec;
1560 local_vec = field_view.local_vec;
1561 dof = field_view.descriptor->dof;
1562
1563 // --- Execute the copy logic based on DoF and Direction ---
1564 if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
1565 PetscReal ***g_array, ***l_array;
1566 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1567 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr); // Use Read for safety
1568
1569 switch (direction) {
1570 case 'i':
1571 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xs] = l_array[k][j][xs-2];
1572 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xe-1] = l_array[k][j][xe+1];
1573 break;
1574 case 'j':
1575 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ys][i] = l_array[k][ys-2][i];
1576 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ye-1][i] = l_array[k][ye+1][i];
1577 break;
1578 case 'k':
1579 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[zs][j][i] = l_array[zs-2][j][i];
1580 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[ze-1][j][i] = l_array[ze+1][j][i];
1581 break;
1582 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1583 }
1584 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1585 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1586
1587 } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
1588 Cmpnts ***g_array, ***l_array;
1589 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1590 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1591
1592 switch (direction) {
1593 case 'i':
1594 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xs] = l_array[k][j][xs-2];
1595 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xe-1] = l_array[k][j][xe+1];
1596 break;
1597 case 'j':
1598 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ys][i] = l_array[k][ys-2][i];
1599 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ye-1][i] = l_array[k][ye+1][i];
1600 break;
1601 case 'k':
1602 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[zs][j][i] = l_array[zs-2][j][i];
1603 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[ze-1][j][i] = l_array[ze+1][j][i];
1604 break;
1605 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1606 }
1607 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1608 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1609 }
1610
1611 PetscFunctionReturn(0);
1612}
FIELD_CAPABILITY_PERIODIC_CELL_SYNC
@ FIELD_CAPABILITY_PERIODIC_CELL_SYNC
Definition field_catalog.h:111
FieldDescriptor::capabilities
unsigned int capabilities
Definition field_catalog.h:129
FieldView::descriptor
const FieldDescriptor * descriptor
Definition field_catalog.h:136
FieldGetView
PetscErrorCode FieldGetView(UserCtx *user, FieldId field_id, FieldView *view)
Resolve the existing DM and global/local vectors for one field.
Definition field_catalog.c:259
FieldView::local_vec
Vec local_vec
Definition field_catalog.h:139
FieldView::dm
DM dm
Definition field_catalog.h:137
FieldDescriptor::dof
PetscInt dof
Definition field_catalog.h:124
FieldView::global_vec
Vec global_vec
Definition field_catalog.h:138
FieldDescriptor::canonical_name
const char * canonical_name
Definition field_catalog.h:121
FieldView
Non-owning runtime objects resolved for one field and UserCtx.
Definition field_catalog.h:135
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SynchronizePeriodicCellFields()

PetscErrorCode SynchronizePeriodicCellFields ( UserCtx *  user,
PetscInt  num_fields,
const FieldId  field_ids[] 
)

Implementation of SynchronizePeriodicCellFields().

Synchronizes periodic endpoint cells for a list of cell-centered fields.

Full API contract is documented with the header declaration in include/Boundaries.h.

Definition at line 1621 of file Boundaries.c.

1622{
1623 PetscErrorCode ierr;
1624 PetscBool periodic_i;
1625 PetscBool periodic_j;
1626 PetscBool periodic_k;
1627
1628 PetscFunctionBeginUser;
1629
1630 PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1631 "Number of cell fields cannot be negative.");
1632 if (num_fields == 0) PetscFunctionReturn(0);
1633 PetscCheck(field_ids != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
1634 "Cell field-ID array cannot be NULL.");
1635
1636 periodic_i =
1637 user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
1638 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC;
1639 periodic_j =
1640 user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
1641 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC;
1642 periodic_k =
1643 user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
1644 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC;
1645
1646 if (!periodic_i && !periodic_j && !periodic_k) PetscFunctionReturn(0);
1647
1648 for (PetscInt field = 0; field < num_fields; field++) {
1649 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1650 }
1651
1652 if (periodic_i) {
1653 for (PetscInt field = 0; field < num_fields; field++) {
1654 ierr = TransferPeriodicFieldByDirection(user, field_ids[field], 'i'); CHKERRQ(ierr);
1655 }
1656 for (PetscInt field = 0; field < num_fields; field++) {
1657 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1658 }
1659 }
1660
1661 if (periodic_j) {
1662 for (PetscInt field = 0; field < num_fields; field++) {
1663 ierr = TransferPeriodicFieldByDirection(user, field_ids[field], 'j'); CHKERRQ(ierr);
1664 }
1665 for (PetscInt field = 0; field < num_fields; field++) {
1666 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1667 }
1668 }
1669
1670 if (periodic_k) {
1671 for (PetscInt field = 0; field < num_fields; field++) {
1672 ierr = TransferPeriodicFieldByDirection(user, field_ids[field], 'k'); CHKERRQ(ierr);
1673 }
1674 for (PetscInt field = 0; field < num_fields; field++) {
1675 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1676 }
1677 }
1678
1679 PetscFunctionReturn(0);
1680}
TransferPeriodicFieldByDirection
static PetscErrorCode TransferPeriodicFieldByDirection(UserCtx *user, FieldId field_id, char direction)
Copies one cell field's wrapped local values onto the owned periodic duplicate plane.
Definition Boundaries.c:1537
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1838
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetPersistentFaceField()

static PetscErrorCode GetPersistentFaceField ( UserCtx *  user,
FieldId  field_id,
char  face_direction,
DM *  dm,
Vec *  global_vec,
Vec *  local_vec,
PetscInt *  dof 
)
static

Resolves one registered persistent single-face-family field.

Definition at line 1687 of file Boundaries.c.

1691{
1692 FieldView field_view;
1693 FieldLayout expected_layout;
1694
1695 PetscFunctionBeginUser;
1696 PetscCheck(face_direction == 'i' || face_direction == 'j' || face_direction == 'k',
1697 PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1698 "Invalid face direction '%c'; expected 'i', 'j', or 'k'.", face_direction);
1699
1700 expected_layout = face_direction == 'i' ? FIELD_LAYOUT_I_FACE :
1701 (face_direction == 'j' ? FIELD_LAYOUT_J_FACE : FIELD_LAYOUT_K_FACE);
1702 PetscCall(FieldGetView(user, field_id, &field_view));
1703 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_FACE_SYNC) != 0u &&
1704 field_view.descriptor->layout == expected_layout,
1705 PETSC_COMM_SELF, PETSC_ERR_SUP,
1706 "Field '%s' is not registered for %c-face periodic synchronization.",
1707 field_view.descriptor->canonical_name, face_direction);
1708
1709 *dm = field_view.dm;
1710 *global_vec = field_view.global_vec;
1711 *local_vec = field_view.local_vec;
1712 *dof = field_view.descriptor->dof;
1713 PetscFunctionReturn(0);
1714}
FieldDescriptor::layout
FieldLayout layout
Definition field_catalog.h:126
FIELD_CAPABILITY_PERIODIC_FACE_SYNC
@ FIELD_CAPABILITY_PERIODIC_FACE_SYNC
Definition field_catalog.h:112
FieldLayout
FieldLayout
Logical storage topology of a field.
Definition field_catalog.h:71
FIELD_LAYOUT_K_FACE
@ FIELD_LAYOUT_K_FACE
Definition field_catalog.h:76
FIELD_LAYOUT_I_FACE
@ FIELD_LAYOUT_I_FACE
Definition field_catalog.h:74
FIELD_LAYOUT_J_FACE
@ FIELD_LAYOUT_J_FACE
Definition field_catalog.h:75
Here is the call graph for this function:
Here is the caller graph for this function:

◆ IsFaceCenterCoordinateField()

static PetscErrorCode IsFaceCenterCoordinateField ( FieldId  field_id,
PetscBool *  is_coordinate 
)
static

Returns whether a registered face field stores physical coordinates.

Definition at line 1719 of file Boundaries.c.

1720{
1721 const FieldDescriptor *descriptor = NULL;
1722
1723 PetscFunctionBeginUser;
1724 PetscCheck(is_coordinate != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
1725 "Face-coordinate result cannot be NULL.");
1726 PetscCall(FieldGetDescriptor(field_id, &descriptor));
1727 *is_coordinate = (PetscBool)((descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_GEOMETRY_SHIFT) != 0u);
1728 PetscFunctionReturn(0);
1729}
FIELD_CAPABILITY_PERIODIC_GEOMETRY_SHIFT
@ FIELD_CAPABILITY_PERIODIC_GEOMETRY_SHIFT
Definition field_catalog.h:114
FieldGetDescriptor
PetscErrorCode FieldGetDescriptor(FieldId field_id, const FieldDescriptor **descriptor)
Return immutable metadata for a valid field identifier.
Definition field_catalog.c:178
FieldDescriptor
Immutable metadata for one field identity.
Definition field_catalog.h:119
Here is the call graph for this function:
Here is the caller graph for this function:

◆ TranslatePeriodicFaceCenterGhosts()

static PetscErrorCode TranslatePeriodicFaceCenterGhosts ( UserCtx *  user,
Vec  local_vec 
)
static

Applies geometric translations to wrapped face-center ghost coordinates.

Definition at line 1734 of file Boundaries.c.

1735{
1736 DMDALocalInfo info;
1737 Cmpnts ***array;
1738 const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
1739 const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
1740
1741 PetscFunctionBeginUser;
1742 PetscCall(DMDAGetLocalInfo(user->fda, &info));
1743 PetscCall(DMDAVecGetArray(user->fda, local_vec, &array));
1744
1745 for (PetscInt axis = 0; axis < 3; axis++) {
1746 const PetscBool active =
1747 user->boundary_faces[negative_faces[axis]].mathematical_type == PERIODIC ||
1748 user->boundary_faces[positive_faces[axis]].mathematical_type == PERIODIC;
1749 Cmpnts translation;
1750
1751 if (!active) continue;
1752 PetscCheck(user->periodic_translation_valid[axis], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
1753 "Periodic face-center synchronization requires validated %c-direction geometry.",
1754 "XYZ"[axis]);
1755 translation = user->periodic_translation[axis];
1756
1757 for (PetscInt k = info.gzs; k < info.gzs + info.gzm; k++) {
1758 for (PetscInt j = info.gys; j < info.gys + info.gym; j++) {
1759 for (PetscInt i = info.gxs; i < info.gxs + info.gxm; i++) {
1760 PetscReal scale = 0.0;
1761 const PetscInt index = axis == 0 ? i : (axis == 1 ? j : k);
1762 const PetscInt size = axis == 0 ? info.mx : (axis == 1 ? info.my : info.mz);
1763 if (index < 0) scale = -1.0;
1764 else if (index >= size) scale = 1.0;
1765 if (scale == 0.0) continue;
1766 array[k][j][i].x += scale * translation.x;
1767 array[k][j][i].y += scale * translation.y;
1768 array[k][j][i].z += scale * translation.z;
1769 }
1770 }
1771 }
1772 }
1773
1774 PetscCall(DMDAVecRestoreArray(user->fda, local_vec, &array));
1775 PetscFunctionReturn(0);
1776}
Cmpnts::z
PetscScalar z
Definition variables.h:103
Cmpnts::y
PetscScalar y
Definition variables.h:103
UserCtx::periodic_translation
Cmpnts periodic_translation[3]
Definition variables.h:927
UserCtx::periodic_translation_valid
PetscBool periodic_translation_valid[3]
Definition variables.h:928
Here is the caller graph for this function:

◆ TransferPeriodicFaceFieldByDirection()

static PetscErrorCode TransferPeriodicFaceFieldByDirection ( UserCtx *  user,
FieldId  field_id,
char  face_direction,
char  periodic_direction 
)
static

Transfers one registered face-family field along one periodic axis.

Definition at line 1781 of file Boundaries.c.

1783{
1784 PetscErrorCode ierr;
1785 DMDALocalInfo info = user->info;
1786 PetscInt xs = info.xs, xe = info.xs + info.xm;
1787 PetscInt ys = info.ys, ye = info.ys + info.ym;
1788 PetscInt zs = info.zs, ze = info.zs + info.zm;
1789 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1790 DM dm;
1791 Vec global_vec, local_vec;
1792 PetscInt dof;
1793
1794 PetscFunctionBeginUser;
1795 PetscCheck(periodic_direction == 'i' || periodic_direction == 'j' || periodic_direction == 'k',
1796 PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1797 "Invalid periodic direction '%c'; expected 'i', 'j', or 'k'.", periodic_direction);
1798 PetscCall(GetPersistentFaceField(user, field_id, face_direction, &dm, &global_vec, &local_vec, &dof));
1799
1800 if (dof == 1) {
1801 PetscReal ***global_array, ***local_array;
1802 ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1803 ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1804
1805 if (periodic_direction == 'i') {
1806 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
1807 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
1808 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
1809 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
1810 } else if (periodic_direction == 'j') {
1811 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
1812 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
1813 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
1814 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
1815 } else {
1816 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
1817 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
1818 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
1819 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
1820 }
1821 ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1822 ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1823 } else {
1824 Cmpnts ***global_array, ***local_array;
1825 ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1826 ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1827
1828 if (periodic_direction == 'i') {
1829 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
1830 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
1831 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
1832 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
1833 } else if (periodic_direction == 'j') {
1834 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
1835 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
1836 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
1837 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
1838 } else {
1839 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
1840 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
1841 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
1842 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
1843 }
1844 ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1845 ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1846 }
1847
1848 PetscFunctionReturn(0);
1849}
GetPersistentFaceField
static PetscErrorCode GetPersistentFaceField(UserCtx *user, FieldId field_id, char face_direction, DM *dm, Vec *global_vec, Vec *local_vec, PetscInt *dof)
Resolves one registered persistent single-face-family field.
Definition Boundaries.c:1687
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SynchronizePeriodicFaceFields()

PetscErrorCode SynchronizePeriodicFaceFields ( UserCtx *  user,
char  face_direction,
PetscInt  num_fields,
const FieldId  field_ids[] 
)

Synchronizes persistent fields belonging to one face family.

The function performs deterministic I/J/K directional passes with an intermediate ghost refresh after each active periodic direction. It updates persistent global seam/dummy values only; face-specific local stencil repair remains a separate operation.

Parameters
userThe main UserCtx struct.
face_directionFace family shared by every field (‘'i’,'j', or'k'`).
[in]num_fieldsCount of registered face fields.
field_idsRegistered persistent face-field identities.
Returns
PetscErrorCode 0 on success.

Definition at line 1855 of file Boundaries.c.

1857{
1858 PetscErrorCode ierr;
1859 const char periodic_directions[3] = {'i', 'j', 'k'};
1860 const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
1861 const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
1862
1863 PetscFunctionBeginUser;
1864 PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1865 "Number of face fields cannot be negative.");
1866 if (num_fields == 0) PetscFunctionReturn(0);
1867 PetscCheck(field_ids != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
1868 "Face field-ID array cannot be NULL.");
1869
1870 for (PetscInt field = 0; field < num_fields; field++) {
1871 DM dm;
1872 Vec global_vec, local_vec;
1873 PetscInt dof;
1874 PetscBool is_coordinate;
1875 PetscCall(GetPersistentFaceField(user, field_ids[field], face_direction,
1876 &dm, &global_vec, &local_vec, &dof));
1877 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1878 ierr = IsFaceCenterCoordinateField(field_ids[field], &is_coordinate); CHKERRQ(ierr);
1879 if (is_coordinate) {
1880 ierr = TranslatePeriodicFaceCenterGhosts(user, local_vec); CHKERRQ(ierr);
1881 }
1882 }
1883
1884 for (PetscInt direction = 0; direction < 3; direction++) {
1885 const PetscBool active =
1886 user->boundary_faces[negative_faces[direction]].mathematical_type == PERIODIC ||
1887 user->boundary_faces[positive_faces[direction]].mathematical_type == PERIODIC;
1888 if (!active) continue;
1889
1890 for (PetscInt field = 0; field < num_fields; field++) {
1891 ierr = TransferPeriodicFaceFieldByDirection(user, field_ids[field], face_direction,
1892 periodic_directions[direction]); CHKERRQ(ierr);
1893 }
1894 for (PetscInt field = 0; field < num_fields; field++) {
1895 DM dm;
1896 Vec global_vec, local_vec;
1897 PetscInt dof;
1898 PetscBool is_coordinate;
1899 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1900 ierr = IsFaceCenterCoordinateField(field_ids[field], &is_coordinate); CHKERRQ(ierr);
1901 if (is_coordinate) {
1902 PetscCall(GetPersistentFaceField(user, field_ids[field], face_direction,
1903 &dm, &global_vec, &local_vec, &dof));
1904 ierr = TranslatePeriodicFaceCenterGhosts(user, local_vec); CHKERRQ(ierr);
1905 }
1906 }
1907 }
1908
1909 PetscFunctionReturn(0);
1910}
TranslatePeriodicFaceCenterGhosts
static PetscErrorCode TranslatePeriodicFaceCenterGhosts(UserCtx *user, Vec local_vec)
Applies geometric translations to wrapped face-center ghost coordinates.
Definition Boundaries.c:1734
IsFaceCenterCoordinateField
static PetscErrorCode IsFaceCenterCoordinateField(FieldId field_id, PetscBool *is_coordinate)
Returns whether a registered face field stores physical coordinates.
Definition Boundaries.c:1719
TransferPeriodicFaceFieldByDirection
static PetscErrorCode TransferPeriodicFaceFieldByDirection(UserCtx *user, FieldId field_id, char face_direction, char periodic_direction)
Transfers one registered face-family field along one periodic axis.
Definition Boundaries.c:1781
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetPersistentStaggeredField()

static PetscErrorCode GetPersistentStaggeredField ( UserCtx *  user,
FieldId  field_id,
DM *  dm,
Vec *  global_vec,
Vec *  local_vec 
)
static

Resolves one registered persistent component-staggered field.

Definition at line 1917 of file Boundaries.c.

1919{
1920 FieldView field_view;
1921
1922 PetscFunctionBeginUser;
1923 PetscCall(FieldGetView(user, field_id, &field_view));
1924 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_STAGGERED_SYNC) != 0u &&
1925 field_view.descriptor->layout == FIELD_LAYOUT_COMPONENT_STAGGERED,
1926 PETSC_COMM_SELF, PETSC_ERR_SUP,
1927 "Field '%s' is not registered for component-staggered periodic synchronization.",
1928 field_view.descriptor->canonical_name);
1929 *dm = field_view.dm;
1930 *global_vec = field_view.global_vec;
1931 *local_vec = field_view.local_vec;
1932 PetscFunctionReturn(0);
1933}
FIELD_CAPABILITY_PERIODIC_STAGGERED_SYNC
@ FIELD_CAPABILITY_PERIODIC_STAGGERED_SYNC
Definition field_catalog.h:113
FIELD_LAYOUT_COMPONENT_STAGGERED
@ FIELD_LAYOUT_COMPONENT_STAGGERED
Definition field_catalog.h:77
Here is the call graph for this function:
Here is the caller graph for this function:

◆ TransferPeriodicStaggeredFieldByDirection()

static PetscErrorCode TransferPeriodicStaggeredFieldByDirection ( UserCtx *  user,
FieldId  field_id,
char  periodic_direction 
)
static

Transfers one component-staggered field along one periodic axis.

Definition at line 1940 of file Boundaries.c.

1942{
1943 PetscErrorCode ierr;
1944 DMDALocalInfo info = user->info;
1945 PetscInt xs = info.xs, xe = info.xs + info.xm;
1946 PetscInt ys = info.ys, ye = info.ys + info.ym;
1947 PetscInt zs = info.zs, ze = info.zs + info.zm;
1948 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1949 DM dm;
1950 Vec global_vec, local_vec;
1951 Cmpnts ***global_array, ***local_array;
1952
1953 PetscFunctionBeginUser;
1954 PetscCheck(periodic_direction == 'i' || periodic_direction == 'j' || periodic_direction == 'k',
1955 PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1956 "Invalid periodic direction '%c'; expected 'i', 'j', or 'k'.", periodic_direction);
1957 PetscCall(GetPersistentStaggeredField(user, field_id, &dm, &global_vec, &local_vec));
1958
1959 ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1960 ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1961
1962 if (periodic_direction == 'i') {
1963 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
1964 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
1965 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
1966 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
1967 } else if (periodic_direction == 'j') {
1968 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
1969 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
1970 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
1971 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
1972 } else {
1973 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
1974 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
1975 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
1976 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
1977 }
1978
1979 ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1980 ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1981 PetscFunctionReturn(0);
1982}
GetPersistentStaggeredField
static PetscErrorCode GetPersistentStaggeredField(UserCtx *user, FieldId field_id, DM *dm, Vec *global_vec, Vec *local_vec)
Resolves one registered persistent component-staggered field.
Definition Boundaries.c:1917
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SynchronizePeriodicStaggeredFields()

PetscErrorCode SynchronizePeriodicStaggeredFields ( UserCtx *  user,
PetscInt  num_fields,
const FieldId  field_ids[] 
)

Implementation of SynchronizePeriodicStaggeredFields().

Synchronizes persistent component-staggered vector fields.

Definition at line 1989 of file Boundaries.c.

1991{
1992 PetscErrorCode ierr;
1993 const char periodic_directions[3] = {'i', 'j', 'k'};
1994 const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
1995 const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
1996
1997 PetscFunctionBeginUser;
1998 PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1999 "Number of staggered fields cannot be negative.");
2000 if (num_fields == 0) PetscFunctionReturn(0);
2001 PetscCheck(field_ids != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
2002 "Staggered field-ID array cannot be NULL.");
2003
2004 for (PetscInt field = 0; field < num_fields; field++) {
2005 DM dm;
2006 Vec global_vec, local_vec;
2007 PetscCall(GetPersistentStaggeredField(user, field_ids[field], &dm, &global_vec, &local_vec));
2008 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
2009 }
2010
2011 for (PetscInt direction = 0; direction < 3; direction++) {
2012 const PetscBool active =
2013 user->boundary_faces[negative_faces[direction]].mathematical_type == PERIODIC ||
2014 user->boundary_faces[positive_faces[direction]].mathematical_type == PERIODIC;
2015 if (!active) continue;
2016
2017 for (PetscInt field = 0; field < num_fields; field++) {
2018 ierr = TransferPeriodicStaggeredFieldByDirection(user, field_ids[field],
2019 periodic_directions[direction]); CHKERRQ(ierr);
2020 }
2021 for (PetscInt field = 0; field < num_fields; field++) {
2022 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
2023 }
2024 }
2025
2026 PetscFunctionReturn(0);
2027}
TransferPeriodicStaggeredFieldByDirection
static PetscErrorCode TransferPeriodicStaggeredFieldByDirection(UserCtx *user, FieldId field_id, char periodic_direction)
Transfers one component-staggered field along one periodic axis.
Definition Boundaries.c:1940
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PreparePeriodicQuickStencilFields()

PetscErrorCode PreparePeriodicQuickStencilFields ( UserCtx *  user,
Vec  local_vector_field,
Vec  local_scalar_field 
)

Implementation of PreparePeriodicQuickStencilFields().

Repairs the outer adjacent periodic ghosts used by QUICK cell stencils.

Definition at line 2034 of file Boundaries.c.

2036{
2037 DMDALocalInfo info = user->info;
2038 Cmpnts ***vector_array;
2039 PetscReal ***scalar_array;
2040 const PetscInt xs = info.xs, xe = info.xs + info.xm;
2041 const PetscInt ys = info.ys, ye = info.ys + info.ym;
2042 const PetscInt zs = info.zs, ze = info.zs + info.zm;
2043 const PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
2044 const PetscInt gys = info.gys, gye = info.gys + info.gym;
2045 const PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
2046
2047 PetscFunctionBeginUser;
2048 PetscCheck(local_vector_field && local_scalar_field, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
2049 "QUICK stencil repair requires both local vector and scalar fields.");
2050 PetscCall(DMDAVecGetArray(user->fda, local_vector_field, &vector_array));
2051 PetscCall(DMDAVecGetArray(user->da, local_scalar_field, &scalar_array));
2052
2053 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
2054 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++) {
2055 vector_array[k][j][-1] = vector_array[k][j][-3];
2056 scalar_array[k][j][-1] = scalar_array[k][j][-3];
2057 }
2058 }
2059 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == info.mx) {
2060 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++) {
2061 vector_array[k][j][info.mx] = vector_array[k][j][info.mx + 2];
2062 scalar_array[k][j][info.mx] = scalar_array[k][j][info.mx + 2];
2063 }
2064 }
2065 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
2066 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++) {
2067 vector_array[k][-1][i] = vector_array[k][-3][i];
2068 scalar_array[k][-1][i] = scalar_array[k][-3][i];
2069 }
2070 }
2071 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == info.my) {
2072 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++) {
2073 vector_array[k][info.my][i] = vector_array[k][info.my + 2][i];
2074 scalar_array[k][info.my][i] = scalar_array[k][info.my + 2][i];
2075 }
2076 }
2077 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
2078 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++) {
2079 vector_array[-1][j][i] = vector_array[-3][j][i];
2080 scalar_array[-1][j][i] = scalar_array[-3][j][i];
2081 }
2082 }
2083 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == info.mz) {
2084 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++) {
2085 vector_array[info.mz][j][i] = vector_array[info.mz + 2][j][i];
2086 scalar_array[info.mz][j][i] = scalar_array[info.mz + 2][j][i];
2087 }
2088 }
2089
2090 PetscCall(DMDAVecRestoreArray(user->da, local_scalar_field, &scalar_array));
2091 PetscCall(DMDAVecRestoreArray(user->fda, local_vector_field, &vector_array));
2092 PetscFunctionReturn(0);
2093}
Here is the caller graph for this function:

◆ SynchronizePeriodicLocalStaggeredField()

PetscErrorCode SynchronizePeriodicLocalStaggeredField ( UserCtx *  user,
Vec  local_field 
)

Implementation of SynchronizePeriodicLocalStaggeredField().

Synchronizes one local-only component-staggered periodic work field.

Definition at line 2100 of file Boundaries.c.

2101{
2102 DMDALocalInfo info = user->info;
2103 Cmpnts ***array;
2104 const PetscInt xs = info.xs, xe = info.xs + info.xm;
2105 const PetscInt ys = info.ys, ye = info.ys + info.ym;
2106 const PetscInt zs = info.zs, ze = info.zs + info.zm;
2107
2108 PetscFunctionBeginUser;
2109 PetscCheck(local_field, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
2110 "Local staggered field cannot be NULL.");
2111 PetscCall(DMLocalToLocalBegin(user->fda, local_field, INSERT_VALUES, local_field));
2112 PetscCall(DMLocalToLocalEnd(user->fda, local_field, INSERT_VALUES, local_field));
2113 PetscCall(DMDAVecGetArray(user->fda, local_field, &array));
2114
2115 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
2116 for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) array[k][j][0].x = array[k][j][-2].x;
2117 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == info.mx)
2118 for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) array[k][j][info.mx - 1].x = array[k][j][info.mx + 1].x;
2119 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
2120 for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) array[k][0][i].y = array[k][-2][i].y;
2121 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == info.my)
2122 for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) array[k][info.my - 1][i].y = array[k][info.my + 1][i].y;
2123 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
2124 for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) array[0][j][i].z = array[-2][j][i].z;
2125 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == info.mz)
2126 for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) array[info.mz - 1][j][i].z = array[info.mz + 1][j][i].z;
2127
2128 PetscCall(DMDAVecRestoreArray(user->fda, local_field, &array));
2129 PetscCall(DMLocalToLocalBegin(user->fda, local_field, INSERT_VALUES, local_field));
2130 PetscCall(DMLocalToLocalEnd(user->fda, local_field, INSERT_VALUES, local_field));
2131 PetscFunctionReturn(0);
2132}
Here is the caller graph for this function:

◆ ApplyMetricsPeriodicBCs()

PetscErrorCode ApplyMetricsPeriodicBCs ( UserCtx *  user)

Internal helper implementation: ApplyMetricsPeriodicBCs().

(Orchestrator) Updates all metric-related fields in the local ghost cell regions for periodic boundaries.

Local to this translation unit.

Definition at line 2140 of file Boundaries.c.

2141{
2142 PetscErrorCode ierr;
2143 PetscFunctionBeginUser;
2144 PROFILE_FUNCTION_BEGIN;
2145
2146 const FieldId cell_fields[] = {FIELD_ID_AJ};
2147 const FieldId i_face_fields[] = {FIELD_ID_CENTX, FIELD_ID_CSI, FIELD_ID_ICSI,
2148 FIELD_ID_IETA, FIELD_ID_IZET, FIELD_ID_IAJ};
2149 const FieldId j_face_fields[] = {FIELD_ID_CENTY, FIELD_ID_ETA, FIELD_ID_JCSI,
2150 FIELD_ID_JETA, FIELD_ID_JZET, FIELD_ID_JAJ};
2151 const FieldId k_face_fields[] = {FIELD_ID_CENTZ, FIELD_ID_ZET, FIELD_ID_KCSI,
2152 FIELD_ID_KETA, FIELD_ID_KZET, FIELD_ID_KAJ};
2153
2154 ierr = SynchronizePeriodicCellFields(user, 1, cell_fields); CHKERRQ(ierr);
2155 ierr = SynchronizePeriodicFaceFields(user, 'i', 6, i_face_fields); CHKERRQ(ierr);
2156 ierr = SynchronizePeriodicFaceFields(user, 'j', 6, j_face_fields); CHKERRQ(ierr);
2157 ierr = SynchronizePeriodicFaceFields(user, 'k', 6, k_face_fields); CHKERRQ(ierr);
2158
2159 PROFILE_FUNCTION_END;
2160 PetscFunctionReturn(0);
2161}
SynchronizePeriodicFaceFields
PetscErrorCode SynchronizePeriodicFaceFields(UserCtx *user, char face_direction, PetscInt num_fields, const FieldId field_ids[])
Synchronizes persistent fields belonging to one face family.
Definition Boundaries.c:1855
SynchronizePeriodicCellFields
PetscErrorCode SynchronizePeriodicCellFields(UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
Implementation of SynchronizePeriodicCellFields().
Definition Boundaries.c:1621
FieldId
FieldId
Compile-time identity for a catalogued Eulerian field.
Definition field_catalog.h:23
FIELD_ID_JETA
@ FIELD_ID_JETA
Definition field_catalog.h:49
FIELD_ID_CENTZ
@ FIELD_ID_CENTZ
Definition field_catalog.h:44
FIELD_ID_CSI
@ FIELD_ID_CSI
Definition field_catalog.h:35
FIELD_ID_IAJ
@ FIELD_ID_IAJ
Definition field_catalog.h:54
FIELD_ID_KETA
@ FIELD_ID_KETA
Definition field_catalog.h:52
FIELD_ID_JAJ
@ FIELD_ID_JAJ
Definition field_catalog.h:55
FIELD_ID_KAJ
@ FIELD_ID_KAJ
Definition field_catalog.h:56
FIELD_ID_AJ
@ FIELD_ID_AJ
Definition field_catalog.h:39
FIELD_ID_CENTY
@ FIELD_ID_CENTY
Definition field_catalog.h:43
FIELD_ID_KZET
@ FIELD_ID_KZET
Definition field_catalog.h:53
FIELD_ID_IETA
@ FIELD_ID_IETA
Definition field_catalog.h:46
FIELD_ID_ICSI
@ FIELD_ID_ICSI
Definition field_catalog.h:45
FIELD_ID_ETA
@ FIELD_ID_ETA
Definition field_catalog.h:36
FIELD_ID_JCSI
@ FIELD_ID_JCSI
Definition field_catalog.h:48
FIELD_ID_JZET
@ FIELD_ID_JZET
Definition field_catalog.h:50
FIELD_ID_IZET
@ FIELD_ID_IZET
Definition field_catalog.h:47
FIELD_ID_ZET
@ FIELD_ID_ZET
Definition field_catalog.h:37
FIELD_ID_KCSI
@ FIELD_ID_KCSI
Definition field_catalog.h:51
FIELD_ID_CENTX
@ FIELD_ID_CENTX
Definition field_catalog.h:42
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ApplyPeriodicBCs()

PetscErrorCode ApplyPeriodicBCs ( UserCtx *  user)

Internal helper implementation: ApplyPeriodicBCs().

Applies periodic boundary conditions by copying data across domain boundaries for all relevant fields.

Local to this translation unit.

Definition at line 2169 of file Boundaries.c.

2170{
2171 PetscErrorCode ierr;
2172 PetscBool is_any_periodic = PETSC_FALSE;
2173
2174 PetscFunctionBeginUser;
2175
2176 PROFILE_FUNCTION_BEGIN;
2177
2178 for (int i = 0; i < 6; i++) {
2179 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
2180 is_any_periodic = PETSC_TRUE;
2181 break;
2182 }
2183 }
2184
2185 if (!is_any_periodic) {
2186 LOG_ALLOW(GLOBAL,LOG_TRACE, "No periodic boundaries defined; skipping ApplyPeriodicBCs.\n");
2187 PROFILE_FUNCTION_END;
2188 PetscFunctionReturn(0);
2189 }
2190
2191 LOG_ALLOW(GLOBAL, LOG_TRACE, "Applying periodic boundary conditions for all fields.\n");
2192
2193 // STEP 1: Synchronize periodic cell-centered fields in deterministic direction order.
2194 const FieldId cell_fields[] = {FIELD_ID_UCAT, FIELD_ID_P, FIELD_ID_NVERT};
2195 ierr = SynchronizePeriodicCellFields(user, 3, cell_fields); CHKERRQ(ierr);
2196
2197 /* A future temperature field must be catalogued before requesting its typed ghost update. */
2198
2199 // STEP 2: Synchronize persistent staggered endpoints and repair local
2200 // component-normal ghosts through UpdateLocalGhosts().
2201 const FieldId staggered_fields[] = {FIELD_ID_UCONT};
2202 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
2203
2204 // FUTURE EXTENSION: Add new cell fields through SynchronizePeriodicCellFields().
2205 /*
2206 if (user->solve_temperature) {
2207 const char *temperature_field[] = {"Temperature"};
2208 ierr = SynchronizePeriodicCellFields(user, 1, temperature_field); CHKERRQ(ierr);
2209 }
2210 */
2211
2212 PROFILE_FUNCTION_END;
2213 PetscFunctionReturn(0);
2214}
SynchronizePeriodicStaggeredFields
PetscErrorCode SynchronizePeriodicStaggeredFields(UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
Implementation of SynchronizePeriodicStaggeredFields().
Definition Boundaries.c:1989
FIELD_ID_NVERT
@ FIELD_ID_NVERT
Definition field_catalog.h:38
FIELD_ID_UCAT
@ FIELD_ID_UCAT
Definition field_catalog.h:26
FIELD_ID_UCONT
@ FIELD_ID_UCONT
Definition field_catalog.h:27
FIELD_ID_P
@ FIELD_ID_P
Definition field_catalog.h:30
Here is the call graph for this function:
Here is the caller graph for this function:

◆ UpdateDummyCells()

PetscErrorCode UpdateDummyCells ( UserCtx *  user)

Internal helper implementation: UpdateDummyCells().

Updates the dummy cells (ghost nodes) on the faces of the local domain for NON-PERIODIC boundaries.

Local to this translation unit.

Definition at line 2222 of file Boundaries.c.

2223{
2224 PetscErrorCode ierr;
2225 DM fda = user->fda;
2226 DMDALocalInfo info = user->info;
2227 PetscInt xs = info.xs, xe = info.xs + info.xm;
2228 PetscInt ys = info.ys, ye = info.ys + info.ym;
2229 PetscInt zs = info.zs, ze = info.zs + info.zm;
2230 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2231
2232 // --- Calculate shrunken loop ranges to avoid edges and corners ---
2233 PetscInt lxs = xs, lxe = xe;
2234 PetscInt lys = ys, lye = ye;
2235 PetscInt lzs = zs, lze = ze;
2236
2237 if (xs == 0) lxs = xs + 1;
2238 if (ys == 0) lys = ys + 1;
2239 if (zs == 0) lzs = zs + 1;
2240
2241 if (xe == mx) lxe = xe - 1;
2242 if (ye == my) lye = ye - 1;
2243 if (ze == mz) lze = ze - 1;
2244
2245 Cmpnts ***ucat, ***ubcs;
2246 PetscFunctionBeginUser;
2247
2248 ierr = DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2249 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2250
2251 // -X Face
2252 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC && xs == 0) {
2253 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2254 ucat[k][j][xs].x = 2.0 * ubcs[k][j][xs].x - ucat[k][j][xs + 1].x;
2255 ucat[k][j][xs].y = 2.0 * ubcs[k][j][xs].y - ucat[k][j][xs + 1].y;
2256 ucat[k][j][xs].z = 2.0 * ubcs[k][j][xs].z - ucat[k][j][xs + 1].z;
2257 }
2258 }
2259 // +X Face
2260 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC && xe == mx) {
2261 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2262 ucat[k][j][xe-1].x = 2.0 * ubcs[k][j][xe-1].x - ucat[k][j][xe - 2].x;
2263 ucat[k][j][xe-1].y = 2.0 * ubcs[k][j][xe-1].y - ucat[k][j][xe - 2].y;
2264 ucat[k][j][xe-1].z = 2.0 * ubcs[k][j][xe-1].z - ucat[k][j][xe - 2].z;
2265 }
2266 }
2267
2268 // -Y Face
2269 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC && ys == 0) {
2270 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2271 ucat[k][ys][i].x = 2.0 * ubcs[k][ys][i].x - ucat[k][ys + 1][i].x;
2272 ucat[k][ys][i].y = 2.0 * ubcs[k][ys][i].y - ucat[k][ys + 1][i].y;
2273 ucat[k][ys][i].z = 2.0 * ubcs[k][ys][i].z - ucat[k][ys + 1][i].z;
2274 }
2275 }
2276 // +Y Face
2277 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC && ye == my) {
2278 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2279 ucat[k][ye-1][i].x = 2.0 * ubcs[k][ye-1][i].x - ucat[k][ye-2][i].x;
2280 ucat[k][ye-1][i].y = 2.0 * ubcs[k][ye-1][i].y - ucat[k][ye-2][i].y;
2281 ucat[k][ye-1][i].z = 2.0 * ubcs[k][ye-1][i].z - ucat[k][ye-2][i].z;
2282 }
2283 }
2284
2285 // -Z Face
2286 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC && zs == 0) {
2287 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2288 ucat[zs][j][i].x = 2.0 * ubcs[zs][j][i].x - ucat[zs + 1][j][i].x;
2289 ucat[zs][j][i].y = 2.0 * ubcs[zs][j][i].y - ucat[zs + 1][j][i].y;
2290 ucat[zs][j][i].z = 2.0 * ubcs[zs][j][i].z - ucat[zs + 1][j][i].z;
2291 }
2292 }
2293 // +Z Face
2294 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC && ze == mz) {
2295 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2296 ucat[ze-1][j][i].x = 2.0 * ubcs[ze-1][j][i].x - ucat[ze-2][j][i].x;
2297 ucat[ze-1][j][i].y = 2.0 * ubcs[ze-1][j][i].y - ucat[ze-2][j][i].y;
2298 ucat[ze-1][j][i].z = 2.0 * ubcs[ze-1][j][i].z - ucat[ze-2][j][i].z;
2299 }
2300 }
2301
2302 ierr = DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2303 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2304
2305 PetscFunctionReturn(0);
2306}
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:123
UserCtx::Bcs
BCS Bcs
Definition variables.h:934
UserCtx::Ucat
Vec Ucat
Definition variables.h:939
Here is the caller graph for this function:

◆ UpdateCornerNodes()

PetscErrorCode UpdateCornerNodes ( UserCtx *  user)

Internal helper implementation: UpdateCornerNodes().

Updates the corner and edge ghost nodes of the local domain by averaging.

Local to this translation unit.

Definition at line 2314 of file Boundaries.c.

2315{
2316 PetscErrorCode ierr;
2317 DM da = user->da, fda = user->fda;
2318 DMDALocalInfo info = user->info;
2319 PetscInt xs = info.xs, xe = info.xs + info.xm;
2320 PetscInt ys = info.ys, ye = info.ys + info.ym;
2321 PetscInt zs = info.zs, ze = info.zs + info.zm;
2322 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2323
2324 Cmpnts ***ucat;
2325 PetscReal ***p;
2326
2327 PetscFunctionBeginUser;
2328
2329 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2330 ierr = DMDAVecGetArray(da, user->P, &p); CHKERRQ(ierr);
2331
2332 // --- Update Edges and Corners by Averaging ---
2333 // The order of these blocks ensures that corners (where 3 faces meet) are
2334 // computed using data from edges (where 2 faces meet), which are computed first.
2335// Edges connected to the -Z face (k=zs)
2336 if (zs == 0) {
2337 if (xs == 0) {
2338 for (PetscInt j = ys; j < ye; j++) {
2339 p[zs][j][xs] = 0.5 * (p[zs+1][j][xs] + p[zs][j][xs+1]);
2340 ucat[zs][j][xs].x = 0.5 * (ucat[zs+1][j][xs].x + ucat[zs][j][xs+1].x);
2341 ucat[zs][j][xs].y = 0.5 * (ucat[zs+1][j][xs].y + ucat[zs][j][xs+1].y);
2342 ucat[zs][j][xs].z = 0.5 * (ucat[zs+1][j][xs].z + ucat[zs][j][xs+1].z);
2343 }
2344 }
2345 if (xe == mx) {
2346 for (PetscInt j = ys; j < ye; j++) {
2347 p[zs][j][mx-1] = 0.5 * (p[zs+1][j][mx-1] + p[zs][j][mx-2]);
2348 ucat[zs][j][mx-1].x = 0.5 * (ucat[zs+1][j][mx-1].x + ucat[zs][j][mx-2].x);
2349 ucat[zs][j][mx-1].y = 0.5 * (ucat[zs+1][j][mx-1].y + ucat[zs][j][mx-2].y);
2350 ucat[zs][j][mx-1].z = 0.5 * (ucat[zs+1][j][mx-1].z + ucat[zs][j][mx-2].z);
2351 }
2352 }
2353 if (ys == 0) {
2354 for (PetscInt i = xs; i < xe; i++) {
2355 p[zs][ys][i] = 0.5 * (p[zs+1][ys][i] + p[zs][ys+1][i]);
2356 ucat[zs][ys][i].x = 0.5 * (ucat[zs+1][ys][i].x + ucat[zs][ys+1][i].x);
2357 ucat[zs][ys][i].y = 0.5 * (ucat[zs+1][ys][i].y + ucat[zs][ys+1][i].y);
2358 ucat[zs][ys][i].z = 0.5 * (ucat[zs+1][ys][i].z + ucat[zs][ys+1][i].z);
2359 }
2360 }
2361 if (ye == my) {
2362 for (PetscInt i = xs; i < xe; i++) {
2363 p[zs][my-1][i] = 0.5 * (p[zs+1][my-1][i] + p[zs][my-2][i]);
2364 ucat[zs][my-1][i].x = 0.5 * (ucat[zs+1][my-1][i].x + ucat[zs][my-2][i].x);
2365 ucat[zs][my-1][i].y = 0.5 * (ucat[zs+1][my-1][i].y + ucat[zs][my-2][i].y);
2366 ucat[zs][my-1][i].z = 0.5 * (ucat[zs+1][my-1][i].z + ucat[zs][my-2][i].z);
2367 }
2368 }
2369 }
2370
2371 // Edges connected to the +Z face (k=ze-1)
2372 if (ze == mz) {
2373 if (xs == 0) {
2374 for (PetscInt j = ys; j < ye; j++) {
2375 p[mz-1][j][xs] = 0.5 * (p[mz-2][j][xs] + p[mz-1][j][xs+1]);
2376 ucat[mz-1][j][xs].x = 0.5 * (ucat[mz-2][j][xs].x + ucat[mz-1][j][xs+1].x);
2377 ucat[mz-1][j][xs].y = 0.5 * (ucat[mz-2][j][xs].y + ucat[mz-1][j][xs+1].y);
2378 ucat[mz-1][j][xs].z = 0.5 * (ucat[mz-2][j][xs].z + ucat[mz-1][j][xs+1].z);
2379 }
2380 }
2381 if (xe == mx) {
2382 for (PetscInt j = ys; j < ye; j++) {
2383 p[mz-1][j][mx-1] = 0.5 * (p[mz-2][j][mx-1] + p[mz-1][j][mx-2]);
2384 ucat[mz-1][j][mx-1].x = 0.5 * (ucat[mz-2][j][mx-1].x + ucat[mz-1][j][mx-2].x);
2385 ucat[mz-1][j][mx-1].y = 0.5 * (ucat[mz-2][j][mx-1].y + ucat[mz-1][j][mx-2].y);
2386 ucat[mz-1][j][mx-1].z = 0.5 * (ucat[mz-2][j][mx-1].z + ucat[mz-1][j][mx-2].z);
2387 }
2388 }
2389 if (ys == 0) {
2390 for (PetscInt i = xs; i < xe; i++) {
2391 p[mz-1][ys][i] = 0.5 * (p[mz-2][ys][i] + p[mz-1][ys+1][i]);
2392 ucat[mz-1][ys][i].x = 0.5 * (ucat[mz-2][ys][i].x + ucat[mz-1][ys+1][i].x);
2393 ucat[mz-1][ys][i].y = 0.5 * (ucat[mz-2][ys][i].y + ucat[mz-1][ys+1][i].y);
2394 ucat[mz-1][ys][i].z = 0.5 * (ucat[mz-2][ys][i].z + ucat[mz-1][ys+1][i].z);
2395 }
2396 }
2397 if (ye == my) {
2398 for (PetscInt i = xs; i < xe; i++) {
2399 p[mz-1][my-1][i] = 0.5 * (p[mz-2][my-1][i] + p[mz-1][my-2][i]);
2400 ucat[mz-1][my-1][i].x = 0.5 * (ucat[mz-2][my-1][i].x + ucat[mz-1][my-2][i].x);
2401 ucat[mz-1][my-1][i].y = 0.5 * (ucat[mz-2][my-1][i].y + ucat[mz-1][my-2][i].y);
2402 ucat[mz-1][my-1][i].z = 0.5 * (ucat[mz-2][my-1][i].z + ucat[mz-1][my-2][i].z);
2403 }
2404 }
2405 }
2406
2407 // Remaining edges on the XY plane (that are not on Z faces)
2408 if (ys == 0) {
2409 if (xs == 0) {
2410 for (PetscInt k = zs; k < ze; k++) {
2411 p[k][ys][xs] = 0.5 * (p[k][ys+1][xs] + p[k][ys][xs+1]);
2412 ucat[k][ys][xs].x = 0.5 * (ucat[k][ys+1][xs].x + ucat[k][ys][xs+1].x);
2413 ucat[k][ys][xs].y = 0.5 * (ucat[k][ys+1][xs].y + ucat[k][ys][xs+1].y);
2414 ucat[k][ys][xs].z = 0.5 * (ucat[k][ys+1][xs].z + ucat[k][ys][xs+1].z);
2415 }
2416 }
2417 if (xe == mx) {
2418 for (PetscInt k = zs; k < ze; k++) {
2419 p[k][ys][mx-1] = 0.5 * (p[k][ys+1][mx-1] + p[k][ys][mx-2]);
2420 ucat[k][ys][mx-1].x = 0.5 * (ucat[k][ys+1][mx-1].x + ucat[k][ys][mx-2].x);
2421 ucat[k][ys][mx-1].y = 0.5 * (ucat[k][ys+1][mx-1].y + ucat[k][ys][mx-2].y);
2422 ucat[k][ys][mx-1].z = 0.5 * (ucat[k][ys+1][mx-1].z + ucat[k][ys][mx-2].z);
2423 }
2424 }
2425 }
2426
2427 if (ye == my) {
2428 if (xs == 0) {
2429 for (PetscInt k = zs; k < ze; k++) {
2430 p[k][my-1][xs] = 0.5 * (p[k][my-2][xs] + p[k][my-1][xs+1]);
2431 ucat[k][my-1][xs].x = 0.5 * (ucat[k][my-2][xs].x + ucat[k][my-1][xs+1].x);
2432 ucat[k][my-1][xs].y = 0.5 * (ucat[k][my-2][xs].y + ucat[k][my-1][xs+1].y);
2433 ucat[k][my-1][xs].z = 0.5 * (ucat[k][my-2][xs].z + ucat[k][my-1][xs+1].z);
2434 }
2435 }
2436 if (xe == mx) {
2437 for (PetscInt k = zs; k < ze; k++) {
2438 p[k][my-1][mx-1] = 0.5 * (p[k][my-2][mx-1] + p[k][my-1][mx-2]);
2439 ucat[k][my-1][mx-1].x = 0.5 * (ucat[k][my-2][mx-1].x + ucat[k][my-1][mx-2].x);
2440 ucat[k][my-1][mx-1].y = 0.5 * (ucat[k][my-2][mx-1].y + ucat[k][my-1][mx-2].y);
2441 ucat[k][my-1][mx-1].z = 0.5 * (ucat[k][my-2][mx-1].z + ucat[k][my-1][mx-2].z);
2442 }
2443 }
2444 }
2445
2446 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2447 ierr = DMDAVecRestoreArray(da, user->P, &p); CHKERRQ(ierr);
2448
2449 PetscFunctionReturn(0);
2450}
Here is the caller graph for this function:

◆ ApplyWallFunction()

PetscErrorCode ApplyWallFunction ( UserCtx *  user)

Internal helper implementation: ApplyWallFunction().

Applies wall function modeling to near-wall velocities for all wall-type boundaries.

Local to this translation unit.

Definition at line 2458 of file Boundaries.c.

2459{
2460 PetscErrorCode ierr;
2461 SimCtx *simCtx = user->simCtx;
2462 DMDALocalInfo *info = &user->info;
2463
2464 PetscFunctionBeginUser;
2465
2466 // =========================================================================
2467 // STEP 0: Early exit if wall functions are disabled
2468 // =========================================================================
2469 if (!simCtx->wallfunction) {
2470 PetscFunctionReturn(0);
2471 }
2472
2473 LOG_ALLOW(LOCAL, LOG_DEBUG, "Processing wall function boundaries.\n");
2474
2475 // =========================================================================
2476 // STEP 1: Get read/write access to all necessary field arrays
2477 // =========================================================================
2478 Cmpnts ***velocity_cartesian; // Cartesian velocity (modified)
2479 Cmpnts ***velocity_contravariant; // Contravariant velocity (set to zero at walls)
2480 Cmpnts ***velocity_boundary; // Boundary condition velocity (kept at zero)
2481 Cmpnts ***csi, ***eta, ***zet; // Metric tensor components (face normals)
2482 PetscReal ***node_vertex_flag; // Fluid/solid indicator (0=fluid, 1=solid)
2483 PetscReal ***cell_jacobian; // Grid Jacobian (1/volume)
2484 PetscReal ***friction_velocity; // u_tau (friction velocity field)
2485
2486 ierr = DMDAVecGetArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2487 ierr = DMDAVecGetArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2488 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2489 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2490 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2491 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2492 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2493 ierr = DMDAVecGetArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2494 ierr = DMDAVecGetArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2495
2496 // =========================================================================
2497 // STEP 2: Define loop bounds (owned portion of the grid for this MPI rank)
2498 // =========================================================================
2499 PetscInt grid_start_i = info->xs, grid_end_i = info->xs + info->xm;
2500 PetscInt grid_start_j = info->ys, grid_end_j = info->ys + info->ym;
2501 PetscInt grid_start_k = info->zs, grid_end_k = info->zs + info->zm;
2502 PetscInt grid_size_i = info->mx, grid_size_j = info->my, grid_size_k = info->mz;
2503
2504 // Shrunken loop bounds: exclude domain edges and corners to avoid double-counting
2505 PetscInt loop_start_i = grid_start_i, loop_end_i = grid_end_i;
2506 PetscInt loop_start_j = grid_start_j, loop_end_j = grid_end_j;
2507 PetscInt loop_start_k = grid_start_k, loop_end_k = grid_end_k;
2508
2509 if (grid_start_i == 0) loop_start_i = grid_start_i + 1;
2510 if (grid_end_i == grid_size_i) loop_end_i = grid_end_i - 1;
2511 if (grid_start_j == 0) loop_start_j = grid_start_j + 1;
2512 if (grid_end_j == grid_size_j) loop_end_j = grid_end_j - 1;
2513 if (grid_start_k == 0) loop_start_k = grid_start_k + 1;
2514 if (grid_end_k == grid_size_k) loop_end_k = grid_end_k - 1;
2515
2516 // Wall roughness parameter (smooth wall by default, configurable via -wall_roughness).
2517 const PetscReal wall_roughness_height = user->simCtx->wall_roughness_height;
2518
2519 // =========================================================================
2520 // STEP 3: Process each of the 6 domain faces
2521 // =========================================================================
2522 for (int face_index = 0; face_index < 6; face_index++) {
2523 BCFace current_face_id = (BCFace)face_index;
2524 BoundaryFaceConfig *face_config = &user->boundary_faces[current_face_id];
2525
2526 // Only process faces that are mathematical walls (applies to no-slip, moving, slip, etc.)
2527 if (face_config->mathematical_type != WALL) {
2528 continue;
2529 }
2530
2531 // Check if this MPI rank owns part of this face
2532 PetscBool rank_owns_this_face;
2533 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM,
2534 current_face_id, &rank_owns_this_face); CHKERRQ(ierr);
2535
2536 if (!rank_owns_this_face) {
2537 continue;
2538 }
2539
2540 LOG_ALLOW(LOCAL, LOG_TRACE, "Processing Face %d (%s)\n",
2541 current_face_id, BCFaceToString(current_face_id));
2542
2543 // =====================================================================
2544 // Process each face with appropriate indexing
2545 // =====================================================================
2546 switch(current_face_id) {
2547
2548 // =================================================================
2549 // NEGATIVE X FACE (i = 0, normal points in +X direction)
2550 // =================================================================
2551 case BC_FACE_NEG_X: {
2552 if (grid_start_i == 0) {
2553 const PetscInt ghost_cell_index = grid_start_i;
2554 const PetscInt first_interior_cell = grid_start_i + 1;
2555 const PetscInt second_interior_cell = grid_start_i + 2;
2556
2557 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2558 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2559
2560 // Skip if this is a solid cell (embedded boundary)
2561 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2562
2563 // Calculate face area from contravariant metric tensor
2564 PetscReal face_area = sqrt(
2565 csi[k][j][ghost_cell_index].x * csi[k][j][ghost_cell_index].x +
2566 csi[k][j][ghost_cell_index].y * csi[k][j][ghost_cell_index].y +
2567 csi[k][j][ghost_cell_index].z * csi[k][j][ghost_cell_index].z
2568 );
2569
2570 // Compute wall-normal distances using cell Jacobians
2571 // sb = distance from wall to first interior cell center
2572 // sc = distance from wall to second interior cell center
2573 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2574 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2575 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2576
2577 // Compute unit normal vector pointing INTO the domain
2578 PetscReal wall_normal[3];
2579 wall_normal[0] = csi[k][j][ghost_cell_index].x / face_area;
2580 wall_normal[1] = csi[k][j][ghost_cell_index].y / face_area;
2581 wall_normal[2] = csi[k][j][ghost_cell_index].z / face_area;
2582
2583 // Define velocities for wall function calculation
2584 Cmpnts wall_velocity; // Ua = velocity at wall (zero for stationary wall)
2585 Cmpnts reference_velocity; // Uc = velocity at second interior cell
2586
2587 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2588 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2589
2590 // Step 1: Linear interpolation (provides initial guess)
2591 noslip(user, distance_to_second_cell, distance_to_first_cell,
2592 wall_velocity, reference_velocity,
2593 &velocity_cartesian[k][j][first_interior_cell],
2594 wall_normal[0], wall_normal[1], wall_normal[2]);
2595
2596 // Step 2: Apply log-law correction (improves near-wall velocity)
2597 wall_function_loglaw(user, wall_roughness_height,
2598 distance_to_second_cell, distance_to_first_cell,
2599 wall_velocity, reference_velocity,
2600 &velocity_cartesian[k][j][first_interior_cell],
2601 &friction_velocity[k][j][first_interior_cell],
2602 wall_normal[0], wall_normal[1], wall_normal[2]);
2603
2604 // Ensure ghost cell BC remains zero (required for proper extrapolation)
2605 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2606 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2607 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2608 velocity_contravariant[k][j][ghost_cell_index].x = 0.0;
2609 }
2610 }
2611 }
2612 }
2613 } break;
2614
2615 // =================================================================
2616 // POSITIVE X FACE (i = mx-1, normal points in -X direction)
2617 // =================================================================
2618 case BC_FACE_POS_X: {
2619 if (grid_end_i == grid_size_i) {
2620 const PetscInt ghost_cell_index = grid_end_i - 1;
2621 const PetscInt first_interior_cell = grid_end_i - 2;
2622 const PetscInt second_interior_cell = grid_end_i - 3;
2623
2624 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2625 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2626
2627 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2628
2629 PetscReal face_area = sqrt(
2630 csi[k][j][first_interior_cell].x * csi[k][j][first_interior_cell].x +
2631 csi[k][j][first_interior_cell].y * csi[k][j][first_interior_cell].y +
2632 csi[k][j][first_interior_cell].z * csi[k][j][first_interior_cell].z
2633 );
2634
2635 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2636 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2637 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2638
2639 // Note: Normal flipped for +X face to point INTO domain
2640 PetscReal wall_normal[3];
2641 wall_normal[0] = -csi[k][j][first_interior_cell].x / face_area;
2642 wall_normal[1] = -csi[k][j][first_interior_cell].y / face_area;
2643 wall_normal[2] = -csi[k][j][first_interior_cell].z / face_area;
2644
2645 Cmpnts wall_velocity, reference_velocity;
2646 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2647 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2648
2649 noslip(user, distance_to_second_cell, distance_to_first_cell,
2650 wall_velocity, reference_velocity,
2651 &velocity_cartesian[k][j][first_interior_cell],
2652 wall_normal[0], wall_normal[1], wall_normal[2]);
2653
2654 wall_function_loglaw(user, wall_roughness_height,
2655 distance_to_second_cell, distance_to_first_cell,
2656 wall_velocity, reference_velocity,
2657 &velocity_cartesian[k][j][first_interior_cell],
2658 &friction_velocity[k][j][first_interior_cell],
2659 wall_normal[0], wall_normal[1], wall_normal[2]);
2660
2661 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2662 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2663 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2664 velocity_contravariant[k][j][first_interior_cell].x = 0.0;
2665 }
2666 }
2667 }
2668 }
2669 } break;
2670
2671 // =================================================================
2672 // NEGATIVE Y FACE (j = 0, normal points in +Y direction)
2673 // =================================================================
2674 case BC_FACE_NEG_Y: {
2675 if (grid_start_j == 0) {
2676 const PetscInt ghost_cell_index = grid_start_j;
2677 const PetscInt first_interior_cell = grid_start_j + 1;
2678 const PetscInt second_interior_cell = grid_start_j + 2;
2679
2680 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2681 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2682
2683 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2684
2685 PetscReal face_area = sqrt(
2686 eta[k][ghost_cell_index][i].x * eta[k][ghost_cell_index][i].x +
2687 eta[k][ghost_cell_index][i].y * eta[k][ghost_cell_index][i].y +
2688 eta[k][ghost_cell_index][i].z * eta[k][ghost_cell_index][i].z
2689 );
2690
2691 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2692 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2693 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2694
2695 PetscReal wall_normal[3];
2696 wall_normal[0] = eta[k][ghost_cell_index][i].x / face_area;
2697 wall_normal[1] = eta[k][ghost_cell_index][i].y / face_area;
2698 wall_normal[2] = eta[k][ghost_cell_index][i].z / face_area;
2699
2700 Cmpnts wall_velocity, reference_velocity;
2701 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2702 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2703
2704 noslip(user, distance_to_second_cell, distance_to_first_cell,
2705 wall_velocity, reference_velocity,
2706 &velocity_cartesian[k][first_interior_cell][i],
2707 wall_normal[0], wall_normal[1], wall_normal[2]);
2708
2709 wall_function_loglaw(user, wall_roughness_height,
2710 distance_to_second_cell, distance_to_first_cell,
2711 wall_velocity, reference_velocity,
2712 &velocity_cartesian[k][first_interior_cell][i],
2713 &friction_velocity[k][first_interior_cell][i],
2714 wall_normal[0], wall_normal[1], wall_normal[2]);
2715
2716 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2717 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2718 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2719 velocity_contravariant[k][ghost_cell_index][i].y = 0.0;
2720 }
2721 }
2722 }
2723 }
2724 } break;
2725
2726 // =================================================================
2727 // POSITIVE Y FACE (j = my-1, normal points in -Y direction)
2728 // =================================================================
2729 case BC_FACE_POS_Y: {
2730 if (grid_end_j == grid_size_j) {
2731 const PetscInt ghost_cell_index = grid_end_j - 1;
2732 const PetscInt first_interior_cell = grid_end_j - 2;
2733 const PetscInt second_interior_cell = grid_end_j - 3;
2734
2735 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2736 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2737
2738 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2739
2740 PetscReal face_area = sqrt(
2741 eta[k][first_interior_cell][i].x * eta[k][first_interior_cell][i].x +
2742 eta[k][first_interior_cell][i].y * eta[k][first_interior_cell][i].y +
2743 eta[k][first_interior_cell][i].z * eta[k][first_interior_cell][i].z
2744 );
2745
2746 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2747 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2748 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2749
2750 PetscReal wall_normal[3];
2751 wall_normal[0] = -eta[k][first_interior_cell][i].x / face_area;
2752 wall_normal[1] = -eta[k][first_interior_cell][i].y / face_area;
2753 wall_normal[2] = -eta[k][first_interior_cell][i].z / face_area;
2754
2755 Cmpnts wall_velocity, reference_velocity;
2756 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2757 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2758
2759 noslip(user, distance_to_second_cell, distance_to_first_cell,
2760 wall_velocity, reference_velocity,
2761 &velocity_cartesian[k][first_interior_cell][i],
2762 wall_normal[0], wall_normal[1], wall_normal[2]);
2763
2764 wall_function_loglaw(user, wall_roughness_height,
2765 distance_to_second_cell, distance_to_first_cell,
2766 wall_velocity, reference_velocity,
2767 &velocity_cartesian[k][first_interior_cell][i],
2768 &friction_velocity[k][first_interior_cell][i],
2769 wall_normal[0], wall_normal[1], wall_normal[2]);
2770
2771 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2772 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2773 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2774 velocity_contravariant[k][first_interior_cell][i].y = 0.0;
2775 }
2776 }
2777 }
2778 }
2779 } break;
2780
2781 // =================================================================
2782 // NEGATIVE Z FACE (k = 0, normal points in +Z direction)
2783 // =================================================================
2784 case BC_FACE_NEG_Z: {
2785 if (grid_start_k == 0) {
2786 const PetscInt ghost_cell_index = grid_start_k;
2787 const PetscInt first_interior_cell = grid_start_k + 1;
2788 const PetscInt second_interior_cell = grid_start_k + 2;
2789
2790 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2791 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2792
2793 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2794
2795 PetscReal face_area = sqrt(
2796 zet[ghost_cell_index][j][i].x * zet[ghost_cell_index][j][i].x +
2797 zet[ghost_cell_index][j][i].y * zet[ghost_cell_index][j][i].y +
2798 zet[ghost_cell_index][j][i].z * zet[ghost_cell_index][j][i].z
2799 );
2800
2801 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2802 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2803 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2804
2805 PetscReal wall_normal[3];
2806 wall_normal[0] = zet[ghost_cell_index][j][i].x / face_area;
2807 wall_normal[1] = zet[ghost_cell_index][j][i].y / face_area;
2808 wall_normal[2] = zet[ghost_cell_index][j][i].z / face_area;
2809
2810 Cmpnts wall_velocity, reference_velocity;
2811 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2812 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2813
2814 noslip(user, distance_to_second_cell, distance_to_first_cell,
2815 wall_velocity, reference_velocity,
2816 &velocity_cartesian[first_interior_cell][j][i],
2817 wall_normal[0], wall_normal[1], wall_normal[2]);
2818
2819 wall_function_loglaw(user, wall_roughness_height,
2820 distance_to_second_cell, distance_to_first_cell,
2821 wall_velocity, reference_velocity,
2822 &velocity_cartesian[first_interior_cell][j][i],
2823 &friction_velocity[first_interior_cell][j][i],
2824 wall_normal[0], wall_normal[1], wall_normal[2]);
2825
2826 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2827 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2828 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2829 velocity_contravariant[ghost_cell_index][j][i].z = 0.0;
2830 }
2831 }
2832 }
2833 }
2834 } break;
2835
2836 // =================================================================
2837 // POSITIVE Z FACE (k = mz-1, normal points in -Z direction)
2838 // =================================================================
2839 case BC_FACE_POS_Z: {
2840 if (grid_end_k == grid_size_k) {
2841 const PetscInt ghost_cell_index = grid_end_k - 1;
2842 const PetscInt first_interior_cell = grid_end_k - 2;
2843 const PetscInt second_interior_cell = grid_end_k - 3;
2844
2845 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2846 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2847
2848 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2849
2850 PetscReal face_area = sqrt(
2851 zet[first_interior_cell][j][i].x * zet[first_interior_cell][j][i].x +
2852 zet[first_interior_cell][j][i].y * zet[first_interior_cell][j][i].y +
2853 zet[first_interior_cell][j][i].z * zet[first_interior_cell][j][i].z
2854 );
2855
2856 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2857 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2858 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2859
2860 PetscReal wall_normal[3];
2861 wall_normal[0] = -zet[first_interior_cell][j][i].x / face_area;
2862 wall_normal[1] = -zet[first_interior_cell][j][i].y / face_area;
2863 wall_normal[2] = -zet[first_interior_cell][j][i].z / face_area;
2864
2865 Cmpnts wall_velocity, reference_velocity;
2866 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2867 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2868
2869 noslip(user, distance_to_second_cell, distance_to_first_cell,
2870 wall_velocity, reference_velocity,
2871 &velocity_cartesian[first_interior_cell][j][i],
2872 wall_normal[0], wall_normal[1], wall_normal[2]);
2873
2874 wall_function_loglaw(user, wall_roughness_height,
2875 distance_to_second_cell, distance_to_first_cell,
2876 wall_velocity, reference_velocity,
2877 &velocity_cartesian[first_interior_cell][j][i],
2878 &friction_velocity[first_interior_cell][j][i],
2879 wall_normal[0], wall_normal[1], wall_normal[2]);
2880
2881 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2882 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2883 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2884 velocity_contravariant[first_interior_cell][j][i].z = 0.0;
2885 }
2886 }
2887 }
2888 }
2889 } break;
2890 }
2891 }
2892
2893 // =========================================================================
2894 // STEP 4: Restore all arrays and release memory
2895 // =========================================================================
2896 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2897 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2898 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2899 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2900 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2901 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2902 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2903 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2904 ierr = DMDAVecRestoreArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2905
2906 LOG_ALLOW(LOCAL, LOG_DEBUG, "Complete.\n");
2907
2908 PetscFunctionReturn(0);
2909}
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)
Implementation of CanRankServiceFace().
Definition Boundaries.c:127
UserCtx::lFriction_Velocity
Vec lFriction_Velocity
Definition variables.h:935
WALL
@ WALL
Definition variables.h:286
UserCtx::lNvert
Vec lNvert
Definition variables.h:939
UserCtx::KM
PetscInt KM
Definition variables.h:920
UserCtx::lZet
Vec lZet
Definition variables.h:974
UserCtx::Ucont
Vec Ucont
Definition variables.h:939
UserCtx::lCsi
Vec lCsi
Definition variables.h:974
SimCtx::wall_roughness_height
PetscReal wall_roughness_height
Definition variables.h:786
UserCtx::JM
PetscInt JM
Definition variables.h:920
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:822
UserCtx::lAj
Vec lAj
Definition variables.h:974
UserCtx::IM
PetscInt IM
Definition variables.h:920
UserCtx::lEta
Vec lEta
Definition variables.h:974
wall_function_loglaw
void wall_function_loglaw(UserCtx *user, double roughness_height, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z)
Applies log-law wall function with roughness correction.
Definition wallfunction.c:706
noslip
void noslip(UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, double normal_x, double normal_y, double normal_z)
Applies no-slip wall boundary condition with linear interpolation.
Definition wallfunction.c:68
Here is the call graph for this function:
Here is the caller graph for this function:

◆ FinalizePostProjectionCellFields()

PetscErrorCode FinalizePostProjectionCellFields ( UserCtx *  user)

Implementation of FinalizePostProjectionCellFields().

Finalizes cell-centered fields after the projection step.

Full API contract is documented with the header declaration in include/Boundaries.h.

Definition at line 2918 of file Boundaries.c.

2919{
2920 PetscErrorCode ierr;
2921 const FieldId cell_fields[] = {FIELD_ID_UCAT, FIELD_ID_P};
2922
2923 PetscFunctionBeginUser;
2924 PROFILE_FUNCTION_BEGIN;
2925
2926 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing post-projection cell-centered fields.\n");
2927
2928 // Ensure flow-dependent Ubcs handlers see the newly reconstructed Ucat.
2929 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2930 ierr = BoundarySystem_RefreshUbcs(user); CHKERRQ(ierr);
2931
2932 // Establish flat non-periodic faces and periodic endpoints before corners.
2933 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
2934 ierr = SynchronizePeriodicCellFields(user, 2, cell_fields); CHKERRQ(ierr);
2935
2936 // Corner averaging can overwrite periodic endpoints, so restore them after.
2937 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
2938 ierr = SynchronizePeriodicCellFields(user, 2, cell_fields); CHKERRQ(ierr);
2939
2940 // Synchronize explicitly because the periodic helper is a no-op when every
2941 // direction is non-periodic.
2942 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2943 ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
2944
2945 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Post-projection cell-centered fields finalized.\n");
2946 PROFILE_FUNCTION_END;
2947 PetscFunctionReturn(0);
2948}
UpdateDummyCells
PetscErrorCode UpdateDummyCells(UserCtx *user)
Internal helper implementation: UpdateDummyCells().
Definition Boundaries.c:2222
BoundarySystem_RefreshUbcs
PetscErrorCode BoundarySystem_RefreshUbcs(UserCtx *user)
Internal helper implementation: BoundarySystem_RefreshUbcs().
Definition Boundaries.c:1440
UpdateCornerNodes
PetscErrorCode UpdateCornerNodes(UserCtx *user)
Internal helper implementation: UpdateCornerNodes().
Definition Boundaries.c:2314
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ApplyBoundaryConditions()

PetscErrorCode ApplyBoundaryConditions ( UserCtx *  user)

Implementation of ApplyBoundaryConditions().

Main boundary-condition orchestrator executed during solver timestepping.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Boundaries.h.

See also
ApplyBoundaryConditions()

Definition at line 2958 of file Boundaries.c.

2959{
2960 PetscErrorCode ierr;
2961 const FieldId staggered_fields[] = {FIELD_ID_UCONT};
2962 PetscFunctionBeginUser;
2963 PROFILE_FUNCTION_BEGIN;
2964
2965 LOG_ALLOW(GLOBAL,LOG_TRACE,"Boundary Condition Application begins.\n");
2966
2967 // STEP 1: Main iteration loop for applying and converging non-periodic BCs.
2968 // The number of iterations (e.g., 3) allows information to propagate
2969 // between coupled boundaries, like an inlet and a conserving outlet.
2970 for (PetscInt iter = 0; iter < 3; iter++) {
2971 // (a) Execute the boundary system. This phase calculates fluxes across
2972 // the domain and then applies the physical logic for each non-periodic
2973 // handler, setting the `ubcs` (boundary value) array.
2974 ierr = BoundarySystem_ExecuteStep(user); CHKERRQ(ierr);
2975
2976 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Boundary Condition Setup Executed.\n");
2977
2978 // (b) Synchronize the updated ghost cells across all processors to ensure
2979 // all ucont values are current before updating the dummy cells.
2980 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
2981
2982 // (c) Convert updated Contravariant velocities to Cartesian velocities.
2983 ierr = Contra2Cart(user); CHKERRQ(ierr);
2984
2985 // (d) Synchronize the updated Cartesian velocities across all processors
2986 // to ensure all ucat values are current before updating the dummy cells.
2987 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2988
2989 // (e) If Wall functions are enabled, apply them now to adjust near-wall velocities.
2990 if(user->simCtx->wallfunction){
2991 // Apply wall function adjustments to the boundary velocities.
2992 ierr = ApplyWallFunction(user); CHKERRQ(ierr);
2993
2994 // Synchronize the updated Cartesian velocities after wall function adjustments.
2995 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2996
2997 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Wall Function Applied at Walls.\n");
2998 }
2999
3000 // (f) Update the first layer of ghost cells for non-periodic faces using
3001 // the newly computed `ubcs` values.
3002 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
3003
3004 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Dummy Cells/Ghost Cells Updated.\n");
3005
3006 // (g) Handle all periodic boundaries. This is a parallel direct copy
3007 // that sets the absolute constraints for the rest of the solve.
3008 // There is a Ghost update happening inside this function.
3009 ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
3010
3011 // (h) Update the corner and edge ghost nodes. This routine calculates
3012 // values for corners/edges by averaging their neighbors, which have been
3013 // finalized in the steps above (both periodic and non-periodic).
3014 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
3015
3016 // (i) Synchronize the updated edge and corner cells across all processors to ensure
3017 // consistency before the next iteration or finalization.
3018 ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
3019 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
3020 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
3021
3022 // (j) Ensure All the corners are synchronized with a well defined protocol in case of Periodic boundary conditions
3023 // To avoid race conditions.
3024 const FieldId all_fields[] = {FIELD_ID_UCAT, FIELD_ID_P, FIELD_ID_NVERT};
3025 ierr = SynchronizePeriodicCellFields(user, 3, all_fields); CHKERRQ(ierr);
3026
3027 }
3028
3029 // STEP 3: Final ghost node synchronization. This ensures all changes made
3030 // to the global vectors are reflected in the local ghost regions of all
3031 // processors, making the state fully consistent before the next solver stage.
3032 ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
3033 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
3034 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
3035
3036 PROFILE_FUNCTION_END;
3037 PetscFunctionReturn(0);
3038}
ApplyPeriodicBCs
PetscErrorCode ApplyPeriodicBCs(UserCtx *user)
Internal helper implementation: ApplyPeriodicBCs().
Definition Boundaries.c:2169
ApplyWallFunction
PetscErrorCode ApplyWallFunction(UserCtx *user)
Internal helper implementation: ApplyWallFunction().
Definition Boundaries.c:2458
BoundarySystem_ExecuteStep
PetscErrorCode BoundarySystem_ExecuteStep(UserCtx *user)
Implementation of BoundarySystem_ExecuteStep().
Definition Boundaries.c:1018
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2649
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8