PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
Metric.h File Reference
#include <petsc.h>
#include "variables.h"
#include "logging.h"
#include <stdlib.h>
#include "io.h"
#include "setup.h"
Include dependency graph for Metric.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

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

Function Documentation

◆ MetricLogicalToPhysical()

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

Public interface for MetricLogicalToPhysical().

Parameters
userPrimary UserCtx input for the operation.
XParameter X passed to MetricLogicalToPhysical().
iParameter i passed to MetricLogicalToPhysical().
jParameter j passed to MetricLogicalToPhysical().
kParameter k passed to MetricLogicalToPhysical().
xiParameter xi passed to MetricLogicalToPhysical().
etaParameter eta passed to MetricLogicalToPhysical().
ztaParameter zta passed to MetricLogicalToPhysical().
XpParameter Xp passed to MetricLogicalToPhysical().
Returns
PetscErrorCode 0 on success.

Public interface for MetricLogicalToPhysical().

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

See also
MetricLogicalToPhysical()

Definition at line 77 of file Metric.c.

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)
Internal helper implementation: TrilinearBlend().
Definition Metric.c:51
MetricGetCellVertices
PetscErrorCode MetricGetCellVertices(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
Implementation of MetricGetCellVertices().
Definition Metric.c:26
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770
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:

◆ MetricGetCellVertices()

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

Public interface for MetricGetCellVertices().

Parameters
userPrimary UserCtx input for the operation.
XParameter X passed to MetricGetCellVertices().
iParameter i passed to MetricGetCellVertices().
jParameter j passed to MetricGetCellVertices().
kParameter k passed to MetricGetCellVertices().
VParameter V passed to MetricGetCellVertices().
Returns
PetscErrorCode 0 on success.

Public interface for MetricGetCellVertices().

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

See also
MetricGetCellVertices()

Definition at line 26 of file Metric.c.

30{
31 (void)user;
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:297
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
Here is the caller graph for this function:

◆ MetricJacobian()

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

Public interface for MetricJacobian().

Parameters
userPrimary UserCtx input for the operation.
XParameter X passed to MetricJacobian().
iParameter i passed to MetricJacobian().
jParameter j passed to MetricJacobian().
kParameter k passed to MetricJacobian().
xiParameter xi passed to MetricJacobian().
etaParameter eta passed to MetricJacobian().
ztaParameter zta passed to MetricJacobian().
JParameter J passed to MetricJacobian().
detJParameter detJ passed to MetricJacobian().
Returns
PetscErrorCode 0 on success.

Public interface for MetricJacobian().

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

See also
MetricJacobian()

Definition at line 106 of file Metric.c.

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

◆ MetricVelocityContravariant()

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

Public interface for MetricVelocityContravariant().

Parameters
JParameter J passed to MetricVelocityContravariant().
detJParameter detJ passed to MetricVelocityContravariant().
uParameter u passed to MetricVelocityContravariant().
ucParameter uc passed to MetricVelocityContravariant().
Returns
PetscErrorCode 0 on success.

Public interface for MetricVelocityContravariant().

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

See also
MetricVelocityContravariant()

Definition at line 166 of file Metric.c.

168{
169 PetscFunctionBeginUser;
170
171 PROFILE_FUNCTION_BEGIN;
172
173 /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
174 PetscReal gxi[3] = { J[1][1]*J[2][2]-J[1][2]*J[2][1],
175 -J[0][1]*J[2][2]+J[0][2]*J[2][1],
176 J[0][1]*J[1][2]-J[0][2]*J[1][1] };
177 PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
178 J[0][0]*J[2][2]-J[0][2]*J[2][0],
179 -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
180 PetscReal gzta[3] = { J[1][0]*J[2][1]-J[1][1]*J[2][0],
181 -J[0][0]*J[2][1]+J[0][1]*J[2][0],
182 J[0][0]*J[1][1]-J[0][1]*J[1][0] };
183
184 PetscReal invDet = 1.0 / detJ;
185 for (int d=0; d<3; ++d) { gxi[d] *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
186
187 uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
188 uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
189 uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
190
191 PROFILE_FUNCTION_END;
192
193 PetscFunctionReturn(0);
194}
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
csiParameter csi passed to CalculateFaceNormalAndArea().
etaParameter eta passed to CalculateFaceNormalAndArea().
zetParameter zet passed to CalculateFaceNormalAndArea().
niOutput:
njOutput:
nkOutput:
AiOutput:
AjOutput:
AkOutput:
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 237 of file Metric.c.

238{
239 PetscFunctionBeginUser;
240 PROFILE_FUNCTION_BEGIN;
241 double g[3][3];
242 double G[3][3];
243
244 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
245 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
246 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
247
248 InvertCovariantMetricTensor(g, G);
249 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
250 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
251 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
252
253 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
254 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
255 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
256
257 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
258 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
259 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
260
261 *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
262 *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
263 *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
264
265 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
266 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
267 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
268
269 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
270 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
271 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
272
273 PROFILE_FUNCTION_END;
274 PetscFunctionReturn(0);
275}
InvertCovariantMetricTensor
PetscErrorCode InvertCovariantMetricTensor(double covariantTensor[3][3], double contravariantTensor[3][3])
Internal helper implementation: InvertCovariantMetricTensor().
Definition Metric.c:203
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 call graph for this function:
Here is the caller graph for this function:

◆ InvertCovariantMetricTensor()

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

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

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

Parameters
covariantTensorInput: A 3x3 matrix representing the covariant metric tensor.
contravariantTensorOutput: A 3x3 matrix where the inverted result is stored.
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 203 of file Metric.c.

204{
205 PetscFunctionBeginUser;
206
207 const double a11=covariantTensor[0][0], a12=covariantTensor[0][1], a13=covariantTensor[0][2];
208 const double a21=covariantTensor[1][0], a22=covariantTensor[1][1], a23=covariantTensor[1][2];
209 const double a31=covariantTensor[2][0], a32=covariantTensor[2][1], a33=covariantTensor[2][2];
210
211 double det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
212
213 if (fabs(det) < 1.0e-12) {
214 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Matrix is singular, determinant is near zero.");
215 }
216
217 contravariantTensor[0][0] = (a33*a22-a32*a23)/det;
218 contravariantTensor[0][1] = -(a33*a12-a32*a13)/det;
219 contravariantTensor[0][2] = (a23*a12-a22*a13)/det;
220 contravariantTensor[1][0] = -(a33*a21-a31*a23)/det;
221 contravariantTensor[1][1] = (a33*a11-a31*a13)/det;
222 contravariantTensor[1][2] = -(a23*a11-a21*a13)/det;
223 contravariantTensor[2][0] = (a32*a21-a31*a22)/det;
224 contravariantTensor[2][1] = -(a32*a11-a31*a12)/det;
225 contravariantTensor[2][2] = (a22*a11-a21*a12)/det;
226
227 PetscFunctionReturn(0);
228}
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
csiParameter csi passed to ComputeCellCharacteristicLengthScale().
etaParameter eta passed to ComputeCellCharacteristicLengthScale().
zetParameter zet passed to ComputeCellCharacteristicLengthScale().
dxOutput:
dyOutput:
dzOutput:
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 283 of file Metric.c.

284{
285 PetscFunctionBeginUser;
286 PROFILE_FUNCTION_BEGIN;
287 double ni[3], nj[3], nk[3];
288 double Li, Lj, Lk;
289 double Ai, Aj, Ak;
290 double vol = 1./ajc;
291
292 CalculateFaceNormalAndArea(csi, eta, zet, ni, nj, nk, &Ai, &Aj, &Ak);
293 Li = vol / Ai;
294 Lj = vol / Aj;
295 Lk = vol / Ak;
296
297 // Length scale vector = di * ni_vector + dj * nj_vector + dk * nk_vector
298 *dx = fabs( Li * ni[0] + Lj * nj[0] + Lk * nk[0] );
299 *dy = fabs( Li * ni[1] + Lj * nj[1] + Lk * nk[1] );
300 *dz = fabs( Li * ni[2] + Lj * nj[2] + Lk * nk[2] );
301
302 PROFILE_FUNCTION_END;
303 PetscFunctionReturn(0);
304}
CalculateFaceNormalAndArea
PetscErrorCode CalculateFaceNormalAndArea(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
Internal helper implementation: CalculateFaceNormalAndArea().
Definition Metric.c:237
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ApplyPeriodicCorrectionsToCellCentersAndSpacing()

PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing ( UserCtx *  user)

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

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

Parameters
userThe UserCtx containing grid and field data.
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 392 of file Metric.c.

393{
394 PetscErrorCode ierr;
395 DMDALocalInfo info = user->info;
396 PetscInt xs = info.xs, xe = info.xs + info.xm;
397 PetscInt ys = info.ys, ye = info.ys + info.ym;
398 PetscInt zs = info.zs, ze = info.zs + info.zm;
399 PetscInt mx = info.mx, my = info.my, mz = info.mz;
400 Cmpnts ***cent, ***lcent, ***gs;
401 PetscReal delta;
402
403 PetscFunctionBeginUser;
404 PROFILE_FUNCTION_BEGIN;
405
406 // Check if any periodic boundaries exist
407 PetscBool has_periodic = PETSC_FALSE;
408 for (int i = 0; i < 6; i++) {
409 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
410 has_periodic = PETSC_TRUE;
411 break;
412 }
413 }
414
415 if (!has_periodic) {
416 LOG_ALLOW(LOCAL, LOG_TRACE, "No periodic boundaries; skipping corrections for Cent/GridSpace.\n");
417 PROFILE_FUNCTION_END;
418 PetscFunctionReturn(0);
419 }
420
421 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying periodic corrections to Cent and GridSpace.\n");
422
423 // Must update ghosts first before applying corrections
424 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
425 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
426
427 // --- X-direction periodic corrections ---
428 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
429 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
430
431 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
432 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
433 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
434
435 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
436 if (user->cgrid) {
437 for (PetscInt k=zs; k<ze; k++) {
438 for (PetscInt j=ys; j<ye; j++) {
439 cent[k][j][0] = lcent[k][j][-2];
440 }
441 }
442 } else {
443 for (PetscInt k=zs; k<ze; k++) {
444 for (PetscInt j=ys; j<ye; j++) {
445 delta = (gs[k][j][1].x + gs[k][j][-2].x) / 2.0;
446 cent[k][j][0].x = cent[k][j][1].x - delta;
447 cent[k][j][0].y = cent[k][j][1].y;
448 cent[k][j][0].z = cent[k][j][1].z;
449 }
450 }
451 }
452 }
453
454 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
455 if (user->cgrid) {
456 for (PetscInt k=zs; k<ze; k++) {
457 for (PetscInt j=ys; j<ye; j++) {
458 cent[k][j][mx-1] = lcent[k][j][mx+1];
459 }
460 }
461 } else {
462 for (PetscInt k=zs; k<ze; k++) {
463 for (PetscInt j=ys; j<ye; j++) {
464 delta = (gs[k][j][mx-2].x + gs[k][j][mx+1].x) / 2.0;
465 cent[k][j][mx-1].x = cent[k][j][mx-2].x + delta;
466 cent[k][j][mx-1].y = cent[k][j][mx-2].y;
467 cent[k][j][mx-1].z = cent[k][j][mx-2].z;
468 }
469 }
470 }
471 }
472
473 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
474 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
475 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
476 }
477 // --- Y-direction periodic corrections ---
478 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
479 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
480
481 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
482 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
483 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
484
485 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
486 if (user->cgrid) {
487 for (PetscInt k=zs; k<ze; k++) {
488 for (PetscInt i=xs; i<xe; i++) {
489 cent[k][0][i] = lcent[k][-2][i];
490 }
491 }
492 } else {
493 for (PetscInt k=zs; k<ze; k++) {
494 for (PetscInt i=xs; i<xe; i++) {
495 delta = (gs[k][1][i].y + gs[k][-2][i].y) / 2.0;
496 cent[k][0][i].x = cent[k][1][i].x;
497 cent[k][0][i].y = cent[k][1][i].y - delta;
498 cent[k][0][i].z = cent[k][1][i].z;
499 }
500 }
501 }
502 }
503
504 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
505 if (user->cgrid) {
506 for (PetscInt k=zs; k<ze; k++) {
507 for (PetscInt i=xs; i<xe; i++) {
508 cent[k][my-1][i] = lcent[k][my+1][i];
509 }
510 }
511 } else {
512 for (PetscInt k=zs; k<ze; k++) {
513 for (PetscInt i=xs; i<xe; i++) {
514 delta = (gs[k][my-2][i].y + gs[k][my+1][i].y) / 2.0;
515 cent[k][my-1][i].x = cent[k][my-2][i].x;
516 cent[k][my-1][i].y = cent[k][my-2][i].y + delta;
517 cent[k][my-1][i].z = cent[k][my-2][i].z;
518 }
519 }
520 }
521 }
522
523 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
524 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
525 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
526
527 }
528
529 // --- Z-direction periodic corrections ---
530 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
531 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
532
533 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
534 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
535 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
536
537 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
538 if (user->cgrid) {
539 for (PetscInt j=ys; j<ye; j++) {
540 for (PetscInt i=xs; i<xe; i++) {
541 cent[0][j][i] = lcent[-2][j][i];
542 }
543 }
544 } else {
545 for (PetscInt j=ys; j<ye; j++) {
546 for (PetscInt i=xs; i<xe; i++) {
547 delta = (gs[1][j][i].z + gs[-2][j][i].z) / 2.0;
548 cent[0][j][i].x = cent[1][j][i].x;
549 cent[0][j][i].y = cent[1][j][i].y;
550 cent[0][j][i].z = cent[1][j][i].z - delta;
551 }
552 }
553 }
554 }
555
556 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
557 if (user->cgrid) {
558 for (PetscInt j=ys; j<ye; j++) {
559 for (PetscInt i=xs; i<xe; i++) {
560 cent[mz-1][j][i] = lcent[mz+1][j][i];
561 }
562 }
563 } else {
564 for (PetscInt j=ys; j<ye; j++) {
565 for (PetscInt i=xs; i<xe; i++) {
566 delta = (gs[mz-2][j][i].z + gs[mz+1][j][i].z) / 2.0;
567 cent[mz-1][j][i].x = cent[mz-2][j][i].x;
568 cent[mz-1][j][i].y = cent[mz-2][j][i].y;
569 cent[mz-1][j][i].z = cent[mz-2][j][i].z + delta;
570 }
571 }
572 }
573 }
574
575 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
576 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
577 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
578
579 }
580
581 PROFILE_FUNCTION_END;
582 PetscFunctionReturn(0);
583}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
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:199
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:32
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1361
UserCtx::lCent
Vec lCent
Definition variables.h:858
PERIODIC
@ PERIODIC
Definition variables.h:260
UserCtx::cgrid
PetscInt cgrid
Definition variables.h:826
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:829
UserCtx::lGridSpace
Vec lGridSpace
Definition variables.h:858
UserCtx::info
DMDALocalInfo info
Definition variables.h:818
UserCtx::Cent
Vec Cent
Definition variables.h:858
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:247
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:246
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ApplyPeriodicCorrectionsToIFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter ( UserCtx *  user)

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

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

Parameters
userThe UserCtx containing grid and field data.
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 591 of file Metric.c.

592{
593 PetscErrorCode ierr;
594 DMDALocalInfo info = user->info;
595 PetscInt xs = info.xs, xe = info.xs + info.xm;
596 PetscInt mx = info.mx;
597 PetscInt gys = info.gys, gye = info.gys + info.gym;
598 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
599 Cmpnts ***centx, ***gs;
600
601 PetscFunctionBeginUser;
602 PROFILE_FUNCTION_BEGIN;
603
604 // Check if X-periodic boundaries exist
605 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC &&
606 user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
607 LOG_ALLOW(LOCAL, LOG_TRACE, "No X-periodic boundaries; skipping Centx corrections.\n");
608 PROFILE_FUNCTION_END;
609 PetscFunctionReturn(0);
610 }
611
612 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying X-periodic corrections to Centx.\n");
613
614 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
615 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
616
617 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
618 if (user->cgrid) {
619 for (PetscInt k=gzs+1; k<gze; k++) {
620 for (PetscInt j=gys+1; j<gye; j++) {
621 centx[k][j][-1].x = centx[k][j][-3].x;
622 centx[k][j][-1].y = centx[k][j][-3].y;
623 centx[k][j][-1].z = centx[k][j][-3].z;
624 }
625 }
626 } else {
627 for (PetscInt k=gzs+1; k<gze; k++) {
628 for (PetscInt j=gys+1; j<gye; j++) {
629 centx[k][j][-1].x = centx[k][j][0].x - gs[k][j][-2].x;
630 centx[k][j][-1].y = centx[k][j][0].y;
631 centx[k][j][-1].z = centx[k][j][0].z;
632 }
633 }
634 }
635 }
636
637 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
638 if (user->cgrid) {
639 for (PetscInt k=gzs+1; k<gze; k++) {
640 for (PetscInt j=gys+1; j<gye; j++) {
641 centx[k][j][mx-1].x = centx[k][j][mx+1].x;
642 centx[k][j][mx-1].y = centx[k][j][mx+1].y;
643 centx[k][j][mx-1].z = centx[k][j][mx+1].z;
644 }
645 }
646 } else {
647 for (PetscInt k=gzs+1; k<gze; k++) {
648 for (PetscInt j=gys+1; j<gye; j++) {
649 centx[k][j][mx-1].x = centx[k][j][mx-2].x + gs[k][j][mx+1].x;
650 centx[k][j][mx-1].y = centx[k][j][mx-2].y;
651 centx[k][j][mx-1].z = centx[k][j][mx-2].z;
652 }
653 }
654 }
655 }
656
657 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
658 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
659
660 PROFILE_FUNCTION_END;
661 PetscFunctionReturn(0);
662}
UserCtx::Centx
Vec Centx
Definition variables.h:859
Here is the caller graph for this function:

◆ ApplyPeriodicCorrectionsToJFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter ( UserCtx *  user)

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

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

Parameters
userThe UserCtx containing grid and field data.
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 671 of file Metric.c.

672{
673 PetscErrorCode ierr;
674 DMDALocalInfo info = user->info;
675 PetscInt ys = info.ys, ye = info.ys + info.ym;
676 PetscInt my = info.my;
677 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
678 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
679 Cmpnts ***centy, ***gs;
680
681 PetscFunctionBeginUser;
682 PROFILE_FUNCTION_BEGIN;
683
684 // Check if Y-periodic boundaries exist
685 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC &&
686 user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
687 LOG_ALLOW(LOCAL, LOG_TRACE, "No Y-periodic boundaries; skipping Centy corrections.\n");
688 PROFILE_FUNCTION_END;
689 PetscFunctionReturn(0);
690 }
691
692 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying Y-periodic corrections to Centy.\n");
693
694 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
695 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
696
697 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
698 if (user->cgrid) {
699 for (PetscInt k=gzs+1; k<gze; k++) {
700 for (PetscInt i=gxs+1; i<gxe; i++) {
701 centy[k][-1][i].x = centy[k][-3][i].x;
702 centy[k][-1][i].y = centy[k][-3][i].y;
703 centy[k][-1][i].z = centy[k][-3][i].z;
704 }
705 }
706 } else {
707 for (PetscInt k=gzs+1; k<gze; k++) {
708 for (PetscInt i=gxs+1; i<gxe; i++) {
709 centy[k][-1][i].x = centy[k][0][i].x;
710 centy[k][-1][i].y = centy[k][0][i].y - gs[k][-2][i].y;
711 centy[k][-1][i].z = centy[k][0][i].z;
712 }
713 }
714 }
715 }
716
717 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
718 if (user->cgrid) {
719 for (PetscInt k=gzs+1; k<gze; k++) {
720 for (PetscInt i=gxs+1; i<gxe; i++) {
721 centy[k][my-1][i].x = centy[k][my+1][i].x;
722 centy[k][my-1][i].y = centy[k][my+1][i].y;
723 centy[k][my-1][i].z = centy[k][my+1][i].z;
724 }
725 }
726 } else {
727 for (PetscInt k=gzs+1; k<gze; k++) {
728 for (PetscInt i=gxs+1; i<gxe; i++) {
729 centy[k][my-1][i].x = centy[k][my-2][i].x;
730 centy[k][my-1][i].y = centy[k][my-2][i].y + gs[k][my+1][i].y;
731 centy[k][my-1][i].z = centy[k][my-2][i].z;
732 }
733 }
734 }
735 }
736
737 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
738 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
739
740 PROFILE_FUNCTION_END;
741 PetscFunctionReturn(0);
742}
UserCtx::Centy
Vec Centy
Definition variables.h:859
Here is the caller graph for this function:

◆ ApplyPeriodicCorrectionsToKFaceCenter()

PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter ( UserCtx *  user)

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

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

Parameters
userThe UserCtx containing grid and field data.
Returns
PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 751 of file Metric.c.

752{
753 PetscErrorCode ierr;
754 DMDALocalInfo info = user->info;
755 PetscInt zs = info.zs, ze = info.zs + info.zm;
756 PetscInt mz = info.mz;
757 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
758 PetscInt gys = info.gys, gye = info.gys + info.gym;
759 Cmpnts ***centz, ***gs;
760
761 PetscFunctionBeginUser;
762 PROFILE_FUNCTION_BEGIN;
763
764 // Check if Z-periodic boundaries exist
765 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC &&
766 user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
767 LOG_ALLOW(LOCAL, LOG_TRACE, "No Z-periodic boundaries; skipping Centz corrections.\n");
768 PROFILE_FUNCTION_END;
769 PetscFunctionReturn(0);
770 }
771
772 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying Z-periodic corrections to Centz.\n");
773
774 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
775 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
776
777 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
778 if (user->cgrid) {
779 for (PetscInt j=gys+1; j<gye; j++) {
780 for (PetscInt i=gxs+1; i<gxe; i++) {
781 centz[-1][j][i].x = centz[-3][j][i].x;
782 centz[-1][j][i].y = centz[-3][j][i].y;
783 centz[-1][j][i].z = centz[-3][j][i].z;
784 }
785 }
786 } else {
787 for (PetscInt j=gys+1; j<gye; j++) {
788 for (PetscInt i=gxs+1; i<gxe; i++) {
789 centz[-1][j][i].x = centz[0][j][i].x;
790 centz[-1][j][i].y = centz[0][j][i].y;
791 centz[-1][j][i].z = centz[0][j][i].z - gs[-2][j][i].z;
792 }
793 }
794 }
795 }
796
797 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
798 if (user->cgrid) {
799 for (PetscInt j=gys+1; j<gye; j++) {
800 for (PetscInt i=gxs+1; i<gxe; i++) {
801 centz[mz-1][j][i].x = centz[mz+1][j][i].x;
802 centz[mz-1][j][i].y = centz[mz+1][j][i].y;
803 centz[mz-1][j][i].z = centz[mz+1][j][i].z;
804 }
805 }
806 } else {
807 for (PetscInt j=gys+1; j<gye; j++) {
808 for (PetscInt i=gxs+1; i<gxe; i++) {
809 centz[mz-1][j][i].x = centz[mz-2][j][i].x;
810 centz[mz-1][j][i].y = centz[mz-2][j][i].y;
811 centz[mz-1][j][i].z = centz[mz-2][j][i].z + gs[mz+1][j][i].z;
812 }
813 }
814 }
815 }
816
817 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
818 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
819
820 PROFILE_FUNCTION_END;
821 PetscFunctionReturn(0);
822}
UserCtx::Centz
Vec Centz
Definition variables.h:859
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.

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

Local to this translation unit.

Definition at line 830 of file Metric.c.

831{
832 PetscErrorCode ierr;
833 DMDALocalInfo info;
834 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
835 Cmpnts ***nodal_coords_arr;
836 Vec localCoords_from_dm;
837
838 PetscFunctionBeginUser;
839
840 PROFILE_FUNCTION_BEGIN;
841
842 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
843
844 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
845
846 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
847 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
848 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
849 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
850
851 // --- 2. Get arrays for output global Vecs from UserCtx ---
852 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
853 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
854 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
855
856 // Define owned node ranges (global indices)
857 PetscInt xs = info.xs, xe = info.xs + info.xm;
858 PetscInt ys = info.ys, ye = info.ys + info.ym;
859 PetscInt zs = info.zs, ze = info.zs + info.zm;
860
861 // Global domain dimensions (total number of nodes)
862 PetscInt mx = info.mx;
863 PetscInt my = info.my;
864 PetscInt mz = info.mz;
865
866 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
867 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
868 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
869 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
870 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
871
872 // These represent the surface area of the curvilinear cell face and the normal rotated such that the direction of increasing coordinate is maintained.
873 // The metric vectors (Csi, Eta, Zet) are defined to point in the direction of their corresponding increasing computational coordinate.
874
875 // Calculate Csi
876 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
877 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
878 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
879
880 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);
881 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);
882 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);
883 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);
884 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);
885 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);
886
887 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
888 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
889 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
890 }
891 }
892 }
893
894 // Calculate Eta
895 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
896 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
897 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
898
899 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);
900 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);
901 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);
902 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);
903 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);
904 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);
905
906 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
907 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
908 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
909 }
910 }
911 }
912
913 // Calculate Zet
914 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
915 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
916 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
917
918 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);
919 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);
920 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);
921 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);
922 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);
923 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);
924
925 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
926 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
927 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
928 }
929 }
930 }
931
932 // --- 4. Boundary Extrapolation ---
933 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
934 PetscInt i_bnd, j_bnd, k_bnd;
935
936 if (xs == 0) { // If this rank owns the global i=0 boundary
937 i_bnd = 0;
938 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
939 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
940 if (i_bnd + 1 < mx) {
941 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
942 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
943 }
944 }
945 }
946 }
947 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
948 i_bnd = mx - 1;
949 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
950 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
951 if (i_bnd - 1 >= 0) {
952 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
953 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
954 }
955 }
956 }
957 }
958 if (ys == 0) {
959 j_bnd = 0;
960 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
961 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
962 if (j_bnd + 1 < my) {
963 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
964 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
965 }
966 }
967 }
968 }
969 if (ye == my) {
970 j_bnd = my - 1;
971 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
972 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
973 if (j_bnd - 1 >= 0) {
974 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
975 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
976 }
977 }
978 }
979 }
980 if (zs == 0) {
981 k_bnd = 0;
982 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
983 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
984 if (k_bnd + 1 < mz) {
985 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
986 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
987 }
988 }
989 }
990 }
991 if (ze == mz) {
992 k_bnd = mz - 1;
993 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
994 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
995 if (k_bnd - 1 >= 0) {
996 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
997 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
998 }
999 }
1000 }
1001 }
1002
1003 if (info.xs==0 && info.ys==0 && info.zs==0) {
1004 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
1005 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
1006 }
1007
1008 // --- 5. Restore all arrays ---
1009 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1010 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
1011 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
1012 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
1013
1014 // --- 6. Assemble Global Vectors ---
1015 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
1016 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
1017 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
1018 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
1019
1020 // --- 7. Update Local Ghosted Versions ---
1021 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
1022 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
1023 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
1024 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
1025
1026 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
1027
1028 PROFILE_FUNCTION_END;
1029
1030 PetscFunctionReturn(0);
1031}
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
UserCtx::Zet
Vec Zet
Definition variables.h:858
UserCtx::Csi
Vec Csi
Definition variables.h:858
UserCtx::Eta
Vec Eta
Definition variables.h:858
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeCellCenteredJacobianInverse()

PetscErrorCode ComputeCellCenteredJacobianInverse ( UserCtx *  user)

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

This version includes boundary extrapolation.

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

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

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

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

See also
ComputeCellCenteredJacobianInverse()

Definition at line 1042 of file Metric.c.

1043{
1044 PetscErrorCode ierr;
1045 DMDALocalInfo info;
1046 PetscScalar ***aj_arr;
1047 Cmpnts ***nodal_coords_arr;
1048 Vec localCoords_from_dm;
1049
1050 PetscFunctionBeginUser;
1051 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
1052
1053 // --- 1. Get Nodal Coordinates and Output Array ---
1054 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
1055 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1056 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1057 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
1058
1059 // Define owned node ranges (global indices)
1060 PetscInt xs = info.xs, xe = info.xs + info.xm;
1061 PetscInt ys = info.ys, ye = info.ys + info.ym;
1062 PetscInt zs = info.zs, ze = info.zs + info.zm;
1063
1064 // Global domain dimensions (total number of nodes)
1065 PetscInt mx = info.mx;
1066 PetscInt my = info.my;
1067 PetscInt mz = info.mz;
1068
1069 // --- 2. Calculate Aj for INTERIOR Stencils ---
1070
1071 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
1072 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
1073 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
1074
1075 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
1076 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
1077 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
1078
1079 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
1080 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
1081 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
1082
1083 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) );
1084
1085 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) );
1086
1087 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) );
1088
1089 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) );
1090
1091 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) );
1092
1093 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) );
1094
1095 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) );
1096
1097 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) );
1098
1099 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) );
1100
1101 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);
1102 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
1103 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
1104 }
1105 }
1106 }
1107
1108 // --- 4. Boundary Extrapolation for Aj ---
1109 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
1110 PetscInt i_bnd, j_bnd, k_bnd;
1111
1112 if (xs == 0) {
1113 i_bnd = 0;
1114 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1115 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1116 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
1117 }
1118 }
1119 }
1120 if (xe == mx) {
1121 i_bnd = mx - 1;
1122 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1123 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1124 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
1125 }
1126 }
1127 }
1128 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
1129 if (ys == 0) {
1130 j_bnd = 0;
1131 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1132 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1133 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
1134 }
1135 }
1136 }
1137 if (ye == my) {
1138 j_bnd = my - 1;
1139 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1140 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1141 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
1142 }
1143 }
1144 }
1145 if (zs == 0) {
1146 k_bnd = 0;
1147 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1148 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1149 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
1150 }
1151 }
1152 }
1153 if (ze == mz) {
1154 k_bnd = mz - 1;
1155 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1156 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1157 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
1158 }
1159 }
1160 }
1161
1162 // --- 5. Restore arrays ---
1163 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1164 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
1165
1166 // --- 6. Assemble Global Vector ---
1167 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
1168 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
1169 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
1170
1171 // --- 7. Update Local Ghosted Version ---
1172 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
1173 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
1174
1175 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
1176 PetscFunctionReturn(0);
1177}
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 metadata intentionally omitted in API docs.

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

Local to this translation unit.

Definition at line 314 of file Metric.c.

315{
316 PetscErrorCode ierr;
317 PetscReal aj_min, aj_max;
318 PetscMPIInt rank;
319
320 PetscFunctionBeginUser;
321
322 PROFILE_FUNCTION_BEGIN;
323
324 /* ---------------- step 1: global extrema of Aj ---------------- */
325 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
326
327 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr); /* already global */
328 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
329
330 LOG_ALLOW(GLOBAL, LOG_INFO,
331 "[orientation] Global Aj range: [%.3e , %.3e]\n",
332 (double)aj_min, (double)aj_max);
333
334 /* ---------------- step 2: detect malformed mesh ---------------- */
335 if (aj_min < 0.0 && aj_max > 0.0)
336 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
337 "Mixed Jacobian signs detected – grid is topologically inconsistent.");
338
339 /* Default: grid is right-handed unless proven otherwise */
340 PetscInt orientation = +1;
341
342 /* ---------------- step 3: repair left-handed mesh -------------- */
343 if (aj_max < 0.0) { /* entire domain has Aj < 0 */
344 orientation = -1;
345
346 if (!rank)
347 LOG_ALLOW(LOCAL, LOG_INFO,
348 "[orientation] Detected left-handed grid – flipping metric vectors\n");
349
350 /* Flip sign of *all* metric vectors and Aj */
351 ierr = VecScale(user->Csi, -1.0); CHKERRQ(ierr);
352 ierr = VecScale(user->Eta, -1.0); CHKERRQ(ierr);
353 ierr = VecScale(user->Zet, -1.0); CHKERRQ(ierr);
354 ierr = VecScale(user->Aj , -1.0); CHKERRQ(ierr);
355
356 /* Local ghost regions now stale – refresh */
357 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
358 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
359 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
360 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
361
362 /* Sanity print: Aj must be > 0 now */
363 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
364 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
365
366 if (aj_min <= 0.0)
367 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
368 "Failed to flip grid orientation – Aj still non-positive.");
369 else if (aj_min && aj_max > 0.0)
370 orientation = +1;
371 }
372
373 /* ---------------- step 4: store result in UserCtx -------------- */
374 user->GridOrientation = orientation;
375
376 if (!rank)
377 LOG_ALLOW(LOCAL, LOG_INFO,
378 "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
379 (orientation>0) ? "right" : "left", orientation);
380
381 PROFILE_FUNCTION_END;
382
383 PetscFunctionReturn(0);
384}
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:824
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.

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

Local to this translation unit.

Definition at line 1186 of file Metric.c.

1187{
1188 PetscErrorCode ierr;
1189 DMDALocalInfo info;
1190 Vec lCoords;
1191 const Cmpnts ***coor;
1192 Cmpnts ***cent, ***gs;
1193 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
1194 PetscInt xs,ys,zs,xe,ye,ze,mx,my,mz;
1195
1196 PetscFunctionBeginUser;
1197
1198 PROFILE_FUNCTION_BEGIN;
1199
1200 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);
1201
1202 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1203 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1204 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1205
1206 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1207 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1208
1209 xs = info.xs; xe = info.xs + info.xm;
1210 ys = info.ys; ye = info.ys + info.ym;
1211 zs = info.zs; ze = info.zs + info.zm;
1212 mx = info.mx; my = info.my; mz = info.mz;
1213
1214 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
1215 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
1216 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
1217
1218 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
1219 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
1220 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
1221
1222 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
1223 for (PetscInt k=k_start_node; k<k_end_node; k++) {
1224 for (PetscInt j=j_start_node; j<j_end_node; j++) {
1225 for (PetscInt i=i_start_node; i<i_end_node; i++) {
1226 // Calculate cell center as the average of its 8 corner nodes
1227 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);
1228 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);
1229 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);
1230
1231 // Calculate Grid Spacing in i-direction (distance between i-face centers)
1232 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);
1233 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);
1234 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);
1235 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);
1236 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);
1237 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);
1238 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1239
1240 // Calculate Grid Spacing in j-direction (distance between j-face centers)
1241 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);
1242 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);
1243 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);
1244 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);
1245 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);
1246 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);
1247 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1248
1249 // Calculate Grid Spacing in k-direction (distance between k-face centers)
1250 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);
1251 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);
1252 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);
1253 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);
1254 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);
1255 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);
1256 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1257 }
1258 }
1259 }
1260
1261 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1262 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1263 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1264
1265 // Assemble and update ghost regions for the new data
1266 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1267 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
1268 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
1269 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
1270
1271 ierr = ApplyPeriodicCorrectionsToCellCentersAndSpacing(user); CHKERRQ(ierr);
1272
1273 // Final assembly and ghost update after corrections
1274 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr);
1275 ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1276 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
1277
1278 PROFILE_FUNCTION_END;
1279
1280 PetscFunctionReturn(0);
1281}
ApplyPeriodicCorrectionsToCellCentersAndSpacing
PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToCellCentersAndSpacing().
Definition Metric.c:392
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:858
SimCtx::rank
PetscMPIInt rank
Definition variables.h:646
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
UserCtx::_this
PetscInt _this
Definition variables.h:824
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:874
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeIFaceMetrics()

PetscErrorCode ComputeIFaceMetrics ( UserCtx *  user)

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

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

Parameters
userThe UserCtx for a specific grid level. Populates user->ICsi, etc.
Returns
PetscErrorCode 0 on success, or a PETSc error code on failure.

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

Local to this translation unit.

Definition at line 1289 of file Metric.c.

1290{
1291 PetscErrorCode ierr;
1292 DMDALocalInfo info;
1293 Vec lCoords;
1294 const Cmpnts ***coor;
1295 Cmpnts ***centx; //***gs;
1296 const Cmpnts ***centx_const;
1297 Cmpnts ***icsi, ***ieta, ***izet;
1298 PetscScalar ***iaj;
1299 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1300
1301 PetscFunctionBeginUser;
1302
1303 PROFILE_FUNCTION_BEGIN;
1304
1305 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1306
1307 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1308 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1309 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1310 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1311 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1312 PetscInt gys = info.gys, gye = info.gys + info.gym;
1313 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1314
1315 PetscInt lxe = xe;
1316 PetscInt lys = ys; PetscInt lye = ye;
1317 PetscInt lzs = zs; PetscInt lze = ze;
1318
1319 if (ys==0) lys = ys+1;
1320 if (zs==0) lzs = zs+1;
1321
1322 if (xe==mx) lxe=xe-1;
1323 if (ye==my) lye=ye-1;
1324 if (ze==mz) lze=ze-1;
1325
1326 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1327 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1328 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1329 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1330 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1331
1332 //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);
1333
1334 // Loop over the ghosted region to calculate all local face centers
1335 // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
1336 PetscInt j_end = (ye == mx)? my - 1:gye;
1337 PetscInt k_end = (ze == mz)? mz - 1:gze;
1338
1339 for (PetscInt k = gzs + 1; k < k_end; k++) {
1340 for (PetscInt j = gys + 1; j < j_end; j++) {
1341 for (PetscInt i = gxs; i < gxe; i++) {
1342 //----- DEBUG ------
1343 //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);
1344 //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,
1345 // coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
1346 // coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
1347 // coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
1348 // coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
1349
1350 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);
1351 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);
1352 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);
1353
1354 //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);
1355 }
1356 }
1357 }
1358
1359 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
1360 /*
1361 if(xs==0){
1362 for(PetscInt k=gzs+1;k < gze; k++){
1363 for(PetscInt j=gys+1;j < gye; j++){
1364 PetscInt i=0;
1365 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
1366 centx[k][j][i-1].y=centx[k][j][i].y;
1367 centx[k][j][i-1].z=centx[k][j][i].z;
1368 }
1369 }
1370 }
1371 if (xe==mx){
1372 for(PetscInt k=gzs+1; k<gze; k++) {
1373 for (PetscInt j=gys+1; j<gye;j++) {
1374 PetscInt i=mx-1;
1375 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
1376 centx[k][j][i].y=centx[k][j][i-1].y;
1377 centx[k][j][i].z=centx[k][j][i-1].z;
1378 }
1379 }
1380 }
1381 */
1382
1383 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1384 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1385
1386 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1387
1388 // For periodic BCs, exchange ghosts between processes
1389 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
1390 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
1391 ierr = DMLocalToLocalBegin(user->fda, user->Centx, INSERT_VALUES, user->Centx); CHKERRQ(ierr);
1392 ierr = DMLocalToLocalEnd(user->fda, user->Centx, INSERT_VALUES, user->Centx); CHKERRQ(ierr);
1393 }
1394
1395 ierr = ApplyPeriodicCorrectionsToIFaceCenter(user); CHKERRQ(ierr);
1396
1397 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1398
1399 // --- Part 2: Calculate metrics using face-centered coordinates ---
1400 ierr = DMDAVecGetArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1401 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1402 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1403 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1404 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1405
1406 // Loop over the OWNED region where we will store the final metrics
1407 for (PetscInt k=lzs; k<lze; k++) {
1408 for (PetscInt j=lys; j<lye; j++) {
1409 for (PetscInt i=xs; i<lxe; i++) {
1410
1411 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1412 if(i == 0 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC){
1413 // Forward difference at the domain's min-i boundary
1414 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
1415 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
1416 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
1417 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1418 // Backward difference at the domain's max-i boundary
1419 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
1420 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
1421 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
1422 } else { // Central difference in the interior (or if PERIODIC BCs)
1423 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
1424 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
1425 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
1426 }
1427
1428 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1429 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1430 // Forward difference
1431 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
1432 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
1433 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
1434 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1435 // Backward difference
1436 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
1437 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
1438 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
1439 } else { // Central difference (interior or PERIODIC)
1440 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
1441 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
1442 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
1443 }
1444
1445 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1446 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1447 // Forward difference
1448 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
1449 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
1450 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
1451 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1452 // Backward difference
1453 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
1454 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
1455 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
1456 } else { // Central difference (Interior + PERIODIC)
1457 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
1458 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
1459 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
1460 }
1461
1462 // --- Metric calculations (identical to legacy FormMetrics) ---
1463 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1464 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1465 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1466
1467 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1468 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1469 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1470
1471 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
1472 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1473 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
1474
1475 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
1476 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
1477 iaj[k][j][i] = 1.0 / iaj[k][j][i];
1478 }
1479 }
1480 }
1481 }
1482
1483 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1484 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1485 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1486 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1487 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1488
1489 // --- Part 3: Assemble global vectors and update local ghosts ---
1490 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
1491 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
1492 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
1493 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
1494
1495 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
1496 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
1497 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
1498 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
1499
1500 PROFILE_FUNCTION_END;
1501
1502 PetscFunctionReturn(0);
1503}
ApplyPeriodicCorrectionsToIFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToIFaceCenter().
Definition Metric.c:591
UserCtx::IZet
Vec IZet
Definition variables.h:860
UserCtx::IEta
Vec IEta
Definition variables.h:860
UserCtx::ICsi
Vec ICsi
Definition variables.h:860
UserCtx::IAj
Vec IAj
Definition variables.h:860
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.

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

Local to this translation unit.

Definition at line 1511 of file Metric.c.

1512{
1513 PetscErrorCode ierr;
1514 DMDALocalInfo info;
1515 Vec lCoords;
1516 const Cmpnts ***coor;
1517 Cmpnts ***centy; //***gs;
1518 const Cmpnts ***centy_const;
1519 Cmpnts ***jcsi, ***jeta, ***jzet;
1520 PetscScalar ***jaj;
1521 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1522
1523 PetscFunctionBeginUser;
1524
1525 PROFILE_FUNCTION_BEGIN;
1526
1527 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1528
1529 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1530 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1531 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1532 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1533 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1534 PetscInt gys = info.gys, gye = info.gys + info.gym;
1535 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1536
1537 PetscInt lxs = xs; PetscInt lxe = xe;
1538 PetscInt lye = ye;
1539 PetscInt lzs = zs; PetscInt lze = ze;
1540
1541 if (xs==0) lxs = xs+1;
1542 if (zs==0) lzs = zs+1;
1543
1544 if (xe==mx) lxe=xe-1;
1545 if (ye==my) lye=ye-1;
1546 if (ze==mz) lze=ze-1;
1547
1548 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1549 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1550 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1551 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1552 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1553
1554 // Loop over the ghosted region to calculate all local face centers
1555 // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
1556 PetscInt k_end = (ze == mz)? mz - 1:gze;
1557 PetscInt i_end = (xe == mx)? mx - 1:gxe;
1558
1559 for (PetscInt k = gzs + 1; k < k_end; k++) {
1560 for (PetscInt j = gys; j < gye; j++) {
1561 for (PetscInt i = gxs + 1; i < i_end; i++) {
1562 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);
1563 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);
1564 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);
1565 }
1566 }
1567 }
1568
1569 /*
1570 if(ys==0){
1571 for(PetscInt k=gzs+1;k < gze; k++){
1572 for(PetscInt i=gxs+1;j < gxe; i++){
1573 PetscInt j=0;
1574 centy[k][j-1][i].x=centy[k][j][i].x;
1575 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1576 centy[k][j-1][i].z=centy[k][j][i].z;
1577 }
1578 }
1579 }
1580 if (ye==my){
1581 for(PetscInt k=gzs+1; k<gze; k++) {
1582 for (PetscInt i=gxs+1; j<gxe;i++) {
1583 PetscInt j=my-1;
1584 centy[k][j][i].x=centy[k][j-1][i].x
1585 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1586 centy[k][j][i].z=centy[k][j-1][i].z;
1587 }
1588 }
1589 }
1590 */
1591
1592 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1593 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1594 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1595
1596 // For periodic BCs, exchange ghosts between processes
1597 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
1598 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
1599 ierr = DMLocalToLocalBegin(user->fda, user->Centy, INSERT_VALUES, user->Centy); CHKERRQ(ierr);
1600 ierr = DMLocalToLocalEnd(user->fda, user->Centy, INSERT_VALUES, user->Centy); CHKERRQ(ierr);
1601 }
1602
1603 ierr = ApplyPeriodicCorrectionsToJFaceCenter(user); CHKERRQ(ierr);
1604
1605 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1606
1607 // --- Part 2: Calculate metrics using face-centered coordinates ---
1608 ierr = DMDAVecGetArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1609 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1610 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1611 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1612 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1613
1614 // Loop over the OWNED region where we will store the final metrics
1615 for (PetscInt k=lzs; k<lze; k++) {
1616 for (PetscInt j=ys; j<lye; j++) {
1617 for (PetscInt i=lxs; i<lxe; i++) {
1618
1619 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1620 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1621 // Forward difference at the domain's min-i boundary
1622 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1623 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1624 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1625 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1626 // Backward difference at the domain's max-i boundary
1627 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1628 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1629 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1630 } else { // Central difference in the interior or PERIODIC
1631 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1632 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1633 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1634 }
1635
1636 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1637 if (j == 0 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1638 // Forward difference
1639 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1640 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1641 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1642 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1643 // Backward difference
1644 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1645 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1646 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1647 } else { // Central difference (interior or PERIODIC)
1648 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1649 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1650 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1651 }
1652
1653 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1654 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1655 // Forward difference
1656 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1657 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1658 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1659 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1660 // Backward difference
1661 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1662 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1663 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1664 } else { // Central difference (Interior or PERIODIC)
1665 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1666 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1667 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1668 }
1669
1670 // --- Metric calculations (identical to legacy FormMetrics) ---
1671 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1672 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1673 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1674
1675 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1676 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1677 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1678
1679 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1680 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1681 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1682
1683 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1684 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1685 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1686 }
1687 }
1688 }
1689 }
1690
1691 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1692 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1693 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1694 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1695 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1696
1697 // --- Part 3: Assemble global vectors and update local ghosts ---
1698 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1699 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1700 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1701 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1702
1703 ierr = UpdateLocalGhosts(user, "JCsi"); CHKERRQ(ierr);
1704 ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
1705 ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
1706 ierr = UpdateLocalGhosts(user, "JAj"); CHKERRQ(ierr);
1707
1708 PROFILE_FUNCTION_END;
1709
1710 PetscFunctionReturn(0);
1711}
ApplyPeriodicCorrectionsToJFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToJFaceCenter().
Definition Metric.c:671
UserCtx::JCsi
Vec JCsi
Definition variables.h:861
UserCtx::JEta
Vec JEta
Definition variables.h:861
UserCtx::JZet
Vec JZet
Definition variables.h:861
UserCtx::JAj
Vec JAj
Definition variables.h:861
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeKFaceMetrics()

PetscErrorCode ComputeKFaceMetrics ( UserCtx *  user)

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

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

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

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

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

Local to this translation unit.

Definition at line 1719 of file Metric.c.

1720{
1721 PetscErrorCode ierr;
1722 DMDALocalInfo info;
1723 Vec lCoords;
1724 const Cmpnts ***coor;
1725 Cmpnts ***centz; //***gs;
1726 const Cmpnts ***centz_const;
1727 Cmpnts ***kcsi, ***keta, ***kzet;
1728 PetscScalar ***kaj;
1729 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1730
1731 PetscFunctionBeginUser;
1732
1733 PROFILE_FUNCTION_BEGIN;
1734
1735 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1736
1737 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1738 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1739 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1740 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1741 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1742 PetscInt gys = info.gys, gye = info.gys + info.gym;
1743 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1744
1745 PetscInt lxs = xs; PetscInt lxe = xe;
1746 PetscInt lys = ys; PetscInt lye = ye;
1747 PetscInt lze = ze;
1748
1749 if (xs==0) lxs = xs+1;
1750 if (ys==0) lys = ys+1;
1751
1752 if (xe==mx) lxe=xe-1;
1753 if (ye==my) lye=ye-1;
1754 if (ze==mz) lze=ze-1;
1755
1756 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1757 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1758 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1759 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1760 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1761
1762 // Loop over the ghosted region to calculate all local face centers
1763 // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
1764 for (PetscInt k = gzs; k < gze; k++) {
1765 for (PetscInt j = gys; j < gye; j++) {
1766 for (PetscInt i = gxs + 1; i < gxe; i++) {
1767 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);
1768 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);
1769 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);
1770 }
1771 }
1772 }
1773
1774 /*
1775 if(zs==0){
1776 for(PetscInt j=gys+1;j < gye; j++){
1777 for(PetscInt i=gxs+1;j < gxe; i++){
1778 PetscInt k=0;
1779 centz[k-1][j][i].x=centz[k][j][i].x;
1780 centz[k-1][j][i].y=centz[k][j][i].y;
1781 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1782 }
1783 }
1784 }
1785 if (ze==mz){
1786 for(PetscInt j=gys+1; j<gye; j++) {
1787 for (PetscInt i=gxs+1; j<gxe;i++) {
1788 PetscInt k=mz-1;
1789 centy[k][j][i].x=centy[k-1][j][i].x
1790 centy[k][j][i].y=centy[k-1][j][i].y;
1791 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1792 }
1793 }
1794 }
1795 */
1796
1797 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1798 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1799 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1800
1801 // For periodic BCs, exchange ghosts between processes
1802 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
1803 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
1804 ierr = DMLocalToLocalBegin(user->fda, user->Centz, INSERT_VALUES, user->Centz); CHKERRQ(ierr);
1805 ierr = DMLocalToLocalEnd(user->fda, user->Centz, INSERT_VALUES, user->Centz); CHKERRQ(ierr);
1806 }
1807
1808 ierr = ApplyPeriodicCorrectionsToKFaceCenter(user); CHKERRQ(ierr);
1809
1810 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1811
1812 // --- Part 2: Calculate metrics using face-centered coordinates ---
1813 ierr = DMDAVecGetArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1814 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1815 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1816 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1817 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1818
1819 // Loop over the OWNED region where we will store the final metrics
1820 for (PetscInt k=zs; k<lze; k++) {
1821 for (PetscInt j=lys; j<lye; j++) {
1822 for (PetscInt i=lxs; i<lxe; i++) {
1823
1824 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1825 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1826 // Forward difference at the domain's min-i boundary
1827 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
1828 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
1829 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
1830 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1831 // Backward difference at the domain's max-i boundary
1832 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
1833 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
1834 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
1835 } else { // Central difference in the interior (or PERIODIC)
1836 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
1837 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
1838 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
1839 }
1840
1841 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1842 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1843 // Forward difference
1844 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
1845 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
1846 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
1847 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1848 // Backward difference
1849 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
1850 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
1851 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
1852 } else { // Central difference (interior or PERIODIC)
1853 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
1854 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
1855 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
1856 }
1857
1858 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1859 if (k == 0 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1860 // Forward difference
1861 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
1862 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
1863 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
1864 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1865 // Backward difference
1866 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
1867 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
1868 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
1869 } else { // Central difference (Interior or PERIODIC)
1870 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
1871 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
1872 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
1873 }
1874
1875 // --- Metric calculations (identical to legacy FormMetrics) ---
1876 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1877 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1878 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1879
1880 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1881 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1882 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1883
1884 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1885 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1886 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1887
1888 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
1889 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
1890 kaj[k][j][i] = 1.0 / kaj[k][j][i];
1891 }
1892 }
1893 }
1894 }
1895
1896 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1897 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1898 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1899 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1900 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1901
1902 // --- Part 3: Assemble global vectors and update local ghosts ---
1903 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
1904 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
1905 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
1906 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
1907
1908 ierr = UpdateLocalGhosts(user, "KCsi"); CHKERRQ(ierr);
1909 ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
1910 ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
1911 ierr = UpdateLocalGhosts(user, "KAj"); CHKERRQ(ierr);
1912
1913 PROFILE_FUNCTION_END;
1914
1915 PetscFunctionReturn(0);
1916}
ApplyPeriodicCorrectionsToKFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToKFaceCenter().
Definition Metric.c:751
UserCtx::KAj
Vec KAj
Definition variables.h:862
UserCtx::KEta
Vec KEta
Definition variables.h:862
UserCtx::KZet
Vec KZet
Definition variables.h:862
UserCtx::KCsi
Vec KCsi
Definition variables.h:862
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeMetricsDivergence()

PetscErrorCode ComputeMetricsDivergence ( UserCtx *  user)

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

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

Parameters
userThe UserCtx for a specific grid level (typically the finest).
Returns
PetscErrorCode 0 on success, or a PETSc error code on failure.

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

Local to this translation unit.

Definition at line 1944 of file Metric.c.

1945{
1946 DM da = user->da, fda = user->fda;
1947 DMDALocalInfo info = user->info;
1948 PetscInt xs = info.xs, xe = info.xs + info.xm;
1949 PetscInt ys = info.ys, ye = info.ys + info.ym;
1950 PetscInt zs = info.zs, ze = info.zs + info.zm;
1951 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1952 PetscInt lxs, lys, lzs, lxe, lye, lze;
1953 PetscInt i, j, k;
1954 Vec Div;
1955 PetscReal ***div, ***aj;
1956 Cmpnts ***csi, ***eta, ***zet;
1957 PetscReal maxdiv;
1958
1959 PetscFunctionBeginUser;
1960
1961 PROFILE_FUNCTION_BEGIN;
1962
1963 lxs = xs; lxe = xe;
1964 lys = ys; lye = ye;
1965 lzs = zs; lze = ze;
1966
1967 if (xs == 0) lxs = xs + 1;
1968 if (ys == 0) lys = ys + 1;
1969 if (zs == 0) lzs = zs + 1;
1970
1971 if (xe == mx) lxe = xe - 1;
1972 if (ye == my) lye = ye - 1;
1973 if (ze == mz) lze = ze - 1;
1974
1975 DMDAVecGetArray(fda, user->lCsi, &csi);
1976 DMDAVecGetArray(fda, user->lEta, &eta);
1977 DMDAVecGetArray(fda, user->lZet, &zet);
1978 DMDAVecGetArray(da, user->lAj, &aj);
1979
1980 VecDuplicate(user->P, &Div);
1981 VecSet(Div, 0.);
1982 DMDAVecGetArray(da, Div, &div);
1983
1984 for (k = lzs; k < lze; k++) {
1985 for (j = lys; j < lye; j++) {
1986 for (i = lxs; i < lxe; i++) {
1987 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1988 eta[k][j][i].x - eta[k][j-1][i].x +
1989 zet[k][j][i].x - zet[k-1][j][i].x +
1990 csi[k][j][i].y - csi[k][j][i-1].y +
1991 eta[k][j][i].y - eta[k][j-1][i].y +
1992 zet[k][j][i].y - zet[k-1][j][i].y +
1993 csi[k][j][i].z - csi[k][j][i-1].z +
1994 eta[k][j][i].z - eta[k][j-1][i].z +
1995 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1996 div[k][j][i] = fabs(divergence);
1997 }
1998 }
1999 }
2000
2001 DMDAVecRestoreArray(da, Div, &div);
2002
2003 PetscInt MaxFlatIndex = -1;
2004 VecMax(Div, &MaxFlatIndex, &maxdiv);
2005 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %" PetscInt_FMT ".\n",maxdiv,MaxFlatIndex);
2006
2007 for (k=zs; k<ze; k++) {
2008 for (j=ys; j<ye; j++) {
2009 for (i=xs; i<xe; i++) {
2010 if (Gidx(i,j,k,user) == MaxFlatIndex) {
2011 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,(int)k,(int)j,(int)i);
2012 }
2013 }
2014 }
2015 }
2016
2017
2018 DMDAVecRestoreArray(fda, user->lCsi, &csi);
2019 DMDAVecRestoreArray(fda, user->lEta, &eta);
2020 DMDAVecRestoreArray(fda, user->lZet, &zet);
2021 DMDAVecRestoreArray(da, user->lAj, &aj);
2022 VecDestroy(&Div);
2023
2024
2025 PROFILE_FUNCTION_END;
2026
2027 PetscFunctionReturn(0);
2028}
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Internal helper implementation: Gidx().
Definition Metric.c:1922
UserCtx::lZet
Vec lZet
Definition variables.h:858
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
UserCtx::lAj
Vec lAj
Definition variables.h:858
UserCtx::lEta
Vec lEta
Definition variables.h:858
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

Computes the max-min values of the grid metrics.

Local to this translation unit.

Definition at line 2036 of file Metric.c.

2037{
2038
2039 DMDALocalInfo info = user->info;
2040 PetscInt xs = info.xs, xe = info.xs + info.xm;
2041 PetscInt ys = info.ys, ye = info.ys + info.ym;
2042 PetscInt zs = info.zs, ze = info.zs + info.zm;
2043 PetscInt i, j, k;
2044
2045 PetscFunctionBeginUser;
2046
2047 PROFILE_FUNCTION_BEGIN;
2048
2049 PetscReal CsiMax, EtaMax, ZetMax;
2050 PetscReal ICsiMax, IEtaMax, IZetMax;
2051 PetscReal JCsiMax, JEtaMax, JZetMax;
2052 PetscReal KCsiMax, KEtaMax, KZetMax;
2053 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
2054
2055 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
2056 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
2057 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
2058 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
2059 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
2060
2061 // Max Values
2062 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
2063 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
2064 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
2065
2066 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
2067 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
2068 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
2069
2070 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
2071 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
2072 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
2073
2074 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
2075 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
2076 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
2077
2078 VecMax(user->lAj,&AjMaxArg,&AjMax);
2079 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
2080 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
2081 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
2082
2083 VecMax(user->lAj,&AjMaxArg,&AjMax);
2084 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
2085 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
2086 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
2087
2088 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
2089
2090 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
2091 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
2092 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
2093 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
2094 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
2095
2096 for (k=zs; k<ze; k++) {
2097 for (j=ys; j<ye; j++) {
2098 for (i=xs; i<xe; i++) {
2099 if (Gidx(i,j,k,user) == CsiMaxArg) {
2100 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
2101 }
2102 if (Gidx(i,j,k,user) == EtaMaxArg) {
2103 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
2104 }
2105 if (Gidx(i,j,k,user) == ZetMaxArg) {
2106 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
2107 }
2108 if (Gidx(i,j,k,user) == ICsiMaxArg) {
2109 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
2110 }
2111 if (Gidx(i,j,k,user) == IEtaMaxArg) {
2112 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
2113 }
2114 if (Gidx(i,j,k,user) == IZetMaxArg) {
2115 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
2116 }
2117 if (Gidx(i,j,k,user) == JCsiMaxArg) {
2118 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
2119 }
2120 if (Gidx(i,j,k,user) == JEtaMaxArg) {
2121 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
2122 }
2123 if (Gidx(i,j,k,user) == JZetMaxArg) {
2124 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
2125 }
2126 if (Gidx(i,j,k,user) == KCsiMaxArg) {
2127 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
2128 }
2129 if (Gidx(i,j,k,user) == KEtaMaxArg) {
2130 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
2131 }
2132 if (Gidx(i,j,k,user) == KZetMaxArg) {
2133 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
2134 }
2135 if (Gidx(i,j,k,user) == AjMaxArg) {
2136 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
2137 }
2138 if (Gidx(i,j,k,user) == IAjMaxArg) {
2139 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
2140 }
2141 if (Gidx(i,j,k,user) == JAjMaxArg) {
2142 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
2143 }
2144 if (Gidx(i,j,k,user) == KAjMaxArg) {
2145 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
2146 }
2147 }
2148 }
2149 }
2150
2151 /*
2152 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
2153 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
2154 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
2155 */
2156
2157 PROFILE_FUNCTION_END;
2158
2159 PetscFunctionReturn(0);
2160}
UserCtx::lIEta
Vec lIEta
Definition variables.h:860
UserCtx::lIZet
Vec lIZet
Definition variables.h:860
UserCtx::lIAj
Vec lIAj
Definition variables.h:860
UserCtx::lKEta
Vec lKEta
Definition variables.h:862
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:861
UserCtx::lKZet
Vec lKZet
Definition variables.h:862
UserCtx::lJEta
Vec lJEta
Definition variables.h:861
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:862
UserCtx::lJZet
Vec lJZet
Definition variables.h:861
UserCtx::lICsi
Vec lICsi
Definition variables.h:860
UserCtx::lJAj
Vec lJAj
Definition variables.h:861
UserCtx::lKAj
Vec lKAj
Definition variables.h:862
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

Orchestrates the calculation of all grid metrics.

Local to this translation unit.

Definition at line 2168 of file Metric.c.

2169{
2170 PetscErrorCode ierr;
2171 UserMG *usermg = &simCtx->usermg;
2172 MGCtx *mgctx = usermg->mgctx;
2173 PetscInt nblk = simCtx->block_number;
2174
2175 PetscFunctionBeginUser;
2176
2177 PROFILE_FUNCTION_BEGIN;
2178
2179 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
2180
2181 // Loop through all levels and all blocks
2182 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
2183 for (PetscInt bi = 0; bi < nblk; bi++) {
2184 UserCtx *user = &mgctx[level].user[bi];
2185 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
2186
2187 // Call the modern, modular helper functions for each UserCtx.
2188 // These functions are self-contained and operate on the data within the provided context.
2189 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
2190 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
2191 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
2192 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
2193 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
2194 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
2195 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
2196
2197 // Apply Periodic Boundary Condition Adjustments if necessary
2198 ierr = ApplyMetricsPeriodicBCs(user); CHKERRQ(ierr);
2199 // Diagnostics
2200 ierr = ComputeMetricNorms(user);
2201 if (level == usermg->mglevels - 1) {
2202 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
2203 }
2204 }
2205 }
2206
2207 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
2208
2209 PROFILE_FUNCTION_END;
2210
2211 PetscFunctionReturn(0);
2212}
ApplyMetricsPeriodicBCs
PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
(Orchestrator) Updates all metric-related fields in the local ghost cell regions for periodic boundar...
Definition Boundaries.c:1895
ComputeMetricNorms
PetscErrorCode ComputeMetricNorms(UserCtx *user)
Internal helper implementation: ComputeMetricNorms().
Definition Metric.c:2036
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Internal helper implementation: CheckAndFixGridOrientation().
Definition Metric.c:314
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Internal helper implementation: ComputeCellCentersAndSpacing().
Definition Metric.c:1186
ComputeJFaceMetrics
PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeJFaceMetrics().
Definition Metric.c:1511
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Internal helper implementation: ComputeMetricsDivergence().
Definition Metric.c:1944
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeFaceMetrics().
Definition Metric.c:830
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Implementation of ComputeCellCenteredJacobianInverse().
Definition Metric.c:1042
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeKFaceMetrics().
Definition Metric.c:1719
ComputeIFaceMetrics
PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeIFaceMetrics().
Definition Metric.c:1289
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:252
MGCtx::user
UserCtx * user
Definition variables.h:528
SimCtx::block_number
PetscInt block_number
Definition variables.h:712
SimCtx::usermg
UserMG usermg
Definition variables.h:764
UserMG::mglevels
PetscInt mglevels
Definition variables.h:535
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:538
MGCtx
Context for Multigrid operations.
Definition variables.h:527
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:811
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:534
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8