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)
	Public interface for MetricLogicalToPhysical().
PetscErrorCode	MetricGetCellVertices (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
	Public interface for MetricGetCellVertices().
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)
	Public interface for MetricJacobian().
PetscErrorCode	MetricVelocityContravariant (const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
	Public interface for MetricVelocityContravariant().
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 covariant metric tensor to obtain the contravariant metric tensor.
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	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)
	Ensure a 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
	)

Public interface for MetricLogicalToPhysical().

Parameters

user	Primary UserCtx input for the operation.
X	Parameter X passed to MetricLogicalToPhysical().
i	Parameter i passed to MetricLogicalToPhysical().
j	Parameter j passed to MetricLogicalToPhysical().
k	Parameter k passed to MetricLogicalToPhysical().
xi	Parameter xi passed to MetricLogicalToPhysical().
eta	Parameter eta passed to MetricLogicalToPhysical().
zta	Parameter zta passed to MetricLogicalToPhysical().
Xp	Parameter Xp passed to MetricLogicalToPhysical().

Returns

PetscErrorCode 0 on success.

Public interface for MetricLogicalToPhysical().

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

See also

MetricLogicalToPhysical()

Definition at line 77 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]
	)

Public interface for MetricGetCellVertices().

Parameters

user	Primary UserCtx input for the operation.
X	Parameter X passed to MetricGetCellVertices().
i	Parameter i passed to MetricGetCellVertices().
j	Parameter j passed to MetricGetCellVertices().
k	Parameter k passed to MetricGetCellVertices().
V	Parameter V passed to MetricGetCellVertices().

Returns

PetscErrorCode 0 on success.

Public interface for MetricGetCellVertices().

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
	)

Public interface for MetricJacobian().

Parameters

user	Primary UserCtx input for the operation.
X	Parameter X passed to MetricJacobian().
i	Parameter i passed to MetricJacobian().
j	Parameter j passed to MetricJacobian().
k	Parameter k passed to MetricJacobian().
xi	Parameter xi passed to MetricJacobian().
eta	Parameter eta passed to MetricJacobian().
zta	Parameter zta passed to MetricJacobian().
J	Parameter J passed to MetricJacobian().
detJ	Parameter detJ passed to MetricJacobian().

Returns

PetscErrorCode 0 on success.

Public interface for MetricJacobian().

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

See also

MetricJacobian()

Definition at line 106 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]
	)

Public interface for MetricVelocityContravariant().

Parameters

J	Parameter J passed to MetricVelocityContravariant().
detJ	Parameter detJ passed to MetricVelocityContravariant().
u	Parameter u passed to MetricVelocityContravariant().
uc	Parameter uc passed to MetricVelocityContravariant().

Returns

PetscErrorCode 0 on success.

Public interface for MetricVelocityContravariant().

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

See also

MetricVelocityContravariant()

Definition at line 166 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	Parameter csi passed to CalculateFaceNormalAndArea().
eta	Parameter eta passed to CalculateFaceNormalAndArea().
zet	Parameter zet passed to CalculateFaceNormalAndArea().
ni	Output:
nj	Output:
nk	Output:
Ai	Output:
Aj	Output:
Ak	Output:

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 237 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;
    
    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 covariant metric tensor to obtain the contravariant metric tensor.

In curvilinear coordinates, the input matrix g contains the dot products of the covariant basis vectors (e.g., g_ij = e_i . e_j). Its inverse, G, is the contravariant metric tensor, which is essential for transforming vectors and tensors between coordinate systems.

Parameters

covariantTensor	Input: A 3x3 matrix representing the covariant metric tensor.
contravariantTensor	Output: A 3x3 matrix where the inverted result is stored.

Returns

PetscErrorCode 0 on success.

Inverts the 3x3 covariant metric tensor to obtain the contravariant metric tensor.

Local to this translation unit.

Definition at line 203 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);
 
    if (fabs(det) < 1.0e-12) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Matrix is singular, determinant is near zero.");
    }
 
    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.

Parameters

ajc	The
csi	Parameter csi passed to ComputeCellCharacteristicLengthScale().
eta	Parameter eta passed to ComputeCellCharacteristicLengthScale().
zet	Parameter zet passed to ComputeCellCharacteristicLengthScale().
dx	Output:
dy	Output:
dz	Output:

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 283 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;
 
        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:

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 392 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, "Cent"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "GridSpace"); 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 591 of file Metric.c.

{
    const char *fields[] = {"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 607 of file Metric.c.

{
    const char *fields[] = {"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 623 of file Metric.c.

{
    const char *fields[] = {"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 638 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, "Csi"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "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 850 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, "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

)

Ensure a 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) → flip Csi, Eta, Zet, Aj & update local ghosts.
4	Store user->orientation = ±1 so BC / IC routines can apply sign-aware logic if they care about inlet direction.

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.

Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain.

Local to this translation unit.

Definition at line 314 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.");
 
    /* Default: grid is right-handed unless proven otherwise */
    PetscInt orientation = +1;
 
    /* ---------------- step 3: repair left-handed mesh -------------- */
    if (aj_max < 0.0) {                  /* entire domain has Aj < 0   */
        orientation = -1;
 
        if (!rank)
            LOG_ALLOW(LOCAL, LOG_INFO,
                "[orientation] Detected left-handed grid – flipping metric vectors\n");
 
        /* Flip sign of *all* metric vectors and Aj                     */
        ierr = VecScale(user->Csi, -1.0); CHKERRQ(ierr);
        ierr = VecScale(user->Eta, -1.0); CHKERRQ(ierr);
        ierr = VecScale(user->Zet, -1.0); CHKERRQ(ierr);
        ierr = VecScale(user->Aj , -1.0); CHKERRQ(ierr);
 
        /* Local ghost regions now stale – refresh                      */
        ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, "Aj");  CHKERRQ(ierr);
 
        /* Sanity print: Aj must be > 0 now                             */
        ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
        ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
 
        if (aj_min <= 0.0)
            SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
                "Failed to flip grid orientation – Aj still non-positive.");
    else if (aj_min && aj_max > 0.0)
      orientation = +1;
    }
 
    /* ---------------- step 4: store result in UserCtx -------------- */
    user->GridOrientation = orientation;
 
    if (!rank)
        LOG_ALLOW(LOCAL, LOG_INFO,
            "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
            (orientation>0) ? "right" : "left", orientation);
 
    PROFILE_FUNCTION_END;        
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

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 994 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, "Cent"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "GridSpace"); 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, "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 1097 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 char *face_centers[] = {"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, "ICsi"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "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 1308 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 char *face_centers[] = {"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, "JCsi"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "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 1503 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 char *face_centers[] = {"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, "KCsi"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "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 1718 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 1810 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 1942 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: