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

Go to the source code of this file.

Macros

#define __FUNCT__   "ComputeDrivenChannelFlowSource"
 

Functions

PetscErrorCode ComputeDrivenChannelFlowSource (UserCtx *user, Vec Rct)
 Internal helper implementation: ComputeDrivenChannelFlowSource().
 

Macro Definition Documentation

◆ __FUNCT__

#define __FUNCT__   "ComputeDrivenChannelFlowSource"

Definition at line 9 of file BodyForces.c.

Function Documentation

◆ ComputeDrivenChannelFlowSource()

PetscErrorCode ComputeDrivenChannelFlowSource ( UserCtx *  user,
Vec  Rct 
)

Internal helper implementation: ComputeDrivenChannelFlowSource().

Applies a momentum source term to drive flow in a periodic channel or pipe.

Local to this translation unit.

Definition at line 14 of file BodyForces.c.

15{
16 PetscErrorCode ierr;
17 SimCtx *simCtx = user->simCtx;
18 PetscFunctionBeginUser;
19
20 // --- Step 1: Discover if and where a driven flow is active ---
21 char drivenDirection = ' '; // Use space as a null/not-found indicator
22 for (int i = 0; i < 6; i++) {
23 BCHandlerType handler_type = user->boundary_faces[i].handler_type;
24 if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
25 handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
26 {
27 switch (user->boundary_faces[i].face_id) {
28 case BC_FACE_NEG_X: case BC_FACE_POS_X: drivenDirection = 'X'; break;
29 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: drivenDirection = 'Y'; break;
30 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: drivenDirection = 'Z'; break;
31 }
32 break; // Found it, no need to check other faces
33 }
34 }
35
36 // --- Step 2: Early exit if no driven flow is configured ---
37 if (drivenDirection == ' ') {
38 PetscFunctionReturn(0);
39 }
40
41 // --- Step 3: Get the control signal and exit if no correction is needed ---
42 PetscReal bulkVelocityCorrection = simCtx->bulkVelocityCorrection;
43 if (PetscAbsReal(bulkVelocityCorrection) < 1.0e-12) {
44 PetscFunctionReturn(0);
45 }
46
47 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Applying driven flow momentum source in '%c' direction.\n",
48 simCtx->rank, user->_this, drivenDirection);
49 LOG_ALLOW(LOCAL, LOG_DEBUG, " - Received Bulk Velocity Correction: %le\n", bulkVelocityCorrection);
50
51 // --- Step 4: Setup for calculation ---
52 DMDALocalInfo info = user->info;
53 PetscInt i, j, k;
54 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
55 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
56 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
57 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
58 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
59 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
60
61 Cmpnts ***rct, ***csi, ***eta, ***zet;
62 PetscReal ***nvert;
63 ierr = DMDAVecGetArray(user->fda, Rct, &rct); CHKERRQ(ierr);
64 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
65 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
66 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
67 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
68
69 // Calculate the driving force magnitude for the current timestep, smoothed
70 // with the value from the previous step for stability.
71 //
72 // ONCE PER PHYSICAL STEP, NOT ONCE PER CALL. This function runs from
73 // ComputeRHS, which executes once per Jameson RK stage under the Picard
74 // solver and once per residual evaluation under Newton-Krylov. The smoothing
75 // below carries state in simCtx across calls, so advancing it every call
76 // would walk the applied force toward its target within a single timestep
77 // (0.5, then 0.75, then 0.875 ... of the way there). The force would then
78 // depend on how many residual evaluations preceded it - history dependence
79 // that MomentumNewtonKrylov_FormResidual() explicitly forbids, and that
80 // breaks the constant-forcing assumption behind the Picard shadow-Jacobian
81 // estimate. Resolve it once for the step and reuse it thereafter.
82 const PetscReal forceScalingFactor = simCtx->forceScalingFactor;
83 PetscReal drivingForceMagnitude;
84
85 if (simCtx->drivingForceStep != simCtx->step) {
86 const PetscReal targetForce = (bulkVelocityCorrection / simCtx->dt / 1.0 * COEF_TIME_ACCURACY); // replaced simCtx->st with 1.0.
87 drivingForceMagnitude = (simCtx->drivingForceMagnitude * 0.5) + (targetForce * 0.5);
88 simCtx->drivingForceMagnitude = drivingForceMagnitude;
89 simCtx->drivingForceStep = simCtx->step;
90 } else {
91 drivingForceMagnitude = simCtx->drivingForceMagnitude;
92 }
93
94 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Previous driving force: %le\n", simCtx->drivingForceMagnitude);
95 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - New smoothed driving force: %le\n", drivingForceMagnitude);
96 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Force scaling factor: %f\n", simCtx->forceScalingFactor);
97
98 PetscBool hasLoggedApplication = PETSC_FALSE; // Flag to log details only once per rank.
99 // --- Step 5: Apply the momentum source to the correct RHS component ---
100 for (k = lzs; k < lze; k++) {
101 for (j = lys; j < lye; j++) {
102 for (i = lxs; i < lxe; i++) {
103 if (nvert[k][j][i] < 0.1) { // Apply only to fluid cells
104 PetscReal faceArea = 0.0;
105 PetscReal momentumSource = 0.0;
106
107 switch (drivenDirection) {
108 case 'X':
109 faceArea = sqrt(csi[k][j][i].x * csi[k][j][i].x + csi[k][j][i].y * csi[k][j][i].y + csi[k][j][i].z * csi[k][j][i].z);
110 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
111 rct[k][j][i].x += momentumSource;
112
113 // Log details for the very first point where force is applied on this rank.
114 if (!hasLoggedApplication) {
115 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
116 hasLoggedApplication = PETSC_TRUE;
117 }
118 break;
119 case 'Y':
120 faceArea = sqrt(eta[k][j][i].x * eta[k][j][i].x + eta[k][j][i].y * eta[k][j][i].y + eta[k][j][i].z * eta[k][j][i].z);
121 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
122 rct[k][j][i].y += momentumSource;
123
124 // Log details for the very first point where force is applied on this rank.
125 if (!hasLoggedApplication) {
126 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
127 hasLoggedApplication = PETSC_TRUE;
128 }
129 break;
130 case 'Z':
131 faceArea = sqrt(zet[k][j][i].x * zet[k][j][i].x + zet[k][j][i].y * zet[k][j][i].y + zet[k][j][i].z * zet[k][j][i].z);
132 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
133 rct[k][j][i].z += momentumSource;
134
135 // Log details for the very first point where force is applied on this rank.
136 if (!hasLoggedApplication) {
137 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
138 hasLoggedApplication = PETSC_TRUE;
139 }
140 break;
141 }
142 }
143 }
144 }
145 }
146
147 // --- Step 6: Restore arrays ---
148 ierr = DMDAVecRestoreArray(user->fda, Rct, &rct); CHKERRQ(ierr);
149 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
150 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
151 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
152 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
153
154 PetscFunctionReturn(0);
155}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
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
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
SimCtx::rank
PetscMPIInt rank
Definition variables.h:698
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:931
UserCtx::lNvert
Vec lNvert
Definition variables.h:939
SimCtx::forceScalingFactor
PetscReal forceScalingFactor
Definition variables.h:801
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
UserCtx::lZet
Vec lZet
Definition variables.h:974
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:303
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
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:369
UserCtx::_this
PetscInt _this
Definition variables.h:924
SimCtx::dt
PetscReal dt
Definition variables.h:710
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:813
Cmpnts::x
PetscScalar x
Definition variables.h:103
UserCtx::lCsi
Vec lCsi
Definition variables.h:974
Cmpnts::z
PetscScalar z
Definition variables.h:103
SimCtx::drivingForceStep
PetscInt drivingForceStep
Definition variables.h:806
SimCtx::step
PetscInt step
Definition variables.h:703
UserCtx::info
DMDALocalInfo info
Definition variables.h:918
Cmpnts::y
PetscScalar y
Definition variables.h:103
UserCtx::lEta
Vec lEta
Definition variables.h:974
COEF_TIME_ACCURACY
#define COEF_TIME_ACCURACY
Coefficient controlling the temporal accuracy scheme (e.g., 1.5 for 2nd Order Backward Difference).
Definition variables.h:57
SimCtx::drivingForceMagnitude
PetscReal drivingForceMagnitude
Definition variables.h:801
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
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102
SimCtx
The master context for the entire simulation.
Definition variables.h:695
Here is the caller graph for this function:
Generated by doxygen 1.9.8