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 Implementation of \ref MetricGetCellVertices().
22 * @details Full API contract (arguments, ownership, side effects) is documented with
23 * the header declaration in `include/Metric.h`.
24 * @see MetricGetCellVertices()
25 */
26PetscErrorCode MetricGetCellVertices(UserCtx *user,
27 const Cmpnts ***X, /* coord array */
28 PetscInt i,PetscInt j,PetscInt k,
29 Cmpnts V[8])
30{
31 (void)user;
32 PetscFunctionBeginUser;
33 for (PetscInt c = 0; c < 8; ++c) {
34 PetscInt ii = i + ((c & 1) ? 1 : 0);
35 PetscInt jj = j + ((c & 2) ? 1 : 0);
36 PetscInt kk = k + ((c & 4) ? 1 : 0);
37 LOG_LOOP_ALLOW(GLOBAL, LOG_VERBOSE,i+j+k,10," ii: %d,jj:%d,kk:%d - Retrieved.\n",ii,jj,kk);
38 V[c] = X[kk][jj][ii];
39 }
40 PetscFunctionReturn(0);
41}
42
43#undef __FUNCT__
44#define __FUNCT__ "TrilinearBlend"
45
46/* ------------------------------------------------------------------------- */
47/**
48 * @brief Blend eight corner values at the supplied trilinear reference coordinates.
49 */
50static inline void TrilinearBlend(const Cmpnts V[8],
51 PetscReal xi,PetscReal eta,PetscReal zta,
52 Cmpnts *Xp)
53{
54 PetscReal x=0,y=0,z=0;
55 for (PetscInt c=0;c<8;++c) {
56 PetscReal N = ((c&1)?xi : 1.0-xi ) *
57 ((c&2)?eta: 1.0-eta) *
58 ((c&4)?zta: 1.0-zta);
59 x += N * V[c].x;
60 y += N * V[c].y;
61 z += N * V[c].z;
62 }
63 Xp->x = x; Xp->y = y; Xp->z = z;
64}
65
66
67#undef __FUNCT__
68#define __FUNCT__ "MetricLogicalToPhysical"
69/* ------------------------------------------------------------------------- */
70/**
71 * @brief Implementation of \ref MetricLogicalToPhysical().
72 * @details Full API contract (arguments, ownership, side effects) is documented with
73 * the header declaration in `include/Metric.h`.
74 * @see MetricLogicalToPhysical()
75 */
76PetscErrorCode MetricLogicalToPhysical(UserCtx *user,
77 const Cmpnts ***X,
78 PetscInt i,PetscInt j,PetscInt k,
79 PetscReal xi,PetscReal eta,PetscReal zta,
80 Cmpnts *Xp)
81{
82 PetscErrorCode ierr;
83 Cmpnts V[8];
84 PetscFunctionBeginUser;
85
86 PROFILE_FUNCTION_BEGIN;
87
88 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
89 TrilinearBlend(V,xi,eta,zta,Xp);
90
91 PROFILE_FUNCTION_END;
92
93 PetscFunctionReturn(0);
94}
95
96#undef __FUNCT__
97#define __FUNCT__ "MetricJacobian"
98/* ------------------------------------------------------------------------- */
99/**
100 * @brief Implementation of \ref MetricJacobian().
101 * @details Full API contract (arguments, ownership, side effects) is documented with
102 * the header declaration in `include/Metric.h`.
103 * @see MetricJacobian()
104 */
105PetscErrorCode MetricJacobian(UserCtx *user,
106 const Cmpnts ***X,
107 PetscInt i,PetscInt j,PetscInt k,
108 PetscReal xi,PetscReal eta,PetscReal zta,
109 PetscReal J[3][3], PetscReal *detJ)
110{
111 PetscErrorCode ierr;
112 Cmpnts V[8];
113 PetscFunctionBeginUser;
114
115 PROFILE_FUNCTION_BEGIN;
116
117 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
118
119 /* derivatives of trilinear shape functions */
120 PetscReal dN_dXi[8], dN_dEta[8], dN_dZta[8];
121 for (PetscInt c=0;c<8;++c) {
122 PetscReal sx = (c & 1) ? 1.0 : -1.0;
123 PetscReal sy = (c & 2) ? 1.0 : -1.0;
124 PetscReal sz = (c & 4) ? 1.0 : -1.0;
125 dN_dXi [c] = 0.125 * sx * ( (c&2?eta:1-eta) ) * ( (c&4?zta:1-zta) );
126 dN_dEta[c] = 0.125 * sy * ( (c&1?xi :1-xi ) ) * ( (c&4?zta:1-zta) );
127 dN_dZta[c] = 0.125 * sz * ( (c&1?xi :1-xi ) ) * ( (c&2?eta:1-eta) );
128 }
129
130 /* assemble Jacobian */
131 PetscReal x_xi=0,y_xi=0,z_xi=0,
132 x_eta=0,y_eta=0,z_eta=0,
133 x_zta=0,y_zta=0,z_zta=0;
134 for (PetscInt c=0;c<8;++c) {
135 x_xi += dN_dXi [c]*V[c].x; y_xi += dN_dXi [c]*V[c].y; z_xi += dN_dXi [c]*V[c].z;
136 x_eta += dN_dEta[c]*V[c].x; y_eta += dN_dEta[c]*V[c].y; z_eta += dN_dEta[c]*V[c].z;
137 x_zta += dN_dZta[c]*V[c].x; y_zta += dN_dZta[c]*V[c].y; z_zta += dN_dZta[c]*V[c].z;
138 }
139
140 J[0][0]=x_xi; J[0][1]=x_eta; J[0][2]=x_zta;
141 J[1][0]=y_xi; J[1][1]=y_eta; J[1][2]=y_zta;
142 J[2][0]=z_xi; J[2][1]=z_eta; J[2][2]=z_zta;
143
144 if (detJ) {
145 *detJ = x_xi*(y_eta*z_zta - y_zta*z_eta)
146 - x_eta*(y_xi*z_zta - y_zta*z_xi)
147 + x_zta*(y_xi*z_eta - y_eta*z_xi);
148 }
149
150 PROFILE_FUNCTION_END;
151
152 PetscFunctionReturn(0);
153}
154
155
156#undef __FUNCT__
157#define __FUNCT__ "MetricVelocityContravariant"
158/* ------------------------------------------------------------------------- */
159/**
160 * @brief Implementation of \ref MetricVelocityContravariant().
161 * @details Full API contract (arguments, ownership, side effects) is documented with
162 * the header declaration in `include/Metric.h`.
163 * @see MetricVelocityContravariant()
164 */
165PetscErrorCode MetricVelocityContravariant(const PetscReal J[3][3], PetscReal detJ,
166 const PetscReal u[3], PetscReal uc[3])
167{
168 PetscFunctionBeginUser;
169
170 PROFILE_FUNCTION_BEGIN;
171
172 /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
173 PetscReal gxi[3] = { J[1][1]*J[2][2]-J[1][2]*J[2][1],
174 -J[0][1]*J[2][2]+J[0][2]*J[2][1],
175 J[0][1]*J[1][2]-J[0][2]*J[1][1] };
176 PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
177 J[0][0]*J[2][2]-J[0][2]*J[2][0],
178 -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
179 PetscReal gzta[3] = { J[1][0]*J[2][1]-J[1][1]*J[2][0],
180 -J[0][0]*J[2][1]+J[0][1]*J[2][0],
181 J[0][0]*J[1][1]-J[0][1]*J[1][0] };
182
183 PetscReal invDet = 1.0 / detJ;
184 for (int d=0; d<3; ++d) { gxi[d] *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
185
186 uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
187 uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
188 uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
189
190 PROFILE_FUNCTION_END;
191
192 PetscFunctionReturn(0);
193}
194
195/////////////////////////////////////////////////////////////////////////////////
196#undef __FUNCT__
197#define __FUNCT__ "InvertCovariantMetricTensor"
198/**
199 * @brief Internal helper implementation: `InvertCovariantMetricTensor()`.
200 * @details Local to this translation unit.
201 */
202PetscErrorCode InvertCovariantMetricTensor(double covariantTensor[3][3], double contravariantTensor[3][3])
203{
204 PetscFunctionBeginUser;
205
206 const double a11=covariantTensor[0][0], a12=covariantTensor[0][1], a13=covariantTensor[0][2];
207 const double a21=covariantTensor[1][0], a22=covariantTensor[1][1], a23=covariantTensor[1][2];
208 const double a31=covariantTensor[2][0], a32=covariantTensor[2][1], a33=covariantTensor[2][2];
209
210 double det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
211
212 /* Degeneracy has to be judged against the rows' own magnitude, not against a fixed
213 number. The rows here are face-area vectors, so the determinant carries their
214 product - for an orthogonal cell it is exactly the square of the cell volume, and
215 it shrinks as the mesh is refined while the cell stays perfectly well formed. An
216 absolute floor therefore rejects fine grids: a wall-resolved cell of 5e-9 volume
217 has a determinant near 3e-17. The ratio below is 1 for an orthogonal cell and
218 approaches 0 only as the rows become coplanar, which is the actual failure. */
219 const double row_scale = sqrt(a11*a11 + a12*a12 + a13*a13) *
220 sqrt(a21*a21 + a22*a22 + a23*a23) *
221 sqrt(a31*a31 + a32*a32 + a33*a33);
222
223 if (row_scale <= 0.0 || fabs(det) <= 1.0e-10 * row_scale) {
224 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT,
225 "Matrix is singular: |det| = %g is degenerate against the rows' own scale "
226 "%g (ratio %g); the rows are coplanar or one of them vanishes.",
227 (double)fabs(det), (double)row_scale,
228 (double)(row_scale > 0.0 ? fabs(det)/row_scale : 0.0));
229 }
230
231 contravariantTensor[0][0] = (a33*a22-a32*a23)/det;
232 contravariantTensor[0][1] = -(a33*a12-a32*a13)/det;
233 contravariantTensor[0][2] = (a23*a12-a22*a13)/det;
234 contravariantTensor[1][0] = -(a33*a21-a31*a23)/det;
235 contravariantTensor[1][1] = (a33*a11-a31*a13)/det;
236 contravariantTensor[1][2] = -(a23*a11-a21*a13)/det;
237 contravariantTensor[2][0] = (a32*a21-a31*a22)/det;
238 contravariantTensor[2][1] = -(a32*a11-a31*a12)/det;
239 contravariantTensor[2][2] = (a22*a11-a21*a12)/det;
240
241 PetscFunctionReturn(0);
242}
243
244
245#undef __FUNCT__
246#define __FUNCT__ "CalculateFaceNormalAndArea"
247/**
248 * @brief Internal helper implementation: `CalculateFaceNormalAndArea()`.
249 * @details Local to this translation unit.
250 */
251PetscErrorCode CalculateFaceNormalAndArea(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
252{
253 PetscFunctionBeginUser;
254 PROFILE_FUNCTION_BEGIN;
255 double g[3][3];
256 double G[3][3];
257
258 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
259 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
260 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
261
262 PetscCall(InvertCovariantMetricTensor(g, G));
263 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
264 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
265 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
266
267 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
268 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
269 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
270
271 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
272 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
273 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
274
275 *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
276 *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
277 *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
278
279 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
280 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
281 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
282
283 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
284 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
285 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
286
287 PROFILE_FUNCTION_END;
288 PetscFunctionReturn(0);
289}
290
291#undef __FUNCT__
292#define __FUNCT__ "ComputeCellCharacteristicLengthScale"
293/**
294 * @brief Internal helper implementation: `ComputeCellCharacteristicLengthScale()`.
295 * @details Local to this translation unit.
296 */
297PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
298{
299 PetscFunctionBeginUser;
300 PROFILE_FUNCTION_BEGIN;
301 double ni[3], nj[3], nk[3];
302 double Li, Lj, Lk;
303 double Ai, Aj, Ak;
304 double vol = 1./ajc;
305
306 PetscCall(CalculateFaceNormalAndArea(csi, eta, zet, ni, nj, nk, &Ai, &Aj, &Ak));
307 Li = vol / Ai;
308 Lj = vol / Aj;
309 Lk = vol / Ak;
310
311 // Length scale vector = di * ni_vector + dj * nj_vector + dk * nk_vector
312 *dx = fabs( Li * ni[0] + Lj * nj[0] + Lk * nk[0] );
313 *dy = fabs( Li * ni[1] + Lj * nj[1] + Lk * nk[1] );
314 *dz = fabs( Li * ni[2] + Lj * nj[2] + Lk * nk[2] );
315
316 PROFILE_FUNCTION_END;
317 PetscFunctionReturn(0);
318}
319
320#undef __FUNCT__
321#define __FUNCT__ "ComputeCellDirectionalExtents"
322/**
323 * @brief Implementation of \ref ComputeCellDirectionalExtents().
324 * @details Full API contract is documented with the header declaration in
325 * `include/Metric.h`.
326 * @see ComputeCellDirectionalExtents()
327 */
328PetscErrorCode ComputeCellDirectionalExtents(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet,
329 double *l_xi, double *l_eta, double *l_zeta)
330{
331 const double area_xi = sqrt(csi.x*csi.x + csi.y*csi.y + csi.z*csi.z);
332 const double area_eta = sqrt(eta.x*eta.x + eta.y*eta.y + eta.z*eta.z);
333 const double area_zeta = sqrt(zet.x*zet.x + zet.y*zet.y + zet.z*zet.z);
334
335 PetscFunctionBeginUser;
336 PetscCheck(ajc > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
337 "Cell Jacobian must be positive to define cell extents; received %g.", (double)ajc);
338 PetscCheck(area_xi > 0.0 && area_eta > 0.0 && area_zeta > 0.0, PETSC_COMM_SELF,
339 PETSC_ERR_ARG_OUTOFRANGE,
340 "A cell face has zero area (%g, %g, %g), so the cell has no extent across it.",
341 area_xi, area_eta, area_zeta);
342
343 /* Volume over the area of a face pair is the distance between those faces: the
344 cell's extent in that grid direction. It uses nothing but the cell's own metrics,
345 so it is the same wherever and however the cell sits in space. */
346 *l_xi = (1.0/ajc)/area_xi;
347 *l_eta = (1.0/ajc)/area_eta;
348 *l_zeta = (1.0/ajc)/area_zeta;
349 PetscFunctionReturn(0);
350}
351
352#undef __FUNCT__
353#define __FUNCT__ "ComputeCellEdgeVectors"
354/**
355 * @brief Implementation of \ref ComputeCellEdgeVectors().
356 * @details Full API contract is documented with the header declaration in
357 * `include/Metric.h`.
358 * @see ComputeCellEdgeVectors()
359 */
360PetscErrorCode ComputeCellEdgeVectors(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet,
361 Cmpnts edges[3])
362{
363 PetscFunctionBeginUser;
364 PetscCheck(ajc > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
365 "Cell Jacobian must be positive to define cell edges; received %g.", (double)ajc);
366 /* The face-area vectors are J grad(xi), J grad(eta), J grad(zeta), with J the cell
367 volume. The covariant basis dx/dxi equals J grad(eta) x grad(zeta), so it is
368 (eta x zet) / J - the cell's edge along xi, direction and length together. */
369 edges[0].x = ajc * (eta.y*zet.z - eta.z*zet.y);
370 edges[0].y = ajc * (eta.z*zet.x - eta.x*zet.z);
371 edges[0].z = ajc * (eta.x*zet.y - eta.y*zet.x);
372 edges[1].x = ajc * (zet.y*csi.z - zet.z*csi.y);
373 edges[1].y = ajc * (zet.z*csi.x - zet.x*csi.z);
374 edges[1].z = ajc * (zet.x*csi.y - zet.y*csi.x);
375 edges[2].x = ajc * (csi.y*eta.z - csi.z*eta.y);
376 edges[2].y = ajc * (csi.z*eta.x - csi.x*eta.z);
377 edges[2].z = ajc * (csi.x*eta.y - csi.y*eta.x);
378 PetscFunctionReturn(0);
379}
380
381
382#undef __FUNCT__
383#define __FUNCT__ "CheckAndFixGridOrientation"
384/* -------------------------------------------------------------------------- */
385/**
386 * @brief Internal helper implementation: `CheckAndFixGridOrientation()`.
387 * @details Local to this translation unit.
388 */
389PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
390{
391 PetscErrorCode ierr;
392 PetscReal aj_min, aj_max;
393 PetscMPIInt rank;
394
395 PetscFunctionBeginUser;
396
397 PROFILE_FUNCTION_BEGIN;
398
399 /* ---------------- step 1: global extrema of Aj ---------------- */
400 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
401
402 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr); /* already global */
403 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
404
405 LOG_ALLOW(GLOBAL, LOG_INFO,
406 "[orientation] Global Aj range: [%.3e , %.3e]\n",
407 (double)aj_min, (double)aj_max);
408
409 /* ---------------- step 2: detect malformed mesh ---------------- */
410 if (aj_min < 0.0 && aj_max > 0.0)
411 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
412 "Mixed Jacobian signs detected – grid is topologically inconsistent.");
413
414 /* A uniformly left-handed grid is refused. The former "repair" negated Csi, Eta,
415 Zet and Aj on this level only, left the face and cell-centre metrics as they
416 were, and then reset the orientation flag to +1; a mirrored duct solved that way
417 reached forty times its bulk velocity. Renumbering one logical axis is exact. */
418 PetscInt orientation = +1;
419 if (aj_max < 0.0) {
420 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER,
421 "The grid is left-handed: every cell Jacobian is negative, as a mirrored or "
422 "oddly permuted grid is. The solver needs right-handed logical axes. Renumber "
423 "one logical axis (grid.gen: append reverse:axis=i|j|k to the transforms); "
424 "that axis's two faces trade names.");
425 }
426
427 /* ---------------- step 4: store result in UserCtx -------------- */
428 user->GridOrientation = orientation;
429
430 if (!rank)
431 LOG_ALLOW(LOCAL, LOG_INFO, "[orientation] Grid confirmed right-handed.\n");
432
433 PROFILE_FUNCTION_END;
434
435 PetscFunctionReturn(0);
436}
437
438#undef __FUNCT__
439#define __FUNCT__ "ApplyPeriodicCorrectionsToCellCentersAndSpacing"
440/**
441 * @brief Internal helper implementation: `ApplyPeriodicCorrectionsToCellCentersAndSpacing()`.
442 * @details Local to this translation unit.
443 */
444PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing(UserCtx *user)
445{
446 PetscErrorCode ierr;
447 DMDALocalInfo info = user->info;
448 PetscInt xs = info.xs, xe = info.xs + info.xm;
449 PetscInt ys = info.ys, ye = info.ys + info.ym;
450 PetscInt zs = info.zs, ze = info.zs + info.zm;
451 PetscInt mx = info.mx, my = info.my, mz = info.mz;
452 Cmpnts ***cent, ***lcent, ***gs;
453 PetscReal delta;
454
455 PetscFunctionBeginUser;
456 PROFILE_FUNCTION_BEGIN;
457
458 // Check if any periodic boundaries exist
459 PetscBool has_periodic = PETSC_FALSE;
460 for (int i = 0; i < 6; i++) {
461 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
462 has_periodic = PETSC_TRUE;
463 break;
464 }
465 }
466
467 if (!has_periodic) {
468 LOG_ALLOW(LOCAL, LOG_TRACE, "No periodic boundaries; skipping corrections for Cent/GridSpace.\n");
469 PROFILE_FUNCTION_END;
470 PetscFunctionReturn(0);
471 }
472
473 LOG_ALLOW(LOCAL, LOG_DEBUG, "Applying periodic corrections to Cent and GridSpace.\n");
474
475 // Must update ghosts first before applying corrections
476 ierr = UpdateLocalGhosts(user, FIELD_ID_CENT); CHKERRQ(ierr);
477 ierr = UpdateLocalGhosts(user, FIELD_ID_GRID_SPACE); CHKERRQ(ierr);
478
479 // --- X-direction periodic corrections ---
480 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
481 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
482
483 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
484 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
485 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
486
487 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
488 if (user->cgrid) {
489 for (PetscInt k=zs; k<ze; k++) {
490 for (PetscInt j=ys; j<ye; j++) {
491 cent[k][j][0] = lcent[k][j][-2];
492 }
493 }
494 } else {
495 for (PetscInt k=zs; k<ze; k++) {
496 for (PetscInt j=ys; j<ye; j++) {
497 delta = (gs[k][j][1].x + gs[k][j][-2].x) / 2.0;
498 cent[k][j][0].x = cent[k][j][1].x - delta;
499 cent[k][j][0].y = cent[k][j][1].y;
500 cent[k][j][0].z = cent[k][j][1].z;
501 }
502 }
503 }
504 }
505
506 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
507 if (user->cgrid) {
508 for (PetscInt k=zs; k<ze; k++) {
509 for (PetscInt j=ys; j<ye; j++) {
510 cent[k][j][mx-1] = lcent[k][j][mx+1];
511 }
512 }
513 } else {
514 for (PetscInt k=zs; k<ze; k++) {
515 for (PetscInt j=ys; j<ye; j++) {
516 delta = (gs[k][j][mx-2].x + gs[k][j][mx+1].x) / 2.0;
517 cent[k][j][mx-1].x = cent[k][j][mx-2].x + delta;
518 cent[k][j][mx-1].y = cent[k][j][mx-2].y;
519 cent[k][j][mx-1].z = cent[k][j][mx-2].z;
520 }
521 }
522 }
523 }
524
525 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
526 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
527 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
528 }
529 // --- Y-direction periodic corrections ---
530 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
531 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
532
533 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
534 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
535 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
536
537 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
538 if (user->cgrid) {
539 for (PetscInt k=zs; k<ze; k++) {
540 for (PetscInt i=xs; i<xe; i++) {
541 cent[k][0][i] = lcent[k][-2][i];
542 }
543 }
544 } else {
545 for (PetscInt k=zs; k<ze; k++) {
546 for (PetscInt i=xs; i<xe; i++) {
547 delta = (gs[k][1][i].y + gs[k][-2][i].y) / 2.0;
548 cent[k][0][i].x = cent[k][1][i].x;
549 cent[k][0][i].y = cent[k][1][i].y - delta;
550 cent[k][0][i].z = cent[k][1][i].z;
551 }
552 }
553 }
554 }
555
556 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
557 if (user->cgrid) {
558 for (PetscInt k=zs; k<ze; k++) {
559 for (PetscInt i=xs; i<xe; i++) {
560 cent[k][my-1][i] = lcent[k][my+1][i];
561 }
562 }
563 } else {
564 for (PetscInt k=zs; k<ze; k++) {
565 for (PetscInt i=xs; i<xe; i++) {
566 delta = (gs[k][my-2][i].y + gs[k][my+1][i].y) / 2.0;
567 cent[k][my-1][i].x = cent[k][my-2][i].x;
568 cent[k][my-1][i].y = cent[k][my-2][i].y + delta;
569 cent[k][my-1][i].z = cent[k][my-2][i].z;
570 }
571 }
572 }
573 }
574
575 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
576 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
577 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
578
579 }
580
581 // --- Z-direction periodic corrections ---
582 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
583 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
584
585 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
586 ierr = DMDAVecGetArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
587 ierr = DMDAVecGetArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
588
589 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
590 if (user->cgrid) {
591 for (PetscInt j=ys; j<ye; j++) {
592 for (PetscInt i=xs; i<xe; i++) {
593 cent[0][j][i] = lcent[-2][j][i];
594 }
595 }
596 } else {
597 for (PetscInt j=ys; j<ye; j++) {
598 for (PetscInt i=xs; i<xe; i++) {
599 delta = (gs[1][j][i].z + gs[-2][j][i].z) / 2.0;
600 cent[0][j][i].x = cent[1][j][i].x;
601 cent[0][j][i].y = cent[1][j][i].y;
602 cent[0][j][i].z = cent[1][j][i].z - delta;
603 }
604 }
605 }
606 }
607
608 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
609 if (user->cgrid) {
610 for (PetscInt j=ys; j<ye; j++) {
611 for (PetscInt i=xs; i<xe; i++) {
612 cent[mz-1][j][i] = lcent[mz+1][j][i];
613 }
614 }
615 } else {
616 for (PetscInt j=ys; j<ye; j++) {
617 for (PetscInt i=xs; i<xe; i++) {
618 delta = (gs[mz-2][j][i].z + gs[mz+1][j][i].z) / 2.0;
619 cent[mz-1][j][i].x = cent[mz-2][j][i].x;
620 cent[mz-1][j][i].y = cent[mz-2][j][i].y;
621 cent[mz-1][j][i].z = cent[mz-2][j][i].z + delta;
622 }
623 }
624 }
625 }
626
627 ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace, &gs); CHKERRQ(ierr);
628 ierr = DMDAVecRestoreArray(user->fda, user->lCent, &lcent); CHKERRQ(ierr);
629 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
630
631 }
632
633 PROFILE_FUNCTION_END;
634 PetscFunctionReturn(0);
635}
636
637#undef __FUNCT__
638#define __FUNCT__ "ApplyPeriodicCorrectionsToIFaceCenter"
639/**
640 * @brief Internal helper implementation: `ApplyPeriodicCorrectionsToIFaceCenter()`.
641 * @details Local to this translation unit.
642 */
643PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter(UserCtx *user)
644{
645 const FieldId fields[] = {FIELD_ID_CENTX};
646
647 PetscFunctionBeginUser;
648 PetscCall(SynchronizePeriodicFaceFields(user, 'i', 1, fields));
649 PetscFunctionReturn(0);
650}
651
652
653#undef __FUNCT__
654#define __FUNCT__ "ApplyPeriodicCorrectionsToJFaceCenter"
655/**
656 * @brief Internal helper implementation: `ApplyPeriodicCorrectionsToJFaceCenter()`.
657 * @details Local to this translation unit.
658 */
659PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
660{
661 const FieldId fields[] = {FIELD_ID_CENTY};
662
663 PetscFunctionBeginUser;
664 PetscCall(SynchronizePeriodicFaceFields(user, 'j', 1, fields));
665 PetscFunctionReturn(0);
666}
667
668
669#undef __FUNCT__
670#define __FUNCT__ "ApplyPeriodicCorrectionsToKFaceCenter"
671/**
672 * @brief Internal helper implementation: `ApplyPeriodicCorrectionsToKFaceCenter()`.
673 * @details Local to this translation unit.
674 */
675PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
676{
677 const FieldId fields[] = {FIELD_ID_CENTZ};
678
679 PetscFunctionBeginUser;
680 PetscCall(SynchronizePeriodicFaceFields(user, 'k', 1, fields));
681 PetscFunctionReturn(0);
682}
683
684#undef __FUNCT__
685#define __FUNCT__ "ComputeFaceMetrics"
686/**
687 * @brief Internal helper implementation: `ComputeFaceMetrics()`.
688 * @details Local to this translation unit.
689 */
690PetscErrorCode ComputeFaceMetrics(UserCtx *user)
691{
692 PetscErrorCode ierr;
693 DMDALocalInfo info;
694 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
695 Cmpnts ***nodal_coords_arr;
696 Vec localCoords_from_dm;
697
698 PetscFunctionBeginUser;
699
700 PROFILE_FUNCTION_BEGIN;
701
702 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
703
704 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
705
706 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
707 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
708 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
709 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
710
711 // --- 2. Get arrays for output global Vecs from UserCtx ---
712 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
713 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
714 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
715
716 // Define owned node ranges (global indices)
717 PetscInt xs = info.xs, xe = info.xs + info.xm;
718 PetscInt ys = info.ys, ye = info.ys + info.ym;
719 PetscInt zs = info.zs, ze = info.zs + info.zm;
720
721 // Global domain dimensions (total number of nodes)
722 PetscInt mx = info.mx;
723 PetscInt my = info.my;
724 PetscInt mz = info.mz;
725
726 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
727 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
728 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
729 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
730 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
731
732 // These represent the surface area of the curvilinear cell face and the normal rotated such that the direction of increasing coordinate is maintained.
733 // The metric vectors (Csi, Eta, Zet) are defined to point in the direction of their corresponding increasing computational coordinate.
734
735 // Calculate Csi
736 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
737 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
738 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
739
740 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);
741 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);
742 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);
743 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);
744 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);
745 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);
746
747 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
748 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
749 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
750 }
751 }
752 }
753
754 // Calculate Eta
755 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
756 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
757 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
758
759 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);
760 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);
761 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);
762 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);
763 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);
764 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);
765
766 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
767 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
768 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
769 }
770 }
771 }
772
773 // Calculate Zet
774 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
775 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
776 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
777
778 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);
779 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);
780 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);
781 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);
782 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);
783 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);
784
785 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
786 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
787 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
788 }
789 }
790 }
791
792 // --- 4. Boundary Extrapolation ---
793 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
794 PetscInt i_bnd, j_bnd, k_bnd;
795
796 if (xs == 0) { // If this rank owns the global i=0 boundary
797 i_bnd = 0;
798 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
799 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
800 if (i_bnd + 1 < mx) {
801 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
802 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
803 }
804 }
805 }
806 }
807 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
808 i_bnd = mx - 1;
809 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
810 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
811 if (i_bnd - 1 >= 0) {
812 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
813 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
814 }
815 }
816 }
817 }
818 if (ys == 0) {
819 j_bnd = 0;
820 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
821 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
822 if (j_bnd + 1 < my) {
823 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
824 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
825 }
826 }
827 }
828 }
829 if (ye == my) {
830 j_bnd = my - 1;
831 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
832 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
833 if (j_bnd - 1 >= 0) {
834 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
835 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
836 }
837 }
838 }
839 }
840 if (zs == 0) {
841 k_bnd = 0;
842 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
843 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
844 if (k_bnd + 1 < mz) {
845 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
846 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
847 }
848 }
849 }
850 }
851 if (ze == mz) {
852 k_bnd = mz - 1;
853 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
854 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
855 if (k_bnd - 1 >= 0) {
856 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
857 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
858 }
859 }
860 }
861 }
862
863 if (info.xs==0 && info.ys==0 && info.zs==0) {
864 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
865 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
866 }
867
868 // --- 5. Restore all arrays ---
869 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
870 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
871 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
872 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
873
874 // --- 6. Assemble Global Vectors ---
875 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
876 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
877 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
878 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
879
880 // --- 7. Update Local Ghosted Versions ---
881 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
882 ierr = UpdateLocalGhosts(user, FIELD_ID_CSI); CHKERRQ(ierr);
883 ierr = UpdateLocalGhosts(user, FIELD_ID_ETA); CHKERRQ(ierr);
884 ierr = UpdateLocalGhosts(user, FIELD_ID_ZET); CHKERRQ(ierr);
885
886 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
887
888 PROFILE_FUNCTION_END;
889
890 PetscFunctionReturn(0);
891}
892
893
894#undef __FUNCT__
895#define __FUNCT__ "ComputeCellCenteredJacobianInverse"
896/**
897 * @brief Implementation of \ref ComputeCellCenteredJacobianInverse().
898 * @details Full API contract (arguments, ownership, side effects) is documented with
899 * the header declaration in `include/Metric.h`.
900 * @see ComputeCellCenteredJacobianInverse()
901 */
902PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
903{
904 PetscErrorCode ierr;
905 DMDALocalInfo info;
906 PetscScalar ***aj_arr;
907 Cmpnts ***nodal_coords_arr;
908 Vec localCoords_from_dm;
909
910 PetscFunctionBeginUser;
911 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
912
913 // --- 1. Get Nodal Coordinates and Output Array ---
914 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
915 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
916 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
917 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
918
919 // Define owned node ranges (global indices)
920 PetscInt xs = info.xs, xe = info.xs + info.xm;
921 PetscInt ys = info.ys, ye = info.ys + info.ym;
922 PetscInt zs = info.zs, ze = info.zs + info.zm;
923
924 // Global domain dimensions (total number of nodes)
925 PetscInt mx = info.mx;
926 PetscInt my = info.my;
927 PetscInt mz = info.mz;
928
929 // --- 2. Calculate Aj for INTERIOR Stencils ---
930
931 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
932 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
933 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
934
935 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
936 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
937 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
938
939 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
940 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
941 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
942
943 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) );
944
945 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) );
946
947 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) );
948
949 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) );
950
951 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) );
952
953 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) );
954
955 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) );
956
957 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) );
958
959 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) );
960
961 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);
962 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
963 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
964 }
965 }
966 }
967
968 // --- 4. Boundary Extrapolation for Aj ---
969 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
970 PetscInt i_bnd, j_bnd, k_bnd;
971
972 if (xs == 0) {
973 i_bnd = 0;
974 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
975 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
976 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
977 }
978 }
979 }
980 if (xe == mx) {
981 i_bnd = mx - 1;
982 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
983 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
984 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
985 }
986 }
987 }
988 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
989 if (ys == 0) {
990 j_bnd = 0;
991 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
992 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
993 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
994 }
995 }
996 }
997 if (ye == my) {
998 j_bnd = my - 1;
999 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
1000 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1001 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
1002 }
1003 }
1004 }
1005 if (zs == 0) {
1006 k_bnd = 0;
1007 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1008 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1009 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
1010 }
1011 }
1012 }
1013 if (ze == mz) {
1014 k_bnd = mz - 1;
1015 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
1016 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
1017 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
1018 }
1019 }
1020 }
1021
1022 // --- 5. Restore arrays ---
1023 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
1024 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
1025
1026 // --- 6. Assemble Global Vector ---
1027 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
1028 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
1029 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
1030
1031 // --- 7. Update Local Ghosted Version ---
1032 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
1033 ierr = UpdateLocalGhosts(user, FIELD_ID_AJ); CHKERRQ(ierr);
1034
1035 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
1036 PetscFunctionReturn(0);
1037}
1038
1039
1040#undef __FUNCT__
1041#define __FUNCT__ "ComputeCellCentersAndSpacing"
1042/**
1043 * @brief Internal helper implementation: `ComputeCellCentersAndSpacing()`.
1044 * @details Local to this translation unit.
1045 */
1046PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
1047{
1048 PetscErrorCode ierr;
1049 DMDALocalInfo info;
1050 Vec lCoords;
1051 const Cmpnts ***coor;
1052 Cmpnts ***cent, ***gs;
1053 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
1054 PetscInt xs,ys,zs,xe,ye,ze,mx,my,mz;
1055
1056 PetscFunctionBeginUser;
1057
1058 PROFILE_FUNCTION_BEGIN;
1059
1060 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);
1061
1062 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1063 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1064 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1065
1066 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1067 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1068
1069 xs = info.xs; xe = info.xs + info.xm;
1070 ys = info.ys; ye = info.ys + info.ym;
1071 zs = info.zs; ze = info.zs + info.zm;
1072 mx = info.mx; my = info.my; mz = info.mz;
1073
1074 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
1075 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
1076 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
1077
1078 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
1079 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
1080 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
1081
1082 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
1083 for (PetscInt k=k_start_node; k<k_end_node; k++) {
1084 for (PetscInt j=j_start_node; j<j_end_node; j++) {
1085 for (PetscInt i=i_start_node; i<i_end_node; i++) {
1086 // Calculate cell center as the average of its 8 corner nodes
1087 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);
1088 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);
1089 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);
1090
1091 // Calculate Grid Spacing in i-direction (distance between i-face centers)
1092 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);
1093 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);
1094 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);
1095 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);
1096 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);
1097 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);
1098 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1099
1100 // Calculate Grid Spacing in j-direction (distance between j-face centers)
1101 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);
1102 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);
1103 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);
1104 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);
1105 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);
1106 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);
1107 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1108
1109 // Calculate Grid Spacing in k-direction (distance between k-face centers)
1110 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);
1111 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);
1112 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);
1113 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);
1114 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);
1115 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);
1116 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1117 }
1118 }
1119 }
1120
1121 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1122 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1123 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1124
1125 // Assemble and update ghost regions for the new data
1126 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1127 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
1128 ierr = UpdateLocalGhosts(user, FIELD_ID_CENT); CHKERRQ(ierr);
1129 ierr = UpdateLocalGhosts(user, FIELD_ID_GRID_SPACE); CHKERRQ(ierr);
1130
1131 ierr = ApplyPeriodicCorrectionsToCellCentersAndSpacing(user); CHKERRQ(ierr);
1132
1133 // Final assembly and ghost update after corrections
1134 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr);
1135 ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1136 ierr = UpdateLocalGhosts(user, FIELD_ID_CENT); CHKERRQ(ierr);
1137
1138 PROFILE_FUNCTION_END;
1139
1140 PetscFunctionReturn(0);
1141}
1142
1143#undef __FUNCT__
1144#define __FUNCT__ "ComputeIFaceMetrics"
1145/**
1146 * @brief Internal helper implementation: `ComputeIFaceMetrics()`.
1147 * @details Local to this translation unit.
1148 */
1149PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
1150{
1151 PetscErrorCode ierr;
1152 DMDALocalInfo info;
1153 Vec lCoords;
1154 const Cmpnts ***coor;
1155 Cmpnts ***centx; //***gs;
1156 const Cmpnts ***centx_const;
1157 Cmpnts ***icsi, ***ieta, ***izet;
1158 PetscScalar ***iaj;
1159 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1160
1161 PetscFunctionBeginUser;
1162
1163 PROFILE_FUNCTION_BEGIN;
1164
1165 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1166
1167 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1168 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1169 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1170 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1171 PetscInt lxe = xe;
1172 PetscInt lys = ys; PetscInt lye = ye;
1173 PetscInt lzs = zs; PetscInt lze = ze;
1174
1175 if (ys==0) lys = ys+1;
1176 if (zs==0) lzs = zs+1;
1177
1178 if (xe==mx) lxe=xe-1;
1179 if (ye==my) lye=ye-1;
1180 if (ze==mz) lze=ze-1;
1181
1182 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1183 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1184 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1185 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1186 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1187
1188 //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);
1189
1190 // Populate only owned physical face centers. Periodic endpoint and ghost
1191 // coordinates are established by the canonical face-field synchronizer.
1192 for (PetscInt k = PetscMax(zs, 1); k < PetscMin(ze, mz - 1); k++) {
1193 for (PetscInt j = PetscMax(ys, 1); j < PetscMin(ye, my - 1); j++) {
1194 for (PetscInt i = xs; i < PetscMin(xe, mx - 1); i++) {
1195 //----- DEBUG ------
1196 //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);
1197 //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,
1198 // coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
1199 // coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
1200 // coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
1201 // coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
1202
1203 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);
1204 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);
1205 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);
1206
1207 //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);
1208 }
1209 }
1210 }
1211
1212 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
1213 /*
1214 if(xs==0){
1215 for(PetscInt k=gzs+1;k < gze; k++){
1216 for(PetscInt j=gys+1;j < gye; j++){
1217 PetscInt i=0;
1218 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
1219 centx[k][j][i-1].y=centx[k][j][i].y;
1220 centx[k][j][i-1].z=centx[k][j][i].z;
1221 }
1222 }
1223 }
1224 if (xe==mx){
1225 for(PetscInt k=gzs+1; k<gze; k++) {
1226 for (PetscInt j=gys+1; j<gye;j++) {
1227 PetscInt i=mx-1;
1228 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
1229 centx[k][j][i].y=centx[k][j][i-1].y;
1230 centx[k][j][i].z=centx[k][j][i-1].z;
1231 }
1232 }
1233 }
1234 */
1235
1236 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1237 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1238
1239 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1240
1241 {
1242 const FieldId face_centers[] = {FIELD_ID_CENTX};
1243 ierr = SynchronizePeriodicFaceFields(user, 'i', 1, face_centers); CHKERRQ(ierr);
1244 }
1245
1246 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1247
1248 // --- Part 2: Calculate metrics using face-centered coordinates ---
1249 ierr = DMDAVecGetArrayRead(user->fda, user->lCentx, &centx_const); CHKERRQ(ierr);
1250 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1251 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1252 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1253 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1254
1255 // Loop over the OWNED region where we will store the final metrics
1256 for (PetscInt k=lzs; k<lze; k++) {
1257 for (PetscInt j=lys; j<lye; j++) {
1258 for (PetscInt i=xs; i<lxe; i++) {
1259
1260 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1261 if(i == 0 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC){
1262 // Forward difference at the domain's min-i boundary
1263 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
1264 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
1265 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
1266 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1267 // Backward difference at the domain's max-i boundary
1268 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
1269 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
1270 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
1271 } else { // Central difference in the interior (or if PERIODIC BCs)
1272 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
1273 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
1274 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
1275 }
1276
1277 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1278 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1279 // Forward difference
1280 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
1281 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
1282 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
1283 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1284 // Backward difference
1285 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
1286 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
1287 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
1288 } else { // Central difference (interior or PERIODIC)
1289 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
1290 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
1291 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
1292 }
1293
1294 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1295 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1296 // Forward difference
1297 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
1298 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
1299 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
1300 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1301 // Backward difference
1302 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
1303 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
1304 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
1305 } else { // Central difference (Interior + PERIODIC)
1306 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
1307 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
1308 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
1309 }
1310
1311 // --- Metric calculations (identical to legacy FormMetrics) ---
1312 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1313 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1314 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1315
1316 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1317 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1318 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1319
1320 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
1321 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1322 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
1323
1324 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
1325 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
1326 iaj[k][j][i] = 1.0 / iaj[k][j][i];
1327 }
1328 }
1329 }
1330 }
1331
1332 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentx, &centx_const); CHKERRQ(ierr);
1333 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1334 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1335 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1336 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1337
1338 // --- Part 3: Assemble global vectors and update local ghosts ---
1339 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
1340 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
1341 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
1342 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
1343
1344 ierr = UpdateLocalGhosts(user, FIELD_ID_ICSI); CHKERRQ(ierr);
1345 ierr = UpdateLocalGhosts(user, FIELD_ID_IETA); CHKERRQ(ierr);
1346 ierr = UpdateLocalGhosts(user, FIELD_ID_IZET); CHKERRQ(ierr);
1347 ierr = UpdateLocalGhosts(user, FIELD_ID_IAJ); CHKERRQ(ierr);
1348
1349 PROFILE_FUNCTION_END;
1350
1351 PetscFunctionReturn(0);
1352}
1353
1354#undef __FUNCT__
1355#define __FUNCT__ "ComputeJFaceMetrics"
1356/**
1357 * @brief Internal helper implementation: `ComputeJFaceMetrics()`.
1358 * @details Local to this translation unit.
1359 */
1360PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
1361{
1362 PetscErrorCode ierr;
1363 DMDALocalInfo info;
1364 Vec lCoords;
1365 const Cmpnts ***coor;
1366 Cmpnts ***centy; //***gs;
1367 const Cmpnts ***centy_const;
1368 Cmpnts ***jcsi, ***jeta, ***jzet;
1369 PetscScalar ***jaj;
1370 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1371
1372 PetscFunctionBeginUser;
1373
1374 PROFILE_FUNCTION_BEGIN;
1375
1376 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1377
1378 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1379 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1380 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1381 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1382 PetscInt lxs = xs; PetscInt lxe = xe;
1383 PetscInt lye = ye;
1384 PetscInt lzs = zs; PetscInt lze = ze;
1385
1386 if (xs==0) lxs = xs+1;
1387 if (zs==0) lzs = zs+1;
1388
1389 if (xe==mx) lxe=xe-1;
1390 if (ye==my) lye=ye-1;
1391 if (ze==mz) lze=ze-1;
1392
1393 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1394 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1395 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1396 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1397 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1398
1399 for (PetscInt k = PetscMax(zs, 1); k < PetscMin(ze, mz - 1); k++) {
1400 for (PetscInt j = ys; j < PetscMin(ye, my - 1); j++) {
1401 for (PetscInt i = PetscMax(xs, 1); i < PetscMin(xe, mx - 1); i++) {
1402 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);
1403 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);
1404 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);
1405 }
1406 }
1407 }
1408
1409 /*
1410 if(ys==0){
1411 for(PetscInt k=gzs+1;k < gze; k++){
1412 for(PetscInt i=gxs+1;j < gxe; i++){
1413 PetscInt j=0;
1414 centy[k][j-1][i].x=centy[k][j][i].x;
1415 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1416 centy[k][j-1][i].z=centy[k][j][i].z;
1417 }
1418 }
1419 }
1420 if (ye==my){
1421 for(PetscInt k=gzs+1; k<gze; k++) {
1422 for (PetscInt i=gxs+1; j<gxe;i++) {
1423 PetscInt j=my-1;
1424 centy[k][j][i].x=centy[k][j-1][i].x
1425 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1426 centy[k][j][i].z=centy[k][j-1][i].z;
1427 }
1428 }
1429 }
1430 */
1431
1432 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1433 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1434 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1435
1436 {
1437 const FieldId face_centers[] = {FIELD_ID_CENTY};
1438 ierr = SynchronizePeriodicFaceFields(user, 'j', 1, face_centers); CHKERRQ(ierr);
1439 }
1440
1441 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-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->lCenty, &centy_const); CHKERRQ(ierr);
1445 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1446 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1447 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1448 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1449
1450 // Loop over the OWNED region where we will store the final metrics
1451 for (PetscInt k=lzs; k<lze; k++) {
1452 for (PetscInt j=ys; 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 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1457 // Forward difference at the domain's min-i boundary
1458 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1459 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1460 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1461 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1462 // Backward difference at the domain's max-i boundary
1463 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1464 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1465 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1466 } else { // Central difference in the interior or PERIODIC
1467 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1468 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1469 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1470 }
1471
1472 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1473 if (j == 0 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1474 // Forward difference
1475 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1476 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1477 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1478 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1479 // Backward difference
1480 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1481 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1482 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1483 } else { // Central difference (interior or PERIODIC)
1484 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1485 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1486 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1487 }
1488
1489 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1490 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1491 // Forward difference
1492 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1493 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1494 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1495 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1496 // Backward difference
1497 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1498 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1499 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1500 } else { // Central difference (Interior or PERIODIC)
1501 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1502 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1503 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1504 }
1505
1506 // --- Metric calculations (identical to legacy FormMetrics) ---
1507 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1508 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1509 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1510
1511 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1512 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1513 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1514
1515 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1516 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1517 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1518
1519 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1520 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1521 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1522 }
1523 }
1524 }
1525 }
1526
1527 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCenty, &centy_const); CHKERRQ(ierr);
1528 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1529 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1530 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1531 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1532
1533 // --- Part 3: Assemble global vectors and update local ghosts ---
1534 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1535 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1536 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1537 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1538
1539 ierr = UpdateLocalGhosts(user, FIELD_ID_JCSI); CHKERRQ(ierr);
1540 ierr = UpdateLocalGhosts(user, FIELD_ID_JETA); CHKERRQ(ierr);
1541 ierr = UpdateLocalGhosts(user, FIELD_ID_JZET); CHKERRQ(ierr);
1542 ierr = UpdateLocalGhosts(user, FIELD_ID_JAJ); CHKERRQ(ierr);
1543
1544 PROFILE_FUNCTION_END;
1545
1546 PetscFunctionReturn(0);
1547}
1548
1549#undef __FUNCT__
1550#define __FUNCT__ "ComputeJFaceMetrics"
1551/**
1552 * @brief Internal helper implementation: `ComputeKFaceMetrics()`.
1553 * @details Local to this translation unit.
1554 */
1555PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
1556{
1557 PetscErrorCode ierr;
1558 DMDALocalInfo info;
1559 Vec lCoords;
1560 const Cmpnts ***coor;
1561 Cmpnts ***centz; //***gs;
1562 const Cmpnts ***centz_const;
1563 Cmpnts ***kcsi, ***keta, ***kzet;
1564 PetscScalar ***kaj;
1565 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1566
1567 PetscFunctionBeginUser;
1568
1569 PROFILE_FUNCTION_BEGIN;
1570
1571 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1572
1573 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1574 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1575 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1576 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1577 PetscInt lxs = xs; PetscInt lxe = xe;
1578 PetscInt lys = ys; PetscInt lye = ye;
1579 PetscInt lze = ze;
1580
1581 if (xs==0) lxs = xs+1;
1582 if (ys==0) lys = ys+1;
1583
1584 if (xe==mx) lxe=xe-1;
1585 if (ye==my) lye=ye-1;
1586 if (ze==mz) lze=ze-1;
1587
1588 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1589 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1590 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1591 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1592 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1593
1594 for (PetscInt k = zs; k < PetscMin(ze, mz - 1); k++) {
1595 for (PetscInt j = PetscMax(ys, 1); j < PetscMin(ye, my - 1); j++) {
1596 for (PetscInt i = PetscMax(xs, 1); i < PetscMin(xe, mx - 1); i++) {
1597 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);
1598 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);
1599 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);
1600 }
1601 }
1602 }
1603
1604 /*
1605 if(zs==0){
1606 for(PetscInt j=gys+1;j < gye; j++){
1607 for(PetscInt i=gxs+1;j < gxe; i++){
1608 PetscInt k=0;
1609 centz[k-1][j][i].x=centz[k][j][i].x;
1610 centz[k-1][j][i].y=centz[k][j][i].y;
1611 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1612 }
1613 }
1614 }
1615 if (ze==mz){
1616 for(PetscInt j=gys+1; j<gye; j++) {
1617 for (PetscInt i=gxs+1; j<gxe;i++) {
1618 PetscInt k=mz-1;
1619 centy[k][j][i].x=centy[k-1][j][i].x
1620 centy[k][j][i].y=centy[k-1][j][i].y;
1621 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1622 }
1623 }
1624 }
1625 */
1626
1627 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1628 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1629 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1630
1631 {
1632 const FieldId face_centers[] = {FIELD_ID_CENTZ};
1633 ierr = SynchronizePeriodicFaceFields(user, 'k', 1, face_centers); CHKERRQ(ierr);
1634 }
1635
1636 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1637
1638 // --- Part 2: Calculate metrics using face-centered coordinates ---
1639 ierr = DMDAVecGetArrayRead(user->fda, user->lCentz, &centz_const); CHKERRQ(ierr);
1640 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1641 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1642 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1643 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1644
1645 // Loop over the OWNED region where we will store the final metrics
1646 for (PetscInt k=zs; k<lze; k++) {
1647 for (PetscInt j=lys; j<lye; j++) {
1648 for (PetscInt i=lxs; i<lxe; i++) {
1649
1650 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1651 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1652 // Forward difference at the domain's min-i boundary
1653 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
1654 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
1655 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
1656 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1657 // Backward difference at the domain's max-i boundary
1658 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
1659 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
1660 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
1661 } else { // Central difference in the interior (or PERIODIC)
1662 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
1663 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
1664 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
1665 }
1666
1667 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1668 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1669 // Forward difference
1670 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
1671 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
1672 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
1673 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1674 // Backward difference
1675 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
1676 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
1677 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
1678 } else { // Central difference (interior or PERIODIC)
1679 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
1680 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
1681 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
1682 }
1683
1684 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1685 if (k == 0 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1686 // Forward difference
1687 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
1688 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
1689 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
1690 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1691 // Backward difference
1692 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
1693 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
1694 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
1695 } else { // Central difference (Interior or PERIODIC)
1696 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
1697 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
1698 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
1699 }
1700
1701 // --- Metric calculations (identical to legacy FormMetrics) ---
1702 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1703 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1704 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1705
1706 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1707 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1708 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1709
1710 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1711 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1712 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1713
1714 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
1715 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
1716 kaj[k][j][i] = 1.0 / kaj[k][j][i];
1717 }
1718 }
1719 }
1720 }
1721
1722 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentz, &centz_const); CHKERRQ(ierr);
1723 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1724 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1725 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1726 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1727
1728 // --- Part 3: Assemble global vectors and update local ghosts ---
1729 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
1730 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
1731 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
1732 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
1733
1734 ierr = UpdateLocalGhosts(user, FIELD_ID_KCSI); CHKERRQ(ierr);
1735 ierr = UpdateLocalGhosts(user, FIELD_ID_KETA); CHKERRQ(ierr);
1736 ierr = UpdateLocalGhosts(user, FIELD_ID_KZET); CHKERRQ(ierr);
1737 ierr = UpdateLocalGhosts(user, FIELD_ID_KAJ); CHKERRQ(ierr);
1738
1739 PROFILE_FUNCTION_END;
1740
1741 PetscFunctionReturn(0);
1742}
1743
1744/**
1745 * @brief Convert logical cell indices into the flattened global cell identifier.
1746 */
1747static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
1748{
1749 PetscInt nidx;
1750 DMDALocalInfo info = user->info;
1751
1752 PetscInt mx = info.mx, my = info.my;
1753
1754 AO ao;
1755 DMDAGetAO(user->da, &ao);
1756 nidx=i+j*mx+k*mx*my;
1757
1758 AOApplicationToPetsc(ao,1,&nidx);
1759
1760 return (nidx);
1761}
1762
1763#undef __FUNCT__
1764#define __FUNCT__ "ComputeMetricsDivergence"
1765/**
1766 * @brief Internal helper implementation: `ComputeMetricsDivergence()`.
1767 * @details Local to this translation unit.
1768 */
1769PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
1770{
1771 DM da = user->da, fda = user->fda;
1772 DMDALocalInfo info = user->info;
1773 PetscInt xs = info.xs, xe = info.xs + info.xm;
1774 PetscInt ys = info.ys, ye = info.ys + info.ym;
1775 PetscInt zs = info.zs, ze = info.zs + info.zm;
1776 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1777 PetscInt lxs, lys, lzs, lxe, lye, lze;
1778 PetscInt i, j, k;
1779 Vec Div;
1780 PetscReal ***div, ***aj;
1781 Cmpnts ***csi, ***eta, ***zet;
1782 PetscReal maxdiv;
1783
1784 PetscFunctionBeginUser;
1785
1786 PROFILE_FUNCTION_BEGIN;
1787
1788 lxs = xs; lxe = xe;
1789 lys = ys; lye = ye;
1790 lzs = zs; lze = ze;
1791
1792 if (xs == 0) lxs = xs + 1;
1793 if (ys == 0) lys = ys + 1;
1794 if (zs == 0) lzs = zs + 1;
1795
1796 if (xe == mx) lxe = xe - 1;
1797 if (ye == my) lye = ye - 1;
1798 if (ze == mz) lze = ze - 1;
1799
1800 DMDAVecGetArray(fda, user->lCsi, &csi);
1801 DMDAVecGetArray(fda, user->lEta, &eta);
1802 DMDAVecGetArray(fda, user->lZet, &zet);
1803 DMDAVecGetArray(da, user->lAj, &aj);
1804
1805 VecDuplicate(user->P, &Div);
1806 VecSet(Div, 0.);
1807 DMDAVecGetArray(da, Div, &div);
1808
1809 for (k = lzs; k < lze; k++) {
1810 for (j = lys; j < lye; j++) {
1811 for (i = lxs; i < lxe; i++) {
1812 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1813 eta[k][j][i].x - eta[k][j-1][i].x +
1814 zet[k][j][i].x - zet[k-1][j][i].x +
1815 csi[k][j][i].y - csi[k][j][i-1].y +
1816 eta[k][j][i].y - eta[k][j-1][i].y +
1817 zet[k][j][i].y - zet[k-1][j][i].y +
1818 csi[k][j][i].z - csi[k][j][i-1].z +
1819 eta[k][j][i].z - eta[k][j-1][i].z +
1820 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1821 div[k][j][i] = fabs(divergence);
1822 }
1823 }
1824 }
1825
1826 DMDAVecRestoreArray(da, Div, &div);
1827
1828 PetscInt MaxFlatIndex = -1;
1829 VecMax(Div, &MaxFlatIndex, &maxdiv);
1830 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %" PetscInt_FMT ".\n",maxdiv,MaxFlatIndex);
1831
1832 for (k=zs; k<ze; k++) {
1833 for (j=ys; j<ye; j++) {
1834 for (i=xs; i<xe; i++) {
1835 if (Gidx(i,j,k,user) == MaxFlatIndex) {
1836 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,(int)k,(int)j,(int)i);
1837 }
1838 }
1839 }
1840 }
1841
1842
1843 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1844 DMDAVecRestoreArray(fda, user->lEta, &eta);
1845 DMDAVecRestoreArray(fda, user->lZet, &zet);
1846 DMDAVecRestoreArray(da, user->lAj, &aj);
1847 VecDestroy(&Div);
1848
1849
1850 PROFILE_FUNCTION_END;
1851
1852 PetscFunctionReturn(0);
1853}
1854
1855#undef __FUNCT__
1856#define __FUNCT__ "ComputeMetricNorms"
1857/**
1858 * @brief Internal helper implementation: `ComputeMetricNorms()`.
1859 * @details Local to this translation unit.
1860 */
1861PetscErrorCode ComputeMetricNorms(UserCtx *user)
1862{
1863
1864 DMDALocalInfo info = user->info;
1865 PetscInt xs = info.xs, xe = info.xs + info.xm;
1866 PetscInt ys = info.ys, ye = info.ys + info.ym;
1867 PetscInt zs = info.zs, ze = info.zs + info.zm;
1868 PetscInt i, j, k;
1869
1870 PetscFunctionBeginUser;
1871
1872 PROFILE_FUNCTION_BEGIN;
1873
1874 PetscReal CsiMax, EtaMax, ZetMax;
1875 PetscReal ICsiMax, IEtaMax, IZetMax;
1876 PetscReal JCsiMax, JEtaMax, JZetMax;
1877 PetscReal KCsiMax, KEtaMax, KZetMax;
1878 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
1879
1880 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
1881 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
1882 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
1883 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
1884 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
1885
1886 // Max Values
1887 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
1888 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
1889 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
1890
1891 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
1892 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
1893 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
1894
1895 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
1896 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
1897 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
1898
1899 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
1900 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
1901 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
1902
1903 VecMax(user->lAj,&AjMaxArg,&AjMax);
1904 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1905 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1906 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1907
1908 VecMax(user->lAj,&AjMaxArg,&AjMax);
1909 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1910 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1911 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1912
1913 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
1914
1915 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
1916 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
1917 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
1918 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
1919 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
1920
1921 for (k=zs; k<ze; k++) {
1922 for (j=ys; j<ye; j++) {
1923 for (i=xs; i<xe; i++) {
1924 if (Gidx(i,j,k,user) == CsiMaxArg) {
1925 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
1926 }
1927 if (Gidx(i,j,k,user) == EtaMaxArg) {
1928 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
1929 }
1930 if (Gidx(i,j,k,user) == ZetMaxArg) {
1931 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
1932 }
1933 if (Gidx(i,j,k,user) == ICsiMaxArg) {
1934 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
1935 }
1936 if (Gidx(i,j,k,user) == IEtaMaxArg) {
1937 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
1938 }
1939 if (Gidx(i,j,k,user) == IZetMaxArg) {
1940 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
1941 }
1942 if (Gidx(i,j,k,user) == JCsiMaxArg) {
1943 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
1944 }
1945 if (Gidx(i,j,k,user) == JEtaMaxArg) {
1946 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
1947 }
1948 if (Gidx(i,j,k,user) == JZetMaxArg) {
1949 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
1950 }
1951 if (Gidx(i,j,k,user) == KCsiMaxArg) {
1952 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
1953 }
1954 if (Gidx(i,j,k,user) == KEtaMaxArg) {
1955 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
1956 }
1957 if (Gidx(i,j,k,user) == KZetMaxArg) {
1958 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
1959 }
1960 if (Gidx(i,j,k,user) == AjMaxArg) {
1961 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
1962 }
1963 if (Gidx(i,j,k,user) == IAjMaxArg) {
1964 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
1965 }
1966 if (Gidx(i,j,k,user) == JAjMaxArg) {
1967 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
1968 }
1969 if (Gidx(i,j,k,user) == KAjMaxArg) {
1970 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
1971 }
1972 }
1973 }
1974 }
1975
1976 /*
1977 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
1978 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
1979 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
1980 */
1981
1982 PROFILE_FUNCTION_END;
1983
1984 PetscFunctionReturn(0);
1985}
1986
1987#undef __FUNCT__
1988#define __FUNCT__ "ComputeAllGridMetrics"
1989/**
1990 * @brief Internal helper implementation: `CalculateAllGridMetrics()`.
1991 * @details Local to this translation unit.
1992 */
1993PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
1994{
1995 PetscErrorCode ierr;
1996 UserMG *usermg = &simCtx->usermg;
1997 MGCtx *mgctx = usermg->mgctx;
1998 PetscInt nblk = simCtx->block_number;
1999
2000 PetscFunctionBeginUser;
2001
2002 PROFILE_FUNCTION_BEGIN;
2003
2004 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
2005
2006 // Loop through all levels and all blocks
2007 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
2008 for (PetscInt bi = 0; bi < nblk; bi++) {
2009 UserCtx *user = &mgctx[level].user[bi];
2010 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
2011
2012 // Call the modern, modular helper functions for each UserCtx.
2013 // These functions are self-contained and operate on the data within the provided context.
2014 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
2015 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
2016 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
2017 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
2018 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
2019 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
2020 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
2021
2022 // Apply Periodic Boundary Condition Adjustments if necessary
2023 ierr = ApplyMetricsPeriodicBCs(user); CHKERRQ(ierr);
2024 // Diagnostics
2025 ierr = ComputeMetricNorms(user);
2026 if (level == usermg->mglevels - 1) {
2027 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
2028 }
2029 }
2030 }
2031
2032 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
2033
2034 PROFILE_FUNCTION_END;
2035
2036 PetscFunctionReturn(0);
2037}
2038/* ------------------------------------------------------------------------- */
2039/* End of Metric.c */
ApplyMetricsPeriodicBCs
PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
(Orchestrator) Updates all metric-related fields in the local ghost cell regions for periodic boundar...
Definition Boundaries.c:2166
SynchronizePeriodicFaceFields
PetscErrorCode SynchronizePeriodicFaceFields(UserCtx *user, char face_direction, PetscInt num_fields, const FieldId field_ids[])
Synchronizes persistent fields belonging to one face family.
Definition Boundaries.c:1881
ComputeCellEdgeVectors
PetscErrorCode ComputeCellEdgeVectors(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, Cmpnts edges[3])
Implementation of ComputeCellEdgeVectors().
Definition Metric.c:360
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Internal helper implementation: CalculateAllGridMetrics().
Definition Metric.c:1993
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)
Implementation of MetricJacobian().
Definition Metric.c:105
InvertCovariantMetricTensor
PetscErrorCode InvertCovariantMetricTensor(double covariantTensor[3][3], double contravariantTensor[3][3])
Internal helper implementation: InvertCovariantMetricTensor().
Definition Metric.c:202
ComputeMetricNorms
PetscErrorCode ComputeMetricNorms(UserCtx *user)
Internal helper implementation: ComputeMetricNorms().
Definition Metric.c:1861
ApplyPeriodicCorrectionsToCellCentersAndSpacing
PetscErrorCode ApplyPeriodicCorrectionsToCellCentersAndSpacing(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToCellCentersAndSpacing().
Definition Metric.c:444
TrilinearBlend
static void TrilinearBlend(const Cmpnts V[8], PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Blend eight corner values at the supplied trilinear reference coordinates.
Definition Metric.c:50
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Convert logical cell indices into the flattened global cell identifier.
Definition Metric.c:1747
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Internal helper implementation: CheckAndFixGridOrientation().
Definition Metric.c:389
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Internal helper implementation: ComputeCellCentersAndSpacing().
Definition Metric.c:1046
MetricGetCellVertices
PetscErrorCode MetricGetCellVertices(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
Implementation of MetricGetCellVertices().
Definition Metric.c:26
ComputeJFaceMetrics
PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeJFaceMetrics().
Definition Metric.c:1360
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Internal helper implementation: ComputeMetricsDivergence().
Definition Metric.c:1769
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeFaceMetrics().
Definition Metric.c:690
ComputeCellCharacteristicLengthScale
PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
Internal helper implementation: ComputeCellCharacteristicLengthScale().
Definition Metric.c:297
CalculateFaceNormalAndArea
PetscErrorCode CalculateFaceNormalAndArea(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
Internal helper implementation: CalculateFaceNormalAndArea().
Definition Metric.c:251
ApplyPeriodicCorrectionsToIFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToIFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToIFaceCenter().
Definition Metric.c:643
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Implementation of ComputeCellCenteredJacobianInverse().
Definition Metric.c:902
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeKFaceMetrics().
Definition Metric.c:1555
ApplyPeriodicCorrectionsToJFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToJFaceCenter().
Definition Metric.c:659
ComputeCellDirectionalExtents
PetscErrorCode ComputeCellDirectionalExtents(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *l_xi, double *l_eta, double *l_zeta)
Implementation of ComputeCellDirectionalExtents().
Definition Metric.c:328
MetricLogicalToPhysical
PetscErrorCode MetricLogicalToPhysical(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Implementation of MetricLogicalToPhysical().
Definition Metric.c:76
ComputeIFaceMetrics
PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeIFaceMetrics().
Definition Metric.c:1149
ApplyPeriodicCorrectionsToKFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToKFaceCenter().
Definition Metric.c:675
MetricVelocityContravariant
PetscErrorCode MetricVelocityContravariant(const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
Implementation of MetricVelocityContravariant().
Definition Metric.c:165
FieldId
FieldId
Compile-time identity for a catalogued Eulerian field.
Definition field_catalog.h:23
FIELD_ID_JETA
@ FIELD_ID_JETA
Definition field_catalog.h:51
FIELD_ID_CENTZ
@ FIELD_ID_CENTZ
Definition field_catalog.h:46
FIELD_ID_CSI
@ FIELD_ID_CSI
Definition field_catalog.h:37
FIELD_ID_IAJ
@ FIELD_ID_IAJ
Definition field_catalog.h:56
FIELD_ID_KETA
@ FIELD_ID_KETA
Definition field_catalog.h:54
FIELD_ID_JAJ
@ FIELD_ID_JAJ
Definition field_catalog.h:57
FIELD_ID_KAJ
@ FIELD_ID_KAJ
Definition field_catalog.h:58
FIELD_ID_AJ
@ FIELD_ID_AJ
Definition field_catalog.h:41
FIELD_ID_CENTY
@ FIELD_ID_CENTY
Definition field_catalog.h:45
FIELD_ID_GRID_SPACE
@ FIELD_ID_GRID_SPACE
Definition field_catalog.h:43
FIELD_ID_KZET
@ FIELD_ID_KZET
Definition field_catalog.h:55
FIELD_ID_IETA
@ FIELD_ID_IETA
Definition field_catalog.h:48
FIELD_ID_ICSI
@ FIELD_ID_ICSI
Definition field_catalog.h:47
FIELD_ID_ETA
@ FIELD_ID_ETA
Definition field_catalog.h:38
FIELD_ID_CENT
@ FIELD_ID_CENT
Definition field_catalog.h:42
FIELD_ID_JCSI
@ FIELD_ID_JCSI
Definition field_catalog.h:50
FIELD_ID_JZET
@ FIELD_ID_JZET
Definition field_catalog.h:52
FIELD_ID_IZET
@ FIELD_ID_IZET
Definition field_catalog.h:49
FIELD_ID_ZET
@ FIELD_ID_ZET
Definition field_catalog.h:39
FIELD_ID_KCSI
@ FIELD_ID_KCSI
Definition field_catalog.h:53
FIELD_ID_CENTX
@ FIELD_ID_CENTX
Definition field_catalog.h:44
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:298
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:253
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:894
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:33
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:34
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:885
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:2489
UserCtx::lCent
Vec lCent
Definition variables.h:1148
PERIODIC
@ PERIODIC
Definition variables.h:318
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:1148
MGCtx::user
UserCtx * user
Definition variables.h:729
SimCtx::rank
PetscMPIInt rank
Definition variables.h:862
UserCtx::cgrid
PetscInt cgrid
Definition variables.h:1094
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:1099
SimCtx::block_number
PetscInt block_number
Definition variables.h:952
UserCtx::lIEta
Vec lIEta
Definition variables.h:1151
UserCtx::lIZet
Vec lIZet
Definition variables.h:1151
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:1077
UserCtx::Centz
Vec Centz
Definition variables.h:1149
SimCtx::usermg
UserMG usermg
Definition variables.h:1015
UserCtx::_this
PetscInt _this
Definition variables.h:1092
UserCtx::lKEta
Vec lKEta
Definition variables.h:1153
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:1152
Cmpnts::x
PetscScalar x
Definition variables.h:122
UserCtx::Centx
Vec Centx
Definition variables.h:1149
UserCtx::lKZet
Vec lKZet
Definition variables.h:1153
UserCtx::lJEta
Vec lJEta
Definition variables.h:1152
UserCtx::lGridSpace
Vec lGridSpace
Definition variables.h:1148
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:1165
Cmpnts::z
PetscScalar z
Definition variables.h:122
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:1153
UserCtx::lCenty
Vec lCenty
Definition variables.h:1150
UserMG::mglevels
PetscInt mglevels
Definition variables.h:736
UserCtx::lJZet
Vec lJZet
Definition variables.h:1152
UserCtx::lCentx
Vec lCentx
Definition variables.h:1150
UserCtx::lICsi
Vec lICsi
Definition variables.h:1151
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:1092
UserCtx::info
DMDALocalInfo info
Definition variables.h:1086
Cmpnts::y
PetscScalar y
Definition variables.h:122
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:739
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:392
UserCtx::Centy
Vec Centy
Definition variables.h:1149
UserCtx::lCentz
Vec lCentz
Definition variables.h:1150
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:288
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:290
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:289
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:290
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:288
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:289
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:121
MGCtx
Context for Multigrid operations.
Definition variables.h:728
SimCtx
The master context for the entire simulation.
Definition variables.h:859
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:1074
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:735
Generated by doxygen 1.9.8