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 wrapper: map (cell index, ξ,η,ζ) to (x,y,z).
PetscErrorCode	MetricGetCellVertices (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
	Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
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)
	Compute Jacobian matrix and its determinant at (xi,eta,zta).
PetscErrorCode	MetricVelocityContravariant (const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
	Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).
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)
	Applies periodic boundary corrections to cell centers (Cent) and grid spacing (GridSpace).
PetscErrorCode	ApplyPeriodicCorrectionsToKFaceCenter (UserCtx *user)
	Applies periodic boundary corrections to i-face centers (Centx).
PetscErrorCode	ApplyPeriodicCorrectionsToJFaceCenter (UserCtx *user)
	Applies periodic boundary corrections to j-face centers (Centy).
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 wrapper: map (cell index, ξ,η,ζ) to (x,y,z).

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

Extract the eight vertex coordinates of the hexahedral cell (i,j,k).

Vertices are returned in the standard trilinear ordering: bit 0 → x-corner, bit 1 → y-corner, bit 2 → z-corner. (000 = origin of the cell, 111 = far corner.)

Definition at line 27 of file Metric.c.

{
  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
	)

Compute Jacobian matrix and its determinant at (xi,eta,zta).

   J = [ x_ξ  x_η  x_ζ ]
       [ y_ξ  y_η  y_ζ ]
       [ z_ξ  z_η  z_ζ ]

This is handy for converting physical velocities (u,v,w) into contravariant components and for volume weighting.

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

Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).

Definition at line 167 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);
}

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,eta,zet	The metric vectors at the cell center.
ni	Output: A 3-element array for the unit normal vector of the i-face.
nj	Output: A 3-element array for the unit normal vector of the j-face.
nk	Output: A 3-element array for the unit normal vector of the k-face.
Ai	Output: Pointer to store the area of the i-face.
Aj	Output: Pointer to store the area of the j-face.
Ak	Output: Pointer to store the area of the k-face.

Returns

PetscErrorCode 0 on success.

Definition at line 257 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.

Definition at line 212 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 Jacobian of the grid transformation (1 / cell volume).
csi,eta,zet	The metric vectors at the cell center.
dx	Output: Pointer to store the characteristic length in the x-direction.
dy	Output: Pointer to store the characteristic length in the y-direction.
dz	Output: Pointer to store the characteristic length in the z-direction.

Returns

PetscErrorCode 0 on success.

Definition at line 313 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

)

Applies periodic boundary corrections to cell centers (Cent) and grid spacing (GridSpace).

This function handles the special logic needed when periodic boundaries are present. For coarse grids (cgrid), it directly copies from ghost regions. For fine grids, it calculates the boundary values using the GridSpace information.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

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

ApplyPeriodicCorrectionsToKFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter

(

UserCtx *

user

)

Applies periodic boundary corrections to i-face centers (Centx).

Applies periodic boundary corrections to k-face centers (Centz).

Only X-direction periodicity affects Centx. This function must be called after Centx has been initially computed but before it's used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Only Z-direction periodicity affects Centz. This function must be called after Centz has been initially computed but before it's used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Applies periodic boundary corrections to i-face centers (Centx).

Only Z-direction periodicity affects Centz. This function must be called after Centz has been initially computed but before it's used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Definition at line 828 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    PetscInt       zs = info.zs, ze = info.zs + info.zm;
    PetscInt       mz = info.mz;
    PetscInt       gxs = info.gxs, gxe = info.gxs + info.gxm;
    PetscInt       gys = info.gys, gye = info.gys + info.gym;
    PetscInt       gzs = info.gzs, gze = info.gzs + info.gzm;
    Cmpnts       ***centz, ***gs;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // Check if Z-periodic boundaries exist
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC && 
        user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
        LOG_ALLOW(LOCAL, LOG_TRACE, "No Z-periodic boundaries; skipping Centz corrections.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying Z-periodic corrections to Centz.\n");
 
    ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); 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=gys+1; j<gye; j++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centz[-1][j][i].x = centz[-3][j][i].x;
                    centz[-1][j][i].y = centz[-3][j][i].y;
                    centz[-1][j][i].z = centz[-3][j][i].z;
                }
            }
        } else {
            for (PetscInt j=gys+1; j<gye; j++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centz[-1][j][i].x = centz[0][j][i].x;
                    centz[-1][j][i].y = centz[0][j][i].y;
                    centz[-1][j][i].z = centz[0][j][i].z - gs[-2][j][i].z;
                }
            }
        }
    }
 
    if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
        if (user->cgrid) {
            for (PetscInt j=gys+1; j<gye; j++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centz[mz-1][j][i].x = centz[mz+1][j][i].x;
                    centz[mz-1][j][i].y = centz[mz+1][j][i].y;
                    centz[mz-1][j][i].z = centz[mz+1][j][i].z;
                }
            }
        } else {
            for (PetscInt j=gys+1; j<gye; j++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centz[mz-1][j][i].x = centz[mz-2][j][i].x;
                    centz[mz-1][j][i].y = centz[mz-2][j][i].y;
                    centz[mz-1][j][i].z = centz[mz-2][j][i].z + gs[mz+1][j][i].z;
                }
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ApplyPeriodicCorrectionsToJFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter

(

UserCtx *

user

)

Applies periodic boundary corrections to j-face centers (Centy).

Only Y-direction periodicity affects Centy. This function must be called after Centy has been initially computed but before it's used for metric calculations.

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Definition at line 742 of file Metric.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    PetscInt       ys = info.ys, ye = info.ys + info.ym;
    PetscInt       my = info.my;
    PetscInt       gxs = info.gxs, gxe = info.gxs + info.gxm;
    PetscInt       gys = info.gys, gye = info.gys + info.gym;
    PetscInt       gzs = info.gzs, gze = info.gzs + info.gzm;
    Cmpnts       ***centy, ***gs;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // Check if Y-periodic boundaries exist
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC && 
        user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
        LOG_ALLOW(LOCAL, LOG_TRACE, "No Y-periodic boundaries; skipping Centy corrections.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying Y-periodic corrections to Centy.\n");
 
    ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); 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=gzs+1; k<gze; k++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centy[k][-1][i].x = centy[k][-3][i].x;
                    centy[k][-1][i].y = centy[k][-3][i].y;
                    centy[k][-1][i].z = centy[k][-3][i].z;
                }
            }
        } else {
            for (PetscInt k=gzs+1; k<gze; k++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centy[k][-1][i].x = centy[k][0][i].x;
                    centy[k][-1][i].y = centy[k][0][i].y - gs[k][-2][i].y;
                    centy[k][-1][i].z = centy[k][0][i].z;
                }
            }
        }
    }
 
    if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
        if (user->cgrid) {
            for (PetscInt k=gzs+1; k<gze; k++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centy[k][my-1][i].x = centy[k][my+1][i].x;
                    centy[k][my-1][i].y = centy[k][my+1][i].y;
                    centy[k][my-1][i].z = centy[k][my+1][i].z;
                }
            }
        } else {
            for (PetscInt k=gzs+1; k<gze; k++) {
                for (PetscInt i=gxs+1; i<gxe; i++) {
                    centy[k][my-1][i].x = centy[k][my-2][i].x;
                    centy[k][my-1][i].y = centy[k][my-2][i].y + gs[k][my+1][i].y;
                    centy[k][my-1][i].z = centy[k][my-2][i].z;
                }
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller 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.

  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.

This function calculates the face area vectors, which are fundamental to the finite volume method on a curvilinear grid. The process is performed in two main stages to ensure robustness and correctness across the entire domain, including physical boundaries.

Stage 1: Interior Face Calculation The core of the function uses centered finite-difference stencils on the nodal coordinates (coor) to compute the metric terms. For example, csi at a node (k,j,i) is calculated using a 2x2 stencil of nodes in the J-K plane.

Crucially, these stencils are only valid for interior nodes/faces where all required neighboring nodes exist. The loops are intentionally constructed (e.g., i_loop_start = 1) to skip the calculation directly on the domain's physical boundaries (e.g., at i=0, j=0, etc.), as the stencil would require out-of-bounds data.

Stage 2: Boundary Face Extrapolation After the interior metrics are computed, the values for the faces lying on the physical boundaries of the domain are populated. This is done via zero-order extrapolation, which simply means copying the values from the nearest valid interior face layer.

For example, on a rank owning the global i=0 boundary, zet_arr[k][j][0] is explicitly set to be equal to zet_arr[k][j][1]. Similarly, on a rank owning the i=mx-1 boundary, zet_arr[k][j][mx-1] is set to zet_arr[k][j][mx-2].

This two-stage approach ensures that every node in the physical domain, including the boundaries, has a valid metric value associated with it.

The function concludes by assembling the global PETSc vectors and updating the local ghosted versions, making the computed metrics ready for use by the solver.

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.

• Csi,Eta,Zeta are face centered quantities which represent the surface area in magnitude and the direction of positive computational coordinate increase in direction.

Definition at line 956 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.

Parameters

[in,out]

user

Pointer to the UserCtx structure.

Returns

PetscErrorCode 0 on success.

Definition at line 1171 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

<your name>

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

vishal kandala

Definition at line 368 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.

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.

Note

Grid Spacing represents the vector connecting two adjacent grid cell faces. It is a vector whose magnitude gives the distance between two faces and direction is from one face center to the other.

Definition at line 1326 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.

This function calculates the metric terms (ICsi, IEta, IZet) and the inverse Jacobian (IAj) located at the geometric center of each constant-i 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->Centx vector.
2. It then uses a boundary-aware, second-order finite difference stencil on the Centx 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->ICsi, user->IEta, user->IZet, and user->IAj vectors.

Returns

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

Note

Icsi,Ieta,Izet represent the computational coordinate gradients scaled by local volume inverse ((Iaj)

Definition at line 1450 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 gxs = info.gxs, gxe = info.gxs + info.gxm;
    PetscInt gys = info.gys, gye = info.gys + info.gym;
    PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
    
    PetscInt lxs = xs; PetscInt lxe = xe;
    PetscInt lys = ys; PetscInt lye = ye;
    PetscInt lzs = zs; PetscInt 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;
 
    // --- 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);
 
    // Loop over the ghosted region to calculate all local face centers
    // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
    PetscInt j_end = (ye == mx)? my - 1:gye;
    PetscInt k_end = (ze == mz)? mz - 1:gze;
 
    for (PetscInt k = gzs + 1; k < k_end; k++) {
        for (PetscInt j = gys + 1; j < j_end; j++) {
            for (PetscInt i = gxs; i < gxe; 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);
 
    // For periodic BCs, exchange ghosts between processes
    if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || 
        user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
        ierr = DMLocalToLocalBegin(user->fda, user->Centx, INSERT_VALUES, user->Centx); CHKERRQ(ierr);
        ierr = DMLocalToLocalEnd(user->fda, user->Centx, INSERT_VALUES, user->Centx); CHKERRQ(ierr);
    }
 
    ierr = ApplyPeriodicCorrectionsToIFaceCenter(user); 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->Centx, &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->Centx, &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.

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.

Note

Jcsi,Jeta,Jzet represent the computational coordinate gradients scaled by local volume inverse ((Jaj)

Definition at line 1692 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 gxs = info.gxs, gxe = info.gxs + info.gxm;
    PetscInt gys = info.gys, gye = info.gys + info.gym;
    PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
    
    PetscInt lxs = xs; PetscInt lxe = xe;
    PetscInt lys = ys; PetscInt lye = ye;
    PetscInt lzs = zs; PetscInt 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;
 
    // --- 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);
 
    // Loop over the ghosted region to calculate all local face centers
    // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
    PetscInt k_end = (ze == mz)? mz - 1:gze;
    PetscInt i_end = (xe == mx)? mx - 1:gxe;
 
    for (PetscInt k = gzs + 1; k < k_end; k++) {
        for (PetscInt j = gys; j < gye; j++) {
            for (PetscInt i = gxs + 1; i < i_end; 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);
 
    // For periodic BCs, exchange ghosts between processes
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || 
        user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
        ierr = DMLocalToLocalBegin(user->fda, user->Centy, INSERT_VALUES, user->Centy); CHKERRQ(ierr);
        ierr = DMLocalToLocalEnd(user->fda, user->Centy, INSERT_VALUES, user->Centy); CHKERRQ(ierr);
    }
 
    ierr = ApplyPeriodicCorrectionsToJFaceCenter(user); 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->Centy, &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->Centy, &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.

This function calculates the metric terms (KCsi, KEta, KZet) and the inverse Jacobian (KAj) located at the geometric center of each constant-k 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.

Note

Kcsi,Keta,Kzet represent the computational coordinate gradients scaled by local volume inverse ((Kaj)

Definition at line 1920 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 gxs = info.gxs, gxe = info.gxs + info.gxm;
    PetscInt gys = info.gys, gye = info.gys + info.gym;
    PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
    
    PetscInt lxs = xs; PetscInt lxe = xe;
    PetscInt lys = ys; PetscInt lye = ye;
    PetscInt lzs = zs; PetscInt 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;
 
    // --- 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);
 
    // Loop over the ghosted region to calculate all local face centers
    // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
    PetscInt i_end = (xe == mx)? mx - 1:gxe;
    PetscInt j_end = (ye == my)? my - 1:gye;
    for (PetscInt k = gzs; k < gze; k++) {
        for (PetscInt j = gys; j < gye; j++) {
            for (PetscInt i = gxs + 1; i < gxe; 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);
 
    // For periodic BCs, exchange ghosts between processes
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || 
        user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
        ierr = DMLocalToLocalBegin(user->fda, user->Centz, INSERT_VALUES, user->Centz); CHKERRQ(ierr);
        ierr = DMLocalToLocalEnd(user->fda, user->Centz, INSERT_VALUES, user->Centz); CHKERRQ(ierr);
    }
 
    ierr = ApplyPeriodicCorrectionsToKFaceCenter(user); 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->Centz, &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->Centz, &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.

This function serves as a diagnostic tool to assess the quality of the grid metrics. It calculates the divergence of the face metrics (Csi, Eta, Zet) and scales it by the inverse of the cell Jacobian. The maximum divergence value is then located, and its grid coordinates are printed to the console, helping to pinpoint areas of potential grid quality issues.

Parameters

user	The UserCtx, containing all necessary grid data.

Returns

PetscErrorCode

Definition at line 2152 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);
 
  PetscReal MaxFlatIndex;
  VecMax(Div, &MaxFlatIndex, &maxdiv);
  LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %d.\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

Definition at line 2249 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

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