Metric.h File Reference

#include <petsc.h>
#include "variables.h"
#include "logging.h"
#include <stdlib.h>
#include "io.h"
#include "setup.h"

Include dependency graph for Metric.h:

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

Go to the source code of this file.

Functions
PetscErrorCode	MetricLogicalToPhysical (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
	Maps a logical point inside one hexahedral cell to physical space.
PetscErrorCode	MetricGetCellVertices (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
	Collects the eight node coordinates of one logical hexahedral cell.
PetscErrorCode	MetricJacobian (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, PetscReal J[3][3], PetscReal *detJ)
	Evaluates the trilinear mapping Jacobian and its determinant in one cell.
PetscErrorCode	MetricVelocityContravariant (const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
	Converts a Cartesian velocity vector to contravariant logical components.
PetscErrorCode	CalculateFaceNormalAndArea (Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
	Computes the unit normal vectors and areas of the three faces of a computational cell.
PetscErrorCode	InvertCovariantMetricTensor (double covariantTensor[3][3], double contravariantTensor[3][3])
	Inverts the 3x3 matrix of face-area vectors to obtain the covariant directions.
PetscErrorCode	ComputeCellCharacteristicLengthScale (PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
	Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.
PetscErrorCode	ComputeCellDirectionalExtents (PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *l_xi, double *l_eta, double *l_zeta)
	Computes a cell's extent along each of its own grid directions.
PetscErrorCode	ComputeCellEdgeVectors (PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, Cmpnts edges[3])
	Computes a cell's edge vectors along its three grid directions.
PetscErrorCode	ApplyPeriodicCorrectionsToCellCentersAndSpacing (UserCtx *user)
	Builds translated periodic images for cell centers and grid spacing.
PetscErrorCode	ApplyPeriodicCorrectionsToIFaceCenter (UserCtx *user)
	Builds translated periodic images for i-face centers (Centx).
PetscErrorCode	ApplyPeriodicCorrectionsToJFaceCenter (UserCtx *user)
	Builds translated periodic images for j-face centers (Centy).
PetscErrorCode	ApplyPeriodicCorrectionsToKFaceCenter (UserCtx *user)
	Builds translated periodic images for k-face centers (Centz).
PetscErrorCode	ComputeFaceMetrics (UserCtx *user)
	Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation, and stores them in the corresponding global Vec members of the UserCtx structure (user->Csi, user->Eta, user->Zet).
PetscErrorCode	ComputeCellCenteredJacobianInverse (UserCtx *user)
	Calculates the cell-centered inverse Jacobian determinant (1/J) for INTERIOR cells and stores it in user->Aj.
PetscErrorCode	CheckAndFixGridOrientation (UserCtx *user)
	Verify and, when consistently inverted, repair the right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain.
PetscErrorCode	ComputeCellCentersAndSpacing (UserCtx *user)
	Computes the physical location of cell centers and the spacing between them.
PetscErrorCode	ComputeIFaceMetrics (UserCtx *user)
	Computes metrics centered on constant-i faces (i-faces).
PetscErrorCode	ComputeJFaceMetrics (UserCtx *user)
	Computes metrics centered on constant-j faces (j-faces).
PetscErrorCode	ComputeKFaceMetrics (UserCtx *user)
	Computes metrics centered on constant-k faces (k-faces).
PetscErrorCode	ComputeMetricsDivergence (UserCtx *user)
	Performs a diagnostic check on the divergence of the face area metric vectors.
PetscErrorCode	ComputeMetricNorms (UserCtx *user)
	Computes the max-min values of the grid metrics.
PetscErrorCode	CalculateAllGridMetrics (SimCtx *simCtx)
	Orchestrates the calculation of all grid metrics.

Function Documentation

MetricLogicalToPhysical()

PetscErrorCode MetricLogicalToPhysical	(	UserCtx *	user,
		const Cmpnts ***	X,
		PetscInt	i,
		PetscInt	j,
		PetscInt	k,
		PetscReal	xi,
		PetscReal	eta,
		PetscReal	zta,
		Cmpnts *	Xp
	)

Maps a logical point inside one hexahedral cell to physical space.

Uses trilinear interpolation of the eight cell vertices. Logical coordinates are measured from the cell's lower logical corner and normally lie in [0, 1].

Parameters

	user	Grid context used for cell indexing.
	X	Ghosted node-coordinate array for the block.
	i	Cell index in the xi direction.
	j	Cell index in the eta direction.
	k	Cell index in the zeta direction.
	xi	Local xi coordinate within the cell.
	eta	Local eta coordinate within the cell.
	zta	Local zeta coordinate within the cell.
[out]	Xp	Physical Cartesian position evaluated at the logical point.

Returns

PetscErrorCode 0 on success.

Maps a logical point inside one hexahedral cell to physical space.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Metric.h.

See also

MetricLogicalToPhysical()

Definition at line 76 of file Metric.c.

{
  PetscErrorCode ierr;
  Cmpnts V[8];
  PetscFunctionBeginUser;
  
  PROFILE_FUNCTION_BEGIN;
 
  ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
  TrilinearBlend(V,xi,eta,zta,Xp);
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MetricGetCellVertices()

PetscErrorCode MetricGetCellVertices	(	UserCtx *	user,
		const Cmpnts ***	X,
		PetscInt	i,
		PetscInt	j,
		PetscInt	k,
		Cmpnts	V[8]
	)

Collects the eight node coordinates of one logical hexahedral cell.

Vertex ordering matches the trilinear metric kernels in this module.

Parameters

	user	Grid context used for cell indexing.
	X	Ghosted node-coordinate array for the block.
	i	Cell index in the xi direction.
	j	Cell index in the eta direction.
	k	Cell index in the zeta direction.
[out]	V	Eight physical vertex coordinates in metric-kernel order.

Returns

PetscErrorCode 0 on success.

Collects the eight node coordinates of one logical hexahedral cell.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Metric.h.

See also

MetricGetCellVertices()

Definition at line 26 of file Metric.c.

{
  (void)user;
  PetscFunctionBeginUser;
  for (PetscInt c = 0; c < 8; ++c) {
    PetscInt ii = i + ((c & 1) ? 1 : 0);
    PetscInt jj = j + ((c & 2) ? 1 : 0);
    PetscInt kk = k + ((c & 4) ? 1 : 0);
    LOG_LOOP_ALLOW(GLOBAL, LOG_VERBOSE,i+j+k,10," ii: %d,jj:%d,kk:%d - Retrieved.\n",ii,jj,kk);
    V[c] = X[kk][jj][ii];
  }
  PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MetricJacobian()

PetscErrorCode MetricJacobian	(	UserCtx *	user,
		const Cmpnts ***	X,
		PetscInt	i,
		PetscInt	j,
		PetscInt	k,
		PetscReal	xi,
		PetscReal	eta,
		PetscReal	zta,
		PetscReal	J[3][3],
		PetscReal *	detJ
	)

Evaluates the trilinear mapping Jacobian and its determinant in one cell.

The returned matrix is \(\partial(x,y,z)/\partial(\xi,\eta,\zeta)\); callers use its determinant to reject degenerate or inverted cells.

Parameters

	user	Grid context used for cell indexing.
	X	Ghosted node-coordinate array for the block.
	i	Cell index in the xi direction.
	j	Cell index in the eta direction.
	k	Cell index in the zeta direction.
	xi	Local xi coordinate within the cell.
	eta	Local eta coordinate within the cell.
	zta	Local zeta coordinate within the cell.
[out]	J	Mapping Jacobian in row-major Cartesian/logical form.
[out]	detJ	Determinant of J.

Returns

PetscErrorCode 0 on success.

Evaluates the trilinear mapping Jacobian and its determinant in one cell.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Metric.h.

See also

MetricJacobian()

Definition at line 105 of file Metric.c.

{
  PetscErrorCode ierr;
  Cmpnts V[8];
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
 
  /* derivatives of trilinear shape functions */
  PetscReal dN_dXi[8], dN_dEta[8], dN_dZta[8];
  for (PetscInt c=0;c<8;++c) {
    PetscReal sx = (c & 1) ?  1.0 : -1.0;
    PetscReal sy = (c & 2) ?  1.0 : -1.0;
    PetscReal sz = (c & 4) ?  1.0 : -1.0;
    dN_dXi [c] = 0.125 * sx * ( (c&2?eta:1-eta) ) * ( (c&4?zta:1-zta) );
    dN_dEta[c] = 0.125 * sy * ( (c&1?xi :1-xi ) ) * ( (c&4?zta:1-zta) );
    dN_dZta[c] = 0.125 * sz * ( (c&1?xi :1-xi ) ) * ( (c&2?eta:1-eta) );
  }
 
  /* assemble Jacobian */
  PetscReal x_xi=0,y_xi=0,z_xi=0,
            x_eta=0,y_eta=0,z_eta=0,
            x_zta=0,y_zta=0,z_zta=0;
  for (PetscInt c=0;c<8;++c) {
    x_xi  += dN_dXi [c]*V[c].x;  y_xi  += dN_dXi [c]*V[c].y;  z_xi  += dN_dXi [c]*V[c].z;
    x_eta += dN_dEta[c]*V[c].x;  y_eta += dN_dEta[c]*V[c].y;  z_eta += dN_dEta[c]*V[c].z;
    x_zta += dN_dZta[c]*V[c].x;  y_zta += dN_dZta[c]*V[c].y;  z_zta += dN_dZta[c]*V[c].z;
  }
 
  J[0][0]=x_xi;  J[0][1]=x_eta;  J[0][2]=x_zta;
  J[1][0]=y_xi;  J[1][1]=y_eta;  J[1][2]=y_zta;
  J[2][0]=z_xi;  J[2][1]=z_eta;  J[2][2]=z_zta;
 
  if (detJ) {
    *detJ = x_xi*(y_eta*z_zta - y_zta*z_eta)
          - x_eta*(y_xi*z_zta - y_zta*z_xi)
          + x_zta*(y_xi*z_eta - y_eta*z_xi);
  }
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

MetricVelocityContravariant()

PetscErrorCode MetricVelocityContravariant	(	const PetscReal	J[3][3],
		PetscReal	detJ,
		const PetscReal	u[3],
		PetscReal	uc[3]
	)

Converts a Cartesian velocity vector to contravariant logical components.

Uses the inverse of the physical-to-logical Jacobian; detJ must describe the same mapping as J and must be nonzero.

Parameters

	J	Mapping Jacobian at the evaluation point.
	detJ	Determinant of J.
	u	Cartesian velocity components.
[out]	uc	Contravariant xi/eta/zeta velocity components.

Returns

PetscErrorCode 0 on success.

Converts a Cartesian velocity vector to contravariant logical components.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Metric.h.

See also

MetricVelocityContravariant()

Definition at line 165 of file Metric.c.

{
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
  PetscReal gxi[3]  = {  J[1][1]*J[2][2]-J[1][2]*J[2][1],
                        -J[0][1]*J[2][2]+J[0][2]*J[2][1],
                         J[0][1]*J[1][2]-J[0][2]*J[1][1] };
  PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
                         J[0][0]*J[2][2]-J[0][2]*J[2][0],
                        -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
  PetscReal gzta[3] = {  J[1][0]*J[2][1]-J[1][1]*J[2][0],
                        -J[0][0]*J[2][1]+J[0][1]*J[2][0],
                         J[0][0]*J[1][1]-J[0][1]*J[1][0] };
 
  PetscReal invDet = 1.0 / detJ;
  for (int d=0; d<3; ++d) { gxi[d]  *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
 
  uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
  uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
  uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the caller graph for this function:

CalculateFaceNormalAndArea()

PetscErrorCode CalculateFaceNormalAndArea	(	Cmpnts	csi,
		Cmpnts	eta,
		Cmpnts	zet,
		double	ni[3],
		double	nj[3],
		double	nk[3],
		double *	Ai,
		double *	Aj,
		double *	Ak
	)

Computes the unit normal vectors and areas of the three faces of a computational cell.

Given the metric vectors (csi, eta, zet), this function calculates the geometric properties of the cell faces aligned with the i, j, and k directions.

Parameters

	csi	Covariant xi metric vector.
	eta	Covariant eta metric vector.
	zet	Covariant zeta metric vector.
[out]	ni	Unit normal of the xi-normal face.
[out]	nj	Unit normal of the eta-normal face.
[out]	nk	Unit normal of the zeta-normal face.
[out]	Ai	Area of the xi-normal face.
[out]	Aj	Area of the eta-normal face.
[out]	Ak	Area of the zeta-normal face.

Deprecated:

Unreachable from any executable: its only caller is ComputeCellCharacteristicLengthScale, which has none. Kept for now with its unit tests; do not add new callers - use ComputeCellDirectionalExtents or ComputeCellEdgeVectors, which do not depend on the cell's orientation.

Returns

PetscErrorCode 0 on success.

Computes the unit normal vectors and areas of the three faces of a computational cell.

Local to this translation unit.

Definition at line 251 of file Metric.c.

{
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    double g[3][3];
    double G[3][3];
    
    g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
    g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
    g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
    
    PetscCall(InvertCovariantMetricTensor(g, G));
    double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
    double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
    double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
          
    double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
    double nx_j = xeta, ny_j = yeta, nz_j = zeta;
    double nx_k = xzet, ny_k = yzet, nz_k = zzet;
          
    double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
    double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
    double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
 
    *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] );  // area
    *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
    *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
        
    nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
    nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
    nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
 
    ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
    nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
    nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InvertCovariantMetricTensor()

PetscErrorCode InvertCovariantMetricTensor	(	double	covariantTensor[3][3],
		double	contravariantTensor[3][3]
	)

Inverts the 3x3 matrix of face-area vectors to obtain the covariant directions.

The caller fills the rows with csi, eta and zet - the face-area vectors, which are the contravariant basis scaled by the cell Jacobian. The inverse's columns are the covariant (tangent) directions dx/dxi, dx/deta, dx/dzeta divided by that Jacobian, which is what CalculateFaceNormalAndArea() normalizes.

Degeneracy is judged relative to the rows' own magnitude, because the determinant of area vectors scales as the square of the cell volume: an absolute floor would reject a well-formed cell purely for being small.

Parameters

	covariantTensor	Input: rows are the three face-area vectors.
[out]	contravariantTensor	Inverse written in place.

Deprecated:

Unreachable from any executable: its only caller is CalculateFaceNormalAndArea, which is reached only from the deprecated ComputeCellCharacteristicLengthScale. Kept for now with its unit tests; do not add new callers - use ComputeCellDirectionalExtents or ComputeCellEdgeVectors, which do not depend on the cell's orientation.

Returns

PetscErrorCode 0 on success; PETSC_ERR_MAT_LU_ZRPVT when the rows are coplanar or one of them vanishes.

Inverts the 3x3 matrix of face-area vectors to obtain the covariant directions.

Local to this translation unit.

Definition at line 202 of file Metric.c.

{
    PetscFunctionBeginUser;
 
    const double a11=covariantTensor[0][0], a12=covariantTensor[0][1], a13=covariantTensor[0][2];
    const double a21=covariantTensor[1][0], a22=covariantTensor[1][1], a23=covariantTensor[1][2];
    const double a31=covariantTensor[2][0], a32=covariantTensor[2][1], a33=covariantTensor[2][2];
 
    double det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
 
    /* Degeneracy has to be judged against the rows' own magnitude, not against a fixed
       number. The rows here are face-area vectors, so the determinant carries their
       product - for an orthogonal cell it is exactly the square of the cell volume, and
       it shrinks as the mesh is refined while the cell stays perfectly well formed. An
       absolute floor therefore rejects fine grids: a wall-resolved cell of 5e-9 volume
       has a determinant near 3e-17. The ratio below is 1 for an orthogonal cell and
       approaches 0 only as the rows become coplanar, which is the actual failure. */
    const double row_scale = sqrt(a11*a11 + a12*a12 + a13*a13) *
                             sqrt(a21*a21 + a22*a22 + a23*a23) *
                             sqrt(a31*a31 + a32*a32 + a33*a33);
 
    if (row_scale <= 0.0 || fabs(det) <= 1.0e-10 * row_scale) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT,
                "Matrix is singular: |det| = %g is degenerate against the rows' own scale "
                "%g (ratio %g); the rows are coplanar or one of them vanishes.",
                (double)fabs(det), (double)row_scale,
                (double)(row_scale > 0.0 ? fabs(det)/row_scale : 0.0));
    }
 
    contravariantTensor[0][0] = (a33*a22-a32*a23)/det;
    contravariantTensor[0][1] = -(a33*a12-a32*a13)/det;
    contravariantTensor[0][2] = (a23*a12-a22*a13)/det;
    contravariantTensor[1][0] = -(a33*a21-a31*a23)/det;
    contravariantTensor[1][1] = (a33*a11-a31*a13)/det;
    contravariantTensor[1][2] = -(a23*a11-a21*a13)/det;
    contravariantTensor[2][0] = (a32*a21-a31*a22)/det;
    contravariantTensor[2][1] = -(a32*a11-a31*a12)/det;
    contravariantTensor[2][2] = (a22*a11-a21*a12)/det;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ComputeCellCharacteristicLengthScale()

PetscErrorCode ComputeCellCharacteristicLengthScale	(	PetscReal	ajc,
		Cmpnts	csi,
		Cmpnts	eta,
		Cmpnts	zet,
		double *	dx,
		double *	dy,
		double *	dz
	)

Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.

For a non-uniform, non-orthogonal cell, there is no single "dx". This function computes an effective length scale in each Cartesian direction based on the cell volume and the areas of its faces: the three outputs are the Cartesian components of the cell diagonal, so they change when an otherwise identical cell is rotated. Use ComputeCellDirectionalExtents for a measure of the cell itself.

Parameters

	ajc	Cell Jacobian/volume metric at the cell center.
	csi	Covariant xi metric vector.
	eta	Covariant eta metric vector.
	zet	Covariant zeta metric vector.
[out]	dx	Effective physical length associated with xi variation.
[out]	dy	Effective physical length associated with eta variation.
[out]	dz	Effective physical length associated with zeta variation.

Deprecated:

No production caller. Its outputs are Cartesian projections of the cell diagonal, so they change when an identical cell is rotated; the LES filter widths and the Clark term stopped using it for that reason. Kept for now with its unit tests; do not add new callers - use ComputeCellDirectionalExtents or ComputeCellEdgeVectors, which do not depend on the cell's orientation.

Returns

PetscErrorCode 0 on success.

Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.

Local to this translation unit.

Definition at line 297 of file Metric.c.

{
        PetscFunctionBeginUser;
        PROFILE_FUNCTION_BEGIN;
        double ni[3], nj[3], nk[3];
        double Li, Lj, Lk;
        double Ai, Aj, Ak;
        double vol = 1./ajc;
 
        PetscCall(CalculateFaceNormalAndArea(csi, eta, zet, ni, nj, nk, &Ai, &Aj, &Ak));
        Li = vol / Ai;
        Lj = vol / Aj;
        Lk = vol / Ak;
 
        // Length scale vector = di * ni_vector + dj * nj_vector + dk * nk_vector
        *dx = fabs( Li * ni[0] + Lj * nj[0] + Lk * nk[0] );
        *dy = fabs( Li * ni[1] + Lj * nj[1] + Lk * nk[1] );
        *dz = fabs( Li * ni[2] + Lj * nj[2] + Lk * nk[2] );
 
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeCellDirectionalExtents()

PetscErrorCode ComputeCellDirectionalExtents	(	PetscReal	ajc,
		Cmpnts	csi,
		Cmpnts	eta,
		Cmpnts	zet,
		double *	l_xi,
		double *	l_eta,
		double *	l_zeta
	)

Computes a cell's extent along each of its own grid directions.

Each extent is the cell volume divided by the area of the face pair it crosses, which is the distance between those two faces. Unlike ComputeCellCharacteristicLengthScale, whose outputs are the Cartesian components of the cell diagonal and so change when the same cell is rotated, these depend on the cell alone: a filter width built from them is the same for a cell in a straight section and for an identical cell turned by a bend.

Parameters

	ajc	Inverse cell volume at the cell center.
	csi	Face-area vector of the xi faces.
	eta	Face-area vector of the eta faces.
	zet	Face-area vector of the zeta faces.
[out]	l_xi	Extent across the xi faces.
[out]	l_eta	Extent across the eta faces.
[out]	l_zeta	Extent across the zeta faces.

Returns

PetscErrorCode 0 on success; PETSC_ERR_ARG_OUTOFRANGE for a non-positive Jacobian or a face of zero area.

Computes a cell's extent along each of its own grid directions.

Full API contract is documented with the header declaration in include/Metric.h.

See also

ComputeCellDirectionalExtents()

Definition at line 328 of file Metric.c.

{
    const double area_xi   = sqrt(csi.x*csi.x + csi.y*csi.y + csi.z*csi.z);
    const double area_eta  = sqrt(eta.x*eta.x + eta.y*eta.y + eta.z*eta.z);
    const double area_zeta = sqrt(zet.x*zet.x + zet.y*zet.y + zet.z*zet.z);
 
    PetscFunctionBeginUser;
    PetscCheck(ajc > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Cell Jacobian must be positive to define cell extents; received %g.", (double)ajc);
    PetscCheck(area_xi > 0.0 && area_eta > 0.0 && area_zeta > 0.0, PETSC_COMM_SELF,
               PETSC_ERR_ARG_OUTOFRANGE,
               "A cell face has zero area (%g, %g, %g), so the cell has no extent across it.",
               area_xi, area_eta, area_zeta);
 
    /* Volume over the area of a face pair is the distance between those faces: the
       cell's extent in that grid direction. It uses nothing but the cell's own metrics,
       so it is the same wherever and however the cell sits in space. */
    *l_xi   = (1.0/ajc)/area_xi;
    *l_eta  = (1.0/ajc)/area_eta;
    *l_zeta = (1.0/ajc)/area_zeta;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ComputeCellEdgeVectors()

PetscErrorCode ComputeCellEdgeVectors	(	PetscReal	ajc,
		Cmpnts	csi,
		Cmpnts	eta,
		Cmpnts	zet,
		Cmpnts	edges[3]
	)

Computes a cell's edge vectors along its three grid directions.

Each is the covariant basis vector dx/dxi_k over one index step: the physical displacement from one cell to the next along that grid line, with both its direction and its length. Projecting a velocity gradient onto it gives the change in velocity across the cell in that direction, which is what a direction-resolved subgrid model needs and what does not depend on how the grid is oriented.

Parameters

	ajc	Inverse cell volume at the cell center.
	csi	Face-area vector of the xi faces.
	eta	Face-area vector of the eta faces.
	zet	Face-area vector of the zeta faces.
[out]	edges	Edge vectors along xi, eta and zeta, in that order.

Returns

PetscErrorCode 0 on success; PETSC_ERR_ARG_OUTOFRANGE for a non-positive Jacobian.

Computes a cell's edge vectors along its three grid directions.

Full API contract is documented with the header declaration in include/Metric.h.

See also

ComputeCellEdgeVectors()

Definition at line 360 of file Metric.c.

{
    PetscFunctionBeginUser;
    PetscCheck(ajc > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Cell Jacobian must be positive to define cell edges; received %g.", (double)ajc);
    /* The face-area vectors are J grad(xi), J grad(eta), J grad(zeta), with J the cell
       volume. The covariant basis dx/dxi equals J grad(eta) x grad(zeta), so it is
       (eta x zet) / J - the cell's edge along xi, direction and length together. */
    edges[0].x = ajc * (eta.y*zet.z - eta.z*zet.y);
    edges[0].y = ajc * (eta.z*zet.x - eta.x*zet.z);
    edges[0].z = ajc * (eta.x*zet.y - eta.y*zet.x);
    edges[1].x = ajc * (zet.y*csi.z - zet.z*csi.y);
    edges[1].y = ajc * (zet.z*csi.x - zet.x*csi.z);
    edges[1].z = ajc * (zet.x*csi.y - zet.y*csi.x);
    edges[2].x = ajc * (csi.y*eta.z - csi.z*eta.y);
    edges[2].y = ajc * (csi.z*eta.x - csi.x*eta.z);
    edges[2].z = ajc * (csi.x*eta.y - csi.y*eta.x);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ApplyPeriodicCorrectionsToCellCentersAndSpacing()

PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing

(

UserCtx *

user

)

Builds translated periodic images for cell centers and grid spacing.

PETSc wraps field indices but does not translate coordinates. This routine applies the validated per-axis geometric translation to wrapped center coordinates and refreshes the associated local grid spacing.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Builds translated periodic images for cell centers and grid spacing.

Local to this translation unit.

Definition at line 444 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    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;
    Cmpnts       ***cent, ***lcent, ***gs;
    PetscReal      delta;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // Check if any periodic boundaries exist
    PetscBool has_periodic = PETSC_FALSE;
    for (int i = 0; i < 6; i++) {
        if (user->boundary_faces[i].mathematical_type == PERIODIC) {
            has_periodic = PETSC_TRUE;
            break;
        }
    }
 
    if (!has_periodic) {
        LOG_ALLOW(LOCAL, LOG_TRACE, "No periodic boundaries; skipping corrections for Cent/GridSpace.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying periodic corrections to Cent and GridSpace.\n");
 
    // Must update ghosts first before applying corrections
    ierr = UpdateLocalGhosts(user, FIELD_ID_CENT); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_GRID_SPACE); CHKERRQ(ierr);
 
    // --- X-direction periodic corrections ---
    if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || 
        user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
        
        ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
 
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
            if (user->cgrid) {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt j=ys; j<ye; j++) {
                        cent[k][j][0] = lcent[k][j][-2];
                    }
                }
            } else {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt j=ys; j<ye; j++) {
                        delta = (gs[k][j][1].x + gs[k][j][-2].x) / 2.0;
                        cent[k][j][0].x = cent[k][j][1].x - delta;
                        cent[k][j][0].y = cent[k][j][1].y;
                        cent[k][j][0].z = cent[k][j][1].z;
                    }
                }
            }
        }
 
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
            if (user->cgrid) {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt j=ys; j<ye; j++) {
                        cent[k][j][mx-1] = lcent[k][j][mx+1];
                    }
                }
            } else {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt j=ys; j<ye; j++) {
                        delta = (gs[k][j][mx-2].x + gs[k][j][mx+1].x) / 2.0;
                        cent[k][j][mx-1].x = cent[k][j][mx-2].x + delta;
                        cent[k][j][mx-1].y = cent[k][j][mx-2].y;
                        cent[k][j][mx-1].z = cent[k][j][mx-2].z;
                    }
                }
            }
        }
 
        ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
    }    
    // --- Y-direction periodic corrections ---
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || 
        user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
        
        ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
 
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
            if (user->cgrid) {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        cent[k][0][i] = lcent[k][-2][i];
                    }
                }
            } else {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        delta = (gs[k][1][i].y + gs[k][-2][i].y) / 2.0;
                        cent[k][0][i].x = cent[k][1][i].x;
                        cent[k][0][i].y = cent[k][1][i].y - delta;
                        cent[k][0][i].z = cent[k][1][i].z;
                    }
                }
            }
        }
 
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
            if (user->cgrid) {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        cent[k][my-1][i] = lcent[k][my+1][i];
                    }
                }
            } else {
                for (PetscInt k=zs; k<ze; k++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        delta = (gs[k][my-2][i].y + gs[k][my+1][i].y) / 2.0;
                        cent[k][my-1][i].x = cent[k][my-2][i].x;
                        cent[k][my-1][i].y = cent[k][my-2][i].y + delta;
                        cent[k][my-1][i].z = cent[k][my-2][i].z;
                    }
                }
            }
        }
 
        ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
 
    }
 
    // --- Z-direction periodic corrections ---
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || 
        user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
        
        ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
 
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
            if (user->cgrid) {
                for (PetscInt j=ys; j<ye; j++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        cent[0][j][i] = lcent[-2][j][i];
                    }
                }
            } else {
                for (PetscInt j=ys; j<ye; j++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        delta = (gs[1][j][i].z + gs[-2][j][i].z) / 2.0;
                        cent[0][j][i].x = cent[1][j][i].x;
                        cent[0][j][i].y = cent[1][j][i].y;
                        cent[0][j][i].z = cent[1][j][i].z - delta;
                    }
                }
            }
        }
 
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
            if (user->cgrid) {
                for (PetscInt j=ys; j<ye; j++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        cent[mz-1][j][i] = lcent[mz+1][j][i];
                    }
                }
            } else {
                for (PetscInt j=ys; j<ye; j++) {
                    for (PetscInt i=xs; i<xe; i++) {
                        delta = (gs[mz-2][j][i].z + gs[mz+1][j][i].z) / 2.0;
                        cent[mz-1][j][i].x = cent[mz-2][j][i].x;
                        cent[mz-1][j][i].y = cent[mz-2][j][i].y;
                        cent[mz-1][j][i].z = cent[mz-2][j][i].z + delta;
                    }
                }
            }
        }
 
        ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
 
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyPeriodicCorrectionsToIFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter

(

UserCtx *

user

)

Builds translated periodic images for i-face centers (Centx).

Every active periodic axis can affect the ghost images of Centx. This must be called after Centx is computed and before it is used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Builds translated periodic images for i-face centers (Centx).

Local to this translation unit.

Definition at line 643 of file Metric.c.

{
    const FieldId fields[] = {FIELD_ID_CENTX};
 
    PetscFunctionBeginUser;
    PetscCall(SynchronizePeriodicFaceFields(user, 'i', 1, fields));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

ApplyPeriodicCorrectionsToJFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter

(

UserCtx *

user

)

Builds translated periodic images for j-face centers (Centy).

Every active periodic axis can affect the ghost images of Centy. This must be called after Centy is computed and before it is used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Builds translated periodic images for j-face centers (Centy).

Local to this translation unit.

Definition at line 659 of file Metric.c.

{
    const FieldId fields[] = {FIELD_ID_CENTY};
 
    PetscFunctionBeginUser;
    PetscCall(SynchronizePeriodicFaceFields(user, 'j', 1, fields));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

ApplyPeriodicCorrectionsToKFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter

(

UserCtx *

user

)

Builds translated periodic images for k-face centers (Centz).

Every active periodic axis can affect the ghost images of Centz. This must be called after Centz is computed and before it is used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Builds translated periodic images for k-face centers (Centz).

Local to this translation unit.

Definition at line 675 of file Metric.c.

{
    const FieldId fields[] = {FIELD_ID_CENTZ};
 
    PetscFunctionBeginUser;
    PetscCall(SynchronizePeriodicFaceFields(user, 'k', 1, fields));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

ComputeFaceMetrics()

PetscErrorCode ComputeFaceMetrics

(

UserCtx *

user

)

Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation, and stores them in the corresponding global Vec members of the UserCtx structure (user->Csi, user->Eta, user->Zet).

This is a self-contained routine that performs the following steps:

1. Obtains local ghosted nodal coordinates using DMGetCoordinatesLocal.
2. Calculates metrics for INTERIOR faces where finite difference stencils are valid.
3. EXTRAPOLATES metrics for faces on the physical domain boundaries by copying from the nearest computed interior face.
4. Assembles the global user->Csi, user->Eta, user->Zet Vecs.
5. Updates the local ghosted user->lCsi, user->lEta, user->lZet Vecs.

Parameters

[in,out]

user

Pointer to the UserCtx structure.

Returns

PetscErrorCode 0 on success.

Note

• This function is a complete "compute and make ready" unit for Csi, Eta, and Zet.
• It's recommended to call VecZeroEntries on user->Csi, Eta, Zet before this if they might contain old data.

Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation, and stores them in the corresponding global Vec members of the UserCtx structure (user->Csi, user->Eta, user->Zet).

Local to this translation unit.

Definition at line 690 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Cmpnts       ***csi_arr, ***eta_arr, ***zet_arr;
    Cmpnts       ***nodal_coords_arr;
    Vec            localCoords_from_dm;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
    ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
    if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
    ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
 
    // --- 2. Get arrays for output global Vecs from UserCtx ---
    ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
 
    // Define owned node ranges (global indices)
    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;
 
    // Global domain dimensions (total number of nodes)
    PetscInt mx = info.mx;
    PetscInt my = info.my;
    PetscInt mz = info.mz;
 
    // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
    // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
    PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
    PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
    PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
 
   // These represent the surface area of the curvilinear cell face and the normal rotated such that the direction of increasing coordinate is maintained.
   // The metric vectors (Csi, Eta, Zet) are defined to point in the direction of their corresponding increasing computational coordinate.
 
    // Calculate Csi
    for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
        for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
            for (PetscInt i_node = xs; i_node < xe; ++i_node) {
          
                PetscReal dx_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node].x - nodal_coords_arr[k_node][j_node-1][i_node].x - nodal_coords_arr[k_node-1][j_node-1][i_node].x);
                PetscReal dy_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node].y - nodal_coords_arr[k_node][j_node-1][i_node].y - nodal_coords_arr[k_node-1][j_node-1][i_node].y);
                PetscReal dz_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node].z - nodal_coords_arr[k_node][j_node-1][i_node].z - nodal_coords_arr[k_node-1][j_node-1][i_node].z);
                PetscReal dx_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node][j_node][i_node].x - nodal_coords_arr[k_node-1][j_node-1][i_node].x - nodal_coords_arr[k_node-1][j_node][i_node].x);
                PetscReal dy_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node][j_node][i_node].y - nodal_coords_arr[k_node-1][j_node-1][i_node].y - nodal_coords_arr[k_node-1][j_node][i_node].y);
                PetscReal dz_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node][j_node][i_node].z - nodal_coords_arr[k_node-1][j_node-1][i_node].z - nodal_coords_arr[k_node-1][j_node][i_node].z);
 
        csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
                csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
                csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
            }
        }
    }
 
    // Calculate Eta
    for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
        for (PetscInt j_node = ys; j_node < ye; ++j_node) {
            for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
          
                PetscReal dx_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node].x - nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node-1][j_node][i_node-1].x);
                PetscReal dy_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node].y - nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node-1][j_node][i_node-1].y);
                PetscReal dz_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node].z - nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node-1][j_node][i_node-1].z);
                PetscReal dx_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node-1][j_node][i_node].x - nodal_coords_arr[k_node-1][j_node][i_node-1].x);
                PetscReal dy_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node-1][j_node][i_node].y - nodal_coords_arr[k_node-1][j_node][i_node-1].y);
                PetscReal dz_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node-1][j_node][i_node].z - nodal_coords_arr[k_node-1][j_node][i_node-1].z);
 
        eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi  - dz_dzeta * dy_dxi;
                eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi  - dx_dzeta * dz_dxi;
                eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi  - dy_dzeta * dx_dxi;
            }
        }
    }
 
    // Calculate Zet
    for (PetscInt k_node = zs; k_node < ze; ++k_node) {
        for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
            for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
          
                PetscReal dx_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node].x - nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node][j_node-1][i_node-1].x);
                PetscReal dy_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node].y - nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node][j_node-1][i_node-1].y);
                PetscReal dz_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node].z - nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node][j_node-1][i_node-1].z);
                PetscReal dx_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node][j_node-1][i_node].x - nodal_coords_arr[k_node][j_node-1][i_node-1].x);
                PetscReal dy_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node][j_node-1][i_node].y - nodal_coords_arr[k_node][j_node-1][i_node-1].y);
                PetscReal dz_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node][j_node-1][i_node].z - nodal_coords_arr[k_node][j_node-1][i_node-1].z);
 
        zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
                zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
                zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
            }
        }
    }
 
    // --- 4. Boundary Extrapolation ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
    PetscInt i_bnd, j_bnd, k_bnd;
 
    if (xs == 0) { // If this rank owns the global i=0 boundary
        i_bnd = 0;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
                if (i_bnd + 1 < mx) {
                    eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
                    zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
                }
            }
        }
    }
    if (xe == mx) { // If this rank owns the global i=mx-1 boundary
        i_bnd = mx - 1;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
                if (i_bnd - 1 >= 0) {
                    eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
                    zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
                }
            }
        }
    }
    if (ys == 0) {
        j_bnd = 0;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                if (j_bnd + 1 < my) {
                    csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
                    zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
                }
            }
        }
    }
    if (ye == my) {
        j_bnd = my - 1;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                if (j_bnd - 1 >= 0) {
                    csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
                    zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
                }
            }
        }
    }
    if (zs == 0) {
        k_bnd = 0;
        for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                if (k_bnd + 1 < mz) {
                    csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
                    eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
                }
            }
        }
    }
    if (ze == mz) {
        k_bnd = mz - 1;
        for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                if (k_bnd - 1 >= 0) {
                    csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
                    eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
                }
            }
        }
    }
 
    if (info.xs==0 && info.ys==0 && info.zs==0) {
      PetscReal dot = zet_arr[0][0][0].z;  /* dot with global +z */
      LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f   (should be >0 for right-handed grid)\n", dot);
    }
    
    // --- 5. Restore all arrays ---
    ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
 
    // --- 6. Assemble Global Vectors ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
    ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
 
    // --- 7. Update Local Ghosted Versions ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
    ierr = UpdateLocalGhosts(user, FIELD_ID_CSI); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_ETA); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_ZET); CHKERRQ(ierr);
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeCellCenteredJacobianInverse()

PetscErrorCode ComputeCellCenteredJacobianInverse

(

UserCtx *

user

)

Calculates the cell-centered inverse Jacobian determinant (1/J) for INTERIOR cells and stores it in user->Aj.

This version includes boundary extrapolation.

Nodal coordinates are obtained internally. Refer to previous Doxygen comments for details on physical locations and storage convention (aj_arr[k_n][j_n][i_n] for cell C(i_n-1,j_n-1,k_n-1)).

Parameters

[in,out]

user

Pointer to the UserCtx structure.

Returns

PetscErrorCode 0 on success.

Calculates the cell-centered inverse Jacobian determinant (1/J) for INTERIOR cells and stores it in user->Aj.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/Metric.h.

See also

ComputeCellCenteredJacobianInverse()

Definition at line 902 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscScalar  ***aj_arr;
    Cmpnts       ***nodal_coords_arr;
    Vec            localCoords_from_dm;
 
    PetscFunctionBeginUser;
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
 
    // --- 1. Get Nodal Coordinates and Output Array ---
    ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
 
    // Define owned node ranges (global indices)
    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;
 
    // Global domain dimensions (total number of nodes)
    PetscInt mx = info.mx;
    PetscInt my = info.my;
    PetscInt mz = info.mz;
 
    // --- 2. Calculate Aj for INTERIOR Stencils ---
    
    PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
    PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
    PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
 
    PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
    PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
    PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
 
    for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
        for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
            for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
          
                PetscReal dx_dxi = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node-1][i_node].x) - (nodal_coords_arr[k_node][j_node][i_node-1].x + nodal_coords_arr[k_node][j_node-1][i_node-1].x + nodal_coords_arr[k_node-1][j_node][i_node-1].x + nodal_coords_arr[k_node-1][j_node-1][i_node-1].x) );
 
        PetscReal dy_dxi = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node-1][i_node].y) - (nodal_coords_arr[k_node][j_node][i_node-1].y + nodal_coords_arr[k_node][j_node-1][i_node-1].y + nodal_coords_arr[k_node-1][j_node][i_node-1].y + nodal_coords_arr[k_node-1][j_node-1][i_node-1].y) );
 
        PetscReal dz_dxi = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node-1][i_node].z) - (nodal_coords_arr[k_node][j_node][i_node-1].z + nodal_coords_arr[k_node][j_node-1][i_node-1].z + nodal_coords_arr[k_node-1][j_node][i_node-1].z + nodal_coords_arr[k_node-1][j_node-1][i_node-1].z) );
 
        PetscReal dx_deta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x + nodal_coords_arr[k_node-1][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node-1].x) - (nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node-1].x + nodal_coords_arr[k_node-1][j_node-1][i_node].x + nodal_coords_arr[k_node-1][j_node-1][i_node-1].x) );
 
        PetscReal dy_deta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y + nodal_coords_arr[k_node-1][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node-1].y) - (nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node-1].y + nodal_coords_arr[k_node-1][j_node-1][i_node].y + nodal_coords_arr[k_node-1][j_node-1][i_node-1].y) );
 
        PetscReal dz_deta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z + nodal_coords_arr[k_node-1][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node-1].z) - (nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node-1].z + nodal_coords_arr[k_node-1][j_node-1][i_node].z + nodal_coords_arr[k_node-1][j_node-1][i_node-1].z) );
 
        PetscReal dx_dzeta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x + nodal_coords_arr[k_node][j_node-1][i_node-1].x) - (nodal_coords_arr[k_node-1][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node-1][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node-1].x + nodal_coords_arr[k_node-1][j_node-1][i_node-1].x) );
 
        PetscReal dy_dzeta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y + nodal_coords_arr[k_node][j_node-1][i_node-1].y) - (nodal_coords_arr[k_node-1][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node-1][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node-1].y + nodal_coords_arr[k_node-1][j_node-1][i_node-1].y) );
 
        PetscReal dz_dzeta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z + nodal_coords_arr[k_node][j_node-1][i_node-1].z) - (nodal_coords_arr[k_node-1][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node-1][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node-1].z + nodal_coords_arr[k_node-1][j_node-1][i_node-1].z) );
 
        PetscReal jacobian_det = dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) - dy_dxi * (dx_deta * dz_dzeta - dz_deta * dx_dzeta) + dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta);
                if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
                aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
            }
        }
    }
 
    // --- 4. Boundary Extrapolation for Aj ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
    PetscInt i_bnd, j_bnd, k_bnd;
 
    if (xs == 0) {
        i_bnd = 0;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
                if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
            }
        }
    }
    if (xe == mx) {
        i_bnd = mx - 1;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
                if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
            }
        }
    }
    // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
    if (ys == 0) {
        j_bnd = 0;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
            }
        }
    }
    if (ye == my) {
        j_bnd = my - 1;
        for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
            }
        }
    }
    if (zs == 0) {
        k_bnd = 0;
        for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
            }
        }
    }
    if (ze == mz) {
        k_bnd = mz - 1;
        for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
            for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
                 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
            }
        }
    }
 
    // --- 5. Restore arrays ---
    ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
 
    // --- 6. Assemble Global Vector ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
    ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
 
    // --- 7. Update Local Ghosted Version ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
    ierr = UpdateLocalGhosts(user, FIELD_ID_AJ); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CheckAndFixGridOrientation()

PetscErrorCode CheckAndFixGridOrientation

(

UserCtx *

user

)

Verify and, when consistently inverted, repair the right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain.

The metric-generation kernels are completely algebraic, so they will happily deliver a left-handed basis if the mesh file enumerates nodes in the opposite ζ-direction.
This routine makes the orientation explicit and—if needed—repairs it once per run:

Step	Action
1	Compute global Aj_min, Aj_max.
2	Mixed signs (Aj_min < 0 && Aj_max > 0) → abort: the mesh is topologically inconsistent.
3	All negative (Aj_max < 0) → abort: the grid is left-handed. Negating the metric vectors is not a valid repair (a mirrored duct solved that way was wrong by a factor of forty), so the grid must be renumbered along one logical axis instead.
4	Store user->GridOrientation = +1, the only orientation the solver accepts.

Parameters

[in,out]

user

Fully initialised UserCtx that already contains Csi, Eta, Zet, Aj, their local ghosts, and valid distributed DMs.

Returns

0 on success or a PETSc error code on failure.

Note

Call immediately after ComputeCellCenteredJacobianInverse() and before any routine that differentiates or applies BCs.

Author metadata intentionally omitted in API docs.

Verify and, when consistently inverted, repair the right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain.

Local to this translation unit.

Definition at line 389 of file Metric.c.

{
    PetscErrorCode ierr;
    PetscReal      aj_min, aj_max;
    PetscMPIInt    rank;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    /* ---------------- step 1: global extrema of Aj ---------------- */
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
 
    ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);  /* already global */
    ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
        "[orientation] Global Aj range: [%.3e , %.3e]\n",
        (double)aj_min, (double)aj_max);
 
    /* ---------------- step 2: detect malformed mesh ---------------- */
    if (aj_min < 0.0 && aj_max > 0.0)
        SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
            "Mixed Jacobian signs detected – grid is topologically inconsistent.");
 
    /* A uniformly left-handed grid is refused. The former "repair" negated Csi, Eta,
       Zet and Aj on this level only, left the face and cell-centre metrics as they
       were, and then reset the orientation flag to +1; a mirrored duct solved that way
       reached forty times its bulk velocity. Renumbering one logical axis is exact. */
    PetscInt orientation = +1;
    if (aj_max < 0.0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER,
                "The grid is left-handed: every cell Jacobian is negative, as a mirrored or "
                "oddly permuted grid is. The solver needs right-handed logical axes. Renumber "
                "one logical axis (grid.gen: append reverse:axis=i|j|k to the transforms); "
                "that axis's two faces trade names.");
    }
 
    /* ---------------- step 4: store result in UserCtx -------------- */
    user->GridOrientation = orientation;
 
    if (!rank)
        LOG_ALLOW(LOCAL, LOG_INFO, "[orientation] Grid confirmed right-handed.\n");
 
    PROFILE_FUNCTION_END;        
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ComputeCellCentersAndSpacing()

PetscErrorCode ComputeCellCentersAndSpacing

(

UserCtx *

user

)

Computes the physical location of cell centers and the spacing between them.

This function calculates two key geometric properties from the nodal coordinates:

1. Cent: A vector field storing the (x,y,z) coordinates of the center of each grid cell.
2. GridSpace: A vector field storing the physical distance between adjacent cell centers in the i, j, and k computational directions.

It is a direct adaptation of the corresponding logic from the legacy FormMetrics.

Parameters

user	The UserCtx for a specific grid level. The function populates user->Cent and user->GridSpace.

Returns

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

Computes the physical location of cell centers and the spacing between them.

Local to this translation unit.

Definition at line 1046 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Vec            lCoords;
    const Cmpnts ***coor;
    Cmpnts       ***cent, ***gs;
    PetscReal    xcp, ycp, zcp, xcm, ycm, zcm;
    PetscInt     xs,ys,zs,xe,ye,ze,mx,my,mz;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing cell centers and spacing for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
 
    ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
 
    xs = info.xs; xe = info.xs + info.xm;
    ys = info.ys; ye = info.ys + info.ym;
    zs = info.zs; ze = info.zs + info.zm;
    mx = info.mx; my =  info.my; mz = info.mz;
 
    PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
    PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
    PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
 
    PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
    PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
    PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;    
    
    // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
    for (PetscInt k=k_start_node; k<k_end_node; k++) {
        for (PetscInt j=j_start_node; j<j_end_node; j++) {
            for (PetscInt i=i_start_node; i<i_end_node; i++) {
                // Calculate cell center as the average of its 8 corner nodes
                cent[k][j][i].x = 0.125 * (coor[k][j][i].x + coor[k][j-1][i].x + coor[k-1][j][i].x + coor[k-1][j-1][i].x + coor[k][j][i-1].x + coor[k][j-1][i-1].x + coor[k-1][j][i-1].x + coor[k-1][j-1][i-1].x);
                cent[k][j][i].y = 0.125 * (coor[k][j][i].y + coor[k][j-1][i].y + coor[k-1][j][i].y + coor[k-1][j-1][i].y + coor[k][j][i-1].y + coor[k][j-1][i-1].y + coor[k-1][j][i-1].y + coor[k-1][j-1][i-1].y);
                cent[k][j][i].z = 0.125 * (coor[k][j][i].z + coor[k][j-1][i].z + coor[k-1][j][i].z + coor[k-1][j-1][i].z + coor[k][j][i-1].z + coor[k][j-1][i-1].z + coor[k-1][j][i-1].z + coor[k-1][j-1][i-1].z);
 
                // Calculate Grid Spacing in i-direction (distance between i-face centers)
                xcp = 0.25 * (coor[k][j][i].x + coor[k][j-1][i].x + coor[k-1][j-1][i].x + coor[k-1][j][i].x);
                ycp = 0.25 * (coor[k][j][i].y + coor[k][j-1][i].y + coor[k-1][j-1][i].y + coor[k-1][j][i].y);
                zcp = 0.25 * (coor[k][j][i].z + coor[k][j-1][i].z + coor[k-1][j-1][i].z + coor[k-1][j][i].z);
                xcm = 0.25 * (coor[k][j][i-1].x + coor[k][j-1][i-1].x + coor[k-1][j-1][i-1].x + coor[k-1][j][i-1].x);
                ycm = 0.25 * (coor[k][j][i-1].y + coor[k][j-1][i-1].y + coor[k-1][j-1][i-1].y + coor[k-1][j][i-1].y);
                zcm = 0.25 * (coor[k][j][i-1].z + coor[k][j-1][i-1].z + coor[k-1][j-1][i-1].z + coor[k-1][j][i-1].z);
                gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
       
                // Calculate Grid Spacing in j-direction (distance between j-face centers)
                xcp = 0.25 * (coor[k][j][i].x + coor[k][j][i-1].x + coor[k-1][j][i].x + coor[k-1][j][i-1].x);
                ycp = 0.25 * (coor[k][j][i].y + coor[k][j][i-1].y + coor[k-1][j][i].y + coor[k-1][j][i-1].y);
                zcp = 0.25 * (coor[k][j][i].z + coor[k][j][i-1].z + coor[k-1][j][i].z + coor[k-1][j][i-1].z);
                xcm = 0.25 * (coor[k][j-1][i].x + coor[k][j-1][i-1].x + coor[k-1][j-1][i].x + coor[k-1][j-1][i-1].x);
                ycm = 0.25 * (coor[k][j-1][i].y + coor[k][j-1][i-1].y + coor[k-1][j-1][i].y + coor[k-1][j-1][i-1].y);
                zcm = 0.25 * (coor[k][j-1][i].z + coor[k][j-1][i-1].z + coor[k-1][j-1][i].z + coor[k-1][j-1][i-1].z);
                gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
 
                // Calculate Grid Spacing in k-direction (distance between k-face centers)
                xcp = 0.25 * (coor[k][j][i].x + coor[k][j][i-1].x + coor[k][j-1][i].x + coor[k][j-1][i-1].x);
                ycp = 0.25 * (coor[k][j][i].y + coor[k][j][i-1].y + coor[k][j-1][i].y + coor[k][j-1][i-1].y);
                zcp = 0.25 * (coor[k][j][i].z + coor[k][j][i-1].z + coor[k][j-1][i].z + coor[k][j-1][i-1].z);
                xcm = 0.25 * (coor[k-1][j][i].x + coor[k-1][j][i-1].x + coor[k-1][j-1][i].x + coor[k-1][j-1][i-1].x);
                ycm = 0.25 * (coor[k-1][j][i].y + coor[k-1][j][i-1].y + coor[k-1][j-1][i].y + coor[k-1][j-1][i-1].y);
                zcm = 0.25 * (coor[k-1][j][i].z + coor[k-1][j-1][i-1].z + coor[k-1][j-1][i].z + coor[k-1][j-1][i-1].z);
                gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
            }
        }
    }
    
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
 
    // Assemble and update ghost regions for the new data
    ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_CENT); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_GRID_SPACE); CHKERRQ(ierr);
 
    ierr = ApplyPeriodicCorrectionsToCellCentersAndSpacing(user); CHKERRQ(ierr);
 
    // Final assembly and ghost update after corrections
    ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); 
    ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_CENT); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeIFaceMetrics()

PetscErrorCode ComputeIFaceMetrics

(

UserCtx *

user

)

Computes metrics centered on constant-i faces (i-faces).

This function calculates the metric terms (ICsi, IEta, IZet) and the inverse Jacobian (IAj) located at the center of each i-face. The stencils use i-face-centered coordinates (Centx) which must be computed first. The logic is a direct adaptation of the legacy FormMetrics function.

Parameters

user	The UserCtx for a specific grid level. Populates user->ICsi, etc.

Returns

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

Computes metrics centered on constant-i faces (i-faces).

Local to this translation unit.

Definition at line 1149 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Vec            lCoords;
    const Cmpnts ***coor;
    Cmpnts       ***centx; //***gs;
    const Cmpnts ***centx_const;
    Cmpnts       ***icsi, ***ieta, ***izet;
    PetscScalar  ***iaj;
    PetscReal      dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
    PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
    PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
    PetscInt lxe = xe;
    PetscInt lys = ys; PetscInt lye = ye;
    PetscInt lzs = zs; PetscInt lze = ze;
 
    if (ys==0) lys = ys+1;
    if (zs==0) lzs = zs+1;
 
    if (xe==mx) lxe=xe-1;
    if (ye==my) lye=ye-1;
    if (ze==mz) lze=ze-1;
 
    // --- Part 1: Calculate the location of i-face centers (Centx) ---
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
    //  ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
 
    //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Calculating i-face centers (Centx) with i[%d,%d], j[%d,%d], k[%d,%d] ...\n", user->simCtx->rank,gxs,gxe,gys,gye,gzs,gze);
 
    // Populate only owned physical face centers. Periodic endpoint and ghost
    // coordinates are established by the canonical face-field synchronizer.
    for (PetscInt k = PetscMax(zs, 1); k < PetscMin(ze, mz - 1); k++) {
        for (PetscInt j = PetscMax(ys, 1); j < PetscMin(ye, my - 1); j++) {
            for (PetscInt i = xs; i < PetscMin(xe, mx - 1); i++) {
                //----- DEBUG ------
                //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:     Calculating i-face center at (k=%d, j=%d, i=%d)\n", user->simCtx->rank, k, j, i);
                //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:       Using corner nodes: (%f,%f,%f), (%f,%f,%f), (%f,%f,%f), (%f,%f,%f)\n", user->simCtx->rank,
                //          coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
                //          coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
                //          coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
                //          coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
 
                centx[k][j][i].x = 0.25 * (coor[k][j][i].x + coor[k-1][j][i].x + coor[k][j-1][i].x + coor[k-1][j-1][i].x);
                centx[k][j][i].y = 0.25 * (coor[k][j][i].y + coor[k-1][j][i].y + coor[k][j-1][i].y + coor[k-1][j-1][i].y);
                centx[k][j][i].z = 0.25 * (coor[k][j][i].z + coor[k-1][j][i].z + coor[k][j-1][i].z + coor[k-1][j-1][i].z);
            
                //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:       Calculated i-face center: (%f,%f,%f)\n", user->simCtx->rank, centx[k][j][i].x, centx[k][j][i].y, centx[k][j][i].z);
            }
        }
    }
 
    //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
    /*
    if(xs==0){
      for(PetscInt k=gzs+1;k < gze; k++){
    for(PetscInt j=gys+1;j < gye; j++){
      PetscInt i=0; 
      centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
      centx[k][j][i-1].y=centx[k][j][i].y;
      centx[k][j][i-1].z=centx[k][j][i].z;
    }
      }
    }
    if (xe==mx){
      for(PetscInt k=gzs+1; k<gze; k++) {
    for (PetscInt j=gys+1; j<gye;j++) {
      PetscInt i=mx-1; 
      centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
      centx[k][j][i].y=centx[k][j][i-1].y;
      centx[k][j][i].z=centx[k][j][i-1].z;
    }
      }
    }
    */
    
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
 
    // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
 
    {
        const FieldId face_centers[] = {FIELD_ID_CENTX};
        ierr = SynchronizePeriodicFaceFields(user, 'i', 1, face_centers); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
 
    // --- Part 2: Calculate metrics using face-centered coordinates ---
    ierr = DMDAVecGetArrayRead(user->fda, user->lCentx, &centx_const); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
 
    // Loop over the OWNED region where we will store the final metrics
    for (PetscInt k=lzs; k<lze; k++) {
        for (PetscInt j=lys; j<lye; j++) {
            for (PetscInt i=xs; i<lxe; i++) {
 
                // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
                if(i == 0 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC){
                        // Forward difference at the domain's min-i boundary
                        dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
                        dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
                        dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
                } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) { 
                    // Backward difference at the domain's max-i boundary
                    dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
                    dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
                    dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
                } else { // Central difference in the interior  (or if PERIODIC BCs)
                    dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
                    dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
                    dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
                }
 
                // --- Stencil Logic for d/deta (derivative in j-direction) ---
                if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) { 
                    // Forward difference
                    dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
                    dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
                    dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
                } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) { 
                    // Backward difference
                    dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
                    dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
                    dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
                } else { // Central difference (interior or PERIODIC)
                    dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
                    dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
                    dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
                }
 
                // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
                if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) { 
                    // Forward difference
                    dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
                    dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
                    dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
                } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) { 
                    // Backward difference
                    dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
                    dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
                    dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
                } else { // Central difference (Interior + PERIODIC)
                    dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
                    dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
                    dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
                }
 
                // --- Metric calculations (identical to legacy FormMetrics) ---
                icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
                icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
                icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
 
                ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
                ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
                ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
 
                izet[k][j][i].x = dydc * dzde - dzdc * dyde;
                izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
                izet[k][j][i].z = dxdc * dyde - dydc * dxde;
 
                iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
                if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
                    iaj[k][j][i] = 1.0 / iaj[k][j][i];
                }
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentx, &centx_const); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
    
    // --- Part 3: Assemble global vectors and update local ghosts ---
    ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
 
    ierr = UpdateLocalGhosts(user, FIELD_ID_ICSI); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_IETA); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_IZET); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_IAJ); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeJFaceMetrics()

PetscErrorCode ComputeJFaceMetrics

(

UserCtx *

user

)

Computes metrics centered on constant-j faces (j-faces).

This function calculates the metric terms (JCsi, JEta, JZet) and the inverse Jacobian (JAj) located at the geometric center of each constant-j face. This is a critical step for staggered-grid finite difference schemes.

The process is a direct and faithful refactoring of the corresponding logic from the legacy FormMetrics function:

1. It first calculates the physical (x,y,z) coordinates of the center of each i-face and stores them in the user->Centy vector.
2. It then uses a boundary-aware, second-order finite difference stencil on the Centy field to compute the derivatives (e.g., d(x)/d(csi)).
   • Central differences are used in the grid interior.
   • One-sided differences are used at the physical domain boundaries.
3. Finally, these derivatives are used to compute the final metric terms and the inverse Jacobian, which are stored in their respective Vec objects.

Parameters

user	The UserCtx for a specific grid level. This function populates the user->JCsi, user->JEta, user->JZet, and user->JAj vectors.

Returns

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

Computes metrics centered on constant-j faces (j-faces).

Local to this translation unit.

Definition at line 1360 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Vec            lCoords;
    const Cmpnts ***coor;
    Cmpnts       ***centy; //***gs;
    const Cmpnts ***centy_const;
    Cmpnts       ***jcsi, ***jeta, ***jzet;
    PetscScalar  ***jaj;
    PetscReal      dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
    PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
    PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
    PetscInt lxs = xs; PetscInt lxe = xe;
    PetscInt lye = ye;
    PetscInt lzs = zs; PetscInt lze = ze;
 
    if (xs==0) lxs = xs+1;
    if (zs==0) lzs = zs+1;
 
    if (xe==mx) lxe=xe-1;
    if (ye==my) lye=ye-1;
    if (ze==mz) lze=ze-1;
 
    // --- Part 1: Calculate the location of i-face centers (Centx) ---
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
    //  ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
 
    for (PetscInt k = PetscMax(zs, 1); k < PetscMin(ze, mz - 1); k++) {
        for (PetscInt j = ys; j < PetscMin(ye, my - 1); j++) {
            for (PetscInt i = PetscMax(xs, 1); i < PetscMin(xe, mx - 1); i++) {
                centy[k][j][i].x = 0.25 * (coor[k][j][i].x + coor[k-1][j][i].x + coor[k][j][i-1].x + coor[k-1][j][i-1].x);
                centy[k][j][i].y = 0.25 * (coor[k][j][i].y + coor[k-1][j][i].y + coor[k][j][i-1].y + coor[k-1][j][i-1].y);
                centy[k][j][i].z = 0.25 * (coor[k][j][i].z + coor[k-1][j][i].z + coor[k][j][i-1].z + coor[k-1][j][i-1].z);
            }
        }
    }
 
    /*
    if(ys==0){
      for(PetscInt k=gzs+1;k < gze; k++){
    for(PetscInt i=gxs+1;j < gxe; i++){
      PetscInt j=0; 
      centy[k][j-1][i].x=centy[k][j][i].x;
      centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
      centy[k][j-1][i].z=centy[k][j][i].z;
    }
      }
    }
    if (ye==my){
      for(PetscInt k=gzs+1; k<gze; k++) {
    for (PetscInt i=gxs+1; j<gxe;i++) {
      PetscInt j=my-1; 
      centy[k][j][i].x=centy[k][j-1][i].x
      centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
      centy[k][j][i].z=centy[k][j-1][i].z;
    }
      }
    }
    */
    
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
    // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
 
    {
        const FieldId face_centers[] = {FIELD_ID_CENTY};
        ierr = SynchronizePeriodicFaceFields(user, 'j', 1, face_centers); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
 
    // --- Part 2: Calculate metrics using face-centered coordinates ---
    ierr = DMDAVecGetArrayRead(user->fda, user->lCenty, &centy_const); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
 
    // Loop over the OWNED region where we will store the final metrics
    for (PetscInt k=lzs; k<lze; k++) {
        for (PetscInt j=ys; j<lye; j++) {
            for (PetscInt i=lxs; i<lxe; i++) {
 
                // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
                if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) { 
                    // Forward difference at the domain's min-i boundary
                    dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
                    dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
                    dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
                } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) { 
                    // Backward difference at the domain's max-i boundary
                    dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
                    dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
                    dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
                } else { // Central difference in the interior or PERIODIC
                    dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
                    dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
                    dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
                }
 
                // --- Stencil Logic for d/deta (derivative in j-direction) ---
                if (j == 0 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) { 
                    // Forward difference
                    dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
                    dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
                    dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
                } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
                    // Backward difference
                    dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
                    dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
                    dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
                } else { // Central difference (interior or PERIODIC)
                    dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
                    dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
                    dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
                }
 
                // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
                if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) { 
                    // Forward difference
                    dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
                    dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
                    dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
                } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) { 
                    // Backward difference
                    dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
                    dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
                    dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
                } else { // Central difference (Interior or PERIODIC)
                    dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
                    dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
                    dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
                }
 
                // --- Metric calculations (identical to legacy FormMetrics) ---
                jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
                jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
                jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
 
                jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
                jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
                jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
 
                jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
                jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
                jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
 
                jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
                if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
                    jaj[k][j][i] = 1.0 / jaj[k][j][i];
                }
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCenty, &centy_const); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
    
    // --- Part 3: Assemble global vectors and update local ghosts ---
    ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
 
    ierr = UpdateLocalGhosts(user, FIELD_ID_JCSI); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_JETA); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_JZET); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_JAJ); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeKFaceMetrics()

PetscErrorCode ComputeKFaceMetrics

(

UserCtx *

user

)

Computes metrics centered on constant-k faces (k-faces).

This function calculates the metric terms (KCsi, KEta, KZet) and the inverse Jacobian (KAj) located at the geometric center of each constant-j face. This is a critical step for staggered-grid finite difference schemes.

The process is a direct and faithful refactoring of the corresponding logic from the legacy FormMetrics function:

1. It first calculates the physical (x,y,z) coordinates of the center of each i-face and stores them in the user->Centz vector.
2. It then uses a boundary-aware, second-order finite difference stencil on the Centz field to compute the derivatives (e.g., d(x)/d(csi)).
   • Central differences are used in the grid interior.
   • One-sided differences are used at the physical domain boundaries.
3. Finally, these derivatives are used to compute the final metric terms and the inverse Jacobian, which are stored in their respective Vec objects.

Parameters

user	The UserCtx for a specific grid level. This function populates the user->KCsi, user->KEta, user->KZet, and user->KAj vectors.

Returns

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

Computes metrics centered on constant-k faces (k-faces).

Local to this translation unit.

Definition at line 1555 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Vec            lCoords;
    const Cmpnts ***coor;
    Cmpnts       ***centz; //***gs;
    const Cmpnts ***centz_const;
    Cmpnts       ***kcsi, ***keta, ***kzet;
    PetscScalar  ***kaj;
    PetscReal      dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
    PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
    PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
    PetscInt lxs = xs; PetscInt lxe = xe;
    PetscInt lys = ys; PetscInt lye = ye;
    PetscInt lze = ze;
 
    if (xs==0) lxs = xs+1;
    if (ys==0) lys = ys+1;
 
    if (xe==mx) lxe=xe-1;
    if (ye==my) lye=ye-1;
    if (ze==mz) lze=ze-1;
 
    // --- Part 1: Calculate the location of i-face centers (Centx) ---
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
    //  ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
 
    for (PetscInt k = zs; k < PetscMin(ze, mz - 1); k++) {
        for (PetscInt j = PetscMax(ys, 1); j < PetscMin(ye, my - 1); j++) {
            for (PetscInt i = PetscMax(xs, 1); i < PetscMin(xe, mx - 1); i++) {
                centz[k][j][i].x = 0.25 * (coor[k][j][i].x + coor[k][j-1][i].x + coor[k][j][i-1].x + coor[k][j-1][i-1].x);
                centz[k][j][i].y = 0.25 * (coor[k][j][i].y + coor[k][j-1][i].y + coor[k][j][i-1].y + coor[k][j-1][i-1].y);
                centz[k][j][i].z = 0.25 * (coor[k][j][i].z + coor[k][j-1][i].z + coor[k][j][i-1].z + coor[k][j-1][i-1].z);
            }
        }
    }
 
    /*
    if(zs==0){
      for(PetscInt j=gys+1;j < gye; j++){
    for(PetscInt i=gxs+1;j < gxe; i++){
      PetscInt k=0; 
      centz[k-1][j][i].x=centz[k][j][i].x;
      centz[k-1][j][i].y=centz[k][j][i].y;
      centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
    }
      }
    }
    if (ze==mz){
      for(PetscInt j=gys+1; j<gye; j++) {
    for (PetscInt i=gxs+1; j<gxe;i++) {
      PetscInt k=mz-1; 
      centy[k][j][i].x=centy[k-1][j][i].x
      centy[k][j][i].y=centy[k-1][j][i].y;
      centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
    }
      }
    }
    */
    
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
    // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
 
    {
        const FieldId face_centers[] = {FIELD_ID_CENTZ};
        ierr = SynchronizePeriodicFaceFields(user, 'k', 1, face_centers); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
 
    // --- Part 2: Calculate metrics using face-centered coordinates ---
    ierr = DMDAVecGetArrayRead(user->fda, user->lCentz, &centz_const); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
 
    // Loop over the OWNED region where we will store the final metrics
    for (PetscInt k=zs; k<lze; k++) {
        for (PetscInt j=lys; j<lye; j++) {
            for (PetscInt i=lxs; i<lxe; i++) {
 
                // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
                if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) { 
                    // Forward difference at the domain's min-i boundary
                    dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
                    dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
                    dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
                } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) { 
                    // Backward difference at the domain's max-i boundary
                    dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
                    dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
                    dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
                } else { // Central difference in the interior (or PERIODIC)
                    dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
                    dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
                    dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
                }
 
                // --- Stencil Logic for d/deta (derivative in j-direction) ---
                if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) { 
                    // Forward difference
                    dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
                    dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
                    dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
                } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) { 
                    // Backward difference
                    dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
                    dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
                    dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
                } else { // Central difference (interior or PERIODIC)
                    dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
                    dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
                    dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
                }
 
                // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
                if (k == 0 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) { 
                    // Forward difference
                    dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
                    dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
                    dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
                } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) { 
                    // Backward difference
                    dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
                    dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
                    dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
                } else { // Central difference (Interior or PERIODIC)
                    dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
                    dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
                    dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
                }
 
                // --- Metric calculations (identical to legacy FormMetrics) ---
                kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
                kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
                kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
 
                keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
                keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
                keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
 
                kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
                kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
                kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
 
                kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
                if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
                    kaj[k][j][i] = 1.0 / kaj[k][j][i];
                }
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentz, &centz_const); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
    
    // --- Part 3: Assemble global vectors and update local ghosts ---
    ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
 
    ierr = UpdateLocalGhosts(user, FIELD_ID_KCSI); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_KETA); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_KZET); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_KAJ); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMetricsDivergence()

PetscErrorCode ComputeMetricsDivergence

(

UserCtx *

user

)

Performs a diagnostic check on the divergence of the face area metric vectors.

For a closed cell, the sum of the face area vectors should be zero (Gauss's divergence theorem). This function computes a measure of this divergence and reports the maximum value over the domain. A small value indicates a well-formed grid. This is a direct adaptation of the legacy function.

Parameters

user	The UserCtx for a specific grid level (typically the finest).

Returns

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

Performs a diagnostic check on the divergence of the face area metric vectors.

Local to this translation unit.

Definition at line 1769 of file Metric.c.

{
  DM            da = user->da, fda = user->fda;
  DMDALocalInfo info = user->info;
  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;
  PetscInt      lxs, lys, lzs, lxe, lye, lze;
  PetscInt      i, j, k;
  Vec           Div;
  PetscReal     ***div, ***aj;
  Cmpnts        ***csi, ***eta, ***zet;
  PetscReal     maxdiv;
 
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs == 0) lxs = xs + 1;
  if (ys == 0) lys = ys + 1;
  if (zs == 0) lzs = zs + 1;
 
  if (xe == mx) lxe = xe - 1;
  if (ye == my) lye = ye - 1;
  if (ze == mz) lze = ze - 1;
 
  DMDAVecGetArray(fda, user->lCsi, &csi);
  DMDAVecGetArray(fda, user->lEta, &eta);
  DMDAVecGetArray(fda, user->lZet, &zet);
  DMDAVecGetArray(da, user->lAj, &aj);
 
  VecDuplicate(user->P, &Div);
  VecSet(Div, 0.);
  DMDAVecGetArray(da, Div, &div);
 
  for (k = lzs; k < lze; k++) {
    for (j = lys; j < lye; j++) {
      for (i = lxs; i < lxe; i++) {
        PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
                                eta[k][j][i].x - eta[k][j-1][i].x +
                                zet[k][j][i].x - zet[k-1][j][i].x +
                                csi[k][j][i].y - csi[k][j][i-1].y +
                                eta[k][j][i].y - eta[k][j-1][i].y +
                                zet[k][j][i].y - zet[k-1][j][i].y +
                                csi[k][j][i].z - csi[k][j][i-1].z +
                                eta[k][j][i].z - eta[k][j-1][i].z +
                                zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
        div[k][j][i] = fabs(divergence);
      }
    }
  }
 
  DMDAVecRestoreArray(da, Div, &div);
 
  PetscInt MaxFlatIndex = -1;
  VecMax(Div, &MaxFlatIndex, &maxdiv);
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %" PetscInt_FMT ".\n",maxdiv,MaxFlatIndex);
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    if (Gidx(i,j,k,user) == MaxFlatIndex) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,(int)k,(int)j,(int)i);
    }
      }
    }
  }
 
  
  DMDAVecRestoreArray(fda, user->lCsi, &csi);
  DMDAVecRestoreArray(fda, user->lEta, &eta);
  DMDAVecRestoreArray(fda, user->lZet, &zet);
  DMDAVecRestoreArray(da, user->lAj, &aj);
  VecDestroy(&Div);
 
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMetricNorms()

PetscErrorCode ComputeMetricNorms

(

UserCtx *

user

)

Computes the max-min values of the grid metrics.

This function serves as a diagnostic tool to assess the quality of the grid metrics. It calculates the bounds of the face metrics (Csi, Eta, Zet).

Parameters

user	The UserCtx, containing all necessary grid data.

Returns

PetscErrorCode

Computes the max-min values of the grid metrics.

Local to this translation unit.

Definition at line 1861 of file Metric.c.

{
  
  DMDALocalInfo info = user->info;
  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      i, j, k;
  
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  PetscReal CsiMax, EtaMax, ZetMax;
  PetscReal ICsiMax, IEtaMax, IZetMax;
  PetscReal JCsiMax, JEtaMax, JZetMax;
  PetscReal KCsiMax, KEtaMax, KZetMax;
  PetscReal AjMax, IAjMax, JAjMax, KAjMax;
 
  PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
  PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
  PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
  PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
  PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;  
 
  // Max Values
  VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
  VecMax(user->lEta,&EtaMaxArg,&EtaMax);
  VecMax(user->lZet,&ZetMaxArg,&ZetMax);
 
  VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
  VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
  VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
 
  VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
  VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
  VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
 
  VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
  VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
  VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
 
  VecMax(user->lAj,&AjMaxArg,&AjMax);
  VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
  VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
  VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
 
  VecMax(user->lAj,&AjMaxArg,&AjMax);
  VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
  VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
  VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
 
  LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
  
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    if (Gidx(i,j,k,user) == CsiMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == EtaMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == ZetMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == ICsiMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == IEtaMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == IZetMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == JCsiMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == JEtaMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == JZetMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == KCsiMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == KEtaMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == KZetMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == AjMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == IAjMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == JAjMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
    }
    if (Gidx(i,j,k,user) == KAjMaxArg) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
    }
      }
    }
  }
 
  /*
  VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
  VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
  VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
  */
  
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateAllGridMetrics()

PetscErrorCode CalculateAllGridMetrics

(

SimCtx *

simCtx

)

Orchestrates the calculation of all grid metrics.

This function iterates through every UserCtx in the multigrid and multi-block hierarchy. For each context, it calls a series of modern, modular helper functions to compute the face metrics (Csi, Eta, Zet), the cell-centered inverse Jacobian (Aj), and to validate the grid's orientation.

Parameters

simCtx

The master SimCtx, containing the configured UserCtx hierarchy.

Returns

PetscErrorCode

Orchestrates the calculation of all grid metrics.

Local to this translation unit.

Definition at line 1993 of file Metric.c.

{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx = usermg->mgctx;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
 
    // Loop through all levels and all blocks
    for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *user = &mgctx[level].user[bi];
            LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
 
            // Call the modern, modular helper functions for each UserCtx.
            // These functions are self-contained and operate on the data within the provided context.
            ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
            ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
            ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
            ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
            ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
            ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
            ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr); 
 
            // Apply Periodic Boundary Condition Adjustments if necessary
            ierr = ApplyMetricsPeriodicBCs(user); CHKERRQ(ierr);
            // Diagnostics
        ierr = ComputeMetricNorms(user);
            if (level == usermg->mglevels - 1) {
                ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
            }
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: