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