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