interpolation.h File Reference

#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscdmcomposite.h>
#include <petscerror.h>
#include <petscsys.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "setup.h"

Include dependency graph for interpolation.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	NUM_WEIGHTS   8
#define	InterpolateFieldFromCornerToCenter(field, centfield, user)
	Generic macro to call the appropriate interpolation function based on the field type.
#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 size variable.
#define	InterpolateCornerToFaceCenter(corner_arr, faceX_arr, faceY_arr, faceZ_arr, user_ctx)
	A type-generic macro that interpolates a field from corner nodes to all face centers.
#define	PieceWiseLinearInterpolation(fieldName, fieldPtr, i, j, k, outPtr)
	Macro that calls either the scalar or vector piecewise interpolation function based on the type of the fieldPtr parameter (3D array).
#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 the fieldPtr parameter (3D array).

Functions
PetscErrorCode	PieceWiseLinearInterpolation_Scalar (const char *fieldName, PetscReal ***fieldScal, PetscInt iCell, PetscInt jCell, PetscInt kCell, PetscReal *val)
	Returns the scalar value from the input cell index without blending.
PetscErrorCode	PieceWiseLinearInterpolation_Vector (const char *fieldName, Cmpnts ***fieldVec, PetscInt iCell, PetscInt jCell, PetscInt kCell, Cmpnts *vec)
	Returns the vector value from the input cell index without blending.
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 interpolated 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 interpolated vector (e.g., velocity) at a given point.
PetscErrorCode	InterpolateEulerFieldToSwarm (UserCtx *user, FieldId source_field_id, ParticleFieldId target_field_id)
	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	InterpolateParticleVelocities (UserCtx *user)
	Interpolates particle velocities using trilinear interpolation.
PetscErrorCode	InterpolateFieldFromCornerToCenter_Scalar (PetscReal ***field_arr, PetscReal ***centfield_arr, UserCtx *user)
	Safely interpolate a scalar field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.
PetscErrorCode	InterpolateFieldFromCornerToCenter_Vector (Cmpnts ***field_arr, Cmpnts ***centfield_arr, UserCtx *user)
	Safely interpolate a vector field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.
PetscErrorCode	TestCornerToCenterInterpolation (UserCtx *user)
	Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Vector (Cmpnts ***centfield_arr, Cmpnts ***corner_arr, UserCtx *user)
	Interpolates a vector field from cell centers to corner nodes.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Scalar (PetscReal ***centfield_arr, PetscReal ***corner_arr, UserCtx *user)
	Interpolates a scalar field from cell centers to corner nodes.
PetscErrorCode	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.
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)
	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.

Macro Definition Documentation

NUM_WEIGHTS

#define NUM_WEIGHTS   8

Definition at line 24 of file interpolation.h.

InterpolateFieldFromCornerToCenter

#define InterpolateFieldFromCornerToCenter	(		field,
			centfield,
			user
	)

Value:

  ( (void)sizeof(char[1 - 2*!!(!__builtin_types_compatible_p(typeof(field), typeof(centfield)))]),\
    _Generic((field),                                                                       \
      PetscReal ***: InterpolateFieldFromCornerToCenter_Scalar,                             \
      Cmpnts ***:    InterpolateFieldFromCornerToCenter_Vector                                \
    )(field, centfield, user) )

Generic macro to call the appropriate interpolation function based on the field type.

This macro will select either the scalar or vector interpolation function based on the type of the 'field' pointer. It also performs a compile-time check that 'centfield' is of the same type as 'field'. If the types are not the same, a compile-time error is produced.

Usage: InterpolateFieldFromCornerToCenter(field, centfield, user);

Definition at line 36 of file interpolation.h.

                   : InterpolateFieldFromCornerToCenter_Scalar,                             \
      Cmpnts ***:    InterpolateFieldFromCornerToCenter_Vector                                \
    )(field, centfield, user) )

InterpolateFieldFromCenterToCorner

#define InterpolateFieldFromCenterToCorner	(		blockSize,
			centfield_ptr,
			corner_ptr,
			user_ctx
	)

Value:

    ( (blockSize) == 1 ? \
        InterpolateFieldFromCenterToCorner_Scalar((PetscReal***)(centfield_ptr), (PetscReal***)(corner_ptr), (user_ctx)) : \
        InterpolateFieldFromCenterToCorner_Vector((Cmpnts***)(centfield_ptr), (Cmpnts***)(corner_ptr), (user_ctx)) \
    )

Macro to dispatch to the correct scalar or vector center-to-corner function based on a runtime block size variable.

This macro uses a ternary operator to inspect the runtime value of 'blockSize' and select the appropriate implementation. It expects the input and output pointers to be void* and handles casting them to the correct, strongly-typed pointers.

Parameters

blockSize	The runtime block size (1 for scalar, 3 for vector).
centfield_ptr	The input void* pointer to the 3D cell-centered data array.
corner_ptr	The output void* pointer to the 3D corner data array.
user_ctx	The UserCtx structure.

Definition at line 56 of file interpolation.h.

                                                                                                                         : \
        InterpolateFieldFromCenterToCorner_Vector((Cmpnts***)(centfield_ptr), (Cmpnts***)(corner_ptr), (user_ctx)) \
    )

InterpolateCornerToFaceCenter

#define InterpolateCornerToFaceCenter	(		corner_arr,
			faceX_arr,
			faceY_arr,
			faceZ_arr,
			user_ctx
	)

Value:

    _Generic((corner_arr),                                                                  \
        PetscReal***: InterpolateCornerToFaceCenter_Scalar,                                 \
        Cmpnts***:    InterpolateCornerToFaceCenter_Vector                                  \
    )(corner_arr, faceX_arr, faceY_arr, faceZ_arr, user_ctx)

A type-generic macro that interpolates a field from corner nodes to all face centers.

This macro uses C11's _Generic feature to dispatch to the appropriate underlying function based on the type of the input corner_arr.

• If corner_arr is of type PetscReal***, it calls InterpolateCornerToFaceCenter_Scalar.
• If corner_arr is of type Cmpnts***, it calls InterpolateCornerToFaceCenter_Vector.

Parameters

corner_arr	Ghosted node-centered array (global indexing). The type of this argument determines which function is called.
faceX_arr	Local array for X-faces.
faceY_arr	Local array for Y-faces.
faceZ_arr	Local array for Z-faces.
user_ctx	User context containing DMDA 'fda'.

Definition at line 79 of file interpolation.h.

                    : InterpolateCornerToFaceCenter_Scalar,                                 \
        Cmpnts***:    InterpolateCornerToFaceCenter_Vector                                  \
    )(corner_arr, faceX_arr, faceY_arr, faceZ_arr, user_ctx)

PieceWiseLinearInterpolation

#define PieceWiseLinearInterpolation	(		fieldName,
			fieldPtr,
			i,
			j,
			k,
			outPtr
	)

Value:

  _Generic((fieldPtr),                                                     \
    PetscReal ***: PieceWiseLinearInterpolation_Scalar,                    \
    Cmpnts      ***: PieceWiseLinearInterpolation_Vector                   \
  )(fieldName, fieldPtr, i, j, k, outPtr)

Macro that calls either the scalar or vector piecewise interpolation function based on the type of the fieldPtr parameter (3D array).

Usage example:

// For scalar: PetscReal ***fieldScal; PetscReal outVal; PieceWiseLinearInterpolation(fieldName, fieldScal, i, j, k, &outVal);

// For vector: Cmpnts ***fieldVec; Cmpnts vec; PieceWiseLinearInterpolation(fieldName, fieldVec, i, j, k, &vec);

Definition at line 101 of file interpolation.h.

                 : PieceWiseLinearInterpolation_Scalar,                    \
    Cmpnts      ***: PieceWiseLinearInterpolation_Vector                   \
  )(fieldName, fieldPtr, i, j, k, outPtr)

TrilinearInterpolation

#define TrilinearInterpolation	(		fieldName,
			fieldPtr,
			i,
			j,
			k,
			a1,
			a2,
			a3,
			outPtr
	)

Value:

  _Generic((fieldPtr),                                                                     \
    PetscReal ***: TrilinearInterpolation_Scalar,                                          \
    Cmpnts      ***: TrilinearInterpolation_Vector                                         \
  )(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 the fieldPtr parameter (3D array).

Usage example:

PetscReal result; Cmpnts vec;

// For scalars: TrilinearInterpolation(fieldName, fieldScal, i, j, k, a1, a2, a3, &result);

// For vectors: TrilinearInterpolation(fieldName, fieldVec, i, j, k, a1, a2, a3, &vec);

Definition at line 166 of file interpolation.h.

                 : TrilinearInterpolation_Scalar,                                          \
    Cmpnts      ***: TrilinearInterpolation_Vector                                         \
  )(fieldName, fieldPtr, i, j, k, a1, a2, a3, outPtr)

Function Documentation

PieceWiseLinearInterpolation_Scalar()

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

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

This is a first-order piecewise-constant lookup helper used in workflows that intentionally disable trilinear blending.

Parameters

[in]	fieldName	Field label used for diagnostics.
[in]	fieldScal	Scalar field array indexed as [k][j][i].
[in]	iCell	Cell i-index.
[in]	jCell	Cell j-index.
[in]	kCell	Cell k-index.
[out]	val	Output scalar value.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 462 of file interpolation.c.

{
  PetscFunctionBegin;
  *val = fieldScal[kCell][jCell][iCell];
  
  // Optional logging
  LOG_ALLOW(LOCAL, LOG_VERBOSE,
      "Field '%s' at (i=%d, j=%d, k=%d) => val=%.6f\n",
      fieldName, iCell, jCell, kCell, *val);
 
  PetscFunctionReturn(0);
}

PieceWiseLinearInterpolation_Vector()

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

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

This is a first-order piecewise-constant lookup helper used in workflows that intentionally disable trilinear blending.

Parameters

[in]	fieldName	Field label used for diagnostics.
[in]	fieldVec	Vector field array indexed as [k][j][i].
[in]	iCell	Cell i-index.
[in]	jCell	Cell j-index.
[in]	kCell	Cell k-index.
[out]	vec	Output vector value.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 487 of file interpolation.c.

{
  PetscFunctionBegin;
  vec->x = fieldVec[kCell][jCell][iCell].x;
  vec->y = fieldVec[kCell][jCell][iCell].y;
  vec->z = fieldVec[kCell][jCell][iCell].z;
 
  // Optional logging
  LOG_ALLOW(LOCAL, LOG_VERBOSE,
      "Field '%s' at (i=%d, j=%d, k=%d) => (x=%.6f, y=%.6f, z=%.6f)\n",
      fieldName, iCell, jCell, kCell, vec->x, vec->y, vec->z);
 
  PetscFunctionReturn(0);
}

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 interpolated scalar at a given point.

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

Parameters

fieldName	A
fieldScal	3D
i	Integral
j	Integral
k	Integral
a1	Normalized
a2	Normalized
a3	Normalized
val	Pointer

Returns

PetscErrorCode 0 on success.

Computes the trilinear interpolated scalar at a given point.

Local to this translation unit.

Definition at line 576 of file interpolation.c.

{
    PetscFunctionBegin; // PETSc macro for error/stack tracing
 
    // Compute the 8 corner weights
    PetscReal wcorner[8];
    ComputeTrilinearWeights(a1, a2, a3, wcorner);
 
    // Offsets for cell corners
    PetscInt i1 = i + 1;
    PetscInt j1 = j + 1;
    PetscInt k1 = k + 1;
 
    // Initialize the output scalar
    PetscReal sum = 0.0;
 
    // Corner 0 => (i, j, k)
    sum += wcorner[0] * fieldScal[k ][j ][i ];
    // Corner 1 => (i+1, j, k)
    sum += wcorner[1] * fieldScal[k ][j ][i1];
    // Corner 2 => (i, j+1, k)
    sum += wcorner[2] * fieldScal[k ][j1][i ];
    // Corner 3 => (i+1, j+1, k)
    sum += wcorner[3] * fieldScal[k ][j1][i1];
    // Corner 4 => (i, j, k+1)
    sum += wcorner[4] * fieldScal[k1][j ][i ];
    // Corner 5 => (i+1, j, k+1)
    sum += wcorner[5] * fieldScal[k1][j ][i1];
    // Corner 6 => (i, j+1, k+1)
    sum += wcorner[6] * fieldScal[k1][j1][i ];
    // Corner 7 => (i+1, j+1, k+1)
    sum += wcorner[7] * fieldScal[k1][j1][i1];
 
    *val = sum;
 
    // Logging (optional)
    LOG_ALLOW(LOCAL, LOG_VERBOSE,
        "Field '%s' at (i=%d, j=%d, k=%d), "
        "a1=%.6f, a2=%.6f, a3=%.6f -> val=%.6f.\n",
        fieldName, i, j, k, a1, a2, a3, *val);
 
    // LOG_ALLOW_SYNC(GLOBAL, LOG_INFO,
    //    "TrilinearInterpolation_Scalar: Completed interpolation for field '%s' across local cells.\n",
      //    fieldName);
 
    PetscFunctionReturn(0);
}

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
	)

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

each cell of type Cmpnts. This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same interface.

Parameters

fieldName	A
fieldVec	3D
i	Integral
j	Integral
k	Integral
a1	Normalized
a2	Normalized
a3	Normalized
vec	Pointer

Returns

PetscErrorCode 0 on success.

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.

{
  PetscFunctionBegin; // PETSc macro for error/stack tracing
 
  // Compute the 8 corner weights
  PetscErrorCode ierr;
  PetscReal wcorner[8];
  PetscMPIInt rank;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  
  LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Computing Trilinear Weights.\n",rank);
  ComputeTrilinearWeights(a1, a2, a3, wcorner);
  
  LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d]  Trilinear weights computed for local cell %d,%d,%d.\n",rank,i,j,k);
 
  // For partial interpolation, we'll keep track of how many corners are valid
  // and how much sum of weights is used. Then we do a final normalization.
  PetscReal sumW = 0.0;
  Cmpnts    accum = {0.0, 0.0, 0.0};
 
  // The eight corner indices, with their weights:
  // corners: (i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), etc.
  // We store them in an array to iterate cleanly.
  const PetscInt cornerOffsets[8][3] = {
    {0, 0, 0},
    {1, 0, 0},
    {0, 1, 0},
    {1, 1, 0},
    {0, 0, 1},
    {1, 0, 1},
    {0, 1, 1},
    {1, 1, 1}
  };
 
  // Weighted partial sum
  for (PetscInt c = 0; c < 8; c++) {
    const PetscInt di = cornerOffsets[c][0];
    const PetscInt dj = cornerOffsets[c][1];
    const PetscInt dk = cornerOffsets[c][2];
    PetscInt iC = i + di;
    PetscInt jC = j + dj;
    PetscInt kC = k + dk;
 
    /*
    // skip if out of domain
    // (Assuming you know global domain is [0..mx), [0..my), [0..mz).)
    if (iC < 0 || iC >= (PetscInt)userGlobalMx ||
        jC < 0 || jC >= (PetscInt)userGlobalMy ||
        kC < 0 || kC >= (PetscInt)userGlobalMz)
    {
      // skip this corner
      continue;
    }
 
    */
 
    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);
    
    // Otherwise, accumulate
    accum.x += wcorner[c] * fieldVec[kC][jC][iC].x;
    accum.y += wcorner[c] * fieldVec[kC][jC][iC].y;
    accum.z += wcorner[c] * fieldVec[kC][jC][iC].z;
    sumW += wcorner[c];
  }
 
  // If sumW=0 => out-of-range or zero weighting => set (0,0,0)
  if (sumW > 1.0e-14) {
    vec->x = accum.x / sumW;
    vec->y = accum.y / sumW;
    vec->z = accum.z / sumW;
  } else {
    vec->x = 0.0;  vec->y = 0.0;  vec->z = 0.0;
  }
 
  PetscFunctionReturn(0);
}

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
	)

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]	user	User context with DMDA, DMSwarm, etc.
[in]	source_field_id	Typed identity of the cell-centered Eulerian source.
[in]	target_field_id	Typed identity of the PETSC_REAL particle destination.

Returns

PetscErrorCode 0 on success.

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.

{
    PetscErrorCode ierr;
    FieldView source_view;
    const ParticleFieldDescriptor *target_descriptor = NULL;
 
    PetscFunctionBegin;
 
    ierr = FieldGetView(user, source_field_id, &source_view); CHKERRQ(ierr);
    ierr = ParticleFieldGetDescriptor(target_field_id, &target_descriptor); CHKERRQ(ierr);
    PetscCheck(source_view.descriptor->layout == FIELD_LAYOUT_CELL_CENTERED,
               PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
               "Eulerian-to-particle interpolation requires a cell-centered source; '%s' uses layout %s.",
               source_view.descriptor->canonical_name,
               FieldLayoutName(source_view.descriptor->layout));
    PetscCheck(target_descriptor->data_type == PETSC_REAL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
               "Eulerian interpolation destination '%s' must use PETSC_REAL storage.",
               target_descriptor->canonical_name);
    PetscCheck(source_view.descriptor->dof == target_descriptor->components,
               PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
               "Eulerian source '%s' has %d components but particle destination '%s' has %d.",
               source_view.descriptor->canonical_name, source_view.descriptor->dof,
               target_descriptor->canonical_name, target_descriptor->components);
 
    if (user->simCtx->interpolationMethod == INTERP_TRILINEAR) {
        ierr = InterpolateEulerFieldFromCenterToSwarm(user, source_view.local_vec,
                                                       source_view.descriptor->canonical_name,
                                                       target_descriptor->canonical_name); CHKERRQ(ierr);
    } else {
        ierr = InterpolateEulerFieldFromCornerToSwarm(user, source_view.local_vec,
                                                       source_view.descriptor->canonical_name,
                                                       target_descriptor->canonical_name); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

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 interpolates:

• 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]

user

Pointer 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.

Interpolates all relevant fields from the DMDA to the DMSwarm.

Local to this translation unit.

Definition at line 1400 of file interpolation.c.

{
  PetscErrorCode ierr;
  PetscMPIInt rank;
  PetscFunctionBegin;
 
  PROFILE_FUNCTION_BEGIN;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
 
  LOG_ALLOW(LOCAL, LOG_DEBUG,
             " Interpolation of ucat to velocity begins on rank %d.\n",rank);
  ierr = InterpolateEulerFieldToSwarm(user, FIELD_ID_UCAT,
                                      PARTICLE_FIELD_ID_VELOCITY); CHKERRQ(ierr);
  ierr = InterpolateEulerFieldToSwarm(user, FIELD_ID_DIFFUSIVITY,
                                      PARTICLE_FIELD_ID_DIFFUSIVITY); CHKERRQ(ierr);
  ierr = InterpolateEulerFieldToSwarm(user, FIELD_ID_DIFFUSIVITY_GRADIENT,
                                      PARTICLE_FIELD_ID_DIFFUSIVITY_GRADIENT); CHKERRQ(ierr);
  /* Add fields as necessary here*/
                                      
  ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
  LOG_ALLOW(LOCAL, LOG_INFO,
        "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InterpolateParticleVelocities()

PetscErrorCode InterpolateParticleVelocities

(

UserCtx *

user

)

Interpolates particle velocities using trilinear interpolation.

Parameters

[in]

user

Pointer to the user-defined context containing grid and swarm information.

Returns

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

InterpolateFieldFromCornerToCenter_Scalar()

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

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

For each cell center in the physical region of the cell-centered DM (fda), this function averages the up to 8 surrounding scalar values from the coordinate DM (da). On boundaries, where fewer corners are available, a partial average is computed.

The coordinate DM (da) is built on corners (IM+1 x JM+1 x KM+1) while the cell-centered DM (fda) covers the physical cells (IM x JM x KM). Index offsets are adjusted via DMDAGetLocalInfo.

Parameters

[in]	field_arr	3D array of corner-based scalar data (from user->da).
[out]	centfield_arr	3D array for interpolated cell-center scalar data (for user->fda).
[in]	user	User context containing: da : DM for the coordinate (corner) data. fda : DM for the cell-centered data.

Returns

PetscErrorCode 0 on success.

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.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    
    // Get local and global grid dimensions
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    // Determine loop bounds to compute for interior cells only, matching the code's convention.
    // Start at index 1 if the process owns the global boundary at index 0.
    PetscInt is = (xs == 0) ? 1 : xs;
    PetscInt js = (ys == 0) ? 1 : ys;
    PetscInt ks = (zs == 0) ? 1 : zs;
 
    // Stop one cell short if the process owns the global boundary at the max index.
    PetscInt ie = (xe == mx) ? xe - 1 : xe;
    PetscInt je = (ye == my) ? ye - 1 : ye;
    PetscInt ke = (ze == mz) ? ze - 1 : ze;
 
    // Loop over the locally owned INTERIOR cells.
    for (PetscInt k = ks; k < ke; k++) {
        for (PetscInt j = js; j < je; j++) {
            for (PetscInt i = is; i < ie; i++) {
                // Calculate cell center value as the average of its 8 corner nodes
                centfield_arr[k][j][i] = 0.125 * (field_arr[k][j][i] + field_arr[k][j-1][i] +
                                                  field_arr[k-1][j][i] + field_arr[k-1][j-1][i] +
                                                  field_arr[k][j][i-1] + field_arr[k][j-1][i-1] +
                                                  field_arr[k-1][j][i-1] + field_arr[k-1][j-1][i-1]);
            }
        }
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

InterpolateFieldFromCornerToCenter_Vector()

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

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

For each cell center in the physical region of the cell-centered DM (fda), this function averages the 8 surrounding corner values from the coordinate DM (da). The coordinate DM (da) is built on corners (IM+1 x JM+1 x KM+1) while the cell-centered DM (fda) covers the physical cells (IM x JM x KM). Index offsets are adjusted using DMDAGetLocalInfo.

Parameters

[in]	field_arr	3D array of corner-based vector data (from user->da).
[out]	centfield_arr	3D array for interpolated cell-center vector data (for user->fda).
[in]	user	User context containing: da : DM for the coordinate (corner) data. fda : DM for the cell-centered data.

Returns

PetscErrorCode 0 on success.

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.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Get local and global grid dimensions
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    // Determine loop bounds to compute for interior cells only, matching the code's convention.
    // Start at index 1 if the process owns the global boundary at index 0.
    PetscInt is = (xs == 0) ? 1 : xs;
    PetscInt js = (ys == 0) ? 1 : ys;
    PetscInt ks = (zs == 0) ? 1 : zs;
 
    // Stop one cell short if the process owns the global boundary at the max index.
    PetscInt ie = (xe == mx) ? xe - 1 : xe;
    PetscInt je = (ye == my) ? ye - 1 : ye;
    PetscInt ke = (ze == mz) ? ze - 1 : ze;
 
    // Loop over the locally owned INTERIOR cells.
    for (PetscInt k = ks; k < ke; k++) {
        for (PetscInt j = js; j < je; j++) {
            for (PetscInt i = is; i < ie; i++) {
                // Calculate cell center value as the average of its 8 corner nodes
                centfield_arr[k][j][i].x = 0.125 * (field_arr[k][j][i].x + field_arr[k][j-1][i].x +
                                                    field_arr[k-1][j][i].x + field_arr[k-1][j-1][i].x +
                                                    field_arr[k][j][i-1].x + field_arr[k][j-1][i-1].x +
                                                    field_arr[k-1][j][i-1].x + field_arr[k-1][j-1][i-1].x);
 
                centfield_arr[k][j][i].y = 0.125 * (field_arr[k][j][i].y + field_arr[k][j-1][i].y +
                                                    field_arr[k-1][j][i].y + field_arr[k-1][j-1][i].y +
                                                    field_arr[k][j][i-1].y + field_arr[k][j-1][i-1].y +
                                                    field_arr[k-1][j][i-1].y + field_arr[k-1][j-1][i-1].y);
 
                centfield_arr[k][j][i].z = 0.125 * (field_arr[k][j][i].z + field_arr[k][j-1][i].z +
                                                    field_arr[k-1][j][i].z + field_arr[k-1][j-1][i].z +
                                                    field_arr[k][j][i-1].z + field_arr[k][j-1][i-1].z +
                                                    field_arr[k-1][j][i-1].z + field_arr[k-1][j-1][i-1].z);
            }
        }
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TestCornerToCenterInterpolation()

PetscErrorCode TestCornerToCenterInterpolation

(

UserCtx *

user

)

Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.

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

1. Takes the corner-centered nodal coordinates (from DMGetCoordinatesLocal) as input.
2. Uses the InterpolateFieldFromCornerToCenter macro to interpolate these coordinates to the cell centers, storing the result in a new temporary vector.
3. Compares this new vector with the user->Cent vector, which is assumed to have been computed by ComputeCellCentersAndSpacing and serves as the ground truth.
4. A 2-norm of the difference is computed. If it is below a small tolerance, the test passes.

Note

This function should be called immediately after ComputeCellCentersAndSpacing has been successfully executed.

Parameters

user	The UserCtx for a specific grid level.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.

Local to this translation unit.

Definition at line 137 of file interpolation.c.

{
    PetscErrorCode ierr;
    Vec            lCoords, TestCent;
    Cmpnts         ***coor_arr, ***test_cent_arr;
    PetscReal      diff_norm;
 
    PetscFunctionBeginUser;
 
    // 1. Create a temporary vector to hold the result of our interpolation.
    //    It must have the same layout and size as the ground-truth user->Cent vector.
    ierr = VecDuplicate(user->Cent, &TestCent); CHKERRQ(ierr);
 
    // 2. Get the input (corner coordinates) and output (our test vector) arrays.
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
 
    // 3. Call the generic interpolation macro.
    //    The macro will see that `coor_arr` is of type `Cmpnts***` and correctly
    //    call the `InterpolateFieldFromCornerToCenter_Vector` function.
    ierr = InterpolateFieldFromCornerToCenter(coor_arr, test_cent_arr, user); CHKERRQ(ierr);
 
    // 4. Restore the arrays.
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
 
    // 5. IMPORTANT: Assemble the vector so its values are communicated across processors
    //    and it's ready for global operations like VecNorm.
    ierr = VecAssemblyBegin(TestCent); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(TestCent); CHKERRQ(ierr);
 
    // 6. Compare the result with the ground truth.
    //    We compute TestCent = -1.0 * user->Cent + 1.0 * TestCent.
    //    This calculates the difference vector: TestCent - user->Cent.
    ierr = VecAXPY(TestCent, -1.0, user->Cent); CHKERRQ(ierr);
 
    //    Now, compute the L2 norm of the difference vector. If the functions are
    //    identical, the norm should be zero (or very close due to floating point).
    ierr = VecNorm(TestCent, NORM_2, &diff_norm); CHKERRQ(ierr);
 
    // 7. Report the result and clean up.
    if (diff_norm < 1.0e-12) {
        LOG_ALLOW(GLOBAL,LOG_DEBUG,"[SUCCESS] Test passed. Norm of difference is %g.\n", (double)diff_norm);
    } else {
        LOG_ALLOW(GLOBAL,LOG_DEBUG, "[FAILURE] Test failed. Norm of difference is %g.\n", (double)diff_norm);
    }
 
    ierr = VecDestroy(&TestCent); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

InterpolateFieldFromCenterToCorner_Vector()

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

Interpolates a vector 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_arr	Input: 3D array (ghosted) of cell-centered data, accessed via GLOBAL indices.
[out]	corner_arr	Output: 3D array (ghosted) where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Interpolates a vector field from cell centers to corner nodes.

Local to this translation unit.

Definition at line 196 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
    PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
    PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
 
    PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
    PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
    PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
 
 
    // Global grid dimensions (used for valid cell check)
    PetscInt IM = info.mx - 1; // Total nodes in i-direction
    PetscInt JM = info.my - 1; // Total nodes in j-direction
    PetscInt KM = info.mz - 1; // Total nodes in k-direction
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
        "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
        rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
 
    // Loop over the GLOBAL indices of the NODES owned by this processor
    for (PetscInt k = zs_node; k < ze_node; k++) {
        for (PetscInt j = ys_node; j < ye_node; j++) {
            for (PetscInt i = xs_node; i < xe_node; i++) {
                Cmpnts sum = {0.0, 0.0, 0.0};
                PetscInt count = 0;
 
                // DEBUG 1 TEST
                /*
                if(rank == 1 && i == 24 && j == 12 && k == 49){
                    PetscInt i_cell_A = i - 1; 
                    PetscInt i_cell_B = i;
 
                    Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
                    Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
 
                    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",
                             rank, k, j, i,
                             k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
                             k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
 
                }
                */
                
                // Skip processing the unused last node in each dimension.
                if(i >= IM || j >= JM || k >= KM){
                    continue;
                }
                // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
                // 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
                for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
                    for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
                        for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
 
                            // These are still GLOBAL cell indices
                            PetscInt global_cell_k = k + dk_offset;
                            PetscInt global_cell_j = j + dj_offset;
                            PetscInt global_cell_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index
                            if (global_cell_i >= 0 && global_cell_i < nCellsX &&
                                global_cell_j >= 0 && global_cell_j < nCellsY &&
                                global_cell_k >= 0 && global_cell_k < nCellsZ)
                            {
                                    Cmpnts cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
 
                                    LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f, %.2f, %.2f)\n",
                                        rank,global_cell_k,global_cell_j,global_cell_i,cell_val.x, cell_val.y, cell_val.z);
                    
                                    sum.x += cell_val.x;
                                    sum.y += cell_val.y;
                                    sum.z += cell_val.z;
                                    count++;                                  
                            }
                        }
                    }
                }
 
                PetscInt i_global_write = i; // Global index in GLOBAL array.
                PetscInt j_global_write = j;
                PetscInt k_global_write = k;
                
                // We write directly into the array using the global loop indices.
                if (count > 0) {
                    corner_arr[k_global_write][j_global_write][i_global_write].x = sum.x / (PetscReal)count;
                    corner_arr[k_global_write][j_global_write][i_global_write].y = sum.y / (PetscReal)count;
                    corner_arr[k_global_write][j_global_write][i_global_write].z = sum.z / (PetscReal)count;
                } else {
                    // This case should ideally not happen for a valid owned node, but as a failsafe:
                    corner_arr[k_global_write][j_global_write][i_global_write] = (Cmpnts){0.0, 0.0, 0.0};
                }
               
                // DEBUG 2
                /*
                if(rank == 1){
                    if(i == 11 && j == 11 && k == 49){
                        Cmpnts ucat_node = corner_arr[k][j][i];
                        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",
                                 rank, k, j, i,
                                 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                    }
                }
 
                if(rank == 0 && i == 24 && j == 12 && k == 0){
                    Cmpnts ucat_node = corner_arr[k][j][i];
                    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",
                             rank, k, j, i,
                             k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                }
                */
              //  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);
            }
        }
    }
    PROFILE_FUNCTION_END;
    return 0;
}

InterpolateFieldFromCenterToCorner_Scalar()

PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar	(	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_arr	Input: 3D array (ghosted) of scalar data at cell centers, accessed via GLOBAL indices.
[out]	corner_arr	Output: 3D array (ghosted) where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Interpolates a scalar field from cell centers to corner nodes.

Local to this translation unit.

Definition at line 331 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
    PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
    PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
 
    PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
    PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
    PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
 
 
    // Global grid dimensions (used for valid cell check)
    PetscInt IM = info.mx - 1; // Total nodes in i-direction
    PetscInt JM = info.my - 1; // Total nodes in j-direction
    PetscInt KM = info.mz - 1; // Total nodes in k-direction
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
        "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
        rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
 
    // Loop over the GLOBAL indices of the NODES owned by this processor
    for (PetscInt k = zs_node; k < ze_node; k++) {
        for (PetscInt j = ys_node; j < ye_node; j++) {
            for (PetscInt i = xs_node; i < xe_node; i++) {
                PetscReal sum = 0.0;
                PetscInt count = 0;
 
                // DEBUG 1 TEST
                /*
                if(rank == 1 && i == 24 && j == 12 && k == 49){
                    PetscInt i_cell_A = i - 1; 
                    PetscInt i_cell_B = i;
 
                    Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
                    Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
 
                    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",
                             rank, k, j, i,
                             k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
                             k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
 
                }
                */
                
                // Skip processing the unused last node in each dimension.
                if(i >= IM || j >= JM || k >= KM){
                    continue;
                }
                // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
                // 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
                for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
                    for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
                        for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
 
                            // These are still GLOBAL cell indices
                            PetscInt global_cell_k = k + dk_offset;
                            PetscInt global_cell_j = j + dj_offset;
                            PetscInt global_cell_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index
                            if (global_cell_i >= 0 && global_cell_i < nCellsX &&
                                global_cell_j >= 0 && global_cell_j < nCellsY &&
                                global_cell_k >= 0 && global_cell_k < nCellsZ)
                            {
                                    PetscReal cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
 
                                    LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f)\n",
                                        rank,global_cell_k,global_cell_j,global_cell_i,cell_val);
                    
                                    sum += cell_val;
                                    count++;                                  
                            }
                        }
                    }
                }
 
                PetscInt i_global_write = i; // Global index in GLOBAL array.
                PetscInt j_global_write = j;
                PetscInt k_global_write = k;
                
                // We write directly into the array using the global loop indices.
                if (count > 0) {
                    corner_arr[k_global_write][j_global_write][i_global_write] = sum / (PetscReal)count;
                } else {
                    // This case should ideally not happen for a valid owned node, but as a failsafe:
                    corner_arr[k_global_write][j_global_write][i_global_write] = 0.0;
                }
               
                // DEBUG 2
                /*
                if(rank == 1){
                    if(i == 11 && j == 11 && k == 49){
                        Cmpnts ucat_node = corner_arr[k][j][i];
                        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",
                                 rank, k, j, i,
                                 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                    }
                }
 
                if(rank == 0 && i == 24 && j == 12 && k == 0){
                    Cmpnts ucat_node = corner_arr[k][j][i];
                    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",
                             rank, k, j, i,
                             k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                }
                */
              //  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);
            }
        }
    }
    PROFILE_FUNCTION_END;
    return 0;
}

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_arr	Ghosted node-centered array (global indexing) from user->fda.
[out]	faceX_arr	Local array for X-faces sized [zm][ym][xm+1].
[out]	faceY_arr	Local array for Y-faces sized [zm][ym+1][xm].
[out]	faceZ_arr	Local array for Z-faces sized [zm+1][ym][xm].
[in]	user	User context containing DMDA 'fda' and GetOwnedCellRange.

Returns

PetscErrorCode 0 on success, non-zero on failure.

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

Local to this translation unit.

Definition at line 1953 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Determine owned-cell ranges based on corner-node ownership
    PetscInt xs, xm, ys, ym, zs, zm;
    ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
 
    // Global exclusive end indices for cells
    PetscInt xe = xs + xm;
    PetscInt ye = ys + ym;
    PetscInt ze = zs + zm;
 
    // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i <= xe; ++i) {
                PetscInt i_loc = i - xs; // 0..xm
                // Average the four corners of the Y-Z face at X = i
                PetscReal sum = corner_arr[k  ][j  ][i]
                              + corner_arr[k+1][j  ][i]
                              + corner_arr[k  ][j+1][i]
                              + corner_arr[k+1][j+1][i];
                faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j <= ye; ++j) {
            PetscInt j_loc = j - ys; // 0..ym
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                // Average the four corners of the X-Z face at Y = j
                PetscReal sum = corner_arr[k  ][j][i  ]
                              + corner_arr[k+1][j][i  ]
                              + corner_arr[k  ][j][i+1]
                              + corner_arr[k+1][j][i+1];
                faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
    for (PetscInt k = zs; k <= ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                // Average the four corners of the X-Y face at Z = k
                PetscReal sum = corner_arr[k][j  ][i  ]
                              + corner_arr[k][j  ][i+1]
                              + corner_arr[k][j+1][i  ]
                              + corner_arr[k][j+1][i+1];
                faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

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_arr	Ghosted 3-component array (global node indices).
[out]	faceX_arr	Local array of Cmpnts for X-faces sized [zm][ym][xm+1].
[out]	faceY_arr	Local array of Cmpnts for Y-faces sized [zm][ym+1][xm].
[out]	faceZ_arr	Local array of Cmpnts for Z-faces sized [zm+1][ym][xm].
[in]	user	User context containing DMDA 'fda'.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 2043 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info;
    PetscMPIInt rank;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
    
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    PetscInt xs, xm, ys, ym, zs, zm;
    ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
 
    PetscInt xe = xs + xm;
    PetscInt ye = ys + ym;
    PetscInt ze = zs + zm;
 
    // X-faces
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i <= xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k  ][j  ][i].x + corner_arr[k+1][j  ][i].x
                      + corner_arr[k  ][j+1][i].x + corner_arr[k+1][j+1][i].x;
                sum.y = corner_arr[k  ][j  ][i].y + corner_arr[k+1][j  ][i].y
                      + corner_arr[k  ][j+1][i].y + corner_arr[k+1][j+1][i].y;
                sum.z = corner_arr[k  ][j  ][i].z + corner_arr[k+1][j  ][i].z
                      + corner_arr[k  ][j+1][i].z + corner_arr[k+1][j+1][i].z;
                faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d x-face Interpolation complete.\n", rank);
    
    // Y-faces
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j <= ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k  ][j][i  ].x + corner_arr[k+1][j][i  ].x
                      + corner_arr[k  ][j][i+1].x + corner_arr[k+1][j][i+1].x;
                sum.y = corner_arr[k  ][j][i  ].y + corner_arr[k+1][j][i  ].y
                      + corner_arr[k  ][j][i+1].y + corner_arr[k+1][j][i+1].y;
                sum.z = corner_arr[k  ][j][i  ].z + corner_arr[k+1][j][i  ].z
                      + corner_arr[k  ][j][i+1].z + corner_arr[k+1][j][i+1].z;
                faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d y-face Interpolation complete.\n", rank);
    
    // Z-faces
    for (PetscInt k = zs; k <= ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k][j  ][i  ].x + corner_arr[k][j  ][i+1].x
                      + corner_arr[k][j+1][i  ].x + corner_arr[k][j+1][i+1].x;
                sum.y = corner_arr[k][j  ][i  ].y + corner_arr[k][j  ][i+1].y
                      + corner_arr[k][j+1][i  ].y + corner_arr[k][j+1][i+1].y;
                sum.z = corner_arr[k][j  ][i  ].z + corner_arr[k][j  ][i+1].z
                      + corner_arr[k][j+1][i  ].z + corner_arr[k][j+1][i+1].z;
                faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d z-face Interpolation complete.\n", rank);
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function: