Metric.h File Reference

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

Include dependency graph for Metric.h:

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

Go to the source code of this file.

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

Function Documentation

MetricLogicalToPhysical()

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

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

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

Parameters

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

Returns

PetscErrorCode 0 on success.

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

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

See also

MetricLogicalToPhysical()

Definition at line 76 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

MetricGetCellVertices()

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

Collects the eight node coordinates of one logical hexahedral cell.

Vertex ordering matches the trilinear metric kernels in this module.

Parameters

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

Returns

PetscErrorCode 0 on success.

Collects the eight node coordinates of one logical hexahedral cell.

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

See also

MetricGetCellVertices()

Definition at line 26 of file Metric.c.

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

Here is the caller graph for this function:

MetricJacobian()

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

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

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

Parameters

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

Returns

PetscErrorCode 0 on success.

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

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

See also

MetricJacobian()

Definition at line 105 of file Metric.c.

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

Here is the call graph for this function:

MetricVelocityContravariant()

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

Converts a Cartesian velocity vector to contravariant logical components.

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

Parameters

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

Returns

PetscErrorCode 0 on success.

Converts a Cartesian velocity vector to contravariant logical components.

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

See also

MetricVelocityContravariant()

Definition at line 165 of file Metric.c.

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

Here is the caller graph for this function:

CalculateFaceNormalAndArea()

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

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

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

Parameters

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

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 236 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.
[out]	contravariantTensor	Inverse metric tensor written in place.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 202 of file Metric.c.

{
    PetscFunctionBeginUser;
 
    const double a11=covariantTensor[0][0], a12=covariantTensor[0][1], a13=covariantTensor[0][2];
    const double a21=covariantTensor[1][0], a22=covariantTensor[1][1], a23=covariantTensor[1][2];
    const double a31=covariantTensor[2][0], a32=covariantTensor[2][1], a33=covariantTensor[2][2];
 
    double det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
 
    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	Cell Jacobian/volume metric at the cell center.
	csi	Covariant xi metric vector.
	eta	Covariant eta metric vector.
	zet	Covariant zeta metric vector.
[out]	dx	Effective physical length associated with xi variation.
[out]	dy	Effective physical length associated with eta variation.
[out]	dz	Effective physical length associated with zeta variation.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 282 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyPeriodicCorrectionsToCellCentersAndSpacing()

PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing

(

UserCtx *

user

)

Builds translated periodic images for cell centers and grid spacing.

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

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

Builds translated periodic images for cell centers and grid spacing.

Local to this translation unit.

Definition at line 391 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyPeriodicCorrectionsToIFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter

(

UserCtx *

user

)

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

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

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 590 of file Metric.c.

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

Here is the call graph for this function:

ApplyPeriodicCorrectionsToJFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter

(

UserCtx *

user

)

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

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

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 606 of file Metric.c.

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

Here is the call graph for this function:

ApplyPeriodicCorrectionsToKFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter

(

UserCtx *

user

)

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

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

Parameters

user	The UserCtx containing grid and field data.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 622 of file Metric.c.

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

Here is the call graph for this function:

ComputeFaceMetrics()

PetscErrorCode ComputeFaceMetrics

(

UserCtx *

user

)

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

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

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

Parameters

[in,out]

user

Pointer to the UserCtx structure.

Returns

PetscErrorCode 0 on success.

Note

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

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

Local to this translation unit.

Definition at line 637 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeCellCenteredJacobianInverse()

PetscErrorCode ComputeCellCenteredJacobianInverse

(

UserCtx *

user

)

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

This version includes boundary extrapolation.

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

Parameters

[in,out]

user

Pointer to the UserCtx structure.

Returns

PetscErrorCode 0 on success.

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

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

See also

ComputeCellCenteredJacobianInverse()

Definition at line 849 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

CheckAndFixGridOrientation()

PetscErrorCode CheckAndFixGridOrientation

(

UserCtx *

user

)

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

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

Step	Action
1	Compute global Aj_min, Aj_max.
2	Mixed signs (Aj_min < 0 && Aj_max > 0) → abort: the mesh is topologically inconsistent.
3	All negative (Aj_max < 0) → flip Csi, Eta, Zet, Aj & update local ghosts.
4	Store user->orientation = ±1 so BC / IC routines can apply sign-aware logic if they care about inlet direction.

Parameters

[in,out]

user

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

Returns

0 on success or a PETSc error code on failure.

Note

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

Author metadata intentionally omitted in API docs.

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

Local to this translation unit.

Definition at line 313 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, FIELD_ID_CSI); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_ETA); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_ZET); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_AJ);  CHKERRQ(ierr);
 
        /* Sanity print: Aj must be > 0 now                             */
        ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
        ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
 
        if (aj_min <= 0.0)
            SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
                "Failed to flip grid orientation – Aj still non-positive.");
    else if (aj_min && aj_max > 0.0)
      orientation = +1;
    }
 
    /* ---------------- step 4: store result in UserCtx -------------- */
    user->GridOrientation = orientation;
 
    if (!rank)
        LOG_ALLOW(LOCAL, LOG_INFO,
            "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
            (orientation>0) ? "right" : "left", orientation);
 
    PROFILE_FUNCTION_END;        
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeCellCentersAndSpacing()

PetscErrorCode ComputeCellCentersAndSpacing

(

UserCtx *

user

)

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

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

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

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

Parameters

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

Returns

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

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

Local to this translation unit.

Definition at line 993 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeIFaceMetrics()

PetscErrorCode ComputeIFaceMetrics

(

UserCtx *

user

)

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

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

Parameters

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

Returns

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

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

Local to this translation unit.

Definition at line 1096 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeJFaceMetrics()

PetscErrorCode ComputeJFaceMetrics

(

UserCtx *

user

)

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

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

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

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

Parameters

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

Returns

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

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

Local to this translation unit.

Definition at line 1307 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeKFaceMetrics()

PetscErrorCode ComputeKFaceMetrics

(

UserCtx *

user

)

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

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

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

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

Parameters

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

Returns

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

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

Local to this translation unit.

Definition at line 1502 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMetricsDivergence()

PetscErrorCode ComputeMetricsDivergence

(

UserCtx *

user

)

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

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

Parameters

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

Returns

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

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

Local to this translation unit.

Definition at line 1716 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMetricNorms()

PetscErrorCode ComputeMetricNorms

(

UserCtx *

user

)

Computes the max-min values of the grid metrics.

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

Parameters

user	The UserCtx, containing all necessary grid data.

Returns

PetscErrorCode

Computes the max-min values of the grid metrics.

Local to this translation unit.

Definition at line 1808 of file Metric.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateAllGridMetrics()

PetscErrorCode CalculateAllGridMetrics

(

SimCtx *

simCtx

)

Orchestrates the calculation of all grid metrics.

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

Parameters

simCtx

The master SimCtx, containing the configured UserCtx hierarchy.

Returns

PetscErrorCode

Orchestrates the calculation of all grid metrics.

Local to this translation unit.

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