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

Go to the source code of this file.

Macros

#define __FUNCT__   "DimensionalizeField"
 
#define __FUNCT__   "DimensionalizeAllLoadedFields"
 
#define __FUNCT__   "ComputeNodalAverage"
 
#define __FUNCT__   "ComputeQCriterion"
 
#define __FUNCT__   "NormalizeRelativeField"
 
#define __FUNCT__   "ComputeSpecificKE"
 
#define __FUNCT__   "ComputeDisplacement"
 

Functions

PetscErrorCode DimensionalizeField (UserCtx *user, const char *field_name)
 Implementation of DimensionalizeField().
 
PetscErrorCode DimensionalizeAllLoadedFields (UserCtx *user)
 Internal helper implementation: DimensionalizeAllLoadedFields().
 
PetscErrorCode ComputeNodalAverage (UserCtx *user, const char *in_field_name, const char *out_field_name)
 Implementation of ComputeNodalAverage().
 
PetscErrorCode ComputeQCriterion (UserCtx *user)
 Implementation of ComputeQCriterion().
 
PetscErrorCode NormalizeRelativeField (UserCtx *user, const char *relative_field_name)
 Implementation of NormalizeRelativeField().
 
PetscErrorCode ComputeSpecificKE (UserCtx *user, const char *velocity_field, const char *ske_field)
 Internal helper implementation: ComputeSpecificKE().
 
PetscErrorCode ComputeDisplacement (UserCtx *user, const char *disp_field)
 Internal helper implementation: ComputeDisplacement().
 

Macro Definition Documentation

◆ __FUNCT__ [1/7]

#define __FUNCT__   "DimensionalizeField"

Definition at line 5 of file postprocessing_kernels.c.

◆ __FUNCT__ [2/7]

#define __FUNCT__   "DimensionalizeAllLoadedFields"

Definition at line 5 of file postprocessing_kernels.c.

◆ __FUNCT__ [3/7]

#define __FUNCT__   "ComputeNodalAverage"

Definition at line 5 of file postprocessing_kernels.c.

◆ __FUNCT__ [4/7]

#define __FUNCT__   "ComputeQCriterion"

Definition at line 5 of file postprocessing_kernels.c.

◆ __FUNCT__ [5/7]

#define __FUNCT__   "NormalizeRelativeField"

Definition at line 5 of file postprocessing_kernels.c.

◆ __FUNCT__ [6/7]

#define __FUNCT__   "ComputeSpecificKE"

Definition at line 5 of file postprocessing_kernels.c.

◆ __FUNCT__ [7/7]

#define __FUNCT__   "ComputeDisplacement"

Definition at line 5 of file postprocessing_kernels.c.

Function Documentation

◆ DimensionalizeField()

PetscErrorCode DimensionalizeField ( UserCtx *  user,
const char *  field_name 
)

Implementation of DimensionalizeField().

Scales a specified field from non-dimensional to dimensional units in-place.

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

See also
DimensionalizeField()

Definition at line 12 of file postprocessing_kernels.c.

13{
14 PetscErrorCode ierr;
15 SimCtx *simCtx = user->simCtx;
16 Vec target_vec = NULL;
17 PetscReal scale_factor = 1.0;
18 char field_type[64] = "Unknown";
19 PetscBool is_swarm_field = PETSC_FALSE; // Flag for special swarm handling
20 const char *swarm_field_name = NULL; // Name of the field within the swarm
21
22 PetscFunctionBeginUser;
23 PROFILE_FUNCTION_BEGIN;
24 if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx is NULL.");
25 if (!field_name) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "field_name is NULL.");
26
27 // --- 1. Identify the target Vec and the correct scaling factor ---
28 if (strcasecmp(field_name, "Ucat") == 0) {
29 target_vec = user->Ucat;
30 scale_factor = simCtx->scaling.U_ref;
31 strcpy(field_type, "Cartesian Velocity (L/T)");
32 } else if (strcasecmp(field_name, "Ucont") == 0) {
33 target_vec = user->Ucont;
34 scale_factor = simCtx->scaling.U_ref * simCtx->scaling.L_ref * simCtx->scaling.L_ref;
35 strcpy(field_type, "Contravariant Volume Flux (L^3/T)");
36 } else if (strcasecmp(field_name, "P") == 0) {
37 target_vec = user->P;
38 scale_factor = simCtx->scaling.P_ref;
39 strcpy(field_type, "Pressure (M L^-1 T^-2)");
40 } else if (strcasecmp(field_name, "Coordinates") == 0) {
41 ierr = DMGetCoordinates(user->da, &target_vec); CHKERRQ(ierr);
42 scale_factor = simCtx->scaling.L_ref;
43 strcpy(field_type, "Grid Coordinates (L)");
44 } else if (strcasecmp(field_name, "ParticlePosition") == 0) {
45 is_swarm_field = PETSC_TRUE;
46 swarm_field_name = "position";
47 scale_factor = simCtx->scaling.L_ref;
48 strcpy(field_type, "Particle Position (L)");
49 } else if (strcasecmp(field_name, "ParticleVelocity") == 0) {
50 is_swarm_field = PETSC_TRUE;
51 swarm_field_name = "velocity";
52 scale_factor = simCtx->scaling.U_ref;
53 strcpy(field_type, "Particle Velocity (L/T)");
54 } else {
55 LOG(GLOBAL, LOG_WARNING, "DimensionalizeField: Unknown or unhandled field_name '%s'. Field will not be scaled.\n", field_name);
56 PROFILE_FUNCTION_END;
57 PetscFunctionReturn(0);
58 }
59
60 // --- 2. Check for trivial scaling ---
61 if (PetscAbsReal(scale_factor - 1.0) < PETSC_MACHINE_EPSILON) {
62 LOG(GLOBAL, LOG_DEBUG, "DimensionalizeField: Scaling factor for '%s' is 1.0. Skipping operation.\n", field_name);
63 PROFILE_FUNCTION_END;
64 PetscFunctionReturn(0);
65 }
66
67 // --- 3. Perform the in-place scaling operation ---
68 LOG(GLOBAL, LOG_INFO, "Scaling '%s' field (%s) by factor %.4e.\n", field_name, field_type, scale_factor);
69
70 if (is_swarm_field) {
71 // Special handling for DMSwarm fields
72 ierr = DMSwarmCreateGlobalVectorFromField(user->swarm, swarm_field_name, &target_vec); CHKERRQ(ierr);
73 ierr = VecScale(target_vec, scale_factor); CHKERRQ(ierr);
74 ierr = DMSwarmDestroyGlobalVectorFromField(user->swarm, swarm_field_name, &target_vec); CHKERRQ(ierr);
75 } else {
76 // Standard handling for PETSc Vecs
77 if (target_vec) {
78 ierr = VecScale(target_vec, scale_factor); CHKERRQ(ierr);
79 } else {
80 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "Target vector for field '%s' was not found or is NULL.", field_name);
81 }
82 }
83
84 // --- 4. Post-scaling updates for special cases ---
85 if (strcasecmp(field_name, "Coordinates") == 0) {
86 Vec l_coords;
87 ierr = DMGetCoordinates(user->da, &target_vec); CHKERRQ(ierr); // Re-get the global vector
88 ierr = DMGetCoordinatesLocal(user->da, &l_coords); CHKERRQ(ierr);
89 ierr = DMGlobalToLocalBegin(user->fda, target_vec, INSERT_VALUES, l_coords); CHKERRQ(ierr);
90 ierr = DMGlobalToLocalEnd(user->fda, target_vec, INSERT_VALUES, l_coords); CHKERRQ(ierr);
91 }
92
93 PROFILE_FUNCTION_END;
94 PetscFunctionReturn(0);
95}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
LOG
#define LOG(scope, level, fmt,...)
Logging macro for PETSc-based applications with scope control.
Definition logging.h:83
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
ScalingCtx::L_ref
PetscReal L_ref
Definition variables.h:625
UserCtx::Ucont
Vec Ucont
Definition variables.h:837
UserCtx::Ucat
Vec Ucat
Definition variables.h:837
SimCtx::scaling
ScalingCtx scaling
Definition variables.h:707
ScalingCtx::P_ref
PetscReal P_ref
Definition variables.h:628
ScalingCtx::U_ref
PetscReal U_ref
Definition variables.h:626
SimCtx
The master context for the entire simulation.
Definition variables.h:643
Here is the caller graph for this function:

◆ DimensionalizeAllLoadedFields()

PetscErrorCode DimensionalizeAllLoadedFields ( UserCtx *  user)

Internal helper implementation: DimensionalizeAllLoadedFields().

Orchestrates the dimensionalization of all relevant fields loaded from a file.

Local to this translation unit.

Definition at line 103 of file postprocessing_kernels.c.

104{
105 PetscErrorCode ierr;
106 SimCtx *simCtx = user->simCtx;
107
108 PetscFunctionBeginUser;
109 PROFILE_FUNCTION_BEGIN;
110
111 LOG(GLOBAL, LOG_INFO, "--- Converting all loaded fields to dimensional units ---\n");
112
113 // Scale the grid itself first
114 ierr = DimensionalizeField(user, "Coordinates"); CHKERRQ(ierr);
115
116 // Scale primary fluid fields
117 ierr = DimensionalizeField(user, "Ucat"); CHKERRQ(ierr);
118 ierr = DimensionalizeField(user, "Ucont"); CHKERRQ(ierr);
119 ierr = DimensionalizeField(user, "P"); CHKERRQ(ierr);
120
121 // If particles are present, scale their fields
122 if (simCtx->np > 0 && user->swarm) {
123 ierr = DimensionalizeField(user, "ParticlePosition"); CHKERRQ(ierr);
124 ierr = DimensionalizeField(user, "ParticleVelocity"); CHKERRQ(ierr);
125 }
126
127 LOG(GLOBAL, LOG_INFO, "--- Field dimensionalization complete ---\n");
128
129 PROFILE_FUNCTION_END;
130 PetscFunctionReturn(0);
131}
DimensionalizeField
PetscErrorCode DimensionalizeField(UserCtx *user, const char *field_name)
Implementation of DimensionalizeField().
Definition postprocessing_kernels.c:12
SimCtx::np
PetscInt np
Definition variables.h:739
Here is the call graph for this function:

◆ ComputeNodalAverage()

PetscErrorCode ComputeNodalAverage ( UserCtx *  user,
const char *  in_field_name,
const char *  out_field_name 
)

Implementation of ComputeNodalAverage().

Interpolates a cell-centered field to nodal locations using local stencil averaging.

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

See also
ComputeNodalAverage()

Definition at line 143 of file postprocessing_kernels.c.

144{
145 PetscErrorCode ierr;
146 Vec in_vec_local = NULL, out_vec_global = NULL;
147 DM dm_in = NULL, dm_out = NULL;
148 PetscInt dof = 0;
149
150 PetscFunctionBeginUser;
151 PROFILE_FUNCTION_BEGIN;
152 LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeNodalAverage on '%s' -> '%s'.\n", in_field_name, out_field_name);
153
154 // --- 1. Map string names to PETSc objects ---
155 if (strcasecmp(in_field_name, "P") == 0) { in_vec_local = user->lP; dm_in = user->da; dof = 1; }
156 else if (strcasecmp(in_field_name, "Ucat") == 0) { in_vec_local = user->lUcat; dm_in = user->fda; dof = 3; }
157 else if (strcasecmp(in_field_name, "Psi") == 0) { in_vec_local = user->lPsi; dm_in = user->da; dof = 1; }
158 // ... (add other fields as needed) ...
159 else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Unknown input field name for nodal averaging: %s", in_field_name);
160
161 if (strcasecmp(out_field_name, "P_nodal") == 0) { out_vec_global = user->P_nodal; dm_out = user->da; }
162 else if (strcasecmp(out_field_name, "Ucat_nodal") == 0) { out_vec_global = user->Ucat_nodal; dm_out = user->fda; }
163 else if (strcasecmp(out_field_name, "Psi_nodal") == 0) { out_vec_global = user->Psi_nodal; dm_out = user->da; }
164 // ... (add other fields as needed) ...
165 else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Unknown output field name for nodal averaging: %s", out_field_name);
166
167 // --- 2. Ensure Input Data Ghosts are Up-to-Date ---
168 ierr = UpdateLocalGhosts(user, in_field_name); CHKERRQ(ierr);
169
170 // --- 3. Get DMDA info and array pointers ---
171 DMDALocalInfo info;
172 ierr = DMDAGetLocalInfo(dm_out, &info); CHKERRQ(ierr);
173
174 if (dof == 1) { // --- Scalar Field Averaging ---
175 const PetscReal ***l_in_arr;
176 PetscReal ***g_out_arr;
177 ierr = DMDAVecGetArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
178 ierr = DMDAVecGetArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
179
180 // Loop over the output NODE locations. The loop bounds match the required
181 // size of the final subsampled grid.
182 for (PetscInt k = info.zs; k < info.zs + info.zm - 1; k++) {
183 for (PetscInt j = info.ys; j < info.ys + info.ym - 1; j++) {
184 for (PetscInt i = info.xs; i < info.xs + info.xm - 1; i++) {
185 g_out_arr[k][j][i] = 0.125 * (l_in_arr[k][j][i] + l_in_arr[k][j][i+1] +
186 l_in_arr[k][j+1][i] + l_in_arr[k][j+1][i+1] +
187 l_in_arr[k+1][j][i] + l_in_arr[k+1][j][i+1] +
188 l_in_arr[k+1][j+1][i] + l_in_arr[k+1][j+1][i+1]);
189 }
190 }
191 }
192 ierr = DMDAVecRestoreArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
193 ierr = DMDAVecRestoreArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
194
195 } else if (dof == 3) { // --- Vector Field Averaging ---
196 const Cmpnts ***l_in_arr;
197 Cmpnts ***g_out_arr;
198 ierr = DMDAVecGetArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
199 ierr = DMDAVecGetArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
200
201 for (PetscInt k = info.zs; k < info.zs + info.zm - 1; k++) {
202 for (PetscInt j = info.ys; j < info.ys + info.ym - 1; j++) {
203 for (PetscInt i = info.xs; i < info.xs + info.xm - 1; i++) {
204 g_out_arr[k][j][i].x = 0.125 * (l_in_arr[k][j][i].x + l_in_arr[k][j][i+1].x +
205 l_in_arr[k][j+1][i].x + l_in_arr[k][j+1][i+1].x +
206 l_in_arr[k+1][j][i].x + l_in_arr[k+1][j][i+1].x +
207 l_in_arr[k+1][j+1][i].x + l_in_arr[k+1][j+1][i+1].x);
208
209 g_out_arr[k][j][i].y = 0.125 * (l_in_arr[k][j][i].y + l_in_arr[k][j][i+1].y +
210 l_in_arr[k][j+1][i].y + l_in_arr[k][j+1][i+1].y +
211 l_in_arr[k+1][j][i].y + l_in_arr[k+1][j][i+1].y +
212 l_in_arr[k+1][j+1][i].y + l_in_arr[k+1][j+1][i+1].y);
213
214 g_out_arr[k][j][i].z = 0.125 * (l_in_arr[k][j][i].z + l_in_arr[k][j][i+1].z +
215 l_in_arr[k][j+1][i].z + l_in_arr[k][j+1][i+1].z +
216 l_in_arr[k+1][j][i].z + l_in_arr[k+1][j][i+1].z +
217 l_in_arr[k+1][j+1][i].z + l_in_arr[k+1][j+1][i+1].z);
218 }
219 }
220 }
221 ierr = DMDAVecRestoreArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
222 ierr = DMDAVecRestoreArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
223 }
224 PROFILE_FUNCTION_END;
225 PetscFunctionReturn(0);
226}
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
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1361
UserCtx::P_nodal
Vec P_nodal
Definition variables.h:887
UserCtx::Ucat_nodal
Vec Ucat_nodal
Definition variables.h:888
UserCtx::lPsi
Vec lPsi
Definition variables.h:883
Cmpnts::x
PetscScalar x
Definition variables.h:101
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Psi_nodal
Vec Psi_nodal
Definition variables.h:890
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
Cmpnts::y
PetscScalar y
Definition variables.h:101
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeQCriterion()

PetscErrorCode ComputeQCriterion ( UserCtx *  user)

Implementation of ComputeQCriterion().

Computes the Q-criterion diagnostic from the local velocity-gradient tensor.

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

See also
ComputeQCriterion()

Definition at line 237 of file postprocessing_kernels.c.

238{
239 PetscErrorCode ierr;
240 DMDALocalInfo info;
241 const Cmpnts ***lucat, ***lcsi, ***leta, ***lzet;
242 const PetscReal***laj, ***lnvert;
243 PetscReal ***gq;
244
245 PetscFunctionBeginUser;
246 PROFILE_FUNCTION_BEGIN;
247 LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeQCriterion.\n");
248
249 // --- 1. Ensure all required ghost values are up-to-date ---
250 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
251 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
252 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
253 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
254 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
255 ierr = UpdateLocalGhosts(user, "Nvert"); CHKERRQ(ierr);
256
257 // --- 2. Get DMDA info and array pointers ---
258 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
259
260 ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, (void*)&lucat); CHKERRQ(ierr);
261 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (void*)&lcsi); CHKERRQ(ierr);
262 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (void*)&leta); CHKERRQ(ierr);
263 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (void*)&lzet); CHKERRQ(ierr);
264 ierr = DMDAVecGetArrayRead(user->da, user->lAj, (void*)&laj); CHKERRQ(ierr);
265 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (void*)&lnvert); CHKERRQ(ierr);
266 ierr = DMDAVecGetArray(user->da, user->Qcrit, (void*)&gq); CHKERRQ(ierr);
267
268 // --- 3. Define Loop Bounds for INTERIOR Cells ---
269 PetscInt i_start = (info.xs == 0) ? 1 : info.xs;
270 PetscInt i_end = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
271 PetscInt j_start = (info.ys == 0) ? 1 : info.ys;
272 PetscInt j_end = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
273 PetscInt k_start = (info.zs == 0) ? 1 : info.zs;
274 PetscInt k_end = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
275
276 // --- 4. Main Computation Loop ---
277 for (PetscInt k = k_start; k < k_end; k++) {
278 for (PetscInt j = j_start; j < j_end; j++) {
279 for (PetscInt i = i_start; i < i_end; i++) {
280
281 // Calculate velocity derivatives in computational space (central differences)
282 PetscReal uc = 0.5 * (lucat[k][j][i+1].x - lucat[k][j][i-1].x);
283 PetscReal vc = 0.5 * (lucat[k][j][i+1].y - lucat[k][j][i-1].y);
284 PetscReal wc = 0.5 * (lucat[k][j][i+1].z - lucat[k][j][i-1].z);
285
286 PetscReal ue = 0.5 * (lucat[k][j+1][i].x - lucat[k][j-1][i].x);
287 PetscReal ve = 0.5 * (lucat[k][j+1][i].y - lucat[k][j-1][i].y);
288 PetscReal we = 0.5 * (lucat[k][j+1][i].z - lucat[k][j-1][i].z);
289
290 PetscReal uz = 0.5 * (lucat[k+1][j][i].x - lucat[k-1][j][i].x);
291 PetscReal vz = 0.5 * (lucat[k+1][j][i].y - lucat[k-1][j][i].y);
292 PetscReal wz = 0.5 * (lucat[k+1][j][i].z - lucat[k-1][j][i].z);
293
294 // Average metrics to the cell center
295 PetscReal csi1 = 0.5 * (lcsi[k][j][i].x + lcsi[k][j][i-1].x) * laj[k][j][i];
296 PetscReal csi2 = 0.5 * (lcsi[k][j][i].y + lcsi[k][j][i-1].y) * laj[k][j][i];
297 PetscReal csi3 = 0.5 * (lcsi[k][j][i].z + lcsi[k][j][i-1].z) * laj[k][j][i];
298
299 PetscReal eta1 = 0.5 * (leta[k][j][i].x + leta[k][j-1][i].x) * laj[k][j][i];
300 PetscReal eta2 = 0.5 * (leta[k][j][i].y + leta[k][j-1][i].y) * laj[k][j][i];
301 PetscReal eta3 = 0.5 * (leta[k][j][i].z + leta[k][j-1][i].z) * laj[k][j][i];
302
303 PetscReal zet1 = 0.5 * (lzet[k][j][i].x + lzet[k-1][j][i].x) * laj[k][j][i];
304 PetscReal zet2 = 0.5 * (lzet[k][j][i].y + lzet[k-1][j][i].y) * laj[k][j][i];
305 PetscReal zet3 = 0.5 * (lzet[k][j][i].z + lzet[k-1][j][i].z) * laj[k][j][i];
306
307 // Calculate velocity gradient tensor components d_ij = du_i/dx_j
308 PetscReal d11 = uc * csi1 + ue * eta1 + uz * zet1;
309 PetscReal d12 = uc * csi2 + ue * eta2 + uz * zet2;
310 PetscReal d13 = uc * csi3 + ue * eta3 + uz * zet3;
311
312 PetscReal d21 = vc * csi1 + ve * eta1 + vz * zet1;
313 PetscReal d22 = vc * csi2 + ve * eta2 + vz * zet2;
314 PetscReal d23 = vc * csi3 + ve * eta3 + vz * zet3;
315
316 PetscReal d31 = wc * csi1 + we * eta1 + wz * zet1;
317 PetscReal d32 = wc * csi2 + we * eta2 + wz * zet2;
318 PetscReal d33 = wc * csi3 + we * eta3 + wz * zet3;
319
320 // Strain-Rate Tensor S_ij = 0.5 * (d_ij + d_ji)
321 PetscReal s11 = d11;
322 PetscReal s12 = 0.5 * (d12 + d21);
323 PetscReal s13 = 0.5 * (d13 + d31);
324 PetscReal s22 = d22;
325 PetscReal s23 = 0.5 * (d23 + d32);
326 PetscReal s33 = d33;
327
328 // Vorticity Tensor Omega_ij = 0.5 * (d_ij - d_ji)
329 PetscReal w12 = 0.5 * (d12 - d21);
330 PetscReal w13 = 0.5 * (d13 - d31);
331 PetscReal w23 = 0.5 * (d23 - d32);
332
333 // Squared norms of the tensors
334 PetscReal s_norm_sq = s11*s11 + s22*s22 + s33*s33 + 2.0*(s12*s12 + s13*s13 + s23*s23);
335 PetscReal w_norm_sq = 2.0 * (w12*w12 + w13*w13 + w23*w23);
336
337 gq[k][j][i] = 0.5 * (w_norm_sq - s_norm_sq);
338
339 if (lnvert[k][j][i] > 0.1) {
340 gq[k][j][i] = 0.0;
341 }
342 }
343 }
344 }
345
346 // --- 5. Restore arrays ---
347 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, (void*)&lucat); CHKERRQ(ierr);
348 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (void*)&lcsi); CHKERRQ(ierr);
349 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (void*)&leta); CHKERRQ(ierr);
350 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (void*)&lzet); CHKERRQ(ierr);
351 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, (void*)&laj); CHKERRQ(ierr);
352 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (void*)&lnvert); CHKERRQ(ierr);
353 ierr = DMDAVecRestoreArray(user->da, user->Qcrit, (void*)&gq); CHKERRQ(ierr);
354
355 PROFILE_FUNCTION_END;
356 PetscFunctionReturn(0);
357}
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
UserCtx::lZet
Vec lZet
Definition variables.h:858
UserCtx::Qcrit
Vec Qcrit
Definition variables.h:889
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
UserCtx::lAj
Vec lAj
Definition variables.h:858
UserCtx::lEta
Vec lEta
Definition variables.h:858
Here is the call graph for this function:
Here is the caller graph for this function:

◆ NormalizeRelativeField()

PetscErrorCode NormalizeRelativeField ( UserCtx *  user,
const char *  relative_field_name 
)

Implementation of NormalizeRelativeField().

Normalizes a relative scalar field using the configured reference pressure scale.

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

See also
NormalizeRelativeField()

Definition at line 367 of file postprocessing_kernels.c.

368{
369 PetscErrorCode ierr;
370 Vec P_vec = NULL;
371 PetscMPIInt rank;
372 PetscInt ip=1, jp=1, kp=1; // Default reference point
373 PetscReal p_ref = 0.0;
374 PetscInt ref_point_global_idx[1];
375 PetscScalar ref_value_local[1];
376 IS is_from, is_to;
377 VecScatter scatter_ctx;
378 Vec ref_vec_seq;
379 PostProcessParams *pps = user->simCtx->pps;
380
381 // Fetch referenc point from pps.
382 ip = pps->reference[0];
383 jp = pps->reference[1];
384 kp = pps->reference[2];
385
386 PetscFunctionBeginUser;
387 PROFILE_FUNCTION_BEGIN;
388 LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running NormalizeRelativeField on '%s'.\n", relative_field_name);
389
390 // --- 1. Map string argument to the PETSc Vec ---
391 if (strcasecmp(relative_field_name, "P") == 0) {
392 P_vec = user->P;
393 } else {
394 SETERRQ(PETSC_COMM_SELF, 1, "NormalizeRelativeField only supports the primary 'P' field , not '%s' currently.", relative_field_name);
395 }
396
397 // --- 2. Get reference point from options and calculate its global DA index ---
398 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
399
400 // Convert the (i,j,k) logical grid coordinates to the global 1D index used by the DMDA vector
401 ref_point_global_idx[0] = kp * (user->IM * user->JM) + jp * user->IM + ip;
402
403 // --- 3. Robustly Scatter the single reference value to rank 0 ---
404 // Create an Index Set (IS) for the source point (from the global vector)
405 ierr = ISCreateGeneral(PETSC_COMM_WORLD, 1, ref_point_global_idx, PETSC_COPY_VALUES, &is_from); CHKERRQ(ierr);
406
407 // Create a sequential vector on rank 0 to hold the result
408 ierr = VecCreateSeq(PETSC_COMM_SELF, 1, &ref_vec_seq); CHKERRQ(ierr);
409
410 // Create an Index Set for the destination point (index 0 of the new sequential vector)
411 PetscInt dest_idx[1] = {0};
412 ierr = ISCreateGeneral(PETSC_COMM_SELF, 1, dest_idx, PETSC_COPY_VALUES, &is_to); CHKERRQ(ierr);
413
414 // Create the scatter context and perform the scatter
415 ierr = VecScatterCreate(P_vec, is_from, ref_vec_seq, is_to, &scatter_ctx); CHKERRQ(ierr);
416 ierr = VecScatterBegin(scatter_ctx, P_vec, ref_vec_seq, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
417 ierr = VecScatterEnd(scatter_ctx, P_vec, ref_vec_seq, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
418
419 // On rank 0, get the value. On other ranks, this will do nothing.
420 if (rank == 0) {
421 ierr = VecGetValues(ref_vec_seq, 1, dest_idx, ref_value_local); CHKERRQ(ierr);
422 p_ref = ref_value_local[0];
423 LOG_ALLOW(LOCAL, LOG_DEBUG, "%s reference point (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") has value %g.\n", relative_field_name, jp, kp, ip, p_ref);
424 }
425
426 // --- 4. Broadcast the reference value from rank 0 to all other processes ---
427 ierr = MPI_Bcast(&p_ref, 1, MPIU_REAL, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
428
429 // --- 5. Perform the normalization (in-place shift) on the full distributed vector ---
430 ierr = VecShift(P_vec, -p_ref); CHKERRQ(ierr);
431 LOG_ALLOW(GLOBAL, LOG_DEBUG, "%s field normalized by subtracting %g.\n", relative_field_name, p_ref);
432
433 // --- 6. Cleanup ---
434 ierr = ISDestroy(&is_from); CHKERRQ(ierr);
435 ierr = ISDestroy(&is_to); CHKERRQ(ierr);
436 ierr = VecScatterDestroy(&scatter_ctx); CHKERRQ(ierr);
437 ierr = VecDestroy(&ref_vec_seq); CHKERRQ(ierr);
438
439 PROFILE_FUNCTION_END;
440 PetscFunctionReturn(0);
441}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
PostProcessParams::reference
PetscInt reference[3]
Definition variables.h:582
UserCtx::JM
PetscInt JM
Definition variables.h:820
SimCtx::pps
PostProcessParams * pps
Definition variables.h:798
UserCtx::IM
PetscInt IM
Definition variables.h:820
PostProcessParams
Holds all configuration parameters for a post-processing run.
Definition variables.h:553
Here is the caller graph for this function:

◆ ComputeSpecificKE()

PetscErrorCode ComputeSpecificKE ( UserCtx *  user,
const char *  velocity_field,
const char *  ske_field 
)

Internal helper implementation: ComputeSpecificKE().

Computes the specific kinetic energy (KE per unit mass) for each particle.

Local to this translation unit.

Definition at line 452 of file postprocessing_kernels.c.

453{
454 PetscErrorCode ierr;
455 PetscInt n_local;
456 const PetscScalar (*vel_arr)[3]; // Access velocity as array of 3-component vectors
457 PetscScalar *ske_arr;
458
459 PetscFunctionBeginUser;
460 PROFILE_FUNCTION_BEGIN;
461 LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeSpecificKE ('%s' -> '%s').\n", velocity_field, ske_field);
462
463 // Get local data arrays from the DMSwarm
464 ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
465 if (n_local == 0) PetscFunctionReturn(0);
466
467 // Get read-only access to velocity and write access to the output field
468 ierr = DMSwarmGetField(user->swarm, velocity_field, NULL, NULL, (void**)&vel_arr); CHKERRQ(ierr);
469 ierr = DMSwarmGetField(user->post_swarm, ske_field, NULL, NULL, (void**)&ske_arr); CHKERRQ(ierr);
470
471 // Main computation loop
472 for (PetscInt p = 0; p < n_local; p++) {
473 const PetscScalar u = vel_arr[p][0];
474 const PetscScalar v = vel_arr[p][1];
475 const PetscScalar w = vel_arr[p][2];
476 const PetscScalar vel_sq = u*u + v*v + w*w;
477 ske_arr[p] = 0.5 * vel_sq;
478 }
479
480 // Restore arrays
481 ierr = DMSwarmRestoreField(user->swarm, velocity_field, NULL, NULL, (void**)&vel_arr); CHKERRQ(ierr);
482 ierr = DMSwarmRestoreField(user->post_swarm, ske_field, NULL, NULL, (void**)&ske_arr); CHKERRQ(ierr);
483
484 PROFILE_FUNCTION_END;
485 PetscFunctionReturn(0);
486}
UserCtx::post_swarm
DM post_swarm
Definition variables.h:886
Here is the caller graph for this function:

◆ ComputeDisplacement()

PetscErrorCode ComputeDisplacement ( UserCtx *  user,
const char *  disp_field 
)

Internal helper implementation: ComputeDisplacement().

Computes the displacement magnitude |r_i - r_0| for each particle (per-particle VTK kernel).

Local to this translation unit.

Definition at line 494 of file postprocessing_kernels.c.

495{
496 PetscErrorCode ierr;
497 PetscInt n_local;
498 const PetscReal (*pos_arr)[3];
499 PetscScalar *disp_out;
500 SimCtx *simCtx = user->simCtx;
501
502 PetscFunctionBeginUser;
503 PROFILE_FUNCTION_BEGIN;
504 LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeDisplacement (-> '%s').\n", disp_field);
505
506 ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
507 if (n_local == 0) PetscFunctionReturn(0);
508
509 const PetscReal x0 = simCtx->psrc_x;
510 const PetscReal y0 = simCtx->psrc_y;
511 const PetscReal z0 = simCtx->psrc_z;
512
513 ierr = DMSwarmGetField(user->swarm, "position", NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
514 ierr = DMSwarmGetField(user->post_swarm, disp_field, NULL, NULL, (void**)&disp_out); CHKERRQ(ierr);
515
516 for (PetscInt p = 0; p < n_local; p++) {
517 const PetscReal dx = pos_arr[p][0] - x0;
518 const PetscReal dy = pos_arr[p][1] - y0;
519 const PetscReal dz = pos_arr[p][2] - z0;
520 disp_out[p] = PetscSqrtReal(dx*dx + dy*dy + dz*dz);
521 }
522
523 ierr = DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
524 ierr = DMSwarmRestoreField(user->post_swarm, disp_field, NULL, NULL, (void**)&disp_out); CHKERRQ(ierr);
525
526 PROFILE_FUNCTION_END;
527 PetscFunctionReturn(0);
528}
SimCtx::psrc_x
PetscReal psrc_x
Definition variables.h:706
SimCtx::psrc_z
PetscReal psrc_z
Point source location for PARTICLE_INIT_POINT_SOURCE.
Definition variables.h:706
SimCtx::psrc_y
PetscReal psrc_y
Definition variables.h:706
Here is the caller graph for this function:
Generated by doxygen 1.9.8