PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
Metric.c File Reference
#include <petsc.h>
#include "Metric.h"
Include dependency graph for Metric.c:

Go to the source code of this file.

Macros

#define __FUNCT__   "MetricGetCellVertices"
 
#define __FUNCT__   "TrilinearBlend"
 
#define __FUNCT__   "MetricLogicalToPhysical"
 
#define __FUNCT__   "MetricJacobian"
 
#define __FUNCT__   "MetricVelocityContravariant"
 
#define __FUNCT__   "InvertCovariantMetricTensor"
 
#define __FUNCT__   "CalculateFaceNormalAndArea"
 
#define __FUNCT__   "ComputeCellCharacteristicLengthScale"
 
#define __FUNCT__   "CheckAndFixGridOrientation"
 
#define __FUNCT__   "ComputeFaceMetrics"
 
#define __FUNCT__   "ComputeCellCenteredJacobianInverse"
 
#define __FUNCT__   "ComputeCellCentersAndSpacing"
 
#define __FUNCT__   "ComputeIFaceMetrics"
 
#define __FUNCT__   "ComputeJFaceMetrics"
 
#define __FUNCT__   "ComputeJFaceMetrics"
 
#define __FUNCT__   "ComputeMetricsDivergence"
 
#define __FUNCT__   "ComputeMetricNorms"
 
#define __FUNCT__   "ComputeAllGridMetrics"
 

Functions

PetscErrorCode MetricGetCellVertices (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
 Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
 
static void TrilinearBlend (const Cmpnts V[8], PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
 Map logical coordinates to physical space using trilinear basis.
 
PetscErrorCode MetricLogicalToPhysical (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
 Public wrapper: map (cell index, ξ,η,ζ) to (x,y,z).
 
PetscErrorCode MetricJacobian (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, PetscReal J[3][3], PetscReal *detJ)
 Compute Jacobian matrix and its determinant at (xi,eta,zta).
 
PetscErrorCode MetricVelocityContravariant (const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
 Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).
 
PetscErrorCode InvertCovariantMetricTensor (double covariantTensor[3][3], double contravariantTensor[3][3])
 Inverts the 3x3 covariant metric tensor to obtain the contravariant metric tensor.
 
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 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 CheckAndFixGridOrientation (UserCtx *user)
 Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain.
 
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), including boundary extrapolation, stores it in user->Aj, assembles user->Aj, and updates user->lAj.
 
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).
 
static PetscInt Gidx (PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
 
PetscErrorCode ComputeMetricsDivergence (UserCtx *user)
 Computes the divergence of the grid metrics and identifies the maximum value.
 
PetscErrorCode ComputeMetricNorms (UserCtx *user)
 Computes the max-min values of the grid metrics.
 
PetscErrorCode CalculateAllGridMetrics (SimCtx *simCtx)
 Orchestrates the calculation of all grid metrics.
 

Macro Definition Documentation

◆ __FUNCT__ [1/18]

#define __FUNCT__   "MetricGetCellVertices"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [2/18]

#define __FUNCT__   "TrilinearBlend"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [3/18]

#define __FUNCT__   "MetricLogicalToPhysical"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [4/18]

#define __FUNCT__   "MetricJacobian"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [5/18]

#define __FUNCT__   "MetricVelocityContravariant"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [6/18]

#define __FUNCT__   "InvertCovariantMetricTensor"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [7/18]

#define __FUNCT__   "CalculateFaceNormalAndArea"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [8/18]

#define __FUNCT__   "ComputeCellCharacteristicLengthScale"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [9/18]

#define __FUNCT__   "CheckAndFixGridOrientation"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [10/18]

#define __FUNCT__   "ComputeFaceMetrics"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [11/18]

#define __FUNCT__   "ComputeCellCenteredJacobianInverse"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [12/18]

#define __FUNCT__   "ComputeCellCentersAndSpacing"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [13/18]

#define __FUNCT__   "ComputeIFaceMetrics"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [14/18]

#define __FUNCT__   "ComputeJFaceMetrics"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [15/18]

#define __FUNCT__   "ComputeJFaceMetrics"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [16/18]

#define __FUNCT__   "ComputeMetricsDivergence"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [17/18]

#define __FUNCT__   "ComputeMetricNorms"

Definition at line 18 of file Metric.c.

◆ __FUNCT__ [18/18]

#define __FUNCT__   "ComputeAllGridMetrics"

Definition at line 18 of file Metric.c.

Function Documentation

◆ MetricGetCellVertices()

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

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

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

Definition at line 27 of file Metric.c.

31{
32 PetscFunctionBeginUser;
33 for (PetscInt c = 0; c < 8; ++c) {
34 PetscInt ii = i + ((c & 1) ? 1 : 0);
35 PetscInt jj = j + ((c & 2) ? 1 : 0);
36 PetscInt kk = k + ((c & 4) ? 1 : 0);
37 LOG_LOOP_ALLOW(GLOBAL, LOG_VERBOSE,i+j+k,10," ii: %d,jj:%d,kk:%d - Retrieved.\n",ii,jj,kk);
38 V[c] = X[kk][jj][ii];
39 }
40 PetscFunctionReturn(0);
41}
LOG_LOOP_ALLOW
#define LOG_LOOP_ALLOW(scope, level, iterVar, interval, fmt,...)
Logs a message inside a loop, but only every interval iterations.
Definition logging.h:313
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
Here is the caller graph for this function:

◆ TrilinearBlend()

static void TrilinearBlend ( const Cmpnts  V[8],
PetscReal  xi,
PetscReal  eta,
PetscReal  zta,
Cmpnts *  Xp 
)
inlinestatic

Map logical coordinates to physical space using trilinear basis.

Parameters
[in]VArray of the eight vertex coordinates (MetricGetCellVertices).
[in]xi,eta,ztaLogical coordinates in [0,1].
[out]XpPhysical coordinate.

Definition at line 54 of file Metric.c.

57{
58 PetscReal x=0,y=0,z=0;
59 for (PetscInt c=0;c<8;++c) {
60 PetscReal N = ((c&1)?xi : 1.0-xi ) *
61 ((c&2)?eta: 1.0-eta) *
62 ((c&4)?zta: 1.0-zta);
63 x += N * V[c].x;
64 y += N * V[c].y;
65 z += N * V[c].z;
66 }
67 Xp->x = x; Xp->y = y; Xp->z = z;
68}
Cmpnts::x
PetscScalar x
Definition variables.h:101
Cmpnts::z
PetscScalar z
Definition variables.h:101
Cmpnts::y
PetscScalar y
Definition variables.h:101
Here is the caller graph for this function:

◆ MetricLogicalToPhysical()

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

Public wrapper: map (cell index, ξ,η,ζ) to (x,y,z).

Definition at line 77 of file Metric.c.

82{
83 PetscErrorCode ierr;
84 Cmpnts V[8];
85 PetscFunctionBeginUser;
86
87 PROFILE_FUNCTION_BEGIN;
88
89 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
90 TrilinearBlend(V,xi,eta,zta,Xp);
91
92 PROFILE_FUNCTION_END;
93
94 PetscFunctionReturn(0);
95}
TrilinearBlend
static void TrilinearBlend(const Cmpnts V[8], PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Map logical coordinates to physical space using trilinear basis.
Definition Metric.c:54
MetricGetCellVertices
PetscErrorCode MetricGetCellVertices(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
Definition Metric.c:27
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
Here is the call graph for this function:
Here is the caller graph for this function:

◆ MetricJacobian()

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

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

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

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

Definition at line 110 of file Metric.c.

115{
116 PetscErrorCode ierr;
117 Cmpnts V[8];
118 PetscFunctionBeginUser;
119
120 PROFILE_FUNCTION_BEGIN;
121
122 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
123
124 /* derivatives of trilinear shape functions */
125 PetscReal dN_dXi[8], dN_dEta[8], dN_dZta[8];
126 for (PetscInt c=0;c<8;++c) {
127 PetscReal sx = (c & 1) ? 1.0 : -1.0;
128 PetscReal sy = (c & 2) ? 1.0 : -1.0;
129 PetscReal sz = (c & 4) ? 1.0 : -1.0;
130 dN_dXi [c] = 0.125 * sx * ( (c&2?eta:1-eta) ) * ( (c&4?zta:1-zta) );
131 dN_dEta[c] = 0.125 * sy * ( (c&1?xi :1-xi ) ) * ( (c&4?zta:1-zta) );
132 dN_dZta[c] = 0.125 * sz * ( (c&1?xi :1-xi ) ) * ( (c&2?eta:1-eta) );
133 }
134
135 /* assemble Jacobian */
136 PetscReal x_xi=0,y_xi=0,z_xi=0,
137 x_eta=0,y_eta=0,z_eta=0,
138 x_zta=0,y_zta=0,z_zta=0;
139 for (PetscInt c=0;c<8;++c) {
140 x_xi += dN_dXi [c]*V[c].x; y_xi += dN_dXi [c]*V[c].y; z_xi += dN_dXi [c]*V[c].z;
141 x_eta += dN_dEta[c]*V[c].x; y_eta += dN_dEta[c]*V[c].y; z_eta += dN_dEta[c]*V[c].z;
142 x_zta += dN_dZta[c]*V[c].x; y_zta += dN_dZta[c]*V[c].y; z_zta += dN_dZta[c]*V[c].z;
143 }
144
145 J[0][0]=x_xi; J[0][1]=x_eta; J[0][2]=x_zta;
146 J[1][0]=y_xi; J[1][1]=y_eta; J[1][2]=y_zta;
147 J[2][0]=z_xi; J[2][1]=z_eta; J[2][2]=z_zta;
148
149 if (detJ) {
150 *detJ = x_xi*(y_eta*z_zta - y_zta*z_eta)
151 - x_eta*(y_xi*z_zta - y_zta*z_xi)
152 + x_zta*(y_xi*z_eta - y_eta*z_xi);
153 }
154
155 PROFILE_FUNCTION_END;
156
157 PetscFunctionReturn(0);
158}
Here is the call graph for this function:

◆ MetricVelocityContravariant()

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

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

Definition at line 167 of file Metric.c.

169{
170 PetscFunctionBeginUser;
171
172 PROFILE_FUNCTION_BEGIN;
173
174 /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
175 PetscReal gxi[3] = { J[1][1]*J[2][2]-J[1][2]*J[2][1],
176 -J[0][1]*J[2][2]+J[0][2]*J[2][1],
177 J[0][1]*J[1][2]-J[0][2]*J[1][1] };
178 PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
179 J[0][0]*J[2][2]-J[0][2]*J[2][0],
180 -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
181 PetscReal gzta[3] = { J[1][0]*J[2][1]-J[1][1]*J[2][0],
182 -J[0][0]*J[2][1]+J[0][1]*J[2][0],
183 J[0][0]*J[1][1]-J[0][1]*J[1][0] };
184
185 PetscReal invDet = 1.0 / detJ;
186 for (int d=0; d<3; ++d) { gxi[d] *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
187
188 uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
189 uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
190 uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
191
192 PROFILE_FUNCTION_END;
193
194 PetscFunctionReturn(0);
195}

◆ 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
covariantTensorInput: A 3x3 matrix representing the covariant metric tensor.
contravariantTensorOutput: A 3x3 matrix where the inverted result is stored.
Returns
PetscErrorCode 0 on success.

Definition at line 212 of file Metric.c.

213{
214 PetscFunctionBeginUser;
215
216 const double a11=covariantTensor[0][0], a12=covariantTensor[0][1], a13=covariantTensor[0][2];
217 const double a21=covariantTensor[1][0], a22=covariantTensor[1][1], a23=covariantTensor[1][2];
218 const double a31=covariantTensor[2][0], a32=covariantTensor[2][1], a33=covariantTensor[2][2];
219
220 double det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
221
222 if (fabs(det) < 1.0e-12) {
223 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Matrix is singular, determinant is near zero.");
224 }
225
226 contravariantTensor[0][0] = (a33*a22-a32*a23)/det;
227 contravariantTensor[0][1] = -(a33*a12-a32*a13)/det;
228 contravariantTensor[0][2] = (a23*a12-a22*a13)/det;
229 contravariantTensor[1][0] = -(a33*a21-a31*a23)/det;
230 contravariantTensor[1][1] = (a33*a11-a31*a13)/det;
231 contravariantTensor[1][2] = -(a23*a11-a21*a13)/det;
232 contravariantTensor[2][0] = (a32*a21-a31*a22)/det;
233 contravariantTensor[2][1] = -(a32*a11-a31*a12)/det;
234 contravariantTensor[2][2] = (a22*a11-a21*a12)/det;
235
236 PetscFunctionReturn(0);
237}
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,eta,zetThe metric vectors at the cell center.
niOutput: A 3-element array for the unit normal vector of the i-face.
njOutput: A 3-element array for the unit normal vector of the j-face.
nkOutput: A 3-element array for the unit normal vector of the k-face.
AiOutput: Pointer to store the area of the i-face.
AjOutput: Pointer to store the area of the j-face.
AkOutput: Pointer to store the area of the k-face.
Returns
PetscErrorCode 0 on success.

Definition at line 257 of file Metric.c.

258{
259 PetscFunctionBeginUser;
260 PROFILE_FUNCTION_BEGIN;
261 double g[3][3];
262 double G[3][3];
263
264 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
265 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
266 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
267
268 InvertCovariantMetricTensor(g, G);
269 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
270 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
271 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
272
273 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
274 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
275 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
276
277 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
278 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
279 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
280
281 *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
282 *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
283 *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
284
285 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
286 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
287 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
288
289 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
290 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
291 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
292
293 PROFILE_FUNCTION_END;
294 PetscFunctionReturn(0);
295}
InvertCovariantMetricTensor
PetscErrorCode InvertCovariantMetricTensor(double covariantTensor[3][3], double contravariantTensor[3][3])
Inverts the 3x3 covariant metric tensor to obtain the contravariant metric tensor.
Definition Metric.c:212
Here is the call graph for this function:
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
ajcThe Jacobian of the grid transformation (1 / cell volume).
csi,eta,zetThe metric vectors at the cell center.
dxOutput: Pointer to store the characteristic length in the x-direction.
dyOutput: Pointer to store the characteristic length in the y-direction.
dzOutput: Pointer to store the characteristic length in the z-direction.
Returns
PetscErrorCode 0 on success.

Definition at line 313 of file Metric.c.

314{
315 PetscFunctionBeginUser;
316 PROFILE_FUNCTION_BEGIN;
317 double ni[3], nj[3], nk[3];
318 double Li, Lj, Lk;
319 double Ai, Aj, Ak;
320 double vol = 1./ajc;
321
322 CalculateFaceNormalAndArea(csi, eta, zet, ni, nj, nk, &Ai, &Aj, &Ak);
323 Li = vol / Ai;
324 Lj = vol / Aj;
325 Lk = vol / Ak;
326
327 // Length scale vector = di * ni_vector + dj * nj_vector + dk * nk_vector
328 *dx = fabs( Li * ni[0] + Lj * nj[0] + Lk * nk[0] );
329 *dy = fabs( Li * ni[1] + Lj * nj[1] + Lk * nk[1] );
330 *dz = fabs( Li * ni[2] + Lj * nj[2] + Lk * nk[2] );
331
332 PROFILE_FUNCTION_END;
333 PetscFunctionReturn(0);
334}
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.
Definition Metric.c:257
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CheckAndFixGridOrientation()

PetscErrorCode CheckAndFixGridOrientation ( UserCtx *  user)

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

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

Step Action
1 Compute global Aj_min, Aj_max.
2 Mixed signs (Aj_min < 0 && Aj_max > 0) → abort: the mesh is topologically inconsistent.
3 All negative (Aj_max < 0) → flip
Csi, Eta, Zet, Aj & update local ghosts.
4 Store user->orientation = ±1 so BC / IC routines can apply sign-aware logic if they care about inlet direction.
Parameters
[in,out]userFully initialised UserCtx that already contains
Csi, Eta, Zet, Aj, their local ghosts, and valid distributed DMs.
Returns
0 on success or a PETSc error code on failure.
Note
Call immediately after ComputeCellCenteredJacobianInverse() and before any routine that differentiates or applies BCs.
Author
vishal kandala

Definition at line 368 of file Metric.c.

369{
370 PetscErrorCode ierr;
371 PetscReal aj_min, aj_max;
372 PetscMPIInt rank;
373
374 PetscFunctionBeginUser;
375
376 PROFILE_FUNCTION_BEGIN;
377
378 /* ---------------- step 1: global extrema of Aj ---------------- */
379 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
380
381 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr); /* already global */
382 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
383
384 LOG_ALLOW(GLOBAL, LOG_INFO,
385 "[orientation] Global Aj range: [%.3e , %.3e]\n",
386 (double)aj_min, (double)aj_max);
387
388 /* ---------------- step 2: detect malformed mesh ---------------- */
389 if (aj_min < 0.0 && aj_max > 0.0)
390 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
391 "Mixed Jacobian signs detected – grid is topologically inconsistent.");
392
393 /* Default: grid is right-handed unless proven otherwise */
394 PetscInt orientation = +1;
395
396 /* ---------------- step 3: repair left-handed mesh -------------- */
397 if (aj_max < 0.0) { /* entire domain has Aj < 0 */
398 orientation = -1;
399
400 if (!rank)
401 LOG_ALLOW(LOCAL, LOG_INFO,
402 "[orientation] Detected left-handed grid – flipping metric vectors\n");
403
404 /* Flip sign of *all* metric vectors and Aj */
405 ierr = VecScale(user->Csi, -1.0); CHKERRQ(ierr);
406 ierr = VecScale(user->Eta, -1.0); CHKERRQ(ierr);
407 ierr = VecScale(user->Zet, -1.0); CHKERRQ(ierr);
408 ierr = VecScale(user->Aj , -1.0); CHKERRQ(ierr);
409
410 /* Local ghost regions now stale – refresh */
411 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
412 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
413 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
414 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
415
416 /* Sanity print: Aj must be > 0 now */
417 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
418 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
419
420 if (aj_min <= 0.0)
421 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
422 "Failed to flip grid orientation – Aj still non-positive.");
423 else if (aj_min && aj_max > 0.0)
424 orientation = +1;
425 }
426
427 /* ---------------- step 4: store result in UserCtx -------------- */
428 user->GridOrientation = orientation;
429
430 if (!rank)
431 LOG_ALLOW(LOCAL, LOG_INFO,
432 "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
433 (orientation>0) ? "right" : "left", orientation);
434
435 PROFILE_FUNCTION_END;
436
437 PetscFunctionReturn(0);
438}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:201
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1091
UserCtx::Zet
Vec Zet
Definition variables.h:705
UserCtx::Csi
Vec Csi
Definition variables.h:705
UserCtx::Eta
Vec Eta
Definition variables.h:705
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:674
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeFaceMetrics()

PetscErrorCode ComputeFaceMetrics ( UserCtx *  user)

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

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

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

Definition at line 464 of file Metric.c.

465{
466 PetscErrorCode ierr;
467 DMDALocalInfo info;
468 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
469 Cmpnts ***nodal_coords_arr;
470 Vec localCoords_from_dm;
471
472 PetscFunctionBeginUser;
473
474 PROFILE_FUNCTION_BEGIN;
475
476 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
477
478 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
479
480 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
481 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
482 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
483 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
484
485 // --- 2. Get arrays for output global Vecs from UserCtx ---
486 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
487 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
488 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
489
490 // Define owned node ranges (global indices)
491 PetscInt xs = info.xs, xe = info.xs + info.xm;
492 PetscInt ys = info.ys, ye = info.ys + info.ym;
493 PetscInt zs = info.zs, ze = info.zs + info.zm;
494
495 // Global domain dimensions (total number of nodes)
496 PetscInt mx = info.mx;
497 PetscInt my = info.my;
498 PetscInt mz = info.mz;
499
500 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
501 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
502 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
503 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
504 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
505
506 // Calculate Csi
507 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
508 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
509 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
510
511 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);
512 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);
513 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);
514 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);
515 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);
516 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);
517
518 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
519 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
520 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
521 }
522 }
523 }
524
525 // Calculate Eta
526 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
527 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
528 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
529
530 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);
531 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);
532 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);
533 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);
534 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);
535 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);
536
537 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
538 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
539 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
540 }
541 }
542 }
543
544 // Calculate Zet
545 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
546 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
547 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
548
549 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);
550 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);
551 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);
552 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);
553 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);
554 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);
555
556 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
557 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
558 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
559 }
560 }
561 }
562
563 // --- 4. Boundary Extrapolation ---
564 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
565 PetscInt i_bnd, j_bnd, k_bnd;
566
567 if (xs == 0) { // If this rank owns the global i=0 boundary
568 i_bnd = 0;
569 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
570 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
571 if (i_bnd + 1 < mx) {
572 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
573 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
574 }
575 }
576 }
577 }
578 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
579 i_bnd = mx - 1;
580 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
581 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
582 if (i_bnd - 1 >= 0) {
583 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
584 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
585 }
586 }
587 }
588 }
589 if (ys == 0) {
590 j_bnd = 0;
591 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
592 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
593 if (j_bnd + 1 < my) {
594 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
595 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
596 }
597 }
598 }
599 }
600 if (ye == my) {
601 j_bnd = my - 1;
602 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
603 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
604 if (j_bnd - 1 >= 0) {
605 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
606 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
607 }
608 }
609 }
610 }
611 if (zs == 0) {
612 k_bnd = 0;
613 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
614 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
615 if (k_bnd + 1 < mz) {
616 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
617 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
618 }
619 }
620 }
621 }
622 if (ze == mz) {
623 k_bnd = mz - 1;
624 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
625 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
626 if (k_bnd - 1 >= 0) {
627 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
628 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
629 }
630 }
631 }
632 }
633
634 if (info.xs==0 && info.ys==0 && info.zs==0) {
635 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
636 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
637 }
638
639 // --- 5. Restore all arrays ---
640 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
641 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
642 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
643 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
644
645 // --- 6. Assemble Global Vectors ---
646 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
647 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
648 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
649 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
650
651 // --- 7. Update Local Ghosted Versions ---
652 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
653 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
654 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
655 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
656
657 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
658
659 PROFILE_FUNCTION_END;
660
661 PetscFunctionReturn(0);
662}
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
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), including boundary extrapolation, stores it in user->Aj, assembles user->Aj, and updates user->lAj.

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

Parameters
[in,out]userPointer to the UserCtx structure.
Returns
PetscErrorCode 0 on success.

Definition at line 676 of file Metric.c.

677{
678 PetscErrorCode ierr;
679 DMDALocalInfo info;
680 PetscScalar ***aj_arr;
681 Cmpnts ***nodal_coords_arr;
682 Vec localCoords_from_dm;
683
684 PetscFunctionBeginUser;
685 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
686
687 // --- 1. Get Nodal Coordinates and Output Array ---
688 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
689 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
690 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
691 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
692
693 // Define owned node ranges (global indices)
694 PetscInt xs = info.xs, xe = info.xs + info.xm;
695 PetscInt ys = info.ys, ye = info.ys + info.ym;
696 PetscInt zs = info.zs, ze = info.zs + info.zm;
697
698 // Global domain dimensions (total number of nodes)
699 PetscInt mx = info.mx;
700 PetscInt my = info.my;
701 PetscInt mz = info.mz;
702
703 // --- 2. Calculate Aj for INTERIOR Stencils ---
704
705 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
706 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
707 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
708
709 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
710 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
711 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
712
713 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
714 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
715 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
716
717 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) );
718
719 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) );
720
721 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) );
722
723 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) );
724
725 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) );
726
727 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) );
728
729 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) );
730
731 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) );
732
733 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) );
734
735 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);
736 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
737 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
738 }
739 }
740 }
741
742 // --- 4. Boundary Extrapolation for Aj ---
743 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
744 PetscInt i_bnd, j_bnd, k_bnd;
745
746 if (xs == 0) {
747 i_bnd = 0;
748 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
749 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
750 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
751 }
752 }
753 }
754 if (xe == mx) {
755 i_bnd = mx - 1;
756 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
757 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
758 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
759 }
760 }
761 }
762 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
763 if (ys == 0) {
764 j_bnd = 0;
765 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
766 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
767 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
768 }
769 }
770 }
771 if (ye == my) {
772 j_bnd = my - 1;
773 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
774 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
775 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
776 }
777 }
778 }
779 if (zs == 0) {
780 k_bnd = 0;
781 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
782 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
783 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
784 }
785 }
786 }
787 if (ze == mz) {
788 k_bnd = mz - 1;
789 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
790 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
791 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
792 }
793 }
794 }
795
796 // --- 5. Restore arrays ---
797 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
798 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
799
800 // --- 6. Assemble Global Vector ---
801 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
802 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
803 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
804
805 // --- 7. Update Local Ghosted Version ---
806 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
807 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
808
809 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
810 PetscFunctionReturn(0);
811}
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
userThe 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.

Definition at line 829 of file Metric.c.

830{
831 PetscErrorCode ierr;
832 DMDALocalInfo info;
833 Vec lCoords;
834 const Cmpnts ***coor;
835 Cmpnts ***cent, ***gs;
836 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
837
838 PetscFunctionBeginUser;
839
840 PROFILE_FUNCTION_BEGIN;
841
842 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);
843
844 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
845 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
846 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
847
848 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
849 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
850
851 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
852 for (PetscInt k=info.zs+1; k<info.zs+info.zm; k++) {
853 for (PetscInt j=info.ys+1; j<info.ys+info.ym; j++) {
854 for (PetscInt i=info.xs+1; i<info.xs+info.xm; i++) {
855 // Calculate cell center as the average of its 8 corner nodes
856 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);
857 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);
858 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);
859
860 // Calculate Grid Spacing in i-direction (distance between i-face centers)
861 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);
862 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);
863 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);
864 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);
865 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);
866 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);
867 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
868
869 // Calculate Grid Spacing in j-direction (distance between j-face centers)
870 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);
871 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);
872 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);
873 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);
874 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);
875 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);
876 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
877
878 // Calculate Grid Spacing in k-direction (distance between k-face centers)
879 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);
880 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);
881 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);
882 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);
883 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);
884 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);
885 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
886 }
887 }
888 }
889
890 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
891 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
892 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
893
894 // Assemble and update ghost regions for the new data
895 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
896 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
897 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
898 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
899
900 PROFILE_FUNCTION_END;
901
902 PetscFunctionReturn(0);
903}
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:705
SimCtx::rank
PetscMPIInt rank
Definition variables.h:541
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
UserCtx::_this
PetscInt _this
Definition variables.h:674
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:721
UserCtx::Cent
Vec Cent
Definition variables.h:705
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 geometric center of each constant-i face. This is a critical step for staggered-grid finite difference schemes.

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

  1. It first calculates the physical (x,y,z) coordinates of the center of each i-face and stores them in the user->Centx vector.
  2. It then uses a boundary-aware, second-order finite difference stencil on the Centx field to compute the derivatives (e.g., d(x)/d(csi)).
    • Central differences are used in the grid interior.
    • One-sided differences are used at the physical domain boundaries.
  3. Finally, these derivatives are used to compute the final metric terms and the inverse Jacobian, which are stored in their respective Vec objects.
Parameters
userThe UserCtx for a specific grid level. This function populates the user->ICsi, user->IEta, user->IZet, and user->IAj vectors.
Returns
PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 929 of file Metric.c.

930{
931 PetscErrorCode ierr;
932 DMDALocalInfo info;
933 Vec lCoords;
934 const Cmpnts ***coor;
935 Cmpnts ***centx; //***gs;
936 const Cmpnts ***centx_const;
937 Cmpnts ***icsi, ***ieta, ***izet;
938 PetscScalar ***iaj;
939 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
940
941 PetscFunctionBeginUser;
942
943 PROFILE_FUNCTION_BEGIN;
944
945 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
946
947 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
948 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
949 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
950 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
951 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
952 PetscInt gys = info.gys, gye = info.gys + info.gym;
953 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
954
955 PetscInt lxs = xs; PetscInt lxe = xe;
956 PetscInt lys = ys; PetscInt lye = ye;
957 PetscInt lzs = zs; PetscInt lze = ze;
958
959 if (xs==0) lxs = xs+1;
960 if (ys==0) lys = ys+1;
961 if (zs==0) lzs = zs+1;
962
963 if (xe==mx) lxe=xe-1;
964 if (ye==my) lye=ye-1;
965 if (ze==mz) lze=ze-1;
966
967 // --- Part 1: Calculate the location of i-face centers (Centx) ---
968 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
969 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
970 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
971 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
972
973 //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);
974
975 // Loop over the ghosted region to calculate all local face centers
976 for (PetscInt k = gzs + 1; k < gze; k++) {
977 for (PetscInt j = gys + 1; j < gye; j++) {
978 for (PetscInt i = gxs; i < gxe; i++) {
979 //----- DEBUG ------
980 //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);
981 //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,
982 // coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
983 // coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
984 // coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
985 // coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
986
987 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);
988 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);
989 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);
990
991 //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);
992 }
993 }
994 }
995
996 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
997 /*
998 if(xs==0){
999 for(PetscInt k=gzs+1;k < gze; k++){
1000 for(PetscInt j=gys+1;j < gye; j++){
1001 PetscInt i=0;
1002 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
1003 centx[k][j][i-1].y=centx[k][j][i].y;
1004 centx[k][j][i-1].z=centx[k][j][i].z;
1005 }
1006 }
1007 }
1008 if (xe==mx){
1009 for(PetscInt k=gzs+1; k<gze; k++) {
1010 for (PetscInt j=gys+1; j<gye;j++) {
1011 PetscInt i=mx-1;
1012 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
1013 centx[k][j][i].y=centx[k][j][i-1].y;
1014 centx[k][j][i].z=centx[k][j][i-1].z;
1015 }
1016 }
1017 }
1018 */
1019
1020 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1021 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1022
1023 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1024
1025 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1026
1027 // --- Part 2: Calculate metrics using face-centered coordinates ---
1028 ierr = DMDAVecGetArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1029 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1030 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1031 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1032 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1033
1034 // Loop over the OWNED region where we will store the final metrics
1035 for (PetscInt k=lzs; k<lze; k++) {
1036 for (PetscInt j=lys; j<lye; j++) {
1037 for (PetscInt i=xs; i<lxe; i++) {
1038
1039 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1040 if (i == 0) { // Forward difference at the domain's min-i boundary
1041 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
1042 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
1043 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
1044 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1045 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
1046 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
1047 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
1048 } else { // Central difference in the interior
1049 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
1050 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
1051 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
1052 }
1053
1054 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1055 if (j == 1) { // Forward difference
1056 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
1057 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
1058 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
1059 } else if (j == my - 2) { // Backward difference
1060 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
1061 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
1062 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
1063 } else { // Central difference
1064 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
1065 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
1066 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
1067 }
1068
1069 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1070 if (k == 1) { // Forward difference
1071 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
1072 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
1073 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
1074 } else if (k == mz - 2) { // Backward difference
1075 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
1076 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
1077 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
1078 } else { // Central difference
1079 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
1080 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
1081 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
1082 }
1083
1084 // --- Metric calculations (identical to legacy FormMetrics) ---
1085 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1086 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1087 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1088
1089 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1090 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1091 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1092
1093 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
1094 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1095 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
1096
1097 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
1098 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
1099 iaj[k][j][i] = 1.0 / iaj[k][j][i];
1100 }
1101 }
1102 }
1103 }
1104
1105 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1106 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1107 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1108 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1109 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1110
1111 // --- Part 3: Assemble global vectors and update local ghosts ---
1112 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
1113 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
1114 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
1115 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
1116
1117 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
1118 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
1119 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
1120 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
1121
1122 PROFILE_FUNCTION_END;
1123
1124 PetscFunctionReturn(0);
1125}
UserCtx::IZet
Vec IZet
Definition variables.h:707
UserCtx::IEta
Vec IEta
Definition variables.h:707
UserCtx::Centx
Vec Centx
Definition variables.h:706
UserCtx::ICsi
Vec ICsi
Definition variables.h:707
UserCtx::IAj
Vec IAj
Definition variables.h:707
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
userThe 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.

Definition at line 1151 of file Metric.c.

1152{
1153 PetscErrorCode ierr;
1154 DMDALocalInfo info;
1155 Vec lCoords;
1156 const Cmpnts ***coor;
1157 Cmpnts ***centy; //***gs;
1158 const Cmpnts ***centy_const;
1159 Cmpnts ***jcsi, ***jeta, ***jzet;
1160 PetscScalar ***jaj;
1161 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1162
1163 PetscFunctionBeginUser;
1164
1165 PROFILE_FUNCTION_BEGIN;
1166
1167 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1168
1169 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1170 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1171 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1172 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1173 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1174 PetscInt gys = info.gys, gye = info.gys + info.gym;
1175 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1176
1177 PetscInt lxs = xs; PetscInt lxe = xe;
1178 PetscInt lys = ys; PetscInt lye = ye;
1179 PetscInt lzs = zs; PetscInt lze = ze;
1180
1181 if (xs==0) lxs = xs+1;
1182 if (ys==0) lys = ys+1;
1183 if (zs==0) lzs = zs+1;
1184
1185 if (xe==mx) lxe=xe-1;
1186 if (ye==my) lye=ye-1;
1187 if (ze==mz) lze=ze-1;
1188
1189 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1190 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1191 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1192 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1193 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1194
1195 // Loop over the ghosted region to calculate all local face centers
1196 for (PetscInt k = gzs + 1; k < gze; k++) {
1197 for (PetscInt j = gys; j < gye; j++) {
1198 for (PetscInt i = gxs + 1; i < gxe; i++) {
1199 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);
1200 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);
1201 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);
1202 }
1203 }
1204 }
1205
1206 /*
1207 if(ys==0){
1208 for(PetscInt k=gzs+1;k < gze; k++){
1209 for(PetscInt i=gxs+1;j < gxe; i++){
1210 PetscInt j=0;
1211 centy[k][j-1][i].x=centy[k][j][i].x;
1212 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1213 centy[k][j-1][i].z=centy[k][j][i].z;
1214 }
1215 }
1216 }
1217 if (ye==my){
1218 for(PetscInt k=gzs+1; k<gze; k++) {
1219 for (PetscInt i=gxs+1; j<gxe;i++) {
1220 PetscInt j=my-1;
1221 centy[k][j][i].x=centy[k][j-1][i].x
1222 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1223 centy[k][j][i].z=centy[k][j-1][i].z;
1224 }
1225 }
1226 }
1227 */
1228
1229 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1230 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1231 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1232
1233 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1234
1235 // --- Part 2: Calculate metrics using face-centered coordinates ---
1236 ierr = DMDAVecGetArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1237 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1238 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1239 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1240 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1241
1242 // Loop over the OWNED region where we will store the final metrics
1243 for (PetscInt k=lzs; k<lze; k++) {
1244 for (PetscInt j=ys; j<lye; j++) {
1245 for (PetscInt i=lxs; i<lxe; i++) {
1246
1247 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1248 if (i == 1) { // Forward difference at the domain's min-i boundary
1249 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1250 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1251 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1252 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1253 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1254 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1255 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1256 } else { // Central difference in the interior
1257 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1258 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1259 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1260 }
1261
1262 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1263 if (j == 0) { // Forward difference
1264 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1265 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1266 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1267 } else if (j == my - 2) { // Backward difference
1268 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1269 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1270 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1271 } else { // Central difference
1272 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1273 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1274 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1275 }
1276
1277 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1278 if (k == 1) { // Forward difference
1279 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1280 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1281 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1282 } else if (k == mz - 2) { // Backward difference
1283 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1284 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1285 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1286 } else { // Central difference
1287 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1288 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1289 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1290 }
1291
1292 // --- Metric calculations (identical to legacy FormMetrics) ---
1293 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1294 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1295 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1296
1297 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1298 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1299 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1300
1301 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1302 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1303 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1304
1305 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1306 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1307 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1308 }
1309 }
1310 }
1311 }
1312
1313 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1314 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1315 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1316 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1317 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1318
1319 // --- Part 3: Assemble global vectors and update local ghosts ---
1320 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1321 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1322 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1323 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1324
1325 ierr = UpdateLocalGhosts(user, "JCsi"); CHKERRQ(ierr);
1326 ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
1327 ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
1328 ierr = UpdateLocalGhosts(user, "JAj"); CHKERRQ(ierr);
1329
1330 PROFILE_FUNCTION_END;
1331
1332 PetscFunctionReturn(0);
1333}
UserCtx::JCsi
Vec JCsi
Definition variables.h:708
UserCtx::JEta
Vec JEta
Definition variables.h:708
UserCtx::JZet
Vec JZet
Definition variables.h:708
UserCtx::JAj
Vec JAj
Definition variables.h:708
UserCtx::Centy
Vec Centy
Definition variables.h:706
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-k face. This is a critical step for staggered-grid finite difference schemes.

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

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

Definition at line 1359 of file Metric.c.

1360{
1361 PetscErrorCode ierr;
1362 DMDALocalInfo info;
1363 Vec lCoords;
1364 const Cmpnts ***coor;
1365 Cmpnts ***centz; //***gs;
1366 const Cmpnts ***centz_const;
1367 Cmpnts ***kcsi, ***keta, ***kzet;
1368 PetscScalar ***kaj;
1369 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1370
1371 PetscFunctionBeginUser;
1372
1373 PROFILE_FUNCTION_BEGIN;
1374
1375 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1376
1377 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1378 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1379 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1380 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1381 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1382 PetscInt gys = info.gys, gye = info.gys + info.gym;
1383 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1384
1385 PetscInt lxs = xs; PetscInt lxe = xe;
1386 PetscInt lys = ys; PetscInt lye = ye;
1387 PetscInt lzs = zs; PetscInt lze = ze;
1388
1389 if (xs==0) lxs = xs+1;
1390 if (ys==0) lys = ys+1;
1391 if (zs==0) lzs = zs+1;
1392
1393 if (xe==mx) lxe=xe-1;
1394 if (ye==my) lye=ye-1;
1395 if (ze==mz) lze=ze-1;
1396
1397 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1398 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1399 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1400 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1401 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1402
1403 // Loop over the ghosted region to calculate all local face centers
1404 for (PetscInt k = gzs; k < gze; k++) {
1405 for (PetscInt j = gys; j < gye; j++) {
1406 for (PetscInt i = gxs + 1; i < gxe; i++) {
1407 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);
1408 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);
1409 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);
1410 }
1411 }
1412 }
1413
1414 /*
1415 if(zs==0){
1416 for(PetscInt j=gys+1;j < gye; j++){
1417 for(PetscInt i=gxs+1;j < gxe; i++){
1418 PetscInt k=0;
1419 centz[k-1][j][i].x=centz[k][j][i].x;
1420 centz[k-1][j][i].y=centz[k][j][i].y;
1421 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1422 }
1423 }
1424 }
1425 if (ze==mz){
1426 for(PetscInt j=gys+1; j<gye; j++) {
1427 for (PetscInt i=gxs+1; j<gxe;i++) {
1428 PetscInt k=mz-1;
1429 centy[k][j][i].x=centy[k-1][j][i].x
1430 centy[k][j][i].y=centy[k-1][j][i].y;
1431 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1432 }
1433 }
1434 }
1435 */
1436
1437 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1438 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1439 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1440
1441 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1442
1443 // --- Part 2: Calculate metrics using face-centered coordinates ---
1444 ierr = DMDAVecGetArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1445 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1446 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1447 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1448 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1449
1450 // Loop over the OWNED region where we will store the final metrics
1451 for (PetscInt k=zs; k<lze; k++) {
1452 for (PetscInt j=lys; j<lye; j++) {
1453 for (PetscInt i=lxs; i<lxe; i++) {
1454
1455 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1456 if (i == 1) { // Forward difference at the domain's min-i boundary
1457 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
1458 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
1459 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
1460 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1461 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
1462 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
1463 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
1464 } else { // Central difference in the interior
1465 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
1466 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
1467 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
1468 }
1469
1470 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1471 if (j == 1) { // Forward difference
1472 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
1473 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
1474 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
1475 } else if (j == my - 2) { // Backward difference
1476 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
1477 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
1478 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
1479 } else { // Central difference
1480 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
1481 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
1482 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
1483 }
1484
1485 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1486 if (k == 0) { // Forward difference
1487 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
1488 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
1489 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
1490 } else if (k == mz - 2) { // Backward difference
1491 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
1492 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
1493 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
1494 } else { // Central difference
1495 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
1496 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
1497 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
1498 }
1499
1500 // --- Metric calculations (identical to legacy FormMetrics) ---
1501 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1502 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1503 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1504
1505 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1506 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1507 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1508
1509 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1510 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1511 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1512
1513 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
1514 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
1515 kaj[k][j][i] = 1.0 / kaj[k][j][i];
1516 }
1517 }
1518 }
1519 }
1520
1521 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1522 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1523 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1524 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1525 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1526
1527 // --- Part 3: Assemble global vectors and update local ghosts ---
1528 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
1529 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
1530 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
1531 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
1532
1533 ierr = UpdateLocalGhosts(user, "KCsi"); CHKERRQ(ierr);
1534 ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
1535 ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
1536 ierr = UpdateLocalGhosts(user, "KAj"); CHKERRQ(ierr);
1537
1538 PROFILE_FUNCTION_END;
1539
1540 PetscFunctionReturn(0);
1541}
UserCtx::KAj
Vec KAj
Definition variables.h:709
UserCtx::Centz
Vec Centz
Definition variables.h:706
UserCtx::KEta
Vec KEta
Definition variables.h:709
UserCtx::KZet
Vec KZet
Definition variables.h:709
UserCtx::KCsi
Vec KCsi
Definition variables.h:709
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Gidx()

static PetscInt Gidx ( PetscInt  i,
PetscInt  j,
PetscInt  k,
UserCtx *  user 
)
static

Definition at line 1543 of file Metric.c.

1544{
1545 PetscInt nidx;
1546 DMDALocalInfo info = user->info;
1547
1548 PetscInt mx = info.mx, my = info.my;
1549
1550 AO ao;
1551 DMDAGetAO(user->da, &ao);
1552 nidx=i+j*mx+k*mx*my;
1553
1554 AOApplicationToPetsc(ao,1,&nidx);
1555
1556 return (nidx);
1557}
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
Here is the caller graph for this function:

◆ ComputeMetricsDivergence()

PetscErrorCode ComputeMetricsDivergence ( UserCtx *  user)

Computes the divergence of the grid metrics and identifies the maximum value.

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

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

Parameters
userThe UserCtx, containing all necessary grid data.
Returns
PetscErrorCode

Definition at line 1573 of file Metric.c.

1574{
1575 DM da = user->da, fda = user->fda;
1576 DMDALocalInfo info = user->info;
1577 PetscInt xs = info.xs, xe = info.xs + info.xm;
1578 PetscInt ys = info.ys, ye = info.ys + info.ym;
1579 PetscInt zs = info.zs, ze = info.zs + info.zm;
1580 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1581 PetscInt lxs, lys, lzs, lxe, lye, lze;
1582 PetscInt i, j, k;
1583 Vec Div;
1584 PetscReal ***div, ***aj;
1585 Cmpnts ***csi, ***eta, ***zet;
1586 PetscReal maxdiv;
1587
1588 PetscFunctionBeginUser;
1589
1590 PROFILE_FUNCTION_BEGIN;
1591
1592 lxs = xs; lxe = xe;
1593 lys = ys; lye = ye;
1594 lzs = zs; lze = ze;
1595
1596 if (xs == 0) lxs = xs + 1;
1597 if (ys == 0) lys = ys + 1;
1598 if (zs == 0) lzs = zs + 1;
1599
1600 if (xe == mx) lxe = xe - 1;
1601 if (ye == my) lye = ye - 1;
1602 if (ze == mz) lze = ze - 1;
1603
1604 DMDAVecGetArray(fda, user->lCsi, &csi);
1605 DMDAVecGetArray(fda, user->lEta, &eta);
1606 DMDAVecGetArray(fda, user->lZet, &zet);
1607 DMDAVecGetArray(da, user->lAj, &aj);
1608
1609 VecDuplicate(user->P, &Div);
1610 VecSet(Div, 0.);
1611 DMDAVecGetArray(da, Div, &div);
1612
1613 for (k = lzs; k < lze; k++) {
1614 for (j = lys; j < lye; j++) {
1615 for (i = lxs; i < lxe; i++) {
1616 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1617 eta[k][j][i].x - eta[k][j-1][i].x +
1618 zet[k][j][i].x - zet[k-1][j][i].x +
1619 csi[k][j][i].y - csi[k][j][i-1].y +
1620 eta[k][j][i].y - eta[k][j-1][i].y +
1621 zet[k][j][i].y - zet[k-1][j][i].y +
1622 csi[k][j][i].z - csi[k][j][i-1].z +
1623 eta[k][j][i].z - eta[k][j-1][i].z +
1624 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1625 div[k][j][i] = fabs(divergence);
1626 }
1627 }
1628 }
1629
1630 DMDAVecRestoreArray(da, Div, &div);
1631
1632 PetscReal MaxFlatIndex;
1633 VecMax(Div, &MaxFlatIndex, &maxdiv);
1634 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %d.\n",maxdiv,MaxFlatIndex);
1635
1636 for (k=zs; k<ze; k++) {
1637 for (j=ys; j<ye; j++) {
1638 for (i=xs; i<xe; i++) {
1639 if (Gidx(i,j,k,user) == MaxFlatIndex) {
1640 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,k,j,i);
1641 }
1642 }
1643 }
1644 }
1645
1646
1647 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1648 DMDAVecRestoreArray(fda, user->lEta, &eta);
1649 DMDAVecRestoreArray(fda, user->lZet, &zet);
1650 DMDAVecRestoreArray(da, user->lAj, &aj);
1651 VecDestroy(&Div);
1652
1653
1654 PROFILE_FUNCTION_END;
1655
1656 PetscFunctionReturn(0);
1657}
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition Metric.c:1543
UserCtx::lZet
Vec lZet
Definition variables.h:705
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::lEta
Vec lEta
Definition variables.h:705
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
userThe UserCtx, containing all necessary grid data.
Returns
PetscErrorCode

Definition at line 1670 of file Metric.c.

1671{
1672
1673 DMDALocalInfo info = user->info;
1674 PetscInt xs = info.xs, xe = info.xs + info.xm;
1675 PetscInt ys = info.ys, ye = info.ys + info.ym;
1676 PetscInt zs = info.zs, ze = info.zs + info.zm;
1677 PetscInt i, j, k;
1678
1679 PetscFunctionBeginUser;
1680
1681 PROFILE_FUNCTION_BEGIN;
1682
1683 PetscReal CsiMax, EtaMax, ZetMax;
1684 PetscReal ICsiMax, IEtaMax, IZetMax;
1685 PetscReal JCsiMax, JEtaMax, JZetMax;
1686 PetscReal KCsiMax, KEtaMax, KZetMax;
1687 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
1688
1689 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
1690 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
1691 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
1692 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
1693 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
1694
1695 // Max Values
1696 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
1697 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
1698 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
1699
1700 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
1701 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
1702 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
1703
1704 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
1705 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
1706 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
1707
1708 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
1709 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
1710 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
1711
1712 VecMax(user->lAj,&AjMaxArg,&AjMax);
1713 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1714 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1715 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1716
1717 VecMax(user->lAj,&AjMaxArg,&AjMax);
1718 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1719 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1720 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1721
1722 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
1723
1724 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
1725 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
1726 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
1727 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
1728 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
1729
1730 for (k=zs; k<ze; k++) {
1731 for (j=ys; j<ye; j++) {
1732 for (i=xs; i<xe; i++) {
1733 if (Gidx(i,j,k,user) == CsiMaxArg) {
1734 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
1735 }
1736 if (Gidx(i,j,k,user) == EtaMaxArg) {
1737 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
1738 }
1739 if (Gidx(i,j,k,user) == ZetMaxArg) {
1740 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
1741 }
1742 if (Gidx(i,j,k,user) == ICsiMaxArg) {
1743 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
1744 }
1745 if (Gidx(i,j,k,user) == IEtaMaxArg) {
1746 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
1747 }
1748 if (Gidx(i,j,k,user) == IZetMaxArg) {
1749 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
1750 }
1751 if (Gidx(i,j,k,user) == JCsiMaxArg) {
1752 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
1753 }
1754 if (Gidx(i,j,k,user) == JEtaMaxArg) {
1755 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
1756 }
1757 if (Gidx(i,j,k,user) == JZetMaxArg) {
1758 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
1759 }
1760 if (Gidx(i,j,k,user) == KCsiMaxArg) {
1761 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
1762 }
1763 if (Gidx(i,j,k,user) == KEtaMaxArg) {
1764 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
1765 }
1766 if (Gidx(i,j,k,user) == KZetMaxArg) {
1767 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
1768 }
1769 if (Gidx(i,j,k,user) == AjMaxArg) {
1770 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
1771 }
1772 if (Gidx(i,j,k,user) == IAjMaxArg) {
1773 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
1774 }
1775 if (Gidx(i,j,k,user) == JAjMaxArg) {
1776 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
1777 }
1778 if (Gidx(i,j,k,user) == KAjMaxArg) {
1779 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
1780 }
1781 }
1782 }
1783 }
1784
1785 /*
1786 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
1787 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
1788 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
1789 */
1790
1791 PROFILE_FUNCTION_END;
1792
1793 PetscFunctionReturn(0);
1794}
UserCtx::lIEta
Vec lIEta
Definition variables.h:707
UserCtx::lIZet
Vec lIZet
Definition variables.h:707
UserCtx::lIAj
Vec lIAj
Definition variables.h:707
UserCtx::lKEta
Vec lKEta
Definition variables.h:709
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:708
UserCtx::lKZet
Vec lKZet
Definition variables.h:709
UserCtx::lJEta
Vec lJEta
Definition variables.h:708
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:709
UserCtx::lJZet
Vec lJZet
Definition variables.h:708
UserCtx::lICsi
Vec lICsi
Definition variables.h:707
UserCtx::lJAj
Vec lJAj
Definition variables.h:708
UserCtx::lKAj
Vec lKAj
Definition variables.h:709
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
simCtxThe master SimCtx, containing the configured UserCtx hierarchy.
Returns
PetscErrorCode

Definition at line 1809 of file Metric.c.

1810{
1811 PetscErrorCode ierr;
1812 UserMG *usermg = &simCtx->usermg;
1813 MGCtx *mgctx = usermg->mgctx;
1814 PetscInt nblk = simCtx->block_number;
1815
1816 PetscFunctionBeginUser;
1817
1818 PROFILE_FUNCTION_BEGIN;
1819
1820 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
1821
1822 // Loop through all levels and all blocks
1823 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
1824 for (PetscInt bi = 0; bi < nblk; bi++) {
1825 UserCtx *user = &mgctx[level].user[bi];
1826 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
1827
1828 // Call the modern, modular helper functions for each UserCtx.
1829 // These functions are self-contained and operate on the data within the provided context.
1830 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
1831 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
1832 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
1833 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
1834 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
1835 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
1836 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
1837
1838 // Diagnostics
1839 ierr = ComputeMetricNorms(user);
1840 if (level == usermg->mglevels - 1) {
1841 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
1842 }
1843 }
1844 }
1845
1846 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
1847
1848 PROFILE_FUNCTION_END;
1849
1850 PetscFunctionReturn(0);
1851}
ComputeMetricNorms
PetscErrorCode ComputeMetricNorms(UserCtx *user)
Computes the max-min values of the grid metrics.
Definition Metric.c:1670
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain...
Definition Metric.c:368
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Computes the physical location of cell centers and the spacing between them.
Definition Metric.c:829
ComputeJFaceMetrics
PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
Computes metrics centered on constant-j faces (j-faces).
Definition Metric.c:1151
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Computes the divergence of the grid metrics and identifies the maximum value.
Definition Metric.c:1573
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation,...
Definition Metric.c:464
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Calculates the cell-centered inverse Jacobian determinant (1/J), including boundary extrapolation,...
Definition Metric.c:676
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Computes metrics centered on constant-k faces (k-faces).
Definition Metric.c:1359
ComputeIFaceMetrics
PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
Computes metrics centered on constant-i faces (i-faces).
Definition Metric.c:929
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:268
MGCtx::user
UserCtx * user
Definition variables.h:427
SimCtx::block_number
PetscInt block_number
Definition variables.h:593
SimCtx::usermg
UserMG usermg
Definition variables.h:631
UserMG::mglevels
PetscInt mglevels
Definition variables.h:434
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:437
MGCtx
Context for Multigrid operations.
Definition variables.h:426
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:661
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:433
Here is the call graph for this function:
Here is the caller graph for this function: