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__ "ComputeFaceMetrics"
442/**
443 * @brief Computes the primary face metric components (Csi, Eta, Zet), including
444 * boundary extrapolation, and stores them in the corresponding global Vec
445 * members of the UserCtx structure (user->Csi, user->Eta, user->Zet).
446 *
447 * This is a self-contained routine that performs the following steps:
448 * 1. Obtains local ghosted nodal coordinates using DMGetCoordinatesLocal.
449 * 2. Calculates metrics for INTERIOR faces where finite difference stencils are valid.
450 * 3. EXTRAPOLATES metrics for faces on the physical domain boundaries by copying
451 * from the nearest computed interior face.
452 * 4. Assembles the global `user->Csi`, `user->Eta`, `user->Zet` Vecs.
453 * 5. Updates the local ghosted `user->lCsi`, `user->lEta`, `user->lZet` Vecs.
454 *
455 * @param[in,out] user Pointer to the UserCtx structure.
456 *
457 * @return PetscErrorCode 0 on success.
458 *
459 * @note
460 * - This function is a complete "compute and make ready" unit for Csi, Eta, and Zet.
461 * - It's recommended to call `VecZeroEntries` on user->Csi, Eta, Zet before this
462 * if they might contain old data.
463 */
464PetscErrorCode ComputeFaceMetrics(UserCtx *user)
465{
466 PetscErrorCode ierr;
467 DMDALocalInfo info;
468 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
469 Cmpnts ***nodal_coords_arr;
470 Vec localCoords_from_dm;
471
472 PetscFunctionBeginUser;
473
474 PROFILE_FUNCTION_BEGIN;
475
476 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
477
478 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
479
480 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
481 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
482 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
483 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
484
485 // --- 2. Get arrays for output global Vecs from UserCtx ---
486 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
487 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
488 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
489
490 // Define owned node ranges (global indices)
491 PetscInt xs = info.xs, xe = info.xs + info.xm;
492 PetscInt ys = info.ys, ye = info.ys + info.ym;
493 PetscInt zs = info.zs, ze = info.zs + info.zm;
494
495 // Global domain dimensions (total number of nodes)
496 PetscInt mx = info.mx;
497 PetscInt my = info.my;
498 PetscInt mz = info.mz;
499
500 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
501 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
502 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
503 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
504 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
505
506 // Calculate Csi
507 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
508 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
509 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
510
511 PetscReal dx_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node].x - nodal_coords_arr[k_node][j_node-1][i_node].x - nodal_coords_arr[k_node-1][j_node-1][i_node].x);
512 PetscReal dy_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node].y - nodal_coords_arr[k_node][j_node-1][i_node].y - nodal_coords_arr[k_node-1][j_node-1][i_node].y);
513 PetscReal dz_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node].z - nodal_coords_arr[k_node][j_node-1][i_node].z - nodal_coords_arr[k_node-1][j_node-1][i_node].z);
514 PetscReal dx_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node][j_node][i_node].x - nodal_coords_arr[k_node-1][j_node-1][i_node].x - nodal_coords_arr[k_node-1][j_node][i_node].x);
515 PetscReal dy_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node][j_node][i_node].y - nodal_coords_arr[k_node-1][j_node-1][i_node].y - nodal_coords_arr[k_node-1][j_node][i_node].y);
516 PetscReal dz_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node][j_node][i_node].z - nodal_coords_arr[k_node-1][j_node-1][i_node].z - nodal_coords_arr[k_node-1][j_node][i_node].z);
517
518 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
519 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
520 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
521 }
522 }
523 }
524
525 // Calculate Eta
526 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
527 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
528 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
529
530 PetscReal dx_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node].x - nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node-1][j_node][i_node-1].x);
531 PetscReal dy_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node].y - nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node-1][j_node][i_node-1].y);
532 PetscReal dz_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node].z - nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node-1][j_node][i_node-1].z);
533 PetscReal dx_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node-1][j_node][i_node].x - nodal_coords_arr[k_node-1][j_node][i_node-1].x);
534 PetscReal dy_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node-1][j_node][i_node].y - nodal_coords_arr[k_node-1][j_node][i_node-1].y);
535 PetscReal dz_dzeta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node-1][j_node][i_node].z - nodal_coords_arr[k_node-1][j_node][i_node-1].z);
536
537 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
538 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
539 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
540 }
541 }
542 }
543
544 // Calculate Zet
545 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
546 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
547 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
548
549 PetscReal dx_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node].x - nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node][j_node-1][i_node-1].x);
550 PetscReal dy_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node].y - nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node][j_node-1][i_node-1].y);
551 PetscReal dz_dxi = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node].z - nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node][j_node-1][i_node-1].z);
552 PetscReal dx_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x - nodal_coords_arr[k_node][j_node-1][i_node].x - nodal_coords_arr[k_node][j_node-1][i_node-1].x);
553 PetscReal dy_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y - nodal_coords_arr[k_node][j_node-1][i_node].y - nodal_coords_arr[k_node][j_node-1][i_node-1].y);
554 PetscReal dz_deta = 0.5 * (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z - nodal_coords_arr[k_node][j_node-1][i_node].z - nodal_coords_arr[k_node][j_node-1][i_node-1].z);
555
556 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
557 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
558 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
559 }
560 }
561 }
562
563 // --- 4. Boundary Extrapolation ---
564 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
565 PetscInt i_bnd, j_bnd, k_bnd;
566
567 if (xs == 0) { // If this rank owns the global i=0 boundary
568 i_bnd = 0;
569 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
570 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
571 if (i_bnd + 1 < mx) {
572 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
573 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
574 }
575 }
576 }
577 }
578 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
579 i_bnd = mx - 1;
580 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
581 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
582 if (i_bnd - 1 >= 0) {
583 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
584 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
585 }
586 }
587 }
588 }
589 if (ys == 0) {
590 j_bnd = 0;
591 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
592 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
593 if (j_bnd + 1 < my) {
594 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
595 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
596 }
597 }
598 }
599 }
600 if (ye == my) {
601 j_bnd = my - 1;
602 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
603 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
604 if (j_bnd - 1 >= 0) {
605 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
606 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
607 }
608 }
609 }
610 }
611 if (zs == 0) {
612 k_bnd = 0;
613 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
614 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
615 if (k_bnd + 1 < mz) {
616 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
617 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
618 }
619 }
620 }
621 }
622 if (ze == mz) {
623 k_bnd = mz - 1;
624 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
625 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
626 if (k_bnd - 1 >= 0) {
627 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
628 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
629 }
630 }
631 }
632 }
633
634 if (info.xs==0 && info.ys==0 && info.zs==0) {
635 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
636 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
637 }
638
639 // --- 5. Restore all arrays ---
640 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
641 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
642 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
643 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
644
645 // --- 6. Assemble Global Vectors ---
646 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
647 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
648 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
649 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
650
651 // --- 7. Update Local Ghosted Versions ---
652 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
653 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
654 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
655 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
656
657 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
658
659 PROFILE_FUNCTION_END;
660
661 PetscFunctionReturn(0);
662}
663
664
665#undef __FUNCT__
666#define __FUNCT__ "ComputeCellCenteredJacobianInverse"
667/**
668 * @brief Calculates the cell-centered inverse Jacobian determinant (1/J), including
669 * boundary extrapolation, stores it in `user->Aj`, assembles `user->Aj`, and
670 * updates `user->lAj`.
671 *
672 * @param[in,out] user Pointer to the UserCtx structure.
673 *
674 * @return PetscErrorCode 0 on success.
675 */
676PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
677{
678 PetscErrorCode ierr;
679 DMDALocalInfo info;
680 PetscScalar ***aj_arr;
681 Cmpnts ***nodal_coords_arr;
682 Vec localCoords_from_dm;
683
684 PetscFunctionBeginUser;
685 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
686
687 // --- 1. Get Nodal Coordinates and Output Array ---
688 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
689 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
690 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
691 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
692
693 // Define owned node ranges (global indices)
694 PetscInt xs = info.xs, xe = info.xs + info.xm;
695 PetscInt ys = info.ys, ye = info.ys + info.ym;
696 PetscInt zs = info.zs, ze = info.zs + info.zm;
697
698 // Global domain dimensions (total number of nodes)
699 PetscInt mx = info.mx;
700 PetscInt my = info.my;
701 PetscInt mz = info.mz;
702
703 // --- 2. Calculate Aj for INTERIOR Stencils ---
704
705 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
706 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
707 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
708
709 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
710 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
711 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
712
713 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
714 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
715 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
716
717 PetscReal dx_dxi = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node-1][i_node].x) - (nodal_coords_arr[k_node][j_node][i_node-1].x + nodal_coords_arr[k_node][j_node-1][i_node-1].x + nodal_coords_arr[k_node-1][j_node][i_node-1].x + nodal_coords_arr[k_node-1][j_node-1][i_node-1].x) );
718
719 PetscReal dy_dxi = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node-1][i_node].y) - (nodal_coords_arr[k_node][j_node][i_node-1].y + nodal_coords_arr[k_node][j_node-1][i_node-1].y + nodal_coords_arr[k_node-1][j_node][i_node-1].y + nodal_coords_arr[k_node-1][j_node-1][i_node-1].y) );
720
721 PetscReal dz_dxi = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node-1][i_node].z) - (nodal_coords_arr[k_node][j_node][i_node-1].z + nodal_coords_arr[k_node][j_node-1][i_node-1].z + nodal_coords_arr[k_node-1][j_node][i_node-1].z + nodal_coords_arr[k_node-1][j_node-1][i_node-1].z) );
722
723 PetscReal dx_deta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x + nodal_coords_arr[k_node-1][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node-1].x) - (nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node-1].x + nodal_coords_arr[k_node-1][j_node-1][i_node].x + nodal_coords_arr[k_node-1][j_node-1][i_node-1].x) );
724
725 PetscReal dy_deta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y + nodal_coords_arr[k_node-1][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node-1].y) - (nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node-1].y + nodal_coords_arr[k_node-1][j_node-1][i_node].y + nodal_coords_arr[k_node-1][j_node-1][i_node-1].y) );
726
727 PetscReal dz_deta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z + nodal_coords_arr[k_node-1][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node-1].z) - (nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node-1].z + nodal_coords_arr[k_node-1][j_node-1][i_node].z + nodal_coords_arr[k_node-1][j_node-1][i_node-1].z) );
728
729 PetscReal dx_dzeta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].x + nodal_coords_arr[k_node][j_node-1][i_node].x + nodal_coords_arr[k_node][j_node][i_node-1].x + nodal_coords_arr[k_node][j_node-1][i_node-1].x) - (nodal_coords_arr[k_node-1][j_node][i_node].x + nodal_coords_arr[k_node-1][j_node-1][i_node].x + nodal_coords_arr[k_node-1][j_node][i_node-1].x + nodal_coords_arr[k_node-1][j_node-1][i_node-1].x) );
730
731 PetscReal dy_dzeta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].y + nodal_coords_arr[k_node][j_node-1][i_node].y + nodal_coords_arr[k_node][j_node][i_node-1].y + nodal_coords_arr[k_node][j_node-1][i_node-1].y) - (nodal_coords_arr[k_node-1][j_node][i_node].y + nodal_coords_arr[k_node-1][j_node-1][i_node].y + nodal_coords_arr[k_node-1][j_node][i_node-1].y + nodal_coords_arr[k_node-1][j_node-1][i_node-1].y) );
732
733 PetscReal dz_dzeta = 0.25 * ( (nodal_coords_arr[k_node][j_node][i_node].z + nodal_coords_arr[k_node][j_node-1][i_node].z + nodal_coords_arr[k_node][j_node][i_node-1].z + nodal_coords_arr[k_node][j_node-1][i_node-1].z) - (nodal_coords_arr[k_node-1][j_node][i_node].z + nodal_coords_arr[k_node-1][j_node-1][i_node].z + nodal_coords_arr[k_node-1][j_node][i_node-1].z + nodal_coords_arr[k_node-1][j_node-1][i_node-1].z) );
734
735 PetscReal jacobian_det = dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) - dy_dxi * (dx_deta * dz_dzeta - dz_deta * dx_dzeta) + dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta);
736 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
737 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
738 }
739 }
740 }
741
742 // --- 4. Boundary Extrapolation for Aj ---
743 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
744 PetscInt i_bnd, j_bnd, k_bnd;
745
746 if (xs == 0) {
747 i_bnd = 0;
748 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
749 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
750 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
751 }
752 }
753 }
754 if (xe == mx) {
755 i_bnd = mx - 1;
756 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
757 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
758 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
759 }
760 }
761 }
762 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
763 if (ys == 0) {
764 j_bnd = 0;
765 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
766 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
767 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
768 }
769 }
770 }
771 if (ye == my) {
772 j_bnd = my - 1;
773 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
774 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
775 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
776 }
777 }
778 }
779 if (zs == 0) {
780 k_bnd = 0;
781 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
782 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
783 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
784 }
785 }
786 }
787 if (ze == mz) {
788 k_bnd = mz - 1;
789 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
790 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
791 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
792 }
793 }
794 }
795
796 // --- 5. Restore arrays ---
797 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
798 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
799
800 // --- 6. Assemble Global Vector ---
801 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
802 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
803 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
804
805 // --- 7. Update Local Ghosted Version ---
806 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
807 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
808
809 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
810 PetscFunctionReturn(0);
811}
812
813
814#undef __FUNCT__
815#define __FUNCT__ "ComputeCellCentersAndSpacing"
816/**
817 * @brief Computes the physical location of cell centers and the spacing between them.
818 *
819 * This function calculates two key geometric properties from the nodal coordinates:
820 * 1. `Cent`: A vector field storing the (x,y,z) coordinates of the center of each grid cell.
821 * 2. `GridSpace`: A vector field storing the physical distance between adjacent
822 * cell centers in the i, j, and k computational directions.
823 *
824 * It is a direct adaptation of the corresponding logic from the legacy `FormMetrics`.
825 *
826 * @param user The UserCtx for a specific grid level. The function populates `user->Cent` and `user->GridSpace`.
827 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
828 */
829PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
830{
831 PetscErrorCode ierr;
832 DMDALocalInfo info;
833 Vec lCoords;
834 const Cmpnts ***coor;
835 Cmpnts ***cent, ***gs;
836 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
837
838 PetscFunctionBeginUser;
839
840 PROFILE_FUNCTION_BEGIN;
841
842 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing cell centers and spacing for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
843
844 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
845 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
846 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
847
848 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
849 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
850
851 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
852 for (PetscInt k=info.zs+1; k<info.zs+info.zm; k++) {
853 for (PetscInt j=info.ys+1; j<info.ys+info.ym; j++) {
854 for (PetscInt i=info.xs+1; i<info.xs+info.xm; i++) {
855 // Calculate cell center as the average of its 8 corner nodes
856 cent[k][j][i].x = 0.125 * (coor[k][j][i].x + coor[k][j-1][i].x + coor[k-1][j][i].x + coor[k-1][j-1][i].x + coor[k][j][i-1].x + coor[k][j-1][i-1].x + coor[k-1][j][i-1].x + coor[k-1][j-1][i-1].x);
857 cent[k][j][i].y = 0.125 * (coor[k][j][i].y + coor[k][j-1][i].y + coor[k-1][j][i].y + coor[k-1][j-1][i].y + coor[k][j][i-1].y + coor[k][j-1][i-1].y + coor[k-1][j][i-1].y + coor[k-1][j-1][i-1].y);
858 cent[k][j][i].z = 0.125 * (coor[k][j][i].z + coor[k][j-1][i].z + coor[k-1][j][i].z + coor[k-1][j-1][i].z + coor[k][j][i-1].z + coor[k][j-1][i-1].z + coor[k-1][j][i-1].z + coor[k-1][j-1][i-1].z);
859
860 // Calculate Grid Spacing in i-direction (distance between i-face centers)
861 xcp = 0.25 * (coor[k][j][i].x + coor[k][j-1][i].x + coor[k-1][j-1][i].x + coor[k-1][j][i].x);
862 ycp = 0.25 * (coor[k][j][i].y + coor[k][j-1][i].y + coor[k-1][j-1][i].y + coor[k-1][j][i].y);
863 zcp = 0.25 * (coor[k][j][i].z + coor[k][j-1][i].z + coor[k-1][j-1][i].z + coor[k-1][j][i].z);
864 xcm = 0.25 * (coor[k][j][i-1].x + coor[k][j-1][i-1].x + coor[k-1][j-1][i-1].x + coor[k-1][j][i-1].x);
865 ycm = 0.25 * (coor[k][j][i-1].y + coor[k][j-1][i-1].y + coor[k-1][j-1][i-1].y + coor[k-1][j][i-1].y);
866 zcm = 0.25 * (coor[k][j][i-1].z + coor[k][j-1][i-1].z + coor[k-1][j-1][i-1].z + coor[k-1][j][i-1].z);
867 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
868
869 // Calculate Grid Spacing in j-direction (distance between j-face centers)
870 xcp = 0.25 * (coor[k][j][i].x + coor[k][j][i-1].x + coor[k-1][j][i].x + coor[k-1][j][i-1].x);
871 ycp = 0.25 * (coor[k][j][i].y + coor[k][j][i-1].y + coor[k-1][j][i].y + coor[k-1][j][i-1].y);
872 zcp = 0.25 * (coor[k][j][i].z + coor[k][j][i-1].z + coor[k-1][j][i].z + coor[k-1][j][i-1].z);
873 xcm = 0.25 * (coor[k][j-1][i].x + coor[k][j-1][i-1].x + coor[k-1][j-1][i].x + coor[k-1][j-1][i-1].x);
874 ycm = 0.25 * (coor[k][j-1][i].y + coor[k][j-1][i-1].y + coor[k-1][j-1][i].y + coor[k-1][j-1][i-1].y);
875 zcm = 0.25 * (coor[k][j-1][i].z + coor[k][j-1][i-1].z + coor[k-1][j-1][i].z + coor[k-1][j-1][i-1].z);
876 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
877
878 // Calculate Grid Spacing in k-direction (distance between k-face centers)
879 xcp = 0.25 * (coor[k][j][i].x + coor[k][j][i-1].x + coor[k][j-1][i].x + coor[k][j-1][i-1].x);
880 ycp = 0.25 * (coor[k][j][i].y + coor[k][j][i-1].y + coor[k][j-1][i].y + coor[k][j-1][i-1].y);
881 zcp = 0.25 * (coor[k][j][i].z + coor[k][j][i-1].z + coor[k][j-1][i].z + coor[k][j-1][i-1].z);
882 xcm = 0.25 * (coor[k-1][j][i].x + coor[k-1][j][i-1].x + coor[k-1][j-1][i].x + coor[k-1][j-1][i-1].x);
883 ycm = 0.25 * (coor[k-1][j][i].y + coor[k-1][j][i-1].y + coor[k-1][j-1][i].y + coor[k-1][j-1][i-1].y);
884 zcm = 0.25 * (coor[k-1][j][i].z + coor[k-1][j-1][i-1].z + coor[k-1][j-1][i].z + coor[k-1][j-1][i-1].z);
885 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
886 }
887 }
888 }
889
890 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
891 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
892 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
893
894 // Assemble and update ghost regions for the new data
895 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
896 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
897 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
898 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
899
900 PROFILE_FUNCTION_END;
901
902 PetscFunctionReturn(0);
903}
904
905#undef __FUNCT__
906#define __FUNCT__ "ComputeIFaceMetrics"
907/**
908 * @brief Computes metrics centered on constant-i faces (i-faces).
909 *
910 * This function calculates the metric terms (`ICsi`, `IEta`, `IZet`) and the
911 * inverse Jacobian (`IAj`) located at the geometric center of each constant-i
912 * face. This is a critical step for staggered-grid finite difference schemes.
913 *
914 * The process is a direct and faithful refactoring of the corresponding logic
915 * from the legacy `FormMetrics` function:
916 * 1. It first calculates the physical (x,y,z) coordinates of the center of
917 * each i-face and stores them in the `user->Centx` vector.
918 * 2. It then uses a boundary-aware, second-order finite difference stencil on
919 * the `Centx` field to compute the derivatives (e.g., d(x)/d(csi)).
920 * - Central differences are used in the grid interior.
921 * - One-sided differences are used at the physical domain boundaries.
922 * 3. Finally, these derivatives are used to compute the final metric terms and
923 * the inverse Jacobian, which are stored in their respective `Vec` objects.
924 *
925 * @param user The UserCtx for a specific grid level. This function populates
926 * the `user->ICsi`, `user->IEta`, `user->IZet`, and `user->IAj` vectors.
927 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
928 */
929PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
930{
931 PetscErrorCode ierr;
932 DMDALocalInfo info;
933 Vec lCoords;
934 const Cmpnts ***coor;
935 Cmpnts ***centx; //***gs;
936 const Cmpnts ***centx_const;
937 Cmpnts ***icsi, ***ieta, ***izet;
938 PetscScalar ***iaj;
939 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
940
941 PetscFunctionBeginUser;
942
943 PROFILE_FUNCTION_BEGIN;
944
945 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
946
947 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
948 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
949 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
950 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
951 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
952 PetscInt gys = info.gys, gye = info.gys + info.gym;
953 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
954
955 PetscInt lxs = xs; PetscInt lxe = xe;
956 PetscInt lys = ys; PetscInt lye = ye;
957 PetscInt lzs = zs; PetscInt lze = ze;
958
959 if (xs==0) lxs = xs+1;
960 if (ys==0) lys = ys+1;
961 if (zs==0) lzs = zs+1;
962
963 if (xe==mx) lxe=xe-1;
964 if (ye==my) lye=ye-1;
965 if (ze==mz) lze=ze-1;
966
967 // --- Part 1: Calculate the location of i-face centers (Centx) ---
968 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
969 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
970 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
971 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
972
973 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calculating i-face centers (Centx) with i[%d,%d], j[%d,%d], k[%d,%d] ...\n", user->simCtx->rank,gxs,gxe,gys,gye,gzs,gze);
974
975 // Loop over the ghosted region to calculate all local face centers
976 for (PetscInt k = gzs + 1; k < gze; k++) {
977 for (PetscInt j = gys + 1; j < gye; j++) {
978 for (PetscInt i = gxs; i < gxe; i++) {
979 //----- DEBUG ------
980 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calculating i-face center at (k=%d, j=%d, i=%d)\n", user->simCtx->rank, k, j, i);
981 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Using corner nodes: (%f,%f,%f), (%f,%f,%f), (%f,%f,%f), (%f,%f,%f)\n", user->simCtx->rank,
982 // coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
983 // coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
984 // coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
985 // coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
986
987 centx[k][j][i].x = 0.25 * (coor[k][j][i].x + coor[k-1][j][i].x + coor[k][j-1][i].x + coor[k-1][j-1][i].x);
988 centx[k][j][i].y = 0.25 * (coor[k][j][i].y + coor[k-1][j][i].y + coor[k][j-1][i].y + coor[k-1][j-1][i].y);
989 centx[k][j][i].z = 0.25 * (coor[k][j][i].z + coor[k-1][j][i].z + coor[k][j-1][i].z + coor[k-1][j-1][i].z);
990
991 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calculated i-face center: (%f,%f,%f)\n", user->simCtx->rank, centx[k][j][i].x, centx[k][j][i].y, centx[k][j][i].z);
992 }
993 }
994 }
995
996 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
997 /*
998 if(xs==0){
999 for(PetscInt k=gzs+1;k < gze; k++){
1000 for(PetscInt j=gys+1;j < gye; j++){
1001 PetscInt i=0;
1002 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
1003 centx[k][j][i-1].y=centx[k][j][i].y;
1004 centx[k][j][i-1].z=centx[k][j][i].z;
1005 }
1006 }
1007 }
1008 if (xe==mx){
1009 for(PetscInt k=gzs+1; k<gze; k++) {
1010 for (PetscInt j=gys+1; j<gye;j++) {
1011 PetscInt i=mx-1;
1012 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
1013 centx[k][j][i].y=centx[k][j][i-1].y;
1014 centx[k][j][i].z=centx[k][j][i-1].z;
1015 }
1016 }
1017 }
1018 */
1019
1020 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1021 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1022
1023 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1024
1025 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1026
1027 // --- Part 2: Calculate metrics using face-centered coordinates ---
1028 ierr = DMDAVecGetArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1029 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1030 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1031 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1032 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1033
1034 // Loop over the OWNED region where we will store the final metrics
1035 for (PetscInt k=lzs; k<lze; k++) {
1036 for (PetscInt j=lys; j<lye; j++) {
1037 for (PetscInt i=xs; i<lxe; i++) {
1038
1039 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1040 if (i == 0) { // Forward difference at the domain's min-i boundary
1041 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
1042 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
1043 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
1044 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1045 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
1046 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
1047 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
1048 } else { // Central difference in the interior
1049 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
1050 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
1051 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
1052 }
1053
1054 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1055 if (j == 1) { // Forward difference
1056 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
1057 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
1058 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
1059 } else if (j == my - 2) { // Backward difference
1060 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
1061 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
1062 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
1063 } else { // Central difference
1064 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
1065 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
1066 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
1067 }
1068
1069 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1070 if (k == 1) { // Forward difference
1071 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
1072 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
1073 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
1074 } else if (k == mz - 2) { // Backward difference
1075 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
1076 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
1077 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
1078 } else { // Central difference
1079 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
1080 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
1081 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
1082 }
1083
1084 // --- Metric calculations (identical to legacy FormMetrics) ---
1085 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1086 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1087 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1088
1089 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1090 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1091 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1092
1093 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
1094 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1095 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
1096
1097 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
1098 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
1099 iaj[k][j][i] = 1.0 / iaj[k][j][i];
1100 }
1101 }
1102 }
1103 }
1104
1105 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
1106 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1107 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1108 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1109 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1110
1111 // --- Part 3: Assemble global vectors and update local ghosts ---
1112 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
1113 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
1114 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
1115 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
1116
1117 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
1118 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
1119 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
1120 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
1121
1122 PROFILE_FUNCTION_END;
1123
1124 PetscFunctionReturn(0);
1125}
1126
1127#undef __FUNCT__
1128#define __FUNCT__ "ComputeJFaceMetrics"
1129/**
1130 * @brief Computes metrics centered on constant-j faces (j-faces).
1131 *
1132 * This function calculates the metric terms (`JCsi`, `JEta`, `JZet`) and the
1133 * inverse Jacobian (`JAj`) located at the geometric center of each constant-j
1134 * face. This is a critical step for staggered-grid finite difference schemes.
1135 *
1136 * The process is a direct and faithful refactoring of the corresponding logic
1137 * from the legacy `FormMetrics` function:
1138 * 1. It first calculates the physical (x,y,z) coordinates of the center of
1139 * each i-face and stores them in the `user->Centy` vector.
1140 * 2. It then uses a boundary-aware, second-order finite difference stencil on
1141 * the `Centy` field to compute the derivatives (e.g., d(x)/d(csi)).
1142 * - Central differences are used in the grid interior.
1143 * - One-sided differences are used at the physical domain boundaries.
1144 * 3. Finally, these derivatives are used to compute the final metric terms and
1145 * the inverse Jacobian, which are stored in their respective `Vec` objects.
1146 *
1147 * @param user The UserCtx for a specific grid level. This function populates
1148 * the `user->JCsi`, `user->JEta`, `user->JZet`, and `user->JAj` vectors.
1149 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1150 */
1151PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
1152{
1153 PetscErrorCode ierr;
1154 DMDALocalInfo info;
1155 Vec lCoords;
1156 const Cmpnts ***coor;
1157 Cmpnts ***centy; //***gs;
1158 const Cmpnts ***centy_const;
1159 Cmpnts ***jcsi, ***jeta, ***jzet;
1160 PetscScalar ***jaj;
1161 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1162
1163 PetscFunctionBeginUser;
1164
1165 PROFILE_FUNCTION_BEGIN;
1166
1167 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1168
1169 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1170 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1171 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1172 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1173 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1174 PetscInt gys = info.gys, gye = info.gys + info.gym;
1175 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1176
1177 PetscInt lxs = xs; PetscInt lxe = xe;
1178 PetscInt lys = ys; PetscInt lye = ye;
1179 PetscInt lzs = zs; PetscInt lze = ze;
1180
1181 if (xs==0) lxs = xs+1;
1182 if (ys==0) lys = ys+1;
1183 if (zs==0) lzs = zs+1;
1184
1185 if (xe==mx) lxe=xe-1;
1186 if (ye==my) lye=ye-1;
1187 if (ze==mz) lze=ze-1;
1188
1189 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1190 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1191 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1192 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1193 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1194
1195 // Loop over the ghosted region to calculate all local face centers
1196 for (PetscInt k = gzs + 1; k < gze; k++) {
1197 for (PetscInt j = gys; j < gye; j++) {
1198 for (PetscInt i = gxs + 1; i < gxe; i++) {
1199 centy[k][j][i].x = 0.25 * (coor[k][j][i].x + coor[k-1][j][i].x + coor[k][j][i-1].x + coor[k-1][j][i-1].x);
1200 centy[k][j][i].y = 0.25 * (coor[k][j][i].y + coor[k-1][j][i].y + coor[k][j][i-1].y + coor[k-1][j][i-1].y);
1201 centy[k][j][i].z = 0.25 * (coor[k][j][i].z + coor[k-1][j][i].z + coor[k][j][i-1].z + coor[k-1][j][i-1].z);
1202 }
1203 }
1204 }
1205
1206 /*
1207 if(ys==0){
1208 for(PetscInt k=gzs+1;k < gze; k++){
1209 for(PetscInt i=gxs+1;j < gxe; i++){
1210 PetscInt j=0;
1211 centy[k][j-1][i].x=centy[k][j][i].x;
1212 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1213 centy[k][j-1][i].z=centy[k][j][i].z;
1214 }
1215 }
1216 }
1217 if (ye==my){
1218 for(PetscInt k=gzs+1; k<gze; k++) {
1219 for (PetscInt i=gxs+1; j<gxe;i++) {
1220 PetscInt j=my-1;
1221 centy[k][j][i].x=centy[k][j-1][i].x
1222 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1223 centy[k][j][i].z=centy[k][j-1][i].z;
1224 }
1225 }
1226 }
1227 */
1228
1229 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1230 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1231 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1232
1233 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1234
1235 // --- Part 2: Calculate metrics using face-centered coordinates ---
1236 ierr = DMDAVecGetArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1237 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1238 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1239 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1240 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1241
1242 // Loop over the OWNED region where we will store the final metrics
1243 for (PetscInt k=lzs; k<lze; k++) {
1244 for (PetscInt j=ys; j<lye; j++) {
1245 for (PetscInt i=lxs; i<lxe; i++) {
1246
1247 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1248 if (i == 1) { // Forward difference at the domain's min-i boundary
1249 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1250 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1251 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1252 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1253 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1254 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1255 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1256 } else { // Central difference in the interior
1257 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1258 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1259 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1260 }
1261
1262 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1263 if (j == 0) { // Forward difference
1264 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1265 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1266 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1267 } else if (j == my - 2) { // Backward difference
1268 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1269 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1270 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1271 } else { // Central difference
1272 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1273 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1274 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1275 }
1276
1277 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1278 if (k == 1) { // Forward difference
1279 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1280 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1281 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1282 } else if (k == mz - 2) { // Backward difference
1283 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1284 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1285 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1286 } else { // Central difference
1287 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1288 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1289 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1290 }
1291
1292 // --- Metric calculations (identical to legacy FormMetrics) ---
1293 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1294 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1295 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1296
1297 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1298 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1299 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1300
1301 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1302 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1303 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1304
1305 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1306 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1307 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1308 }
1309 }
1310 }
1311 }
1312
1313 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1314 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1315 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1316 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1317 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1318
1319 // --- Part 3: Assemble global vectors and update local ghosts ---
1320 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1321 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1322 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1323 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1324
1325 ierr = UpdateLocalGhosts(user, "JCsi"); CHKERRQ(ierr);
1326 ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
1327 ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
1328 ierr = UpdateLocalGhosts(user, "JAj"); CHKERRQ(ierr);
1329
1330 PROFILE_FUNCTION_END;
1331
1332 PetscFunctionReturn(0);
1333}
1334
1335#undef __FUNCT__
1336#define __FUNCT__ "ComputeJFaceMetrics"
1337/**
1338 * @brief Computes metrics centered on constant-k faces (k-faces).
1339 *
1340 * This function calculates the metric terms (`KCsi`, `KEta`, `KZet`) and the
1341 * inverse Jacobian (`KAj`) located at the geometric center of each constant-k
1342 * face. This is a critical step for staggered-grid finite difference schemes.
1343 *
1344 * The process is a direct and faithful refactoring of the corresponding logic
1345 * from the legacy `FormMetrics` function:
1346 * 1. It first calculates the physical (x,y,z) coordinates of the center of
1347 * each i-face and stores them in the `user->Centz` vector.
1348 * 2. It then uses a boundary-aware, second-order finite difference stencil on
1349 * the `Centz` field to compute the derivatives (e.g., d(x)/d(csi)).
1350 * - Central differences are used in the grid interior.
1351 * - One-sided differences are used at the physical domain boundaries.
1352 * 3. Finally, these derivatives are used to compute the final metric terms and
1353 * the inverse Jacobian, which are stored in their respective `Vec` objects.
1354 *
1355 * @param user The UserCtx for a specific grid level. This function populates
1356 * the `user->KCsi`, `user->KEta`, `user->KZet`, and `user->KAj` vectors.
1357 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1358 */
1359PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
1360{
1361 PetscErrorCode ierr;
1362 DMDALocalInfo info;
1363 Vec lCoords;
1364 const Cmpnts ***coor;
1365 Cmpnts ***centz; //***gs;
1366 const Cmpnts ***centz_const;
1367 Cmpnts ***kcsi, ***keta, ***kzet;
1368 PetscScalar ***kaj;
1369 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1370
1371 PetscFunctionBeginUser;
1372
1373 PROFILE_FUNCTION_BEGIN;
1374
1375 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1376
1377 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1378 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1379 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1380 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1381 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1382 PetscInt gys = info.gys, gye = info.gys + info.gym;
1383 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1384
1385 PetscInt lxs = xs; PetscInt lxe = xe;
1386 PetscInt lys = ys; PetscInt lye = ye;
1387 PetscInt lzs = zs; PetscInt lze = ze;
1388
1389 if (xs==0) lxs = xs+1;
1390 if (ys==0) lys = ys+1;
1391 if (zs==0) lzs = zs+1;
1392
1393 if (xe==mx) lxe=xe-1;
1394 if (ye==my) lye=ye-1;
1395 if (ze==mz) lze=ze-1;
1396
1397 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1398 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1399 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1400 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1401 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1402
1403 // Loop over the ghosted region to calculate all local face centers
1404 for (PetscInt k = gzs; k < gze; k++) {
1405 for (PetscInt j = gys; j < gye; j++) {
1406 for (PetscInt i = gxs + 1; i < gxe; i++) {
1407 centz[k][j][i].x = 0.25 * (coor[k][j][i].x + coor[k][j-1][i].x + coor[k][j][i-1].x + coor[k][j-1][i-1].x);
1408 centz[k][j][i].y = 0.25 * (coor[k][j][i].y + coor[k][j-1][i].y + coor[k][j][i-1].y + coor[k][j-1][i-1].y);
1409 centz[k][j][i].z = 0.25 * (coor[k][j][i].z + coor[k][j-1][i].z + coor[k][j][i-1].z + coor[k][j-1][i-1].z);
1410 }
1411 }
1412 }
1413
1414 /*
1415 if(zs==0){
1416 for(PetscInt j=gys+1;j < gye; j++){
1417 for(PetscInt i=gxs+1;j < gxe; i++){
1418 PetscInt k=0;
1419 centz[k-1][j][i].x=centz[k][j][i].x;
1420 centz[k-1][j][i].y=centz[k][j][i].y;
1421 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1422 }
1423 }
1424 }
1425 if (ze==mz){
1426 for(PetscInt j=gys+1; j<gye; j++) {
1427 for (PetscInt i=gxs+1; j<gxe;i++) {
1428 PetscInt k=mz-1;
1429 centy[k][j][i].x=centy[k-1][j][i].x
1430 centy[k][j][i].y=centy[k-1][j][i].y;
1431 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1432 }
1433 }
1434 }
1435 */
1436
1437 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1438 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1439 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1440
1441 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1442
1443 // --- Part 2: Calculate metrics using face-centered coordinates ---
1444 ierr = DMDAVecGetArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1445 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1446 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1447 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1448 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1449
1450 // Loop over the OWNED region where we will store the final metrics
1451 for (PetscInt k=zs; k<lze; k++) {
1452 for (PetscInt j=lys; j<lye; j++) {
1453 for (PetscInt i=lxs; i<lxe; i++) {
1454
1455 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1456 if (i == 1) { // Forward difference at the domain's min-i boundary
1457 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
1458 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
1459 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
1460 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1461 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
1462 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
1463 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
1464 } else { // Central difference in the interior
1465 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
1466 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
1467 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
1468 }
1469
1470 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1471 if (j == 1) { // Forward difference
1472 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
1473 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
1474 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
1475 } else if (j == my - 2) { // Backward difference
1476 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
1477 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
1478 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
1479 } else { // Central difference
1480 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
1481 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
1482 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
1483 }
1484
1485 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1486 if (k == 0) { // Forward difference
1487 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
1488 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
1489 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
1490 } else if (k == mz - 2) { // Backward difference
1491 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
1492 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
1493 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
1494 } else { // Central difference
1495 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
1496 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
1497 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
1498 }
1499
1500 // --- Metric calculations (identical to legacy FormMetrics) ---
1501 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1502 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1503 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1504
1505 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1506 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1507 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1508
1509 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1510 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1511 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1512
1513 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
1514 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
1515 kaj[k][j][i] = 1.0 / kaj[k][j][i];
1516 }
1517 }
1518 }
1519 }
1520
1521 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1522 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1523 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1524 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1525 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1526
1527 // --- Part 3: Assemble global vectors and update local ghosts ---
1528 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
1529 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
1530 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
1531 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
1532
1533 ierr = UpdateLocalGhosts(user, "KCsi"); CHKERRQ(ierr);
1534 ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
1535 ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
1536 ierr = UpdateLocalGhosts(user, "KAj"); CHKERRQ(ierr);
1537
1538 PROFILE_FUNCTION_END;
1539
1540 PetscFunctionReturn(0);
1541}
1542
1543static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
1544{
1545 PetscInt nidx;
1546 DMDALocalInfo info = user->info;
1547
1548 PetscInt mx = info.mx, my = info.my;
1549
1550 AO ao;
1551 DMDAGetAO(user->da, &ao);
1552 nidx=i+j*mx+k*mx*my;
1553
1554 AOApplicationToPetsc(ao,1,&nidx);
1555
1556 return (nidx);
1557}
1558
1559#undef __FUNCT__
1560#define __FUNCT__ "ComputeMetricsDivergence"
1561/**
1562 * @brief Computes the divergence of the grid metrics and identifies the maximum value.
1563 *
1564 * This function serves as a diagnostic tool to assess the quality of the grid
1565 * metrics. It calculates the divergence of the face metrics (Csi, Eta, Zet)
1566 * and scales it by the inverse of the cell Jacobian. The maximum divergence
1567 * value is then located, and its grid coordinates are printed to the console,
1568 * helping to pinpoint areas of potential grid quality issues.
1569 *
1570 * @param user The UserCtx, containing all necessary grid data.
1571 * @return PetscErrorCode
1572 */
1573PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
1574{
1575 DM da = user->da, fda = user->fda;
1576 DMDALocalInfo info = user->info;
1577 PetscInt xs = info.xs, xe = info.xs + info.xm;
1578 PetscInt ys = info.ys, ye = info.ys + info.ym;
1579 PetscInt zs = info.zs, ze = info.zs + info.zm;
1580 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1581 PetscInt lxs, lys, lzs, lxe, lye, lze;
1582 PetscInt i, j, k;
1583 Vec Div;
1584 PetscReal ***div, ***aj;
1585 Cmpnts ***csi, ***eta, ***zet;
1586 PetscReal maxdiv;
1587
1588 PetscFunctionBeginUser;
1589
1590 PROFILE_FUNCTION_BEGIN;
1591
1592 lxs = xs; lxe = xe;
1593 lys = ys; lye = ye;
1594 lzs = zs; lze = ze;
1595
1596 if (xs == 0) lxs = xs + 1;
1597 if (ys == 0) lys = ys + 1;
1598 if (zs == 0) lzs = zs + 1;
1599
1600 if (xe == mx) lxe = xe - 1;
1601 if (ye == my) lye = ye - 1;
1602 if (ze == mz) lze = ze - 1;
1603
1604 DMDAVecGetArray(fda, user->lCsi, &csi);
1605 DMDAVecGetArray(fda, user->lEta, &eta);
1606 DMDAVecGetArray(fda, user->lZet, &zet);
1607 DMDAVecGetArray(da, user->lAj, &aj);
1608
1609 VecDuplicate(user->P, &Div);
1610 VecSet(Div, 0.);
1611 DMDAVecGetArray(da, Div, &div);
1612
1613 for (k = lzs; k < lze; k++) {
1614 for (j = lys; j < lye; j++) {
1615 for (i = lxs; i < lxe; i++) {
1616 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1617 eta[k][j][i].x - eta[k][j-1][i].x +
1618 zet[k][j][i].x - zet[k-1][j][i].x +
1619 csi[k][j][i].y - csi[k][j][i-1].y +
1620 eta[k][j][i].y - eta[k][j-1][i].y +
1621 zet[k][j][i].y - zet[k-1][j][i].y +
1622 csi[k][j][i].z - csi[k][j][i-1].z +
1623 eta[k][j][i].z - eta[k][j-1][i].z +
1624 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1625 div[k][j][i] = fabs(divergence);
1626 }
1627 }
1628 }
1629
1630 DMDAVecRestoreArray(da, Div, &div);
1631
1632 PetscReal MaxFlatIndex;
1633 VecMax(Div, &MaxFlatIndex, &maxdiv);
1634 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %d.\n",maxdiv,MaxFlatIndex);
1635
1636 for (k=zs; k<ze; k++) {
1637 for (j=ys; j<ye; j++) {
1638 for (i=xs; i<xe; i++) {
1639 if (Gidx(i,j,k,user) == MaxFlatIndex) {
1640 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,k,j,i);
1641 }
1642 }
1643 }
1644 }
1645
1646
1647 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1648 DMDAVecRestoreArray(fda, user->lEta, &eta);
1649 DMDAVecRestoreArray(fda, user->lZet, &zet);
1650 DMDAVecRestoreArray(da, user->lAj, &aj);
1651 VecDestroy(&Div);
1652
1653
1654 PROFILE_FUNCTION_END;
1655
1656 PetscFunctionReturn(0);
1657}
1658
1659#undef __FUNCT__
1660#define __FUNCT__ "ComputeMetricNorms"
1661/**
1662 * @brief Computes the max-min values of the grid metrics.
1663 *
1664 * This function serves as a diagnostic tool to assess the quality of the grid
1665 * metrics. It calculates the bounds of the face metrics (Csi, Eta, Zet).
1666 *
1667 * @param user The UserCtx, containing all necessary grid data.
1668 * @return PetscErrorCode
1669 */
1670PetscErrorCode ComputeMetricNorms(UserCtx *user)
1671{
1672
1673 DMDALocalInfo info = user->info;
1674 PetscInt xs = info.xs, xe = info.xs + info.xm;
1675 PetscInt ys = info.ys, ye = info.ys + info.ym;
1676 PetscInt zs = info.zs, ze = info.zs + info.zm;
1677 PetscInt i, j, k;
1678
1679 PetscFunctionBeginUser;
1680
1681 PROFILE_FUNCTION_BEGIN;
1682
1683 PetscReal CsiMax, EtaMax, ZetMax;
1684 PetscReal ICsiMax, IEtaMax, IZetMax;
1685 PetscReal JCsiMax, JEtaMax, JZetMax;
1686 PetscReal KCsiMax, KEtaMax, KZetMax;
1687 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
1688
1689 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
1690 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
1691 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
1692 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
1693 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
1694
1695 // Max Values
1696 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
1697 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
1698 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
1699
1700 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
1701 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
1702 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
1703
1704 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
1705 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
1706 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
1707
1708 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
1709 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
1710 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
1711
1712 VecMax(user->lAj,&AjMaxArg,&AjMax);
1713 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1714 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1715 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1716
1717 VecMax(user->lAj,&AjMaxArg,&AjMax);
1718 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1719 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1720 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1721
1722 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
1723
1724 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
1725 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
1726 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
1727 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
1728 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
1729
1730 for (k=zs; k<ze; k++) {
1731 for (j=ys; j<ye; j++) {
1732 for (i=xs; i<xe; i++) {
1733 if (Gidx(i,j,k,user) == CsiMaxArg) {
1734 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
1735 }
1736 if (Gidx(i,j,k,user) == EtaMaxArg) {
1737 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
1738 }
1739 if (Gidx(i,j,k,user) == ZetMaxArg) {
1740 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
1741 }
1742 if (Gidx(i,j,k,user) == ICsiMaxArg) {
1743 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
1744 }
1745 if (Gidx(i,j,k,user) == IEtaMaxArg) {
1746 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
1747 }
1748 if (Gidx(i,j,k,user) == IZetMaxArg) {
1749 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
1750 }
1751 if (Gidx(i,j,k,user) == JCsiMaxArg) {
1752 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
1753 }
1754 if (Gidx(i,j,k,user) == JEtaMaxArg) {
1755 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
1756 }
1757 if (Gidx(i,j,k,user) == JZetMaxArg) {
1758 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
1759 }
1760 if (Gidx(i,j,k,user) == KCsiMaxArg) {
1761 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
1762 }
1763 if (Gidx(i,j,k,user) == KEtaMaxArg) {
1764 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
1765 }
1766 if (Gidx(i,j,k,user) == KZetMaxArg) {
1767 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
1768 }
1769 if (Gidx(i,j,k,user) == AjMaxArg) {
1770 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
1771 }
1772 if (Gidx(i,j,k,user) == IAjMaxArg) {
1773 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
1774 }
1775 if (Gidx(i,j,k,user) == JAjMaxArg) {
1776 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
1777 }
1778 if (Gidx(i,j,k,user) == KAjMaxArg) {
1779 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
1780 }
1781 }
1782 }
1783 }
1784
1785 /*
1786 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
1787 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
1788 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
1789 */
1790
1791 PROFILE_FUNCTION_END;
1792
1793 PetscFunctionReturn(0);
1794}
1795
1796#undef __FUNCT__
1797#define __FUNCT__ "ComputeAllGridMetrics"
1798/**
1799 * @brief Orchestrates the calculation of all grid metrics.
1800 *
1801 * This function iterates through every UserCtx in the multigrid and multi-block
1802 * hierarchy. For each context, it calls a series of modern, modular helper
1803 * functions to compute the face metrics (Csi, Eta, Zet), the cell-centered
1804 * inverse Jacobian (Aj), and to validate the grid's orientation.
1805 *
1806 * @param simCtx The master SimCtx, containing the configured UserCtx hierarchy.
1807 * @return PetscErrorCode
1808 */
1809PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
1810{
1811 PetscErrorCode ierr;
1812 UserMG *usermg = &simCtx->usermg;
1813 MGCtx *mgctx = usermg->mgctx;
1814 PetscInt nblk = simCtx->block_number;
1815
1816 PetscFunctionBeginUser;
1817
1818 PROFILE_FUNCTION_BEGIN;
1819
1820 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
1821
1822 // Loop through all levels and all blocks
1823 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
1824 for (PetscInt bi = 0; bi < nblk; bi++) {
1825 UserCtx *user = &mgctx[level].user[bi];
1826 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
1827
1828 // Call the modern, modular helper functions for each UserCtx.
1829 // These functions are self-contained and operate on the data within the provided context.
1830 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
1831 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
1832 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
1833 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
1834 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
1835 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
1836 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
1837
1838 // Diagnostics
1839 ierr = ComputeMetricNorms(user);
1840 if (level == usermg->mglevels - 1) {
1841 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
1842 }
1843 }
1844 }
1845
1846 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
1847
1848 PROFILE_FUNCTION_END;
1849
1850 PetscFunctionReturn(0);
1851}
1852/* ------------------------------------------------------------------------- */
1853/* End of Metric.c */
1854
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Orchestrates the calculation of all grid metrics.
Definition Metric.c:1809
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:1670
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:1543
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain...
Definition Metric.c:368
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Computes the physical location of cell centers and the spacing between them.
Definition Metric.c:829
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:1151
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Computes the divergence of the grid metrics and identifies the maximum value.
Definition Metric.c:1573
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation,...
Definition Metric.c:464
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
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Calculates the cell-centered inverse Jacobian determinant (1/J), including boundary extrapolation,...
Definition Metric.c:676
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Computes metrics centered on constant-k faces (k-faces).
Definition Metric.c:1359
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:929
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:313
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:268
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:201
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1091
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:705
UserCtx::JCsi
Vec JCsi
Definition variables.h:708
UserCtx::KAj
Vec KAj
Definition variables.h:709
MGCtx::user
UserCtx * user
Definition variables.h:427
UserCtx::JEta
Vec JEta
Definition variables.h:708
UserCtx::Zet
Vec Zet
Definition variables.h:705
SimCtx::rank
PetscMPIInt rank
Definition variables.h:541
SimCtx::block_number
PetscInt block_number
Definition variables.h:593
UserCtx::lIEta
Vec lIEta
Definition variables.h:707
UserCtx::lIZet
Vec lIZet
Definition variables.h:707
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
UserCtx::IZet
Vec IZet
Definition variables.h:707
UserCtx::Centz
Vec Centz
Definition variables.h:706
UserCtx::IEta
Vec IEta
Definition variables.h:707
UserCtx::lZet
Vec lZet
Definition variables.h:705
SimCtx::usermg
UserMG usermg
Definition variables.h:631
UserCtx::Csi
Vec Csi
Definition variables.h:705
UserCtx::lIAj
Vec lIAj
Definition variables.h:707
UserCtx::_this
PetscInt _this
Definition variables.h:674
UserCtx::lKEta
Vec lKEta
Definition variables.h:709
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:708
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::JZet
Vec JZet
Definition variables.h:708
UserCtx::Centx
Vec Centx
Definition variables.h:706
UserCtx::lKZet
Vec lKZet
Definition variables.h:709
UserCtx::Eta
Vec Eta
Definition variables.h:705
UserCtx::lJEta
Vec lJEta
Definition variables.h:708
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:721
UserCtx::ICsi
Vec ICsi
Definition variables.h:707
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:709
UserMG::mglevels
PetscInt mglevels
Definition variables.h:434
UserCtx::lJZet
Vec lJZet
Definition variables.h:708
UserCtx::IAj
Vec IAj
Definition variables.h:707
UserCtx::JAj
Vec JAj
Definition variables.h:708
UserCtx::KEta
Vec KEta
Definition variables.h:709
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::lICsi
Vec lICsi
Definition variables.h:707
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:674
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:705
UserCtx::KZet
Vec KZet
Definition variables.h:709
UserCtx::Cent
Vec Cent
Definition variables.h:705
UserCtx::KCsi
Vec KCsi
Definition variables.h:709
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:437
UserCtx::Centy
Vec Centy
Definition variables.h:706
UserCtx::lJAj
Vec lJAj
Definition variables.h:708
UserCtx::lKAj
Vec lKAj
Definition variables.h:709
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
MGCtx
Context for Multigrid operations.
Definition variables.h:426
SimCtx
The master context for the entire simulation.
Definition variables.h:538
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:661
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:433