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

Detailed Description

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

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

Definition in file interpolation.c.

Macro Definition Documentation

◆ NUM_WEIGHTS

#define NUM_WEIGHTS   8

Definition at line 13 of file interpolation.c.

◆ ERROR_MSG_BUFFER_SIZE

#define ERROR_MSG_BUFFER_SIZE   256

Definition at line 16 of file interpolation.c.

◆ __FUNCT [1/6]

#define __FUNCT   "PiecWiseLinearInterpolation_Scalar"

Definition at line 488 of file interpolation.c.

◆ __FUNCT [2/6]

#define __FUNCT   "PiecWiseLinearInterpolation_Vector"

Definition at line 488 of file interpolation.c.

◆ __FUNCT [3/6]

#define __FUNCT   "ComputeTrilinearWeights"

Definition at line 488 of file interpolation.c.

◆ __FUNCT [4/6]

#define __FUNCT   "TrilinearInterpolation_Scalar"

Definition at line 488 of file interpolation.c.

◆ __FUNCT [5/6]

#define __FUNCT   "TrilinearInterpolation_Vector"

Definition at line 488 of file interpolation.c.

◆ __FUNCT [6/6]

#define __FUNCT   "InterpolateEulerFieldToSwarmForParticle"

Definition at line 488 of file interpolation.c.

Function Documentation

◆ InterpolateFieldFromCornerToCenter_Vector()

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

Interpolates a vector field from corner nodes to cell centers.

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

This function estimates the value of a vector field at the center of each computational cell (whose origin node is owned by this rank) by averaging the values from the 8 corner nodes defining that cell.

Parameters
[in]field_arrInput: 3D array view (ghosted) of vector data at nodes, obtained from a Vec associated with user->fda. Accessed via GLOBAL node indices. Ghost values must be up-to-date.
[out]centfield_arrOutput: 3D local array (0-indexed) where interpolated cell center values are stored. Sized [zm_cell_local][ym_cell_local][xm_cell_local] by the caller. centfield_arr[k_loc][j_loc][i_loc] stores the value for the cell corresponding to the (i_loc, j_loc, k_loc)-th owned cell origin.
[in]userUser context containing DMDA information (fda, IM, JM, KM).
Returns
PetscErrorCode 0 on success.

Definition at line 39 of file interpolation.c.

43{
44 PetscErrorCode ierr;
45 PetscMPIInt rank;
46 PetscFunctionBeginUser;
47 PROFILE_FUNCTION_BEGIN;
48
49 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
50 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
51 "Rank %d starting InterpolateFieldFromCornerToCenter_Vector.\n", rank);
52
53 // Get local info based on the NODE-based DMDA (user->fda)
54 // This DM defines the ownership of the nodes from which we interpolate,
55 // and thus defines the cells whose centers we will compute.
56 DMDALocalInfo info_fda;
57 ierr = DMDAGetLocalInfo(user->fda, &info_fda); CHKERRQ(ierr);
58
59 // Determine owned CELL ranges based on owning their origin NODES (from user->fda)
60 PetscInt xs_cell_global_i, xm_cell_local_i; // Global start cell index, local count of owned cells
61 PetscInt ys_cell_global_j, ym_cell_local_j;
62 PetscInt zs_cell_global_k, zm_cell_local_k;
63
64
65 ierr = GetOwnedCellRange(&info_fda, 0, &xs_cell_global_i, &xm_cell_local_i); CHKERRQ(ierr);
66 ierr = GetOwnedCellRange(&info_fda, 1, &ys_cell_global_j, &ym_cell_local_j); CHKERRQ(ierr);
67 ierr = GetOwnedCellRange(&info_fda, 2, &zs_cell_global_k, &zm_cell_local_k); CHKERRQ(ierr);
68
69 // Exclusive end global cell indices
70 PetscInt xe_cell_global_i_excl = xs_cell_global_i + xm_cell_local_i;
71 PetscInt ye_cell_global_j_excl = ys_cell_global_j + ym_cell_local_j;
72 PetscInt ze_cell_global_k_excl = zs_cell_global_k + zm_cell_local_k;
73
74 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d : Processing Owned Cells (global indices) k=%d..%d (count %d), j=%d..%d (count %d), i=%d..%d (count %d)\n",
75 rank,
76 zs_cell_global_k, ze_cell_global_k_excl -1, zm_cell_local_k,
77 ys_cell_global_j, ye_cell_global_j_excl -1, ym_cell_local_j,
78 xs_cell_global_i, xe_cell_global_i_excl -1, xm_cell_local_i);
79
80 // Loop over the GLOBAL indices of the CELLS whose origin nodes are owned by this processor
81 // (according to user->fda partitioning)
82 for (PetscInt k_glob_cell = zs_cell_global_k; k_glob_cell < ze_cell_global_k_excl; k_glob_cell++) {
83 PetscInt k_local = k_glob_cell - zs_cell_global_k; // 0-based local index for centfield_arr
84
85 for (PetscInt j_glob_cell = ys_cell_global_j; j_glob_cell < ye_cell_global_j_excl; j_glob_cell++) {
86 PetscInt j_local = j_glob_cell - ys_cell_global_j; // 0-based local index
87
88 for (PetscInt i_glob_cell = xs_cell_global_i; i_glob_cell < xe_cell_global_i_excl; i_glob_cell++) {
89 PetscInt i_local = i_glob_cell - xs_cell_global_i; // 0-based local index
90
91 Cmpnts sum = {0.0, 0.0, 0.0};
92 // Count is always 8 for a hexahedral cell
93 // PetscInt count = 0; // Not strictly needed if we assume 8 corners
94
95 // Loop over the 8 corner NODES that define cell C(i_glob_cell, j_glob_cell, k_glob_cell)
96 for (PetscInt dk = 0; dk < 2; dk++) {
97 for (PetscInt dj = 0; dj < 2; dj++) {
98 for (PetscInt di = 0; di < 2; di++) {
99 PetscInt ni_glob = i_glob_cell + di; // Global NODE index i or i+1
100 PetscInt nj_glob = j_glob_cell + dj; // Global NODE index j or j+1
101 PetscInt nk_glob = k_glob_cell + dk; // Global NODE index k or k+1
102
103 // Access the input 'field_arr' (node values) using GLOBAL node indices.
104 // Ghosts in field_arr must be valid.
105 // DMDAVecGetArrayRead handles mapping global indices to local memory.
106 sum.x += field_arr[nk_glob][nj_glob][ni_glob].x;
107 sum.y += field_arr[nk_glob][nj_glob][ni_glob].y;
108 sum.z += field_arr[nk_glob][nj_glob][ni_glob].z;
109 // count++;
110 }
111 }
112 }
113
114 // Calculate average value representing the cell center
115 Cmpnts center_value;
116 center_value.x = sum.x / 8.0;
117 center_value.y = sum.y / 8.0;
118 center_value.z = sum.z / 8.0;
119
120 // Store result into the local, 0-indexed output array 'centfield_arr'
121 // Defensive check (should not be strictly necessary if caller allocated centfield_arr
122 // with dimensions xm_cell_local_i, ym_cell_local_j, zm_cell_local_k)
123 if (i_local < 0 || i_local >= xm_cell_local_i ||
124 j_local < 0 || j_local >= ym_cell_local_j ||
125 k_local < 0 || k_local >= zm_cell_local_k) {
126 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local index out of bounds for centfield_arr write!");
127 }
128 centfield_arr[k_local][j_local][i_local] = center_value;
129
130 } // End loop i_glob_cell
131 } // End loop j_glob_cell
132 } // End loop k_glob_cell
133
134 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
135 "Rank %d completed InterpolateFieldFromCornerToCenter_Vector.\n", rank);
136
137 PROFILE_FUNCTION_END;
138
139 PetscFunctionReturn(0);
140}
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:268
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
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:201
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:34
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
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:731
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global, PetscInt *xm_cell_local)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:1663
Cmpnts::x
PetscScalar x
Definition variables.h:101
Cmpnts::z
PetscScalar z
Definition variables.h:101
Cmpnts::y
PetscScalar y
Definition variables.h:101
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
Here is the call graph for this function:

◆ InterpolateFieldFromCornerToCenter_Scalar()

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

Interpolates a scalar field from corner nodes to cell centers.

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

This function estimates the value of a scalar field at the center of each computational cell (whose origin node is owned by this rank) by averaging the values from the 8 corner nodes defining that cell.

Parameters
[in]field_arrInput: 3D array view (ghosted) of scalar data at nodes, obtained from a Vec associated with user->da. Accessed via GLOBAL node indices. Ghost values must be up-to-date.
[out]centfield_arrOutput: 3D local array (0-indexed) where interpolated cell center values are stored. Sized [zm_cell_local][ym_cell_local][xm_cell_local] by the caller. centfield_arr[k_loc][j_loc][i_loc] stores the value for the cell corresponding to the (i_loc, j_loc, k_loc)-th owned cell origin.
[in]userUser context containing DMDA information (da, IM, JM, KM).
Returns
PetscErrorCode 0 on success.

Definition at line 166 of file interpolation.c.

170{
171 PetscErrorCode ierr;
172 PetscMPIInt rank;
173 PetscFunctionBeginUser;
174 PROFILE_FUNCTION_BEGIN;
175 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
176 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
177 "Rank %d starting InterpolateFieldFromCornerToCenter_Scalar.\n", rank);
178
179 // Get local info based on user->da (which is node-based but defines cell origins)
180 DMDALocalInfo info_da_nodes;
181 ierr = DMDAGetLocalInfo(user->da, &info_da_nodes); CHKERRQ(ierr);
182
183 // Determine owned CELL ranges based on owning their origin NODES (from user->da)
184 PetscInt xs_cell_global_i, xm_cell_local_i; // Global start cell index, local count of owned cells
185 PetscInt ys_cell_global_j, ym_cell_local_j;
186 PetscInt zs_cell_global_k, zm_cell_local_k;
187
188 ierr = GetOwnedCellRange(&info_da_nodes, 0, &xs_cell_global_i, &xm_cell_local_i); CHKERRQ(ierr);
189 ierr = GetOwnedCellRange(&info_da_nodes, 1, &ys_cell_global_j, &ym_cell_local_j); CHKERRQ(ierr);
190 ierr = GetOwnedCellRange(&info_da_nodes, 2, &zs_cell_global_k, &zm_cell_local_k); CHKERRQ(ierr);
191
192 // Exclusive end global cell indices
193 PetscInt xe_cell_global_i_excl = xs_cell_global_i + xm_cell_local_i;
194 PetscInt ye_cell_global_j_excl = ys_cell_global_j + ym_cell_local_j;
195 PetscInt ze_cell_global_k_excl = zs_cell_global_k + zm_cell_local_k;
196
197 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d: Processing Owned Cells (global indices) k=%d..%d (count %d), j=%d..%d (count %d), i=%d..%d (count %d)\n",
198 rank,
199 zs_cell_global_k, ze_cell_global_k_excl -1, zm_cell_local_k,
200 ys_cell_global_j, ye_cell_global_j_excl -1, ym_cell_local_j,
201 xs_cell_global_i, xe_cell_global_i_excl -1, xm_cell_local_i);
202
203 // Loop over the GLOBAL indices of the CELLS whose origin nodes are owned by this processor
204 // (according to user->da partitioning)
205 for (PetscInt k_glob_cell = zs_cell_global_k; k_glob_cell < ze_cell_global_k_excl; k_glob_cell++) {
206 PetscInt k_local = k_glob_cell - zs_cell_global_k; // 0-based local index for centfield_arr
207
208 for (PetscInt j_glob_cell = ys_cell_global_j; j_glob_cell < ye_cell_global_j_excl; j_glob_cell++) {
209 PetscInt j_local = j_glob_cell - ys_cell_global_j; // 0-based local index
210
211 for (PetscInt i_glob_cell = xs_cell_global_i; i_glob_cell < xe_cell_global_i_excl; i_glob_cell++) {
212 PetscInt i_local = i_glob_cell - xs_cell_global_i; // 0-based local index
213
214 PetscReal sum = 0.0;
215 // Count is always 8 for a hexahedral cell if all nodes are accessible
216 // PetscInt count = 0;
217
218 // Loop over the 8 NODE indices that define cell C(i_glob_cell, j_glob_cell, k_glob_cell)
219 for (PetscInt dk = 0; dk < 2; dk++) {
220 for (PetscInt dj = 0; dj < 2; dj++) {
221 for (PetscInt di = 0; di < 2; di++) {
222 PetscInt ni_glob = i_glob_cell + di; // Global NODE index for corner
223 PetscInt nj_glob = j_glob_cell + dj; // Global NODE index for corner
224 PetscInt nk_glob = k_glob_cell + dk; // Global NODE index for corner
225
226 // Access the input 'field_arr' (node values) using GLOBAL node indices.
227 sum += field_arr[nk_glob][nj_glob][ni_glob];
228 // count++;
229 }
230 }
231 }
232
233 PetscReal center_value = sum / 8.0;
234
235 // Store the calculated cell center value into the output array 'centfield_arr'
236 // using the LOCAL 0-based cell index.
237 // Defensive check (should not be strictly necessary if caller allocated centfield_arr correctly)
238 if (i_local < 0 || i_local >= xm_cell_local_i ||
239 j_local < 0 || j_local >= ym_cell_local_j ||
240 k_local < 0 || k_local >= zm_cell_local_k) {
241 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local index out of bounds for centfield_arr write in Scalar version!");
242 }
243 centfield_arr[k_local][j_local][i_local] = center_value;
244
245 } // End loop i_glob_cell
246 } // End loop j_glob_cell
247 } // End loop k_glob_cell
248
249 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
250 "Rank %d completed InterpolateFieldFromCornerToCenter_Scalar.\n", rank);
251
252 PROFILE_FUNCTION_END;
253
254 PetscFunctionReturn(0);
255}
Here is the call graph for this function:

◆ InterpolateFieldFromCenterToCorner_Vector_Petsc()

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

Interpolates a vector field from cell centers to corner nodes.

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

Definition at line 270 of file interpolation.c.

274{
275 PetscErrorCode ierr;
276 DMDALocalInfo info;
277 PetscMPIInt rank;
278 PROFILE_FUNCTION_BEGIN;
279 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
280 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
281
282 // Node ownership range (GLOBAL indices)
283 PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
284 PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
285 PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
286
287 PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
288 PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
289 PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
290
291
292 // Global grid dimensions (used for valid cell check)
293 PetscInt IM = info.mx - 1; // Total nodes in i-direction
294 PetscInt JM = info.my - 1; // Total nodes in j-direction
295 PetscInt KM = info.mz - 1; // Total nodes in k-direction
296
297 // Ghosted array starting indices (from the same info struct)
298 PetscInt gxs = info.gxs;
299 PetscInt gys = info.gys;
300 PetscInt gzs = info.gzs;
301
302 PetscInt xs = info.xs;
303 PetscInt ys = info.ys;
304 PetscInt zs = info.zs;
305
306 LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
307 "[Rank %d] Starting InterpolateFieldFromCenterToCorner_Vector_Petsc: Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
308 rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
309
310 // Loop over the GLOBAL indices of the NODES owned by this processor
311 for (PetscInt k = zs_node; k < ze_node; k++) {
312 for (PetscInt j = ys_node; j < ye_node; j++) {
313 for (PetscInt i = xs_node; i < xe_node; i++) {
314 Cmpnts sum = {0.0, 0.0, 0.0};
315 PetscInt count = 0;
316
317 // DEBUG 1 TEST
318 /*
319 if(rank == 1 && i == 24 && j == 12 && k == 49){
320 PetscInt i_cell_A = i - 1;
321 PetscInt i_cell_B = i;
322
323 Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
324 Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
325
326 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",
327 rank, k, j, i,
328 k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
329 k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
330
331 }
332 */
333
334 // Skip processing the unused last node in each dimension.
335 if(i >= IM || j >= JM || k >= KM){
336 continue;
337 }
338 // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
339 // 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
340 for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
341 for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
342 for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
343
344 // These are still GLOBAL cell indices
345 PetscInt global_cell_k = k + dk_offset;
346 PetscInt global_cell_j = j + dj_offset;
347 PetscInt global_cell_i = i + di_offset;
348
349 // Check if this corresponds to a valid GLOBAL cell index
350 if (global_cell_i >= 0 && global_cell_i < nCellsX &&
351 global_cell_j >= 0 && global_cell_j < nCellsY &&
352 global_cell_k >= 0 && global_cell_k < nCellsZ)
353 {
354 Cmpnts cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
355
356 LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f, %.2f, %.2f)\n",
357 rank,global_cell_k,global_cell_j,global_cell_i,cell_val.x, cell_val.y, cell_val.z);
358
359 sum.x += cell_val.x;
360 sum.y += cell_val.y;
361 sum.z += cell_val.z;
362 count++;
363 }
364 }
365 }
366 }
367
368 PetscInt i_global_write = i; // Global index in GLOBAL array.
369 PetscInt j_global_write = j;
370 PetscInt k_global_write = k;
371
372 // We write directly into the array using the global loop indices.
373 if (count > 0) {
374 corner_arr[k_global_write][j_global_write][i_global_write].x = sum.x / (PetscReal)count;
375 corner_arr[k_global_write][j_global_write][i_global_write].y = sum.y / (PetscReal)count;
376 corner_arr[k_global_write][j_global_write][i_global_write].z = sum.z / (PetscReal)count;
377 } else {
378 // This case should ideally not happen for a valid owned node, but as a failsafe:
379 corner_arr[k_global_write][j_global_write][i_global_write] = (Cmpnts){0.0, 0.0, 0.0};
380 }
381
382 // DEBUG 2
383 /*
384 if(rank == 1){
385 if(i == 11 && j == 11 && k == 49){
386 Cmpnts ucat_node = corner_arr[k][j][i];
387 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",
388 rank, k, j, i,
389 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
390 }
391 }
392
393 if(rank == 0 && i == 24 && j == 12 && k == 0){
394 Cmpnts ucat_node = corner_arr[k][j][i];
395 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",
396 rank, k, j, i,
397 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
398 }
399 */
400 // 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);
401 }
402 }
403 }
404 PROFILE_FUNCTION_END;
405 return 0;
406}
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:349
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35

◆ InterpolateFieldFromCenterToCorner_Scalar_Petsc()

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

Interpolates a scalar field from cell centers to corner nodes.

This version is adapted to write directly into a ghosted local array obtained from DMDAVecGetArray(), which allows using GLOBAL indices for writing to the OWNED portion of the array.

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

Definition at line 424 of file interpolation.c.

428{
429 PetscErrorCode ierr;
430 DMDALocalInfo info;
431 PROFILE_FUNCTION_BEGIN;
432 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
433
434 // Node ownership range (GLOBAL indices)
435 PetscInt xs_node = info.xs, xe_node = info.xs + info.xm;
436 PetscInt ys_node = info.ys, ye_node = info.ys + info.ym;
437 PetscInt zs_node = info.zs, ze_node = info.zs + info.zm;
438
439 // Global grid dimensions (number of nodes)
440 PetscInt IM = info.mx-1;
441 PetscInt JM = info.my-1;
442 PetscInt KM = info.mz-1;
443
444 // Loop over the GLOBAL indices of the NODES owned by this processor
445 for (PetscInt k = zs_node; k < ze_node; k++) {
446 for (PetscInt j = ys_node; j < ye_node; j++) {
447 for (PetscInt i = xs_node; i < xe_node; i++) {
448 PetscReal sum = 0.0;
449 PetscInt count = 0;
450
451 // Loop over the 8 potential cells surrounding node N(k,j,i)
452 for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
453 for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
454 for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
455 // This is the global index of the cell-centered data point
456 PetscInt cell_idx_k = k + dk_offset;
457 PetscInt cell_idx_j = j + dj_offset;
458 PetscInt cell_idx_i = i + di_offset;
459
460 // Check if this corresponds to a valid GLOBAL cell index.
461 // Cell indices run from 0 to NumberOfNodes-2.
462 if (cell_idx_i >= 0 && cell_idx_i < IM - 1 &&
463 cell_idx_j >= 0 && cell_idx_j < JM - 1 &&
464 cell_idx_k >= 0 && cell_idx_k < KM - 1)
465 {
466 sum += centfield_arr[cell_idx_k][cell_idx_j][cell_idx_i];
467 count++;
468 }
469 }
470 }
471 }
472
473 // Calculate average and write to the output array using GLOBAL indices
474 if (count > 0) {
475 corner_arr[k][j][i] = sum / (PetscReal)count;
476 } else {
477 corner_arr[k][j][i] = 0.0; // Failsafe
478 }
479 }
480 }
481 }
482 PROFILE_FUNCTION_END;
483 return 0;
484}

◆ PieceWiseLinearInterpolation_Scalar()

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

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

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

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

Definition at line 502 of file interpolation.c.

509{
510 PetscFunctionBegin;
511 *val = fieldScal[kCell][jCell][iCell];
512
513 // Optional logging
514 LOG_ALLOW(LOCAL, LOG_VERBOSE,
515 "Field '%s' at (i=%d, j=%d, k=%d) => val=%.6f\n",
516 fieldName, iCell, jCell, kCell, *val);
517
518 PetscFunctionReturn(0);
519}

◆ PieceWiseLinearInterpolation_Vector()

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

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

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

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

Definition at line 539 of file interpolation.c.

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

◆ ComputeTrilinearWeights()

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

Computes the trilinear interpolation weights from the interpolation coefficients.

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

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

Definition at line 579 of file interpolation.c.

579 {
580 LOG_ALLOW(GLOBAL, LOG_VERBOSE, "Computing weights for a1=%f, a2=%f, a3=%f.\n", a1, a2, a3);
581
582 // Ensure a1, a2, a3 are within [0,1]
583 a1 = PetscMax(0.0, PetscMin(1.0, a1));
584 a2 = PetscMax(0.0, PetscMin(1.0, a2));
585 a3 = PetscMax(0.0, PetscMin(1.0, a3));
586
587 const PetscReal oa1 = 1.0 - a1;
588 const PetscReal oa2 = 1.0 - a2;
589 const PetscReal oa3 = 1.0 - a3;
590
591 w[0] = oa1 * oa2 * oa3; /* cornerOffsets[0] => (0,0,0) */
592 w[1] = a1 * oa2 * oa3; /* cornerOffsets[1] => (1,0,0) */
593 w[2] = oa1 * a2 * oa3; /* cornerOffsets[2] => (0,1,0) */
594 w[3] = a1 * a2 * oa3; /* cornerOffsets[3] => (1,1,0) */
595 w[4] = oa1 * oa2 * a3; /* cornerOffsets[4] => (0,0,1) */
596 w[5] = a1 * oa2 * a3; /* cornerOffsets[5] => (1,0,1) */
597 w[6] = oa1 * a2 * a3; /* cornerOffsets[6] => (0,1,1) */
598 w[7] = a1 * a2 * a3; /* cornerOffsets[7] => (1,1,1) */
599
600 // Log the computed weights for debugging
601 LOG_ALLOW(LOCAL,LOG_VERBOSE, "Weights computed - "
602 "w0=%f, w1=%f, w2=%f, w3=%f, w4=%f, w5=%f, w6=%f, w7=%f. \n",
603 w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7]);
604}
Here is the caller graph for this function:

◆ TrilinearInterpolation_Scalar()

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

Computes the trilinear Interpolateted scalar at a given point.

Computes the trilinear interpolated scalar at a given point.

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

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

Definition at line 622 of file interpolation.c.

632{
633 PetscFunctionBegin; // PETSc macro for error/stack tracing
634
635 // Compute the 8 corner weights
636 PetscReal wcorner[8];
637 ComputeTrilinearWeights(a1, a2, a3, wcorner);
638
639 // Offsets for cell corners
640 PetscInt i1 = i + 1;
641 PetscInt j1 = j + 1;
642 PetscInt k1 = k + 1;
643
644 // Initialize the output scalar
645 PetscReal sum = 0.0;
646
647 // Corner 0 => (i, j, k)
648 sum += wcorner[0] * fieldScal[k ][j ][i ];
649 // Corner 1 => (i+1, j, k)
650 sum += wcorner[1] * fieldScal[k ][j ][i1];
651 // Corner 2 => (i, j+1, k)
652 sum += wcorner[2] * fieldScal[k ][j1][i ];
653 // Corner 3 => (i+1, j+1, k)
654 sum += wcorner[3] * fieldScal[k ][j1][i1];
655 // Corner 4 => (i, j, k+1)
656 sum += wcorner[4] * fieldScal[k1][j ][i ];
657 // Corner 5 => (i+1, j, k+1)
658 sum += wcorner[5] * fieldScal[k1][j ][i1];
659 // Corner 6 => (i, j+1, k+1)
660 sum += wcorner[6] * fieldScal[k1][j1][i ];
661 // Corner 7 => (i+1, j+1, k+1)
662 sum += wcorner[7] * fieldScal[k1][j1][i1];
663
664 *val = sum;
665
666 // Logging (optional)
667 LOG_ALLOW(LOCAL, LOG_VERBOSE,
668 "Field '%s' at (i=%d, j=%d, k=%d), "
669 "a1=%.6f, a2=%.6f, a3=%.6f -> val=%.6f.\n",
670 fieldName, i, j, k, a1, a2, a3, *val);
671
672 // LOG_ALLOW_SYNC(GLOBAL, LOG_INFO,
673 // "TrilinearInterpolation_Scalar: Completed interpolation for field '%s' across local cells.\n",
674 // fieldName);
675
676 PetscFunctionReturn(0);
677}
ComputeTrilinearWeights
static void ComputeTrilinearWeights(PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w)
Computes the trilinear interpolation weights from the interpolation coefficients.
Definition interpolation.c:579
Here is the call graph for this function:

◆ TrilinearInterpolation_Vector()

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

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

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

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

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

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

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

Definition at line 702 of file interpolation.c.

712{
713 PetscFunctionBegin; // PETSc macro for error/stack tracing
714
715 // Compute the 8 corner weights
716 PetscErrorCode ierr;
717 PetscReal wcorner[8];
718 PetscMPIInt rank;
719
720 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
721
722 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Computing Trilinear Weights.\n",rank);
723 ComputeTrilinearWeights(a1, a2, a3, wcorner);
724
725 LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Trilinear weights computed for local cell %d,%d,%d.\n",rank,i,j,k);
726
727 // For partial interpolation, we'll keep track of how many corners are valid
728 // and how much sum of weights is used. Then we do a final normalization.
729 PetscReal sumW = 0.0;
730 Cmpnts accum = {0.0, 0.0, 0.0};
731
732 // The eight corner indices, with their weights:
733 // corners: (i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), etc.
734 // We store them in an array to iterate cleanly.
735 const PetscInt cornerOffsets[8][3] = {
736 {0, 0, 0},
737 {1, 0, 0},
738 {0, 1, 0},
739 {1, 1, 0},
740 {0, 0, 1},
741 {1, 0, 1},
742 {0, 1, 1},
743 {1, 1, 1}
744 };
745
746 // Weighted partial sum
747 for (PetscInt c = 0; c < 8; c++) {
748 const PetscInt di = cornerOffsets[c][0];
749 const PetscInt dj = cornerOffsets[c][1];
750 const PetscInt dk = cornerOffsets[c][2];
751 PetscInt iC = i + di;
752 PetscInt jC = j + dj;
753 PetscInt kC = k + dk;
754
755 /*
756 // skip if out of domain
757 // (Assuming you know global domain is [0..mx), [0..my), [0..mz).)
758 if (iC < 0 || iC >= (PetscInt)userGlobalMx ||
759 jC < 0 || jC >= (PetscInt)userGlobalMy ||
760 kC < 0 || kC >= (PetscInt)userGlobalMz)
761 {
762 // skip this corner
763 continue;
764 }
765
766 */
767
768 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);
769
770 // Otherwise, accumulate
771 accum.x += wcorner[c] * fieldVec[kC][jC][iC].x;
772 accum.y += wcorner[c] * fieldVec[kC][jC][iC].y;
773 accum.z += wcorner[c] * fieldVec[kC][jC][iC].z;
774 sumW += wcorner[c];
775 }
776
777 // If sumW=0 => out-of-range or zero weighting => set (0,0,0)
778 if (sumW > 1.0e-14) {
779 vec->x = accum.x / sumW;
780 vec->y = accum.y / sumW;
781 vec->z = accum.z / sumW;
782 } else {
783 vec->x = 0.0; vec->y = 0.0; vec->z = 0.0;
784 }
785
786 PetscFunctionReturn(0);
787}
Here is the call graph for this function:

◆ InterpolateEulerFieldToSwarmForParticle()

static PetscErrorCode InterpolateEulerFieldToSwarmForParticle ( const char *  fieldName,
void *  fieldPtr,
PetscInt  iCell,
PetscInt  jCell,
PetscInt  kCell,
PetscReal  a1,
PetscReal  a2,
PetscReal  a3,
void *  swarmOut,
PetscInt  p,
PetscInt  blockSize 
)
inlinestatic

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

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

Parameters
[in]fieldNameA string identifying the field (e.g. "velocity", "temperature", etc.).
[in]fieldPtrA Pointer to the local DMDA array for the field:
  • (PetscReal ***) if scalar (blockSize = 1)
  • (Cmpnts ***) if vector (blockSize = 3)
[in]iCell,jCell,kCellThe cell indices for this particle (already clamped if needed).
[in]a1,a2,a3Interpolation coefficients in [0..1]. If using PiecewiseLinearInterpolation, these are currently ignored.
[out]swarmOutA Pointer to the DMSwarm output array:
  • (PetscReal*) if scalar dof=1
  • (PetscReal*) if vector dof=3 (storing x,y,z in consecutive reals)
  • or (Cmpnts*) if you store the result as a Cmpnts struct
[in]pThe local particle index (used to compute the correct offset in swarmOut).
[in]blockSizeThe number of degrees of freedom (1 => scalar, 3 => vector).

This routine demonstrates a switch between:

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

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

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

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

Definition at line 822 of file interpolation.c.

834{
835 PetscErrorCode ierr;
836 PetscFunctionBegin;
837
838 PROFILE_FUNCTION_BEGIN;
839
840 // Optional logging at start
841 LOG_ALLOW(LOCAL, LOG_VERBOSE,
842 "field='%s', blockSize=%d, "
843 "cell IDs=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
844 fieldName, blockSize, iCell, jCell, kCell, a1, a2, a3);
845
846 /*
847 If blockSize=1, we PetscInterpret the fieldPtr as a 3D array of PetscReal (scalar).
848 If blockSize=3, we PetscInterpret the fieldPtr as a 3D array of Cmpnts (vector).
849 */
850 if (blockSize == 1) {
851 /* Scalar field: Cast fieldPtr to (PetscReal ***). */
852 PetscReal ***fieldScal = (PetscReal ***) fieldPtr;
853 PetscReal val;
854
855 // Currently using trilinear.
856 ierr = TrilinearInterpolation(fieldName, fieldScal,
857 iCell, jCell, kCell,
858 a1, a2, a3,
859 &val);
860 CHKERRQ(ierr);
861
862 // Alternative (commented) call to PiecewiseLinearInterpolation (zeroth order) :
863 // PetscErrorCode ierr = PiecewiseLinearInterpolation(fieldName,
864 // fieldScal,
865 // iCell, jCell, kCell,
866 // &val);
867 // CHKERRQ(ierr);
868
869 // Write the scalar result to the swarm output at index [p].
870 ((PetscReal*)swarmOut)[p] = val;
871
872 LOG_ALLOW(LOCAL, LOG_VERBOSE,
873 "field='%s', result=%.6f "
874 "stored at swarmOut index p=%d.\n", fieldName, val, (PetscInt)p);
875 }
876 else if (blockSize == 3) {
877 /* Vector field: Cast fieldPtr to (Cmpnts ***). */
878 Cmpnts ***fieldVec = (Cmpnts ***) fieldPtr;
879 Cmpnts vec;
880
881
882
883
884 // Piecewise interpolation (zeroth order).
885 // PetscErrorCode ierr = PieceWiseLinearInterpolation(fieldName,
886 // fieldVec,
887 // iCell, jCell, kCell,
888 // &vec);
889 // CHKERRQ(ierr);
890
891 // Alternative (commented) call to trilinear:
892 ierr = TrilinearInterpolation(fieldName, fieldVec,
893 iCell, jCell, kCell,
894 a1, a2, a3,
895 &vec);
896 CHKERRQ(ierr);
897
898 // If swarmOut is an array of 3 reals per particle:
899 ((PetscReal*)swarmOut)[3*p + 0] = vec.x;
900 ((PetscReal*)swarmOut)[3*p + 1] = vec.y;
901 ((PetscReal*)swarmOut)[3*p + 2] = vec.z;
902
903 LOG_ALLOW(LOCAL, LOG_VERBOSE,
904 "field='%s', result=(%.6f,%.6f,%.6f) "
905 "stored at swarmOut[3p..3p+2], p=%d.\n",
906 fieldName, vec.x, vec.y, vec.z, (PetscInt)p);
907
908 /*
909 If you store the vector result as a Cmpnts in the swarm, do instead:
910 ((Cmpnts*)swarmOut)[p] = vec;
911 but ensure your DMSwarm field is sized for a Cmpnts struct.
912 */
913 }
914 else {
915 /* If blockSize isn't 1 or 3, we raise an error. */
916 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
917 "InterpolateEulerFieldToSwarmForParticle: only blockSize=1 or 3 supported, got %d.",
918 (PetscInt)blockSize);
919 }
920
921 PROFILE_FUNCTION_END;
922 PetscFunctionReturn(0);
923}
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:130
Here is the caller graph for this function:

◆ InterpolateEulerFieldToSwarm()

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

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

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

Workflow:

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

Definition at line 958 of file interpolation.c.

963{
964 PetscErrorCode ierr;
965 DM fda = user->fda;
966 DM swarm = user->swarm;
967 PetscInt bs;
968 DMDALocalInfo info;
969 PetscMPIInt rank;
970
971 // Generic pointers to the raw data arrays
972 void *cellCenterPtr_read;
973 void *cornerPtr_read_with_ghosts;
974
975 // Swarm-related pointers
976 PetscInt *cellIDs = NULL;
977 PetscReal *weights = NULL;
978 PetscInt *pids = NULL; // For logging particle IDs in warnings
979 void *swarmOut = NULL;
980 PetscInt nLocal;
981
982 PetscFunctionBegin;
983 PROFILE_FUNCTION_BEGIN;
984 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
985 ierr = DMDAGetLocalInfo(fda, &info); CHKERRQ(ierr);
986 ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
987 if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
988
989 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting for field '%s'... \n", fieldName);
990
991 /* STAGE 1: Center-to-Corner Calculation with Communication */
992
993 // (A) Create the corner vectors.
994 if(user->CellFieldAtCorner == NULL){
995 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new global corner vector for field '%s'.\n", rank, fieldName);
996 ierr = DMCreateGlobalVector(fda, &user->CellFieldAtCorner); CHKERRQ(ierr);
997 }else{
998 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing global corner vector for field '%s'.\n", rank, fieldName);
999 }
1000
1001 ierr = VecSet(user->CellFieldAtCorner,0.0); CHKERRQ(ierr);
1002 Vec cornerGlobal = user->CellFieldAtCorner;
1003
1004 if(user->lCellFieldAtCorner == NULL){
1005 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new local corner vector for field '%s'.\n", rank, fieldName);
1006 ierr = DMCreateLocalVector(fda, &user->lCellFieldAtCorner); CHKERRQ(ierr);
1007 }else{
1008 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing local corner vector for field '%s'.\n", rank, fieldName);
1009 }
1010
1011 ierr = VecSet(user->lCellFieldAtCorner,0.0); CHKERRQ(ierr);
1012 Vec cornerLocal = user->lCellFieldAtCorner;
1013
1014 // (B) Get a read-only array from the input cell-centered vector
1015 ierr = DMDAVecGetArrayRead(fda, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
1016
1017 PetscInt size;
1018 ierr = VecGetSize(cornerGlobal, &size); CHKERRQ(ierr);
1019 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1020
1021 size = 0;
1022
1023 ierr = VecGetSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
1024 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1025
1026 size = 0;
1027
1028 ierr = VecGetLocalSize(cornerGlobal, &size); CHKERRQ(ierr);
1029 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector LOCAL size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1030
1031 size = 0;
1032
1033 ierr = VecGetLocalSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
1034 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector LOCAL size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1035
1036 size = 0;
1037
1038 ierr = VecGetSize(cornerLocal, &size); CHKERRQ(ierr);
1039 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
1040
1041 PetscInt xs,ys,zs,gxs,gys,gzs;
1042 ierr = DMDAGetCorners(fda,&xs,&ys,&zs,NULL,NULL,NULL); CHKERRQ(ierr);
1043 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);
1044
1045 ierr = DMDAGetGhostCorners(fda,&gxs,&gys,&gzs,NULL,NULL,NULL); CHKERRQ(ierr);
1046 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);
1047
1048 // DEBUG: Inspect Ucat ghost cells before interpolation
1049 /*
1050 if (bs == 3) {
1051 Cmpnts ***ucat_array = (Cmpnts***)cellCenterPtr_read;
1052 DMDALocalInfo info_debug;
1053 ierr = DMDAGetLocalInfo(fda, &info_debug); CHKERRQ(ierr);
1054
1055 // Only print on rank 0 to avoid clutter
1056 if (rank == 0) {
1057 PetscPrintf(PETSC_COMM_SELF, "\nDEBUG: Inspecting Ucat Ghost Cells...\n");
1058 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n");
1059
1060 // --- Check the MINIMUM-SIDE (-Xi) ---
1061 // The first physical cell is at index info.xs (local index).
1062 // The first ghost cell is at info.xs - 1.
1063 PetscInt i_first_phys = info_debug.xs;
1064 PetscInt i_first_ghost = info_debug.xs - 1;
1065 PetscInt j_mid = info_debug.ys + info_debug.ym / 2;
1066 PetscInt k_mid = info_debug.zs + info_debug.zm / 2;
1067
1068 PetscPrintf(PETSC_COMM_SELF, "MIN-SIDE (-Xi):\n");
1069 PetscPrintf(PETSC_COMM_SELF, " Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1070 k_mid, j_mid, i_first_ghost,
1071 ucat_array[k_mid][j_mid][i_first_ghost].x,
1072 ucat_array[k_mid][j_mid][i_first_ghost].y,
1073 ucat_array[k_mid][j_mid][i_first_ghost].z);
1074 PetscPrintf(PETSC_COMM_SELF, " First Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1075 k_mid, j_mid, i_first_phys,
1076 ucat_array[k_mid][j_mid][i_first_phys].x,
1077 ucat_array[k_mid][j_mid][i_first_phys].y,
1078 ucat_array[k_mid][j_mid][i_first_phys].z);
1079
1080
1081 // --- Check the MAXIMUM-SIDE (+Xi) ---
1082 // The last physical cell is at index info.xs + info.xm - 1.
1083 // The first ghost cell on the max side is at info.xs + info.xm.
1084 PetscInt i_last_phys = info_debug.xs + info_debug.xm - 1;
1085 PetscInt i_last_ghost = info_debug.xs + info_debug.xm;
1086
1087 PetscPrintf(PETSC_COMM_SELF, "MAX-SIDE (+Xi):\n");
1088 PetscPrintf(PETSC_COMM_SELF, " Last Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1089 k_mid, j_mid, i_last_phys,
1090 ucat_array[k_mid][j_mid][i_last_phys].x,
1091 ucat_array[k_mid][j_mid][i_last_phys].y,
1092 ucat_array[k_mid][j_mid][i_last_phys].z);
1093 PetscPrintf(PETSC_COMM_SELF, " Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
1094 k_mid, j_mid, i_last_ghost,
1095 ucat_array[k_mid][j_mid][i_last_ghost].x,
1096 ucat_array[k_mid][j_mid][i_last_ghost].y,
1097 ucat_array[k_mid][j_mid][i_last_ghost].z);
1098 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n\n");
1099 }
1100 }
1101 */
1102 // DEBUG
1103
1104 // (C) Perform the center-to-corner interpolation directly into the global corner vector
1105 void *cornerPtr_write = NULL;
1106 ierr = DMDAVecGetArray(fda, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
1107
1108 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Starting center-to-corner interpolation for '%s'.\n", rank, fieldName);
1109
1110 // SINGLE, CLEAN CALL SITE: The macro handles the runtime dispatch based on 'bs'.
1111 ierr = InterpolateFieldFromCenterToCorner(bs, cellCenterPtr_read, cornerPtr_write, user); CHKERRQ(ierr);
1112
1113 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Finished center-to-corner interpolation for '%s'.\n", rank, fieldName);
1114 ierr = DMDAVecRestoreArray(fda, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
1115
1116 ierr = DMDAVecRestoreArrayRead(fda, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
1117
1118 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1119
1120 //////////// DEBUG
1121 /*
1122 // Log a synchronized header message from all ranks.
1123 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerGlobal BEFORE ghost exchange...\n");
1124
1125 // VecView prints to a PETSc viewer. PETSC_VIEWER_STDOUT_WORLD is a global, synchronized viewer.
1126 // We only need to call it if logging is active for this function.
1127 if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
1128 ierr = VecView(cornerGlobal, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1129 }
1130
1131 // Log a synchronized footer message.
1132 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerGlobal.\n");
1133 */
1134 ////////// DEBUG
1135
1136 // ierr = PetscBarrier((PetscObject)cornerGlobal); CHKERRQ(ierr);
1137
1138 // (D) CRITICAL STEP: Communicate the newly computed corner data to fill ghost regions
1139 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Beginning ghost exchange for corner data...\n", rank);
1140
1141 ierr = DMGlobalToLocalBegin(fda, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
1142 ierr = DMGlobalToLocalEnd(fda, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
1143
1144 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Ghost exchange for corner data complete.\n", rank);
1145
1146 LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE, "getting array from FDA on all ranks.\n");
1147 if (is_function_allowed(__func__) && (int)(LOG_VERBOSE) <= (int)get_log_level()) {
1148 // DEBUG: Inspect cornerLocal after ghost exchange
1149
1150 ierr = LOG_FIELD_ANATOMY(user,"CornerField", "After Corner Velocity Interpolated"); CHKERRQ(ierr);
1151
1152 // DEBUG
1153 //ierr = DMView(user->fda, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1154 }
1155
1156 /////// DEBUG
1157 /*
1158 // Log a synchronized header message from all ranks.
1159 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerLocal AFTER ghost exchange...\n");
1160
1161 // Here, we want each rank to print its own local vector.
1162 // PETSC_VIEWER_STDOUT_SELF is the correct tool for this. The LOG_ALLOW_SYNC
1163 // wrapper will ensure the output from each rank is printed sequentially.
1164 if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
1165 PetscMPIInt rank_d;
1166 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_d); CHKERRQ(ierr);
1167
1168 // Use a synchronized print to create a clean header for each rank's output
1169 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n--- cornerLocal (Rank %d) ---\n", rank_d);
1170 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
1171
1172 // Print the local vector's contents for this rank
1173 ierr = VecView(cornerLocal, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
1174
1175 // Use a synchronized print to create a clean footer
1176 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--- End cornerLocal (Rank %d) ---\n", rank_d);
1177 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
1178 }
1179
1180 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerLocal.\n");
1181 */
1182 ///// DEBUG
1183
1184 /* STAGE 2: Particle Interpolation using Ghosted Corner Data */
1185
1186 // (E) Get the local, GHOSTED array of corner data for the final interpolation step
1187 ierr = DMDAVecGetArrayRead(fda, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
1188
1189 // (F) Retrieve swarm fields for the particle loop
1190 ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
1191 ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1192 ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1193 ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1194 ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1195
1196 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);
1197
1198 // (G) Loop over each local particle
1199 for (PetscInt p = 0; p < nLocal; p++) {
1200 PetscInt iCell_global = cellIDs[3 * p + 0];
1201 PetscInt jCell_global = cellIDs[3 * p + 1];
1202 PetscInt kCell_global = cellIDs[3 * p + 2];
1203
1204
1205 LOG_ALLOW(LOCAL, LOG_VERBOSE,
1206 "[Rank %d] Particle PID %lld: global cell=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
1207 rank, (long long)pids[p], iCell_global, jCell_global, kCell_global,
1208 weights[3 * p + 0], weights[3 * p + 1], weights[3 * p + 2]);
1209 // Safety check: Ensure the entire 8-node stencil is within the valid memory
1210 // region (owned + ghosts) of the local array.
1211
1212 if (iCell_global < info.gxs || iCell_global >= info.gxs + info.gxm - 1 ||
1213 jCell_global < info.gxs || jCell_global >= info.gys + info.gym - 1 ||
1214 kCell_global < info.gxs || kCell_global >= info.gzs + info.gzm - 1)
1215 {
1216 LOG_ALLOW(LOCAL, LOG_WARNING,
1217 "[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",
1218 rank, (long long)pids[p], iCell_global, jCell_global, kCell_global, fieldName);
1219
1220 if (bs == 3) {
1221 ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
1222 ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
1223 ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
1224 } else {
1225 ((PetscReal*)swarmOut)[p] = 0.0;
1226 }
1227 continue;
1228 }
1229
1230 PetscReal a1 = weights[3 * p + 0];
1231 PetscReal a2 = weights[3 * p + 1];
1232 PetscReal a3 = weights[3 * p + 2];
1233
1234 ierr = InterpolateEulerFieldToSwarmForParticle(
1235 fieldName,
1236 cornerPtr_read_with_ghosts,
1237 iCell_global, jCell_global, kCell_global,
1238 a1, a2, a3,
1239 swarmOut, p, bs);
1240 CHKERRQ(ierr);
1241 }
1242
1243 /* (H) Restore all retrieved arrays and vectors */
1244 ierr = DMDAVecRestoreArrayRead(fda, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
1245 ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
1246 ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
1247 ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
1248 ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
1249
1250 LOG_ALLOW(GLOBAL, LOG_INFO, "Finished for field '%s'.\n", fieldName);
1251 PROFILE_FUNCTION_END;
1252 PetscFunctionReturn(0);
1253}
InterpolateEulerFieldToSwarmForParticle
static PetscErrorCode InterpolateEulerFieldToSwarmForParticle(const char *fieldName, void *fieldPtr, PetscInt iCell, PetscInt jCell, PetscInt kCell, PetscReal a1, PetscReal a2, PetscReal a3, void *swarmOut, PetscInt p, PetscInt blockSize)
Interpolates a single field (scalar or vector) for one particle, storing the result in a swarm array.
Definition interpolation.c:822
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:64
is_function_allowed
PetscBool is_function_allowed(const char *functionName)
Checks if a given function is in the allow-list.
Definition logging.c:157
LOG_FIELD_ANATOMY
PetscErrorCode LOG_FIELD_ANATOMY(UserCtx *user, const char *field_name, const char *stage_name)
Logs the anatomy of a specified field at key boundary locations, respecting the solver's specific gri...
Definition logging.c:1435
get_log_level
LogLevel get_log_level()
Retrieves the current logging level from the environment variable LOG_LEVEL.
Definition logging.c:39
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
UserCtx::lCellFieldAtCorner
Vec lCellFieldAtCorner
Definition variables.h:693
UserCtx::CellFieldAtCorner
Vec CellFieldAtCorner
Definition variables.h:693
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateAllFieldsToSwarm()

PetscErrorCode InterpolateAllFieldsToSwarm ( UserCtx *  user)

Interpolates all relevant fields from the DMDA to the DMSwarm.

Currently, it Interpolatetes:

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

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

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

Definition at line 1276 of file interpolation.c.

1277{
1278 PetscErrorCode ierr;
1279 PetscMPIInt rank;
1280 PetscFunctionBegin;
1281
1282 PROFILE_FUNCTION_BEGIN;
1283
1284 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
1285
1286 /*
1287 1) Interpolate the 'velocity' field (user->Ucat) into DMSwarm's "swarmVelocity".
1288 - The function InterpolateOneFieldOverSwarm() loops over all local particles,
1289 retrieves iCell, jCell, kCell, (a1,a2,a3) from the DMSwarm,
1290 and calls the appropriate trilinear interpolation routine.
1291
1292 - "velocity" => just a label used in logging or debugging.
1293 - "swarmVelocity" => the DMSwarm field name where we store the result
1294 (assume dof=3 if user->Ucat has blockSize=3).
1295 */
1296
1297 LOG_ALLOW(LOCAL, LOG_DEBUG,
1298 " Interpolation of ucat to velocity begins on rank %d.\n",rank);
1299 // Make sure to pass the *GLOBAL* Vector to the function below!
1300 ierr = InterpolateEulerFieldToSwarm(user, user->lUcat,
1301 "Ucat",
1302 "velocity"); CHKERRQ(ierr);
1303 /*
1304 2) (OPTIONAL) If you have more fields, you can Interpolatete them similarly.
1305
1306 For example, if user->Tcat is a scalar Vec for "temperature" and the DMSwarm
1307 has a field "Temperature" (with dof=1), you could do:
1308
1309 ierr = InterpolateOneFieldOverSwarm(user, user->Tcat,
1310 "temperature",
1311 "swarmTemperature"); CHKERRQ(ierr);
1312
1313 For pressure:
1314 ierr = InterpolateOneFieldOverSwarm(user, user->Pcat,
1315 "pressure",
1316 "swarmPressure"); CHKERRQ(ierr);
1317
1318 ...and so forth.
1319 */
1320
1321 /*
1322 3) Optionally, synchronize or log that all fields are done
1323 */
1324 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
1325 LOG_ALLOW(LOCAL, LOG_INFO,
1326 "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
1327
1328 PROFILE_FUNCTION_END;
1329
1330 PetscFunctionReturn(0);
1331}
InterpolateEulerFieldToSwarm
PetscErrorCode InterpolateEulerFieldToSwarm(UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
Interpolates a cell-centered field (scalar or vector) onto DMSwarm particles, using a robust,...
Definition interpolation.c:958
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InterpolateCornerToFaceCenter_Scalar()

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

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

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

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

Definition at line 2004 of file interpolation.c.

2010{
2011 PetscErrorCode ierr;
2012 DMDALocalInfo info;
2013
2014 PetscFunctionBeginUser;
2015
2016 PROFILE_FUNCTION_BEGIN;
2017
2018 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2019
2020 // Determine owned-cell ranges based on corner-node ownership
2021 PetscInt xs, xm, ys, ym, zs, zm;
2022 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2023 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2024 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2025
2026 // Global exclusive end indices for cells
2027 PetscInt xe = xs + xm;
2028 PetscInt ye = ys + ym;
2029 PetscInt ze = zs + zm;
2030
2031 // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
2032 for (PetscInt k = zs; k < ze; ++k) {
2033 PetscInt k_loc = k - zs;
2034 for (PetscInt j = ys; j < ye; ++j) {
2035 PetscInt j_loc = j - ys;
2036 for (PetscInt i = xs; i <= xe; ++i) {
2037 PetscInt i_loc = i - xs; // 0..xm
2038 // Average the four corners of the Y-Z face at X = i
2039 PetscReal sum = corner_arr[k ][j ][i]
2040 + corner_arr[k+1][j ][i]
2041 + corner_arr[k ][j+1][i]
2042 + corner_arr[k+1][j+1][i];
2043 faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2044 }
2045 }
2046 }
2047
2048 // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
2049 for (PetscInt k = zs; k < ze; ++k) {
2050 PetscInt k_loc = k - zs;
2051 for (PetscInt j = ys; j <= ye; ++j) {
2052 PetscInt j_loc = j - ys; // 0..ym
2053 for (PetscInt i = xs; i < xe; ++i) {
2054 PetscInt i_loc = i - xs;
2055 // Average the four corners of the X-Z face at Y = j
2056 PetscReal sum = corner_arr[k ][j][i ]
2057 + corner_arr[k+1][j][i ]
2058 + corner_arr[k ][j][i+1]
2059 + corner_arr[k+1][j][i+1];
2060 faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2061 }
2062 }
2063 }
2064
2065 // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
2066 for (PetscInt k = zs; k <= ze; ++k) {
2067 PetscInt k_loc = k - zs;
2068 for (PetscInt j = ys; j < ye; ++j) {
2069 PetscInt j_loc = j - ys;
2070 for (PetscInt i = xs; i < xe; ++i) {
2071 PetscInt i_loc = i - xs;
2072 // Average the four corners of the X-Y face at Z = k
2073 PetscReal sum = corner_arr[k][j ][i ]
2074 + corner_arr[k][j ][i+1]
2075 + corner_arr[k][j+1][i ]
2076 + corner_arr[k][j+1][i+1];
2077 faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
2078 }
2079 }
2080 }
2081
2082 PROFILE_FUNCTION_END;
2083
2084 PetscFunctionReturn(0);
2085}
Here is the call graph for this function:

◆ InterpolateCornerToFaceCenter_Vector()

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

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

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

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

Definition at line 2104 of file interpolation.c.

2110{
2111 PetscErrorCode ierr;
2112 DMDALocalInfo info;
2113 PetscMPIInt rank;
2114
2115 PetscFunctionBeginUser;
2116
2117 PROFILE_FUNCTION_BEGIN;
2118
2119 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2120 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2121 "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
2122
2123 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2124
2125 PetscInt xs, xm, ys, ym, zs, zm;
2126 ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
2127 ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
2128 ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
2129
2130 PetscInt xe = xs + xm;
2131 PetscInt ye = ys + ym;
2132 PetscInt ze = zs + zm;
2133
2134 // X-faces
2135 for (PetscInt k = zs; k < ze; ++k) {
2136 PetscInt k_loc = k - zs;
2137 for (PetscInt j = ys; j < ye; ++j) {
2138 PetscInt j_loc = j - ys;
2139 for (PetscInt i = xs; i <= xe; ++i) {
2140 PetscInt i_loc = i - xs;
2141 Cmpnts sum = {0,0,0};
2142 sum.x = corner_arr[k ][j ][i].x + corner_arr[k+1][j ][i].x
2143 + corner_arr[k ][j+1][i].x + corner_arr[k+1][j+1][i].x;
2144 sum.y = corner_arr[k ][j ][i].y + corner_arr[k+1][j ][i].y
2145 + corner_arr[k ][j+1][i].y + corner_arr[k+1][j+1][i].y;
2146 sum.z = corner_arr[k ][j ][i].z + corner_arr[k+1][j ][i].z
2147 + corner_arr[k ][j+1][i].z + corner_arr[k+1][j+1][i].z;
2148 faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2149 faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2150 faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2151 }
2152 }
2153 }
2154
2155 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2156 "Rank %d x-face Interpolation complete.\n", rank);
2157
2158 // Y-faces
2159 for (PetscInt k = zs; k < ze; ++k) {
2160 PetscInt k_loc = k - zs;
2161 for (PetscInt j = ys; j <= ye; ++j) {
2162 PetscInt j_loc = j - ys;
2163 for (PetscInt i = xs; i < xe; ++i) {
2164 PetscInt i_loc = i - xs;
2165 Cmpnts sum = {0,0,0};
2166 sum.x = corner_arr[k ][j][i ].x + corner_arr[k+1][j][i ].x
2167 + corner_arr[k ][j][i+1].x + corner_arr[k+1][j][i+1].x;
2168 sum.y = corner_arr[k ][j][i ].y + corner_arr[k+1][j][i ].y
2169 + corner_arr[k ][j][i+1].y + corner_arr[k+1][j][i+1].y;
2170 sum.z = corner_arr[k ][j][i ].z + corner_arr[k+1][j][i ].z
2171 + corner_arr[k ][j][i+1].z + corner_arr[k+1][j][i+1].z;
2172 faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2173 faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2174 faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2175 }
2176 }
2177 }
2178
2179 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2180 "Rank %d y-face Interpolation complete.\n", rank);
2181
2182 // Z-faces
2183 for (PetscInt k = zs; k <= ze; ++k) {
2184 PetscInt k_loc = k - zs;
2185 for (PetscInt j = ys; j < ye; ++j) {
2186 PetscInt j_loc = j - ys;
2187 for (PetscInt i = xs; i < xe; ++i) {
2188 PetscInt i_loc = i - xs;
2189 Cmpnts sum = {0,0,0};
2190 sum.x = corner_arr[k][j ][i ].x + corner_arr[k][j ][i+1].x
2191 + corner_arr[k][j+1][i ].x + corner_arr[k][j+1][i+1].x;
2192 sum.y = corner_arr[k][j ][i ].y + corner_arr[k][j ][i+1].y
2193 + corner_arr[k][j+1][i ].y + corner_arr[k][j+1][i+1].y;
2194 sum.z = corner_arr[k][j ][i ].z + corner_arr[k][j ][i+1].z
2195 + corner_arr[k][j+1][i ].z + corner_arr[k][j+1][i+1].z;
2196 faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
2197 faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
2198 faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
2199 }
2200 }
2201 }
2202
2203 LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
2204 "Rank %d z-face Interpolation complete.\n", rank);
2205
2206 PROFILE_FUNCTION_END;
2207 PetscFunctionReturn(0);
2208}
Here is the call graph for this function: