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"
 
#define __FUNCT__   "InterpolateFieldFromCornerToCenter_Scalar"
 
#define __FUNCT__   "TestCornerToCenterInterpolation"
 
#define __FUNCT__   "InterpolateFieldFromCenterToCorner_Vector"
 
#define __FUNCT__   "InterpolateFieldFromCenterToCorner_Scalar"
 
#define __FUNCT   "PiecWiseLinearInterpolation_Scalar"
 
#define __FUNCT   "PiecWiseLinearInterpolation_Vector"
 
#define __FUNCT   "ComputeTrilinearWeights"
 
#define __FUNCT   "TrilinearInterpolation_Scalar"
 
#define __FUNCT   "TrilinearInterpolation_Vector"
 
#define __FUNCT   "InterpolateEulerFieldToSwarmForParticle"
 
#define __FUNCT__   "InterpolateEulerFieldFromCenterToSwarm"
 
#define __FUNCT__   "InterpolateEulerFieldFromCornerToSwarm"
 
#define __FUNCT__   "InterpolateEulerFieldToSwarm"
 
#define __FUNCT__   "InterpolateAllFieldsToSwarm"
 
#define __FUNCT__   "GetScatterTargetInfo"
 
#define __FUNCT__   "GetPersistentLocalVector"
 
#define __FUNCT__   "AccumulateParticleField"
 
#define __FUNCT__   "NormalizeGridVectorByCount"
 
#define __FUNCT__   "ScatterParticleFieldToEulerField_Internal"
 
#define __FUNCT__   "ScatterParticleFieldToEulerField"
 
#define __FUNCT__   "ScatterAllParticleFieldsToEulerFields"
 
#define __FUNCT__   "InterpolateCornerToFaceCenter_Scalar"
 
#define __FUNCT__   "InterpolateCornerToFaceCenter_Vector"
 

Functions

PetscErrorCode InterpolateFieldFromCornerToCenter_Vector (Cmpnts ***field_arr, Cmpnts ***centfield_arr, UserCtx *user)
 Internal helper implementation: InterpolateFieldFromCornerToCenter_Vector().
 
PetscErrorCode InterpolateFieldFromCornerToCenter_Scalar (PetscReal ***field_arr, PetscReal ***centfield_arr, UserCtx *user)
 Internal helper implementation: InterpolateFieldFromCornerToCenter_Scalar().
 
PetscErrorCode TestCornerToCenterInterpolation (UserCtx *user)
 Internal helper implementation: TestCornerToCenterInterpolation().
 
PetscErrorCode InterpolateFieldFromCenterToCorner_Vector (Cmpnts ***centfield_arr, Cmpnts ***corner_arr, UserCtx *user)
 Internal helper implementation: InterpolateFieldFromCenterToCorner_Vector().
 
PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar (PetscReal ***centfield_arr, PetscReal ***corner_arr, UserCtx *user)
 Internal helper implementation: InterpolateFieldFromCenterToCorner_Scalar().
 
PetscErrorCode PieceWiseLinearInterpolation_Scalar (const char *fieldName, PetscReal ***fieldScal, PetscInt iCell, PetscInt jCell, PetscInt kCell, PetscReal *val)
 Internal helper implementation: PieceWiseLinearInterpolation_Scalar().
 
PetscErrorCode PieceWiseLinearInterpolation_Vector (const char *fieldName, Cmpnts ***fieldVec, PetscInt iCell, PetscInt jCell, PetscInt kCell, Cmpnts *vec)
 Internal helper implementation: PieceWiseLinearInterpolation_Vector().
 
static void ComputeTrilinearWeights (PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
 Internal helper implementation: ComputeTrilinearWeights().
 
static void ComputeTrilinearWeightsUnclamped (PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
 Unclamped trilinear weights for boundary extrapolation.
 
PetscErrorCode TrilinearInterpolation_Scalar (const char *fieldName, PetscReal ***fieldScal, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *val)
 Internal helper implementation: TrilinearInterpolation_Scalar().
 
PetscErrorCode TrilinearInterpolation_Vector (const char *fieldName, Cmpnts ***fieldVec, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, Cmpnts *vec)
 Internal helper implementation: TrilinearInterpolation_Vector().
 
static PetscErrorCode InterpolateEulerFieldToSwarmForParticle (const char *fieldName, void *fieldPtr, Particle *particle, void *swarmOut, PetscInt p, PetscInt blockSize)
 Internal helper implementation: InterpolateEulerFieldToSwarmForParticle().
 
static PetscErrorCode InterpolateEulerFieldFromCenterToSwarm (UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
 Direct cell-center trilinear interpolation (second-order on curvilinear grids).
 
static PetscErrorCode InterpolateEulerFieldFromCornerToSwarm (UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
 Corner-averaged interpolation path (legacy).
 
PetscErrorCode InterpolateEulerFieldToSwarm (UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
 Dispatches grid-to-particle interpolation to the method selected in the control file.
 
PetscErrorCode InterpolateAllFieldsToSwarm (UserCtx *user)
 Internal helper implementation: InterpolateAllFieldsToSwarm().
 
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)
 Internal helper implementation: GetPersistentLocalVector().
 
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 implementation: ScatterParticleFieldToEulerField_Internal().
 
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)
 Internal helper implementation: InterpolateCornerToFaceCenter_Scalar().
 
PetscErrorCode InterpolateCornerToFaceCenter_Vector (Cmpnts ***corner_arr, Cmpnts ***faceX_arr, Cmpnts ***faceY_arr, Cmpnts ***faceZ_arr, UserCtx *user)
 Internal helper implementation: InterpolateCornerToFaceCenter_Vector().
 

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 457 of file interpolation.c.

◆ __FUNCT [2/6]

#define __FUNCT   "PiecWiseLinearInterpolation_Vector"

Definition at line 457 of file interpolation.c.

◆ __FUNCT [3/6]

#define __FUNCT   "ComputeTrilinearWeights"

Definition at line 457 of file interpolation.c.

◆ __FUNCT [4/6]

#define __FUNCT   "TrilinearInterpolation_Scalar"

Definition at line 457 of file interpolation.c.

◆ __FUNCT [5/6]

#define __FUNCT   "TrilinearInterpolation_Vector"

Definition at line 457 of file interpolation.c.

◆ __FUNCT [6/6]

#define __FUNCT   "InterpolateEulerFieldToSwarmForParticle"

Definition at line 457 of file interpolation.c.

Function Documentation

◆ InterpolateFieldFromCornerToCenter_Vector()

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

Internal helper implementation: InterpolateFieldFromCornerToCenter_Vector().

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

Local to this translation unit.

Definition at line 25 of file interpolation.c.

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

◆ InterpolateFieldFromCornerToCenter_Scalar()

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

Internal helper implementation: InterpolateFieldFromCornerToCenter_Scalar().

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

Local to this translation unit.

Definition at line 86 of file interpolation.c.

90{
91 PetscErrorCode ierr;
92 DMDALocalInfo info;
93
94 PetscFunctionBeginUser;
95
96 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
97
98 // Get local and global grid dimensions
99 PetscInt xs = info.xs, xe = info.xs + info.xm;
100 PetscInt ys = info.ys, ye = info.ys + info.ym;
101 PetscInt zs = info.zs, ze = info.zs + info.zm;
102 PetscInt mx = info.mx, my = info.my, mz = info.mz;
103
104 // Determine loop bounds to compute for interior cells only, matching the code's convention.
105 // Start at index 1 if the process owns the global boundary at index 0.
106 PetscInt is = (xs == 0) ? 1 : xs;
107 PetscInt js = (ys == 0) ? 1 : ys;
108 PetscInt ks = (zs == 0) ? 1 : zs;
109
110 // Stop one cell short if the process owns the global boundary at the max index.
111 PetscInt ie = (xe == mx) ? xe - 1 : xe;
112 PetscInt je = (ye == my) ? ye - 1 : ye;
113 PetscInt ke = (ze == mz) ? ze - 1 : ze;
114
115 // Loop over the locally owned INTERIOR cells.
116 for (PetscInt k = ks; k < ke; k++) {
117 for (PetscInt j = js; j < je; j++) {
118 for (PetscInt i = is; i < ie; i++) {
119 // Calculate cell center value as the average of its 8 corner nodes
120 centfield_arr[k][j][i] = 0.125 * (field_arr[k][j][i] + field_arr[k][j-1][i] +
121 field_arr[k-1][j][i] + field_arr[k-1][j-1][i] +
122 field_arr[k][j][i-1] + field_arr[k][j-1][i-1] +
123 field_arr[k-1][j][i-1] + field_arr[k-1][j-1][i-1]);
124 }
125 }
126 }
127
128 PetscFunctionReturn(0);
129}
Here is the caller graph for this function:

◆ TestCornerToCenterInterpolation()

PetscErrorCode TestCornerToCenterInterpolation ( UserCtx *  user)

Internal helper implementation: TestCornerToCenterInterpolation().

Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.

Local to this translation unit.

Definition at line 137 of file interpolation.c.

138{
139 PetscErrorCode ierr;
140 Vec lCoords, TestCent;
141 Cmpnts ***coor_arr, ***test_cent_arr;
142 PetscReal diff_norm;
143
144 PetscFunctionBeginUser;
145
146 // 1. Create a temporary vector to hold the result of our interpolation.
147 // It must have the same layout and size as the ground-truth user->Cent vector.
148 ierr = VecDuplicate(user->Cent, &TestCent); CHKERRQ(ierr);
149
150 // 2. Get the input (corner coordinates) and output (our test vector) arrays.
151 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
152 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
153 ierr = DMDAVecGetArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
154
155 // 3. Call the generic interpolation macro.
156 // The macro will see that `coor_arr` is of type `Cmpnts***` and correctly
157 // call the `InterpolateFieldFromCornerToCenter_Vector` function.
158 ierr = InterpolateFieldFromCornerToCenter(coor_arr, test_cent_arr, user); CHKERRQ(ierr);
159
160 // 4. Restore the arrays.
161 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
162 ierr = DMDAVecRestoreArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
163
164 // 5. IMPORTANT: Assemble the vector so its values are communicated across processors
165 // and it's ready for global operations like VecNorm.
166 ierr = VecAssemblyBegin(TestCent); CHKERRQ(ierr);
167 ierr = VecAssemblyEnd(TestCent); CHKERRQ(ierr);
168
169 // 6. Compare the result with the ground truth.
170 // We compute TestCent = -1.0 * user->Cent + 1.0 * TestCent.
171 // This calculates the difference vector: TestCent - user->Cent.
172 ierr = VecAXPY(TestCent, -1.0, user->Cent); CHKERRQ(ierr);
173
174 // Now, compute the L2 norm of the difference vector. If the functions are
175 // identical, the norm should be zero (or very close due to floating point).
176 ierr = VecNorm(TestCent, NORM_2, &diff_norm); CHKERRQ(ierr);
177
178 // 7. Report the result and clean up.
179 if (diff_norm < 1.0e-12) {
180 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[SUCCESS] Test passed. Norm of difference is %g.\n", (double)diff_norm);
181 } else {
182 LOG_ALLOW(GLOBAL,LOG_DEBUG, "[FAILURE] Test failed. Norm of difference is %g.\n", (double)diff_norm);
183 }
184
185 ierr = VecDestroy(&TestCent); CHKERRQ(ierr);
186
187 PetscFunctionReturn(0);
188}
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:45
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
UserCtx::Cent
Vec Cent
Definition variables.h:858
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 
)

Internal helper implementation: InterpolateFieldFromCenterToCorner_Vector().

Interpolates a vector field from cell centers to corner nodes.

Local to this translation unit.

Definition at line 196 of file interpolation.c.

200{
201 PetscErrorCode ierr;
202 DMDALocalInfo info;
203 PetscMPIInt rank;
204 PROFILE_FUNCTION_BEGIN;
205 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
206 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
207
208 // Node ownership range (GLOBAL indices)
209 PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
210 PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
211 PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
212
213 PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
214 PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
215 PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
216
217
218 // Global grid dimensions (used for valid cell check)
219 PetscInt IM = info.mx - 1; // Total nodes in i-direction
220 PetscInt JM = info.my - 1; // Total nodes in j-direction
221 PetscInt KM = info.mz - 1; // Total nodes in k-direction
222
223 LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
224 "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
225 rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
226
227 // Loop over the GLOBAL indices of the NODES owned by this processor
228 for (PetscInt k = zs_node; k < ze_node; k++) {
229 for (PetscInt j = ys_node; j < ye_node; j++) {
230 for (PetscInt i = xs_node; i < xe_node; i++) {
231 Cmpnts sum = {0.0, 0.0, 0.0};
232 PetscInt count = 0;
233
234 // DEBUG 1 TEST
235 /*
236 if(rank == 1 && i == 24 && j == 12 && k == 49){
237 PetscInt i_cell_A = i - 1;
238 PetscInt i_cell_B = i;
239
240 Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
241 Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
242
243 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",
244 rank, k, j, i,
245 k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
246 k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
247
248 }
249 */
250
251 // Skip processing the unused last node in each dimension.
252 if(i >= IM || j >= JM || k >= KM){
253 continue;
254 }
255 // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
256 // 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
257 for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
258 for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
259 for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
260
261 // These are still GLOBAL cell indices
262 PetscInt global_cell_k = k + dk_offset;
263 PetscInt global_cell_j = j + dj_offset;
264 PetscInt global_cell_i = i + di_offset;
265
266 // Check if this corresponds to a valid GLOBAL cell index
267 if (global_cell_i >= 0 && global_cell_i < nCellsX &&
268 global_cell_j >= 0 && global_cell_j < nCellsY &&
269 global_cell_k >= 0 && global_cell_k < nCellsZ)
270 {
271 Cmpnts cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
272
273 LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f, %.2f, %.2f)\n",
274 rank,global_cell_k,global_cell_j,global_cell_i,cell_val.x, cell_val.y, cell_val.z);
275
276 sum.x += cell_val.x;
277 sum.y += cell_val.y;
278 sum.z += cell_val.z;
279 count++;
280 }
281 }
282 }
283 }
284
285 PetscInt i_global_write = i; // Global index in GLOBAL array.
286 PetscInt j_global_write = j;
287 PetscInt k_global_write = k;
288
289 // We write directly into the array using the global loop indices.
290 if (count > 0) {
291 corner_arr[k_global_write][j_global_write][i_global_write].x = sum.x / (PetscReal)count;
292 corner_arr[k_global_write][j_global_write][i_global_write].y = sum.y / (PetscReal)count;
293 corner_arr[k_global_write][j_global_write][i_global_write].z = sum.z / (PetscReal)count;
294 } else {
295 // This case should ideally not happen for a valid owned node, but as a failsafe:
296 corner_arr[k_global_write][j_global_write][i_global_write] = (Cmpnts){0.0, 0.0, 0.0};
297 }
298
299 // DEBUG 2
300 /*
301 if(rank == 1){
302 if(i == 11 && j == 11 && k == 49){
303 Cmpnts ucat_node = corner_arr[k][j][i];
304 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",
305 rank, k, j, i,
306 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
307 }
308 }
309
310 if(rank == 0 && i == 24 && j == 12 && k == 0){
311 Cmpnts ucat_node = corner_arr[k][j][i];
312 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",
313 rank, k, j, i,
314 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
315 }
316 */
317 // 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);
318 }
319 }
320 }
321 PROFILE_FUNCTION_END;
322 return 0;
323}
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:252
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
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:334
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770

◆ InterpolateFieldFromCenterToCorner_Scalar()

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

Internal helper implementation: InterpolateFieldFromCenterToCorner_Scalar().

Interpolates a scalar field from cell centers to corner nodes.

Local to this translation unit.

Definition at line 331 of file interpolation.c.

335{
336 PetscErrorCode ierr;
337 DMDALocalInfo info;
338 PetscMPIInt rank;
339 PROFILE_FUNCTION_BEGIN;
340 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
341 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
342
343 // Node ownership range (GLOBAL indices)
344 PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
345 PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
346 PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
347
348 PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
349 PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
350 PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
351
352
353 // Global grid dimensions (used for valid cell check)
354 PetscInt IM = info.mx - 1; // Total nodes in i-direction
355 PetscInt JM = info.my - 1; // Total nodes in j-direction
356 PetscInt KM = info.mz - 1; // Total nodes in k-direction
357
358 LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
359 "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
360 rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
361
362 // Loop over the GLOBAL indices of the NODES owned by this processor
363 for (PetscInt k = zs_node; k < ze_node; k++) {
364 for (PetscInt j = ys_node; j < ye_node; j++) {
365 for (PetscInt i = xs_node; i < xe_node; i++) {
366 PetscReal sum = 0.0;
367 PetscInt count = 0;
368
369 // DEBUG 1 TEST
370 /*
371 if(rank == 1 && i == 24 && j == 12 && k == 49){
372 PetscInt i_cell_A = i - 1;
373 PetscInt i_cell_B = i;
374
375 Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
376 Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
377
378 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",
379 rank, k, j, i,
380 k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
381 k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
382
383 }
384 */
385
386 // Skip processing the unused last node in each dimension.
387 if(i >= IM || j >= JM || k >= KM){
388 continue;
389 }
390 // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
391 // 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
392 for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
393 for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
394 for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
395
396 // These are still GLOBAL cell indices
397 PetscInt global_cell_k = k + dk_offset;
398 PetscInt global_cell_j = j + dj_offset;
399 PetscInt global_cell_i = i + di_offset;
400
401 // Check if this corresponds to a valid GLOBAL cell index
402 if (global_cell_i >= 0 && global_cell_i < nCellsX &&
403 global_cell_j >= 0 && global_cell_j < nCellsY &&
404 global_cell_k >= 0 && global_cell_k < nCellsZ)
405 {
406 PetscReal cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
407
408 LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f)\n",
409 rank,global_cell_k,global_cell_j,global_cell_i,cell_val);
410
411 sum += cell_val;
412 count++;
413 }
414 }
415 }
416 }
417
418 PetscInt i_global_write = i; // Global index in GLOBAL array.
419 PetscInt j_global_write = j;
420 PetscInt k_global_write = k;
421
422 // We write directly into the array using the global loop indices.
423 if (count > 0) {
424 corner_arr[k_global_write][j_global_write][i_global_write] = sum / (PetscReal)count;
425 } else {
426 // This case should ideally not happen for a valid owned node, but as a failsafe:
427 corner_arr[k_global_write][j_global_write][i_global_write] = 0.0;
428 }
429
430 // DEBUG 2
431 /*
432 if(rank == 1){
433 if(i == 11 && j == 11 && k == 49){
434 Cmpnts ucat_node = corner_arr[k][j][i];
435 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",
436 rank, k, j, i,
437 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
438 }
439 }
440
441 if(rank == 0 && i == 24 && j == 12 && k == 0){
442 Cmpnts ucat_node = corner_arr[k][j][i];
443 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",
444 rank, k, j, i,
445 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
446 }
447 */
448 // 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);
449 }
450 }
451 }
452 PROFILE_FUNCTION_END;
453 return 0;
454}

◆ PieceWiseLinearInterpolation_Scalar()

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

Internal helper implementation: PieceWiseLinearInterpolation_Scalar().

Returns the scalar value from the input cell index without blending.

Local to this translation unit.

Definition at line 462 of file interpolation.c.

469{
470 PetscFunctionBegin;
471 *val = fieldScal[kCell][jCell][iCell];
472
473 // Optional logging
474 LOG_ALLOW(LOCAL, LOG_VERBOSE,
475 "Field '%s' at (i=%d, j=%d, k=%d) => val=%.6f\n",
476 fieldName, iCell, jCell, kCell, *val);
477
478 PetscFunctionReturn(0);
479}

◆ PieceWiseLinearInterpolation_Vector()

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

Internal helper implementation: PieceWiseLinearInterpolation_Vector().

Returns the vector value from the input cell index without blending.

Local to this translation unit.

Definition at line 487 of file interpolation.c.

494{
495 PetscFunctionBegin;
496 vec->x = fieldVec[kCell][jCell][iCell].x;
497 vec->y = fieldVec[kCell][jCell][iCell].y;
498 vec->z = fieldVec[kCell][jCell][iCell].z;
499
500 // Optional logging
501 LOG_ALLOW(LOCAL, LOG_VERBOSE,
502 "Field '%s' at (i=%d, j=%d, k=%d) => (x=%.6f, y=%.6f, z=%.6f)\n",
503 fieldName, iCell, jCell, kCell, vec->x, vec->y, vec->z);
504
505 PetscFunctionReturn(0);
506}

◆ ComputeTrilinearWeights()

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

Internal helper implementation: ComputeTrilinearWeights().

Local to this translation unit.

Definition at line 516 of file interpolation.c.

516 {
517 LOG_ALLOW(GLOBAL, LOG_VERBOSE, "Computing weights for a1=%f, a2=%f, a3=%f.\n", a1, a2, a3);
518
519 // Ensure a1, a2, a3 are within [0,1]
520 a1 = PetscMax(0.0, PetscMin(1.0, a1));
521 a2 = PetscMax(0.0, PetscMin(1.0, a2));
522 a3 = PetscMax(0.0, PetscMin(1.0, a3));
523
524 const PetscReal oa1 = 1.0 - a1;
525 const PetscReal oa2 = 1.0 - a2;
526 const PetscReal oa3 = 1.0 - a3;
527
528 w[0] = oa1 * oa2 * oa3; /* cornerOffsets[0] => (0,0,0) */
529 w[1] = a1 * oa2 * oa3; /* cornerOffsets[1] => (1,0,0) */
530 w[2] = oa1 * a2 * oa3; /* cornerOffsets[2] => (0,1,0) */
531 w[3] = a1 * a2 * oa3; /* cornerOffsets[3] => (1,1,0) */
532 w[4] = oa1 * oa2 * a3; /* cornerOffsets[4] => (0,0,1) */
533 w[5] = a1 * oa2 * a3; /* cornerOffsets[5] => (1,0,1) */
534 w[6] = oa1 * a2 * a3; /* cornerOffsets[6] => (0,1,1) */
535 w[7] = a1 * a2 * a3; /* cornerOffsets[7] => (1,1,1) */
536
537 // Log the computed weights for debugging
538 LOG_ALLOW(LOCAL,LOG_VERBOSE, "Weights computed - "
539 "w0=%f, w1=%f, w2=%f, w3=%f, w4=%f, w5=%f, w6=%f, w7=%f. \n",
540 w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7]);
541}
Here is the caller graph for this function:

◆ ComputeTrilinearWeightsUnclamped()

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

Unclamped trilinear weights for boundary extrapolation.

Identical to ComputeTrilinearWeights() but without clamping a1, a2, a3 to [0,1]. Allows weights outside [0,1] for linear extrapolation at non-periodic boundaries. Weights still sum to 1.0.

Definition at line 549 of file interpolation.c.

549 {
550 LOG_ALLOW(GLOBAL, LOG_VERBOSE, "Computing unclamped weights for a1=%f, a2=%f, a3=%f.\n", a1, a2, a3);
551
552 const PetscReal oa1 = 1.0 - a1;
553 const PetscReal oa2 = 1.0 - a2;
554 const PetscReal oa3 = 1.0 - a3;
555
556 w[0] = oa1 * oa2 * oa3;
557 w[1] = a1 * oa2 * oa3;
558 w[2] = oa1 * a2 * oa3;
559 w[3] = a1 * a2 * oa3;
560 w[4] = oa1 * oa2 * a3;
561 w[5] = a1 * oa2 * a3;
562 w[6] = oa1 * a2 * a3;
563 w[7] = a1 * a2 * a3;
564
565 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Unclamped weights - "
566 "w0=%f, w1=%f, w2=%f, w3=%f, w4=%f, w5=%f, w6=%f, w7=%f.\n",
567 w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7]);
568}
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 
)

Internal helper implementation: TrilinearInterpolation_Scalar().

Computes the trilinear interpolated scalar at a given point.

Local to this translation unit.

Definition at line 577 of file interpolation.c.

587{
588 PetscFunctionBegin; // PETSc macro for error/stack tracing
589
590 // Compute the 8 corner weights
591 PetscReal wcorner[8];
592 ComputeTrilinearWeights(a1, a2, a3, wcorner);
593
594 // Offsets for cell corners
595 PetscInt i1 = i + 1;
596 PetscInt j1 = j + 1;
597 PetscInt k1 = k + 1;
598
599 // Initialize the output scalar
600 PetscReal sum = 0.0;
601
602 // Corner 0 => (i, j, k)
603 sum += wcorner[0] * fieldScal[k ][j ][i ];
604 // Corner 1 => (i+1, j, k)
605 sum += wcorner[1] * fieldScal[k ][j ][i1];
606 // Corner 2 => (i, j+1, k)
607 sum += wcorner[2] * fieldScal[k ][j1][i ];
608 // Corner 3 => (i+1, j+1, k)
609 sum += wcorner[3] * fieldScal[k ][j1][i1];
610 // Corner 4 => (i, j, k+1)
611 sum += wcorner[4] * fieldScal[k1][j ][i ];
612 // Corner 5 => (i+1, j, k+1)
613 sum += wcorner[5] * fieldScal[k1][j ][i1];
614 // Corner 6 => (i, j+1, k+1)
615 sum += wcorner[6] * fieldScal[k1][j1][i ];
616 // Corner 7 => (i+1, j+1, k+1)
617 sum += wcorner[7] * fieldScal[k1][j1][i1];
618
619 *val = sum;
620
621 // Logging (optional)
622 LOG_ALLOW(LOCAL, LOG_VERBOSE,
623 "Field '%s' at (i=%d, j=%d, k=%d), "
624 "a1=%.6f, a2=%.6f, a3=%.6f -> val=%.6f.\n",
625 fieldName, i, j, k, a1, a2, a3, *val);
626
627 // LOG_ALLOW_SYNC(GLOBAL, LOG_INFO,
628 // "TrilinearInterpolation_Scalar: Completed interpolation for field '%s' across local cells.\n",
629 // fieldName);
630
631 PetscFunctionReturn(0);
632}
ComputeTrilinearWeights
static void ComputeTrilinearWeights(PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
Internal helper implementation: ComputeTrilinearWeights().
Definition interpolation.c:516
Here is the call graph for this function:
Here is the caller 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 
)

Internal helper implementation: TrilinearInterpolation_Vector().

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

Local to this translation unit.

Definition at line 641 of file interpolation.c.

651{
652 PetscFunctionBegin; // PETSc macro for error/stack tracing
653
654 // Compute the 8 corner weights
655 PetscErrorCode ierr;
656 PetscReal wcorner[8];
657 PetscMPIInt rank;
658
659 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
660
661 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Computing Trilinear Weights.\n",rank);
662 ComputeTrilinearWeights(a1, a2, a3, wcorner);
663
664 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Trilinear weights computed for local cell %d,%d,%d.\n",rank,i,j,k);
665
666 // For partial interpolation, we'll keep track of how many corners are valid
667 // and how much sum of weights is used. Then we do a final normalization.
668 PetscReal sumW = 0.0;
669 Cmpnts accum = {0.0, 0.0, 0.0};
670
671 // The eight corner indices, with their weights:
672 // corners: (i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), etc.
673 // We store them in an array to iterate cleanly.
674 const PetscInt cornerOffsets[8][3] = {
675 {0, 0, 0},
676 {1, 0, 0},
677 {0, 1, 0},
678 {1, 1, 0},
679 {0, 0, 1},
680 {1, 0, 1},
681 {0, 1, 1},
682 {1, 1, 1}
683 };
684
685 // Weighted partial sum
686 for (PetscInt c = 0; c < 8; c++) {
687 const PetscInt di = cornerOffsets[c][0];
688 const PetscInt dj = cornerOffsets[c][1];
689 const PetscInt dk = cornerOffsets[c][2];
690 PetscInt iC = i + di;
691 PetscInt jC = j + dj;
692 PetscInt kC = k + dk;
693
694 /*
695 // skip if out of domain
696 // (Assuming you know global domain is [0..mx), [0..my), [0..mz).)
697 if (iC < 0 || iC >= (PetscInt)userGlobalMx ||
698 jC < 0 || jC >= (PetscInt)userGlobalMy ||
699 kC < 0 || kC >= (PetscInt)userGlobalMz)
700 {
701 // skip this corner
702 continue;
703 }
704
705 */
706
707 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);
708
709 // Otherwise, accumulate
710 accum.x += wcorner[c] * fieldVec[kC][jC][iC].x;
711 accum.y += wcorner[c] * fieldVec[kC][jC][iC].y;
712 accum.z += wcorner[c] * fieldVec[kC][jC][iC].z;
713 sumW += wcorner[c];
714 }
715
716 // If sumW=0 => out-of-range or zero weighting => set (0,0,0)
717 if (sumW > 1.0e-14) {
718 vec->x = accum.x / sumW;
719 vec->y = accum.y / sumW;
720 vec->z = accum.z / sumW;
721 } else {
722 vec->x = 0.0; vec->y = 0.0; vec->z = 0.0;
723 }
724
725 PetscFunctionReturn(0);
726}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateEulerFieldToSwarmForParticle()

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

Internal helper implementation: InterpolateEulerFieldToSwarmForParticle().

Local to this translation unit.

Definition at line 735 of file interpolation.c.

742{
743 PetscErrorCode ierr;
744 PetscFunctionBegin;
745
746 PetscInt iCell = particle->cell[0];
747 PetscInt jCell = particle->cell[1];
748 PetscInt kCell = particle->cell[2];
749 PetscReal a1 = particle->weights.x;
750 PetscReal a2 = particle->weights.y;
751 PetscReal a3 = particle->weights.z;
752
753 PROFILE_FUNCTION_BEGIN;
754
755 // Optional logging at start
756 LOG_ALLOW(LOCAL, LOG_VERBOSE,
757 "field='%s', blockSize=%d, "
758 "cell IDs=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
759 fieldName, blockSize, iCell, jCell, kCell, a1, a2, a3);
760
761 /*
762 If blockSize=1, we PetscInterpret the fieldPtr as a 3D array of PetscReal (scalar).
763 If blockSize=3, we PetscInterpret the fieldPtr as a 3D array of Cmpnts (vector).
764 */
765 if (blockSize == 1) {
766 /* Scalar field: Cast fieldPtr to (PetscReal ***). */
767 PetscReal ***fieldScal = (PetscReal ***) fieldPtr;
768 PetscReal val;
769
770 // Currently using trilinear.
771 ierr = TrilinearInterpolation(fieldName, fieldScal,
772 iCell, jCell, kCell,
773 a1, a2, a3,
774 &val);
775 CHKERRQ(ierr);
776
777 // Alternative (commented) call to PiecewiseLinearInterpolation (zeroth order) :
778 // PetscErrorCode ierr = PiecewiseLinearInterpolation(fieldName,
779 // fieldScal,
780 // iCell, jCell, kCell,
781 // &val);
782 // CHKERRQ(ierr);
783
784 // Write the scalar result to the swarm output at index [p].
785 ((PetscReal*)swarmOut)[p] = val;
786
787 LOG_ALLOW(LOCAL, LOG_VERBOSE,
788 "field='%s', result=%.6f "
789 "stored at swarmOut index p=%d.\n", fieldName, val, (PetscInt)p);
790 }
791 else if (blockSize == 3) {
792 /* Vector field: Cast fieldPtr to (Cmpnts ***). */
793 Cmpnts ***fieldVec = (Cmpnts ***) fieldPtr;
794 Cmpnts vec;
795
796
797
798
799 // Piecewise interpolation (zeroth order).
800 // PetscErrorCode ierr = PieceWiseLinearInterpolation(fieldName,
801 // fieldVec,
802 // iCell, jCell, kCell,
803 // &vec);
804 // CHKERRQ(ierr);
805
806 // Alternative (commented) call to trilinear:
807 ierr = TrilinearInterpolation(fieldName, fieldVec,
808 iCell, jCell, kCell,
809 a1, a2, a3,
810 &vec);
811 CHKERRQ(ierr);
812
813 // If swarmOut is an array of 3 reals per particle:
814 ((PetscReal*)swarmOut)[3*p + 0] = vec.x;
815 ((PetscReal*)swarmOut)[3*p + 1] = vec.y;
816 ((PetscReal*)swarmOut)[3*p + 2] = vec.z;
817
818 LOG_ALLOW(LOCAL, LOG_VERBOSE,
819 "field='%s', result=(%.6f,%.6f,%.6f) "
820 "stored at swarmOut[3p..3p+2], p=%d.\n",
821 fieldName, vec.x, vec.y, vec.z, (PetscInt)p);
822
823 /*
824 If you store the vector result as a Cmpnts in the swarm, do instead:
825 ((Cmpnts*)swarmOut)[p] = vec;
826 but ensure your DMSwarm field is sized for a Cmpnts struct.
827 */
828 }
829 else {
830 /* If blockSize isn't 1 or 3, we raise an error. */
831 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
832 "InterpolateEulerFieldToSwarmForParticle: only blockSize=1 or 3 supported, got %d.",
833 (PetscInt)blockSize);
834 }
835
836 PROFILE_FUNCTION_END;
837 PetscFunctionReturn(0);
838}
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:166
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:

◆ InterpolateEulerFieldFromCenterToSwarm()

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

Direct cell-center trilinear interpolation (second-order on curvilinear grids).

For each particle, determines the 8 nearest cell centers via octant detection, constructs a dual cell, computes face-distance-based trilinear weights, and interpolates directly from cell-centered data. At non-periodic boundaries where the dual cell cannot be fully formed, octant clamping with unclamped trilinear extrapolation preserves second-order accuracy. Periodic boundaries use ghost cell data directly. No intermediate corner staging or extra ghost exchange is needed.

Definition at line 853 of file interpolation.c.

858{
859 PetscErrorCode ierr;
860 DM swarm = user->swarm;
861 PetscInt bs;
862 DMDALocalInfo info;
863 PetscMPIInt rank;
864 void *fieldPtr = NULL;
865 Cmpnts ***cent = NULL;
866
867 PetscInt *cellIDs = NULL;
868 PetscReal *weights = NULL;
869 PetscInt64 *pids = NULL;
870 void *swarmOut = NULL;
871 PetscReal *pos = NULL;
872 PetscInt *status = NULL;
873 PetscInt nLocal;
874
875 PetscFunctionBegin;
876 PROFILE_FUNCTION_BEGIN;
877 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
878 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
879 ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
880 if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
881
882 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting direct cell-center interpolation for field '%s'...\n", fieldName);
883
884 /* Number of physical cells in each direction */
885 PetscInt nCellsX = info.mx - 2;
886 PetscInt nCellsY = info.my - 2;
887 PetscInt nCellsZ = info.mz - 2;
888
889 /* Periodic flags per direction */
890 PetscBool x_periodic = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
891 PetscBool y_periodic = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
892 PetscBool z_periodic = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
893
894 /* Get read-only arrays for the cell-centered field and cell center coordinates */
895 DM dm_field = (bs == 3) ? user->fda : user->da;
896 ierr = DMDAVecGetArrayRead(dm_field, fieldLocal_cellCentered, &fieldPtr); CHKERRQ(ierr);
897 ierr = DMDAVecGetArrayRead(user->fda, user->lCent, (void *)&cent); CHKERRQ(ierr);
898
899 /* Retrieve swarm fields */
900 ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
901 ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
902 ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
903 ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
904 ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
905 ierr = DMSwarmGetField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
906 ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
907
908 /* Loop over each local particle */
909 for (PetscInt p = 0; p < nLocal; p++) {
910
911 Particle particle;
912 ierr = UnpackSwarmFields(p, pids, weights, pos, cellIDs, NULL, status, NULL, NULL, NULL, &particle); CHKERRQ(ierr);
913
914 /* Step 1: Get host cell indices and existing weights */
915 PetscInt ci = particle.cell[0];
916 PetscInt cj = particle.cell[1];
917 PetscInt ck = particle.cell[2];
918 PetscReal a1 = particle.weights.x;
919 PetscReal a2 = particle.weights.y;
920 PetscReal a3 = particle.weights.z;
921
922 /* Step 2: Determine octant */
923 PetscInt oi = (a1 < 0.5) ? -1 : 0;
924 PetscInt oj = (a2 < 0.5) ? -1 : 0;
925 PetscInt ok = (a3 < 0.5) ? -1 : 0;
926 PetscInt bi = ci + oi;
927 PetscInt bj = cj + oj;
928 PetscInt bk = ck + ok;
929
930 /* Step 3: Clamp octant at non-periodic boundaries */
931 PetscBool needs_unclamped = PETSC_FALSE;
932
933 if (!x_periodic) {
934 if (bi < 0) { bi = 0; needs_unclamped = PETSC_TRUE; }
935 if (bi + 1 >= nCellsX) { bi = nCellsX - 2; needs_unclamped = PETSC_TRUE; }
936 }
937 if (!y_periodic) {
938 if (bj < 0) { bj = 0; needs_unclamped = PETSC_TRUE; }
939 if (bj + 1 >= nCellsY) { bj = nCellsY - 2; needs_unclamped = PETSC_TRUE; }
940 }
941 if (!z_periodic) {
942 if (bk < 0) { bk = 0; needs_unclamped = PETSC_TRUE; }
943 if (bk + 1 >= nCellsZ) { bk = nCellsZ - 2; needs_unclamped = PETSC_TRUE; }
944 }
945
946 /* Step 4: Bounds check — dual cell stencil must fit in ghosted region */
947 if (bi + 1 < info.gxs || bi + 1 >= info.gxs + info.gxm - 1 ||
948 bj + 1 < info.gys || bj + 1 >= info.gys + info.gym - 1 ||
949 bk + 1 < info.gzs || bk + 1 >= info.gzs + info.gzm - 1)
950 {
951 LOG_ALLOW(LOCAL, LOG_WARNING,
952 "[Rank %d] Particle PID %lld: dual cell (%d,%d,%d) out of ghosted region. Zeroing '%s'.\n",
953 rank, (long long)particle.PID, bi, bj, bk, fieldName);
954 if (bs == 3) {
955 ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
956 ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
957 ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
958 } else {
959 ((PetscReal*)swarmOut)[p] = 0.0;
960 }
961 continue;
962 }
963
964 /* Step 5: Construct dual cell from lCent (shifted index: cent[k+1][j+1][i+1] for cell (i,j,k))
965 * Vertex ordering follows GetCellVerticesFromGrid() convention (walkingsearch.c:423-430). */
966 Cell dual_cell;
967 dual_cell.vertices[0] = cent[bk + 1][bj + 1][bi + 1]; /* (bi, bj, bk ) */
968 dual_cell.vertices[1] = cent[bk + 1][bj + 1][bi + 2]; /* (bi+1, bj, bk ) */
969 dual_cell.vertices[2] = cent[bk + 1][bj + 2][bi + 2]; /* (bi+1, bj+1, bk ) */
970 dual_cell.vertices[3] = cent[bk + 1][bj + 2][bi + 1]; /* (bi, bj+1, bk ) */
971 dual_cell.vertices[4] = cent[bk + 2][bj + 2][bi + 1]; /* (bi, bj+1, bk+1) */
972 dual_cell.vertices[5] = cent[bk + 2][bj + 2][bi + 2]; /* (bi+1, bj+1, bk+1) */
973 dual_cell.vertices[6] = cent[bk + 2][bj + 1][bi + 2]; /* (bi+1, bj, bk+1) */
974 dual_cell.vertices[7] = cent[bk + 2][bj + 1][bi + 1]; /* (bi, bj, bk+1) */
975
976 /* Step 6: Compute face distances from particle to dual cell faces */
977 PetscReal d[NUM_FACES];
978 ierr = CalculateDistancesToCellFaces(particle.loc, &dual_cell, d, 1e-11); CHKERRQ(ierr);
979
980 /* Step 7: Compute trilinear parametric coordinates from face distances */
981 PetscReal a1_new = d[LEFT] / (d[LEFT] + d[RIGHT]);
982 PetscReal a2_new = d[BOTTOM] / (d[BOTTOM] + d[TOP]);
983 PetscReal a3_new = d[BACK] / (d[FRONT] + d[BACK]);
984
985 LOG_ALLOW(LOCAL, LOG_VERBOSE,
986 "[Rank %d] PID %lld: dual base=(%d,%d,%d), new weights=(%.4f,%.4f,%.4f), unclamped=%d\n",
987 rank, (long long)particle.PID, bi, bj, bk, a1_new, a2_new, a3_new, (int)needs_unclamped);
988
989 /* Step 8: Call TrilinearInterpolation directly with shifted dual-cell indices.
990 * Shifted index: pass (bi+1, bj+1, bk+1) so the kernel accesses
991 * field[bk+1..bk+2][bj+1..bj+2][bi+1..bi+2] = cells (bi..bi+1, bj..bj+1, bk..bk+1). */
992 if (bs == 1) {
993 PetscReal ***fieldScal = (PetscReal ***)fieldPtr;
994 PetscReal val;
995 if (needs_unclamped) {
996 PetscReal w[8];
997 ComputeTrilinearWeightsUnclamped(a1_new, a2_new, a3_new, w);
998 val = w[0] * fieldScal[bk+1][bj+1][bi+1] + w[1] * fieldScal[bk+1][bj+1][bi+2]
999 + w[2] * fieldScal[bk+1][bj+2][bi+1] + w[3] * fieldScal[bk+1][bj+2][bi+2]
1000 + w[4] * fieldScal[bk+2][bj+1][bi+1] + w[5] * fieldScal[bk+2][bj+1][bi+2]
1001 + w[6] * fieldScal[bk+2][bj+2][bi+1] + w[7] * fieldScal[bk+2][bj+2][bi+2];
1002 } else {
1003 ierr = TrilinearInterpolation_Scalar(fieldName, fieldScal,
1004 bi + 1, bj + 1, bk + 1, a1_new, a2_new, a3_new, &val); CHKERRQ(ierr);
1005 }
1006 ((PetscReal*)swarmOut)[p] = val;
1007 } else {
1008 Cmpnts ***fieldVec = (Cmpnts ***)fieldPtr;
1009 Cmpnts vec;
1010 if (needs_unclamped) {
1011 PetscReal w[8];
1012 ComputeTrilinearWeightsUnclamped(a1_new, a2_new, a3_new, w);
1013 vec.x = w[0]*fieldVec[bk+1][bj+1][bi+1].x + w[1]*fieldVec[bk+1][bj+1][bi+2].x
1014 + w[2]*fieldVec[bk+1][bj+2][bi+1].x + w[3]*fieldVec[bk+1][bj+2][bi+2].x
1015 + w[4]*fieldVec[bk+2][bj+1][bi+1].x + w[5]*fieldVec[bk+2][bj+1][bi+2].x
1016 + w[6]*fieldVec[bk+2][bj+2][bi+1].x + w[7]*fieldVec[bk+2][bj+2][bi+2].x;
1017 vec.y = w[0]*fieldVec[bk+1][bj+1][bi+1].y + w[1]*fieldVec[bk+1][bj+1][bi+2].y
1018 + w[2]*fieldVec[bk+1][bj+2][bi+1].y + w[3]*fieldVec[bk+1][bj+2][bi+2].y
1019 + w[4]*fieldVec[bk+2][bj+1][bi+1].y + w[5]*fieldVec[bk+2][bj+1][bi+2].y
1020 + w[6]*fieldVec[bk+2][bj+2][bi+1].y + w[7]*fieldVec[bk+2][bj+2][bi+2].y;
1021 vec.z = w[0]*fieldVec[bk+1][bj+1][bi+1].z + w[1]*fieldVec[bk+1][bj+1][bi+2].z
1022 + w[2]*fieldVec[bk+1][bj+2][bi+1].z + w[3]*fieldVec[bk+1][bj+2][bi+2].z
1023 + w[4]*fieldVec[bk+2][bj+1][bi+1].z + w[5]*fieldVec[bk+2][bj+1][bi+2].z
1024 + w[6]*fieldVec[bk+2][bj+2][bi+1].z + w[7]*fieldVec[bk+2][bj+2][bi+2].z;
1025 } else {
1026 ierr = TrilinearInterpolation_Vector(fieldName, fieldVec,
1027 bi + 1, bj + 1, bk + 1, a1_new, a2_new, a3_new, &vec); CHKERRQ(ierr);
1028 }
1029 ((PetscReal*)swarmOut)[3*p + 0] = vec.x;
1030 ((PetscReal*)swarmOut)[3*p + 1] = vec.y;
1031 ((PetscReal*)swarmOut)[3*p + 2] = vec.z;
1032 }
1033 }
1034
1035 /* Restore arrays and swarm fields */
1036 ierr = DMDAVecRestoreArrayRead(dm_field, fieldLocal_cellCentered, &fieldPtr); CHKERRQ(ierr);
1037 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCent, (void *)&cent); CHKERRQ(ierr);
1038 ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1039 ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1040 ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1041 ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1042 ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
1043 ierr = DMSwarmRestoreField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
1044
1045 LOG_ALLOW(GLOBAL, LOG_INFO, "Finished direct cell-center interpolation for field '%s'.\n", fieldName);
1046 PROFILE_FUNCTION_END;
1047 PetscFunctionReturn(0);
1048}
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:722
TrilinearInterpolation_Vector
PetscErrorCode TrilinearInterpolation_Vector(const char *fieldName, Cmpnts ***fieldVec, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, Cmpnts *vec)
Internal helper implementation: TrilinearInterpolation_Vector().
Definition interpolation.c:641
TrilinearInterpolation_Scalar
PetscErrorCode TrilinearInterpolation_Scalar(const char *fieldName, PetscReal ***fieldScal, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *val)
Internal helper implementation: TrilinearInterpolation_Scalar().
Definition interpolation.c:577
ComputeTrilinearWeightsUnclamped
static void ComputeTrilinearWeightsUnclamped(PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
Unclamped trilinear weights for boundary extrapolation.
Definition interpolation.c:549
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
UserCtx::lCent
Vec lCent
Definition variables.h:858
PERIODIC
@ PERIODIC
Definition variables.h:260
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:829
Particle::loc
Cmpnts loc
Definition variables.h:168
TOP
@ TOP
Definition variables.h:145
FRONT
@ FRONT
Definition variables.h:145
BOTTOM
@ BOTTOM
Definition variables.h:145
BACK
@ BACK
Definition variables.h:145
LEFT
@ LEFT
Definition variables.h:145
NUM_FACES
@ NUM_FACES
Definition variables.h:145
RIGHT
@ RIGHT
Definition variables.h:145
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
Particle::PID
PetscInt64 PID
Definition variables.h:166
Cell::vertices
Cmpnts vertices[8]
Coordinates of the eight vertices of the cell.
Definition variables.h:161
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
Cell
Defines the vertices of a single hexahedral grid cell.
Definition variables.h:160
Particle
Defines a particle's core properties for Lagrangian tracking.
Definition variables.h:165
CalculateDistancesToCellFaces
PetscErrorCode CalculateDistancesToCellFaces(const Cmpnts p, const Cell *cell, PetscReal *d, const PetscReal threshold)
Computes the signed distances from a point to each face of a cubic cell.
Definition walkingsearch.c:206
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateEulerFieldFromCornerToSwarm()

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

Corner-averaged interpolation path (legacy).

Stages cell-centered data to corners via unweighted averaging, then trilinear interpolation from corners to particles. Local to this translation unit.

Definition at line 1059 of file interpolation.c.

1064{
1065 PetscErrorCode ierr;
1066 DM fda = user->fda;
1067 DM swarm = user->swarm;
1068 PetscInt bs;
1069 DMDALocalInfo info;
1070 PetscMPIInt rank;
1071
1072 // Generic pointers to the raw data arrays
1073 void *cellCenterPtr_read;
1074 void *cornerPtr_read_with_ghosts;
1075
1076 // Swarm-related pointers
1077 PetscInt *cellIDs = NULL;
1078 PetscReal *weights = NULL;
1079 PetscInt64 *pids = NULL; // For logging particle IDs in warnings
1080 void *swarmOut = NULL;
1081 PetscReal *pos = NULL;
1082 PetscInt *status = NULL;
1083 PetscInt nLocal;
1084
1085 PetscFunctionBegin;
1086 PROFILE_FUNCTION_BEGIN;
1087 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1088 ierr = DMDAGetLocalInfo(fda, &info); CHKERRQ(ierr);
1089 ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
1090 if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
1091
1092 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting for field '%s'... \n", fieldName);
1093 // Select the appropriate DM for intermediate corner data based on block size
1094 DM dm_corner = (bs == 3) ? user->fda : user->da;
1095
1096 /* STAGE 1: Center-to-Corner Calculation with Communication */
1097
1098 // (A) Create or reuse the corner vectors.
1099 if(user->CellFieldAtCorner){
1100 PetscInt cached_bs;
1101 ierr = VecGetBlockSize(user->CellFieldAtCorner, &cached_bs); CHKERRQ(ierr);
1102 if(cached_bs != bs){
1103 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);
1104 ierr = VecDestroy(&user->CellFieldAtCorner); CHKERRQ(ierr);
1105 ierr = VecDestroy(&user->lCellFieldAtCorner); CHKERRQ(ierr);
1106 user->CellFieldAtCorner = NULL;
1107 user->lCellFieldAtCorner = NULL;
1108 }
1109 }
1110
1111 if(user->CellFieldAtCorner == NULL){
1112 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new global corner vector for field '%s'.\n", rank, fieldName);
1113 ierr = DMCreateGlobalVector(dm_corner, &user->CellFieldAtCorner); CHKERRQ(ierr);
1114 }else{
1115 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing global corner vector for field '%s'.\n", rank, fieldName);
1116 }
1117
1118 ierr = VecSet(user->CellFieldAtCorner,0.0); CHKERRQ(ierr);
1119 Vec cornerGlobal = user->CellFieldAtCorner;
1120
1121 if(user->lCellFieldAtCorner == NULL){
1122 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new local corner vector for field '%s'.\n", rank, fieldName);
1123 ierr = DMCreateLocalVector(dm_corner, &user->lCellFieldAtCorner); CHKERRQ(ierr);
1124 }else{
1125 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing local corner vector for field '%s'.\n", rank, fieldName);
1126 }
1127
1128 ierr = VecSet(user->lCellFieldAtCorner,0.0); CHKERRQ(ierr);
1129 Vec cornerLocal = user->lCellFieldAtCorner;
1130
1131 // (B) Get a read-only array from the input cell-centered vector
1132 ierr = DMDAVecGetArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
1133
1134 PetscInt size;
1135 ierr = VecGetSize(cornerGlobal, &size); CHKERRQ(ierr);
1136 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1137
1138 size = 0;
1139
1140 ierr = VecGetSize(cornerLocal, &size); CHKERRQ(ierr);
1141 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1142
1143 size = 0;
1144
1145 ierr = VecGetSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
1146 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1147
1148 PetscInt xs,ys,zs,gxs,gys,gzs;
1149
1150 ierr = DMDAGetCorners(dm_corner,&xs,&ys,&zs,NULL,NULL,NULL); CHKERRQ(ierr);
1151 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);
1152
1153 ierr = DMDAGetGhostCorners(dm_corner,&gxs,&gys,&gzs,NULL,NULL,NULL); CHKERRQ(ierr);
1154 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);
1155
1156 // DEBUG: Inspect Ucat ghost cells before interpolation
1157 /*
1158 if (bs == 3) {
1159 Cmpnts ***ucat_array = (Cmpnts***)cellCenterPtr_read;
1160 DMDALocalInfo info_debug;
1161 ierr = DMDAGetLocalInfo(fda, &info_debug); CHKERRQ(ierr);
1162
1163 // Only print on rank 0 to avoid clutter
1164 if (rank == 0) {
1165 PetscPrintf(PETSC_COMM_SELF, "\nDEBUG: Inspecting Ucat Ghost Cells...\n");
1166 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n");
1167
1168 // --- Check the MINIMUM-SIDE (-Xi) ---
1169 // The first physical cell is at index info.xs (local index).
1170 // The first ghost cell is at info.xs - 1.
1171 PetscInt i_first_phys = info_debug.xs;
1172 PetscInt i_first_ghost = info_debug.xs - 1;
1173 PetscInt j_mid = info_debug.ys + info_debug.ym / 2;
1174 PetscInt k_mid = info_debug.zs + info_debug.zm / 2;
1175
1176 PetscPrintf(PETSC_COMM_SELF, "MIN-SIDE (-Xi):\n");
1177 PetscPrintf(PETSC_COMM_SELF, " Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1178 k_mid, j_mid, i_first_ghost,
1179 ucat_array[k_mid][j_mid][i_first_ghost].x,
1180 ucat_array[k_mid][j_mid][i_first_ghost].y,
1181 ucat_array[k_mid][j_mid][i_first_ghost].z);
1182 PetscPrintf(PETSC_COMM_SELF, " First Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1183 k_mid, j_mid, i_first_phys,
1184 ucat_array[k_mid][j_mid][i_first_phys].x,
1185 ucat_array[k_mid][j_mid][i_first_phys].y,
1186 ucat_array[k_mid][j_mid][i_first_phys].z);
1187
1188
1189 // --- Check the MAXIMUM-SIDE (+Xi) ---
1190 // The last physical cell is at index info.xs + info.xm - 1.
1191 // The first ghost cell on the max side is at info.xs + info.xm.
1192 PetscInt i_last_phys = info_debug.xs + info_debug.xm - 1;
1193 PetscInt i_last_ghost = info_debug.xs + info_debug.xm;
1194
1195 PetscPrintf(PETSC_COMM_SELF, "MAX-SIDE (+Xi):\n");
1196 PetscPrintf(PETSC_COMM_SELF, " Last Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1197 k_mid, j_mid, i_last_phys,
1198 ucat_array[k_mid][j_mid][i_last_phys].x,
1199 ucat_array[k_mid][j_mid][i_last_phys].y,
1200 ucat_array[k_mid][j_mid][i_last_phys].z);
1201 PetscPrintf(PETSC_COMM_SELF, " Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1202 k_mid, j_mid, i_last_ghost,
1203 ucat_array[k_mid][j_mid][i_last_ghost].x,
1204 ucat_array[k_mid][j_mid][i_last_ghost].y,
1205 ucat_array[k_mid][j_mid][i_last_ghost].z);
1206 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n\n");
1207 }
1208 }
1209 */
1210 // DEBUG
1211
1212 // (C) Perform the center-to-corner interpolation directly into the global corner vector
1213 void *cornerPtr_write = NULL;
1214 ierr = DMDAVecGetArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
1215
1216 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Starting center-to-corner interpolation for '%s'.\n", rank, fieldName);
1217
1218 // SINGLE, CLEAN CALL SITE: The macro handles the runtime dispatch based on 'bs'.
1219 ierr = InterpolateFieldFromCenterToCorner(bs, cellCenterPtr_read, cornerPtr_write, user); CHKERRQ(ierr);
1220
1221 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Finished center-to-corner interpolation for '%s'.\n", rank, fieldName);
1222 ierr = DMDAVecRestoreArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
1223
1224 ierr = DMDAVecRestoreArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
1225
1226 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1227
1228 //////////// DEBUG
1229 /*
1230 // Log a synchronized header message from all ranks.
1231 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerGlobal BEFORE ghost exchange...\n");
1232
1233 // VecView prints to a PETSc viewer. PETSC_VIEWER_STDOUT_WORLD is a global, synchronized viewer.
1234 // We only need to call it if logging is active for this function.
1235 if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
1236 ierr = VecView(cornerGlobal, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1237 }
1238
1239 // Log a synchronized footer message.
1240 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerGlobal.\n");
1241 */
1242 ////////// DEBUG
1243
1244 // ierr = PetscBarrier((PetscObject)cornerGlobal); CHKERRQ(ierr);
1245
1246 // (D) CRITICAL STEP: Communicate the newly computed corner data to fill ghost regions
1247 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Beginning ghost exchange for corner data...\n", rank);
1248
1249 ierr = DMGlobalToLocalBegin(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
1250 ierr = DMGlobalToLocalEnd(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
1251
1252 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Ghost exchange for corner data complete.\n", rank);
1253
1254 LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE, "getting array from %s on all ranks.\n", dm_corner == user->fda ? "fda" : "da");
1255 if (is_function_allowed(__func__) && (int)(LOG_VERBOSE) <= (int)get_log_level()) {
1256 // DEBUG: Inspect cornerLocal after ghost exchange
1257
1258 ierr = LOG_FIELD_ANATOMY(user,"CornerField", "After Corner Velocity Interpolated"); CHKERRQ(ierr);
1259
1260 // DEBUG
1261 //ierr = DMView(user->fda, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1262 }
1263
1264 /////// DEBUG
1265 /*
1266 // Log a synchronized header message from all ranks.
1267 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerLocal AFTER ghost exchange...\n");
1268
1269 // Here, we want each rank to print its own local vector.
1270 // PETSC_VIEWER_STDOUT_SELF is the correct tool for this. The LOG_ALLOW_SYNC
1271 // wrapper will ensure the output from each rank is printed sequentially.
1272 if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
1273 PetscMPIInt rank_d;
1274 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_d); CHKERRQ(ierr);
1275
1276 // Use a synchronized print to create a clean header for each rank's output
1277 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n--- cornerLocal (Rank %d) ---\n", rank_d);
1278 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
1279
1280 // Print the local vector's contents for this rank
1281 ierr = VecView(cornerLocal, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
1282
1283 // Use a synchronized print to create a clean footer
1284 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--- End cornerLocal (Rank %d) ---\n", rank_d);
1285 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
1286 }
1287
1288 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerLocal.\n");
1289 */
1290 ///// DEBUG
1291
1292 // STAGE 2: Particle Interpolation using Ghosted Corner Data */
1293
1294 // (E) Get the local, GHOSTED array of corner data for the final interpolation step
1295 ierr = DMDAVecGetArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
1296
1297 // (F) Retrieve swarm fields for the particle loop
1298 ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
1299 ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1300 ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1301 ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1302 ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
1303 ierr = DMSwarmGetField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
1304 ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1305
1306 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);
1307
1308 // (G) Loop over each local particle
1309 for (PetscInt p = 0; p < nLocal; p++) {
1310
1311 Particle particle;
1312
1313 ierr = UnpackSwarmFields(p,pids,weights,pos,cellIDs,NULL,status,NULL,NULL,NULL,&particle); CHKERRQ(ierr);
1314
1315 LOG_ALLOW(LOCAL, LOG_VERBOSE,
1316 "[Rank %d] Particle PID %lld: global cell=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
1317 rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2],
1318 particle.weights.x, particle.weights.y, particle.weights.z);
1319 // Safety check: Ensure the entire 8-node stencil is within the valid memory
1320 // region (owned + ghosts) of the local array.
1321
1322 if (particle.cell[0] < info.gxs || particle.cell[0] >= info.gxs + info.gxm - 1 ||
1323 particle.cell[1] < info.gys || particle.cell[1] >= info.gys + info.gym - 1 ||
1324 particle.cell[2] < info.gzs || particle.cell[2] >= info.gzs + info.gzm - 1)
1325 {
1326 LOG_ALLOW(LOCAL, LOG_WARNING,
1327 "[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",
1328 rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2], fieldName);
1329 if (bs == 3) {
1330 ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
1331 ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
1332 ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
1333 } else {
1334 ((PetscReal*)swarmOut)[p] = 0.0;
1335 }
1336 continue;
1337 }
1338
1339 ierr = InterpolateEulerFieldToSwarmForParticle(
1340 fieldName,
1341 cornerPtr_read_with_ghosts,
1342 &particle,
1343 swarmOut, p, bs);
1344 CHKERRQ(ierr);
1345 }
1346
1347 // (H) Restore all retrieved arrays and vectors */
1348 ierr = DMDAVecRestoreArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
1349
1350 // (I) Restore swarm fields
1351 ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1352 ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1353 ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1354 ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1355 ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
1356 ierr = DMSwarmRestoreField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
1357 LOG_ALLOW(GLOBAL, LOG_INFO, "Finished for field '%s'.\n", fieldName);
1358 PROFILE_FUNCTION_END;
1359 PetscFunctionReturn(0);
1360}
InterpolateEulerFieldToSwarmForParticle
static PetscErrorCode InterpolateEulerFieldToSwarmForParticle(const char *fieldName, void *fieldPtr, Particle *particle, void *swarmOut, PetscInt p, PetscInt blockSize)
Internal helper implementation: InterpolateEulerFieldToSwarmForParticle().
Definition interpolation.c:735
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:183
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:1515
get_log_level
LogLevel get_log_level()
Retrieves the current logging level from the environment variable LOG_LEVEL.
Definition logging.c:84
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:32
UserCtx::lCellFieldAtCorner
Vec lCellFieldAtCorner
Definition variables.h:846
UserCtx::CellFieldAtCorner
Vec CellFieldAtCorner
Definition variables.h:846
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateEulerFieldToSwarm()

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

Dispatches grid-to-particle interpolation to the method selected in the control file.

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

Routes to InterpolateEulerFieldFromCenterToSwarm (direct trilinear, second-order) or InterpolateEulerFieldFromCornerToSwarm (corner-averaged, legacy) based on user->simCtx->interpolationMethod.

Definition at line 1371 of file interpolation.c.

1376{
1377 PetscErrorCode ierr;
1378 PetscFunctionBegin;
1379
1380 if (user->simCtx->interpolationMethod == INTERP_TRILINEAR) {
1381 ierr = InterpolateEulerFieldFromCenterToSwarm(user, fieldLocal_cellCentered,
1382 fieldName, swarmOutFieldName); CHKERRQ(ierr);
1383 } else {
1384 ierr = InterpolateEulerFieldFromCornerToSwarm(user, fieldLocal_cellCentered,
1385 fieldName, swarmOutFieldName); CHKERRQ(ierr);
1386 }
1387
1388 PetscFunctionReturn(0);
1389}
InterpolateEulerFieldFromCornerToSwarm
static PetscErrorCode InterpolateEulerFieldFromCornerToSwarm(UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
Corner-averaged interpolation path (legacy).
Definition interpolation.c:1059
InterpolateEulerFieldFromCenterToSwarm
static PetscErrorCode InterpolateEulerFieldFromCenterToSwarm(UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
Direct cell-center trilinear interpolation (second-order on curvilinear grids).
Definition interpolation.c:853
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
SimCtx::interpolationMethod
InterpolationMethod interpolationMethod
Definition variables.h:744
INTERP_TRILINEAR
@ INTERP_TRILINEAR
Definition variables.h:522
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateAllFieldsToSwarm()

PetscErrorCode InterpolateAllFieldsToSwarm ( UserCtx *  user)

Internal helper implementation: InterpolateAllFieldsToSwarm().

Interpolates all relevant fields from the DMDA to the DMSwarm.

Local to this translation unit.

Definition at line 1397 of file interpolation.c.

1398{
1399 PetscErrorCode ierr;
1400 PetscMPIInt rank;
1401 PetscFunctionBegin;
1402
1403 PROFILE_FUNCTION_BEGIN;
1404
1405 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
1406
1407 /*
1408 1) Interpolate the 'velocity' field (user->Ucat) into DMSwarm's "swarmVelocity".
1409 - The function InterpolateOneFieldOverSwarm() loops over all local particles,
1410 retrieves iCell, jCell, kCell, (a1,a2,a3) from the DMSwarm,
1411 and calls the appropriate trilinear interpolation routine.
1412
1413 - "velocity" => just a label used in logging or debugging.
1414 - "swarmVelocity" => the DMSwarm field name where we store the result
1415 (assume dof=3 if user->Ucat has blockSize=3).
1416 */
1417
1418 LOG_ALLOW(LOCAL, LOG_DEBUG,
1419 " Interpolation of ucat to velocity begins on rank %d.\n",rank);
1420 // Make sure to pass the *LOCAL* Vector to the function below!
1421 ierr = InterpolateEulerFieldToSwarm(user, user->lUcat,
1422 "Ucat",
1423 "velocity"); CHKERRQ(ierr);
1424 ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivity,
1425 "Diffusivity", "Diffusivity"); CHKERRQ(ierr);
1426
1427 ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivityGradient,
1428 "DiffusivityGradient", "DiffusivityGradient"); CHKERRQ(ierr);
1429 /* Add fields as necessary here*/
1430
1431 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1432 LOG_ALLOW(LOCAL, LOG_INFO,
1433 "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
1434
1435 PROFILE_FUNCTION_END;
1436
1437 PetscFunctionReturn(0);
1438}
InterpolateEulerFieldToSwarm
PetscErrorCode InterpolateEulerFieldToSwarm(UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
Dispatches grid-to-particle interpolation to the method selected in the control file.
Definition interpolation.c:1371
UserCtx::lDiffusivityGradient
Vec lDiffusivityGradient
Definition variables.h:841
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:840
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 
)

Internal helper implementation: InterpolateCornerToFaceCenter_Scalar().

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

Local to this translation unit.

Definition at line 2069 of file interpolation.c.

2075{
2076 PetscErrorCode ierr;
2077 DMDALocalInfo info;
2078
2079 PetscFunctionBeginUser;
2080
2081 PROFILE_FUNCTION_BEGIN;
2082
2083 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2084
2085 // Determine owned-cell ranges based on corner-node ownership
2086 PetscInt xs, xm, ys, ym, zs, zm;
2087 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2088 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2089 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2090
2091 // Global exclusive end indices for cells
2092 PetscInt xe = xs + xm;
2093 PetscInt ye = ys + ym;
2094 PetscInt ze = zs + zm;
2095
2096 // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
2097 for (PetscInt k = zs; k < ze; ++k) {
2098 PetscInt k_loc = k - zs;
2099 for (PetscInt j = ys; j < ye; ++j) {
2100 PetscInt j_loc = j - ys;
2101 for (PetscInt i = xs; i <= xe; ++i) {
2102 PetscInt i_loc = i - xs; // 0..xm
2103 // Average the four corners of the Y-Z face at X = i
2104 PetscReal sum = corner_arr[k ][j ][i]
2105 + corner_arr[k+1][j ][i]
2106 + corner_arr[k ][j+1][i]
2107 + corner_arr[k+1][j+1][i];
2108 faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2109 }
2110 }
2111 }
2112
2113 // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
2114 for (PetscInt k = zs; k < ze; ++k) {
2115 PetscInt k_loc = k - zs;
2116 for (PetscInt j = ys; j <= ye; ++j) {
2117 PetscInt j_loc = j - ys; // 0..ym
2118 for (PetscInt i = xs; i < xe; ++i) {
2119 PetscInt i_loc = i - xs;
2120 // Average the four corners of the X-Z face at Y = j
2121 PetscReal sum = corner_arr[k ][j][i ]
2122 + corner_arr[k+1][j][i ]
2123 + corner_arr[k ][j][i+1]
2124 + corner_arr[k+1][j][i+1];
2125 faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2126 }
2127 }
2128 }
2129
2130 // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
2131 for (PetscInt k = zs; k <= ze; ++k) {
2132 PetscInt k_loc = k - zs;
2133 for (PetscInt j = ys; j < ye; ++j) {
2134 PetscInt j_loc = j - ys;
2135 for (PetscInt i = xs; i < xe; ++i) {
2136 PetscInt i_loc = i - xs;
2137 // Average the four corners of the X-Y face at Z = k
2138 PetscReal sum = corner_arr[k][j ][i ]
2139 + corner_arr[k][j ][i+1]
2140 + corner_arr[k][j+1][i ]
2141 + corner_arr[k][j+1][i+1];
2142 faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2143 }
2144 }
2145 }
2146
2147 PROFILE_FUNCTION_END;
2148
2149 PetscFunctionReturn(0);
2150}
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:1883
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 
)

Internal helper implementation: InterpolateCornerToFaceCenter_Vector().

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

Local to this translation unit.

Definition at line 2159 of file interpolation.c.

2165{
2166 PetscErrorCode ierr;
2167 DMDALocalInfo info;
2168 PetscMPIInt rank;
2169
2170 PetscFunctionBeginUser;
2171
2172 PROFILE_FUNCTION_BEGIN;
2173
2174 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2175 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2176 "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
2177
2178 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2179
2180 PetscInt xs, xm, ys, ym, zs, zm;
2181 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2182 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2183 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2184
2185 PetscInt xe = xs + xm;
2186 PetscInt ye = ys + ym;
2187 PetscInt ze = zs + zm;
2188
2189 // X-faces
2190 for (PetscInt k = zs; k < ze; ++k) {
2191 PetscInt k_loc = k - zs;
2192 for (PetscInt j = ys; j < ye; ++j) {
2193 PetscInt j_loc = j - ys;
2194 for (PetscInt i = xs; i <= xe; ++i) {
2195 PetscInt i_loc = i - xs;
2196 Cmpnts sum = {0,0,0};
2197 sum.x = corner_arr[k ][j ][i].x + corner_arr[k+1][j ][i].x
2198 + corner_arr[k ][j+1][i].x + corner_arr[k+1][j+1][i].x;
2199 sum.y = corner_arr[k ][j ][i].y + corner_arr[k+1][j ][i].y
2200 + corner_arr[k ][j+1][i].y + corner_arr[k+1][j+1][i].y;
2201 sum.z = corner_arr[k ][j ][i].z + corner_arr[k+1][j ][i].z
2202 + corner_arr[k ][j+1][i].z + corner_arr[k+1][j+1][i].z;
2203 faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2204 faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2205 faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2206 }
2207 }
2208 }
2209
2210 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2211 "Rank %d x-face Interpolation complete.\n", rank);
2212
2213 // Y-faces
2214 for (PetscInt k = zs; k < ze; ++k) {
2215 PetscInt k_loc = k - zs;
2216 for (PetscInt j = ys; j <= ye; ++j) {
2217 PetscInt j_loc = j - ys;
2218 for (PetscInt i = xs; i < xe; ++i) {
2219 PetscInt i_loc = i - xs;
2220 Cmpnts sum = {0,0,0};
2221 sum.x = corner_arr[k ][j][i ].x + corner_arr[k+1][j][i ].x
2222 + corner_arr[k ][j][i+1].x + corner_arr[k+1][j][i+1].x;
2223 sum.y = corner_arr[k ][j][i ].y + corner_arr[k+1][j][i ].y
2224 + corner_arr[k ][j][i+1].y + corner_arr[k+1][j][i+1].y;
2225 sum.z = corner_arr[k ][j][i ].z + corner_arr[k+1][j][i ].z
2226 + corner_arr[k ][j][i+1].z + corner_arr[k+1][j][i+1].z;
2227 faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2228 faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2229 faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2230 }
2231 }
2232 }
2233
2234 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2235 "Rank %d y-face Interpolation complete.\n", rank);
2236
2237 // Z-faces
2238 for (PetscInt k = zs; k <= ze; ++k) {
2239 PetscInt k_loc = k - zs;
2240 for (PetscInt j = ys; j < ye; ++j) {
2241 PetscInt j_loc = j - ys;
2242 for (PetscInt i = xs; i < xe; ++i) {
2243 PetscInt i_loc = i - xs;
2244 Cmpnts sum = {0,0,0};
2245 sum.x = corner_arr[k][j ][i ].x + corner_arr[k][j ][i+1].x
2246 + corner_arr[k][j+1][i ].x + corner_arr[k][j+1][i+1].x;
2247 sum.y = corner_arr[k][j ][i ].y + corner_arr[k][j ][i+1].y
2248 + corner_arr[k][j+1][i ].y + corner_arr[k][j+1][i+1].y;
2249 sum.z = corner_arr[k][j ][i ].z + corner_arr[k][j ][i+1].z
2250 + corner_arr[k][j+1][i ].z + corner_arr[k][j+1][i+1].z;
2251 faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2252 faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2253 faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2254 }
2255 }
2256 }
2257
2258 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2259 "Rank %d z-face Interpolation complete.\n", rank);
2260
2261 PROFILE_FUNCTION_END;
2262 PetscFunctionReturn(0);
2263}
Here is the call graph for this function:
Generated by doxygen 1.9.8