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__   "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)
 Compute the eight trilinear interpolation weights for a particle's local coordinates.
 
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)
 Interpolate one Eulerian field to a single located swarm particle.
 
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, FieldId source_field_id, ParticleFieldId target_field_id)
 Dispatches grid-to-particle interpolation to the method selected in the control file.
 
PetscErrorCode InterpolateAllFieldsToSwarm (UserCtx *user)
 Internal helper implementation: InterpolateAllFieldsToSwarm().
 
PetscErrorCode AccumulateParticleField (DM swarm, ParticleFieldId particle_field_id, 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, ParticleFieldId particle_field_id, DM targetDM, PetscInt expected_dof, Vec eulerFieldAverageVec)
 Accumulate one particle field onto the Eulerian grid using the selected scatter stencil.
 
PetscErrorCode ScatterParticleFieldToEulerField (UserCtx *user, ParticleFieldId particle_field_id, 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:103
Cmpnts::z
PetscScalar z
Definition variables.h:103
Cmpnts::y
PetscScalar y
Definition variables.h:103
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: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:974
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102

◆ 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:253
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:859
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:335
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:850

◆ 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

Compute the eight trilinear interpolation weights for a particle's local coordinates.

Definition at line 515 of file interpolation.c.

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

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

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

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

Interpolate one Eulerian field to a single located swarm particle.

Definition at line 733 of file interpolation.c.

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

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

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

◆ InterpolateEulerFieldToSwarm()

PetscErrorCode InterpolateEulerFieldToSwarm ( UserCtx *  user,
FieldId  source_field_id,
ParticleFieldId  target_field_id 
)

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

1357{
1358 PetscErrorCode ierr;
1359 FieldView source_view;
1360 const ParticleFieldDescriptor *target_descriptor = NULL;
1361
1362 PetscFunctionBegin;
1363
1364 ierr = FieldGetView(user, source_field_id, &source_view); CHKERRQ(ierr);
1365 ierr = ParticleFieldGetDescriptor(target_field_id, &target_descriptor); CHKERRQ(ierr);
1366 PetscCheck(source_view.descriptor->layout == FIELD_LAYOUT_CELL_CENTERED,
1367 PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
1368 "Eulerian-to-particle interpolation requires a cell-centered source; '%s' uses layout %s.",
1369 source_view.descriptor->canonical_name,
1370 FieldLayoutName(source_view.descriptor->layout));
1371 PetscCheck(target_descriptor->data_type == PETSC_REAL,
1372 PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
1373 "Eulerian interpolation destination '%s' must use PETSC_REAL storage.",
1374 target_descriptor->canonical_name);
1375 PetscCheck(source_view.descriptor->dof == target_descriptor->components,
1376 PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
1377 "Eulerian source '%s' has %d components but particle destination '%s' has %d.",
1378 source_view.descriptor->canonical_name, source_view.descriptor->dof,
1379 target_descriptor->canonical_name, target_descriptor->components);
1380
1381 if (user->simCtx->interpolationMethod == INTERP_TRILINEAR) {
1382 ierr = InterpolateEulerFieldFromCenterToSwarm(user, source_view.local_vec,
1383 source_view.descriptor->canonical_name,
1384 target_descriptor->canonical_name); CHKERRQ(ierr);
1385 } else {
1386 ierr = InterpolateEulerFieldFromCornerToSwarm(user, source_view.local_vec,
1387 source_view.descriptor->canonical_name,
1388 target_descriptor->canonical_name); CHKERRQ(ierr);
1389 }
1390
1391 PetscFunctionReturn(0);
1392}
FieldDescriptor::layout
FieldLayout layout
Definition field_catalog.h:126
FieldDescriptor::dof
PetscInt dof
Definition field_catalog.h:124
FIELD_LAYOUT_CELL_CENTERED
@ FIELD_LAYOUT_CELL_CENTERED
Definition field_catalog.h:73
FieldLayoutName
const char * FieldLayoutName(FieldLayout layout)
Return a stable printable label for a field layout.
Definition field_catalog.c:240
InterpolateEulerFieldFromCornerToSwarm
static PetscErrorCode InterpolateEulerFieldFromCornerToSwarm(UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
Corner-averaged interpolation path (legacy).
Definition interpolation.c:1057
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:851
ParticleFieldDescriptor::components
PetscInt components
Definition particle_field_catalog.h:59
ParticleFieldGetDescriptor
PetscErrorCode ParticleFieldGetDescriptor(ParticleFieldId field_id, const ParticleFieldDescriptor **descriptor)
Return immutable metadata for a valid particle field ID.
Definition particle_field_catalog.c:59
ParticleFieldDescriptor::data_type
PetscDataType data_type
Definition particle_field_catalog.h:60
ParticleFieldDescriptor::canonical_name
const char * canonical_name
Definition particle_field_catalog.h:56
ParticleFieldDescriptor
Immutable metadata for one persistent particle field.
Definition particle_field_catalog.h:54
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
SimCtx::interpolationMethod
InterpolationMethod interpolationMethod
Definition variables.h:832
INTERP_TRILINEAR
@ INTERP_TRILINEAR
Definition variables.h:565
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 1400 of file interpolation.c.

1401{
1402 PetscErrorCode ierr;
1403 PetscMPIInt rank;
1404 PetscFunctionBegin;
1405
1406 PROFILE_FUNCTION_BEGIN;
1407
1408 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
1409
1410 LOG_ALLOW(LOCAL, LOG_DEBUG,
1411 " Interpolation of ucat to velocity begins on rank %d.\n",rank);
1412 ierr = InterpolateEulerFieldToSwarm(user, FIELD_ID_UCAT,
1413 PARTICLE_FIELD_ID_VELOCITY); CHKERRQ(ierr);
1414 ierr = InterpolateEulerFieldToSwarm(user, FIELD_ID_DIFFUSIVITY,
1415 PARTICLE_FIELD_ID_DIFFUSIVITY); CHKERRQ(ierr);
1416 ierr = InterpolateEulerFieldToSwarm(user, FIELD_ID_DIFFUSIVITY_GRADIENT,
1417 PARTICLE_FIELD_ID_DIFFUSIVITY_GRADIENT); CHKERRQ(ierr);
1418 /* Add fields as necessary here*/
1419
1420 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1421 LOG_ALLOW(LOCAL, LOG_INFO,
1422 "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
1423
1424 PROFILE_FUNCTION_END;
1425
1426 PetscFunctionReturn(0);
1427}
FIELD_ID_UCAT
@ FIELD_ID_UCAT
Definition field_catalog.h:26
FIELD_ID_DIFFUSIVITY_GRADIENT
@ FIELD_ID_DIFFUSIVITY_GRADIENT
Definition field_catalog.h:34
FIELD_ID_DIFFUSIVITY
@ FIELD_ID_DIFFUSIVITY
Definition field_catalog.h:33
InterpolateEulerFieldToSwarm
PetscErrorCode InterpolateEulerFieldToSwarm(UserCtx *user, FieldId source_field_id, ParticleFieldId target_field_id)
Dispatches grid-to-particle interpolation to the method selected in the control file.
Definition interpolation.c:1353
PARTICLE_FIELD_ID_DIFFUSIVITY_GRADIENT
@ PARTICLE_FIELD_ID_DIFFUSIVITY_GRADIENT
Definition particle_field_catalog.h:30
PARTICLE_FIELD_ID_DIFFUSIVITY
@ PARTICLE_FIELD_ID_DIFFUSIVITY
Definition particle_field_catalog.h:29
PARTICLE_FIELD_ID_VELOCITY
@ PARTICLE_FIELD_ID_VELOCITY
Definition particle_field_catalog.h:26
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 1953 of file interpolation.c.

1959{
1960 PetscErrorCode ierr;
1961 DMDALocalInfo info;
1962
1963 PetscFunctionBeginUser;
1964
1965 PROFILE_FUNCTION_BEGIN;
1966
1967 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
1968
1969 // Determine owned-cell ranges based on corner-node ownership
1970 PetscInt xs, xm, ys, ym, zs, zm;
1971 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
1972 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
1973 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
1974
1975 // Global exclusive end indices for cells
1976 PetscInt xe = xs + xm;
1977 PetscInt ye = ys + ym;
1978 PetscInt ze = zs + zm;
1979
1980 // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
1981 for (PetscInt k = zs; k < ze; ++k) {
1982 PetscInt k_loc = k - zs;
1983 for (PetscInt j = ys; j < ye; ++j) {
1984 PetscInt j_loc = j - ys;
1985 for (PetscInt i = xs; i <= xe; ++i) {
1986 PetscInt i_loc = i - xs; // 0..xm
1987 // Average the four corners of the Y-Z face at X = i
1988 PetscReal sum = corner_arr[k ][j ][i]
1989 + corner_arr[k+1][j ][i]
1990 + corner_arr[k ][j+1][i]
1991 + corner_arr[k+1][j+1][i];
1992 faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
1993 }
1994 }
1995 }
1996
1997 // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
1998 for (PetscInt k = zs; k < ze; ++k) {
1999 PetscInt k_loc = k - zs;
2000 for (PetscInt j = ys; j <= ye; ++j) {
2001 PetscInt j_loc = j - ys; // 0..ym
2002 for (PetscInt i = xs; i < xe; ++i) {
2003 PetscInt i_loc = i - xs;
2004 // Average the four corners of the X-Z face at Y = j
2005 PetscReal sum = corner_arr[k ][j][i ]
2006 + corner_arr[k+1][j][i ]
2007 + corner_arr[k ][j][i+1]
2008 + corner_arr[k+1][j][i+1];
2009 faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2010 }
2011 }
2012 }
2013
2014 // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
2015 for (PetscInt k = zs; k <= ze; ++k) {
2016 PetscInt k_loc = k - zs;
2017 for (PetscInt j = ys; j < ye; ++j) {
2018 PetscInt j_loc = j - ys;
2019 for (PetscInt i = xs; i < xe; ++i) {
2020 PetscInt i_loc = i - xs;
2021 // Average the four corners of the X-Y face at Z = k
2022 PetscReal sum = corner_arr[k][j ][i ]
2023 + corner_arr[k][j ][i+1]
2024 + corner_arr[k][j+1][i ]
2025 + corner_arr[k][j+1][i+1];
2026 faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2027 }
2028 }
2029 }
2030
2031 PROFILE_FUNCTION_END;
2032
2033 PetscFunctionReturn(0);
2034}
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:2285
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 2043 of file interpolation.c.

2049{
2050 PetscErrorCode ierr;
2051 DMDALocalInfo info;
2052 PetscMPIInt rank;
2053
2054 PetscFunctionBeginUser;
2055
2056 PROFILE_FUNCTION_BEGIN;
2057
2058 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2059 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2060 "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
2061
2062 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2063
2064 PetscInt xs, xm, ys, ym, zs, zm;
2065 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2066 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2067 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2068
2069 PetscInt xe = xs + xm;
2070 PetscInt ye = ys + ym;
2071 PetscInt ze = zs + zm;
2072
2073 // X-faces
2074 for (PetscInt k = zs; k < ze; ++k) {
2075 PetscInt k_loc = k - zs;
2076 for (PetscInt j = ys; j < ye; ++j) {
2077 PetscInt j_loc = j - ys;
2078 for (PetscInt i = xs; i <= xe; ++i) {
2079 PetscInt i_loc = i - xs;
2080 Cmpnts sum = {0,0,0};
2081 sum.x = corner_arr[k ][j ][i].x + corner_arr[k+1][j ][i].x
2082 + corner_arr[k ][j+1][i].x + corner_arr[k+1][j+1][i].x;
2083 sum.y = corner_arr[k ][j ][i].y + corner_arr[k+1][j ][i].y
2084 + corner_arr[k ][j+1][i].y + corner_arr[k+1][j+1][i].y;
2085 sum.z = corner_arr[k ][j ][i].z + corner_arr[k+1][j ][i].z
2086 + corner_arr[k ][j+1][i].z + corner_arr[k+1][j+1][i].z;
2087 faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2088 faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2089 faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2090 }
2091 }
2092 }
2093
2094 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2095 "Rank %d x-face Interpolation complete.\n", rank);
2096
2097 // Y-faces
2098 for (PetscInt k = zs; k < ze; ++k) {
2099 PetscInt k_loc = k - zs;
2100 for (PetscInt j = ys; j <= ye; ++j) {
2101 PetscInt j_loc = j - ys;
2102 for (PetscInt i = xs; i < xe; ++i) {
2103 PetscInt i_loc = i - xs;
2104 Cmpnts sum = {0,0,0};
2105 sum.x = corner_arr[k ][j][i ].x + corner_arr[k+1][j][i ].x
2106 + corner_arr[k ][j][i+1].x + corner_arr[k+1][j][i+1].x;
2107 sum.y = corner_arr[k ][j][i ].y + corner_arr[k+1][j][i ].y
2108 + corner_arr[k ][j][i+1].y + corner_arr[k+1][j][i+1].y;
2109 sum.z = corner_arr[k ][j][i ].z + corner_arr[k+1][j][i ].z
2110 + corner_arr[k ][j][i+1].z + corner_arr[k+1][j][i+1].z;
2111 faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2112 faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2113 faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2114 }
2115 }
2116 }
2117
2118 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2119 "Rank %d y-face Interpolation complete.\n", rank);
2120
2121 // Z-faces
2122 for (PetscInt k = zs; k <= ze; ++k) {
2123 PetscInt k_loc = k - zs;
2124 for (PetscInt j = ys; j < ye; ++j) {
2125 PetscInt j_loc = j - ys;
2126 for (PetscInt i = xs; i < xe; ++i) {
2127 PetscInt i_loc = i - xs;
2128 Cmpnts sum = {0,0,0};
2129 sum.x = corner_arr[k][j ][i ].x + corner_arr[k][j ][i+1].x
2130 + corner_arr[k][j+1][i ].x + corner_arr[k][j+1][i+1].x;
2131 sum.y = corner_arr[k][j ][i ].y + corner_arr[k][j ][i+1].y
2132 + corner_arr[k][j+1][i ].y + corner_arr[k][j+1][i+1].y;
2133 sum.z = corner_arr[k][j ][i ].z + corner_arr[k][j ][i+1].z
2134 + corner_arr[k][j+1][i ].z + corner_arr[k][j+1][i+1].z;
2135 faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2136 faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2137 faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2138 }
2139 }
2140 }
2141
2142 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2143 "Rank %d z-face Interpolation complete.\n", rank);
2144
2145 PROFILE_FUNCTION_END;
2146 PetscFunctionReturn(0);
2147}
Here is the call graph for this function:
Generated by doxygen 1.9.8