PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Metric.c
Go to the documentation of this file.
1/* Metric.c ------------------------------------------------------------------
2 * Utility routines for curvilinear-grid metric operations used by the
3 * particle-swarm module.
4 *
5 * – Logical (xi,eta,zta) → physical (x,y,z) mapping via trilinear blend.
6 * – Jacobian matrix and determinant for contravariant velocity conversion.
7 *
8 * The only data this file needs from the application is the DMDA that stores
9 * vertex coordinates in the usual PETSc coordinate DM (user->da) and the
10 * coordinate array type `Cmpnts` (three-component struct {x,y,z}).
11 * ---------------------------------------------------------------------------*/
12
13#include <petsc.h>
14#include "Metric.h" /* forward declarations + Cmpnts + UserCtx */
15
16
17#undef __FUNCT__
18#define __FUNCT__ "MetricGetCellVertices"
19/* ------------------------------------------------------------------------- */
20/**
21 * @brief Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
22 *
23 * Vertices are returned in the standard trilinear ordering: bit 0 → x-corner,
24 * bit 1 → y-corner, bit 2 → z-corner. (000 = origin of the cell, 111 = far
25 * corner.)
26 */
27PetscErrorCode MetricGetCellVertices(UserCtx *user,
28 const Cmpnts ***X, /* coord array */
29 PetscInt i,PetscInt j,PetscInt k,
30 Cmpnts V[8])
31{
32 PetscFunctionBeginUser;
33 for (PetscInt c = 0; c < 8; ++c) {
34 PetscInt ii = i + ((c & 1) ? 1 : 0);
35 PetscInt jj = j + ((c & 2) ? 1 : 0);
36 PetscInt kk = k + ((c & 4) ? 1 : 0);
37 LOG_LOOP_ALLOW(GLOBAL, LOG_VERBOSE,i+j+k,10," ii: %d,jj:%d,kk:%d - Retrieved.\n",ii,jj,kk);
38 V[c] = X[kk][jj][ii];
39 }
40 PetscFunctionReturn(0);
41}
42
43#undef __FUNCT__
44#define __FUNCT__ "TrilinearBlend"
45
46/* ------------------------------------------------------------------------- */
47/**
48 * @brief Map logical coordinates to physical space using trilinear basis.
49 *
50 * @param[in] V Array of the eight vertex coordinates (MetricGetCellVertices).
51 * @param[in] xi,eta,zta Logical coordinates in [0,1].
52 * @param[out] Xp Physical coordinate.
53 */
54static inline void TrilinearBlend(const Cmpnts V[8],
55 PetscReal xi,PetscReal eta,PetscReal zta,
56 Cmpnts *Xp)
57{
58 PetscReal x=0,y=0,z=0;
59 for (PetscInt c=0;c<8;++c) {
60 PetscReal N = ((c&1)?xi : 1.0-xi ) *
61 ((c&2)?eta: 1.0-eta) *
62 ((c&4)?zta: 1.0-zta);
63 x += N * V[c].x;
64 y += N * V[c].y;
65 z += N * V[c].z;
66 }
67 Xp->x = x; Xp->y = y; Xp->z = z;
68}
69
70
71#undef _FUNCT__
72#define __FUNCT__ "MetricLogicalToPhysical"
73/* ------------------------------------------------------------------------- */
74/**
75 * @brief Public wrapper: map (cell index, ξ,η,ζ) to (x,y,z).
76 */
77PetscErrorCode MetricLogicalToPhysical(UserCtx *user,
78 const Cmpnts ***X,
79 PetscInt i,PetscInt j,PetscInt k,
80 PetscReal xi,PetscReal eta,PetscReal zta,
81 Cmpnts *Xp)
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}
96
97#undef __FUNCT__
98#define __FUNCT__ "MetricJacobian"
99/* ------------------------------------------------------------------------- */
100/**
101 * @brief Compute Jacobian matrix and its determinant at (xi,eta,zta).
102 *
103 * J = [ x_ξ x_η x_ζ ]
104 * [ y_ξ y_η y_ζ ]
105 * [ z_ξ z_η z_ζ ]
106 *
107 * This is handy for converting physical velocities (u,v,w) into contravariant
108 * components and for volume weighting.
109 */
110PetscErrorCode MetricJacobian(UserCtx *user,
111 const Cmpnts ***X,
112 PetscInt i,PetscInt j,PetscInt k,
113 PetscReal xi,PetscReal eta,PetscReal zta,
114 PetscReal J[3][3], PetscReal *detJ)
115{
116 PetscErrorCode ierr;
117 Cmpnts V[8];
118 PetscFunctionBeginUser;
119
120 PROFILE_FUNCTION_BEGIN;
121
122 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
123
124 /* derivatives of trilinear shape functions */
125 PetscReal dN_dXi[8], dN_dEta[8], dN_dZta[8];
126 for (PetscInt c=0;c<8;++c) {
127 PetscReal sx = (c & 1) ? 1.0 : -1.0;
128 PetscReal sy = (c & 2) ? 1.0 : -1.0;
129 PetscReal sz = (c & 4) ? 1.0 : -1.0;
130 dN_dXi [c] = 0.125 * sx * ( (c&2?eta:1-eta) ) * ( (c&4?zta:1-zta) );
131 dN_dEta[c] = 0.125 * sy * ( (c&1?xi :1-xi ) ) * ( (c&4?zta:1-zta) );
132 dN_dZta[c] = 0.125 * sz * ( (c&1?xi :1-xi ) ) * ( (c&2?eta:1-eta) );
133 }
134
135 /* assemble Jacobian */
136 PetscReal x_xi=0,y_xi=0,z_xi=0,
137 x_eta=0,y_eta=0,z_eta=0,
138 x_zta=0,y_zta=0,z_zta=0;
139 for (PetscInt c=0;c<8;++c) {
140 x_xi += dN_dXi [c]*V[c].x; y_xi += dN_dXi [c]*V[c].y; z_xi += dN_dXi [c]*V[c].z;
141 x_eta += dN_dEta[c]*V[c].x; y_eta += dN_dEta[c]*V[c].y; z_eta += dN_dEta[c]*V[c].z;
142 x_zta += dN_dZta[c]*V[c].x; y_zta += dN_dZta[c]*V[c].y; z_zta += dN_dZta[c]*V[c].z;
143 }
144
145 J[0][0]=x_xi; J[0][1]=x_eta; J[0][2]=x_zta;
146 J[1][0]=y_xi; J[1][1]=y_eta; J[1][2]=y_zta;
147 J[2][0]=z_xi; J[2][1]=z_eta; J[2][2]=z_zta;
148
149 if (detJ) {
150 *detJ = x_xi*(y_eta*z_zta - y_zta*z_eta)
151 - x_eta*(y_xi*z_zta - y_zta*z_xi)
152 + x_zta*(y_xi*z_eta - y_eta*z_xi);
153 }
154
155 PROFILE_FUNCTION_END;
156
157 PetscFunctionReturn(0);
158}
159
160
161#undef __FUNCT__
162#define __FUNCT__ "MetricVelocityContravariant"
163/* ------------------------------------------------------------------------- */
164/**
165 * @brief Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).
166 */
167PetscErrorCode MetricVelocityContravariant(const PetscReal J[3][3], PetscReal detJ,
168 const PetscReal u[3], PetscReal uc[3])
169{
170 PetscFunctionBeginUser;
171
172 PROFILE_FUNCTION_BEGIN;
173
174 /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
175 PetscReal gxi[3] = { J[1][1]*J[2][2]-J[1][2]*J[2][1],
176 -J[0][1]*J[2][2]+J[0][2]*J[2][1],
177 J[0][1]*J[1][2]-J[0][2]*J[1][1] };
178 PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
179 J[0][0]*J[2][2]-J[0][2]*J[2][0],
180 -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
181 PetscReal gzta[3] = { J[1][0]*J[2][1]-J[1][1]*J[2][0],
182 -J[0][0]*J[2][1]+J[0][1]*J[2][0],
183 J[0][0]*J[1][1]-J[0][1]*J[1][0] };
184
185 PetscReal invDet = 1.0 / detJ;
186 for (int d=0; d<3; ++d) { gxi[d] *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
187
188 uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
189 uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
190 uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
191
192 PROFILE_FUNCTION_END;
193
194 PetscFunctionReturn(0);
195}
196
197/////////////////////////////////////////////////////////////////////////////////
198#undef __FUNCT__
199#define __FUNCT__ "InvertCovariantMetricTensor"
200/**
201 * @brief Inverts the 3x3 covariant metric tensor to obtain the contravariant metric tensor.
202 *
203 * In curvilinear coordinates, the input matrix `g` contains the dot products of the
204 * covariant basis vectors (e.g., g_ij = e_i . e_j). Its inverse, `G`, is the
205 * contravariant metric tensor, which is essential for transforming vectors and tensors
206 * between coordinate systems.
207 *
208 * @param covariantTensor Input: A 3x3 matrix representing the covariant metric tensor.
209 * @param contravariantTensor Output: A 3x3 matrix where the inverted result is stored.
210 * @return PetscErrorCode 0 on success.
211 */
212PetscErrorCode InvertCovariantMetricTensor(double covariantTensor[3][3], double contravariantTensor[3][3])
213{
214 PetscFunctionBeginUser;
215
216 const double a11=covariantTensor[0][0], a12=covariantTensor[0][1], a13=covariantTensor[0][2];
217 const double a21=covariantTensor[1][0], a22=covariantTensor[1][1], a23=covariantTensor[1][2];
218 const double a31=covariantTensor[2][0], a32=covariantTensor[2][1], a33=covariantTensor[2][2];
219
220 double det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
221
222 if (fabs(det) < 1.0e-12) {
223 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Matrix is singular, determinant is near zero.");
224 }
225
226 contravariantTensor[0][0] = (a33*a22-a32*a23)/det;
227 contravariantTensor[0][1] = -(a33*a12-a32*a13)/det;
228 contravariantTensor[0][2] = (a23*a12-a22*a13)/det;
229 contravariantTensor[1][0] = -(a33*a21-a31*a23)/det;
230 contravariantTensor[1][1] = (a33*a11-a31*a13)/det;
231 contravariantTensor[1][2] = -(a23*a11-a21*a13)/det;
232 contravariantTensor[2][0] = (a32*a21-a31*a22)/det;
233 contravariantTensor[2][1] = -(a32*a11-a31*a12)/det;
234 contravariantTensor[2][2] = (a22*a11-a21*a12)/det;
235
236 PetscFunctionReturn(0);
237}
238
239
240#undef __FUNCT__
241#define __FUNCT__ "CalculateFaceNormalAndArea"
242/**
243 * @brief Computes the unit normal vectors and areas of the three faces of a computational cell.
244 *
245 * Given the metric vectors (csi, eta, zet), this function calculates the geometric
246 * properties of the cell faces aligned with the i, j, and k directions.
247 *
248 * @param csi, eta, zet The metric vectors at the cell center.
249 * @param ni Output: A 3-element array for the unit normal vector of the i-face.
250 * @param nj Output: A 3-element array for the unit normal vector of the j-face.
251 * @param nk Output: A 3-element array for the unit normal vector of the k-face.
252 * @param Ai Output: Pointer to store the area of the i-face.
253 * @param Aj Output: Pointer to store the area of the j-face.
254 * @param Ak Output: Pointer to store the area of the k-face.
255 * @return PetscErrorCode 0 on success.
256 */
257PetscErrorCode CalculateFaceNormalAndArea(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
258{
259 PetscFunctionBeginUser;
260 PROFILE_FUNCTION_BEGIN;
261 double g[3][3];
262 double G[3][3];
263
264 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
265 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
266 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
267
268 InvertCovariantMetricTensor(g, G);
269 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
270 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
271 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
272
273 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
274 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
275 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
276
277 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
278 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
279 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
280
281 *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
282 *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
283 *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
284
285 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
286 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
287 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
288
289 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
290 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
291 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
292
293 PROFILE_FUNCTION_END;
294 PetscFunctionReturn(0);
295}
296
297#undef __FUNCT__
298#define __FUNCT__ "ComputeCellCharacteristicLengthScale"
299/**
300 * @brief Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.
301 *
302 * For a non-uniform, non-orthogonal cell, there is no single "dx". This function
303 * computes an effective length scale in each Cartesian direction based on the cell
304 * volume and the areas of its faces.
305 *
306 * @param ajc The Jacobian of the grid transformation (1 / cell volume).
307 * @param csi, eta, zet The metric vectors at the cell center.
308 * @param dx Output: Pointer to store the characteristic length in the x-direction.
309 * @param dy Output: Pointer to store the characteristic length in the y-direction.
310 * @param dz Output: Pointer to store the characteristic length in the z-direction.
311 * @return PetscErrorCode 0 on success.
312 */
313PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
314{
315 PetscFunctionBeginUser;
316 PROFILE_FUNCTION_BEGIN;
317 double ni[3], nj[3], nk[3];
318 double Li, Lj, Lk;
319 double Ai, Aj, Ak;
320 double vol = 1./ajc;
321
322 CalculateFaceNormalAndArea(csi, eta, zet, ni, nj, nk, &Ai, &Aj, &Ak);
323 Li = vol / Ai;
324 Lj = vol / Aj;
325 Lk = vol / Ak;
326
327 // Length scale vector = di * ni_vector + dj * nj_vector + dk * nk_vector
328 *dx = fabs( Li * ni[0] + Lj * nj[0] + Lk * nk[0] );
329 *dy = fabs( Li * ni[1] + Lj * nj[1] + Lk * nk[1] );
330 *dz = fabs( Li * ni[2] + Lj * nj[2] + Lk * nk[2] );
331
332 PROFILE_FUNCTION_END;
333 PetscFunctionReturn(0);
334}
335
336
337#undef __FUNCT__
338#define __FUNCT__ "CheckAndFixGridOrientation"
339/* -------------------------------------------------------------------------- */
340/**
341 * @brief Ensure a **right-handed** metric basis (`Csi`, `Eta`, `Zet`) and a
342 * **positive Jacobian** (`Aj`) over the whole domain.
343 *
344 * The metric-generation kernels are completely algebraic, so they will happily
345 * deliver a *left-handed* basis if the mesh file enumerates nodes in the
346 * opposite ζ-direction.
347 * This routine makes the orientation explicit and—if needed—repairs it
348 * **once per run**:
349 *
350 * | Step | Action |
351 * |------|--------|
352 * | 1 | Compute global `Aj_min`, `Aj_max`. |
353 * | 2 | **Mixed signs** (`Aj_min < 0 && Aj_max > 0`) &rarr; abort: the mesh is topologically inconsistent. |
354 * | 3 | **All negative** (`Aj_max < 0`) &rarr; flip <br>`Csi`, `Eta`, `Zet`, `Aj` & update local ghosts. |
355 * | 4 | Store `user->orientation = ±1` so BC / IC routines can apply sign-aware logic if they care about inlet direction. |
356 *
357 * @param[in,out] user Fully initialised #UserCtx that already contains
358 * `Csi`, `Eta`, `Zet`, `Aj`, their **local** ghosts, and
359 * valid distributed DMs.
360 *
361 * @return `0` on success or a PETSc error code on failure.
362 *
363 * @note Call **immediately after** `ComputeCellCenteredJacobianInverse()` and
364 * before any routine that differentiates or applies BCs.
365 *
366 * @author vishal kandala
367 */
368PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
369{
370 PetscErrorCode ierr;
371 PetscReal aj_min, aj_max;
372 PetscMPIInt rank;
373
374 PetscFunctionBeginUser;
375
376 PROFILE_FUNCTION_BEGIN;
377
378 /* ---------------- step 1: global extrema of Aj ---------------- */
379 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
380
381 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr); /* already global */
382 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
383
384 LOG_ALLOW(GLOBAL, LOG_INFO,
385 "[orientation] Global Aj range: [%.3e , %.3e]\n",
386 (double)aj_min, (double)aj_max);
387
388 /* ---------------- step 2: detect malformed mesh ---------------- */
389 if (aj_min < 0.0 && aj_max > 0.0)
390 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
391 "Mixed Jacobian signs detected – grid is topologically inconsistent.");
392
393 /* Default: grid is right-handed unless proven otherwise */
394 PetscInt orientation = +1;
395
396 /* ---------------- step 3: repair left-handed mesh -------------- */
397 if (aj_max < 0.0) { /* entire domain has Aj < 0 */
398 orientation = -1;
399
400 if (!rank)
401 LOG_ALLOW(LOCAL, LOG_INFO,
402 "[orientation] Detected left-handed grid – flipping metric vectors\n");
403
404 /* Flip sign of *all* metric vectors and Aj */
405 ierr = VecScale(user->Csi, -1.0); CHKERRQ(ierr);
406 ierr = VecScale(user->Eta, -1.0); CHKERRQ(ierr);
407 ierr = VecScale(user->Zet, -1.0); CHKERRQ(ierr);
408 ierr = VecScale(user->Aj , -1.0); CHKERRQ(ierr);
409
410 /* Local ghost regions now stale – refresh */
411 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
412 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
413 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
414 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
415
416 /* Sanity print: Aj must be > 0 now */
417 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
418 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
419
420 if (aj_min <= 0.0)
421 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
422 "Failed to flip grid orientation – Aj still non-positive.");
423 else if (aj_min && aj_max > 0.0)
424 orientation = +1;
425 }
426
427 /* ---------------- step 4: store result in UserCtx -------------- */
428 user->GridOrientation = orientation;
429
430 if (!rank)
431 LOG_ALLOW(LOCAL, LOG_INFO,
432 "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
433 (orientation>0) ? "right" : "left", orientation);
434
435 PROFILE_FUNCTION_END;
436
437 PetscFunctionReturn(0);
438}
439
440#undef __FUNCT__
441#define __FUNCT__ "ApplyPeriodicCorrectionsToCellCentersAndSpacing"
442/**
443 * @brief Applies periodic boundary corrections to cell centers (Cent) and grid spacing (GridSpace).
444 *
445 * This function handles the special logic needed when periodic boundaries are present.
446 * For coarse grids (cgrid), it directly copies from ghost regions. For fine grids,
447 * it calculates the boundary values using the GridSpace information.
448 *
449 * @param user The UserCtx containing grid and field data.
450 * @return PetscErrorCode 0 on success.
451 */
452PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing(UserCtx *user)
453{
454 PetscErrorCode ierr;
455 DMDALocalInfo info = user->info;
456 PetscInt xs = info.xs, xe = info.xs + info.xm;
457 PetscInt ys = info.ys, ye = info.ys + info.ym;
458 PetscInt zs = info.zs, ze = info.zs + info.zm;
459 PetscInt mx = info.mx, my = info.my, mz = info.mz;
460 Cmpnts ***cent, ***lcent, ***gs;
461 PetscReal delta;
462
463 PetscFunctionBeginUser;
464 PROFILE_FUNCTION_BEGIN;
465
466 // Check if any periodic boundaries exist
467 PetscBool has_periodic = PETSC_FALSE;
468 for (int i = 0; i < 6; i++) {
469 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
470 has_periodic = PETSC_TRUE;
471 break;
472 }
473 }
474
475 if (!has_periodic) {
476 LOG_ALLOW(LOCAL, LOG_TRACE, "No periodic boundaries; skipping corrections for Cent/GridSpace.\n");
477 PROFILE_FUNCTION_END;
478 PetscFunctionReturn(0);
479 }
480
481 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying periodic corrections to Cent and GridSpace.\n");
482
483 // Must update ghosts first before applying corrections
484 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
485 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
486
487 // --- X-direction periodic corrections ---
488 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
489 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
490
491 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
492 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
493 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
494
495 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
496 if (user->cgrid) {
497 for (PetscInt k=zs; k<ze; k++) {
498 for (PetscInt j=ys; j<ye; j++) {
499 cent[k][j][0] = lcent[k][j][-2];
500 }
501 }
502 } else {
503 for (PetscInt k=zs; k<ze; k++) {
504 for (PetscInt j=ys; j<ye; j++) {
505 delta = (gs[k][j][1].x + gs[k][j][-2].x) / 2.0;
506 cent[k][j][0].x = cent[k][j][1].x - delta;
507 cent[k][j][0].y = cent[k][j][1].y;
508 cent[k][j][0].z = cent[k][j][1].z;
509 }
510 }
511 }
512 }
513
514 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
515 if (user->cgrid) {
516 for (PetscInt k=zs; k<ze; k++) {
517 for (PetscInt j=ys; j<ye; j++) {
518 cent[k][j][mx-1] = lcent[k][j][mx+1];
519 }
520 }
521 } else {
522 for (PetscInt k=zs; k<ze; k++) {
523 for (PetscInt j=ys; j<ye; j++) {
524 delta = (gs[k][j][mx-2].x + gs[k][j][mx+1].x) / 2.0;
525 cent[k][j][mx-1].x = cent[k][j][mx-2].x + delta;
526 cent[k][j][mx-1].y = cent[k][j][mx-2].y;
527 cent[k][j][mx-1].z = cent[k][j][mx-2].z;
528 }
529 }
530 }
531 }
532
533 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
534 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
535 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
536 }
537 // --- Y-direction periodic corrections ---
538 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
539 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
540
541 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
542 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
543 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
544
545 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
546 if (user->cgrid) {
547 for (PetscInt k=zs; k<ze; k++) {
548 for (PetscInt i=xs; i<xe; i++) {
549 cent[k][0][i] = lcent[k][-2][i];
550 }
551 }
552 } else {
553 for (PetscInt k=zs; k<ze; k++) {
554 for (PetscInt i=xs; i<xe; i++) {
555 delta = (gs[k][1][i].y + gs[k][-2][i].y) / 2.0;
556 cent[k][0][i].x = cent[k][1][i].x;
557 cent[k][0][i].y = cent[k][1][i].y - delta;
558 cent[k][0][i].z = cent[k][1][i].z;
559 }
560 }
561 }
562 }
563
564 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
565 if (user->cgrid) {
566 for (PetscInt k=zs; k<ze; k++) {
567 for (PetscInt i=xs; i<xe; i++) {
568 cent[k][my-1][i] = lcent[k][my+1][i];
569 }
570 }
571 } else {
572 for (PetscInt k=zs; k<ze; k++) {
573 for (PetscInt i=xs; i<xe; i++) {
574 delta = (gs[k][my-2][i].y + gs[k][my+1][i].y) / 2.0;
575 cent[k][my-1][i].x = cent[k][my-2][i].x;
576 cent[k][my-1][i].y = cent[k][my-2][i].y + delta;
577 cent[k][my-1][i].z = cent[k][my-2][i].z;
578 }
579 }
580 }
581 }
582
583 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
584 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
585 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
586
587 }
588
589 // --- Z-direction periodic corrections ---
590 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
591 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
592
593 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
594 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
595 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
596
597 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
598 if (user->cgrid) {
599 for (PetscInt j=ys; j<ye; j++) {
600 for (PetscInt i=xs; i<xe; i++) {
601 cent[0][j][i] = lcent[-2][j][i];
602 }
603 }
604 } else {
605 for (PetscInt j=ys; j<ye; j++) {
606 for (PetscInt i=xs; i<xe; i++) {
607 delta = (gs[1][j][i].z + gs[-2][j][i].z) / 2.0;
608 cent[0][j][i].x = cent[1][j][i].x;
609 cent[0][j][i].y = cent[1][j][i].y;
610 cent[0][j][i].z = cent[1][j][i].z - delta;
611 }
612 }
613 }
614 }
615
616 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
617 if (user->cgrid) {
618 for (PetscInt j=ys; j<ye; j++) {
619 for (PetscInt i=xs; i<xe; i++) {
620 cent[mz-1][j][i] = lcent[mz+1][j][i];
621 }
622 }
623 } else {
624 for (PetscInt j=ys; j<ye; j++) {
625 for (PetscInt i=xs; i<xe; i++) {
626 delta = (gs[mz-2][j][i].z + gs[mz+1][j][i].z) / 2.0;
627 cent[mz-1][j][i].x = cent[mz-2][j][i].x;
628 cent[mz-1][j][i].y = cent[mz-2][j][i].y;
629 cent[mz-1][j][i].z = cent[mz-2][j][i].z + delta;
630 }
631 }
632 }
633 }
634
635 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
636 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
637 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
638
639 }
640
641 PROFILE_FUNCTION_END;
642 PetscFunctionReturn(0);
643}
644
645#undef __FUNCT__
646#define __FUNCT__ "ApplyPeriodicCorrectionsToIFaceCenter"
647/**
648 * @brief Applies periodic boundary corrections to i-face centers (Centx).
649 *
650 * Only X-direction periodicity affects Centx. This function must be called
651 * after Centx has been initially computed but before it's used for metric calculations.
652 *
653 * @param user The UserCtx containing grid and field data.
654 * @return PetscErrorCode 0 on success.
655 */
656PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter(UserCtx *user)
657{
658 PetscErrorCode ierr;
659 DMDALocalInfo info = user->info;
660 PetscInt xs = info.xs, xe = info.xs + info.xm;
661 PetscInt mx = info.mx;
662 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
663 PetscInt gys = info.gys, gye = info.gys + info.gym;
664 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
665 Cmpnts ***centx, ***gs;
666
667 PetscFunctionBeginUser;
668 PROFILE_FUNCTION_BEGIN;
669
670 // Check if X-periodic boundaries exist
671 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC &&
672 user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
673 LOG_ALLOW(LOCAL, LOG_TRACE, "No X-periodic boundaries; skipping Centx corrections.\n");
674 PROFILE_FUNCTION_END;
675 PetscFunctionReturn(0);
676 }
677
678 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying X-periodic corrections to Centx.\n");
679
680 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
681 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
682
683 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
684 if (user->cgrid) {
685 for (PetscInt k=gzs+1; k<gze; k++) {
686 for (PetscInt j=gys+1; j<gye; j++) {
687 centx[k][j][-1].x = centx[k][j][-3].x;
688 centx[k][j][-1].y = centx[k][j][-3].y;
689 centx[k][j][-1].z = centx[k][j][-3].z;
690 }
691 }
692 } else {
693 for (PetscInt k=gzs+1; k<gze; k++) {
694 for (PetscInt j=gys+1; j<gye; j++) {
695 centx[k][j][-1].x = centx[k][j][0].x - gs[k][j][-2].x;
696 centx[k][j][-1].y = centx[k][j][0].y;
697 centx[k][j][-1].z = centx[k][j][0].z;
698 }
699 }
700 }
701 }
702
703 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
704 if (user->cgrid) {
705 for (PetscInt k=gzs+1; k<gze; k++) {
706 for (PetscInt j=gys+1; j<gye; j++) {
707 centx[k][j][mx-1].x = centx[k][j][mx+1].x;
708 centx[k][j][mx-1].y = centx[k][j][mx+1].y;
709 centx[k][j][mx-1].z = centx[k][j][mx+1].z;
710 }
711 }
712 } else {
713 for (PetscInt k=gzs+1; k<gze; k++) {
714 for (PetscInt j=gys+1; j<gye; j++) {
715 centx[k][j][mx-1].x = centx[k][j][mx-2].x + gs[k][j][mx+1].x;
716 centx[k][j][mx-1].y = centx[k][j][mx-2].y;
717 centx[k][j][mx-1].z = centx[k][j][mx-2].z;
718 }
719 }
720 }
721 }
722
723 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
724 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
725
726 PROFILE_FUNCTION_END;
727 PetscFunctionReturn(0);
728}
729
730
731#undef __FUNCT__
732#define __FUNCT__ "ApplyPeriodicCorrectionsToJFaceCenter"
733/**
734 * @brief Applies periodic boundary corrections to j-face centers (Centy).
735 *
736 * Only Y-direction periodicity affects Centy. This function must be called
737 * after Centy has been initially computed but before it's used for metric calculations.
738 *
739 * @param user The UserCtx containing grid and field data.
740 * @return PetscErrorCode 0 on success.
741 */
742PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
743{
744 PetscErrorCode ierr;
745 DMDALocalInfo info = user->info;
746 PetscInt ys = info.ys, ye = info.ys + info.ym;
747 PetscInt my = info.my;
748 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
749 PetscInt gys = info.gys, gye = info.gys + info.gym;
750 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
751 Cmpnts ***centy, ***gs;
752
753 PetscFunctionBeginUser;
754 PROFILE_FUNCTION_BEGIN;
755
756 // Check if Y-periodic boundaries exist
757 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC &&
758 user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
759 LOG_ALLOW(LOCAL, LOG_TRACE, "No Y-periodic boundaries; skipping Centy corrections.\n");
760 PROFILE_FUNCTION_END;
761 PetscFunctionReturn(0);
762 }
763
764 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying Y-periodic corrections to Centy.\n");
765
766 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
767 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
768
769 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
770 if (user->cgrid) {
771 for (PetscInt k=gzs+1; k<gze; k++) {
772 for (PetscInt i=gxs+1; i<gxe; i++) {
773 centy[k][-1][i].x = centy[k][-3][i].x;
774 centy[k][-1][i].y = centy[k][-3][i].y;
775 centy[k][-1][i].z = centy[k][-3][i].z;
776 }
777 }
778 } else {
779 for (PetscInt k=gzs+1; k<gze; k++) {
780 for (PetscInt i=gxs+1; i<gxe; i++) {
781 centy[k][-1][i].x = centy[k][0][i].x;
782 centy[k][-1][i].y = centy[k][0][i].y - gs[k][-2][i].y;
783 centy[k][-1][i].z = centy[k][0][i].z;
784 }
785 }
786 }
787 }
788
789 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
790 if (user->cgrid) {
791 for (PetscInt k=gzs+1; k<gze; k++) {
792 for (PetscInt i=gxs+1; i<gxe; i++) {
793 centy[k][my-1][i].x = centy[k][my+1][i].x;
794 centy[k][my-1][i].y = centy[k][my+1][i].y;
795 centy[k][my-1][i].z = centy[k][my+1][i].z;
796 }
797 }
798 } else {
799 for (PetscInt k=gzs+1; k<gze; k++) {
800 for (PetscInt i=gxs+1; i<gxe; i++) {
801 centy[k][my-1][i].x = centy[k][my-2][i].x;
802 centy[k][my-1][i].y = centy[k][my-2][i].y + gs[k][my+1][i].y;
803 centy[k][my-1][i].z = centy[k][my-2][i].z;
804 }
805 }
806 }
807 }
808
809 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
810 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
811
812 PROFILE_FUNCTION_END;
813 PetscFunctionReturn(0);
814}
815
816
817#undef __FUNCT__
818#define __FUNCT__ "ApplyPeriodicCorrectionsToKFaceCenter"
819/**
820 * @brief Applies periodic boundary corrections to k-face centers (Centz).
821 *
822 * Only Z-direction periodicity affects Centz. This function must be called
823 * after Centz has been initially computed but before it's used for metric calculations.
824 *
825 * @param user The UserCtx containing grid and field data.
826 * @return PetscErrorCode 0 on success.
827 */
828PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
829{
830 PetscErrorCode ierr;
831 DMDALocalInfo info = user->info;
832 PetscInt zs = info.zs, ze = info.zs + info.zm;
833 PetscInt mz = info.mz;
834 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
835 PetscInt gys = info.gys, gye = info.gys + info.gym;
836 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
837 Cmpnts ***centz, ***gs;
838
839 PetscFunctionBeginUser;
840 PROFILE_FUNCTION_BEGIN;
841
842 // Check if Z-periodic boundaries exist
843 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC &&
844 user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
845 LOG_ALLOW(LOCAL, LOG_TRACE, "No Z-periodic boundaries; skipping Centz corrections.\n");
846 PROFILE_FUNCTION_END;
847 PetscFunctionReturn(0);
848 }
849
850 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying Z-periodic corrections to Centz.\n");
851
852 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
853 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
854
855 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
856 if (user->cgrid) {
857 for (PetscInt j=gys+1; j<gye; j++) {
858 for (PetscInt i=gxs+1; i<gxe; i++) {
859 centz[-1][j][i].x = centz[-3][j][i].x;
860 centz[-1][j][i].y = centz[-3][j][i].y;
861 centz[-1][j][i].z = centz[-3][j][i].z;
862 }
863 }
864 } else {
865 for (PetscInt j=gys+1; j<gye; j++) {
866 for (PetscInt i=gxs+1; i<gxe; i++) {
867 centz[-1][j][i].x = centz[0][j][i].x;
868 centz[-1][j][i].y = centz[0][j][i].y;
869 centz[-1][j][i].z = centz[0][j][i].z - gs[-2][j][i].z;
870 }
871 }
872 }
873 }
874
875 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
876 if (user->cgrid) {
877 for (PetscInt j=gys+1; j<gye; j++) {
878 for (PetscInt i=gxs+1; i<gxe; i++) {
879 centz[mz-1][j][i].x = centz[mz+1][j][i].x;
880 centz[mz-1][j][i].y = centz[mz+1][j][i].y;
881 centz[mz-1][j][i].z = centz[mz+1][j][i].z;
882 }
883 }
884 } else {
885 for (PetscInt j=gys+1; j<gye; j++) {
886 for (PetscInt i=gxs+1; i<gxe; i++) {
887 centz[mz-1][j][i].x = centz[mz-2][j][i].x;
888 centz[mz-1][j][i].y = centz[mz-2][j][i].y;
889 centz[mz-1][j][i].z = centz[mz-2][j][i].z + gs[mz+1][j][i].z;
890 }
891 }
892 }
893 }
894
895 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
896 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
897
898 PROFILE_FUNCTION_END;
899 PetscFunctionReturn(0);
900}
901
902#undef __FUNCT__
903#define __FUNCT__ "ComputeFaceMetrics"
904/**
905 * @brief Computes the primary face metric components (Csi, Eta, Zet), including
906 * boundary extrapolation, and stores them in the corresponding global Vec
907 * members of the UserCtx structure (user->Csi, user->Eta, user->Zet).
908 *
909 * This is a self-contained routine that performs the following steps:
910 * 1. Obtains local ghosted nodal coordinates using DMGetCoordinatesLocal.
911 * 2. Calculates metrics for INTERIOR faces where finite difference stencils are valid.
912 * 3. EXTRAPOLATES metrics for faces on the physical domain boundaries by copying
913 * from the nearest computed interior face.
914 * 4. Assembles the global `user->Csi`, `user->Eta`, `user->Zet` Vecs.
915 * 5. Updates the local ghosted `user->lCsi`, `user->lEta`, `user->lZet` Vecs.
916 * * @details
917 * This function calculates the face area vectors, which are fundamental to the finite volume
918 * method on a curvilinear grid. The process is performed in two main stages to ensure
919 * robustness and correctness across the entire domain, including physical boundaries.
920 *
921 * **Stage 1: Interior Face Calculation**
922 * The core of the function uses centered finite-difference stencils on the nodal coordinates (`coor`)
923 * to compute the metric terms. For example, `csi` at a node `(k,j,i)` is calculated using a
924 * 2x2 stencil of nodes in the J-K plane.
925 *
926 * Crucially, these stencils are only valid for **interior nodes/faces** where all required neighboring
927 * nodes exist. The loops are intentionally constructed (e.g., `i_loop_start = 1`) to skip the
928 * calculation directly on the domain's physical boundaries (e.g., at i=0, j=0, etc.), as the
929 * stencil would require out-of-bounds data.
930 *
931 * **Stage 2: Boundary Face Extrapolation**
932 * After the interior metrics are computed, the values for the faces lying on the physical
933 * boundaries of the domain are populated. This is done via zero-order extrapolation, which
934 * simply means **copying the values from the nearest valid interior face layer.**
935 *
936 * For example, on a rank owning the global `i=0` boundary, `zet_arr[k][j][0]` is explicitly set to
937 * be equal to `zet_arr[k][j][1]`. Similarly, on a rank owning the `i=mx-1` boundary,
938 * `zet_arr[k][j][mx-1]` is set to `zet_arr[k][j][mx-2]`.
939 *
940 * This two-stage approach ensures that every node in the physical domain, including the
941 * boundaries, has a valid metric value associated with it.
942 *
943 * The function concludes by assembling the global PETSc vectors and updating the local ghosted
944 * versions, making the computed metrics ready for use by the solver.
945 *
946 * @param[in,out] user Pointer to the UserCtx structure.
947 *
948 * @return PetscErrorCode 0 on success.
949 *
950 * @note
951 * - This function is a complete "compute and make ready" unit for Csi, Eta, and Zet.
952 * - It's recommended to call `VecZeroEntries` on user->Csi, Eta, Zet before this
953 * if they might contain old data.
954 * - Csi,Eta,Zeta are face centered quantities which represent the surface area in magnitude and the direction of positive computational coordinate increase in direction.
955 */
956PetscErrorCode ComputeFaceMetrics(UserCtx *user)
957{
958 PetscErrorCode ierr;
959 DMDALocalInfo info;
960 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
961 Cmpnts ***nodal_coords_arr;
962 Vec localCoords_from_dm;
963
964 PetscFunctionBeginUser;
965
966 PROFILE_FUNCTION_BEGIN;
967
968 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
969
970 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
971
972 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
973 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
974 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
975 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
976
977 // --- 2. Get arrays for output global Vecs from UserCtx ---
978 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
979 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
980 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
981
982 // Define owned node ranges (global indices)
983 PetscInt xs = info.xs, xe = info.xs + info.xm;
984 PetscInt ys = info.ys, ye = info.ys + info.ym;
985 PetscInt zs = info.zs, ze = info.zs + info.zm;
986
987 // Global domain dimensions (total number of nodes)
988 PetscInt mx = info.mx;
989 PetscInt my = info.my;
990 PetscInt mz = info.mz;
991
992 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
993 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
994 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
995 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
996 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
997
998 // These represent the surface area of the curvilinear cell face and the normal rotated such that the direction of increasing coordinate is maintained.
999 // The metric vectors (Csi, Eta, Zet) are defined to point in the direction of their corresponding increasing computational coordinate.
1000
1001 // Calculate Csi
1002 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
1003 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
1004 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
1005
1006 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);
1007 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);
1008 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);
1009 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);
1010 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);
1011 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);
1012
1013 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
1014 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
1015 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
1016 }
1017 }
1018 }
1019
1020 // Calculate Eta
1021 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
1022 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
1023 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
1024
1025 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);
1026 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);
1027 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);
1028 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);
1029 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);
1030 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);
1031
1032 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
1033 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
1034 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
1035 }
1036 }
1037 }
1038
1039 // Calculate Zet
1040 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
1041 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
1042 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
1043
1044 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);
1045 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);
1046 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);
1047 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);
1048 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);
1049 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);
1050
1051 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
1052 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
1053 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
1054 }
1055 }
1056 }
1057
1058 // --- 4. Boundary Extrapolation ---
1059 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
1060 PetscInt i_bnd, j_bnd, k_bnd;
1061
1062 if (xs == 0) { // If this rank owns the global i=0 boundary
1063 i_bnd = 0;
1064 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1065 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1066 if (i_bnd + 1 < mx) {
1067 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
1068 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
1069 }
1070 }
1071 }
1072 }
1073 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
1074 i_bnd = mx - 1;
1075 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1076 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1077 if (i_bnd - 1 >= 0) {
1078 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
1079 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
1080 }
1081 }
1082 }
1083 }
1084 if (ys == 0) {
1085 j_bnd = 0;
1086 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1087 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1088 if (j_bnd + 1 < my) {
1089 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
1090 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
1091 }
1092 }
1093 }
1094 }
1095 if (ye == my) {
1096 j_bnd = my - 1;
1097 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1098 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1099 if (j_bnd - 1 >= 0) {
1100 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
1101 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
1102 }
1103 }
1104 }
1105 }
1106 if (zs == 0) {
1107 k_bnd = 0;
1108 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1109 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1110 if (k_bnd + 1 < mz) {
1111 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
1112 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
1113 }
1114 }
1115 }
1116 }
1117 if (ze == mz) {
1118 k_bnd = mz - 1;
1119 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1120 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1121 if (k_bnd - 1 >= 0) {
1122 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
1123 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
1124 }
1125 }
1126 }
1127 }
1128
1129 if (info.xs==0 && info.ys==0 && info.zs==0) {
1130 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
1131 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
1132 }
1133
1134 // --- 5. Restore all arrays ---
1135 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1136 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
1137 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
1138 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
1139
1140 // --- 6. Assemble Global Vectors ---
1141 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
1142 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
1143 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
1144 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
1145
1146 // --- 7. Update Local Ghosted Versions ---
1147 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
1148 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
1149 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
1150 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
1151
1152 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
1153
1154 PROFILE_FUNCTION_END;
1155
1156 PetscFunctionReturn(0);
1157}
1158
1159
1160#undef __FUNCT__
1161#define __FUNCT__ "ComputeCellCenteredJacobianInverse"
1162/**
1163 * @brief Calculates the cell-centered inverse Jacobian determinant (1/J), including
1164 * boundary extrapolation, stores it in `user->Aj`, assembles `user->Aj`, and
1165 * updates `user->lAj`.
1166 *
1167 * @param[in,out] user Pointer to the UserCtx structure.
1168 *
1169 * @return PetscErrorCode 0 on success.
1170 */
1171PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
1172{
1173 PetscErrorCode ierr;
1174 DMDALocalInfo info;
1175 PetscScalar ***aj_arr;
1176 Cmpnts ***nodal_coords_arr;
1177 Vec localCoords_from_dm;
1178
1179 PetscFunctionBeginUser;
1180 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
1181
1182 // --- 1. Get Nodal Coordinates and Output Array ---
1183 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
1184 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1185 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1186 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
1187
1188 // Define owned node ranges (global indices)
1189 PetscInt xs = info.xs, xe = info.xs + info.xm;
1190 PetscInt ys = info.ys, ye = info.ys + info.ym;
1191 PetscInt zs = info.zs, ze = info.zs + info.zm;
1192
1193 // Global domain dimensions (total number of nodes)
1194 PetscInt mx = info.mx;
1195 PetscInt my = info.my;
1196 PetscInt mz = info.mz;
1197
1198 // --- 2. Calculate Aj for INTERIOR Stencils ---
1199
1200 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
1201 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
1202 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
1203
1204 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
1205 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
1206 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
1207
1208 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
1209 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
1210 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
1211
1212 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) );
1213
1214 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) );
1215
1216 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) );
1217
1218 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) );
1219
1220 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) );
1221
1222 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) );
1223
1224 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) );
1225
1226 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) );
1227
1228 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) );
1229
1230 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);
1231 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
1232 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
1233 }
1234 }
1235 }
1236
1237 // --- 4. Boundary Extrapolation for Aj ---
1238 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
1239 PetscInt i_bnd, j_bnd, k_bnd;
1240
1241 if (xs == 0) {
1242 i_bnd = 0;
1243 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1244 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1245 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
1246 }
1247 }
1248 }
1249 if (xe == mx) {
1250 i_bnd = mx - 1;
1251 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1252 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1253 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
1254 }
1255 }
1256 }
1257 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
1258 if (ys == 0) {
1259 j_bnd = 0;
1260 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1261 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1262 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
1263 }
1264 }
1265 }
1266 if (ye == my) {
1267 j_bnd = my - 1;
1268 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1269 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1270 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
1271 }
1272 }
1273 }
1274 if (zs == 0) {
1275 k_bnd = 0;
1276 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1277 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1278 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
1279 }
1280 }
1281 }
1282 if (ze == mz) {
1283 k_bnd = mz - 1;
1284 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1285 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1286 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
1287 }
1288 }
1289 }
1290
1291 // --- 5. Restore arrays ---
1292 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1293 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
1294
1295 // --- 6. Assemble Global Vector ---
1296 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
1297 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
1298 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
1299
1300 // --- 7. Update Local Ghosted Version ---
1301 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
1302 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
1303
1304 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
1305 PetscFunctionReturn(0);
1306}
1307
1308
1309#undef __FUNCT__
1310#define __FUNCT__ "ComputeCellCentersAndSpacing"
1311/**
1312 * @brief Computes the physical location of cell centers and the spacing between them.
1313 *
1314 * This function calculates two key geometric properties from the nodal coordinates:
1315 * 1. `Cent`: A vector field storing the (x,y,z) coordinates of the center of each grid cell.
1316 * 2. `GridSpace`: A vector field storing the physical distance between adjacent
1317 * cell centers in the i, j, and k computational directions.
1318 *
1319 * It is a direct adaptation of the corresponding logic from the legacy `FormMetrics`.
1320 *
1321 * @param user The UserCtx for a specific grid level. The function populates `user->Cent` and `user->GridSpace`.
1322 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1323 *
1324 * @note Grid Spacing represents the vector connecting two adjacent grid cell faces. It is a vector whose magnitude gives the distance between two faces and direction is from one face center to the other.
1325 */
1326PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
1327{
1328 PetscErrorCode ierr;
1329 DMDALocalInfo info;
1330 Vec lCoords;
1331 const Cmpnts ***coor;
1332 Cmpnts ***cent, ***gs;
1333 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
1334 PetscInt xs,ys,zs,xe,ye,ze,mx,my,mz;
1335
1336 PetscFunctionBeginUser;
1337
1338 PROFILE_FUNCTION_BEGIN;
1339
1340 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);
1341
1342 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1343 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1344 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1345
1346 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1347 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1348
1349 xs = info.xs; xe = info.xs + info.xm;
1350 ys = info.ys; ye = info.ys + info.ym;
1351 zs = info.zs; ze = info.zs + info.zm;
1352 mx = info.mx; my = info.my; mz = info.mz;
1353
1354 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
1355 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
1356 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
1357
1358 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
1359 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
1360 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
1361
1362 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
1363 for (PetscInt k=k_start_node; k<k_end_node; k++) {
1364 for (PetscInt j=j_start_node; j<j_end_node; j++) {
1365 for (PetscInt i=i_start_node; i<i_end_node; i++) {
1366 // Calculate cell center as the average of its 8 corner nodes
1367 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);
1368 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);
1369 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);
1370
1371 // Calculate Grid Spacing in i-direction (distance between i-face centers)
1372 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);
1373 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);
1374 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);
1375 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);
1376 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);
1377 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);
1378 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1379
1380 // Calculate Grid Spacing in j-direction (distance between j-face centers)
1381 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);
1382 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);
1383 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);
1384 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);
1385 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);
1386 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);
1387 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1388
1389 // Calculate Grid Spacing in k-direction (distance between k-face centers)
1390 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);
1391 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);
1392 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);
1393 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);
1394 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);
1395 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);
1396 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1397 }
1398 }
1399 }
1400
1401 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1402 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1403 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1404
1405 // Assemble and update ghost regions for the new data
1406 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1407 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
1408 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
1409 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
1410
1411 ierr = ApplyPeriodicCorrectionsToCellCentersAndSpacing(user); CHKERRQ(ierr);
1412
1413 // Final assembly and ghost update after corrections
1414 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr);
1415 ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1416 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
1417
1418 PROFILE_FUNCTION_END;
1419
1420 PetscFunctionReturn(0);
1421}
1422
1423#undef __FUNCT__
1424#define __FUNCT__ "ComputeIFaceMetrics"
1425/**
1426 * @brief Computes metrics centered on constant-i faces (i-faces).
1427 *
1428 * This function calculates the metric terms (`ICsi`, `IEta`, `IZet`) and the
1429 * inverse Jacobian (`IAj`) located at the geometric center of each constant-i
1430 * face. This is a critical step for staggered-grid finite difference schemes.
1431 *
1432 * The process is a direct and faithful refactoring of the corresponding logic
1433 * from the legacy `FormMetrics` function:
1434 * 1. It first calculates the physical (x,y,z) coordinates of the center of
1435 * each i-face and stores them in the `user->Centx` vector.
1436 * 2. It then uses a boundary-aware, second-order finite difference stencil on
1437 * the `Centx` field to compute the derivatives (e.g., d(x)/d(csi)).
1438 * - Central differences are used in the grid interior.
1439 * - One-sided differences are used at the physical domain boundaries.
1440 * 3. Finally, these derivatives are used to compute the final metric terms and
1441 * the inverse Jacobian, which are stored in their respective `Vec` objects.
1442 *
1443 * @param user The UserCtx for a specific grid level. This function populates
1444 * the `user->ICsi`, `user->IEta`, `user->IZet`, and `user->IAj` vectors.
1445 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1446 *
1447 * @note Icsi,Ieta,Izet represent the computational coordinate gradients scaled by local volume inverse ((Iaj)
1448 *
1449 */
1450PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
1451{
1452 PetscErrorCode ierr;
1453 DMDALocalInfo info;
1454 Vec lCoords;
1455 const Cmpnts ***coor;
1456 Cmpnts ***centx; //***gs;
1457 const Cmpnts ***centx_const;
1458 Cmpnts ***icsi, ***ieta, ***izet;
1459 PetscScalar ***iaj;
1460 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1461
1462 PetscFunctionBeginUser;
1463
1464 PROFILE_FUNCTION_BEGIN;
1465
1466 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1467
1468 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1469 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1470 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1471 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1472 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1473 PetscInt gys = info.gys, gye = info.gys + info.gym;
1474 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1475
1476 PetscInt lxs = xs; PetscInt lxe = xe;
1477 PetscInt lys = ys; PetscInt lye = ye;
1478 PetscInt lzs = zs; PetscInt lze = ze;
1479
1480 if (xs==0) lxs = xs+1;
1481 if (ys==0) lys = ys+1;
1482 if (zs==0) lzs = zs+1;
1483
1484 if (xe==mx) lxe=xe-1;
1485 if (ye==my) lye=ye-1;
1486 if (ze==mz) lze=ze-1;
1487
1488 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1489 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1490 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1491 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1492 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1493
1494 //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);
1495
1496 // Loop over the ghosted region to calculate all local face centers
1497 // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
1498 PetscInt j_end = (ye == mx)? my - 1:gye;
1499 PetscInt k_end = (ze == mz)? mz - 1:gze;
1500
1501 for (PetscInt k = gzs + 1; k < k_end; k++) {
1502 for (PetscInt j = gys + 1; j < j_end; j++) {
1503 for (PetscInt i = gxs; i < gxe; i++) {
1504 //----- DEBUG ------
1505 //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);
1506 //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,
1507 // coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
1508 // coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
1509 // coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
1510 // coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
1511
1512 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);
1513 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);
1514 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);
1515
1516 //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);
1517 }
1518 }
1519 }
1520
1521 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
1522 /*
1523 if(xs==0){
1524 for(PetscInt k=gzs+1;k < gze; k++){
1525 for(PetscInt j=gys+1;j < gye; j++){
1526 PetscInt i=0;
1527 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
1528 centx[k][j][i-1].y=centx[k][j][i].y;
1529 centx[k][j][i-1].z=centx[k][j][i].z;
1530 }
1531 }
1532 }
1533 if (xe==mx){
1534 for(PetscInt k=gzs+1; k<gze; k++) {
1535 for (PetscInt j=gys+1; j<gye;j++) {
1536 PetscInt i=mx-1;
1537 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
1538 centx[k][j][i].y=centx[k][j][i-1].y;
1539 centx[k][j][i].z=centx[k][j][i-1].z;
1540 }
1541 }
1542 }
1543 */
1544
1545 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1546 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1547
1548 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1549
1550 // For periodic BCs, exchange ghosts between processes
1551 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
1552 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
1553 ierr = DMLocalToLocalBegin(user->fda, user->Centx, INSERT_VALUES, user->Centx); CHKERRQ(ierr);
1554 ierr = DMLocalToLocalEnd(user->fda, user->Centx, INSERT_VALUES, user->Centx); CHKERRQ(ierr);
1555 }
1556
1557 ierr = ApplyPeriodicCorrectionsToIFaceCenter(user); CHKERRQ(ierr);
1558
1559 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1560
1561 // --- Part 2: Calculate metrics using face-centered coordinates ---
1562 ierr = DMDAVecGetArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1563 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1564 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1565 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1566 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1567
1568 // Loop over the OWNED region where we will store the final metrics
1569 for (PetscInt k=lzs; k<lze; k++) {
1570 for (PetscInt j=lys; j<lye; j++) {
1571 for (PetscInt i=xs; i<lxe; i++) {
1572
1573 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1574 if(i == 0 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC){
1575 // Forward difference at the domain's min-i boundary
1576 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
1577 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
1578 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
1579 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1580 // Backward difference at the domain's max-i boundary
1581 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
1582 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
1583 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
1584 } else { // Central difference in the interior (or if PERIODIC BCs)
1585 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
1586 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
1587 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
1588 }
1589
1590 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1591 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1592 // Forward difference
1593 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
1594 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
1595 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
1596 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1597 // Backward difference
1598 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
1599 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
1600 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
1601 } else { // Central difference (interior or PERIODIC)
1602 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
1603 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
1604 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
1605 }
1606
1607 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1608 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1609 // Forward difference
1610 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
1611 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
1612 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
1613 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1614 // Backward difference
1615 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
1616 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
1617 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
1618 } else { // Central difference (Interior + PERIODIC)
1619 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
1620 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
1621 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
1622 }
1623
1624 // --- Metric calculations (identical to legacy FormMetrics) ---
1625 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1626 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1627 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1628
1629 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1630 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1631 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1632
1633 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
1634 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1635 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
1636
1637 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
1638 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
1639 iaj[k][j][i] = 1.0 / iaj[k][j][i];
1640 }
1641 }
1642 }
1643 }
1644
1645 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1646 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1647 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1648 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1649 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1650
1651 // --- Part 3: Assemble global vectors and update local ghosts ---
1652 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
1653 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
1654 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
1655 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
1656
1657 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
1658 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
1659 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
1660 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
1661
1662 PROFILE_FUNCTION_END;
1663
1664 PetscFunctionReturn(0);
1665}
1666
1667#undef __FUNCT__
1668#define __FUNCT__ "ComputeJFaceMetrics"
1669/**
1670 * @brief Computes metrics centered on constant-j faces (j-faces).
1671 *
1672 * This function calculates the metric terms (`JCsi`, `JEta`, `JZet`) and the
1673 * inverse Jacobian (`JAj`) located at the geometric center of each constant-j
1674 * face. This is a critical step for staggered-grid finite difference schemes.
1675 *
1676 * The process is a direct and faithful refactoring of the corresponding logic
1677 * from the legacy `FormMetrics` function:
1678 * 1. It first calculates the physical (x,y,z) coordinates of the center of
1679 * each i-face and stores them in the `user->Centy` vector.
1680 * 2. It then uses a boundary-aware, second-order finite difference stencil on
1681 * the `Centy` field to compute the derivatives (e.g., d(x)/d(csi)).
1682 * - Central differences are used in the grid interior.
1683 * - One-sided differences are used at the physical domain boundaries.
1684 * 3. Finally, these derivatives are used to compute the final metric terms and
1685 * the inverse Jacobian, which are stored in their respective `Vec` objects.
1686 *
1687 * @param user The UserCtx for a specific grid level. This function populates
1688 * the `user->JCsi`, `user->JEta`, `user->JZet`, and `user->JAj` vectors.
1689 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1690 * @note Jcsi,Jeta,Jzet represent the computational coordinate gradients scaled by local volume inverse ((Jaj)
1691 */
1692PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
1693{
1694 PetscErrorCode ierr;
1695 DMDALocalInfo info;
1696 Vec lCoords;
1697 const Cmpnts ***coor;
1698 Cmpnts ***centy; //***gs;
1699 const Cmpnts ***centy_const;
1700 Cmpnts ***jcsi, ***jeta, ***jzet;
1701 PetscScalar ***jaj;
1702 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1703
1704 PetscFunctionBeginUser;
1705
1706 PROFILE_FUNCTION_BEGIN;
1707
1708 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1709
1710 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1711 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1712 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1713 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1714 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1715 PetscInt gys = info.gys, gye = info.gys + info.gym;
1716 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1717
1718 PetscInt lxs = xs; PetscInt lxe = xe;
1719 PetscInt lys = ys; PetscInt lye = ye;
1720 PetscInt lzs = zs; PetscInt lze = ze;
1721
1722 if (xs==0) lxs = xs+1;
1723 if (ys==0) lys = ys+1;
1724 if (zs==0) lzs = zs+1;
1725
1726 if (xe==mx) lxe=xe-1;
1727 if (ye==my) lye=ye-1;
1728 if (ze==mz) lze=ze-1;
1729
1730 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1731 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1732 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1733 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1734 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1735
1736 // Loop over the ghosted region to calculate all local face centers
1737 // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
1738 PetscInt k_end = (ze == mz)? mz - 1:gze;
1739 PetscInt i_end = (xe == mx)? mx - 1:gxe;
1740
1741 for (PetscInt k = gzs + 1; k < k_end; k++) {
1742 for (PetscInt j = gys; j < gye; j++) {
1743 for (PetscInt i = gxs + 1; i < i_end; i++) {
1744 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);
1745 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);
1746 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);
1747 }
1748 }
1749 }
1750
1751 /*
1752 if(ys==0){
1753 for(PetscInt k=gzs+1;k < gze; k++){
1754 for(PetscInt i=gxs+1;j < gxe; i++){
1755 PetscInt j=0;
1756 centy[k][j-1][i].x=centy[k][j][i].x;
1757 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1758 centy[k][j-1][i].z=centy[k][j][i].z;
1759 }
1760 }
1761 }
1762 if (ye==my){
1763 for(PetscInt k=gzs+1; k<gze; k++) {
1764 for (PetscInt i=gxs+1; j<gxe;i++) {
1765 PetscInt j=my-1;
1766 centy[k][j][i].x=centy[k][j-1][i].x
1767 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1768 centy[k][j][i].z=centy[k][j-1][i].z;
1769 }
1770 }
1771 }
1772 */
1773
1774 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1775 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1776 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1777
1778 // For periodic BCs, exchange ghosts between processes
1779 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
1780 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
1781 ierr = DMLocalToLocalBegin(user->fda, user->Centy, INSERT_VALUES, user->Centy); CHKERRQ(ierr);
1782 ierr = DMLocalToLocalEnd(user->fda, user->Centy, INSERT_VALUES, user->Centy); CHKERRQ(ierr);
1783 }
1784
1785 ierr = ApplyPeriodicCorrectionsToJFaceCenter(user); CHKERRQ(ierr);
1786
1787 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1788
1789 // --- Part 2: Calculate metrics using face-centered coordinates ---
1790 ierr = DMDAVecGetArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1791 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1792 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1793 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1794 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1795
1796 // Loop over the OWNED region where we will store the final metrics
1797 for (PetscInt k=lzs; k<lze; k++) {
1798 for (PetscInt j=ys; j<lye; j++) {
1799 for (PetscInt i=lxs; i<lxe; i++) {
1800
1801 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1802 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1803 // Forward difference at the domain's min-i boundary
1804 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1805 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1806 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1807 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1808 // Backward difference at the domain's max-i boundary
1809 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1810 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1811 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1812 } else { // Central difference in the interior or PERIODIC
1813 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1814 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1815 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1816 }
1817
1818 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1819 if (j == 0 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1820 // Forward difference
1821 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1822 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1823 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1824 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1825 // Backward difference
1826 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1827 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1828 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1829 } else { // Central difference (interior or PERIODIC)
1830 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1831 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1832 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1833 }
1834
1835 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1836 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1837 // Forward difference
1838 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1839 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1840 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1841 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1842 // Backward difference
1843 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1844 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1845 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1846 } else { // Central difference (Interior or PERIODIC)
1847 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1848 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1849 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1850 }
1851
1852 // --- Metric calculations (identical to legacy FormMetrics) ---
1853 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1854 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1855 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1856
1857 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1858 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1859 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1860
1861 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1862 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1863 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1864
1865 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1866 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1867 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1868 }
1869 }
1870 }
1871 }
1872
1873 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1874 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1875 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1876 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1877 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1878
1879 // --- Part 3: Assemble global vectors and update local ghosts ---
1880 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1881 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1882 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1883 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1884
1885 ierr = UpdateLocalGhosts(user, "JCsi"); CHKERRQ(ierr);
1886 ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
1887 ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
1888 ierr = UpdateLocalGhosts(user, "JAj"); CHKERRQ(ierr);
1889
1890 PROFILE_FUNCTION_END;
1891
1892 PetscFunctionReturn(0);
1893}
1894
1895#undef __FUNCT__
1896#define __FUNCT__ "ComputeJFaceMetrics"
1897/**
1898 * @brief Computes metrics centered on constant-k faces (k-faces).
1899 *
1900 * This function calculates the metric terms (`KCsi`, `KEta`, `KZet`) and the
1901 * inverse Jacobian (`KAj`) located at the geometric center of each constant-k
1902 * face. This is a critical step for staggered-grid finite difference schemes.
1903 *
1904 * The process is a direct and faithful refactoring of the corresponding logic
1905 * from the legacy `FormMetrics` function:
1906 * 1. It first calculates the physical (x,y,z) coordinates of the center of
1907 * each i-face and stores them in the `user->Centz` vector.
1908 * 2. It then uses a boundary-aware, second-order finite difference stencil on
1909 * the `Centz` field to compute the derivatives (e.g., d(x)/d(csi)).
1910 * - Central differences are used in the grid interior.
1911 * - One-sided differences are used at the physical domain boundaries.
1912 * 3. Finally, these derivatives are used to compute the final metric terms and
1913 * the inverse Jacobian, which are stored in their respective `Vec` objects.
1914 *
1915 * @param user The UserCtx for a specific grid level. This function populates
1916 * the `user->KCsi`, `user->KEta`, `user->KZet`, and `user->KAj` vectors.
1917 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1918 * @note Kcsi,Keta,Kzet represent the computational coordinate gradients scaled by local volume inverse ((Kaj)
1919 */
1920PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
1921{
1922 PetscErrorCode ierr;
1923 DMDALocalInfo info;
1924 Vec lCoords;
1925 const Cmpnts ***coor;
1926 Cmpnts ***centz; //***gs;
1927 const Cmpnts ***centz_const;
1928 Cmpnts ***kcsi, ***keta, ***kzet;
1929 PetscScalar ***kaj;
1930 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1931
1932 PetscFunctionBeginUser;
1933
1934 PROFILE_FUNCTION_BEGIN;
1935
1936 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1937
1938 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1939 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1940 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1941 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1942 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1943 PetscInt gys = info.gys, gye = info.gys + info.gym;
1944 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1945
1946 PetscInt lxs = xs; PetscInt lxe = xe;
1947 PetscInt lys = ys; PetscInt lye = ye;
1948 PetscInt lzs = zs; PetscInt lze = ze;
1949
1950 if (xs==0) lxs = xs+1;
1951 if (ys==0) lys = ys+1;
1952 if (zs==0) lzs = zs+1;
1953
1954 if (xe==mx) lxe=xe-1;
1955 if (ye==my) lye=ye-1;
1956 if (ze==mz) lze=ze-1;
1957
1958 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1959 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1960 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1961 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1962 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1963
1964 // Loop over the ghosted region to calculate all local face centers
1965 // To ensure we don't mistakenly calculate unused/dummy elements along non-dominant directions.
1966 PetscInt i_end = (xe == mx)? mx - 1:gxe;
1967 PetscInt j_end = (ye == my)? my - 1:gye;
1968 for (PetscInt k = gzs; k < gze; k++) {
1969 for (PetscInt j = gys; j < gye; j++) {
1970 for (PetscInt i = gxs + 1; i < gxe; i++) {
1971 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);
1972 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);
1973 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);
1974 }
1975 }
1976 }
1977
1978 /*
1979 if(zs==0){
1980 for(PetscInt j=gys+1;j < gye; j++){
1981 for(PetscInt i=gxs+1;j < gxe; i++){
1982 PetscInt k=0;
1983 centz[k-1][j][i].x=centz[k][j][i].x;
1984 centz[k-1][j][i].y=centz[k][j][i].y;
1985 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1986 }
1987 }
1988 }
1989 if (ze==mz){
1990 for(PetscInt j=gys+1; j<gye; j++) {
1991 for (PetscInt i=gxs+1; j<gxe;i++) {
1992 PetscInt k=mz-1;
1993 centy[k][j][i].x=centy[k-1][j][i].x
1994 centy[k][j][i].y=centy[k-1][j][i].y;
1995 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1996 }
1997 }
1998 }
1999 */
2000
2001 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
2002 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
2003 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
2004
2005 // For periodic BCs, exchange ghosts between processes
2006 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
2007 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
2008 ierr = DMLocalToLocalBegin(user->fda, user->Centz, INSERT_VALUES, user->Centz); CHKERRQ(ierr);
2009 ierr = DMLocalToLocalEnd(user->fda, user->Centz, INSERT_VALUES, user->Centz); CHKERRQ(ierr);
2010 }
2011
2012 ierr = ApplyPeriodicCorrectionsToKFaceCenter(user); CHKERRQ(ierr);
2013
2014 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
2015
2016 // --- Part 2: Calculate metrics using face-centered coordinates ---
2017 ierr = DMDAVecGetArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
2018 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
2019 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
2020 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
2021 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
2022
2023 // Loop over the OWNED region where we will store the final metrics
2024 for (PetscInt k=zs; k<lze; k++) {
2025 for (PetscInt j=lys; j<lye; j++) {
2026 for (PetscInt i=lxs; i<lxe; i++) {
2027
2028 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
2029 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
2030 // Forward difference at the domain's min-i boundary
2031 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
2032 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
2033 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
2034 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
2035 // Backward difference at the domain's max-i boundary
2036 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
2037 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
2038 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
2039 } else { // Central difference in the interior (or PERIODIC)
2040 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
2041 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
2042 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
2043 }
2044
2045 // --- Stencil Logic for d/deta (derivative in j-direction) ---
2046 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
2047 // Forward difference
2048 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
2049 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
2050 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
2051 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
2052 // Backward difference
2053 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
2054 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
2055 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
2056 } else { // Central difference (interior or PERIODIC)
2057 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
2058 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
2059 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
2060 }
2061
2062 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
2063 if (k == 0 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
2064 // Forward difference
2065 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
2066 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
2067 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
2068 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
2069 // Backward difference
2070 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
2071 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
2072 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
2073 } else { // Central difference (Interior or PERIODIC)
2074 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
2075 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
2076 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
2077 }
2078
2079 // --- Metric calculations (identical to legacy FormMetrics) ---
2080 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
2081 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
2082 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
2083
2084 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
2085 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
2086 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
2087
2088 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
2089 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
2090 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
2091
2092 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
2093 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
2094 kaj[k][j][i] = 1.0 / kaj[k][j][i];
2095 }
2096 }
2097 }
2098 }
2099
2100 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
2101 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
2102 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
2103 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
2104 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
2105
2106 // --- Part 3: Assemble global vectors and update local ghosts ---
2107 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
2108 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
2109 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
2110 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
2111
2112 ierr = UpdateLocalGhosts(user, "KCsi"); CHKERRQ(ierr);
2113 ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
2114 ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
2115 ierr = UpdateLocalGhosts(user, "KAj"); CHKERRQ(ierr);
2116
2117 PROFILE_FUNCTION_END;
2118
2119 PetscFunctionReturn(0);
2120}
2121
2122static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
2123{
2124 PetscInt nidx;
2125 DMDALocalInfo info = user->info;
2126
2127 PetscInt mx = info.mx, my = info.my;
2128
2129 AO ao;
2130 DMDAGetAO(user->da, &ao);
2131 nidx=i+j*mx+k*mx*my;
2132
2133 AOApplicationToPetsc(ao,1,&nidx);
2134
2135 return (nidx);
2136}
2137
2138#undef __FUNCT__
2139#define __FUNCT__ "ComputeMetricsDivergence"
2140/**
2141 * @brief Computes the divergence of the grid metrics and identifies the maximum value.
2142 *
2143 * This function serves as a diagnostic tool to assess the quality of the grid
2144 * metrics. It calculates the divergence of the face metrics (Csi, Eta, Zet)
2145 * and scales it by the inverse of the cell Jacobian. The maximum divergence
2146 * value is then located, and its grid coordinates are printed to the console,
2147 * helping to pinpoint areas of potential grid quality issues.
2148 *
2149 * @param user The UserCtx, containing all necessary grid data.
2150 * @return PetscErrorCode
2151 */
2152PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
2153{
2154 DM da = user->da, fda = user->fda;
2155 DMDALocalInfo info = user->info;
2156 PetscInt xs = info.xs, xe = info.xs + info.xm;
2157 PetscInt ys = info.ys, ye = info.ys + info.ym;
2158 PetscInt zs = info.zs, ze = info.zs + info.zm;
2159 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2160 PetscInt lxs, lys, lzs, lxe, lye, lze;
2161 PetscInt i, j, k;
2162 Vec Div;
2163 PetscReal ***div, ***aj;
2164 Cmpnts ***csi, ***eta, ***zet;
2165 PetscReal maxdiv;
2166
2167 PetscFunctionBeginUser;
2168
2169 PROFILE_FUNCTION_BEGIN;
2170
2171 lxs = xs; lxe = xe;
2172 lys = ys; lye = ye;
2173 lzs = zs; lze = ze;
2174
2175 if (xs == 0) lxs = xs + 1;
2176 if (ys == 0) lys = ys + 1;
2177 if (zs == 0) lzs = zs + 1;
2178
2179 if (xe == mx) lxe = xe - 1;
2180 if (ye == my) lye = ye - 1;
2181 if (ze == mz) lze = ze - 1;
2182
2183 DMDAVecGetArray(fda, user->lCsi, &csi);
2184 DMDAVecGetArray(fda, user->lEta, &eta);
2185 DMDAVecGetArray(fda, user->lZet, &zet);
2186 DMDAVecGetArray(da, user->lAj, &aj);
2187
2188 VecDuplicate(user->P, &Div);
2189 VecSet(Div, 0.);
2190 DMDAVecGetArray(da, Div, &div);
2191
2192 for (k = lzs; k < lze; k++) {
2193 for (j = lys; j < lye; j++) {
2194 for (i = lxs; i < lxe; i++) {
2195 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
2196 eta[k][j][i].x - eta[k][j-1][i].x +
2197 zet[k][j][i].x - zet[k-1][j][i].x +
2198 csi[k][j][i].y - csi[k][j][i-1].y +
2199 eta[k][j][i].y - eta[k][j-1][i].y +
2200 zet[k][j][i].y - zet[k-1][j][i].y +
2201 csi[k][j][i].z - csi[k][j][i-1].z +
2202 eta[k][j][i].z - eta[k][j-1][i].z +
2203 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
2204 div[k][j][i] = fabs(divergence);
2205 }
2206 }
2207 }
2208
2209 DMDAVecRestoreArray(da, Div, &div);
2210
2211 PetscReal MaxFlatIndex;
2212 VecMax(Div, &MaxFlatIndex, &maxdiv);
2213 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %d.\n",maxdiv,MaxFlatIndex);
2214
2215 for (k=zs; k<ze; k++) {
2216 for (j=ys; j<ye; j++) {
2217 for (i=xs; i<xe; i++) {
2218 if (Gidx(i,j,k,user) == MaxFlatIndex) {
2219 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,(int)k,(int)j,(int)i);
2220 }
2221 }
2222 }
2223 }
2224
2225
2226 DMDAVecRestoreArray(fda, user->lCsi, &csi);
2227 DMDAVecRestoreArray(fda, user->lEta, &eta);
2228 DMDAVecRestoreArray(fda, user->lZet, &zet);
2229 DMDAVecRestoreArray(da, user->lAj, &aj);
2230 VecDestroy(&Div);
2231
2232
2233 PROFILE_FUNCTION_END;
2234
2235 PetscFunctionReturn(0);
2236}
2237
2238#undef __FUNCT__
2239#define __FUNCT__ "ComputeMetricNorms"
2240/**
2241 * @brief Computes the max-min values of the grid metrics.
2242 *
2243 * This function serves as a diagnostic tool to assess the quality of the grid
2244 * metrics. It calculates the bounds of the face metrics (Csi, Eta, Zet).
2245 *
2246 * @param user The UserCtx, containing all necessary grid data.
2247 * @return PetscErrorCode
2248 */
2249PetscErrorCode ComputeMetricNorms(UserCtx *user)
2250{
2251
2252 DMDALocalInfo info = user->info;
2253 PetscInt xs = info.xs, xe = info.xs + info.xm;
2254 PetscInt ys = info.ys, ye = info.ys + info.ym;
2255 PetscInt zs = info.zs, ze = info.zs + info.zm;
2256 PetscInt i, j, k;
2257
2258 PetscFunctionBeginUser;
2259
2260 PROFILE_FUNCTION_BEGIN;
2261
2262 PetscReal CsiMax, EtaMax, ZetMax;
2263 PetscReal ICsiMax, IEtaMax, IZetMax;
2264 PetscReal JCsiMax, JEtaMax, JZetMax;
2265 PetscReal KCsiMax, KEtaMax, KZetMax;
2266 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
2267
2268 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
2269 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
2270 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
2271 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
2272 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
2273
2274 // Max Values
2275 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
2276 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
2277 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
2278
2279 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
2280 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
2281 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
2282
2283 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
2284 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
2285 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
2286
2287 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
2288 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
2289 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
2290
2291 VecMax(user->lAj,&AjMaxArg,&AjMax);
2292 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
2293 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
2294 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
2295
2296 VecMax(user->lAj,&AjMaxArg,&AjMax);
2297 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
2298 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
2299 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
2300
2301 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
2302
2303 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
2304 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
2305 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
2306 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
2307 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
2308
2309 for (k=zs; k<ze; k++) {
2310 for (j=ys; j<ye; j++) {
2311 for (i=xs; i<xe; i++) {
2312 if (Gidx(i,j,k,user) == CsiMaxArg) {
2313 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
2314 }
2315 if (Gidx(i,j,k,user) == EtaMaxArg) {
2316 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
2317 }
2318 if (Gidx(i,j,k,user) == ZetMaxArg) {
2319 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
2320 }
2321 if (Gidx(i,j,k,user) == ICsiMaxArg) {
2322 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
2323 }
2324 if (Gidx(i,j,k,user) == IEtaMaxArg) {
2325 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
2326 }
2327 if (Gidx(i,j,k,user) == IZetMaxArg) {
2328 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
2329 }
2330 if (Gidx(i,j,k,user) == JCsiMaxArg) {
2331 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
2332 }
2333 if (Gidx(i,j,k,user) == JEtaMaxArg) {
2334 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
2335 }
2336 if (Gidx(i,j,k,user) == JZetMaxArg) {
2337 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
2338 }
2339 if (Gidx(i,j,k,user) == KCsiMaxArg) {
2340 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
2341 }
2342 if (Gidx(i,j,k,user) == KEtaMaxArg) {
2343 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
2344 }
2345 if (Gidx(i,j,k,user) == KZetMaxArg) {
2346 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
2347 }
2348 if (Gidx(i,j,k,user) == AjMaxArg) {
2349 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
2350 }
2351 if (Gidx(i,j,k,user) == IAjMaxArg) {
2352 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
2353 }
2354 if (Gidx(i,j,k,user) == JAjMaxArg) {
2355 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
2356 }
2357 if (Gidx(i,j,k,user) == KAjMaxArg) {
2358 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
2359 }
2360 }
2361 }
2362 }
2363
2364 /*
2365 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
2366 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
2367 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
2368 */
2369
2370 PROFILE_FUNCTION_END;
2371
2372 PetscFunctionReturn(0);
2373}
2374
2375#undef __FUNCT__
2376#define __FUNCT__ "ComputeAllGridMetrics"
2377/**
2378 * @brief Orchestrates the calculation of all grid metrics.
2379 *
2380 * This function iterates through every UserCtx in the multigrid and multi-block
2381 * hierarchy. For each context, it calls a series of modern, modular helper
2382 * functions to compute the face metrics (Csi, Eta, Zet), the cell-centered
2383 * inverse Jacobian (Aj), and to validate the grid's orientation.
2384 *
2385 * @param simCtx The master SimCtx, containing the configured UserCtx hierarchy.
2386 * @return PetscErrorCode
2387 */
2388PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
2389{
2390 PetscErrorCode ierr;
2391 UserMG *usermg = &simCtx->usermg;
2392 MGCtx *mgctx = usermg->mgctx;
2393 PetscInt nblk = simCtx->block_number;
2394
2395 PetscFunctionBeginUser;
2396
2397 PROFILE_FUNCTION_BEGIN;
2398
2399 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
2400
2401 // Loop through all levels and all blocks
2402 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
2403 for (PetscInt bi = 0; bi < nblk; bi++) {
2404 UserCtx *user = &mgctx[level].user[bi];
2405 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
2406
2407 // Call the modern, modular helper functions for each UserCtx.
2408 // These functions are self-contained and operate on the data within the provided context.
2409 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
2410 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
2411 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
2412 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
2413 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
2414 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
2415 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
2416
2417 // Apply Periodic Boundary Condition Adjustments if necessary
2418 ierr = ApplyMetricsPeriodicBCs(user); CHKERRQ(ierr);
2419 // Diagnostics
2420 ierr = ComputeMetricNorms(user);
2421 if (level == usermg->mglevels - 1) {
2422 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
2423 }
2424 }
2425 }
2426
2427 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
2428
2429 PROFILE_FUNCTION_END;
2430
2431 PetscFunctionReturn(0);
2432}
2433/* ------------------------------------------------------------------------- */
2434/* End of Metric.c */
2435
ApplyMetricsPeriodicBCs
PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
(Orchestrator) Updates all metric-related fields in the local ghost cell regions for periodic boundar...
Definition Boundaries.c:2099
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Orchestrates the calculation of all grid metrics.
Definition Metric.c:2388
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).
Definition Metric.c:110
InvertCovariantMetricTensor
PetscErrorCode InvertCovariantMetricTensor(double covariantTensor[3][3], double contravariantTensor[3][3])
Inverts the 3x3 covariant metric tensor to obtain the contravariant metric tensor.
Definition Metric.c:212
ComputeMetricNorms
PetscErrorCode ComputeMetricNorms(UserCtx *user)
Computes the max-min values of the grid metrics.
Definition Metric.c:2249
ApplyPeriodicCorrectionsToCellCentersAndSpacing
PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing(UserCtx *user)
Applies periodic boundary corrections to cell centers (Cent) and grid spacing (GridSpace).
Definition Metric.c:452
TrilinearBlend
static void TrilinearBlend(const Cmpnts V[8], PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Map logical coordinates to physical space using trilinear basis.
Definition Metric.c:54
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition Metric.c:2122
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain...
Definition Metric.c:368
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Computes the physical location of cell centers and the spacing between them.
Definition Metric.c:1326
MetricGetCellVertices
PetscErrorCode MetricGetCellVertices(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
Definition Metric.c:27
ComputeJFaceMetrics
PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
Computes metrics centered on constant-j faces (j-faces).
Definition Metric.c:1692
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Computes the divergence of the grid metrics and identifies the maximum value.
Definition Metric.c:2152
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation,...
Definition Metric.c:956
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.
Definition Metric.c:313
CalculateFaceNormalAndArea
PetscErrorCode CalculateFaceNormalAndArea(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
Computes the unit normal vectors and areas of the three faces of a computational cell.
Definition Metric.c:257
ApplyPeriodicCorrectionsToIFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter(UserCtx *user)
Applies periodic boundary corrections to i-face centers (Centx).
Definition Metric.c:656
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Calculates the cell-centered inverse Jacobian determinant (1/J), including boundary extrapolation,...
Definition Metric.c:1171
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Computes metrics centered on constant-k faces (k-faces).
Definition Metric.c:1920
ApplyPeriodicCorrectionsToJFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
Applies periodic boundary corrections to j-face centers (Centy).
Definition Metric.c:742
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 Metric.c:77
ComputeIFaceMetrics
PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
Computes metrics centered on constant-i faces (i-faces).
Definition Metric.c:1450
ApplyPeriodicCorrectionsToKFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
Applies periodic boundary corrections to k-face centers (Centz).
Definition Metric.c:828
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 Metric.c:167
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:312
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:267
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:746
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:33
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:34
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:737
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1157
UserCtx::lCent
Vec lCent
Definition variables.h:776
PERIODIC
@ PERIODIC
Definition variables.h:219
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:776
UserCtx::JCsi
Vec JCsi
Definition variables.h:779
UserCtx::KAj
Vec KAj
Definition variables.h:780
MGCtx::user
UserCtx * user
Definition variables.h:474
UserCtx::JEta
Vec JEta
Definition variables.h:779
UserCtx::Zet
Vec Zet
Definition variables.h:776
SimCtx::rank
PetscMPIInt rank
Definition variables.h:588
UserCtx::cgrid
PetscInt cgrid
Definition variables.h:743
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:746
SimCtx::block_number
PetscInt block_number
Definition variables.h:649
UserCtx::lIEta
Vec lIEta
Definition variables.h:778
UserCtx::lIZet
Vec lIZet
Definition variables.h:778
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:731
UserCtx::IZet
Vec IZet
Definition variables.h:778
UserCtx::Centz
Vec Centz
Definition variables.h:777
UserCtx::IEta
Vec IEta
Definition variables.h:778
UserCtx::lZet
Vec lZet
Definition variables.h:776
SimCtx::usermg
UserMG usermg
Definition variables.h:698
UserCtx::Csi
Vec Csi
Definition variables.h:776
UserCtx::lIAj
Vec lIAj
Definition variables.h:778
UserCtx::_this
PetscInt _this
Definition variables.h:741
UserCtx::lKEta
Vec lKEta
Definition variables.h:780
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:779
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::JZet
Vec JZet
Definition variables.h:779
UserCtx::Centx
Vec Centx
Definition variables.h:777
UserCtx::lKZet
Vec lKZet
Definition variables.h:780
UserCtx::Eta
Vec Eta
Definition variables.h:776
UserCtx::lJEta
Vec lJEta
Definition variables.h:779
UserCtx::lCsi
Vec lCsi
Definition variables.h:776
UserCtx::lGridSpace
Vec lGridSpace
Definition variables.h:776
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:792
UserCtx::ICsi
Vec ICsi
Definition variables.h:778
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:780
UserMG::mglevels
PetscInt mglevels
Definition variables.h:481
UserCtx::lJZet
Vec lJZet
Definition variables.h:779
UserCtx::IAj
Vec IAj
Definition variables.h:778
UserCtx::JAj
Vec JAj
Definition variables.h:779
UserCtx::KEta
Vec KEta
Definition variables.h:780
UserCtx::lAj
Vec lAj
Definition variables.h:776
UserCtx::lICsi
Vec lICsi
Definition variables.h:778
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:741
UserCtx::info
DMDALocalInfo info
Definition variables.h:735
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:776
UserCtx::KZet
Vec KZet
Definition variables.h:780
UserCtx::Cent
Vec Cent
Definition variables.h:776
UserCtx::KCsi
Vec KCsi
Definition variables.h:780
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:484
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:295
UserCtx::Centy
Vec Centy
Definition variables.h:777
UserCtx::lJAj
Vec lJAj
Definition variables.h:779
UserCtx::lKAj
Vec lKAj
Definition variables.h:780
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:204
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:206
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:205
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:206
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:204
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:205
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
MGCtx
Context for Multigrid operations.
Definition variables.h:473
SimCtx
The master context for the entire simulation.
Definition variables.h:585
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:728
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:480