PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
BodyForces.h File Reference
#include "variables.h"
#include "logging.h"
#include "Metric.h"
Include dependency graph for BodyForces.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

PetscErrorCode ComputeDrivenChannelFlowSource (UserCtx *user, Vec Rct)
 Applies a momentum source term to drive flow in a periodic channel or pipe.
 

Function Documentation

◆ ComputeDrivenChannelFlowSource()

PetscErrorCode ComputeDrivenChannelFlowSource ( UserCtx *  user,
Vec  Rct 
)

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

This function is the "engine" of the driven flow control system. It operates by:

  1. Introspecting the boundary condition handlers to see if a DRIVEN_ flow handler is active on any face. This determines if a driven flow is enabled and in which direction ('X', 'Y', or 'Z').
  2. If a driven flow is active, it reads the bulkVelocityCorrection value that was computed by the handler's PreStep method and stored in the SimCtx.
  3. It translates this velocity correction into a momentum source term.
  4. It adds this source term to the appropriate component of the contravariant RHS vector (Rct) for all fluid cells in the domain.

If no driven flow handler is found, this function does nothing.

Parameters
userThe UserCtx containing the simulation state for a single block.
RctThe PETSc Vec for the contravariant RHS, which will be modified in-place.
Returns
PetscErrorCode 0 on success.

Definition at line 29 of file BodyForces.c.

30{
31 PetscErrorCode ierr;
32 SimCtx *simCtx = user->simCtx;
33 PetscFunctionBeginUser;
34
35 // --- Step 1: Discover if and where a driven flow is active ---
36 char drivenDirection = ' '; // Use space as a null/not-found indicator
37 for (int i = 0; i < 6; i++) {
38 BCHandlerType handler_type = user->boundary_faces[i].handler_type;
39 if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
40 handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
41 {
42 switch (user->boundary_faces[i].face_id) {
43 case BC_FACE_NEG_X: case BC_FACE_POS_X: drivenDirection = 'X'; break;
44 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: drivenDirection = 'Y'; break;
45 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: drivenDirection = 'Z'; break;
46 }
47 break; // Found it, no need to check other faces
48 }
49 }
50
51 // --- Step 2: Early exit if no driven flow is configured ---
52 if (drivenDirection == ' ') {
53 PetscFunctionReturn(0);
54 }
55
56 // --- Step 3: Get the control signal and exit if no correction is needed ---
57 PetscReal bulkVelocityCorrection = simCtx->bulkVelocityCorrection;
58 if (PetscAbsReal(bulkVelocityCorrection) < 1.0e-12) {
59 PetscFunctionReturn(0);
60 }
61
62 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Applying driven flow momentum source in '%c' direction.\n",
63 simCtx->rank, user->_this, drivenDirection);
64 LOG_ALLOW(LOCAL, LOG_DEBUG, " - Received Bulk Velocity Correction: %le\n", bulkVelocityCorrection);
65
66 // --- Step 4: Setup for calculation ---
67 DMDALocalInfo info = user->info;
68 PetscInt i, j, k;
69 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
70 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
71 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
72 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
73 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
74 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
75
76 Cmpnts ***rct, ***csi, ***eta, ***zet;
77 PetscReal ***nvert;
78 ierr = DMDAVecGetArray(user->fda, Rct, &rct); CHKERRQ(ierr);
79 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
80 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
81 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
82 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
83
84 // Calculate the driving force magnitude for the current timestep, smoothed
85 // with the value from the previous step for stability.
86 const PetscReal forceScalingFactor = simCtx->forceScalingFactor;
87 PetscReal drivingForceMagnitude = (bulkVelocityCorrection / simCtx->dt / 1.0 * COEF_TIME_ACCURACY); // replaced simCtx->st with 1.0.
88 drivingForceMagnitude = (simCtx->drivingForceMagnitude * 0.5) + (drivingForceMagnitude * 0.5);
89
90 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Previous driving force: %le\n", simCtx->drivingForceMagnitude);
91 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - New smoothed driving force: %le\n", drivingForceMagnitude);
92 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Force scaling factor: %f\n", simCtx->forceScalingFactor);
93
94 PetscBool hasLoggedApplication = PETSC_FALSE; // Flag to log details only once per rank.
95 // --- Step 5: Apply the momentum source to the correct RHS component ---
96 for (k = lzs; k < lze; k++) {
97 for (j = lys; j < lye; j++) {
98 for (i = lxs; i < lxe; i++) {
99 if (nvert[k][j][i] < 0.1) { // Apply only to fluid cells
100 PetscReal faceArea = 0.0;
101 PetscReal momentumSource = 0.0;
102
103 switch (drivenDirection) {
104 case 'X':
105 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);
106 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
107 rct[k][j][i].x += momentumSource;
108
109 // Log details for the very first point where force is applied on this rank.
110 if (!hasLoggedApplication) {
111 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
112 hasLoggedApplication = PETSC_TRUE;
113 }
114 break;
115 case 'Y':
116 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);
117 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
118 rct[k][j][i].y += momentumSource;
119
120 // Log details for the very first point where force is applied on this rank.
121 if (!hasLoggedApplication) {
122 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
123 hasLoggedApplication = PETSC_TRUE;
124 }
125 break;
126 case 'Z':
127 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);
128 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
129 rct[k][j][i].z += momentumSource;
130
131 // Log details for the very first point where force is applied on this rank.
132 if (!hasLoggedApplication) {
133 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
134 hasLoggedApplication = PETSC_TRUE;
135 }
136 break;
137 }
138 }
139 }
140 }
141 }
142
143 // Update simCtx with latest Driving Force
144 simCtx->drivingForceMagnitude = drivingForceMagnitude;
145
146 // --- Step 6: Restore arrays ---
147 ierr = DMDAVecRestoreArray(user->fda, Rct, &rct); CHKERRQ(ierr);
148 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
149 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
150 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
151 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
152
153 PetscFunctionReturn(0);
154}
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:588
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:746
UserCtx::lNvert
Vec lNvert
Definition variables.h:755
SimCtx::forceScalingFactor
PetscReal forceScalingFactor
Definition variables.h:660
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:731
UserCtx::lZet
Vec lZet
Definition variables.h:776
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:230
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:246
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:245
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:296
UserCtx::_this
PetscInt _this
Definition variables.h:741
SimCtx::dt
PetscReal dt
Definition variables.h:599
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:662
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::lCsi
Vec lCsi
Definition variables.h:776
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::info
DMDALocalInfo info
Definition variables.h:735
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:776
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:660
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:204
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:206
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:205
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:206
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:204
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:205
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:585
Here is the caller graph for this function: