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.

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 const PetscReal forceScalingFactor = simCtx->forceScalingFactor;
72 PetscReal drivingForceMagnitude = (bulkVelocityCorrection / simCtx->dt / 1.0 * COEF_TIME_ACCURACY); // replaced simCtx->st with 1.0.
73 drivingForceMagnitude = (simCtx->drivingForceMagnitude * 0.5) + (drivingForceMagnitude * 0.5);
74
75 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Previous driving force: %le\n", simCtx->drivingForceMagnitude);
76 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - New smoothed driving force: %le\n", drivingForceMagnitude);
77 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Force scaling factor: %f\n", simCtx->forceScalingFactor);
78
79 PetscBool hasLoggedApplication = PETSC_FALSE; // Flag to log details only once per rank.
80 // --- Step 5: Apply the momentum source to the correct RHS component ---
81 for (k = lzs; k < lze; k++) {
82 for (j = lys; j < lye; j++) {
83 for (i = lxs; i < lxe; i++) {
84 if (nvert[k][j][i] < 0.1) { // Apply only to fluid cells
85 PetscReal faceArea = 0.0;
86 PetscReal momentumSource = 0.0;
87
88 switch (drivenDirection) {
89 case 'X':
90 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);
91 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
92 rct[k][j][i].x += momentumSource;
93
94 // Log details for the very first point where force is applied on this rank.
95 if (!hasLoggedApplication) {
96 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
97 hasLoggedApplication = PETSC_TRUE;
98 }
99 break;
100 case 'Y':
101 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);
102 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
103 rct[k][j][i].y += momentumSource;
104
105 // Log details for the very first point where force is applied on this rank.
106 if (!hasLoggedApplication) {
107 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
108 hasLoggedApplication = PETSC_TRUE;
109 }
110 break;
111 case 'Z':
112 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);
113 momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
114 rct[k][j][i].z += momentumSource;
115
116 // Log details for the very first point where force is applied on this rank.
117 if (!hasLoggedApplication) {
118 LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
119 hasLoggedApplication = PETSC_TRUE;
120 }
121 break;
122 }
123 }
124 }
125 }
126 }
127
128 // Update simCtx with latest Driving Force
129 simCtx->drivingForceMagnitude = drivingForceMagnitude;
130
131 // --- Step 6: Restore arrays ---
132 ierr = DMDAVecRestoreArray(user->fda, Rct, &rct); CHKERRQ(ierr);
133 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
134 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
135 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
136 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
137
138 PetscFunctionReturn(0);
139}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
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:199
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
SimCtx::rank
PetscMPIInt rank
Definition variables.h:646
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:829
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
SimCtx::forceScalingFactor
PetscReal forceScalingFactor
Definition variables.h:723
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
UserCtx::lZet
Vec lZet
Definition variables.h:858
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:271
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:287
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:286
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:337
UserCtx::_this
PetscInt _this
Definition variables.h:824
SimCtx::dt
PetscReal dt
Definition variables.h:658
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:725
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::info
DMDALocalInfo info
Definition variables.h:818
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:858
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:723
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:247
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:246
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:643
Here is the caller graph for this function:
Generated by doxygen 1.9.8