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 const char *fields[] = {"Centx"};
594
595 PetscFunctionBeginUser;
596 PetscCall(SynchronizePeriodicFaceFields(user, 'i', 1, fields));
597 PetscFunctionReturn(0);
598}
599
600
601#undef __FUNCT__
602#define __FUNCT__ "ApplyPeriodicCorrectionsToJFaceCenter"
603/**
604 * @brief Internal helper implementation: `ApplyPeriodicCorrectionsToJFaceCenter()`.
605 * @details Local to this translation unit.
606 */
607PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
608{
609 const char *fields[] = {"Centy"};
610
611 PetscFunctionBeginUser;
612 PetscCall(SynchronizePeriodicFaceFields(user, 'j', 1, fields));
613 PetscFunctionReturn(0);
614}
615
616
617#undef __FUNCT__
618#define __FUNCT__ "ApplyPeriodicCorrectionsToKFaceCenter"
619/**
620 * @brief Internal helper implementation: `ApplyPeriodicCorrectionsToKFaceCenter()`.
621 * @details Local to this translation unit.
622 */
623PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
624{
625 const char *fields[] = {"Centz"};
626
627 PetscFunctionBeginUser;
628 PetscCall(SynchronizePeriodicFaceFields(user, 'k', 1, fields));
629 PetscFunctionReturn(0);
630}
631
632#undef __FUNCT__
633#define __FUNCT__ "ComputeFaceMetrics"
634/**
635 * @brief Internal helper implementation: `ComputeFaceMetrics()`.
636 * @details Local to this translation unit.
637 */
638PetscErrorCode ComputeFaceMetrics(UserCtx *user)
639{
640 PetscErrorCode ierr;
641 DMDALocalInfo info;
642 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
643 Cmpnts ***nodal_coords_arr;
644 Vec localCoords_from_dm;
645
646 PetscFunctionBeginUser;
647
648 PROFILE_FUNCTION_BEGIN;
649
650 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
651
652 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
653
654 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
655 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
656 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
657 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
658
659 // --- 2. Get arrays for output global Vecs from UserCtx ---
660 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
661 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
662 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
663
664 // Define owned node ranges (global indices)
665 PetscInt xs = info.xs, xe = info.xs + info.xm;
666 PetscInt ys = info.ys, ye = info.ys + info.ym;
667 PetscInt zs = info.zs, ze = info.zs + info.zm;
668
669 // Global domain dimensions (total number of nodes)
670 PetscInt mx = info.mx;
671 PetscInt my = info.my;
672 PetscInt mz = info.mz;
673
674 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
675 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
676 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
677 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
678 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
679
680 // These represent the surface area of the curvilinear cell face and the normal rotated such that the direction of increasing coordinate is maintained.
681 // The metric vectors (Csi, Eta, Zet) are defined to point in the direction of their corresponding increasing computational coordinate.
682
683 // Calculate Csi
684 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
685 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
686 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
687
688 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);
689 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);
690 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);
691 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);
692 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);
693 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);
694
695 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
696 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
697 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
698 }
699 }
700 }
701
702 // Calculate Eta
703 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
704 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
705 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
706
707 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);
708 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);
709 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);
710 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);
711 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);
712 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);
713
714 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
715 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
716 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
717 }
718 }
719 }
720
721 // Calculate Zet
722 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
723 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
724 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
725
726 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);
727 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);
728 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);
729 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);
730 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);
731 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);
732
733 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
734 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
735 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
736 }
737 }
738 }
739
740 // --- 4. Boundary Extrapolation ---
741 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
742 PetscInt i_bnd, j_bnd, k_bnd;
743
744 if (xs == 0) { // If this rank owns the global i=0 boundary
745 i_bnd = 0;
746 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
747 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
748 if (i_bnd + 1 < mx) {
749 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
750 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
751 }
752 }
753 }
754 }
755 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
756 i_bnd = mx - 1;
757 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
758 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
759 if (i_bnd - 1 >= 0) {
760 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
761 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
762 }
763 }
764 }
765 }
766 if (ys == 0) {
767 j_bnd = 0;
768 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
769 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
770 if (j_bnd + 1 < my) {
771 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
772 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
773 }
774 }
775 }
776 }
777 if (ye == my) {
778 j_bnd = my - 1;
779 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
780 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
781 if (j_bnd - 1 >= 0) {
782 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
783 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
784 }
785 }
786 }
787 }
788 if (zs == 0) {
789 k_bnd = 0;
790 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
791 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
792 if (k_bnd + 1 < mz) {
793 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
794 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
795 }
796 }
797 }
798 }
799 if (ze == mz) {
800 k_bnd = mz - 1;
801 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
802 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
803 if (k_bnd - 1 >= 0) {
804 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
805 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
806 }
807 }
808 }
809 }
810
811 if (info.xs==0 && info.ys==0 && info.zs==0) {
812 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
813 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
814 }
815
816 // --- 5. Restore all arrays ---
817 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
818 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
819 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
820 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
821
822 // --- 6. Assemble Global Vectors ---
823 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
824 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
825 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
826 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
827
828 // --- 7. Update Local Ghosted Versions ---
829 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
830 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
831 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
832 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
833
834 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
835
836 PROFILE_FUNCTION_END;
837
838 PetscFunctionReturn(0);
839}
840
841
842#undef __FUNCT__
843#define __FUNCT__ "ComputeCellCenteredJacobianInverse"
844/**
845 * @brief Implementation of \ref ComputeCellCenteredJacobianInverse().
846 * @details Full API contract (arguments, ownership, side effects) is documented with
847 * the header declaration in `include/Metric.h`.
848 * @see ComputeCellCenteredJacobianInverse()
849 */
850PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
851{
852 PetscErrorCode ierr;
853 DMDALocalInfo info;
854 PetscScalar ***aj_arr;
855 Cmpnts ***nodal_coords_arr;
856 Vec localCoords_from_dm;
857
858 PetscFunctionBeginUser;
859 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
860
861 // --- 1. Get Nodal Coordinates and Output Array ---
862 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
863 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
864 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
865 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
866
867 // Define owned node ranges (global indices)
868 PetscInt xs = info.xs, xe = info.xs + info.xm;
869 PetscInt ys = info.ys, ye = info.ys + info.ym;
870 PetscInt zs = info.zs, ze = info.zs + info.zm;
871
872 // Global domain dimensions (total number of nodes)
873 PetscInt mx = info.mx;
874 PetscInt my = info.my;
875 PetscInt mz = info.mz;
876
877 // --- 2. Calculate Aj for INTERIOR Stencils ---
878
879 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
880 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
881 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
882
883 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
884 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
885 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
886
887 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
888 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
889 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
890
891 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) );
892
893 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) );
894
895 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) );
896
897 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) );
898
899 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) );
900
901 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) );
902
903 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) );
904
905 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) );
906
907 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) );
908
909 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);
910 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
911 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
912 }
913 }
914 }
915
916 // --- 4. Boundary Extrapolation for Aj ---
917 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
918 PetscInt i_bnd, j_bnd, k_bnd;
919
920 if (xs == 0) {
921 i_bnd = 0;
922 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
923 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
924 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
925 }
926 }
927 }
928 if (xe == mx) {
929 i_bnd = mx - 1;
930 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
931 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
932 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
933 }
934 }
935 }
936 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
937 if (ys == 0) {
938 j_bnd = 0;
939 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
940 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
941 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
942 }
943 }
944 }
945 if (ye == my) {
946 j_bnd = my - 1;
947 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
948 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
949 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
950 }
951 }
952 }
953 if (zs == 0) {
954 k_bnd = 0;
955 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
956 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
957 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
958 }
959 }
960 }
961 if (ze == mz) {
962 k_bnd = mz - 1;
963 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
964 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
965 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
966 }
967 }
968 }
969
970 // --- 5. Restore arrays ---
971 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
972 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
973
974 // --- 6. Assemble Global Vector ---
975 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
976 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
977 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
978
979 // --- 7. Update Local Ghosted Version ---
980 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
981 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
982
983 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
984 PetscFunctionReturn(0);
985}
986
987
988#undef __FUNCT__
989#define __FUNCT__ "ComputeCellCentersAndSpacing"
990/**
991 * @brief Internal helper implementation: `ComputeCellCentersAndSpacing()`.
992 * @details Local to this translation unit.
993 */
994PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
995{
996 PetscErrorCode ierr;
997 DMDALocalInfo info;
998 Vec lCoords;
999 const Cmpnts ***coor;
1000 Cmpnts ***cent, ***gs;
1001 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
1002 PetscInt xs,ys,zs,xe,ye,ze,mx,my,mz;
1003
1004 PetscFunctionBeginUser;
1005
1006 PROFILE_FUNCTION_BEGIN;
1007
1008 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);
1009
1010 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1011 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1012 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1013
1014 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1015 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1016
1017 xs = info.xs; xe = info.xs + info.xm;
1018 ys = info.ys; ye = info.ys + info.ym;
1019 zs = info.zs; ze = info.zs + info.zm;
1020 mx = info.mx; my = info.my; mz = info.mz;
1021
1022 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
1023 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
1024 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
1025
1026 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
1027 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
1028 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
1029
1030 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
1031 for (PetscInt k=k_start_node; k<k_end_node; k++) {
1032 for (PetscInt j=j_start_node; j<j_end_node; j++) {
1033 for (PetscInt i=i_start_node; i<i_end_node; i++) {
1034 // Calculate cell center as the average of its 8 corner nodes
1035 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);
1036 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);
1037 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);
1038
1039 // Calculate Grid Spacing in i-direction (distance between i-face centers)
1040 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);
1041 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);
1042 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);
1043 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);
1044 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);
1045 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);
1046 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1047
1048 // Calculate Grid Spacing in j-direction (distance between j-face centers)
1049 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);
1050 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);
1051 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);
1052 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);
1053 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);
1054 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);
1055 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1056
1057 // Calculate Grid Spacing in k-direction (distance between k-face centers)
1058 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);
1059 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);
1060 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);
1061 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);
1062 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);
1063 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);
1064 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
1065 }
1066 }
1067 }
1068
1069 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1070 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
1071 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
1072
1073 // Assemble and update ghost regions for the new data
1074 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1075 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
1076 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
1077 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
1078
1079 ierr = ApplyPeriodicCorrectionsToCellCentersAndSpacing(user); CHKERRQ(ierr);
1080
1081 // Final assembly and ghost update after corrections
1082 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr);
1083 ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
1084 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
1085
1086 PROFILE_FUNCTION_END;
1087
1088 PetscFunctionReturn(0);
1089}
1090
1091#undef __FUNCT__
1092#define __FUNCT__ "ComputeIFaceMetrics"
1093/**
1094 * @brief Internal helper implementation: `ComputeIFaceMetrics()`.
1095 * @details Local to this translation unit.
1096 */
1097PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
1098{
1099 PetscErrorCode ierr;
1100 DMDALocalInfo info;
1101 Vec lCoords;
1102 const Cmpnts ***coor;
1103 Cmpnts ***centx; //***gs;
1104 const Cmpnts ***centx_const;
1105 Cmpnts ***icsi, ***ieta, ***izet;
1106 PetscScalar ***iaj;
1107 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1108
1109 PetscFunctionBeginUser;
1110
1111 PROFILE_FUNCTION_BEGIN;
1112
1113 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1114
1115 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1116 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1117 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1118 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1119 PetscInt lxe = xe;
1120 PetscInt lys = ys; PetscInt lye = ye;
1121 PetscInt lzs = zs; PetscInt lze = ze;
1122
1123 if (ys==0) lys = ys+1;
1124 if (zs==0) lzs = zs+1;
1125
1126 if (xe==mx) lxe=xe-1;
1127 if (ye==my) lye=ye-1;
1128 if (ze==mz) lze=ze-1;
1129
1130 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1131 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1132 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1133 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1134 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1135
1136 //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);
1137
1138 // Populate only owned physical face centers. Periodic endpoint and ghost
1139 // coordinates are established by the canonical face-field synchronizer.
1140 for (PetscInt k = PetscMax(zs, 1); k < PetscMin(ze, mz - 1); k++) {
1141 for (PetscInt j = PetscMax(ys, 1); j < PetscMin(ye, my - 1); j++) {
1142 for (PetscInt i = xs; i < PetscMin(xe, mx - 1); i++) {
1143 //----- DEBUG ------
1144 //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);
1145 //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,
1146 // coor[k][j][i].x, coor[k][j][i].y, coor[k][j][i].z,
1147 // coor[k-1][j][i].x, coor[k-1][j][i].y, coor[k-1][j][i].z,
1148 // coor[k][j-1][i].x, coor[k][j-1][i].y, coor[k][j-1][i].z,
1149 // coor[k-1][j-1][i].x, coor[k-1][j-1][i].y, coor[k-1][j-1][i].z);
1150
1151 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);
1152 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);
1153 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);
1154
1155 //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);
1156 }
1157 }
1158 }
1159
1160 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face center coordinates calculated. \n", user->simCtx->rank);
1161 /*
1162 if(xs==0){
1163 for(PetscInt k=gzs+1;k < gze; k++){
1164 for(PetscInt j=gys+1;j < gye; j++){
1165 PetscInt i=0;
1166 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
1167 centx[k][j][i-1].y=centx[k][j][i].y;
1168 centx[k][j][i-1].z=centx[k][j][i].z;
1169 }
1170 }
1171 }
1172 if (xe==mx){
1173 for(PetscInt k=gzs+1; k<gze; k++) {
1174 for (PetscInt j=gys+1; j<gye;j++) {
1175 PetscInt i=mx-1;
1176 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
1177 centx[k][j][i].y=centx[k][j][i-1].y;
1178 centx[k][j][i].z=centx[k][j][i-1].z;
1179 }
1180 }
1181 }
1182 */
1183
1184 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1185 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
1186
1187 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1188
1189 {
1190 const char *face_centers[] = {"Centx"};
1191 ierr = SynchronizePeriodicFaceFields(user, 'i', 1, face_centers); CHKERRQ(ierr);
1192 }
1193
1194 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1195
1196 // --- Part 2: Calculate metrics using face-centered coordinates ---
1197 ierr = DMDAVecGetArrayRead(user->fda, user->lCentx, &centx_const); CHKERRQ(ierr);
1198 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1199 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1200 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1201 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1202
1203 // Loop over the OWNED region where we will store the final metrics
1204 for (PetscInt k=lzs; k<lze; k++) {
1205 for (PetscInt j=lys; j<lye; j++) {
1206 for (PetscInt i=xs; i<lxe; i++) {
1207
1208 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1209 if(i == 0 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC){
1210 // Forward difference at the domain's min-i boundary
1211 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
1212 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
1213 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
1214 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1215 // Backward difference at the domain's max-i boundary
1216 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
1217 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
1218 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
1219 } else { // Central difference in the interior (or if PERIODIC BCs)
1220 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
1221 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
1222 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
1223 }
1224
1225 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1226 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1227 // Forward difference
1228 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
1229 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
1230 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
1231 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1232 // Backward difference
1233 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
1234 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
1235 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
1236 } else { // Central difference (interior or PERIODIC)
1237 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
1238 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
1239 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
1240 }
1241
1242 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1243 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1244 // Forward difference
1245 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
1246 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
1247 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
1248 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1249 // Backward difference
1250 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
1251 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
1252 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
1253 } else { // Central difference (Interior + PERIODIC)
1254 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
1255 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
1256 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
1257 }
1258
1259 // --- Metric calculations (identical to legacy FormMetrics) ---
1260 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1261 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1262 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1263
1264 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1265 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1266 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1267
1268 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
1269 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1270 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
1271
1272 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
1273 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
1274 iaj[k][j][i] = 1.0 / iaj[k][j][i];
1275 }
1276 }
1277 }
1278 }
1279
1280 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentx, &centx_const); CHKERRQ(ierr);
1281 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
1282 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
1283 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
1284 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
1285
1286 // --- Part 3: Assemble global vectors and update local ghosts ---
1287 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
1288 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
1289 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
1290 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
1291
1292 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
1293 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
1294 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
1295 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
1296
1297 PROFILE_FUNCTION_END;
1298
1299 PetscFunctionReturn(0);
1300}
1301
1302#undef __FUNCT__
1303#define __FUNCT__ "ComputeJFaceMetrics"
1304/**
1305 * @brief Internal helper implementation: `ComputeJFaceMetrics()`.
1306 * @details Local to this translation unit.
1307 */
1308PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
1309{
1310 PetscErrorCode ierr;
1311 DMDALocalInfo info;
1312 Vec lCoords;
1313 const Cmpnts ***coor;
1314 Cmpnts ***centy; //***gs;
1315 const Cmpnts ***centy_const;
1316 Cmpnts ***jcsi, ***jeta, ***jzet;
1317 PetscScalar ***jaj;
1318 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1319
1320 PetscFunctionBeginUser;
1321
1322 PROFILE_FUNCTION_BEGIN;
1323
1324 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1325
1326 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1327 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1328 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1329 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1330 PetscInt lxs = xs; PetscInt lxe = xe;
1331 PetscInt lye = ye;
1332 PetscInt lzs = zs; PetscInt lze = ze;
1333
1334 if (xs==0) lxs = xs+1;
1335 if (zs==0) lzs = zs+1;
1336
1337 if (xe==mx) lxe=xe-1;
1338 if (ye==my) lye=ye-1;
1339 if (ze==mz) lze=ze-1;
1340
1341 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1342 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1343 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1344 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1345 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1346
1347 for (PetscInt k = PetscMax(zs, 1); k < PetscMin(ze, mz - 1); k++) {
1348 for (PetscInt j = ys; j < PetscMin(ye, my - 1); j++) {
1349 for (PetscInt i = PetscMax(xs, 1); i < PetscMin(xe, mx - 1); i++) {
1350 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);
1351 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);
1352 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);
1353 }
1354 }
1355 }
1356
1357 /*
1358 if(ys==0){
1359 for(PetscInt k=gzs+1;k < gze; k++){
1360 for(PetscInt i=gxs+1;j < gxe; i++){
1361 PetscInt j=0;
1362 centy[k][j-1][i].x=centy[k][j][i].x;
1363 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1364 centy[k][j-1][i].z=centy[k][j][i].z;
1365 }
1366 }
1367 }
1368 if (ye==my){
1369 for(PetscInt k=gzs+1; k<gze; k++) {
1370 for (PetscInt i=gxs+1; j<gxe;i++) {
1371 PetscInt j=my-1;
1372 centy[k][j][i].x=centy[k][j-1][i].x
1373 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1374 centy[k][j][i].z=centy[k][j-1][i].z;
1375 }
1376 }
1377 }
1378 */
1379
1380 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1381 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1382 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1383
1384 {
1385 const char *face_centers[] = {"Centy"};
1386 ierr = SynchronizePeriodicFaceFields(user, 'j', 1, face_centers); CHKERRQ(ierr);
1387 }
1388
1389 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1390
1391 // --- Part 2: Calculate metrics using face-centered coordinates ---
1392 ierr = DMDAVecGetArrayRead(user->fda, user->lCenty, &centy_const); CHKERRQ(ierr);
1393 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1394 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1395 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1396 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1397
1398 // Loop over the OWNED region where we will store the final metrics
1399 for (PetscInt k=lzs; k<lze; k++) {
1400 for (PetscInt j=ys; j<lye; j++) {
1401 for (PetscInt i=lxs; i<lxe; i++) {
1402
1403 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1404 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1405 // Forward difference at the domain's min-i boundary
1406 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1407 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1408 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1409 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1410 // Backward difference at the domain's max-i boundary
1411 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1412 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1413 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1414 } else { // Central difference in the interior or PERIODIC
1415 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1416 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1417 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1418 }
1419
1420 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1421 if (j == 0 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1422 // Forward difference
1423 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1424 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1425 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1426 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1427 // Backward difference
1428 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1429 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1430 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1431 } else { // Central difference (interior or PERIODIC)
1432 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1433 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1434 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1435 }
1436
1437 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1438 if (k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1439 // Forward difference
1440 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1441 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1442 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1443 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1444 // Backward difference
1445 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1446 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1447 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1448 } else { // Central difference (Interior or PERIODIC)
1449 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1450 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1451 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1452 }
1453
1454 // --- Metric calculations (identical to legacy FormMetrics) ---
1455 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1456 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1457 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1458
1459 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1460 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1461 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1462
1463 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1464 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1465 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1466
1467 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1468 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1469 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1470 }
1471 }
1472 }
1473 }
1474
1475 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCenty, &centy_const); CHKERRQ(ierr);
1476 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1477 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1478 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1479 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1480
1481 // --- Part 3: Assemble global vectors and update local ghosts ---
1482 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1483 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1484 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1485 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1486
1487 ierr = UpdateLocalGhosts(user, "JCsi"); CHKERRQ(ierr);
1488 ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
1489 ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
1490 ierr = UpdateLocalGhosts(user, "JAj"); CHKERRQ(ierr);
1491
1492 PROFILE_FUNCTION_END;
1493
1494 PetscFunctionReturn(0);
1495}
1496
1497#undef __FUNCT__
1498#define __FUNCT__ "ComputeJFaceMetrics"
1499/**
1500 * @brief Internal helper implementation: `ComputeKFaceMetrics()`.
1501 * @details Local to this translation unit.
1502 */
1503PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
1504{
1505 PetscErrorCode ierr;
1506 DMDALocalInfo info;
1507 Vec lCoords;
1508 const Cmpnts ***coor;
1509 Cmpnts ***centz; //***gs;
1510 const Cmpnts ***centz_const;
1511 Cmpnts ***kcsi, ***keta, ***kzet;
1512 PetscScalar ***kaj;
1513 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1514
1515 PetscFunctionBeginUser;
1516
1517 PROFILE_FUNCTION_BEGIN;
1518
1519 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1520
1521 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1522 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1523 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1524 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1525 PetscInt lxs = xs; PetscInt lxe = xe;
1526 PetscInt lys = ys; PetscInt lye = ye;
1527 PetscInt lze = ze;
1528
1529 if (xs==0) lxs = xs+1;
1530 if (ys==0) lys = ys+1;
1531
1532 if (xe==mx) lxe=xe-1;
1533 if (ye==my) lye=ye-1;
1534 if (ze==mz) lze=ze-1;
1535
1536 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1537 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1538 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1539 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1540 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1541
1542 for (PetscInt k = zs; k < PetscMin(ze, mz - 1); k++) {
1543 for (PetscInt j = PetscMax(ys, 1); j < PetscMin(ye, my - 1); j++) {
1544 for (PetscInt i = PetscMax(xs, 1); i < PetscMin(xe, mx - 1); i++) {
1545 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);
1546 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);
1547 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);
1548 }
1549 }
1550 }
1551
1552 /*
1553 if(zs==0){
1554 for(PetscInt j=gys+1;j < gye; j++){
1555 for(PetscInt i=gxs+1;j < gxe; i++){
1556 PetscInt k=0;
1557 centz[k-1][j][i].x=centz[k][j][i].x;
1558 centz[k-1][j][i].y=centz[k][j][i].y;
1559 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1560 }
1561 }
1562 }
1563 if (ze==mz){
1564 for(PetscInt j=gys+1; j<gye; j++) {
1565 for (PetscInt i=gxs+1; j<gxe;i++) {
1566 PetscInt k=mz-1;
1567 centy[k][j][i].x=centy[k-1][j][i].x
1568 centy[k][j][i].y=centy[k-1][j][i].y;
1569 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1570 }
1571 }
1572 }
1573 */
1574
1575 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1576 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1577 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1578
1579 {
1580 const char *face_centers[] = {"Centz"};
1581 ierr = SynchronizePeriodicFaceFields(user, 'k', 1, face_centers); CHKERRQ(ierr);
1582 }
1583
1584 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1585
1586 // --- Part 2: Calculate metrics using face-centered coordinates ---
1587 ierr = DMDAVecGetArrayRead(user->fda, user->lCentz, &centz_const); CHKERRQ(ierr);
1588 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1589 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1590 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1591 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1592
1593 // Loop over the OWNED region where we will store the final metrics
1594 for (PetscInt k=zs; k<lze; k++) {
1595 for (PetscInt j=lys; j<lye; j++) {
1596 for (PetscInt i=lxs; i<lxe; i++) {
1597
1598 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1599 if (i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
1600 // Forward difference at the domain's min-i boundary
1601 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
1602 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
1603 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
1604 } else if (i == mx - 2 && user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
1605 // Backward difference at the domain's max-i boundary
1606 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
1607 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
1608 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
1609 } else { // Central difference in the interior (or PERIODIC)
1610 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
1611 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
1612 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
1613 }
1614
1615 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1616 if (j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
1617 // Forward difference
1618 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
1619 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
1620 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
1621 } else if (j == my - 2 && user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
1622 // Backward difference
1623 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
1624 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
1625 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
1626 } else { // Central difference (interior or PERIODIC)
1627 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
1628 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
1629 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
1630 }
1631
1632 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1633 if (k == 0 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
1634 // Forward difference
1635 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
1636 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
1637 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
1638 } else if (k == mz - 2 && user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
1639 // Backward difference
1640 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
1641 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
1642 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
1643 } else { // Central difference (Interior or PERIODIC)
1644 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
1645 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
1646 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
1647 }
1648
1649 // --- Metric calculations (identical to legacy FormMetrics) ---
1650 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1651 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1652 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1653
1654 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1655 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1656 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1657
1658 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1659 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1660 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1661
1662 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
1663 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
1664 kaj[k][j][i] = 1.0 / kaj[k][j][i];
1665 }
1666 }
1667 }
1668 }
1669
1670 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentz, &centz_const); CHKERRQ(ierr);
1671 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1672 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1673 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1674 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1675
1676 // --- Part 3: Assemble global vectors and update local ghosts ---
1677 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
1678 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
1679 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
1680 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
1681
1682 ierr = UpdateLocalGhosts(user, "KCsi"); CHKERRQ(ierr);
1683 ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
1684 ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
1685 ierr = UpdateLocalGhosts(user, "KAj"); CHKERRQ(ierr);
1686
1687 PROFILE_FUNCTION_END;
1688
1689 PetscFunctionReturn(0);
1690}
1691
1692/**
1693 * @brief Internal helper implementation: `Gidx()`.
1694 * @details Local to this translation unit.
1695 */
1696static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
1697{
1698 PetscInt nidx;
1699 DMDALocalInfo info = user->info;
1700
1701 PetscInt mx = info.mx, my = info.my;
1702
1703 AO ao;
1704 DMDAGetAO(user->da, &ao);
1705 nidx=i+j*mx+k*mx*my;
1706
1707 AOApplicationToPetsc(ao,1,&nidx);
1708
1709 return (nidx);
1710}
1711
1712#undef __FUNCT__
1713#define __FUNCT__ "ComputeMetricsDivergence"
1714/**
1715 * @brief Internal helper implementation: `ComputeMetricsDivergence()`.
1716 * @details Local to this translation unit.
1717 */
1718PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
1719{
1720 DM da = user->da, fda = user->fda;
1721 DMDALocalInfo info = user->info;
1722 PetscInt xs = info.xs, xe = info.xs + info.xm;
1723 PetscInt ys = info.ys, ye = info.ys + info.ym;
1724 PetscInt zs = info.zs, ze = info.zs + info.zm;
1725 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1726 PetscInt lxs, lys, lzs, lxe, lye, lze;
1727 PetscInt i, j, k;
1728 Vec Div;
1729 PetscReal ***div, ***aj;
1730 Cmpnts ***csi, ***eta, ***zet;
1731 PetscReal maxdiv;
1732
1733 PetscFunctionBeginUser;
1734
1735 PROFILE_FUNCTION_BEGIN;
1736
1737 lxs = xs; lxe = xe;
1738 lys = ys; lye = ye;
1739 lzs = zs; lze = ze;
1740
1741 if (xs == 0) lxs = xs + 1;
1742 if (ys == 0) lys = ys + 1;
1743 if (zs == 0) lzs = zs + 1;
1744
1745 if (xe == mx) lxe = xe - 1;
1746 if (ye == my) lye = ye - 1;
1747 if (ze == mz) lze = ze - 1;
1748
1749 DMDAVecGetArray(fda, user->lCsi, &csi);
1750 DMDAVecGetArray(fda, user->lEta, &eta);
1751 DMDAVecGetArray(fda, user->lZet, &zet);
1752 DMDAVecGetArray(da, user->lAj, &aj);
1753
1754 VecDuplicate(user->P, &Div);
1755 VecSet(Div, 0.);
1756 DMDAVecGetArray(da, Div, &div);
1757
1758 for (k = lzs; k < lze; k++) {
1759 for (j = lys; j < lye; j++) {
1760 for (i = lxs; i < lxe; i++) {
1761 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1762 eta[k][j][i].x - eta[k][j-1][i].x +
1763 zet[k][j][i].x - zet[k-1][j][i].x +
1764 csi[k][j][i].y - csi[k][j][i-1].y +
1765 eta[k][j][i].y - eta[k][j-1][i].y +
1766 zet[k][j][i].y - zet[k-1][j][i].y +
1767 csi[k][j][i].z - csi[k][j][i-1].z +
1768 eta[k][j][i].z - eta[k][j-1][i].z +
1769 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1770 div[k][j][i] = fabs(divergence);
1771 }
1772 }
1773 }
1774
1775 DMDAVecRestoreArray(da, Div, &div);
1776
1777 PetscInt MaxFlatIndex = -1;
1778 VecMax(Div, &MaxFlatIndex, &maxdiv);
1779 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %" PetscInt_FMT ".\n",maxdiv,MaxFlatIndex);
1780
1781 for (k=zs; k<ze; k++) {
1782 for (j=ys; j<ye; j++) {
1783 for (i=xs; i<xe; i++) {
1784 if (Gidx(i,j,k,user) == MaxFlatIndex) {
1785 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,(int)k,(int)j,(int)i);
1786 }
1787 }
1788 }
1789 }
1790
1791
1792 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1793 DMDAVecRestoreArray(fda, user->lEta, &eta);
1794 DMDAVecRestoreArray(fda, user->lZet, &zet);
1795 DMDAVecRestoreArray(da, user->lAj, &aj);
1796 VecDestroy(&Div);
1797
1798
1799 PROFILE_FUNCTION_END;
1800
1801 PetscFunctionReturn(0);
1802}
1803
1804#undef __FUNCT__
1805#define __FUNCT__ "ComputeMetricNorms"
1806/**
1807 * @brief Internal helper implementation: `ComputeMetricNorms()`.
1808 * @details Local to this translation unit.
1809 */
1810PetscErrorCode ComputeMetricNorms(UserCtx *user)
1811{
1812
1813 DMDALocalInfo info = user->info;
1814 PetscInt xs = info.xs, xe = info.xs + info.xm;
1815 PetscInt ys = info.ys, ye = info.ys + info.ym;
1816 PetscInt zs = info.zs, ze = info.zs + info.zm;
1817 PetscInt i, j, k;
1818
1819 PetscFunctionBeginUser;
1820
1821 PROFILE_FUNCTION_BEGIN;
1822
1823 PetscReal CsiMax, EtaMax, ZetMax;
1824 PetscReal ICsiMax, IEtaMax, IZetMax;
1825 PetscReal JCsiMax, JEtaMax, JZetMax;
1826 PetscReal KCsiMax, KEtaMax, KZetMax;
1827 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
1828
1829 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
1830 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
1831 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
1832 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
1833 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
1834
1835 // Max Values
1836 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
1837 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
1838 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
1839
1840 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
1841 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
1842 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
1843
1844 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
1845 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
1846 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
1847
1848 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
1849 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
1850 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
1851
1852 VecMax(user->lAj,&AjMaxArg,&AjMax);
1853 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1854 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1855 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1856
1857 VecMax(user->lAj,&AjMaxArg,&AjMax);
1858 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1859 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1860 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1861
1862 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
1863
1864 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
1865 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
1866 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
1867 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
1868 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
1869
1870 for (k=zs; k<ze; k++) {
1871 for (j=ys; j<ye; j++) {
1872 for (i=xs; i<xe; i++) {
1873 if (Gidx(i,j,k,user) == CsiMaxArg) {
1874 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
1875 }
1876 if (Gidx(i,j,k,user) == EtaMaxArg) {
1877 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
1878 }
1879 if (Gidx(i,j,k,user) == ZetMaxArg) {
1880 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
1881 }
1882 if (Gidx(i,j,k,user) == ICsiMaxArg) {
1883 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
1884 }
1885 if (Gidx(i,j,k,user) == IEtaMaxArg) {
1886 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
1887 }
1888 if (Gidx(i,j,k,user) == IZetMaxArg) {
1889 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
1890 }
1891 if (Gidx(i,j,k,user) == JCsiMaxArg) {
1892 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
1893 }
1894 if (Gidx(i,j,k,user) == JEtaMaxArg) {
1895 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
1896 }
1897 if (Gidx(i,j,k,user) == JZetMaxArg) {
1898 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
1899 }
1900 if (Gidx(i,j,k,user) == KCsiMaxArg) {
1901 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
1902 }
1903 if (Gidx(i,j,k,user) == KEtaMaxArg) {
1904 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
1905 }
1906 if (Gidx(i,j,k,user) == KZetMaxArg) {
1907 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
1908 }
1909 if (Gidx(i,j,k,user) == AjMaxArg) {
1910 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
1911 }
1912 if (Gidx(i,j,k,user) == IAjMaxArg) {
1913 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
1914 }
1915 if (Gidx(i,j,k,user) == JAjMaxArg) {
1916 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
1917 }
1918 if (Gidx(i,j,k,user) == KAjMaxArg) {
1919 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
1920 }
1921 }
1922 }
1923 }
1924
1925 /*
1926 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
1927 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
1928 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
1929 */
1930
1931 PROFILE_FUNCTION_END;
1932
1933 PetscFunctionReturn(0);
1934}
1935
1936#undef __FUNCT__
1937#define __FUNCT__ "ComputeAllGridMetrics"
1938/**
1939 * @brief Internal helper implementation: `CalculateAllGridMetrics()`.
1940 * @details Local to this translation unit.
1941 */
1942PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
1943{
1944 PetscErrorCode ierr;
1945 UserMG *usermg = &simCtx->usermg;
1946 MGCtx *mgctx = usermg->mgctx;
1947 PetscInt nblk = simCtx->block_number;
1948
1949 PetscFunctionBeginUser;
1950
1951 PROFILE_FUNCTION_BEGIN;
1952
1953 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
1954
1955 // Loop through all levels and all blocks
1956 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
1957 for (PetscInt bi = 0; bi < nblk; bi++) {
1958 UserCtx *user = &mgctx[level].user[bi];
1959 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
1960
1961 // Call the modern, modular helper functions for each UserCtx.
1962 // These functions are self-contained and operate on the data within the provided context.
1963 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
1964 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
1965 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
1966 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
1967 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
1968 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
1969 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
1970
1971 // Apply Periodic Boundary Condition Adjustments if necessary
1972 ierr = ApplyMetricsPeriodicBCs(user); CHKERRQ(ierr);
1973 // Diagnostics
1974 ierr = ComputeMetricNorms(user);
1975 if (level == usermg->mglevels - 1) {
1976 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
1977 }
1978 }
1979 }
1980
1981 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
1982
1983 PROFILE_FUNCTION_END;
1984
1985 PetscFunctionReturn(0);
1986}
1987/* ------------------------------------------------------------------------- */
1988/* 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:2439
SynchronizePeriodicFaceFields
PetscErrorCode SynchronizePeriodicFaceFields(UserCtx *user, char face_direction, PetscInt num_fields, const char *field_names[])
Synchronizes persistent fields belonging to one face family.
Definition Boundaries.c:1916
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Internal helper implementation: CalculateAllGridMetrics().
Definition Metric.c:1942
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:1810
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:1696
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:994
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:1308
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Internal helper implementation: ComputeMetricsDivergence().
Definition Metric.c:1718
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeFaceMetrics().
Definition Metric.c:638
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:850
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Internal helper implementation: ComputeKFaceMetrics().
Definition Metric.c:1503
ApplyPeriodicCorrectionsToJFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToJFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToJFaceCenter().
Definition Metric.c:607
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:1097
ApplyPeriodicCorrectionsToKFaceCenter
PetscErrorCode ApplyPeriodicCorrectionsToKFaceCenter(UserCtx *user)
Internal helper implementation: ApplyPeriodicCorrectionsToKFaceCenter().
Definition Metric.c:623
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:827
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:818
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1755
UserCtx::lCent
Vec lCent
Definition variables.h:927
PERIODIC
@ PERIODIC
Definition variables.h:290
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:927
UserCtx::JCsi
Vec JCsi
Definition variables.h:931
UserCtx::KAj
Vec KAj
Definition variables.h:932
MGCtx::user
UserCtx * user
Definition variables.h:569
UserCtx::JEta
Vec JEta
Definition variables.h:931
UserCtx::Zet
Vec Zet
Definition variables.h:927
SimCtx::rank
PetscMPIInt rank
Definition variables.h:687
UserCtx::cgrid
PetscInt cgrid
Definition variables.h:891
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:896
SimCtx::block_number
PetscInt block_number
Definition variables.h:768
UserCtx::lIEta
Vec lIEta
Definition variables.h:930
UserCtx::lIZet
Vec lIZet
Definition variables.h:930
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:879
UserCtx::IZet
Vec IZet
Definition variables.h:930
UserCtx::Centz
Vec Centz
Definition variables.h:928
UserCtx::IEta
Vec IEta
Definition variables.h:930
UserCtx::lZet
Vec lZet
Definition variables.h:927
SimCtx::usermg
UserMG usermg
Definition variables.h:821
UserCtx::Csi
Vec Csi
Definition variables.h:927
UserCtx::lIAj
Vec lIAj
Definition variables.h:930
UserCtx::_this
PetscInt _this
Definition variables.h:889
UserCtx::lKEta
Vec lKEta
Definition variables.h:932
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:931
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::JZet
Vec JZet
Definition variables.h:931
UserCtx::Centx
Vec Centx
Definition variables.h:928
UserCtx::lKZet
Vec lKZet
Definition variables.h:932
UserCtx::Eta
Vec Eta
Definition variables.h:927
UserCtx::lJEta
Vec lJEta
Definition variables.h:931
UserCtx::lCsi
Vec lCsi
Definition variables.h:927
UserCtx::lGridSpace
Vec lGridSpace
Definition variables.h:927
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:944
UserCtx::ICsi
Vec ICsi
Definition variables.h:930
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:932
UserCtx::lCenty
Vec lCenty
Definition variables.h:929
UserMG::mglevels
PetscInt mglevels
Definition variables.h:576
UserCtx::lJZet
Vec lJZet
Definition variables.h:931
UserCtx::IAj
Vec IAj
Definition variables.h:930
UserCtx::JAj
Vec JAj
Definition variables.h:931
UserCtx::KEta
Vec KEta
Definition variables.h:932
UserCtx::lCentx
Vec lCentx
Definition variables.h:929
UserCtx::lAj
Vec lAj
Definition variables.h:927
UserCtx::lICsi
Vec lICsi
Definition variables.h:930
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:889
UserCtx::info
DMDALocalInfo info
Definition variables.h:883
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:927
UserCtx::KZet
Vec KZet
Definition variables.h:932
UserCtx::Cent
Vec Cent
Definition variables.h:927
UserCtx::KCsi
Vec KCsi
Definition variables.h:932
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:579
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:366
UserCtx::Centy
Vec Centy
Definition variables.h:928
UserCtx::lCentz
Vec lCentz
Definition variables.h:929
UserCtx::lJAj
Vec lJAj
Definition variables.h:931
UserCtx::lKAj
Vec lKAj
Definition variables.h:932
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:260
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:262
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:261
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:262
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:260
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:261
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
MGCtx
Context for Multigrid operations.
Definition variables.h:568
SimCtx
The master context for the entire simulation.
Definition variables.h:684
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:876
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:575
Generated by doxygen 1.9.8