PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
interpolation.c File Reference

Main program for DMSwarm interpolation using the fdf-curvIB method. More...

#include "interpolation.h"
Include dependency graph for interpolation.c:

Go to the source code of this file.

Macros

#define NUM_WEIGHTS   8
 
#define ERROR_MSG_BUFFER_SIZE   256
 
#define __FUNCT__   "InterpolateFieldFromCornerToCenter_Vector"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "InterpolateFieldFromCornerToCenter_Scalar"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "TestCornerToCenterInterpolation"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "InterpolateFieldFromCenterToCorner_Vector"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "InterpolateFieldFromCenterToCorner_Scalar"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT   "PiecWiseLinearInterpolation_Scalar"
 
#define __FUNCT   "PiecWiseLinearInterpolation_Vector"
 
#define __FUNCT   "ComputeTrilinearWeights"
 
#define __FUNCT   "TrilinearInterpolation_Scalar"
 
#define __FUNCT   "TrilinearInterpolation_Vector"
 
#define __FUNCT   "InterpolateEulerFieldToSwarmForParticle"
 
#define __FUNCT__   "InterpolateEulerFieldToSwarm"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "InterpolateAllFieldsToSwarm"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "GetScatterTargetInfo"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "GetPersistentLocalVector"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "AccumulateParticleField"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "NormalizeGridVectorByCount"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "ScatterParticleFieldToEulerField_Internal"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "ScatterParticleFieldToEulerField"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "ScatterAllParticleFieldsToEulerFields"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "InterpolateCornerToFaceCenter_Scalar"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
#define __FUNCT__   "InterpolateCornerToFaceCenter_Vector"
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 

Functions

PetscErrorCode InterpolateFieldFromCornerToCenter_Vector (Cmpnts ***field_arr, Cmpnts ***centfield_arr, UserCtx *user)
 Interpolates a vector field from corner nodes to cell centers, respecting application conventions.
 
PetscErrorCode InterpolateFieldFromCornerToCenter_Scalar (PetscReal ***field_arr, PetscReal ***centfield_arr, UserCtx *user)
 Interpolates a scalar field from corner nodes to cell centers, respecting application conventions.
 
PetscErrorCode TestCornerToCenterInterpolation (UserCtx *user)
 Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.
 
PetscErrorCode InterpolateFieldFromCenterToCorner_Vector (Cmpnts ***centfield_arr, Cmpnts ***corner_arr, UserCtx *user)
 Interpolates a vector field from cell centers to corner nodes.
 
PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar (PetscReal ***centfield_arr, PetscReal ***corner_arr, UserCtx *user)
 Interpolates a scalar field from cell centers to corner nodes.
 
PetscErrorCode PieceWiseLinearInterpolation_Scalar (const char *fieldName, PetscReal ***fieldScal, PetscInt iCell, PetscInt jCell, PetscInt kCell, PetscReal *val)
 Retrieves the scalar value at the cell (iCell, jCell, kCell).
 
PetscErrorCode PieceWiseLinearInterpolation_Vector (const char *fieldName, Cmpnts ***fieldVec, PetscInt iCell, PetscInt jCell, PetscInt kCell, Cmpnts *vec)
 Retrieves the vector (Cmpnts) at the cell (iCell, jCell, kCell).
 
static void ComputeTrilinearWeights (PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
 Computes the trilinear interpolation weights from the interpolation coefficients.
 
PetscErrorCode TrilinearInterpolation_Scalar (const char *fieldName, PetscReal ***fieldScal, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *val)
 Computes the trilinear Interpolateted scalar at a given point.
 
PetscErrorCode TrilinearInterpolation_Vector (const char *fieldName, Cmpnts ***fieldVec, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, Cmpnts *vec)
 Computes the trilinear Interpolateted vector (e.g., velocity) at a given point.
 
static PetscErrorCode InterpolateEulerFieldToSwarmForParticle (const char *fieldName, void *fieldPtr, Particle *particle, void *swarmOut, PetscInt p, PetscInt blockSize)
 Interpolates a single field (scalar or vector) for one particle, storing the result in a swarm array.
 
PetscErrorCode InterpolateEulerFieldToSwarm (UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
 Interpolates a cell-centered field (scalar or vector) onto DMSwarm particles, using a robust, PETSc-idiomatic two-stage process.
 
PetscErrorCode InterpolateAllFieldsToSwarm (UserCtx *user)
 Interpolates all relevant fields from the DMDA to the DMSwarm.
 
PetscErrorCode GetScatterTargetInfo (UserCtx *user, const char *particleFieldName, DM *targetDM, PetscInt *expected_dof)
 Determines the target Eulerian DM and expected DOF for scattering a given particle field.
 
PetscErrorCode GetPersistentLocalVector (UserCtx *user, const char *fieldName, Vec *localVec)
 Retrieves the persistent local vector (e.g., lPsi, lUcat) for a given field name.
 
PetscErrorCode AccumulateParticleField (DM swarm, const char *particleFieldName, DM gridSumDM, Vec gridSumVec)
 Accumulates a particle field (scalar or vector) into a target grid sum vector.
 
PetscErrorCode NormalizeGridVectorByCount (DM countDM, Vec countVec, DM dataDM, Vec sumVec, Vec avgVec)
 Normalizes a grid vector of sums by a grid vector of counts to produce an average.
 
static PetscErrorCode ScatterParticleFieldToEulerField_Internal (UserCtx *user, const char *particleFieldName, DM targetDM, PetscInt expected_dof, Vec eulerFieldAverageVec)
 Internal helper function to orchestrate the scatter operation (accumulate + normalize).
 
PetscErrorCode ScatterParticleFieldToEulerField (UserCtx *user, const char *particleFieldName, Vec eulerFieldAverageVec)
 Scatters a particle field (scalar or vector) to the corresponding Eulerian field average.
 
PetscErrorCode ScatterAllParticleFieldsToEulerFields (UserCtx *user)
 Scatters a predefined set of particle fields to their corresponding Eulerian fields.
 
PetscErrorCode InterpolateCornerToFaceCenter_Scalar (PetscReal ***corner_arr, PetscReal ***faceX_arr, PetscReal ***faceY_arr, PetscReal ***faceZ_arr, UserCtx *user)
 Interpolates a scalar field from corner nodes to all face centers.
 
PetscErrorCode InterpolateCornerToFaceCenter_Vector (Cmpnts ***corner_arr, Cmpnts ***faceX_arr, Cmpnts ***faceY_arr, Cmpnts ***faceZ_arr, UserCtx *user)
 Interpolates a vector field from corner nodes to all face centers.
 

Detailed Description

Main program for DMSwarm interpolation using the fdf-curvIB method.

Provides routines for interpolation between corner-based and center-based fields in the cell-centered DM (fda), plus partial usage examples for DMSwarm-based field sampling.

Definition in file interpolation.c.

Macro Definition Documentation

◆ NUM_WEIGHTS

#define NUM_WEIGHTS   8

Definition at line 13 of file interpolation.c.

◆ ERROR_MSG_BUFFER_SIZE

#define ERROR_MSG_BUFFER_SIZE   256

Definition at line 16 of file interpolation.c.

◆ __FUNCT [1/6]

#define __FUNCT   "PiecWiseLinearInterpolation_Scalar"

Definition at line 527 of file interpolation.c.

◆ __FUNCT [2/6]

#define __FUNCT   "PiecWiseLinearInterpolation_Vector"

Definition at line 527 of file interpolation.c.

◆ __FUNCT [3/6]

#define __FUNCT   "ComputeTrilinearWeights"

Definition at line 527 of file interpolation.c.

◆ __FUNCT [4/6]

#define __FUNCT   "TrilinearInterpolation_Scalar"

Definition at line 527 of file interpolation.c.

◆ __FUNCT [5/6]

#define __FUNCT   "TrilinearInterpolation_Vector"

Definition at line 527 of file interpolation.c.

◆ __FUNCT [6/6]

#define __FUNCT   "InterpolateEulerFieldToSwarmForParticle"

Definition at line 527 of file interpolation.c.

Function Documentation

◆ InterpolateFieldFromCornerToCenter_Vector()

PetscErrorCode InterpolateFieldFromCornerToCenter_Vector ( Cmpnts ***  field_arr,
Cmpnts ***  centfield_arr,
UserCtx *  user 
)

Interpolates a vector field from corner nodes to cell centers, respecting application conventions.

Safely interpolate a vector field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.

For each interior computational cell (k,j,i), this function calculates the cell-centered value as the arithmetic average of the values at its 8 surrounding corner nodes.

This function adheres to the convention that the valid computational domain for cell-centered quantities excludes the first and last indices in each direction (e.g., cells at i=0 and i=IM). The loop bounds are constructed to match this, making it a direct counterpart to functions like ComputeCellCentersAndSpacing.

Parameters
[in]field_arrInput: Ghosted 3D array of corner-node data, accessed via GLOBAL indices.
[out]centfield_arrOutput: 3D array for interpolated cell-center values, accessed via GLOBAL indices.
[in]userUser context containing DMDA information.
Returns
PetscErrorCode 0 on success.

Definition at line 38 of file interpolation.c.

42{
43 PetscErrorCode ierr;
44 DMDALocalInfo info;
45
46 PetscFunctionBeginUser;
47
48 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
49
50 // Get local and global grid dimensions
51 PetscInt xs = info.xs, xe = info.xs + info.xm;
52 PetscInt ys = info.ys, ye = info.ys + info.ym;
53 PetscInt zs = info.zs, ze = info.zs + info.zm;
54 PetscInt mx = info.mx, my = info.my, mz = info.mz;
55
56 // Determine loop bounds to compute for interior cells only, matching the code's convention.
57 // Start at index 1 if the process owns the global boundary at index 0.
58 PetscInt is = (xs == 0) ? 1 : xs;
59 PetscInt js = (ys == 0) ? 1 : ys;
60 PetscInt ks = (zs == 0) ? 1 : zs;
61
62 // Stop one cell short if the process owns the global boundary at the max index.
63 PetscInt ie = (xe == mx) ? xe - 1 : xe;
64 PetscInt je = (ye == my) ? ye - 1 : ye;
65 PetscInt ke = (ze == mz) ? ze - 1 : ze;
66
67 // Loop over the locally owned INTERIOR cells.
68 for (PetscInt k = ks; k < ke; k++) {
69 for (PetscInt j = js; j < je; j++) {
70 for (PetscInt i = is; i < ie; i++) {
71 // Calculate cell center value as the average of its 8 corner nodes
72 centfield_arr[k][j][i].x = 0.125 * (field_arr[k][j][i].x + field_arr[k][j-1][i].x +
73 field_arr[k-1][j][i].x + field_arr[k-1][j-1][i].x +
74 field_arr[k][j][i-1].x + field_arr[k][j-1][i-1].x +
75 field_arr[k-1][j][i-1].x + field_arr[k-1][j-1][i-1].x);
76
77 centfield_arr[k][j][i].y = 0.125 * (field_arr[k][j][i].y + field_arr[k][j-1][i].y +
78 field_arr[k-1][j][i].y + field_arr[k-1][j-1][i].y +
79 field_arr[k][j][i-1].y + field_arr[k][j-1][i-1].y +
80 field_arr[k-1][j][i-1].y + field_arr[k-1][j-1][i-1].y);
81
82 centfield_arr[k][j][i].z = 0.125 * (field_arr[k][j][i].z + field_arr[k][j-1][i].z +
83 field_arr[k-1][j][i].z + field_arr[k-1][j-1][i].z +
84 field_arr[k][j][i-1].z + field_arr[k][j-1][i-1].z +
85 field_arr[k-1][j][i-1].z + field_arr[k-1][j-1][i-1].z);
86 }
87 }
88 }
89
90 PetscFunctionReturn(0);
91}
Cmpnts::x
PetscScalar x
Definition variables.h:101
Cmpnts::z
PetscScalar z
Definition variables.h:101
Cmpnts::y
PetscScalar y
Definition variables.h:101

◆ InterpolateFieldFromCornerToCenter_Scalar()

PetscErrorCode InterpolateFieldFromCornerToCenter_Scalar ( PetscReal ***  field_arr,
PetscReal ***  centfield_arr,
UserCtx *  user 
)

Interpolates a scalar field from corner nodes to cell centers, respecting application conventions.

Safely interpolate a scalar field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.

For each interior computational cell (k,j,i), this function calculates the cell-centered value as the arithmetic average of the values at its 8 surrounding corner nodes.

This function adheres to the convention that the valid computational domain for cell-centered quantities excludes the first and last indices in each direction (e.g., cells at i=0 and i=IM). The loop bounds are constructed to match this.

Parameters
[in]field_arrInput: Ghosted 3D array of corner-node data, accessed via GLOBAL indices.
[out]centfield_arrOutput: 3D array for interpolated cell-center values, accessed via GLOBAL indices.
[in]userUser context containing DMDA information.
Returns
PetscErrorCode 0 on success.

Definition at line 111 of file interpolation.c.

115{
116 PetscErrorCode ierr;
117 DMDALocalInfo info;
118
119 PetscFunctionBeginUser;
120
121 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
122
123 // Get local and global grid dimensions
124 PetscInt xs = info.xs, xe = info.xs + info.xm;
125 PetscInt ys = info.ys, ye = info.ys + info.ym;
126 PetscInt zs = info.zs, ze = info.zs + info.zm;
127 PetscInt mx = info.mx, my = info.my, mz = info.mz;
128
129 // Determine loop bounds to compute for interior cells only, matching the code's convention.
130 // Start at index 1 if the process owns the global boundary at index 0.
131 PetscInt is = (xs == 0) ? 1 : xs;
132 PetscInt js = (ys == 0) ? 1 : ys;
133 PetscInt ks = (zs == 0) ? 1 : zs;
134
135 // Stop one cell short if the process owns the global boundary at the max index.
136 PetscInt ie = (xe == mx) ? xe - 1 : xe;
137 PetscInt je = (ye == my) ? ye - 1 : ye;
138 PetscInt ke = (ze == mz) ? ze - 1 : ze;
139
140 // Loop over the locally owned INTERIOR cells.
141 for (PetscInt k = ks; k < ke; k++) {
142 for (PetscInt j = js; j < je; j++) {
143 for (PetscInt i = is; i < ie; i++) {
144 // Calculate cell center value as the average of its 8 corner nodes
145 centfield_arr[k][j][i] = 0.125 * (field_arr[k][j][i] + field_arr[k][j-1][i] +
146 field_arr[k-1][j][i] + field_arr[k-1][j-1][i] +
147 field_arr[k][j][i-1] + field_arr[k][j-1][i-1] +
148 field_arr[k-1][j][i-1] + field_arr[k-1][j-1][i-1]);
149 }
150 }
151 }
152
153 PetscFunctionReturn(0);
154}

◆ TestCornerToCenterInterpolation()

PetscErrorCode TestCornerToCenterInterpolation ( UserCtx *  user)

Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.

This function serves as a unit test. It performs the following steps:

  1. Takes the corner-centered nodal coordinates (from DMGetCoordinatesLocal) as input.
  2. Uses the InterpolateFieldFromCornerToCenter macro to interpolate these coordinates to the cell centers, storing the result in a new temporary vector.
  3. Compares this new vector with the user->Cent vector, which is assumed to have been computed by ComputeCellCentersAndSpacing and serves as the ground truth.
  4. A 2-norm of the difference is computed. If it is below a small tolerance, the test passes.
Note
This function should be called immediately after ComputeCellCentersAndSpacing has been successfully executed.
Parameters
userThe UserCtx for a specific grid level.
Returns
PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 175 of file interpolation.c.

176{
177 PetscErrorCode ierr;
178 Vec lCoords, TestCent;
179 Cmpnts ***coor_arr, ***test_cent_arr;
180 PetscReal diff_norm;
181
182 PetscFunctionBeginUser;
183
184 // 1. Create a temporary vector to hold the result of our interpolation.
185 // It must have the same layout and size as the ground-truth user->Cent vector.
186 ierr = VecDuplicate(user->Cent, &TestCent); CHKERRQ(ierr);
187
188 // 2. Get the input (corner coordinates) and output (our test vector) arrays.
189 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
190 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
191 ierr = DMDAVecGetArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
192
193 // 3. Call the generic interpolation macro.
194 // The macro will see that `coor_arr` is of type `Cmpnts***` and correctly
195 // call the `InterpolateFieldFromCornerToCenter_Vector` function.
196 ierr = InterpolateFieldFromCornerToCenter(coor_arr, test_cent_arr, user); CHKERRQ(ierr);
197
198 // 4. Restore the arrays.
199 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
200 ierr = DMDAVecRestoreArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
201
202 // 5. IMPORTANT: Assemble the vector so its values are communicated across processors
203 // and it's ready for global operations like VecNorm.
204 ierr = VecAssemblyBegin(TestCent); CHKERRQ(ierr);
205 ierr = VecAssemblyEnd(TestCent); CHKERRQ(ierr);
206
207 // 6. Compare the result with the ground truth.
208 // We compute TestCent = -1.0 * user->Cent + 1.0 * TestCent.
209 // This calculates the difference vector: TestCent - user->Cent.
210 ierr = VecAXPY(TestCent, -1.0, user->Cent); CHKERRQ(ierr);
211
212 // Now, compute the L2 norm of the difference vector. If the functions are
213 // identical, the norm should be zero (or very close due to floating point).
214 ierr = VecNorm(TestCent, NORM_2, &diff_norm); CHKERRQ(ierr);
215
216 // 7. Report the result and clean up.
217 if (diff_norm < 1.0e-12) {
218 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[SUCCESS] Test passed. Norm of difference is %g.\n", (double)diff_norm);
219 } else {
220 LOG_ALLOW(GLOBAL,LOG_DEBUG, "[FAILURE] Test failed. Norm of difference is %g.\n", (double)diff_norm);
221 }
222
223 ierr = VecDestroy(&TestCent); CHKERRQ(ierr);
224
225 PetscFunctionReturn(0);
226}
InterpolateFieldFromCornerToCenter
#define InterpolateFieldFromCornerToCenter(field, centfield, user)
Generic macro to call the appropriate interpolation function based on the field type.
Definition interpolation.h:36
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
UserCtx::Cent
Vec Cent
Definition variables.h:776
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100

◆ InterpolateFieldFromCenterToCorner_Vector()

PetscErrorCode InterpolateFieldFromCenterToCorner_Vector ( Cmpnts ***  centfield_arr,
Cmpnts ***  corner_arr,
UserCtx *  user 
)

Interpolates a vector field from cell centers to corner nodes.

Parameters
[in]centfield_arrInput: 3D array (ghosted) of cell-centered data, accessed via GLOBAL indices.
[out]corner_arrOutput: 3D array where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]userUser context containing DMDA information.
Returns
PetscErrorCode 0 on success.

Definition at line 241 of file interpolation.c.

245{
246 PetscErrorCode ierr;
247 DMDALocalInfo info;
248 PetscMPIInt rank;
249 PROFILE_FUNCTION_BEGIN;
250 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
251 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
252
253 // Node ownership range (GLOBAL indices)
254 PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
255 PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
256 PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
257
258 PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
259 PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
260 PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
261
262
263 // Global grid dimensions (used for valid cell check)
264 PetscInt IM = info.mx - 1; // Total nodes in i-direction
265 PetscInt JM = info.my - 1; // Total nodes in j-direction
266 PetscInt KM = info.mz - 1; // Total nodes in k-direction
267
268 // Ghosted array starting indices (from the same info struct)
269 PetscInt gxs = info.gxs;
270 PetscInt gys = info.gys;
271 PetscInt gzs = info.gzs;
272
273 PetscInt xs = info.xs;
274 PetscInt ys = info.ys;
275 PetscInt zs = info.zs;
276
277 LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
278 "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
279 rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
280
281 // Loop over the GLOBAL indices of the NODES owned by this processor
282 for (PetscInt k = zs_node; k < ze_node; k++) {
283 for (PetscInt j = ys_node; j < ye_node; j++) {
284 for (PetscInt i = xs_node; i < xe_node; i++) {
285 Cmpnts sum = {0.0, 0.0, 0.0};
286 PetscInt count = 0;
287
288 // DEBUG 1 TEST
289 /*
290 if(rank == 1 && i == 24 && j == 12 && k == 49){
291 PetscInt i_cell_A = i - 1;
292 PetscInt i_cell_B = i;
293
294 Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
295 Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
296
297 PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Read { Valid Cell }Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f) and {Out-of-Bounds Cell}Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
298 rank, k, j, i,
299 k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
300 k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
301
302 }
303 */
304
305 // Skip processing the unused last node in each dimension.
306 if(i >= IM || j >= JM || k >= KM){
307 continue;
308 }
309 // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
310 // The index offsets correspond to the relative position of the indices(shifted) that represent cell-centered field values of cells that share the node N(k,j,i) as a corner
311 for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
312 for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
313 for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
314
315 // These are still GLOBAL cell indices
316 PetscInt global_cell_k = k + dk_offset;
317 PetscInt global_cell_j = j + dj_offset;
318 PetscInt global_cell_i = i + di_offset;
319
320 // Check if this corresponds to a valid GLOBAL cell index
321 if (global_cell_i >= 0 && global_cell_i < nCellsX &&
322 global_cell_j >= 0 && global_cell_j < nCellsY &&
323 global_cell_k >= 0 && global_cell_k < nCellsZ)
324 {
325 Cmpnts cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
326
327 LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f, %.2f, %.2f)\n",
328 rank,global_cell_k,global_cell_j,global_cell_i,cell_val.x, cell_val.y, cell_val.z);
329
330 sum.x += cell_val.x;
331 sum.y += cell_val.y;
332 sum.z += cell_val.z;
333 count++;
334 }
335 }
336 }
337 }
338
339 PetscInt i_global_write = i; // Global index in GLOBAL array.
340 PetscInt j_global_write = j;
341 PetscInt k_global_write = k;
342
343 // We write directly into the array using the global loop indices.
344 if (count > 0) {
345 corner_arr[k_global_write][j_global_write][i_global_write].x = sum.x / (PetscReal)count;
346 corner_arr[k_global_write][j_global_write][i_global_write].y = sum.y / (PetscReal)count;
347 corner_arr[k_global_write][j_global_write][i_global_write].z = sum.z / (PetscReal)count;
348 } else {
349 // This case should ideally not happen for a valid owned node, but as a failsafe:
350 corner_arr[k_global_write][j_global_write][i_global_write] = (Cmpnts){0.0, 0.0, 0.0};
351 }
352
353 // DEBUG 2
354 /*
355 if(rank == 1){
356 if(i == 11 && j == 11 && k == 49){
357 Cmpnts ucat_node = corner_arr[k][j][i];
358 PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
359 rank, k, j, i,
360 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
361 }
362 }
363
364 if(rank == 0 && i == 24 && j == 12 && k == 0){
365 Cmpnts ucat_node = corner_arr[k][j][i];
366 PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
367 rank, k, j, i,
368 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
369 }
370 */
371 // LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,48,"[Rank %d] Node(k,j,i)=%d,%d,%d finished loops and write.\n", rank, k, j, i);
372 }
373 }
374 }
375 PROFILE_FUNCTION_END;
376 return 0;
377}
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:267
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:746
LOG_LOOP_ALLOW_EXACT
#define LOG_LOOP_ALLOW_EXACT(scope, level, var, val, fmt,...)
Logs a custom message if a variable equals a specific value.
Definition logging.h:348
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:34
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:737

◆ InterpolateFieldFromCenterToCorner_Scalar()

PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar ( PetscReal ***  centfield_arr,
PetscReal ***  corner_arr,
UserCtx *  user 
)

Interpolates a scalar field from cell centers to corner nodes.

Parameters
[in]centfield_arrInput: 3D array (ghosted) of cell-centered data, accessed via GLOBAL indices.
[out]corner_arrOutput: 3D array where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]userUser context containing DMDA information.
Returns
PetscErrorCode 0 on success.

Definition at line 392 of file interpolation.c.

396{
397 PetscErrorCode ierr;
398 DMDALocalInfo info;
399 PetscMPIInt rank;
400 PROFILE_FUNCTION_BEGIN;
401 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
402 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
403
404 // Node ownership range (GLOBAL indices)
405 PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
406 PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
407 PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
408
409 PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
410 PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
411 PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
412
413
414 // Global grid dimensions (used for valid cell check)
415 PetscInt IM = info.mx - 1; // Total nodes in i-direction
416 PetscInt JM = info.my - 1; // Total nodes in j-direction
417 PetscInt KM = info.mz - 1; // Total nodes in k-direction
418
419 // Ghosted array starting indices (from the same info struct)
420 PetscInt gxs = info.gxs;
421 PetscInt gys = info.gys;
422 PetscInt gzs = info.gzs;
423
424 PetscInt xs = info.xs;
425 PetscInt ys = info.ys;
426 PetscInt zs = info.zs;
427
428 LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
429 "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
430 rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
431
432 // Loop over the GLOBAL indices of the NODES owned by this processor
433 for (PetscInt k = zs_node; k < ze_node; k++) {
434 for (PetscInt j = ys_node; j < ye_node; j++) {
435 for (PetscInt i = xs_node; i < xe_node; i++) {
436 PetscReal sum = 0.0;
437 PetscInt count = 0;
438
439 // DEBUG 1 TEST
440 /*
441 if(rank == 1 && i == 24 && j == 12 && k == 49){
442 PetscInt i_cell_A = i - 1;
443 PetscInt i_cell_B = i;
444
445 Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
446 Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
447
448 PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Read { Valid Cell }Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f) and {Out-of-Bounds Cell}Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
449 rank, k, j, i,
450 k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
451 k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
452
453 }
454 */
455
456 // Skip processing the unused last node in each dimension.
457 if(i >= IM || j >= JM || k >= KM){
458 continue;
459 }
460 // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
461 // The index offsets correspond to the relative position of the indices(shifted) that represent cell-centered field values of cells that share the node N(k,j,i) as a corner
462 for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
463 for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
464 for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
465
466 // These are still GLOBAL cell indices
467 PetscInt global_cell_k = k + dk_offset;
468 PetscInt global_cell_j = j + dj_offset;
469 PetscInt global_cell_i = i + di_offset;
470
471 // Check if this corresponds to a valid GLOBAL cell index
472 if (global_cell_i >= 0 && global_cell_i < nCellsX &&
473 global_cell_j >= 0 && global_cell_j < nCellsY &&
474 global_cell_k >= 0 && global_cell_k < nCellsZ)
475 {
476 PetscReal cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
477
478 LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f)\n",
479 rank,global_cell_k,global_cell_j,global_cell_i,cell_val);
480
481 sum += cell_val;
482 count++;
483 }
484 }
485 }
486 }
487
488 PetscInt i_global_write = i; // Global index in GLOBAL array.
489 PetscInt j_global_write = j;
490 PetscInt k_global_write = k;
491
492 // We write directly into the array using the global loop indices.
493 if (count > 0) {
494 corner_arr[k_global_write][j_global_write][i_global_write] = sum / (PetscReal)count;
495 } else {
496 // This case should ideally not happen for a valid owned node, but as a failsafe:
497 corner_arr[k_global_write][j_global_write][i_global_write] = 0.0;
498 }
499
500 // DEBUG 2
501 /*
502 if(rank == 1){
503 if(i == 11 && j == 11 && k == 49){
504 Cmpnts ucat_node = corner_arr[k][j][i];
505 PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
506 rank, k, j, i,
507 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
508 }
509 }
510
511 if(rank == 0 && i == 24 && j == 12 && k == 0){
512 Cmpnts ucat_node = corner_arr[k][j][i];
513 PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
514 rank, k, j, i,
515 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
516 }
517 */
518 // LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,48,"[Rank %d] Node(k,j,i)=%d,%d,%d finished loops and write.\n", rank, k, j, i);
519 }
520 }
521 }
522 PROFILE_FUNCTION_END;
523 return 0;
524}

◆ PieceWiseLinearInterpolation_Scalar()

PetscErrorCode PieceWiseLinearInterpolation_Scalar ( const char *  fieldName,
PetscReal ***  fieldScal,
PetscInt  iCell,
PetscInt  jCell,
PetscInt  kCell,
PetscReal *  val 
)

Retrieves the scalar value at the cell (iCell, jCell, kCell).

This function does a simple piecewise (zeroth‐order) interpolation by returning fieldScal[kCell][jCell][iCell], ignoring any fractional coordinates.

Parameters
[in]fieldNameA string identifying the field (e.g. "temperature").
[in]fieldScal3D array of PetscReal (indexed as [k][j][i]).
[in]iCell,jCell,kCellIntegral cell indices to sample.
[out]valPointer to a PetscReal that receives the sampled scalar.
Returns
PetscErrorCode 0 on success

Definition at line 541 of file interpolation.c.

548{
549 PetscFunctionBegin;
550 *val = fieldScal[kCell][jCell][iCell];
551
552 // Optional logging
553 LOG_ALLOW(LOCAL, LOG_VERBOSE,
554 "Field '%s' at (i=%d, j=%d, k=%d) => val=%.6f\n",
555 fieldName, iCell, jCell, kCell, *val);
556
557 PetscFunctionReturn(0);
558}

◆ PieceWiseLinearInterpolation_Vector()

PetscErrorCode PieceWiseLinearInterpolation_Vector ( const char *  fieldName,
Cmpnts ***  fieldVec,
PetscInt  iCell,
PetscInt  jCell,
PetscInt  kCell,
Cmpnts *  vec 
)

Retrieves the vector (Cmpnts) at the cell (iCell, jCell, kCell).

This function simply sets: vec->x = fieldVec[kCell][jCell][iCell].x vec->y = fieldVec[kCell][jCell][iCell].y vec->z = fieldVec[kCell][jCell][iCell].z effectively a nearest‐cell or piecewise approach.

Parameters
[in]fieldNameString identifying the field (e.g. "velocity").
[in]fieldVec3D array of Cmpnts (indexed as [k][j][i]).
[in]iCell,jCell,kCellIntegral cell indices to sample.
[out]vecPointer to a Cmpnts struct that receives the sampled vector.
Returns
PetscErrorCode 0 on success

Definition at line 578 of file interpolation.c.

585{
586 PetscFunctionBegin;
587 vec->x = fieldVec[kCell][jCell][iCell].x;
588 vec->y = fieldVec[kCell][jCell][iCell].y;
589 vec->z = fieldVec[kCell][jCell][iCell].z;
590
591 // Optional logging
592 LOG_ALLOW(LOCAL, LOG_VERBOSE,
593 "Field '%s' at (i=%d, j=%d, k=%d) => (x=%.6f, y=%.6f, z=%.6f)\n",
594 fieldName, iCell, jCell, kCell, vec->x, vec->y, vec->z);
595
596 PetscFunctionReturn(0);
597}

◆ ComputeTrilinearWeights()

static void ComputeTrilinearWeights ( PetscReal  a1,
PetscReal  a2,
PetscReal  a3,
PetscReal *  w 
)
inlinestatic

Computes the trilinear interpolation weights from the interpolation coefficients.

This function computes the weights for trilinear interpolation at the eight corners of a cell using the interpolation coefficients provided along the x, y, and z directions.

Parameters
[in]a1Interpolation coefficient along the x-direction (normalized coordinate within the cell).
[in]a2Interpolation coefficient along the y-direction (normalized coordinate within the cell).
[in]a3Interpolation coefficient along the z-direction (normalized coordinate within the cell).
[out]wArray of 8 weights, each corresponding to one corner of the cell.
Note
  • The coefficients a1, a2, and a3 should be in the range [0, 1].
  • The order of weights corresponds to the eight corners of a hexahedral cell.

Definition at line 618 of file interpolation.c.

618 {
619 LOG_ALLOW(GLOBAL, LOG_VERBOSE, "Computing weights for a1=%f, a2=%f, a3=%f.\n", a1, a2, a3);
620
621 // Ensure a1, a2, a3 are within [0,1]
622 a1 = PetscMax(0.0, PetscMin(1.0, a1));
623 a2 = PetscMax(0.0, PetscMin(1.0, a2));
624 a3 = PetscMax(0.0, PetscMin(1.0, a3));
625
626 const PetscReal oa1 = 1.0 - a1;
627 const PetscReal oa2 = 1.0 - a2;
628 const PetscReal oa3 = 1.0 - a3;
629
630 w[0] = oa1 * oa2 * oa3; /* cornerOffsets[0] => (0,0,0) */
631 w[1] = a1 * oa2 * oa3; /* cornerOffsets[1] => (1,0,0) */
632 w[2] = oa1 * a2 * oa3; /* cornerOffsets[2] => (0,1,0) */
633 w[3] = a1 * a2 * oa3; /* cornerOffsets[3] => (1,1,0) */
634 w[4] = oa1 * oa2 * a3; /* cornerOffsets[4] => (0,0,1) */
635 w[5] = a1 * oa2 * a3; /* cornerOffsets[5] => (1,0,1) */
636 w[6] = oa1 * a2 * a3; /* cornerOffsets[6] => (0,1,1) */
637 w[7] = a1 * a2 * a3; /* cornerOffsets[7] => (1,1,1) */
638
639 // Log the computed weights for debugging
640 LOG_ALLOW(LOCAL,LOG_VERBOSE, "Weights computed - "
641 "w0=%f, w1=%f, w2=%f, w3=%f, w4=%f, w5=%f, w6=%f, w7=%f. \n",
642 w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7]);
643}
Here is the caller graph for this function:

◆ TrilinearInterpolation_Scalar()

PetscErrorCode TrilinearInterpolation_Scalar ( const char *  fieldName,
PetscReal ***  fieldScal,
PetscInt  i,
PetscInt  j,
PetscInt  k,
PetscReal  a1,
PetscReal  a2,
PetscReal  a3,
PetscReal *  val 
)

Computes the trilinear Interpolateted scalar at a given point.

Computes the trilinear interpolated scalar at a given point.

Parameters
[in]fieldNameA string representing the name of the scalar field (e.g., "temperature").
[in]fieldScal3D array of the field from a DMDA (indexed as [k][j][i]), each cell a PetscReal.
[in]i,j,kIntegral cell indices (the "lower" corner in each dimension).
[in]a1,a2,a3Normalized coordinates within the cell ([0,1] range).
[out]valPointer to a PetscReal that will store the Interpolateted scalar.

This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same PetscInterface.

Definition at line 661 of file interpolation.c.

671{
672 PetscFunctionBegin; // PETSc macro for error/stack tracing
673
674 // Compute the 8 corner weights
675 PetscReal wcorner[8];
676 ComputeTrilinearWeights(a1, a2, a3, wcorner);
677
678 // Offsets for cell corners
679 PetscInt i1 = i + 1;
680 PetscInt j1 = j + 1;
681 PetscInt k1 = k + 1;
682
683 // Initialize the output scalar
684 PetscReal sum = 0.0;
685
686 // Corner 0 => (i, j, k)
687 sum += wcorner[0] * fieldScal[k ][j ][i ];
688 // Corner 1 => (i+1, j, k)
689 sum += wcorner[1] * fieldScal[k ][j ][i1];
690 // Corner 2 => (i, j+1, k)
691 sum += wcorner[2] * fieldScal[k ][j1][i ];
692 // Corner 3 => (i+1, j+1, k)
693 sum += wcorner[3] * fieldScal[k ][j1][i1];
694 // Corner 4 => (i, j, k+1)
695 sum += wcorner[4] * fieldScal[k1][j ][i ];
696 // Corner 5 => (i+1, j, k+1)
697 sum += wcorner[5] * fieldScal[k1][j ][i1];
698 // Corner 6 => (i, j+1, k+1)
699 sum += wcorner[6] * fieldScal[k1][j1][i ];
700 // Corner 7 => (i+1, j+1, k+1)
701 sum += wcorner[7] * fieldScal[k1][j1][i1];
702
703 *val = sum;
704
705 // Logging (optional)
706 LOG_ALLOW(LOCAL, LOG_VERBOSE,
707 "Field '%s' at (i=%d, j=%d, k=%d), "
708 "a1=%.6f, a2=%.6f, a3=%.6f -> val=%.6f.\n",
709 fieldName, i, j, k, a1, a2, a3, *val);
710
711 // LOG_ALLOW_SYNC(GLOBAL, LOG_INFO,
712 // "TrilinearInterpolation_Scalar: Completed interpolation for field '%s' across local cells.\n",
713 // fieldName);
714
715 PetscFunctionReturn(0);
716}
ComputeTrilinearWeights
static void ComputeTrilinearWeights(PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
Computes the trilinear interpolation weights from the interpolation coefficients.
Definition interpolation.c:618
Here is the call graph for this function:

◆ TrilinearInterpolation_Vector()

PetscErrorCode TrilinearInterpolation_Vector ( const char *  fieldName,
Cmpnts ***  fieldVec,
PetscInt  i,
PetscInt  j,
PetscInt  k,
PetscReal  a1,
PetscReal  a2,
PetscReal  a3,
Cmpnts *  vec 
)

Computes the trilinear Interpolateted vector (e.g., velocity) at a given point.

Computes the trilinear interpolated vector (e.g., velocity) at a given point.

Parameters
[in]fieldNameA string representing the name of the vector field (e.g., "velocity").
[in]fieldVec3D array of the field from a DMDA (indexed as [k][j][i]), each cell of type Cmpnts.
[in]i,j,kIntegral cell indices (the "lower" corner in each dimension).
[in]a1,a2,a3Normalized coordinates within the cell ([0,1] range).
[out]vecPointer to a Cmpnts struct that will store the Interpolateted vector (x, y, z).

This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same PetscInterface.

Computes the trilinear Interpolateted vector at a given point, with partial weighting if corners go out of range.

If any of the 8 corners (i or i+1, j or j+1, k or k+1) is out of [0..mx), [0..my), [0..mz), that corner is skipped and the corresponding weight is omitted. The total is normalized by the sum of used weights, so we get a partial interpolation near boundaries.

Definition at line 741 of file interpolation.c.

751{
752 PetscFunctionBegin; // PETSc macro for error/stack tracing
753
754 // Compute the 8 corner weights
755 PetscErrorCode ierr;
756 PetscReal wcorner[8];
757 PetscMPIInt rank;
758
759 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
760
761 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Computing Trilinear Weights.\n",rank);
762 ComputeTrilinearWeights(a1, a2, a3, wcorner);
763
764 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Trilinear weights computed for local cell %d,%d,%d.\n",rank,i,j,k);
765
766 // For partial interpolation, we'll keep track of how many corners are valid
767 // and how much sum of weights is used. Then we do a final normalization.
768 PetscReal sumW = 0.0;
769 Cmpnts accum = {0.0, 0.0, 0.0};
770
771 // The eight corner indices, with their weights:
772 // corners: (i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), etc.
773 // We store them in an array to iterate cleanly.
774 const PetscInt cornerOffsets[8][3] = {
775 {0, 0, 0},
776 {1, 0, 0},
777 {0, 1, 0},
778 {1, 1, 0},
779 {0, 0, 1},
780 {1, 0, 1},
781 {0, 1, 1},
782 {1, 1, 1}
783 };
784
785 // Weighted partial sum
786 for (PetscInt c = 0; c < 8; c++) {
787 const PetscInt di = cornerOffsets[c][0];
788 const PetscInt dj = cornerOffsets[c][1];
789 const PetscInt dk = cornerOffsets[c][2];
790 PetscInt iC = i + di;
791 PetscInt jC = j + dj;
792 PetscInt kC = k + dk;
793
794 /*
795 // skip if out of domain
796 // (Assuming you know global domain is [0..mx), [0..my), [0..mz).)
797 if (iC < 0 || iC >= (PetscInt)userGlobalMx ||
798 jC < 0 || jC >= (PetscInt)userGlobalMy ||
799 kC < 0 || kC >= (PetscInt)userGlobalMz)
800 {
801 // skip this corner
802 continue;
803 }
804
805 */
806
807 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] %s[%d][%d][%d] = (%.4f,%.4f,%.4f).\n",rank,fieldName,kC,jC,iC,fieldVec[kC][jC][iC].x,fieldVec[kC][jC][iC].y,fieldVec[kC][jC][iC].z);
808
809 // Otherwise, accumulate
810 accum.x += wcorner[c] * fieldVec[kC][jC][iC].x;
811 accum.y += wcorner[c] * fieldVec[kC][jC][iC].y;
812 accum.z += wcorner[c] * fieldVec[kC][jC][iC].z;
813 sumW += wcorner[c];
814 }
815
816 // If sumW=0 => out-of-range or zero weighting => set (0,0,0)
817 if (sumW > 1.0e-14) {
818 vec->x = accum.x / sumW;
819 vec->y = accum.y / sumW;
820 vec->z = accum.z / sumW;
821 } else {
822 vec->x = 0.0; vec->y = 0.0; vec->z = 0.0;
823 }
824
825 PetscFunctionReturn(0);
826}
Here is the call graph for this function:

◆ InterpolateEulerFieldToSwarmForParticle()

static PetscErrorCode InterpolateEulerFieldToSwarmForParticle ( const char *  fieldName,
void *  fieldPtr,
Particle *  particle,
void *  swarmOut,
PetscInt  p,
PetscInt  blockSize 
)
inlinestatic

Interpolates a single field (scalar or vector) for one particle, storing the result in a swarm array.

This helper function is used inside a larger loop over local particles.

Parameters
[in]fieldNameA string identifying the field (e.g. "velocity", "temperature", etc.).
[in]fieldPtrA Pointer to the local DMDA array for the field:
  • (PetscReal ***) if scalar (blockSize = 1)
  • (Cmpnts ***) if vector (blockSize = 3)
[in]particleA Pointer to the Particle struct containing cell ID and weight information.
[out]swarmOutA Pointer to the DMSwarm output array:
  • (PetscReal*) if scalar dof=1
  • (PetscReal*) if vector dof=3 (storing x,y,z in consecutive reals)
  • or (Cmpnts*) if you store the result as a Cmpnts struct
[in]pThe local particle index (used to compute the correct offset in swarmOut).
[in]blockSizeThe number of degrees of freedom (1 => scalar, 3 => vector).

This routine demonstrates a switch between:

  • PiecewiseLinearInterpolation (zeroth order / nearest cell)
  • TrilinearInterpolation (8-corner weighted).

By default, PiecewiseLinearInterpolation is active, while the TrilinearInterpolation calls are commented out. To switch to trilinear, simply comment/uncomment appropriately.

Logging calls are provided (LOG_ALLOW) that you can adapt to your existing logging system.

Returns
PetscErrorCode 0 on success, non-zero on error.

Definition at line 859 of file interpolation.c.

866{
867 PetscErrorCode ierr;
868 PetscFunctionBegin;
869
870 PetscInt iCell = particle->cell[0];
871 PetscInt jCell = particle->cell[1];
872 PetscInt kCell = particle->cell[2];
873 PetscReal a1 = particle->weights.x;
874 PetscReal a2 = particle->weights.y;
875 PetscReal a3 = particle->weights.z;
876
877 PROFILE_FUNCTION_BEGIN;
878
879 // Optional logging at start
880 LOG_ALLOW(LOCAL, LOG_VERBOSE,
881 "field='%s', blockSize=%d, "
882 "cell IDs=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
883 fieldName, blockSize, iCell, jCell, kCell, a1, a2, a3);
884
885 /*
886 If blockSize=1, we PetscInterpret the fieldPtr as a 3D array of PetscReal (scalar).
887 If blockSize=3, we PetscInterpret the fieldPtr as a 3D array of Cmpnts (vector).
888 */
889 if (blockSize == 1) {
890 /* Scalar field: Cast fieldPtr to (PetscReal ***). */
891 PetscReal ***fieldScal = (PetscReal ***) fieldPtr;
892 PetscReal val;
893
894 // Currently using trilinear.
895 ierr = TrilinearInterpolation(fieldName, fieldScal,
896 iCell, jCell, kCell,
897 a1, a2, a3,
898 &val);
899 CHKERRQ(ierr);
900
901 // Alternative (commented) call to PiecewiseLinearInterpolation (zeroth order) :
902 // PetscErrorCode ierr = PiecewiseLinearInterpolation(fieldName,
903 // fieldScal,
904 // iCell, jCell, kCell,
905 // &val);
906 // CHKERRQ(ierr);
907
908 // Write the scalar result to the swarm output at index [p].
909 ((PetscReal*)swarmOut)[p] = val;
910
911 LOG_ALLOW(LOCAL, LOG_VERBOSE,
912 "field='%s', result=%.6f "
913 "stored at swarmOut index p=%d.\n", fieldName, val, (PetscInt)p);
914 }
915 else if (blockSize == 3) {
916 /* Vector field: Cast fieldPtr to (Cmpnts ***). */
917 Cmpnts ***fieldVec = (Cmpnts ***) fieldPtr;
918 Cmpnts vec;
919
920
921
922
923 // Piecewise interpolation (zeroth order).
924 // PetscErrorCode ierr = PieceWiseLinearInterpolation(fieldName,
925 // fieldVec,
926 // iCell, jCell, kCell,
927 // &vec);
928 // CHKERRQ(ierr);
929
930 // Alternative (commented) call to trilinear:
931 ierr = TrilinearInterpolation(fieldName, fieldVec,
932 iCell, jCell, kCell,
933 a1, a2, a3,
934 &vec);
935 CHKERRQ(ierr);
936
937 // If swarmOut is an array of 3 reals per particle:
938 ((PetscReal*)swarmOut)[3*p + 0] = vec.x;
939 ((PetscReal*)swarmOut)[3*p + 1] = vec.y;
940 ((PetscReal*)swarmOut)[3*p + 2] = vec.z;
941
942 LOG_ALLOW(LOCAL, LOG_VERBOSE,
943 "field='%s', result=(%.6f,%.6f,%.6f) "
944 "stored at swarmOut[3p..3p+2], p=%d.\n",
945 fieldName, vec.x, vec.y, vec.z, (PetscInt)p);
946
947 /*
948 If you store the vector result as a Cmpnts in the swarm, do instead:
949 ((Cmpnts*)swarmOut)[p] = vec;
950 but ensure your DMSwarm field is sized for a Cmpnts struct.
951 */
952 }
953 else {
954 /* If blockSize isn't 1 or 3, we raise an error. */
955 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
956 "InterpolateEulerFieldToSwarmForParticle: only blockSize=1 or 3 supported, got %d.",
957 (PetscInt)blockSize);
958 }
959
960 PROFILE_FUNCTION_END;
961 PetscFunctionReturn(0);
962}
TrilinearInterpolation
#define TrilinearInterpolation(fieldName, fieldPtr, i, j, k, a1, a2, a3, outPtr)
Macro that calls either the scalar or vector trilinear interpolation function based on the type of th...
Definition interpolation.h:122
Particle::cell
PetscInt cell[3]
Definition variables.h:167
Particle::weights
Cmpnts weights
Definition variables.h:170
Here is the caller graph for this function:

◆ InterpolateEulerFieldToSwarm()

PetscErrorCode InterpolateEulerFieldToSwarm ( UserCtx *  user,
Vec  fieldLocal_cellCentered,
const char *  fieldName,
const char *  swarmOutFieldName 
)

Interpolates a cell-centered field (scalar or vector) onto DMSwarm particles, using a robust, PETSc-idiomatic two-stage process.

This function first converts the cell-centered input data to corner-node data, storing this intermediate result in a PETSc Vec to correctly handle the communication of ghost-point information across parallel ranks. It then performs a final trilinear interpolation from the ghosted corner data to each particle's location.

Workflow:

  1. Create temporary PETSc Vecs (cornerGlobal, cornerLocal) to manage the intermediate corner-node data, using the existing nodal DMDA (user->fda).
  2. Call a dispatch macro that uses the runtime block size (bs) to select the correct underlying center-to-corner function, writing results into cornerGlobal.
  3. Perform a ghost-point exchange (DMGlobalToLocal) to transfer the boundary data from cornerGlobal into the ghost regions of cornerLocal.
  4. Loop over all local particles. For each particle: a. Convert its global cell index to a local index relative to the ghosted array. b. Check if the particle's interpolation stencil is fully contained within the owned+ghost region. If not, log a warning and set the result to zero. c. Perform the final trilinear interpolation using the ghosted cornerLocal data.
  5. Restore all PETSc objects to prevent memory leaks.
Parameters
[in]userUser context with DMDA, DMSwarm, etc.
[in]fieldlocal_cellCenteredVec (from fda) with cell-centered data (e.g., Ucat).
[in]fieldNameHuman-readable field name for logging (e.g., "Ucat").
[in]swarmOutFieldNameName of the DMSwarm field where interpolation results go.
Returns
PetscErrorCode 0 on success.

Definition at line 997 of file interpolation.c.

1002{
1003 PetscErrorCode ierr;
1004 DM fda = user->fda;
1005 DM swarm = user->swarm;
1006 PetscInt bs;
1007 DMDALocalInfo info;
1008 PetscMPIInt rank;
1009
1010 // Generic pointers to the raw data arrays
1011 void *cellCenterPtr_read;
1012 void *cornerPtr_read_with_ghosts;
1013
1014 // Swarm-related pointers
1015 PetscInt *cellIDs = NULL;
1016 PetscReal *weights = NULL;
1017 PetscInt64 *pids = NULL; // For logging particle IDs in warnings
1018 void *swarmOut = NULL;
1019 PetscReal *pos = NULL;
1020 PetscInt *status = NULL;
1021 PetscInt nLocal;
1022
1023 PetscFunctionBegin;
1024 PROFILE_FUNCTION_BEGIN;
1025 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1026 ierr = DMDAGetLocalInfo(fda, &info); CHKERRQ(ierr);
1027 ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
1028 if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
1029
1030 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting for field '%s'... \n", fieldName);
1031 // Select the appropriate DM for intermediate corner data based on block size
1032 DM dm_corner = (bs == 3) ? user->fda : user->da;
1033
1034 /* STAGE 1: Center-to-Corner Calculation with Communication */
1035
1036 // (A) Create or reuse the corner vectors.
1037 if(user->CellFieldAtCorner){
1038 PetscInt cached_bs;
1039 ierr = VecGetBlockSize(user->CellFieldAtCorner, &cached_bs); CHKERRQ(ierr);
1040 if(cached_bs != bs){
1041 LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d] Block size mismatch for existing corner vector and current field %s(Cached: %d, New: %d). Replacing.\n", rank, fieldName, cached_bs, bs);
1042 ierr = VecDestroy(&user->CellFieldAtCorner); CHKERRQ(ierr);
1043 ierr = VecDestroy(&user->lCellFieldAtCorner); CHKERRQ(ierr);
1044 user->CellFieldAtCorner = NULL;
1045 user->lCellFieldAtCorner = NULL;
1046 }
1047 }
1048
1049 if(user->CellFieldAtCorner == NULL){
1050 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new global corner vector for field '%s'.\n", rank, fieldName);
1051 ierr = DMCreateGlobalVector(dm_corner, &user->CellFieldAtCorner); CHKERRQ(ierr);
1052 }else{
1053 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing global corner vector for field '%s'.\n", rank, fieldName);
1054 }
1055
1056 ierr = VecSet(user->CellFieldAtCorner,0.0); CHKERRQ(ierr);
1057 Vec cornerGlobal = user->CellFieldAtCorner;
1058
1059 if(user->lCellFieldAtCorner == NULL){
1060 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new local corner vector for field '%s'.\n", rank, fieldName);
1061 ierr = DMCreateLocalVector(dm_corner, &user->lCellFieldAtCorner); CHKERRQ(ierr);
1062 }else{
1063 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing local corner vector for field '%s'.\n", rank, fieldName);
1064 }
1065
1066 ierr = VecSet(user->lCellFieldAtCorner,0.0); CHKERRQ(ierr);
1067 Vec cornerLocal = user->lCellFieldAtCorner;
1068
1069 // (B) Get a read-only array from the input cell-centered vector
1070 ierr = DMDAVecGetArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
1071
1072 PetscInt size;
1073 ierr = VecGetSize(cornerGlobal, &size); CHKERRQ(ierr);
1074 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1075
1076 size = 0;
1077
1078 ierr = VecGetSize(cornerLocal, &size); CHKERRQ(ierr);
1079 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1080
1081 size = 0;
1082
1083 ierr = VecGetSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
1084 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1085
1086 PetscInt xs,ys,zs,gxs,gys,gzs;
1087
1088 ierr = DMDAGetCorners(dm_corner,&xs,&ys,&zs,NULL,NULL,NULL); CHKERRQ(ierr);
1089 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] DMDAGetCorners for field '%s': xs=%d, ys=%d, zs=%d.\n", rank, fieldName, (PetscInt)xs, (PetscInt)ys, (PetscInt)zs);
1090
1091 ierr = DMDAGetGhostCorners(dm_corner,&gxs,&gys,&gzs,NULL,NULL,NULL); CHKERRQ(ierr);
1092 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] DMDAGetGhostCorners for field '%s': gxs=%d, gys=%d, gzs=%d.\n", rank, fieldName, (PetscInt)gxs, (PetscInt)gys, (PetscInt)gzs);
1093
1094 // DEBUG: Inspect Ucat ghost cells before interpolation
1095 /*
1096 if (bs == 3) {
1097 Cmpnts ***ucat_array = (Cmpnts***)cellCenterPtr_read;
1098 DMDALocalInfo info_debug;
1099 ierr = DMDAGetLocalInfo(fda, &info_debug); CHKERRQ(ierr);
1100
1101 // Only print on rank 0 to avoid clutter
1102 if (rank == 0) {
1103 PetscPrintf(PETSC_COMM_SELF, "\nDEBUG: Inspecting Ucat Ghost Cells...\n");
1104 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n");
1105
1106 // --- Check the MINIMUM-SIDE (-Xi) ---
1107 // The first physical cell is at index info.xs (local index).
1108 // The first ghost cell is at info.xs - 1.
1109 PetscInt i_first_phys = info_debug.xs;
1110 PetscInt i_first_ghost = info_debug.xs - 1;
1111 PetscInt j_mid = info_debug.ys + info_debug.ym / 2;
1112 PetscInt k_mid = info_debug.zs + info_debug.zm / 2;
1113
1114 PetscPrintf(PETSC_COMM_SELF, "MIN-SIDE (-Xi):\n");
1115 PetscPrintf(PETSC_COMM_SELF, " Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1116 k_mid, j_mid, i_first_ghost,
1117 ucat_array[k_mid][j_mid][i_first_ghost].x,
1118 ucat_array[k_mid][j_mid][i_first_ghost].y,
1119 ucat_array[k_mid][j_mid][i_first_ghost].z);
1120 PetscPrintf(PETSC_COMM_SELF, " First Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1121 k_mid, j_mid, i_first_phys,
1122 ucat_array[k_mid][j_mid][i_first_phys].x,
1123 ucat_array[k_mid][j_mid][i_first_phys].y,
1124 ucat_array[k_mid][j_mid][i_first_phys].z);
1125
1126
1127 // --- Check the MAXIMUM-SIDE (+Xi) ---
1128 // The last physical cell is at index info.xs + info.xm - 1.
1129 // The first ghost cell on the max side is at info.xs + info.xm.
1130 PetscInt i_last_phys = info_debug.xs + info_debug.xm - 1;
1131 PetscInt i_last_ghost = info_debug.xs + info_debug.xm;
1132
1133 PetscPrintf(PETSC_COMM_SELF, "MAX-SIDE (+Xi):\n");
1134 PetscPrintf(PETSC_COMM_SELF, " Last Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1135 k_mid, j_mid, i_last_phys,
1136 ucat_array[k_mid][j_mid][i_last_phys].x,
1137 ucat_array[k_mid][j_mid][i_last_phys].y,
1138 ucat_array[k_mid][j_mid][i_last_phys].z);
1139 PetscPrintf(PETSC_COMM_SELF, " Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1140 k_mid, j_mid, i_last_ghost,
1141 ucat_array[k_mid][j_mid][i_last_ghost].x,
1142 ucat_array[k_mid][j_mid][i_last_ghost].y,
1143 ucat_array[k_mid][j_mid][i_last_ghost].z);
1144 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n\n");
1145 }
1146 }
1147 */
1148 // DEBUG
1149
1150 // (C) Perform the center-to-corner interpolation directly into the global corner vector
1151 void *cornerPtr_write = NULL;
1152 ierr = DMDAVecGetArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
1153
1154 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Starting center-to-corner interpolation for '%s'.\n", rank, fieldName);
1155
1156 // SINGLE, CLEAN CALL SITE: The macro handles the runtime dispatch based on 'bs'.
1157 ierr = InterpolateFieldFromCenterToCorner(bs, cellCenterPtr_read, cornerPtr_write, user); CHKERRQ(ierr);
1158
1159 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Finished center-to-corner interpolation for '%s'.\n", rank, fieldName);
1160 ierr = DMDAVecRestoreArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
1161
1162 ierr = DMDAVecRestoreArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
1163
1164 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1165
1166 //////////// DEBUG
1167 /*
1168 // Log a synchronized header message from all ranks.
1169 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerGlobal BEFORE ghost exchange...\n");
1170
1171 // VecView prints to a PETSc viewer. PETSC_VIEWER_STDOUT_WORLD is a global, synchronized viewer.
1172 // We only need to call it if logging is active for this function.
1173 if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
1174 ierr = VecView(cornerGlobal, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1175 }
1176
1177 // Log a synchronized footer message.
1178 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerGlobal.\n");
1179 */
1180 ////////// DEBUG
1181
1182 // ierr = PetscBarrier((PetscObject)cornerGlobal); CHKERRQ(ierr);
1183
1184 // (D) CRITICAL STEP: Communicate the newly computed corner data to fill ghost regions
1185 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Beginning ghost exchange for corner data...\n", rank);
1186
1187 ierr = DMGlobalToLocalBegin(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
1188 ierr = DMGlobalToLocalEnd(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
1189
1190 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Ghost exchange for corner data complete.\n", rank);
1191
1192 LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE, "getting array from %s on all ranks.\n", dm_corner == user->fda ? "fda" : "da");
1193 if (is_function_allowed(__func__) && (int)(LOG_VERBOSE) <= (int)get_log_level()) {
1194 // DEBUG: Inspect cornerLocal after ghost exchange
1195
1196 ierr = LOG_FIELD_ANATOMY(user,"CornerField", "After Corner Velocity Interpolated"); CHKERRQ(ierr);
1197
1198 // DEBUG
1199 //ierr = DMView(user->fda, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1200 }
1201
1202 /////// DEBUG
1203 /*
1204 // Log a synchronized header message from all ranks.
1205 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerLocal AFTER ghost exchange...\n");
1206
1207 // Here, we want each rank to print its own local vector.
1208 // PETSC_VIEWER_STDOUT_SELF is the correct tool for this. The LOG_ALLOW_SYNC
1209 // wrapper will ensure the output from each rank is printed sequentially.
1210 if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
1211 PetscMPIInt rank_d;
1212 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_d); CHKERRQ(ierr);
1213
1214 // Use a synchronized print to create a clean header for each rank's output
1215 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n--- cornerLocal (Rank %d) ---\n", rank_d);
1216 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
1217
1218 // Print the local vector's contents for this rank
1219 ierr = VecView(cornerLocal, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
1220
1221 // Use a synchronized print to create a clean footer
1222 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--- End cornerLocal (Rank %d) ---\n", rank_d);
1223 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
1224 }
1225
1226 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerLocal.\n");
1227 */
1228 ///// DEBUG
1229
1230 // STAGE 2: Particle Interpolation using Ghosted Corner Data */
1231
1232 // (E) Get the local, GHOSTED array of corner data for the final interpolation step
1233 ierr = DMDAVecGetArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
1234
1235 // (F) Retrieve swarm fields for the particle loop
1236 ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
1237 ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1238 ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1239 ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1240 ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
1241 ierr = DMSwarmGetField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
1242 ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1243
1244 LOG_ALLOW(LOCAL,LOG_TRACE," Rank %d holds data upto & including %d,%d,%d.\n",rank,info.gxs + info.gxm,info.gys+info.gym,info.gzs+info.gzm);
1245
1246 // (G) Loop over each local particle
1247 for (PetscInt p = 0; p < nLocal; p++) {
1248
1249 Particle particle;
1250
1251 ierr = UnpackSwarmFields(p,pids,weights,pos,cellIDs,NULL,status,NULL,NULL,NULL,&particle); CHKERRQ(ierr);
1252
1253 LOG_ALLOW(LOCAL, LOG_VERBOSE,
1254 "[Rank %d] Particle PID %lld: global cell=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
1255 rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2],
1256 particle.weights.x, particle.weights.y, particle.weights.z);
1257 // Safety check: Ensure the entire 8-node stencil is within the valid memory
1258 // region (owned + ghosts) of the local array.
1259
1260 if (particle.cell[0] < info.gxs || particle.cell[0] >= info.gxs + info.gxm - 1 ||
1261 particle.cell[1] < info.gys || particle.cell[1] >= info.gys + info.gym - 1 ||
1262 particle.cell[2] < info.gzs || particle.cell[2] >= info.gzs + info.gzm - 1)
1263 {
1264 LOG_ALLOW(LOCAL, LOG_WARNING,
1265 "[Rank %d] Particle PID %lld in global cell (%d,%d,%d) is in an un-interpolatable region (requires ghosts of ghosts or is out of bounds). Zeroing field '%s'.\n",
1266 rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2], fieldName);
1267 if (bs == 3) {
1268 ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
1269 ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
1270 ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
1271 } else {
1272 ((PetscReal*)swarmOut)[p] = 0.0;
1273 }
1274 continue;
1275 }
1276
1277 ierr = InterpolateEulerFieldToSwarmForParticle(
1278 fieldName,
1279 cornerPtr_read_with_ghosts,
1280 &particle,
1281 swarmOut, p, bs);
1282 CHKERRQ(ierr);
1283 }
1284
1285 // (H) Restore all retrieved arrays and vectors */
1286 ierr = DMDAVecRestoreArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
1287
1288 // (I) Restore swarm fields
1289 ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1290 ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1291 ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1292 ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1293 ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
1294 ierr = DMSwarmRestoreField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
1295 LOG_ALLOW(GLOBAL, LOG_INFO, "Finished for field '%s'.\n", fieldName);
1296 PROFILE_FUNCTION_END;
1297 PetscFunctionReturn(0);
1298}
UnpackSwarmFields
PetscErrorCode UnpackSwarmFields(PetscInt i, const PetscInt64 *PIDs, const PetscReal *weights, const PetscReal *positions, const PetscInt *cellIndices, PetscReal *velocities, PetscInt *LocStatus, PetscReal *diffusivity, Cmpnts *diffusivitygradient, PetscReal *psi, Particle *particle)
Initializes a Particle struct with data from DMSwarm fields.
Definition ParticleSwarm.c:879
InterpolateEulerFieldToSwarmForParticle
static PetscErrorCode InterpolateEulerFieldToSwarmForParticle(const char *fieldName, void *fieldPtr, Particle *particle, void *swarmOut, PetscInt p, PetscInt blockSize)
Interpolates a single field (scalar or vector) for one particle, storing the result in a swarm array.
Definition interpolation.c:859
InterpolateFieldFromCenterToCorner
#define InterpolateFieldFromCenterToCorner(blockSize, centfield_ptr, corner_ptr, user_ctx)
Macro to dispatch to the correct scalar or vector center-to-corner function based on a runtime block ...
Definition interpolation.h:56
is_function_allowed
PetscBool is_function_allowed(const char *functionName)
Checks if a given function is in the allow-list.
Definition logging.c:157
LOG_FIELD_ANATOMY
PetscErrorCode LOG_FIELD_ANATOMY(UserCtx *user, const char *field_name, const char *stage_name)
Logs the anatomy of a specified field at key boundary locations, respecting the solver's specific gri...
Definition logging.c:1460
get_log_level
LogLevel get_log_level()
Retrieves the current logging level from the environment variable LOG_LEVEL.
Definition logging.c:39
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:29
UserCtx::lCellFieldAtCorner
Vec lCellFieldAtCorner
Definition variables.h:764
UserCtx::CellFieldAtCorner
Vec CellFieldAtCorner
Definition variables.h:764
Particle::PID
PetscInt64 PID
Definition variables.h:166
Particle
Defines a particle's core properties for Lagrangian tracking.
Definition variables.h:165
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateAllFieldsToSwarm()

PetscErrorCode InterpolateAllFieldsToSwarm ( UserCtx *  user)

Interpolates all relevant fields from the DMDA to the DMSwarm.

Currently, it Interpolatetes:

  • user->Ucat (vector field) into the DMSwarm field "swarmVelocity".

To add more fields, duplicate the call to InterpolateOneFieldOverSwarm and provide:

  • The global Vec for that field (e.g. user->Tcat for temperature),
  • A human-readable field name (for logging),
  • A DMSwarm output field name (e.g. "swarmTemperature").
Parameters
[in,out]userPointer to a UserCtx containing:
  • user->da (DM for the grid),
  • user->swarm (DMSwarm for particles),.
  • user->Ucat (Vec for the vector field),
  • possibly more fields like user->Tcat, user->Pcat, etc.
Returns
PetscErrorCode Returns 0 on success, non-zero on failure.

Definition at line 1321 of file interpolation.c.

1322{
1323 PetscErrorCode ierr;
1324 PetscMPIInt rank;
1325 PetscFunctionBegin;
1326
1327 PROFILE_FUNCTION_BEGIN;
1328
1329 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
1330
1331 /*
1332 1) Interpolate the 'velocity' field (user->Ucat) into DMSwarm's "swarmVelocity".
1333 - The function InterpolateOneFieldOverSwarm() loops over all local particles,
1334 retrieves iCell, jCell, kCell, (a1,a2,a3) from the DMSwarm,
1335 and calls the appropriate trilinear interpolation routine.
1336
1337 - "velocity" => just a label used in logging or debugging.
1338 - "swarmVelocity" => the DMSwarm field name where we store the result
1339 (assume dof=3 if user->Ucat has blockSize=3).
1340 */
1341
1342 LOG_ALLOW(LOCAL, LOG_DEBUG,
1343 " Interpolation of ucat to velocity begins on rank %d.\n",rank);
1344 // Make sure to pass the *LOCAL* Vector to the function below!
1345 ierr = InterpolateEulerFieldToSwarm(user, user->lUcat,
1346 "Ucat",
1347 "velocity"); CHKERRQ(ierr);
1348 ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivity,
1349 "Diffusivity", "Diffusivity"); CHKERRQ(ierr);
1350
1351 ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivityGradient,
1352 "DiffusivityGradient", "DiffusivityGradient"); CHKERRQ(ierr);
1353 /* Add fields as necessary here*/
1354
1355 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1356 LOG_ALLOW(LOCAL, LOG_INFO,
1357 "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
1358
1359 PROFILE_FUNCTION_END;
1360
1361 PetscFunctionReturn(0);
1362}
InterpolateEulerFieldToSwarm
PetscErrorCode InterpolateEulerFieldToSwarm(UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
Interpolates a cell-centered field (scalar or vector) onto DMSwarm particles, using a robust,...
Definition interpolation.c:997
UserCtx::lDiffusivityGradient
Vec lDiffusivityGradient
Definition variables.h:759
UserCtx::lUcat
Vec lUcat
Definition variables.h:755
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:758
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateCornerToFaceCenter_Scalar()

PetscErrorCode InterpolateCornerToFaceCenter_Scalar ( PetscReal ***  corner_arr,
PetscReal ***  faceX_arr,
PetscReal ***  faceY_arr,
PetscReal ***  faceZ_arr,
UserCtx *  user 
)

Interpolates a scalar field from corner nodes to all face centers.

This routine computes the average of the four corner-node values defining each face of a hexahedral cell:

  • X-faces (perpendicular to X): face between (i-1,i) in X-dir
  • Y-faces (perpendicular to Y): face between (j-1,j) in Y-dir
  • Z-faces (perpendicular to Z): face between (k-1,k) in Z-dir
Parameters
[in]corner_arrGhosted node-centered array (global indexing) from user->fda.
[out]faceX_arrLocal array for X-faces sized [zm][ym][xm+1].
[out]faceY_arrLocal array for Y-faces sized [zm][ym+1][xm].
[out]faceZ_arrLocal array for Z-faces sized [zm+1][ym][xm].
[in]userUser context containing DMDA 'fda' and GetOwnedCellRange.
Returns
PetscErrorCode 0 on success, non-zero on failure.

Definition at line 2228 of file interpolation.c.

2234{
2235 PetscErrorCode ierr;
2236 DMDALocalInfo info;
2237
2238 PetscFunctionBeginUser;
2239
2240 PROFILE_FUNCTION_BEGIN;
2241
2242 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2243
2244 // Determine owned-cell ranges based on corner-node ownership
2245 PetscInt xs, xm, ys, ym, zs, zm;
2246 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2247 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2248 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2249
2250 // Global exclusive end indices for cells
2251 PetscInt xe = xs + xm;
2252 PetscInt ye = ys + ym;
2253 PetscInt ze = zs + zm;
2254
2255 // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
2256 for (PetscInt k = zs; k < ze; ++k) {
2257 PetscInt k_loc = k - zs;
2258 for (PetscInt j = ys; j < ye; ++j) {
2259 PetscInt j_loc = j - ys;
2260 for (PetscInt i = xs; i <= xe; ++i) {
2261 PetscInt i_loc = i - xs; // 0..xm
2262 // Average the four corners of the Y-Z face at X = i
2263 PetscReal sum = corner_arr[k ][j ][i]
2264 + corner_arr[k+1][j ][i]
2265 + corner_arr[k ][j+1][i]
2266 + corner_arr[k+1][j+1][i];
2267 faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2268 }
2269 }
2270 }
2271
2272 // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
2273 for (PetscInt k = zs; k < ze; ++k) {
2274 PetscInt k_loc = k - zs;
2275 for (PetscInt j = ys; j <= ye; ++j) {
2276 PetscInt j_loc = j - ys; // 0..ym
2277 for (PetscInt i = xs; i < xe; ++i) {
2278 PetscInt i_loc = i - xs;
2279 // Average the four corners of the X-Z face at Y = j
2280 PetscReal sum = corner_arr[k ][j][i ]
2281 + corner_arr[k+1][j][i ]
2282 + corner_arr[k ][j][i+1]
2283 + corner_arr[k+1][j][i+1];
2284 faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2285 }
2286 }
2287 }
2288
2289 // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
2290 for (PetscInt k = zs; k <= ze; ++k) {
2291 PetscInt k_loc = k - zs;
2292 for (PetscInt j = ys; j < ye; ++j) {
2293 PetscInt j_loc = j - ys;
2294 for (PetscInt i = xs; i < xe; ++i) {
2295 PetscInt i_loc = i - xs;
2296 // Average the four corners of the X-Y face at Z = k
2297 PetscReal sum = corner_arr[k][j ][i ]
2298 + corner_arr[k][j ][i+1]
2299 + corner_arr[k][j+1][i ]
2300 + corner_arr[k][j+1][i+1];
2301 faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2302 }
2303 }
2304 }
2305
2306 PROFILE_FUNCTION_END;
2307
2308 PetscFunctionReturn(0);
2309}
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global, PetscInt *xm_cell_local)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:1756
Here is the call graph for this function:

◆ InterpolateCornerToFaceCenter_Vector()

PetscErrorCode InterpolateCornerToFaceCenter_Vector ( Cmpnts ***  corner_arr,
Cmpnts ***  faceX_arr,
Cmpnts ***  faceY_arr,
Cmpnts ***  faceZ_arr,
UserCtx *  user 
)

Interpolates a vector field from corner nodes to all face centers.

Identical to the scalar version, except it averages each component of the Cmpnts struct at the four corner-nodes per face.

Parameters
[in]corner_arrGhosted 3-component array (global node indices).
[out]faceX_arrLocal array of Cmpnts for X-faces sized [zm][ym][xm+1].
[out]faceY_arrLocal array of Cmpnts for Y-faces sized [zm][ym+1][xm].
[out]faceZ_arrLocal array of Cmpnts for Z-faces sized [zm+1][ym][xm].
[in]userUser context containing DMDA 'fda'.
Returns
PetscErrorCode 0 on success.

Definition at line 2328 of file interpolation.c.

2334{
2335 PetscErrorCode ierr;
2336 DMDALocalInfo info;
2337 PetscMPIInt rank;
2338
2339 PetscFunctionBeginUser;
2340
2341 PROFILE_FUNCTION_BEGIN;
2342
2343 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2344 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2345 "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
2346
2347 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2348
2349 PetscInt xs, xm, ys, ym, zs, zm;
2350 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2351 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2352 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2353
2354 PetscInt xe = xs + xm;
2355 PetscInt ye = ys + ym;
2356 PetscInt ze = zs + zm;
2357
2358 // X-faces
2359 for (PetscInt k = zs; k < ze; ++k) {
2360 PetscInt k_loc = k - zs;
2361 for (PetscInt j = ys; j < ye; ++j) {
2362 PetscInt j_loc = j - ys;
2363 for (PetscInt i = xs; i <= xe; ++i) {
2364 PetscInt i_loc = i - xs;
2365 Cmpnts sum = {0,0,0};
2366 sum.x = corner_arr[k ][j ][i].x + corner_arr[k+1][j ][i].x
2367 + corner_arr[k ][j+1][i].x + corner_arr[k+1][j+1][i].x;
2368 sum.y = corner_arr[k ][j ][i].y + corner_arr[k+1][j ][i].y
2369 + corner_arr[k ][j+1][i].y + corner_arr[k+1][j+1][i].y;
2370 sum.z = corner_arr[k ][j ][i].z + corner_arr[k+1][j ][i].z
2371 + corner_arr[k ][j+1][i].z + corner_arr[k+1][j+1][i].z;
2372 faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2373 faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2374 faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2375 }
2376 }
2377 }
2378
2379 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2380 "Rank %d x-face Interpolation complete.\n", rank);
2381
2382 // Y-faces
2383 for (PetscInt k = zs; k < ze; ++k) {
2384 PetscInt k_loc = k - zs;
2385 for (PetscInt j = ys; j <= ye; ++j) {
2386 PetscInt j_loc = j - ys;
2387 for (PetscInt i = xs; i < xe; ++i) {
2388 PetscInt i_loc = i - xs;
2389 Cmpnts sum = {0,0,0};
2390 sum.x = corner_arr[k ][j][i ].x + corner_arr[k+1][j][i ].x
2391 + corner_arr[k ][j][i+1].x + corner_arr[k+1][j][i+1].x;
2392 sum.y = corner_arr[k ][j][i ].y + corner_arr[k+1][j][i ].y
2393 + corner_arr[k ][j][i+1].y + corner_arr[k+1][j][i+1].y;
2394 sum.z = corner_arr[k ][j][i ].z + corner_arr[k+1][j][i ].z
2395 + corner_arr[k ][j][i+1].z + corner_arr[k+1][j][i+1].z;
2396 faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2397 faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2398 faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2399 }
2400 }
2401 }
2402
2403 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2404 "Rank %d y-face Interpolation complete.\n", rank);
2405
2406 // Z-faces
2407 for (PetscInt k = zs; k <= ze; ++k) {
2408 PetscInt k_loc = k - zs;
2409 for (PetscInt j = ys; j < ye; ++j) {
2410 PetscInt j_loc = j - ys;
2411 for (PetscInt i = xs; i < xe; ++i) {
2412 PetscInt i_loc = i - xs;
2413 Cmpnts sum = {0,0,0};
2414 sum.x = corner_arr[k][j ][i ].x + corner_arr[k][j ][i+1].x
2415 + corner_arr[k][j+1][i ].x + corner_arr[k][j+1][i+1].x;
2416 sum.y = corner_arr[k][j ][i ].y + corner_arr[k][j ][i+1].y
2417 + corner_arr[k][j+1][i ].y + corner_arr[k][j+1][i+1].y;
2418 sum.z = corner_arr[k][j ][i ].z + corner_arr[k][j ][i+1].z
2419 + corner_arr[k][j+1][i ].z + corner_arr[k][j+1][i+1].z;
2420 faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2421 faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2422 faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2423 }
2424 }
2425 }
2426
2427 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2428 "Rank %d z-face Interpolation complete.\n", rank);
2429
2430 PROFILE_FUNCTION_END;
2431 PetscFunctionReturn(0);
2432}
Here is the call graph for this function: