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 wrapper: map (cell index, ξ,η,ζ) to (x,y,z).
 
PetscErrorCode MetricGetCellVertices (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
 Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
 
PetscErrorCode MetricJacobian (UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, PetscReal J[3][3], PetscReal *detJ)
 Compute Jacobian matrix and its determinant at (xi,eta,zta).
 
PetscErrorCode MetricVelocityContravariant (const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
 Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).
 
PetscErrorCode ComputeFaceMetrics (UserCtx *user)
 Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation, and stores them in the corresponding global Vec members of the UserCtx structure (user->Csi, user->Eta, user->Zet).
 
PetscErrorCode ComputeCellCenteredJacobianInverse (UserCtx *user)
 Calculates the cell-centered inverse Jacobian determinant (1/J) for INTERIOR cells and stores it in user->Aj.
 
PetscErrorCode CheckAndFixGridOrientation (UserCtx *user)
 Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain.
 
PetscErrorCode ComputeCellCentersAndSpacing (UserCtx *user)
 Computes the physical location of cell centers and the spacing between them.
 
PetscErrorCode ComputeIFaceMetrics (UserCtx *user)
 Computes metrics centered on constant-i faces (i-faces).
 
PetscErrorCode ComputeJFaceMetrics (UserCtx *user)
 Computes metrics centered on constant-j faces (j-faces).
 
PetscErrorCode ComputeKFaceMetrics (UserCtx *user)
 Computes metrics centered on constant-k faces (k-faces).
 
PetscErrorCode ComputeMetricsDivergence (UserCtx *user)
 Performs a diagnostic check on the divergence of the face area metric vectors.
 
PetscErrorCode ComputeMetricNorms (UserCtx *user)
 Computes the max-min values of the grid metrics.
 
PetscErrorCode CalculateAllGridMetrics (SimCtx *simCtx)
 Orchestrates the calculation of all grid metrics.
 

Function Documentation

◆ MetricLogicalToPhysical()

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

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

Definition at line 68 of file Metric.c.

73{
74 PetscErrorCode ierr;
75 Cmpnts V[8];
76 PetscFunctionBeginUser;
77 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
78 TrilinearBlend(V,xi,eta,zta,Xp);
79 PetscFunctionReturn(0);
80}
TrilinearBlend
static void TrilinearBlend(const Cmpnts V[8], PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Map logical coordinates to physical space using trilinear basis.
Definition Metric.c:48
MetricGetCellVertices
PetscErrorCode MetricGetCellVertices(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
Definition Metric.c:24
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:99
Here is the call graph for this function:
Here is the caller graph for this function:

◆ MetricGetCellVertices()

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

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

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

Definition at line 24 of file Metric.c.

28{
29 PetscFunctionBeginUser;
30 for (PetscInt c = 0; c < 8; ++c) {
31 PetscInt ii = i + ((c & 1) ? 1 : 0);
32 PetscInt jj = j + ((c & 2) ? 1 : 0);
33 PetscInt kk = k + ((c & 4) ? 1 : 0);
34 LOG_LOOP_ALLOW(GLOBAL, LOG_DEBUG,i+j+k,10," ii: %d,jj:%d,kk:%d - Retrieved.\n",ii,jj,kk);
35 V[c] = X[kk][jj][ii];
36 }
37 PetscFunctionReturn(0);
38}
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:319
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
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 
)

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

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

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

Definition at line 93 of file Metric.c.

98{
99 PetscErrorCode ierr;
100 Cmpnts V[8];
101 PetscFunctionBeginUser;
102 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
103
104 /* derivatives of trilinear shape functions */
105 PetscReal dN_dXi[8], dN_dEta[8], dN_dZta[8];
106 for (PetscInt c=0;c<8;++c) {
107 PetscReal sx = (c & 1) ? 1.0 : -1.0;
108 PetscReal sy = (c & 2) ? 1.0 : -1.0;
109 PetscReal sz = (c & 4) ? 1.0 : -1.0;
110 dN_dXi [c] = 0.125 * sx * ( (c&2?eta:1-eta) ) * ( (c&4?zta:1-zta) );
111 dN_dEta[c] = 0.125 * sy * ( (c&1?xi :1-xi ) ) * ( (c&4?zta:1-zta) );
112 dN_dZta[c] = 0.125 * sz * ( (c&1?xi :1-xi ) ) * ( (c&2?eta:1-eta) );
113 }
114
115 /* assemble Jacobian */
116 PetscReal x_xi=0,y_xi=0,z_xi=0,
117 x_eta=0,y_eta=0,z_eta=0,
118 x_zta=0,y_zta=0,z_zta=0;
119 for (PetscInt c=0;c<8;++c) {
120 x_xi += dN_dXi [c]*V[c].x; y_xi += dN_dXi [c]*V[c].y; z_xi += dN_dXi [c]*V[c].z;
121 x_eta += dN_dEta[c]*V[c].x; y_eta += dN_dEta[c]*V[c].y; z_eta += dN_dEta[c]*V[c].z;
122 x_zta += dN_dZta[c]*V[c].x; y_zta += dN_dZta[c]*V[c].y; z_zta += dN_dZta[c]*V[c].z;
123 }
124
125 J[0][0]=x_xi; J[0][1]=x_eta; J[0][2]=x_zta;
126 J[1][0]=y_xi; J[1][1]=y_eta; J[1][2]=y_zta;
127 J[2][0]=z_xi; J[2][1]=z_eta; J[2][2]=z_zta;
128
129 if (detJ) {
130 *detJ = x_xi*(y_eta*z_zta - y_zta*z_eta)
131 - x_eta*(y_xi*z_zta - y_zta*z_xi)
132 + x_zta*(y_xi*z_eta - y_eta*z_xi);
133 }
134 PetscFunctionReturn(0);
135}
Here is the call graph for this function:

◆ MetricVelocityContravariant()

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

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

Definition at line 141 of file Metric.c.

143{
144 PetscFunctionBeginUser;
145 /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
146 PetscReal gxi[3] = { J[1][1]*J[2][2]-J[1][2]*J[2][1],
147 -J[0][1]*J[2][2]+J[0][2]*J[2][1],
148 J[0][1]*J[1][2]-J[0][2]*J[1][1] };
149 PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
150 J[0][0]*J[2][2]-J[0][2]*J[2][0],
151 -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
152 PetscReal gzta[3] = { J[1][0]*J[2][1]-J[1][1]*J[2][0],
153 -J[0][0]*J[2][1]+J[0][1]*J[2][0],
154 J[0][0]*J[1][1]-J[0][1]*J[1][0] };
155
156 PetscReal invDet = 1.0 / detJ;
157 for (int d=0; d<3; ++d) { gxi[d] *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
158
159 uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
160 uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
161 uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
162 PetscFunctionReturn(0);
163}

◆ ComputeFaceMetrics()

PetscErrorCode ComputeFaceMetrics ( UserCtx *  user)

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

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

  1. Obtains local ghosted nodal coordinates using DMGetCoordinatesLocal.
  2. Calculates metrics for INTERIOR faces where finite difference stencils are valid.
  3. EXTRAPOLATES metrics for faces on the physical domain boundaries by copying from the nearest computed interior face.
  4. Assembles the global user->Csi, user->Eta, user->Zet Vecs.
  5. Updates the local ghosted user->lCsi, user->lEta, user->lZet Vecs.
Parameters
[in,out]userPointer to the UserCtx structure.
Returns
PetscErrorCode 0 on success.
Note
  • This function is a complete "compute and make ready" unit for Csi, Eta, and Zet.
  • It's recommended to call VecZeroEntries on user->Csi, Eta, Zet before this if they might contain old data.

Definition at line 286 of file Metric.c.

287{
288 PetscErrorCode ierr;
289 DMDALocalInfo info;
290 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
291 Cmpnts ***nodal_coords_arr;
292 Vec localCoords_from_dm;
293
294 PetscFunctionBeginUser;
295 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
296
297 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
298
299 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
300 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
301 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
302 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
303
304 // --- 2. Get arrays for output global Vecs from UserCtx ---
305 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
306 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
307 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
308
309 // Define owned node ranges (global indices)
310 PetscInt xs = info.xs, xe = info.xs + info.xm;
311 PetscInt ys = info.ys, ye = info.ys + info.ym;
312 PetscInt zs = info.zs, ze = info.zs + info.zm;
313
314 // Global domain dimensions (total number of nodes)
315 PetscInt mx = info.mx;
316 PetscInt my = info.my;
317 PetscInt mz = info.mz;
318
319 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
320 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
321 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
322 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
323 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
324
325 // Calculate Csi
326 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
327 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
328 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
329
330 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);
331 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);
332 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);
333 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);
334 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);
335 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);
336
337 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
338 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
339 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
340 }
341 }
342 }
343
344 // Calculate Eta
345 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
346 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
347 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
348
349 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);
350 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);
351 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);
352 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);
353 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);
354 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);
355
356 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
357 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
358 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
359 }
360 }
361 }
362
363 // Calculate Zet
364 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
365 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
366 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
367
368 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);
369 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);
370 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);
371 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);
372 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);
373 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);
374
375 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
376 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
377 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
378 }
379 }
380 }
381
382 // --- 4. Boundary Extrapolation ---
383 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
384 PetscInt i_bnd, j_bnd, k_bnd;
385
386 if (xs == 0) { // If this rank owns the global i=0 boundary
387 i_bnd = 0;
388 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
389 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
390 if (i_bnd + 1 < mx) {
391 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
392 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
393 }
394 }
395 }
396 }
397 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
398 i_bnd = mx - 1;
399 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
400 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
401 if (i_bnd - 1 >= 0) {
402 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
403 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
404 }
405 }
406 }
407 }
408 if (ys == 0) {
409 j_bnd = 0;
410 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
411 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
412 if (j_bnd + 1 < my) {
413 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
414 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
415 }
416 }
417 }
418 }
419 if (ye == my) {
420 j_bnd = my - 1;
421 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
422 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
423 if (j_bnd - 1 >= 0) {
424 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
425 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
426 }
427 }
428 }
429 }
430 if (zs == 0) {
431 k_bnd = 0;
432 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
433 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
434 if (k_bnd + 1 < mz) {
435 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
436 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
437 }
438 }
439 }
440 }
441 if (ze == mz) {
442 k_bnd = mz - 1;
443 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
444 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
445 if (k_bnd - 1 >= 0) {
446 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
447 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
448 }
449 }
450 }
451 }
452
453 if (info.xs==0 && info.ys==0 && info.zs==0) {
454 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
455 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
456 }
457
458 // --- 5. Restore all arrays ---
459 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
460 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
461 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
462 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
463
464 // --- 6. Assemble Global Vectors ---
465 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
466 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
467 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
468 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
469
470 // --- 7. Update Local Ghosted Versions ---
471 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
472 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
473 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
474 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
475
476 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
477 PetscFunctionReturn(0);
478}
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:207
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:779
UserCtx::Zet
Vec Zet
Definition variables.h:673
UserCtx::Csi
Vec Csi
Definition variables.h:673
Cmpnts::x
PetscScalar x
Definition variables.h:100
UserCtx::Eta
Vec Eta
Definition variables.h:673
Cmpnts::z
PetscScalar z
Definition variables.h:100
Cmpnts::y
PetscScalar y
Definition variables.h:100
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.

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

Definition at line 491 of file Metric.c.

492{
493 PetscErrorCode ierr;
494 DMDALocalInfo info;
495 PetscScalar ***aj_arr;
496 Cmpnts ***nodal_coords_arr;
497 Vec localCoords_from_dm;
498
499 PetscFunctionBeginUser;
500 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
501
502 // --- 1. Get Nodal Coordinates and Output Array ---
503 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
504 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
505 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
506 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
507
508 // Define owned node ranges (global indices)
509 PetscInt xs = info.xs, xe = info.xs + info.xm;
510 PetscInt ys = info.ys, ye = info.ys + info.ym;
511 PetscInt zs = info.zs, ze = info.zs + info.zm;
512
513 // Global domain dimensions (total number of nodes)
514 PetscInt mx = info.mx;
515 PetscInt my = info.my;
516 PetscInt mz = info.mz;
517
518 // --- 2. Calculate Aj for INTERIOR Stencils ---
519
520 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
521 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
522 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
523
524 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
525 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
526 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
527
528 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
529 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
530 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
531
532 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) );
533
534 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) );
535
536 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) );
537
538 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) );
539
540 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) );
541
542 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) );
543
544 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) );
545
546 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) );
547
548 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) );
549
550 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);
551 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
552 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
553 }
554 }
555 }
556
557 // --- 4. Boundary Extrapolation for Aj ---
558 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
559 PetscInt i_bnd, j_bnd, k_bnd;
560
561 if (xs == 0) {
562 i_bnd = 0;
563 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
564 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
565 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
566 }
567 }
568 }
569 if (xe == mx) {
570 i_bnd = mx - 1;
571 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
572 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
573 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
574 }
575 }
576 }
577 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
578 if (ys == 0) {
579 j_bnd = 0;
580 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
581 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
582 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
583 }
584 }
585 }
586 if (ye == my) {
587 j_bnd = my - 1;
588 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
589 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
590 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
591 }
592 }
593 }
594 if (zs == 0) {
595 k_bnd = 0;
596 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
597 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
598 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
599 }
600 }
601 }
602 if (ze == mz) {
603 k_bnd = mz - 1;
604 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
605 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
606 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
607 }
608 }
609 }
610
611 // --- 5. Restore arrays ---
612 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
613 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
614
615 // --- 6. Assemble Global Vector ---
616 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
617 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
618 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
619
620 // --- 7. Update Local Ghosted Version ---
621 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
622 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
623
624 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
625 PetscFunctionReturn(0);
626}
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
<your name>

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

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

Definition at line 195 of file Metric.c.

196{
197 PetscErrorCode ierr;
198 PetscReal aj_min, aj_max;
199 PetscMPIInt rank;
200
201 PetscFunctionBeginUser;
202
203 /* ---------------- step 1: global extrema of Aj ---------------- */
204 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
205
206 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr); /* already global */
207 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
208
209 LOG_ALLOW(GLOBAL, LOG_INFO,
210 "[orientation] Global Aj range: [%.3e , %.3e]\n",
211 (double)aj_min, (double)aj_max);
212
213 /* ---------------- step 2: detect malformed mesh ---------------- */
214 if (aj_min < 0.0 && aj_max > 0.0)
215 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
216 "Mixed Jacobian signs detected – grid is topologically inconsistent.");
217
218 /* Default: grid is right-handed unless proven otherwise */
219 PetscInt orientation = +1;
220
221 /* ---------------- step 3: repair left-handed mesh -------------- */
222 if (aj_max < 0.0) { /* entire domain has Aj < 0 */
223 orientation = -1;
224
225 if (!rank)
226 LOG_ALLOW(LOCAL, LOG_INFO,
227 "[orientation] Detected left-handed grid – flipping metric vectors\n");
228
229 /* Flip sign of *all* metric vectors and Aj */
230 ierr = VecScale(user->Csi, -1.0); CHKERRQ(ierr);
231 ierr = VecScale(user->Eta, -1.0); CHKERRQ(ierr);
232 ierr = VecScale(user->Zet, -1.0); CHKERRQ(ierr);
233 ierr = VecScale(user->Aj , -1.0); CHKERRQ(ierr);
234
235 /* Local ghost regions now stale – refresh */
236 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
237 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
238 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
239 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
240
241 /* Sanity print: Aj must be > 0 now */
242 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
243 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
244
245 if (aj_min <= 0.0)
246 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
247 "Failed to flip grid orientation – Aj still non-positive.");
248 else if (aj_min && aj_max > 0.0)
249 orientation = +1;
250 }
251
252 /* ---------------- step 4: store result in UserCtx -------------- */
253 user->GridOrientation = orientation;
254
255 if (!rank)
256 LOG_ALLOW(LOCAL, LOG_INFO,
257 "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
258 (orientation>0) ? "right" : "left", orientation);
259
260 PetscFunctionReturn(0);
261}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:643
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeCellCentersAndSpacing()

PetscErrorCode ComputeCellCentersAndSpacing ( UserCtx *  user)

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

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

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

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

Parameters
userThe UserCtx for a specific grid level. The function populates user->Cent and user->GridSpace.
Returns
PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 644 of file Metric.c.

645{
646 PetscErrorCode ierr;
647 DMDALocalInfo info;
648 Vec lCoords;
649 const Cmpnts ***coor;
650 Cmpnts ***cent, ***gs;
651 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
652
653 PetscFunctionBeginUser;
654 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);
655
656 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
657 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
658 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
659
660 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
661 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
662
663 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
664 for (PetscInt k=info.zs+1; k<info.zs+info.zm; k++) {
665 for (PetscInt j=info.ys+1; j<info.ys+info.ym; j++) {
666 for (PetscInt i=info.xs+1; i<info.xs+info.xm; i++) {
667 // Calculate cell center as the average of its 8 corner nodes
668 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);
669 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);
670 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);
671
672 // Calculate Grid Spacing in i-direction (distance between i-face centers)
673 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);
674 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);
675 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);
676 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);
677 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);
678 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);
679 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
680
681 // Calculate Grid Spacing in j-direction (distance between j-face centers)
682 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);
683 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);
684 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);
685 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);
686 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);
687 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);
688 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
689
690 // Calculate Grid Spacing in k-direction (distance between k-face centers)
691 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);
692 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);
693 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);
694 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);
695 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);
696 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);
697 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
698 }
699 }
700 }
701
702 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
703 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
704 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
705
706 // Assemble and update ghost regions for the new data
707 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
708 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
709 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
710 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
711
712 PetscFunctionReturn(0);
713}
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:673
SimCtx::rank
PetscMPIInt rank
Definition variables.h:516
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:633
UserCtx::_this
PetscInt _this
Definition variables.h:643
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:689
UserCtx::Cent
Vec Cent
Definition variables.h:673
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.

This function calculates the metric terms (ICsi, IEta, IZet) and the inverse Jacobian (IAj) located at the geometric center of each constant-i face. This is a critical step for staggered-grid finite difference schemes.

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

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

Definition at line 739 of file Metric.c.

740{
741 PetscErrorCode ierr;
742 DMDALocalInfo info;
743 Vec lCoords;
744 const Cmpnts ***coor;
745 Cmpnts ***centx; //***gs;
746 const Cmpnts ***centx_const;
747 Cmpnts ***icsi, ***ieta, ***izet;
748 PetscScalar ***iaj;
749 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
750
751 PetscFunctionBeginUser;
752 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
753
754 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
755 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
756 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
757 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
758 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
759 PetscInt gys = info.gys, gye = info.gys + info.gym;
760 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
761
762 PetscInt lxs = xs; PetscInt lxe = xe;
763 PetscInt lys = ys; PetscInt lye = ye;
764 PetscInt lzs = zs; PetscInt lze = ze;
765
766 if (xs==0) lxs = xs+1;
767 if (ys==0) lys = ys+1;
768 if (zs==0) lzs = zs+1;
769
770 if (xe==mx) lxe=xe-1;
771 if (ye==my) lye=ye-1;
772 if (ze==mz) lze=ze-1;
773
774 // --- Part 1: Calculate the location of i-face centers (Centx) ---
775 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
776 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
777 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
778 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
779
780 // Loop over the ghosted region to calculate all local face centers
781 for (PetscInt k = gzs + 1; k < gze; k++) {
782 for (PetscInt j = gys + 1; j < gye; j++) {
783 for (PetscInt i = gxs; i < gxe; i++) {
784 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);
785 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);
786 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);
787 }
788 }
789 }
790
791 /*
792 if(xs==0){
793 for(PetscInt k=gzs+1;k < gze; k++){
794 for(PetscInt j=gys+1;j < gye; j++){
795 PetscInt i=0;
796 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
797 centx[k][j][i-1].y=centx[k][j][i].y;
798 centx[k][j][i-1].z=centx[k][j][i].z;
799 }
800 }
801 }
802 if (xe==mx){
803 for(PetscInt k=gzs+1; k<gze; k++) {
804 for (PetscInt j=gys+1; j<gye;j++) {
805 PetscInt i=mx-1;
806 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
807 centx[k][j][i].y=centx[k][j][i-1].y;
808 centx[k][j][i].z=centx[k][j][i-1].z;
809 }
810 }
811 }
812 */
813
814 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
815 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
816 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
817
818 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
819
820 // --- Part 2: Calculate metrics using face-centered coordinates ---
821 ierr = DMDAVecGetArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
822 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
823 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
824 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
825 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
826
827 // Loop over the OWNED region where we will store the final metrics
828 for (PetscInt k=lzs; k<lze; k++) {
829 for (PetscInt j=lys; j<lye; j++) {
830 for (PetscInt i=xs; i<lxe; i++) {
831
832 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
833 if (i == 0) { // Forward difference at the domain's min-i boundary
834 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
835 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
836 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
837 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
838 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
839 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
840 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
841 } else { // Central difference in the interior
842 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
843 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
844 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
845 }
846
847 // --- Stencil Logic for d/deta (derivative in j-direction) ---
848 if (j == 1) { // Forward difference
849 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
850 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
851 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
852 } else if (j == my - 2) { // Backward difference
853 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
854 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
855 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
856 } else { // Central difference
857 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
858 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
859 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
860 }
861
862 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
863 if (k == 1) { // Forward difference
864 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
865 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
866 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
867 } else if (k == mz - 2) { // Backward difference
868 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
869 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
870 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
871 } else { // Central difference
872 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
873 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
874 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
875 }
876
877 // --- Metric calculations (identical to legacy FormMetrics) ---
878 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
879 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
880 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
881
882 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
883 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
884 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
885
886 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
887 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
888 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
889
890 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
891 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
892 iaj[k][j][i] = 1.0 / iaj[k][j][i];
893 }
894 }
895 }
896 }
897
898 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
899 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
900 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
901 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
902 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
903
904 // --- Part 3: Assemble global vectors and update local ghosts ---
905 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
906 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
907 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
908 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
909
910 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
911 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
912 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
913 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
914
915 PetscFunctionReturn(0);
916}
UserCtx::IZet
Vec IZet
Definition variables.h:675
UserCtx::IEta
Vec IEta
Definition variables.h:675
UserCtx::Centx
Vec Centx
Definition variables.h:674
UserCtx::ICsi
Vec ICsi
Definition variables.h:675
UserCtx::IAj
Vec IAj
Definition variables.h:675
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeJFaceMetrics()

PetscErrorCode ComputeJFaceMetrics ( UserCtx *  user)

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

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

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

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

Definition at line 942 of file Metric.c.

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

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

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

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

Definition at line 1145 of file Metric.c.

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

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

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

Definition at line 1416 of file Metric.c.

1417{
1418 DM da = user->da, fda = user->fda;
1419 DMDALocalInfo info = user->info;
1420 PetscInt xs = info.xs, xe = info.xs + info.xm;
1421 PetscInt ys = info.ys, ye = info.ys + info.ym;
1422 PetscInt zs = info.zs, ze = info.zs + info.zm;
1423 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1424 PetscInt lxs, lys, lzs, lxe, lye, lze;
1425 PetscInt i, j, k;
1426 Vec Div;
1427 PetscReal ***div, ***aj;
1428 Cmpnts ***csi, ***eta, ***zet;
1429 PetscReal maxdiv;
1430
1431 PetscFunctionBeginUser;
1432
1433 lxs = xs; lxe = xe;
1434 lys = ys; lye = ye;
1435 lzs = zs; lze = ze;
1436
1437 if (xs == 0) lxs = xs + 1;
1438 if (ys == 0) lys = ys + 1;
1439 if (zs == 0) lzs = zs + 1;
1440
1441 if (xe == mx) lxe = xe - 1;
1442 if (ye == my) lye = ye - 1;
1443 if (ze == mz) lze = ze - 1;
1444
1445 DMDAVecGetArray(fda, user->lCsi, &csi);
1446 DMDAVecGetArray(fda, user->lEta, &eta);
1447 DMDAVecGetArray(fda, user->lZet, &zet);
1448 DMDAVecGetArray(da, user->lAj, &aj);
1449
1450 VecDuplicate(user->P, &Div);
1451 VecSet(Div, 0.);
1452 DMDAVecGetArray(da, Div, &div);
1453
1454 for (k = lzs; k < lze; k++) {
1455 for (j = lys; j < lye; j++) {
1456 for (i = lxs; i < lxe; i++) {
1457 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1458 eta[k][j][i].x - eta[k][j-1][i].x +
1459 zet[k][j][i].x - zet[k-1][j][i].x +
1460 csi[k][j][i].y - csi[k][j][i-1].y +
1461 eta[k][j][i].y - eta[k][j-1][i].y +
1462 zet[k][j][i].y - zet[k-1][j][i].y +
1463 csi[k][j][i].z - csi[k][j][i-1].z +
1464 eta[k][j][i].z - eta[k][j-1][i].z +
1465 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1466 div[k][j][i] = fabs(divergence);
1467 }
1468 }
1469 }
1470
1471 DMDAVecRestoreArray(da, Div, &div);
1472
1473 PetscReal MaxFlatIndex;
1474 VecMax(Div, &MaxFlatIndex, &maxdiv);
1475 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %d.\n",maxdiv,MaxFlatIndex);
1476
1477 for (k=zs; k<ze; k++) {
1478 for (j=ys; j<ye; j++) {
1479 for (i=xs; i<xe; i++) {
1480 if (Gidx(i,j,k,user) == MaxFlatIndex) {
1481 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,k,j,i);
1482 }
1483 }
1484 }
1485 }
1486
1487
1488 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1489 DMDAVecRestoreArray(fda, user->lEta, &eta);
1490 DMDAVecRestoreArray(fda, user->lZet, &zet);
1491 DMDAVecRestoreArray(da, user->lAj, &aj);
1492 VecDestroy(&Div);
1493
1494 PetscFunctionReturn(0);
1495}
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition Metric.c:1387
UserCtx::lZet
Vec lZet
Definition variables.h:673
UserCtx::lCsi
Vec lCsi
Definition variables.h:673
UserCtx::lAj
Vec lAj
Definition variables.h:673
UserCtx::info
DMDALocalInfo info
Definition variables.h:637
UserCtx::lEta
Vec lEta
Definition variables.h:673
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeMetricNorms()

PetscErrorCode ComputeMetricNorms ( UserCtx *  user)

Computes the max-min values of the grid metrics.

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

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

Definition at line 1506 of file Metric.c.

1507{
1508
1509 DMDALocalInfo info = user->info;
1510 PetscInt xs = info.xs, xe = info.xs + info.xm;
1511 PetscInt ys = info.ys, ye = info.ys + info.ym;
1512 PetscInt zs = info.zs, ze = info.zs + info.zm;
1513 PetscInt i, j, k;
1514
1515 PetscFunctionBeginUser;
1516
1517 PetscReal CsiMax, EtaMax, ZetMax;
1518 PetscReal ICsiMax, IEtaMax, IZetMax;
1519 PetscReal JCsiMax, JEtaMax, JZetMax;
1520 PetscReal KCsiMax, KEtaMax, KZetMax;
1521 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
1522
1523 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
1524 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
1525 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
1526 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
1527 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
1528
1529 // Max Values
1530 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
1531 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
1532 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
1533
1534 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
1535 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
1536 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
1537
1538 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
1539 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
1540 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
1541
1542 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
1543 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
1544 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
1545
1546 VecMax(user->lAj,&AjMaxArg,&AjMax);
1547 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1548 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1549 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1550
1551 VecMax(user->lAj,&AjMaxArg,&AjMax);
1552 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1553 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1554 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1555
1556 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
1557
1558 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
1559 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
1560 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
1561 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
1562 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
1563
1564 for (k=zs; k<ze; k++) {
1565 for (j=ys; j<ye; j++) {
1566 for (i=xs; i<xe; i++) {
1567 if (Gidx(i,j,k,user) == CsiMaxArg) {
1568 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
1569 }
1570 if (Gidx(i,j,k,user) == EtaMaxArg) {
1571 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
1572 }
1573 if (Gidx(i,j,k,user) == ZetMaxArg) {
1574 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
1575 }
1576 if (Gidx(i,j,k,user) == ICsiMaxArg) {
1577 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
1578 }
1579 if (Gidx(i,j,k,user) == IEtaMaxArg) {
1580 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
1581 }
1582 if (Gidx(i,j,k,user) == IZetMaxArg) {
1583 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
1584 }
1585 if (Gidx(i,j,k,user) == JCsiMaxArg) {
1586 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
1587 }
1588 if (Gidx(i,j,k,user) == JEtaMaxArg) {
1589 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
1590 }
1591 if (Gidx(i,j,k,user) == JZetMaxArg) {
1592 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
1593 }
1594 if (Gidx(i,j,k,user) == KCsiMaxArg) {
1595 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
1596 }
1597 if (Gidx(i,j,k,user) == KEtaMaxArg) {
1598 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
1599 }
1600 if (Gidx(i,j,k,user) == KZetMaxArg) {
1601 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
1602 }
1603 if (Gidx(i,j,k,user) == AjMaxArg) {
1604 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
1605 }
1606 if (Gidx(i,j,k,user) == IAjMaxArg) {
1607 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
1608 }
1609 if (Gidx(i,j,k,user) == JAjMaxArg) {
1610 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
1611 }
1612 if (Gidx(i,j,k,user) == KAjMaxArg) {
1613 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
1614 }
1615 }
1616 }
1617 }
1618
1619 /*
1620 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
1621 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
1622 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
1623 */
1624
1625 PetscFunctionReturn(0);
1626}
UserCtx::lIEta
Vec lIEta
Definition variables.h:675
UserCtx::lIZet
Vec lIZet
Definition variables.h:675
UserCtx::lIAj
Vec lIAj
Definition variables.h:675
UserCtx::lKEta
Vec lKEta
Definition variables.h:677
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:676
UserCtx::lKZet
Vec lKZet
Definition variables.h:677
UserCtx::lJEta
Vec lJEta
Definition variables.h:676
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:677
UserCtx::lJZet
Vec lJZet
Definition variables.h:676
UserCtx::lICsi
Vec lICsi
Definition variables.h:675
UserCtx::lJAj
Vec lJAj
Definition variables.h:676
UserCtx::lKAj
Vec lKAj
Definition variables.h:677
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculateAllGridMetrics()

PetscErrorCode CalculateAllGridMetrics ( SimCtx *  simCtx)

Orchestrates the calculation of all grid metrics.

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

Parameters
simCtxThe master SimCtx, containing the configured UserCtx hierarchy.
Returns
PetscErrorCode

Definition at line 1639 of file Metric.c.

1640{
1641 PetscErrorCode ierr;
1642 UserMG *usermg = &simCtx->usermg;
1643 MGCtx *mgctx = usermg->mgctx;
1644 PetscInt nblk = simCtx->block_number;
1645
1646 PetscFunctionBeginUser;
1647 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
1648
1649 // Loop through all levels and all blocks
1650 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
1651 for (PetscInt bi = 0; bi < nblk; bi++) {
1652 UserCtx *user = &mgctx[level].user[bi];
1653 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
1654
1655 // Call the modern, modular helper functions for each UserCtx.
1656 // These functions are self-contained and operate on the data within the provided context.
1657 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
1658 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
1659 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
1660 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
1661 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
1662 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
1663 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
1664
1665 // Diagnostics
1666 ierr = ComputeMetricNorms(user);
1667 if (level == usermg->mglevels - 1) {
1668 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
1669 }
1670 }
1671 }
1672
1673 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
1674 PetscFunctionReturn(0);
1675}
ComputeMetricNorms
PetscErrorCode ComputeMetricNorms(UserCtx *user)
Computes the max-min values of the grid metrics.
Definition Metric.c:1506
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain...
Definition Metric.c:195
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Computes the physical location of cell centers and the spacing between them.
Definition Metric.c:644
ComputeJFaceMetrics
PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
Computes metrics centered on constant-j faces (j-faces).
Definition Metric.c:942
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Computes the divergence of the grid metrics and identifies the maximum value.
Definition Metric.c:1416
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation,...
Definition Metric.c:286
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Calculates the cell-centered inverse Jacobian determinant (1/J), including boundary extrapolation,...
Definition Metric.c:491
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Computes metrics centered on constant-k faces (k-faces).
Definition Metric.c:1145
ComputeIFaceMetrics
PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
Computes metrics centered on constant-i faces (i-faces).
Definition Metric.c:739
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:274
MGCtx::user
UserCtx * user
Definition variables.h:418
SimCtx::block_number
PetscInt block_number
Definition variables.h:562
SimCtx::usermg
UserMG usermg
Definition variables.h:599
UserMG::mglevels
PetscInt mglevels
Definition variables.h:425
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:428
MGCtx
Context for Multigrid operations.
Definition variables.h:417
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:630
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:424
Here is the call graph for this function:
Here is the caller graph for this function: