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/**
18 * @brief Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
19 *
20 * Vertices are returned in the standard trilinear ordering: bit 0 → x-corner,
21 * bit 1 → y-corner, bit 2 → z-corner. (000 = origin of the cell, 111 = far
22 * corner.)
23 */
24PetscErrorCode MetricGetCellVertices(UserCtx *user,
25 const Cmpnts ***X, /* coord array */
26 PetscInt i,PetscInt j,PetscInt k,
27 Cmpnts V[8])
28{
29 PetscFunctionBeginUser;
30 for (PetscInt c = 0; c < 8; ++c) {
31 PetscInt ii = i + ((c & 1) ? 1 : 0);
32 PetscInt jj = j + ((c & 2) ? 1 : 0);
33 PetscInt kk = k + ((c & 4) ? 1 : 0);
34 LOG_LOOP_ALLOW(GLOBAL, LOG_DEBUG,i+j+k,10," ii: %d,jj:%d,kk:%d - Retrieved.\n",ii,jj,kk);
35 V[c] = X[kk][jj][ii];
36 }
37 PetscFunctionReturn(0);
38}
39
40/* ------------------------------------------------------------------------- */
41/**
42 * @brief Map logical coordinates to physical space using trilinear basis.
43 *
44 * @param[in] V Array of the eight vertex coordinates (MetricGetCellVertices).
45 * @param[in] xi,eta,zta Logical coordinates in [0,1].
46 * @param[out] Xp Physical coordinate.
47 */
48static inline void TrilinearBlend(const Cmpnts V[8],
49 PetscReal xi,PetscReal eta,PetscReal zta,
50 Cmpnts *Xp)
51{
52 PetscReal x=0,y=0,z=0;
53 for (PetscInt c=0;c<8;++c) {
54 PetscReal N = ((c&1)?xi : 1.0-xi ) *
55 ((c&2)?eta: 1.0-eta) *
56 ((c&4)?zta: 1.0-zta);
57 x += N * V[c].x;
58 y += N * V[c].y;
59 z += N * V[c].z;
60 }
61 Xp->x = x; Xp->y = y; Xp->z = z;
62}
63
64/* ------------------------------------------------------------------------- */
65/**
66 * @brief Public wrapper: map (cell index, ξ,η,ζ) to (x,y,z).
67 */
68PetscErrorCode MetricLogicalToPhysical(UserCtx *user,
69 const Cmpnts ***X,
70 PetscInt i,PetscInt j,PetscInt k,
71 PetscReal xi,PetscReal eta,PetscReal zta,
72 Cmpnts *Xp)
73{
74 PetscErrorCode ierr;
75 Cmpnts V[8];
76 PetscFunctionBeginUser;
77 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
78 TrilinearBlend(V,xi,eta,zta,Xp);
79 PetscFunctionReturn(0);
80}
81
82/* ------------------------------------------------------------------------- */
83/**
84 * @brief Compute Jacobian matrix and its determinant at (xi,eta,zta).
85 *
86 * J = [ x_ξ x_η x_ζ ]
87 * [ y_ξ y_η y_ζ ]
88 * [ z_ξ z_η z_ζ ]
89 *
90 * This is handy for converting physical velocities (u,v,w) into contravariant
91 * components and for volume weighting.
92 */
93PetscErrorCode MetricJacobian(UserCtx *user,
94 const Cmpnts ***X,
95 PetscInt i,PetscInt j,PetscInt k,
96 PetscReal xi,PetscReal eta,PetscReal zta,
97 PetscReal J[3][3], PetscReal *detJ)
98{
99 PetscErrorCode ierr;
100 Cmpnts V[8];
101 PetscFunctionBeginUser;
102 ierr = MetricGetCellVertices(user,X,i,j,k,V); CHKERRQ(ierr);
103
104 /* derivatives of trilinear shape functions */
105 PetscReal dN_dXi[8], dN_dEta[8], dN_dZta[8];
106 for (PetscInt c=0;c<8;++c) {
107 PetscReal sx = (c & 1) ? 1.0 : -1.0;
108 PetscReal sy = (c & 2) ? 1.0 : -1.0;
109 PetscReal sz = (c & 4) ? 1.0 : -1.0;
110 dN_dXi [c] = 0.125 * sx * ( (c&2?eta:1-eta) ) * ( (c&4?zta:1-zta) );
111 dN_dEta[c] = 0.125 * sy * ( (c&1?xi :1-xi ) ) * ( (c&4?zta:1-zta) );
112 dN_dZta[c] = 0.125 * sz * ( (c&1?xi :1-xi ) ) * ( (c&2?eta:1-eta) );
113 }
114
115 /* assemble Jacobian */
116 PetscReal x_xi=0,y_xi=0,z_xi=0,
117 x_eta=0,y_eta=0,z_eta=0,
118 x_zta=0,y_zta=0,z_zta=0;
119 for (PetscInt c=0;c<8;++c) {
120 x_xi += dN_dXi [c]*V[c].x; y_xi += dN_dXi [c]*V[c].y; z_xi += dN_dXi [c]*V[c].z;
121 x_eta += dN_dEta[c]*V[c].x; y_eta += dN_dEta[c]*V[c].y; z_eta += dN_dEta[c]*V[c].z;
122 x_zta += dN_dZta[c]*V[c].x; y_zta += dN_dZta[c]*V[c].y; z_zta += dN_dZta[c]*V[c].z;
123 }
124
125 J[0][0]=x_xi; J[0][1]=x_eta; J[0][2]=x_zta;
126 J[1][0]=y_xi; J[1][1]=y_eta; J[1][2]=y_zta;
127 J[2][0]=z_xi; J[2][1]=z_eta; J[2][2]=z_zta;
128
129 if (detJ) {
130 *detJ = x_xi*(y_eta*z_zta - y_zta*z_eta)
131 - x_eta*(y_xi*z_zta - y_zta*z_xi)
132 + x_zta*(y_xi*z_eta - y_eta*z_xi);
133 }
134 PetscFunctionReturn(0);
135}
136
137/* ------------------------------------------------------------------------- */
138/**
139 * @brief Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).
140 */
141PetscErrorCode MetricVelocityContravariant(const PetscReal J[3][3], PetscReal detJ,
142 const PetscReal u[3], PetscReal uc[3])
143{
144 PetscFunctionBeginUser;
145 /* contravariant basis vectors (row of adjugate(J)) divided by detJ */
146 PetscReal gxi[3] = { J[1][1]*J[2][2]-J[1][2]*J[2][1],
147 -J[0][1]*J[2][2]+J[0][2]*J[2][1],
148 J[0][1]*J[1][2]-J[0][2]*J[1][1] };
149 PetscReal geta[3] = { -J[1][0]*J[2][2]+J[1][2]*J[2][0],
150 J[0][0]*J[2][2]-J[0][2]*J[2][0],
151 -J[0][0]*J[1][2]+J[0][2]*J[1][0] };
152 PetscReal gzta[3] = { J[1][0]*J[2][1]-J[1][1]*J[2][0],
153 -J[0][0]*J[2][1]+J[0][1]*J[2][0],
154 J[0][0]*J[1][1]-J[0][1]*J[1][0] };
155
156 PetscReal invDet = 1.0 / detJ;
157 for (int d=0; d<3; ++d) { gxi[d] *= invDet; geta[d] *= invDet; gzta[d] *= invDet; }
158
159 uc[0] = gxi [0]*u[0] + gxi [1]*u[1] + gxi [2]*u[2];
160 uc[1] = geta[0]*u[0] + geta[1]*u[1] + geta[2]*u[2];
161 uc[2] = gzta[0]*u[0] + gzta[1]*u[1] + gzta[2]*u[2];
162 PetscFunctionReturn(0);
163}
164
165
166/* -------------------------------------------------------------------------- */
167/**
168 * @brief Ensure a **right-handed** metric basis (`Csi`, `Eta`, `Zet`) and a
169 * **positive Jacobian** (`Aj`) over the whole domain.
170 *
171 * The metric-generation kernels are completely algebraic, so they will happily
172 * deliver a *left-handed* basis if the mesh file enumerates nodes in the
173 * opposite ζ-direction.
174 * This routine makes the orientation explicit and—if needed—repairs it
175 * **once per run**:
176 *
177 * | Step | Action |
178 * |------|--------|
179 * | 1 | Compute global `Aj_min`, `Aj_max`. |
180 * | 2 | **Mixed signs** (`Aj_min < 0 && Aj_max > 0`) &rarr; abort: the mesh is topologically inconsistent. |
181 * | 3 | **All negative** (`Aj_max < 0`) &rarr; flip <br>`Csi`, `Eta`, `Zet`, `Aj` & update local ghosts. |
182 * | 4 | Store `user->orientation = ±1` so BC / IC routines can apply sign-aware logic if they care about inlet direction. |
183 *
184 * @param[in,out] user Fully initialised #UserCtx that already contains
185 * `Csi`, `Eta`, `Zet`, `Aj`, their **local** ghosts, and
186 * valid distributed DMs.
187 *
188 * @return `0` on success or a PETSc error code on failure.
189 *
190 * @note Call **immediately after** `ComputeCellCenteredJacobianInverse()` and
191 * before any routine that differentiates or applies BCs.
192 *
193 * @author vishal kandala
194 */
195PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
196{
197 PetscErrorCode ierr;
198 PetscReal aj_min, aj_max;
199 PetscMPIInt rank;
200
201 PetscFunctionBeginUser;
202
203 /* ---------------- step 1: global extrema of Aj ---------------- */
204 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
205
206 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr); /* already global */
207 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
208
209 LOG_ALLOW(GLOBAL, LOG_INFO,
210 "[orientation] Global Aj range: [%.3e , %.3e]\n",
211 (double)aj_min, (double)aj_max);
212
213 /* ---------------- step 2: detect malformed mesh ---------------- */
214 if (aj_min < 0.0 && aj_max > 0.0)
215 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
216 "Mixed Jacobian signs detected – grid is topologically inconsistent.");
217
218 /* Default: grid is right-handed unless proven otherwise */
219 PetscInt orientation = +1;
220
221 /* ---------------- step 3: repair left-handed mesh -------------- */
222 if (aj_max < 0.0) { /* entire domain has Aj < 0 */
223 orientation = -1;
224
225 if (!rank)
226 LOG_ALLOW(LOCAL, LOG_INFO,
227 "[orientation] Detected left-handed grid – flipping metric vectors\n");
228
229 /* Flip sign of *all* metric vectors and Aj */
230 ierr = VecScale(user->Csi, -1.0); CHKERRQ(ierr);
231 ierr = VecScale(user->Eta, -1.0); CHKERRQ(ierr);
232 ierr = VecScale(user->Zet, -1.0); CHKERRQ(ierr);
233 ierr = VecScale(user->Aj , -1.0); CHKERRQ(ierr);
234
235 /* Local ghost regions now stale – refresh */
236 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
237 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
238 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
239 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
240
241 /* Sanity print: Aj must be > 0 now */
242 ierr = VecMin(user->Aj, NULL, &aj_min); CHKERRQ(ierr);
243 ierr = VecMax(user->Aj, NULL, &aj_max); CHKERRQ(ierr);
244
245 if (aj_min <= 0.0)
246 SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_USER,
247 "Failed to flip grid orientation – Aj still non-positive.");
248 else if (aj_min && aj_max > 0.0)
249 orientation = +1;
250 }
251
252 /* ---------------- step 4: store result in UserCtx -------------- */
253 user->GridOrientation = orientation;
254
255 if (!rank)
256 LOG_ALLOW(LOCAL, LOG_INFO,
257 "[orientation] Grid confirmed %s-handed after flip (orientation=%+d)\n",
258 (orientation>0) ? "right" : "left", orientation);
259
260 PetscFunctionReturn(0);
261}
262
263
264/**
265 * @brief Computes the primary face metric components (Csi, Eta, Zet), including
266 * boundary extrapolation, and stores them in the corresponding global Vec
267 * members of the UserCtx structure (user->Csi, user->Eta, user->Zet).
268 *
269 * This is a self-contained routine that performs the following steps:
270 * 1. Obtains local ghosted nodal coordinates using DMGetCoordinatesLocal.
271 * 2. Calculates metrics for INTERIOR faces where finite difference stencils are valid.
272 * 3. EXTRAPOLATES metrics for faces on the physical domain boundaries by copying
273 * from the nearest computed interior face.
274 * 4. Assembles the global `user->Csi`, `user->Eta`, `user->Zet` Vecs.
275 * 5. Updates the local ghosted `user->lCsi`, `user->lEta`, `user->lZet` Vecs.
276 *
277 * @param[in,out] user Pointer to the UserCtx structure.
278 *
279 * @return PetscErrorCode 0 on success.
280 *
281 * @note
282 * - This function is a complete "compute and make ready" unit for Csi, Eta, and Zet.
283 * - It's recommended to call `VecZeroEntries` on user->Csi, Eta, Zet before this
284 * if they might contain old data.
285 */
286PetscErrorCode ComputeFaceMetrics(UserCtx *user)
287{
288 PetscErrorCode ierr;
289 DMDALocalInfo info;
290 Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
291 Cmpnts ***nodal_coords_arr;
292 Vec localCoords_from_dm;
293
294 PetscFunctionBeginUser;
295 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation and update for Csi, Eta, Zet.\n");
296
297 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
298
299 // --- 1. Get Nodal Physical Coordinates (Local Ghosted Array directly) ---
300 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
301 if (!localCoords_from_dm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocal failed to return a coordinate vector. \n");
302 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
303
304 // --- 2. Get arrays for output global Vecs from UserCtx ---
305 ierr = DMDAVecGetArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
306 ierr = DMDAVecGetArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
307 ierr = DMDAVecGetArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
308
309 // Define owned node ranges (global indices)
310 PetscInt xs = info.xs, xe = info.xs + info.xm;
311 PetscInt ys = info.ys, ye = info.ys + info.ym;
312 PetscInt zs = info.zs, ze = info.zs + info.zm;
313
314 // Global domain dimensions (total number of nodes)
315 PetscInt mx = info.mx;
316 PetscInt my = info.my;
317 PetscInt mz = info.mz;
318
319 // --- 3. Calculate Csi, Eta, Zet for INTERIOR Stencils ---
320 // Start loops from 1 if at global boundary 0 to ensure k_node-1 etc. are valid.
321 PetscInt k_loop_start = (zs == 0) ? zs + 1 : zs;
322 PetscInt j_loop_start = (ys == 0) ? ys + 1 : ys;
323 PetscInt i_loop_start = (xs == 0) ? xs + 1 : xs;
324
325 // Calculate Csi
326 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
327 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
328 for (PetscInt i_node = xs; i_node < xe; ++i_node) {
329
330 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);
331 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);
332 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);
333 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);
334 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);
335 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);
336
337 csi_arr[k_node][j_node][i_node].x = dy_deta * dz_dzeta - dz_deta * dy_dzeta;
338 csi_arr[k_node][j_node][i_node].y = dz_deta * dx_dzeta - dx_deta * dz_dzeta;
339 csi_arr[k_node][j_node][i_node].z = dx_deta * dy_dzeta - dy_deta * dx_dzeta;
340 }
341 }
342 }
343
344 // Calculate Eta
345 for (PetscInt k_node = k_loop_start; k_node < ze; ++k_node) {
346 for (PetscInt j_node = ys; j_node < ye; ++j_node) {
347 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
348
349 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);
350 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);
351 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);
352 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);
353 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);
354 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);
355
356 eta_arr[k_node][j_node][i_node].x = dy_dzeta * dz_dxi - dz_dzeta * dy_dxi;
357 eta_arr[k_node][j_node][i_node].y = dz_dzeta * dx_dxi - dx_dzeta * dz_dxi;
358 eta_arr[k_node][j_node][i_node].z = dx_dzeta * dy_dxi - dy_dzeta * dx_dxi;
359 }
360 }
361 }
362
363 // Calculate Zet
364 for (PetscInt k_node = zs; k_node < ze; ++k_node) {
365 for (PetscInt j_node = j_loop_start; j_node < ye; ++j_node) {
366 for (PetscInt i_node = i_loop_start; i_node < xe; ++i_node) {
367
368 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);
369 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);
370 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);
371 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);
372 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);
373 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);
374
375 zet_arr[k_node][j_node][i_node].x = dy_dxi * dz_deta - dz_dxi * dy_deta;
376 zet_arr[k_node][j_node][i_node].y = dz_dxi * dx_deta - dx_dxi * dz_deta;
377 zet_arr[k_node][j_node][i_node].z = dx_dxi * dy_deta - dy_dxi * dx_deta;
378 }
379 }
380 }
381
382 // --- 4. Boundary Extrapolation ---
383 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Csi, Eta, Zet.\n");
384 PetscInt i_bnd, j_bnd, k_bnd;
385
386 if (xs == 0) { // If this rank owns the global i=0 boundary
387 i_bnd = 0;
388 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
389 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
390 if (i_bnd + 1 < mx) {
391 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd+1];
392 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd+1];
393 }
394 }
395 }
396 }
397 if (xe == mx) { // If this rank owns the global i=mx-1 boundary
398 i_bnd = mx - 1;
399 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
400 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
401 if (i_bnd - 1 >= 0) {
402 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd][j_bnd][i_bnd-1];
403 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd][i_bnd-1];
404 }
405 }
406 }
407 }
408 if (ys == 0) {
409 j_bnd = 0;
410 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
411 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
412 if (j_bnd + 1 < my) {
413 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd+1][i_bnd];
414 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd+1][i_bnd];
415 }
416 }
417 }
418 }
419 if (ye == my) {
420 j_bnd = my - 1;
421 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
422 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
423 if (j_bnd - 1 >= 0) {
424 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd][j_bnd-1][i_bnd];
425 zet_arr[k_bnd][j_bnd][i_bnd] = zet_arr[k_bnd][j_bnd-1][i_bnd];
426 }
427 }
428 }
429 }
430 if (zs == 0) {
431 k_bnd = 0;
432 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
433 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
434 if (k_bnd + 1 < mz) {
435 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd+1][j_bnd][i_bnd];
436 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd+1][j_bnd][i_bnd];
437 }
438 }
439 }
440 }
441 if (ze == mz) {
442 k_bnd = mz - 1;
443 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
444 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
445 if (k_bnd - 1 >= 0) {
446 csi_arr[k_bnd][j_bnd][i_bnd] = csi_arr[k_bnd-1][j_bnd][i_bnd];
447 eta_arr[k_bnd][j_bnd][i_bnd] = eta_arr[k_bnd-1][j_bnd][i_bnd];
448 }
449 }
450 }
451 }
452
453 if (info.xs==0 && info.ys==0 && info.zs==0) {
454 PetscReal dot = zet_arr[0][0][0].z; /* dot with global +z */
455 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Zet(k=0)·ez = %.3f (should be >0 for right-handed grid)\n", dot);
456 }
457
458 // --- 5. Restore all arrays ---
459 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
460 ierr = DMDAVecRestoreArray(user->fda, user->Csi, &csi_arr); CHKERRQ(ierr);
461 ierr = DMDAVecRestoreArray(user->fda, user->Eta, &eta_arr); CHKERRQ(ierr);
462 ierr = DMDAVecRestoreArray(user->fda, user->Zet, &zet_arr); CHKERRQ(ierr);
463
464 // --- 6. Assemble Global Vectors ---
465 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Csi, Eta, Zet.\n");
466 ierr = VecAssemblyBegin(user->Csi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Csi); CHKERRQ(ierr);
467 ierr = VecAssemblyBegin(user->Eta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Eta); CHKERRQ(ierr);
468 ierr = VecAssemblyBegin(user->Zet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Zet); CHKERRQ(ierr);
469
470 // --- 7. Update Local Ghosted Versions ---
471 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lCsi, lEta, lZet.\n");
472 ierr = UpdateLocalGhosts(user, "Csi"); CHKERRQ(ierr);
473 ierr = UpdateLocalGhosts(user, "Eta"); CHKERRQ(ierr);
474 ierr = UpdateLocalGhosts(user, "Zet"); CHKERRQ(ierr);
475
476 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Csi, Eta, Zet.\n");
477 PetscFunctionReturn(0);
478}
479
480
481
482/**
483 * @brief Calculates the cell-centered inverse Jacobian determinant (1/J), including
484 * boundary extrapolation, stores it in `user->Aj`, assembles `user->Aj`, and
485 * updates `user->lAj`.
486 *
487 * @param[in,out] user Pointer to the UserCtx structure.
488 *
489 * @return PetscErrorCode 0 on success.
490 */
491PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
492{
493 PetscErrorCode ierr;
494 DMDALocalInfo info;
495 PetscScalar ***aj_arr;
496 Cmpnts ***nodal_coords_arr;
497 Vec localCoords_from_dm;
498
499 PetscFunctionBeginUser;
500 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting calculation, extrapolation, and update for Aj.\n");
501
502 // --- 1. Get Nodal Coordinates and Output Array ---
503 ierr = DMGetCoordinatesLocal(user->da, &localCoords_from_dm); CHKERRQ(ierr);
504 ierr = DMDAVecGetArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
505 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
506 ierr = DMDAVecGetArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
507
508 // Define owned node ranges (global indices)
509 PetscInt xs = info.xs, xe = info.xs + info.xm;
510 PetscInt ys = info.ys, ye = info.ys + info.ym;
511 PetscInt zs = info.zs, ze = info.zs + info.zm;
512
513 // Global domain dimensions (total number of nodes)
514 PetscInt mx = info.mx;
515 PetscInt my = info.my;
516 PetscInt mz = info.mz;
517
518 // --- 2. Calculate Aj for INTERIOR Stencils ---
519
520 PetscInt k_start_node = (zs == 0) ? zs + 1 : zs;
521 PetscInt j_start_node = (ys == 0) ? ys + 1 : ys;
522 PetscInt i_start_node = (xs == 0) ? xs + 1 : xs;
523
524 PetscInt k_end_node = (ze == mz) ? ze - 1 : ze;
525 PetscInt j_end_node = (ye == my) ? ye - 1 : ye;
526 PetscInt i_end_node = (xe == mx) ? xe - 1 : xe;
527
528 for (PetscInt k_node = k_start_node; k_node < k_end_node; ++k_node) {
529 for (PetscInt j_node = j_start_node; j_node < j_end_node; ++j_node) {
530 for (PetscInt i_node = i_start_node; i_node < i_end_node; ++i_node) {
531
532 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) );
533
534 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) );
535
536 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) );
537
538 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) );
539
540 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) );
541
542 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) );
543
544 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) );
545
546 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) );
547
548 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) );
549
550 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);
551 if (PetscAbsReal(jacobian_det) < 1.0e-18) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Jacobian is near zero..."); }
552 aj_arr[k_node][j_node][i_node] = 1.0 / jacobian_det;
553 }
554 }
555 }
556
557 // --- 4. Boundary Extrapolation for Aj ---
558 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Extrapolating boundary values for Aj. \n");
559 PetscInt i_bnd, j_bnd, k_bnd;
560
561 if (xs == 0) {
562 i_bnd = 0;
563 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
564 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
565 if (i_bnd + 1 < mx) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd+1];
566 }
567 }
568 }
569 if (xe == mx) {
570 i_bnd = mx - 1;
571 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
572 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
573 if (i_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd][i_bnd-1];
574 }
575 }
576 }
577 // (Similar extrapolation blocks for Y and Z boundaries for aj_arr)
578 if (ys == 0) {
579 j_bnd = 0;
580 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
581 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
582 if (j_bnd + 1 < my) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd+1][i_bnd];
583 }
584 }
585 }
586 if (ye == my) {
587 j_bnd = my - 1;
588 for (k_bnd = zs; k_bnd < ze; ++k_bnd) {
589 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
590 if (j_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd][j_bnd-1][i_bnd];
591 }
592 }
593 }
594 if (zs == 0) {
595 k_bnd = 0;
596 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
597 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
598 if (k_bnd + 1 < mz) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd+1][j_bnd][i_bnd];
599 }
600 }
601 }
602 if (ze == mz) {
603 k_bnd = mz - 1;
604 for (j_bnd = ys; j_bnd < ye; ++j_bnd) {
605 for (i_bnd = xs; i_bnd < xe; ++i_bnd) {
606 if (k_bnd - 1 >= 0) aj_arr[k_bnd][j_bnd][i_bnd] = aj_arr[k_bnd-1][j_bnd][i_bnd];
607 }
608 }
609 }
610
611 // --- 5. Restore arrays ---
612 ierr = DMDAVecRestoreArrayRead(user->fda, localCoords_from_dm, &nodal_coords_arr); CHKERRQ(ierr);
613 ierr = DMDAVecRestoreArray(user->da, user->Aj, &aj_arr); CHKERRQ(ierr);
614
615 // --- 6. Assemble Global Vector ---
616 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global Aj.\n");
617 ierr = VecAssemblyBegin(user->Aj); CHKERRQ(ierr);
618 ierr = VecAssemblyEnd(user->Aj); CHKERRQ(ierr);
619
620 // --- 7. Update Local Ghosted Version ---
621 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating local lAj.\n");
622 ierr = UpdateLocalGhosts(user, "Aj"); CHKERRQ(ierr);
623
624 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed calculation, extrapolation, and update for Aj.\n");
625 PetscFunctionReturn(0);
626}
627
628
629#undef __FUNCT__
630#define __FUNCT__ "ComputeCellCentersAndSpacing"
631/**
632 * @brief Computes the physical location of cell centers and the spacing between them.
633 *
634 * This function calculates two key geometric properties from the nodal coordinates:
635 * 1. `Cent`: A vector field storing the (x,y,z) coordinates of the center of each grid cell.
636 * 2. `GridSpace`: A vector field storing the physical distance between adjacent
637 * cell centers in the i, j, and k computational directions.
638 *
639 * It is a direct adaptation of the corresponding logic from the legacy `FormMetrics`.
640 *
641 * @param user The UserCtx for a specific grid level. The function populates `user->Cent` and `user->GridSpace`.
642 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
643 */
644PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
645{
646 PetscErrorCode ierr;
647 DMDALocalInfo info;
648 Vec lCoords;
649 const Cmpnts ***coor;
650 Cmpnts ***cent, ***gs;
651 PetscReal xcp, ycp, zcp, xcm, ycm, zcm;
652
653 PetscFunctionBeginUser;
654 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);
655
656 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
657 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
658 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
659
660 ierr = DMDAVecGetArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
661 ierr = DMDAVecGetArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
662
663 // Loop over the interior OWNED cells (stencil requires i-1, j-1, k-1)
664 for (PetscInt k=info.zs+1; k<info.zs+info.zm; k++) {
665 for (PetscInt j=info.ys+1; j<info.ys+info.ym; j++) {
666 for (PetscInt i=info.xs+1; i<info.xs+info.xm; i++) {
667 // Calculate cell center as the average of its 8 corner nodes
668 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);
669 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);
670 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);
671
672 // Calculate Grid Spacing in i-direction (distance between i-face centers)
673 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);
674 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);
675 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);
676 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);
677 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);
678 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);
679 gs[k][j][i].x = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
680
681 // Calculate Grid Spacing in j-direction (distance between j-face centers)
682 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);
683 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);
684 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);
685 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);
686 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);
687 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);
688 gs[k][j][i].y = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
689
690 // Calculate Grid Spacing in k-direction (distance between k-face centers)
691 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);
692 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);
693 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);
694 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);
695 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);
696 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);
697 gs[k][j][i].z = PetscSqrtReal(PetscSqr(xcp-xcm) + PetscSqr(ycp-ycm) + PetscSqr(zcp-zcm));
698 }
699 }
700 }
701
702 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
703 ierr = DMDAVecRestoreArray(user->fda, user->Cent, &cent); CHKERRQ(ierr);
704 ierr = DMDAVecRestoreArray(user->fda, user->GridSpace, &gs); CHKERRQ(ierr);
705
706 // Assemble and update ghost regions for the new data
707 ierr = VecAssemblyBegin(user->Cent); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->Cent); CHKERRQ(ierr);
708 ierr = VecAssemblyBegin(user->GridSpace); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->GridSpace); CHKERRQ(ierr);
709 ierr = UpdateLocalGhosts(user, "Cent"); CHKERRQ(ierr);
710 ierr = UpdateLocalGhosts(user, "GridSpace"); CHKERRQ(ierr);
711
712 PetscFunctionReturn(0);
713}
714
715#undef __FUNCT__
716#define __FUNCT__ "ComputeIFaceMetrics"
717/**
718 * @brief Computes metrics centered on constant-i faces (i-faces).
719 *
720 * This function calculates the metric terms (`ICsi`, `IEta`, `IZet`) and the
721 * inverse Jacobian (`IAj`) located at the geometric center of each constant-i
722 * face. This is a critical step for staggered-grid finite difference schemes.
723 *
724 * The process is a direct and faithful refactoring of the corresponding logic
725 * from the legacy `FormMetrics` function:
726 * 1. It first calculates the physical (x,y,z) coordinates of the center of
727 * each i-face and stores them in the `user->Centx` vector.
728 * 2. It then uses a boundary-aware, second-order finite difference stencil on
729 * the `Centx` field to compute the derivatives (e.g., d(x)/d(csi)).
730 * - Central differences are used in the grid interior.
731 * - One-sided differences are used at the physical domain boundaries.
732 * 3. Finally, these derivatives are used to compute the final metric terms and
733 * the inverse Jacobian, which are stored in their respective `Vec` objects.
734 *
735 * @param user The UserCtx for a specific grid level. This function populates
736 * the `user->ICsi`, `user->IEta`, `user->IZet`, and `user->IAj` vectors.
737 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
738 */
739PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
740{
741 PetscErrorCode ierr;
742 DMDALocalInfo info;
743 Vec lCoords;
744 const Cmpnts ***coor;
745 Cmpnts ***centx; //***gs;
746 const Cmpnts ***centx_const;
747 Cmpnts ***icsi, ***ieta, ***izet;
748 PetscScalar ***iaj;
749 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
750
751 PetscFunctionBeginUser;
752 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing i-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
753
754 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
755 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
756 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
757 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
758 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
759 PetscInt gys = info.gys, gye = info.gys + info.gym;
760 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
761
762 PetscInt lxs = xs; PetscInt lxe = xe;
763 PetscInt lys = ys; PetscInt lye = ye;
764 PetscInt lzs = zs; PetscInt lze = ze;
765
766 if (xs==0) lxs = xs+1;
767 if (ys==0) lys = ys+1;
768 if (zs==0) lzs = zs+1;
769
770 if (xe==mx) lxe=xe-1;
771 if (ye==my) lye=ye-1;
772 if (ze==mz) lze=ze-1;
773
774 // --- Part 1: Calculate the location of i-face centers (Centx) ---
775 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
776 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
777 ierr = DMDAVecGetArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
778 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
779
780 // Loop over the ghosted region to calculate all local face centers
781 for (PetscInt k = gzs + 1; k < gze; k++) {
782 for (PetscInt j = gys + 1; j < gye; j++) {
783 for (PetscInt i = gxs; i < gxe; i++) {
784 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);
785 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);
786 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);
787 }
788 }
789 }
790
791 /*
792 if(xs==0){
793 for(PetscInt k=gzs+1;k < gze; k++){
794 for(PetscInt j=gys+1;j < gye; j++){
795 PetscInt i=0;
796 centx[k][j][i-1].x=centx[k][j][i].x-gs[k][j][i-2].x;
797 centx[k][j][i-1].y=centx[k][j][i].y;
798 centx[k][j][i-1].z=centx[k][j][i].z;
799 }
800 }
801 }
802 if (xe==mx){
803 for(PetscInt k=gzs+1; k<gze; k++) {
804 for (PetscInt j=gys+1; j<gye;j++) {
805 PetscInt i=mx-1;
806 centx[k][j][i].x=centx[k][j][i-1].x+gs[k][j][i+2].x;
807 centx[k][j][i].y=centx[k][j][i-1].y;
808 centx[k][j][i].z=centx[k][j][i-1].z;
809 }
810 }
811 }
812 */
813
814 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
815 ierr = DMDAVecRestoreArray(user->fda, user->Centx, &centx); CHKERRQ(ierr);
816 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
817
818 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: i-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
819
820 // --- Part 2: Calculate metrics using face-centered coordinates ---
821 ierr = DMDAVecGetArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
822 ierr = DMDAVecGetArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
823 ierr = DMDAVecGetArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
824 ierr = DMDAVecGetArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
825 ierr = DMDAVecGetArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
826
827 // Loop over the OWNED region where we will store the final metrics
828 for (PetscInt k=lzs; k<lze; k++) {
829 for (PetscInt j=lys; j<lye; j++) {
830 for (PetscInt i=xs; i<lxe; i++) {
831
832 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
833 if (i == 0) { // Forward difference at the domain's min-i boundary
834 dxdc = centx_const[k][j][i+1].x - centx_const[k][j][i].x;
835 dydc = centx_const[k][j][i+1].y - centx_const[k][j][i].y;
836 dzdc = centx_const[k][j][i+1].z - centx_const[k][j][i].z;
837 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
838 dxdc = centx_const[k][j][i].x - centx_const[k][j][i-1].x;
839 dydc = centx_const[k][j][i].y - centx_const[k][j][i-1].y;
840 dzdc = centx_const[k][j][i].z - centx_const[k][j][i-1].z;
841 } else { // Central difference in the interior
842 dxdc = 0.5 * (centx_const[k][j][i+1].x - centx_const[k][j][i-1].x);
843 dydc = 0.5 * (centx_const[k][j][i+1].y - centx_const[k][j][i-1].y);
844 dzdc = 0.5 * (centx_const[k][j][i+1].z - centx_const[k][j][i-1].z);
845 }
846
847 // --- Stencil Logic for d/deta (derivative in j-direction) ---
848 if (j == 1) { // Forward difference
849 dxde = centx_const[k][j+1][i].x - centx_const[k][j][i].x;
850 dyde = centx_const[k][j+1][i].y - centx_const[k][j][i].y;
851 dzde = centx_const[k][j+1][i].z - centx_const[k][j][i].z;
852 } else if (j == my - 2) { // Backward difference
853 dxde = centx_const[k][j][i].x - centx_const[k][j-1][i].x;
854 dyde = centx_const[k][j][i].y - centx_const[k][j-1][i].y;
855 dzde = centx_const[k][j][i].z - centx_const[k][j-1][i].z;
856 } else { // Central difference
857 dxde = 0.5 * (centx_const[k][j+1][i].x - centx_const[k][j-1][i].x);
858 dyde = 0.5 * (centx_const[k][j+1][i].y - centx_const[k][j-1][i].y);
859 dzde = 0.5 * (centx_const[k][j+1][i].z - centx_const[k][j-1][i].z);
860 }
861
862 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
863 if (k == 1) { // Forward difference
864 dxdz = centx_const[k+1][j][i].x - centx_const[k][j][i].x;
865 dydz = centx_const[k+1][j][i].y - centx_const[k][j][i].y;
866 dzdz = centx_const[k+1][j][i].z - centx_const[k][j][i].z;
867 } else if (k == mz - 2) { // Backward difference
868 dxdz = centx_const[k][j][i].x - centx_const[k-1][j][i].x;
869 dydz = centx_const[k][j][i].y - centx_const[k-1][j][i].y;
870 dzdz = centx_const[k][j][i].z - centx_const[k-1][j][i].z;
871 } else { // Central difference
872 dxdz = 0.5 * (centx_const[k+1][j][i].x - centx_const[k-1][j][i].x);
873 dydz = 0.5 * (centx_const[k+1][j][i].y - centx_const[k-1][j][i].y);
874 dzdz = 0.5 * (centx_const[k+1][j][i].z - centx_const[k-1][j][i].z);
875 }
876
877 // --- Metric calculations (identical to legacy FormMetrics) ---
878 icsi[k][j][i].x = dyde * dzdz - dzde * dydz;
879 icsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
880 icsi[k][j][i].z = dxde * dydz - dyde * dxdz;
881
882 ieta[k][j][i].x = dydz * dzdc - dzdz * dydc;
883 ieta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
884 ieta[k][j][i].z = dxdz * dydc - dydz * dxdc;
885
886 izet[k][j][i].x = dydc * dzde - dzdc * dyde;
887 izet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
888 izet[k][j][i].z = dxdc * dyde - dydc * dxde;
889
890 iaj[k][j][i] = dxdc * icsi[k][j][i].x + dydc * icsi[k][j][i].y + dzdc * icsi[k][j][i].z;
891 if (PetscAbsScalar(iaj[k][j][i]) > 1e-12) {
892 iaj[k][j][i] = 1.0 / iaj[k][j][i];
893 }
894 }
895 }
896 }
897
898 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centx, &centx_const); CHKERRQ(ierr);
899 ierr = DMDAVecRestoreArray(user->fda, user->ICsi, &icsi); CHKERRQ(ierr);
900 ierr = DMDAVecRestoreArray(user->fda, user->IEta, &ieta); CHKERRQ(ierr);
901 ierr = DMDAVecRestoreArray(user->fda, user->IZet, &izet); CHKERRQ(ierr);
902 ierr = DMDAVecRestoreArray(user->da, user->IAj, &iaj); CHKERRQ(ierr);
903
904 // --- Part 3: Assemble global vectors and update local ghosts ---
905 ierr = VecAssemblyBegin(user->ICsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->ICsi); CHKERRQ(ierr);
906 ierr = VecAssemblyBegin(user->IEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IEta); CHKERRQ(ierr);
907 ierr = VecAssemblyBegin(user->IZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IZet); CHKERRQ(ierr);
908 ierr = VecAssemblyBegin(user->IAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->IAj); CHKERRQ(ierr);
909
910 ierr = UpdateLocalGhosts(user, "ICsi"); CHKERRQ(ierr);
911 ierr = UpdateLocalGhosts(user, "IEta"); CHKERRQ(ierr);
912 ierr = UpdateLocalGhosts(user, "IZet"); CHKERRQ(ierr);
913 ierr = UpdateLocalGhosts(user, "IAj"); CHKERRQ(ierr);
914
915 PetscFunctionReturn(0);
916}
917
918#undef __FUNCT__
919#define __FUNCT__ "ComputeJFaceMetrics"
920/**
921 * @brief Computes metrics centered on constant-j faces (j-faces).
922 *
923 * This function calculates the metric terms (`JCsi`, `JEta`, `JZet`) and the
924 * inverse Jacobian (`JAj`) located at the geometric center of each constant-j
925 * face. This is a critical step for staggered-grid finite difference schemes.
926 *
927 * The process is a direct and faithful refactoring of the corresponding logic
928 * from the legacy `FormMetrics` function:
929 * 1. It first calculates the physical (x,y,z) coordinates of the center of
930 * each i-face and stores them in the `user->Centy` vector.
931 * 2. It then uses a boundary-aware, second-order finite difference stencil on
932 * the `Centy` field to compute the derivatives (e.g., d(x)/d(csi)).
933 * - Central differences are used in the grid interior.
934 * - One-sided differences are used at the physical domain boundaries.
935 * 3. Finally, these derivatives are used to compute the final metric terms and
936 * the inverse Jacobian, which are stored in their respective `Vec` objects.
937 *
938 * @param user The UserCtx for a specific grid level. This function populates
939 * the `user->JCsi`, `user->JEta`, `user->JZet`, and `user->JAj` vectors.
940 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
941 */
942PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
943{
944 PetscErrorCode ierr;
945 DMDALocalInfo info;
946 Vec lCoords;
947 const Cmpnts ***coor;
948 Cmpnts ***centy; //***gs;
949 const Cmpnts ***centy_const;
950 Cmpnts ***jcsi, ***jeta, ***jzet;
951 PetscScalar ***jaj;
952 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
953
954 PetscFunctionBeginUser;
955 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing j-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
956
957 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
958 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
959 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
960 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
961 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
962 PetscInt gys = info.gys, gye = info.gys + info.gym;
963 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
964
965 PetscInt lxs = xs; PetscInt lxe = xe;
966 PetscInt lys = ys; PetscInt lye = ye;
967 PetscInt lzs = zs; PetscInt lze = ze;
968
969 if (xs==0) lxs = xs+1;
970 if (ys==0) lys = ys+1;
971 if (zs==0) lzs = zs+1;
972
973 if (xe==mx) lxe=xe-1;
974 if (ye==my) lye=ye-1;
975 if (ze==mz) lze=ze-1;
976
977 // --- Part 1: Calculate the location of i-face centers (Centx) ---
978 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
979 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
980 ierr = DMDAVecGetArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
981 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
982
983 // Loop over the ghosted region to calculate all local face centers
984 for (PetscInt k = gzs + 1; k < gze; k++) {
985 for (PetscInt j = gys; j < gye; j++) {
986 for (PetscInt i = gxs + 1; i < gxe; i++) {
987 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);
988 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);
989 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);
990 }
991 }
992 }
993
994 /*
995 if(ys==0){
996 for(PetscInt k=gzs+1;k < gze; k++){
997 for(PetscInt i=gxs+1;j < gxe; i++){
998 PetscInt j=0;
999 centy[k][j-1][i].x=centy[k][j][i].x;
1000 centy[k][j-1][i].y=centy[k][j][i].y-gs[k][j-2][i].y;
1001 centy[k][j-1][i].z=centy[k][j][i].z;
1002 }
1003 }
1004 }
1005 if (ye==my){
1006 for(PetscInt k=gzs+1; k<gze; k++) {
1007 for (PetscInt i=gxs+1; j<gxe;i++) {
1008 PetscInt j=my-1;
1009 centy[k][j][i].x=centy[k][j-1][i].x
1010 centy[k][j][i].y=centy[k][j-1][i].y+gs[k][j+2][i].y;
1011 centy[k][j][i].z=centy[k][j-1][i].z;
1012 }
1013 }
1014 }
1015 */
1016
1017 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1018 ierr = DMDAVecRestoreArray(user->fda, user->Centy, &centy); CHKERRQ(ierr);
1019 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1020
1021 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: j-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1022
1023 // --- Part 2: Calculate metrics using face-centered coordinates ---
1024 ierr = DMDAVecGetArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1025 ierr = DMDAVecGetArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1026 ierr = DMDAVecGetArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1027 ierr = DMDAVecGetArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1028 ierr = DMDAVecGetArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1029
1030 // Loop over the OWNED region where we will store the final metrics
1031 for (PetscInt k=lzs; k<lze; k++) {
1032 for (PetscInt j=ys; j<lye; j++) {
1033 for (PetscInt i=lxs; i<lxe; i++) {
1034
1035 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1036 if (i == 1) { // Forward difference at the domain's min-i boundary
1037 dxdc = centy_const[k][j][i+1].x - centy_const[k][j][i].x;
1038 dydc = centy_const[k][j][i+1].y - centy_const[k][j][i].y;
1039 dzdc = centy_const[k][j][i+1].z - centy_const[k][j][i].z;
1040 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1041 dxdc = centy_const[k][j][i].x - centy_const[k][j][i-1].x;
1042 dydc = centy_const[k][j][i].y - centy_const[k][j][i-1].y;
1043 dzdc = centy_const[k][j][i].z - centy_const[k][j][i-1].z;
1044 } else { // Central difference in the interior
1045 dxdc = 0.5 * (centy_const[k][j][i+1].x - centy_const[k][j][i-1].x);
1046 dydc = 0.5 * (centy_const[k][j][i+1].y - centy_const[k][j][i-1].y);
1047 dzdc = 0.5 * (centy_const[k][j][i+1].z - centy_const[k][j][i-1].z);
1048 }
1049
1050 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1051 if (j == 0) { // Forward difference
1052 dxde = centy_const[k][j+1][i].x - centy_const[k][j][i].x;
1053 dyde = centy_const[k][j+1][i].y - centy_const[k][j][i].y;
1054 dzde = centy_const[k][j+1][i].z - centy_const[k][j][i].z;
1055 } else if (j == my - 2) { // Backward difference
1056 dxde = centy_const[k][j][i].x - centy_const[k][j-1][i].x;
1057 dyde = centy_const[k][j][i].y - centy_const[k][j-1][i].y;
1058 dzde = centy_const[k][j][i].z - centy_const[k][j-1][i].z;
1059 } else { // Central difference
1060 dxde = 0.5 * (centy_const[k][j+1][i].x - centy_const[k][j-1][i].x);
1061 dyde = 0.5 * (centy_const[k][j+1][i].y - centy_const[k][j-1][i].y);
1062 dzde = 0.5 * (centy_const[k][j+1][i].z - centy_const[k][j-1][i].z);
1063 }
1064
1065 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1066 if (k == 1) { // Forward difference
1067 dxdz = centy_const[k+1][j][i].x - centy_const[k][j][i].x;
1068 dydz = centy_const[k+1][j][i].y - centy_const[k][j][i].y;
1069 dzdz = centy_const[k+1][j][i].z - centy_const[k][j][i].z;
1070 } else if (k == mz - 2) { // Backward difference
1071 dxdz = centy_const[k][j][i].x - centy_const[k-1][j][i].x;
1072 dydz = centy_const[k][j][i].y - centy_const[k-1][j][i].y;
1073 dzdz = centy_const[k][j][i].z - centy_const[k-1][j][i].z;
1074 } else { // Central difference
1075 dxdz = 0.5 * (centy_const[k+1][j][i].x - centy_const[k-1][j][i].x);
1076 dydz = 0.5 * (centy_const[k+1][j][i].y - centy_const[k-1][j][i].y);
1077 dzdz = 0.5 * (centy_const[k+1][j][i].z - centy_const[k-1][j][i].z);
1078 }
1079
1080 // --- Metric calculations (identical to legacy FormMetrics) ---
1081 jcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1082 jcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1083 jcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1084
1085 jeta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1086 jeta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1087 jeta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1088
1089 jzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1090 jzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1091 jzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1092
1093 jaj[k][j][i] = dxdc * jcsi[k][j][i].x + dydc * jcsi[k][j][i].y + dzdc * jcsi[k][j][i].z;
1094 if (PetscAbsScalar(jaj[k][j][i]) > 1e-12) {
1095 jaj[k][j][i] = 1.0 / jaj[k][j][i];
1096 }
1097 }
1098 }
1099 }
1100
1101 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centy, &centy_const); CHKERRQ(ierr);
1102 ierr = DMDAVecRestoreArray(user->fda, user->JCsi, &jcsi); CHKERRQ(ierr);
1103 ierr = DMDAVecRestoreArray(user->fda, user->JEta, &jeta); CHKERRQ(ierr);
1104 ierr = DMDAVecRestoreArray(user->fda, user->JZet, &jzet); CHKERRQ(ierr);
1105 ierr = DMDAVecRestoreArray(user->da, user->JAj, &jaj); CHKERRQ(ierr);
1106
1107 // --- Part 3: Assemble global vectors and update local ghosts ---
1108 ierr = VecAssemblyBegin(user->JCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JCsi); CHKERRQ(ierr);
1109 ierr = VecAssemblyBegin(user->JEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JEta); CHKERRQ(ierr);
1110 ierr = VecAssemblyBegin(user->JZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JZet); CHKERRQ(ierr);
1111 ierr = VecAssemblyBegin(user->JAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->JAj); CHKERRQ(ierr);
1112
1113 ierr = UpdateLocalGhosts(user, "JCsi"); CHKERRQ(ierr);
1114 ierr = UpdateLocalGhosts(user, "JEta"); CHKERRQ(ierr);
1115 ierr = UpdateLocalGhosts(user, "JZet"); CHKERRQ(ierr);
1116 ierr = UpdateLocalGhosts(user, "JAj"); CHKERRQ(ierr);
1117
1118 PetscFunctionReturn(0);
1119}
1120
1121#undef __FUNCT__
1122#define __FUNCT__ "ComputeJFaceMetrics"
1123/**
1124 * @brief Computes metrics centered on constant-k faces (k-faces).
1125 *
1126 * This function calculates the metric terms (`KCsi`, `KEta`, `KZet`) and the
1127 * inverse Jacobian (`KAj`) located at the geometric center of each constant-k
1128 * face. This is a critical step for staggered-grid finite difference schemes.
1129 *
1130 * The process is a direct and faithful refactoring of the corresponding logic
1131 * from the legacy `FormMetrics` function:
1132 * 1. It first calculates the physical (x,y,z) coordinates of the center of
1133 * each i-face and stores them in the `user->Centz` vector.
1134 * 2. It then uses a boundary-aware, second-order finite difference stencil on
1135 * the `Centz` field to compute the derivatives (e.g., d(x)/d(csi)).
1136 * - Central differences are used in the grid interior.
1137 * - One-sided differences are used at the physical domain boundaries.
1138 * 3. Finally, these derivatives are used to compute the final metric terms and
1139 * the inverse Jacobian, which are stored in their respective `Vec` objects.
1140 *
1141 * @param user The UserCtx for a specific grid level. This function populates
1142 * the `user->KCsi`, `user->KEta`, `user->KZet`, and `user->KAj` vectors.
1143 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1144 */
1145PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
1146{
1147 PetscErrorCode ierr;
1148 DMDALocalInfo info;
1149 Vec lCoords;
1150 const Cmpnts ***coor;
1151 Cmpnts ***centz; //***gs;
1152 const Cmpnts ***centz_const;
1153 Cmpnts ***kcsi, ***keta, ***kzet;
1154 PetscScalar ***kaj;
1155 PetscReal dxdc, dydc, dzdc, dxde, dyde, dzde, dxdz, dydz, dzdz;
1156
1157 PetscFunctionBeginUser;
1158 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Computing k-face metrics for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1159
1160 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1161 PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
1162 PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
1163 PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
1164 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1165 PetscInt gys = info.gys, gye = info.gys + info.gym;
1166 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1167
1168 PetscInt lxs = xs; PetscInt lxe = xe;
1169 PetscInt lys = ys; PetscInt lye = ye;
1170 PetscInt lzs = zs; PetscInt lze = ze;
1171
1172 if (xs==0) lxs = xs+1;
1173 if (ys==0) lys = ys+1;
1174 if (zs==0) lzs = zs+1;
1175
1176 if (xe==mx) lxe=xe-1;
1177 if (ye==my) lye=ye-1;
1178 if (ze==mz) lze=ze-1;
1179
1180 // --- Part 1: Calculate the location of i-face centers (Centx) ---
1181 ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
1182 ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1183 ierr = DMDAVecGetArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1184 // ierr = DMDAVecGetArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1185
1186 // Loop over the ghosted region to calculate all local face centers
1187 for (PetscInt k = gzs; k < gze; k++) {
1188 for (PetscInt j = gys; j < gye; j++) {
1189 for (PetscInt i = gxs + 1; i < gxe; i++) {
1190 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);
1191 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);
1192 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);
1193 }
1194 }
1195 }
1196
1197 /*
1198 if(zs==0){
1199 for(PetscInt j=gys+1;j < gye; j++){
1200 for(PetscInt i=gxs+1;j < gxe; i++){
1201 PetscInt k=0;
1202 centz[k-1][j][i].x=centz[k][j][i].x;
1203 centz[k-1][j][i].y=centz[k][j][i].y;
1204 centz[k-1][j][i].z=centz[k][j][i].z-gs[k-2][j][i].z;
1205 }
1206 }
1207 }
1208 if (ze==mz){
1209 for(PetscInt j=gys+1; j<gye; j++) {
1210 for (PetscInt i=gxs+1; j<gxe;i++) {
1211 PetscInt k=mz-1;
1212 centy[k][j][i].x=centy[k-1][j][i].x
1213 centy[k][j][i].y=centy[k-1][j][i].y;
1214 centz[k][j][i].z=centz[k-1][j][i].z+gs[k+2][j][1].z;
1215 }
1216 }
1217 }
1218 */
1219
1220 ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor); CHKERRQ(ierr);
1221 ierr = DMDAVecRestoreArray(user->fda, user->Centz, &centz); CHKERRQ(ierr);
1222 // ierr = DMDAVecRestoreArray(user->fda, user->lGridSpace,&gs); CHKERRQ(ierr);
1223
1224 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: k-face centers (Centx) calculated and ghosts updated.\n", user->simCtx->rank);
1225
1226 // --- Part 2: Calculate metrics using face-centered coordinates ---
1227 ierr = DMDAVecGetArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1228 ierr = DMDAVecGetArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1229 ierr = DMDAVecGetArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1230 ierr = DMDAVecGetArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1231 ierr = DMDAVecGetArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1232
1233 // Loop over the OWNED region where we will store the final metrics
1234 for (PetscInt k=zs; k<lze; k++) {
1235 for (PetscInt j=lys; j<lye; j++) {
1236 for (PetscInt i=lxs; i<lxe; i++) {
1237
1238 // --- Stencil Logic for d/dcsi (derivative in i-direction) ---
1239 if (i == 1) { // Forward difference at the domain's min-i boundary
1240 dxdc = centz_const[k][j][i+1].x - centz_const[k][j][i].x;
1241 dydc = centz_const[k][j][i+1].y - centz_const[k][j][i].y;
1242 dzdc = centz_const[k][j][i+1].z - centz_const[k][j][i].z;
1243 } else if (i == mx - 2) { // Backward difference at the domain's max-i boundary
1244 dxdc = centz_const[k][j][i].x - centz_const[k][j][i-1].x;
1245 dydc = centz_const[k][j][i].y - centz_const[k][j][i-1].y;
1246 dzdc = centz_const[k][j][i].z - centz_const[k][j][i-1].z;
1247 } else { // Central difference in the interior
1248 dxdc = 0.5 * (centz_const[k][j][i+1].x - centz_const[k][j][i-1].x);
1249 dydc = 0.5 * (centz_const[k][j][i+1].y - centz_const[k][j][i-1].y);
1250 dzdc = 0.5 * (centz_const[k][j][i+1].z - centz_const[k][j][i-1].z);
1251 }
1252
1253 // --- Stencil Logic for d/deta (derivative in j-direction) ---
1254 if (j == 1) { // Forward difference
1255 dxde = centz_const[k][j+1][i].x - centz_const[k][j][i].x;
1256 dyde = centz_const[k][j+1][i].y - centz_const[k][j][i].y;
1257 dzde = centz_const[k][j+1][i].z - centz_const[k][j][i].z;
1258 } else if (j == my - 2) { // Backward difference
1259 dxde = centz_const[k][j][i].x - centz_const[k][j-1][i].x;
1260 dyde = centz_const[k][j][i].y - centz_const[k][j-1][i].y;
1261 dzde = centz_const[k][j][i].z - centz_const[k][j-1][i].z;
1262 } else { // Central difference
1263 dxde = 0.5 * (centz_const[k][j+1][i].x - centz_const[k][j-1][i].x);
1264 dyde = 0.5 * (centz_const[k][j+1][i].y - centz_const[k][j-1][i].y);
1265 dzde = 0.5 * (centz_const[k][j+1][i].z - centz_const[k][j-1][i].z);
1266 }
1267
1268 // --- Stencil Logic for d/dzeta (derivative in k-direction) ---
1269 if (k == 0) { // Forward difference
1270 dxdz = centz_const[k+1][j][i].x - centz_const[k][j][i].x;
1271 dydz = centz_const[k+1][j][i].y - centz_const[k][j][i].y;
1272 dzdz = centz_const[k+1][j][i].z - centz_const[k][j][i].z;
1273 } else if (k == mz - 2) { // Backward difference
1274 dxdz = centz_const[k][j][i].x - centz_const[k-1][j][i].x;
1275 dydz = centz_const[k][j][i].y - centz_const[k-1][j][i].y;
1276 dzdz = centz_const[k][j][i].z - centz_const[k-1][j][i].z;
1277 } else { // Central difference
1278 dxdz = 0.5 * (centz_const[k+1][j][i].x - centz_const[k-1][j][i].x);
1279 dydz = 0.5 * (centz_const[k+1][j][i].y - centz_const[k-1][j][i].y);
1280 dzdz = 0.5 * (centz_const[k+1][j][i].z - centz_const[k-1][j][i].z);
1281 }
1282
1283 // --- Metric calculations (identical to legacy FormMetrics) ---
1284 kcsi[k][j][i].x = dyde * dzdz - dzde * dydz;
1285 kcsi[k][j][i].y = -dxde * dzdz + dzde * dxdz;
1286 kcsi[k][j][i].z = dxde * dydz - dyde * dxdz;
1287
1288 keta[k][j][i].x = dydz * dzdc - dzdz * dydc;
1289 keta[k][j][i].y = -dxdz * dzdc + dzdz * dxdc;
1290 keta[k][j][i].z = dxdz * dydc - dydz * dxdc;
1291
1292 kzet[k][j][i].x = dydc * dzde - dzdc * dyde;
1293 kzet[k][j][i].y = -dxdc * dzde + dzdc * dxde;
1294 kzet[k][j][i].z = dxdc * dyde - dydc * dxde;
1295
1296 kaj[k][j][i] = dxdc * kcsi[k][j][i].x + dydc * kcsi[k][j][i].y + dzdc * kcsi[k][j][i].z;
1297 if (PetscAbsScalar(kaj[k][j][i]) > 1e-12) {
1298 kaj[k][j][i] = 1.0 / kaj[k][j][i];
1299 }
1300 }
1301 }
1302 }
1303
1304 ierr = DMDAVecRestoreArrayRead(user->fda, user->Centz, &centz_const); CHKERRQ(ierr);
1305 ierr = DMDAVecRestoreArray(user->fda, user->KCsi, &kcsi); CHKERRQ(ierr);
1306 ierr = DMDAVecRestoreArray(user->fda, user->KEta, &keta); CHKERRQ(ierr);
1307 ierr = DMDAVecRestoreArray(user->fda, user->KZet, &kzet); CHKERRQ(ierr);
1308 ierr = DMDAVecRestoreArray(user->da, user->KAj, &kaj); CHKERRQ(ierr);
1309
1310 // --- Part 3: Assemble global vectors and update local ghosts ---
1311 ierr = VecAssemblyBegin(user->KCsi); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KCsi); CHKERRQ(ierr);
1312 ierr = VecAssemblyBegin(user->KEta); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KEta); CHKERRQ(ierr);
1313 ierr = VecAssemblyBegin(user->KZet); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KZet); CHKERRQ(ierr);
1314 ierr = VecAssemblyBegin(user->KAj); CHKERRQ(ierr); ierr = VecAssemblyEnd(user->KAj); CHKERRQ(ierr);
1315
1316 ierr = UpdateLocalGhosts(user, "KCsi"); CHKERRQ(ierr);
1317 ierr = UpdateLocalGhosts(user, "KEta"); CHKERRQ(ierr);
1318 ierr = UpdateLocalGhosts(user, "KZet"); CHKERRQ(ierr);
1319 ierr = UpdateLocalGhosts(user, "KAj"); CHKERRQ(ierr);
1320
1321 PetscFunctionReturn(0);
1322}
1323
1324#undef __FUNCT__
1325#define __FUNCT__ "MetricsDivergence"
1326/**
1327 * @brief Performs a diagnostic check on the divergence of the face area metric vectors.
1328 *
1329 * For a closed cell, the sum of the face area vectors should be zero (Gauss's
1330 * divergence theorem). This function computes a measure of this divergence and
1331 * reports the maximum value over the domain. A small value indicates a
1332 * well-formed grid. This is a direct adaptation of the legacy function.
1333 *
1334 * @param user The UserCtx for a specific grid level (typically the finest).
1335 * @return PetscErrorCode 0 on success, or a PETSc error code on failure.
1336 */
1337PetscErrorCode MetricsDivergence(UserCtx *user)
1338{
1339 PetscErrorCode ierr;
1340 DMDALocalInfo info;
1341 Vec Div;
1342 PetscReal ***div_arr;
1343 const Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
1344 const PetscScalar ***aj_arr;
1345 PetscReal maxdiv;
1346 PetscInt max_idx;
1347
1348 PetscFunctionBeginUser;
1349 LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Performing metric divergence check for level %d block %d...\n", user->simCtx->rank, user->thislevel, user->_this);
1350
1351 ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
1352
1353 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
1354 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
1355 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
1356 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &aj_arr); CHKERRQ(ierr);
1357
1358 ierr = VecDuplicate(user->P, &Div); CHKERRQ(ierr);
1359 ierr = VecSet(Div, 0.0); CHKERRQ(ierr);
1360 ierr = DMDAVecGetArray(user->da, Div, &div_arr); CHKERRQ(ierr);
1361
1362 // Loop over the interior OWNED cells
1363 for (PetscInt k=info.xs+1; k<info.xs+info.xm; k++) {
1364 for (PetscInt j=info.ys+1; j<info.ys+info.ym; j++) {
1365 for (PetscInt i=info.xs+1; i<info.xs+info.xm; i++) {
1366 PetscReal div_x = (csi_arr[k][j][i].x - csi_arr[k][j][i-1].x) + (eta_arr[k][j][i].x - eta_arr[k][j-1][i].x) + (zet_arr[k][j][i].x - zet_arr[k-1][j][i].x);
1367 PetscReal div_y = (csi_arr[k][j][i].y - csi_arr[k][j][i-1].y) + (eta_arr[k][j][i].y - eta_arr[k][j-1][i].y) + (zet_arr[k][j][i].y - zet_arr[k-1][j][i].y);
1368 PetscReal div_z = (csi_arr[k][j][i].z - csi_arr[k][j][i-1].z) + (eta_arr[k][j][i].z - eta_arr[k][j-1][i].z) + (zet_arr[k][j][i].z - zet_arr[k-1][j][i].z);
1369 div_arr[k][j][i] = PetscAbsScalar((div_x + div_y + div_z) * aj_arr[k][j][i]);
1370 }
1371 }
1372 }
1373
1374 ierr = DMDAVecRestoreArray(user->da, Div, &div_arr); CHKERRQ(ierr);
1375 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
1376 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
1377 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
1378 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &aj_arr); CHKERRQ(ierr);
1379
1380 ierr = VecMax(Div, &max_idx, &maxdiv); CHKERRQ(ierr);
1381 LOG_ALLOW(GLOBAL, LOG_INFO, "Maximum metric divergence: %e\n", maxdiv);
1382
1383 ierr = VecDestroy(&Div); CHKERRQ(ierr);
1384 PetscFunctionReturn(0);
1385}
1386
1387static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
1388
1389{
1390 PetscInt nidx;
1391 DMDALocalInfo info = user->info;
1392
1393 PetscInt mx = info.mx, my = info.my;
1394
1395 AO ao;
1396 DMDAGetAO(user->da, &ao);
1397 nidx=i+j*mx+k*mx*my;
1398
1399 AOApplicationToPetsc(ao,1,&nidx);
1400
1401 return (nidx);
1402}
1403
1404/**
1405 * @brief Computes the divergence of the grid metrics and identifies the maximum value.
1406 *
1407 * This function serves as a diagnostic tool to assess the quality of the grid
1408 * metrics. It calculates the divergence of the face metrics (Csi, Eta, Zet)
1409 * and scales it by the inverse of the cell Jacobian. The maximum divergence
1410 * value is then located, and its grid coordinates are printed to the console,
1411 * helping to pinpoint areas of potential grid quality issues.
1412 *
1413 * @param user The UserCtx, containing all necessary grid data.
1414 * @return PetscErrorCode
1415 */
1416PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
1417{
1418 DM da = user->da, fda = user->fda;
1419 DMDALocalInfo info = user->info;
1420 PetscInt xs = info.xs, xe = info.xs + info.xm;
1421 PetscInt ys = info.ys, ye = info.ys + info.ym;
1422 PetscInt zs = info.zs, ze = info.zs + info.zm;
1423 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1424 PetscInt lxs, lys, lzs, lxe, lye, lze;
1425 PetscInt i, j, k;
1426 Vec Div;
1427 PetscReal ***div, ***aj;
1428 Cmpnts ***csi, ***eta, ***zet;
1429 PetscReal maxdiv;
1430
1431 PetscFunctionBeginUser;
1432
1433 lxs = xs; lxe = xe;
1434 lys = ys; lye = ye;
1435 lzs = zs; lze = ze;
1436
1437 if (xs == 0) lxs = xs + 1;
1438 if (ys == 0) lys = ys + 1;
1439 if (zs == 0) lzs = zs + 1;
1440
1441 if (xe == mx) lxe = xe - 1;
1442 if (ye == my) lye = ye - 1;
1443 if (ze == mz) lze = ze - 1;
1444
1445 DMDAVecGetArray(fda, user->lCsi, &csi);
1446 DMDAVecGetArray(fda, user->lEta, &eta);
1447 DMDAVecGetArray(fda, user->lZet, &zet);
1448 DMDAVecGetArray(da, user->lAj, &aj);
1449
1450 VecDuplicate(user->P, &Div);
1451 VecSet(Div, 0.);
1452 DMDAVecGetArray(da, Div, &div);
1453
1454 for (k = lzs; k < lze; k++) {
1455 for (j = lys; j < lye; j++) {
1456 for (i = lxs; i < lxe; i++) {
1457 PetscReal divergence = (csi[k][j][i].x - csi[k][j][i-1].x +
1458 eta[k][j][i].x - eta[k][j-1][i].x +
1459 zet[k][j][i].x - zet[k-1][j][i].x +
1460 csi[k][j][i].y - csi[k][j][i-1].y +
1461 eta[k][j][i].y - eta[k][j-1][i].y +
1462 zet[k][j][i].y - zet[k-1][j][i].y +
1463 csi[k][j][i].z - csi[k][j][i-1].z +
1464 eta[k][j][i].z - eta[k][j-1][i].z +
1465 zet[k][j][i].z - zet[k-1][j][i].z) * aj[k][j][i];
1466 div[k][j][i] = fabs(divergence);
1467 }
1468 }
1469 }
1470
1471 DMDAVecRestoreArray(da, Div, &div);
1472
1473 PetscReal MaxFlatIndex;
1474 VecMax(Div, &MaxFlatIndex, &maxdiv);
1475 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence is %e at flat index %d.\n",maxdiv,MaxFlatIndex);
1476
1477 for (k=zs; k<ze; k++) {
1478 for (j=ys; j<ye; j++) {
1479 for (i=xs; i<xe; i++) {
1480 if (Gidx(i,j,k,user) == MaxFlatIndex) {
1481 LOG_ALLOW(GLOBAL,LOG_INFO,"The Maximum Metric Divergence(%e) is at location [%d][%d][%d]. \n", maxdiv,k,j,i);
1482 }
1483 }
1484 }
1485 }
1486
1487
1488 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1489 DMDAVecRestoreArray(fda, user->lEta, &eta);
1490 DMDAVecRestoreArray(fda, user->lZet, &zet);
1491 DMDAVecRestoreArray(da, user->lAj, &aj);
1492 VecDestroy(&Div);
1493
1494 PetscFunctionReturn(0);
1495}
1496
1497/**
1498 * @brief Computes the max-min values of the grid metrics.
1499 *
1500 * This function serves as a diagnostic tool to assess the quality of the grid
1501 * metrics. It calculates the bounds of the face metrics (Csi, Eta, Zet).
1502 *
1503 * @param user The UserCtx, containing all necessary grid data.
1504 * @return PetscErrorCode
1505 */
1506PetscErrorCode ComputeMetricNorms(UserCtx *user)
1507{
1508
1509 DMDALocalInfo info = user->info;
1510 PetscInt xs = info.xs, xe = info.xs + info.xm;
1511 PetscInt ys = info.ys, ye = info.ys + info.ym;
1512 PetscInt zs = info.zs, ze = info.zs + info.zm;
1513 PetscInt i, j, k;
1514
1515 PetscFunctionBeginUser;
1516
1517 PetscReal CsiMax, EtaMax, ZetMax;
1518 PetscReal ICsiMax, IEtaMax, IZetMax;
1519 PetscReal JCsiMax, JEtaMax, JZetMax;
1520 PetscReal KCsiMax, KEtaMax, KZetMax;
1521 PetscReal AjMax, IAjMax, JAjMax, KAjMax;
1522
1523 PetscInt CsiMaxArg, EtaMaxArg, ZetMaxArg;
1524 PetscInt ICsiMaxArg, IEtaMaxArg, IZetMaxArg;
1525 PetscInt JCsiMaxArg, JEtaMaxArg, JZetMaxArg;
1526 PetscInt KCsiMaxArg, KEtaMaxArg, KZetMaxArg;
1527 PetscInt AjMaxArg, IAjMaxArg, JAjMaxArg, KAjMaxArg;
1528
1529 // Max Values
1530 VecMax(user->lCsi,&CsiMaxArg,&CsiMax);
1531 VecMax(user->lEta,&EtaMaxArg,&EtaMax);
1532 VecMax(user->lZet,&ZetMaxArg,&ZetMax);
1533
1534 VecMax(user->lICsi,&ICsiMaxArg,&ICsiMax);
1535 VecMax(user->lIEta,&IEtaMaxArg,&IEtaMax);
1536 VecMax(user->lIZet,&IZetMaxArg,&IZetMax);
1537
1538 VecMax(user->lJCsi,&JCsiMaxArg,&JCsiMax);
1539 VecMax(user->lJEta,&JEtaMaxArg,&JEtaMax);
1540 VecMax(user->lJZet,&JZetMaxArg,&JZetMax);
1541
1542 VecMax(user->lKCsi,&KCsiMaxArg,&KCsiMax);
1543 VecMax(user->lKEta,&KEtaMaxArg,&KEtaMax);
1544 VecMax(user->lKZet,&KZetMaxArg,&KZetMax);
1545
1546 VecMax(user->lAj,&AjMaxArg,&AjMax);
1547 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1548 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1549 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1550
1551 VecMax(user->lAj,&AjMaxArg,&AjMax);
1552 VecMax(user->lIAj,&IAjMaxArg,&IAjMax);
1553 VecMax(user->lJAj,&JAjMaxArg,&JAjMax);
1554 VecMax(user->lKAj,&KAjMaxArg,&KAjMax);
1555
1556 LOG_ALLOW(GLOBAL,LOG_INFO," Metric Norms for MG level %d .\n",user->thislevel);
1557
1558 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: CsiMax = %le, EtaMax = %le, ZetMax = %le.\n",CsiMax,EtaMax,ZetMax);
1559 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: ICsiMax = %le, IEtaMax = %le, IZetMax = %le.\n",ICsiMax,IEtaMax,IZetMax);
1560 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: JCsiMax = %le, JEtaMax = %le, JZetMax = %le.\n",JCsiMax,JEtaMax,JZetMax);
1561 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Metric Values are: KCsiMax = %le, KEtaMax = %le, KZetMax = %le.\n",KCsiMax,KEtaMax,KZetMax);
1562 LOG_ALLOW(GLOBAL,LOG_INFO,"The Max Volumes(Inverse) are: Aj = %le, IAj = %le, JAj = %le, KAj = %le.\n",AjMax,IAjMax,JAjMax,KAjMax);
1563
1564 for (k=zs; k<ze; k++) {
1565 for (j=ys; j<ye; j++) {
1566 for (i=xs; i<xe; i++) {
1567 if (Gidx(i,j,k,user) == CsiMaxArg) {
1568 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Csi = %le is at [%d][%d][%d] \n", CsiMax,k,j,i);
1569 }
1570 if (Gidx(i,j,k,user) == EtaMaxArg) {
1571 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Eta = %le is at [%d][%d][%d] \n", EtaMax,k,j,i);
1572 }
1573 if (Gidx(i,j,k,user) == ZetMaxArg) {
1574 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Zet = %le is at [%d][%d][%d] \n", ZetMax,k,j,i);
1575 }
1576 if (Gidx(i,j,k,user) == ICsiMaxArg) {
1577 LOG_ALLOW(GLOBAL,LOG_INFO,"Max ICsi = %le is at [%d][%d][%d] \n", ICsiMax,k,j,i);
1578 }
1579 if (Gidx(i,j,k,user) == IEtaMaxArg) {
1580 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IEta = %le is at [%d][%d][%d] \n", IEtaMax,k,j,i);
1581 }
1582 if (Gidx(i,j,k,user) == IZetMaxArg) {
1583 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IZet = %le is at [%d][%d][%d] \n", IZetMax,k,j,i);
1584 }
1585 if (Gidx(i,j,k,user) == JCsiMaxArg) {
1586 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JCsi = %le is at [%d][%d][%d] \n", JCsiMax,k,j,i);
1587 }
1588 if (Gidx(i,j,k,user) == JEtaMaxArg) {
1589 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JEta = %le is at [%d][%d][%d] \n", JEtaMax,k,j,i);
1590 }
1591 if (Gidx(i,j,k,user) == JZetMaxArg) {
1592 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JZet = %le is at [%d][%d][%d] \n", JZetMax,k,j,i);
1593 }
1594 if (Gidx(i,j,k,user) == KCsiMaxArg) {
1595 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KCsi = %le is at [%d][%d][%d] \n", KCsiMax,k,j,i);
1596 }
1597 if (Gidx(i,j,k,user) == KEtaMaxArg) {
1598 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KEta = %le is at [%d][%d][%d] \n", KEtaMax,k,j,i);
1599 }
1600 if (Gidx(i,j,k,user) == KZetMaxArg) {
1601 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KZet = %le is at [%d][%d][%d] \n", KZetMax,k,j,i);
1602 }
1603 if (Gidx(i,j,k,user) == AjMaxArg) {
1604 LOG_ALLOW(GLOBAL,LOG_INFO,"Max Aj = %le is at [%d][%d][%d] \n", AjMax,k,j,i);
1605 }
1606 if (Gidx(i,j,k,user) == IAjMaxArg) {
1607 LOG_ALLOW(GLOBAL,LOG_INFO,"Max IAj = %le is at [%d][%d][%d] \n", IAjMax,k,j,i);
1608 }
1609 if (Gidx(i,j,k,user) == JAjMaxArg) {
1610 LOG_ALLOW(GLOBAL,LOG_INFO,"Max JAj = %le is at [%d][%d][%d] \n", JAjMax,k,j,i);
1611 }
1612 if (Gidx(i,j,k,user) == KAjMaxArg) {
1613 LOG_ALLOW(GLOBAL,LOG_INFO,"Max KAj = %le is at [%d][%d][%d] \n", KAjMax,k,j,i);
1614 }
1615 }
1616 }
1617 }
1618
1619 /*
1620 VecView(user->lCsi,PETSC_VIEWER_STDOUT_WORLD);
1621 VecView(user->lEta,PETSC_VIEWER_STDOUT_WORLD);
1622 VecView(user->lZet,PETSC_VIEWER_STDOUT_WORLD);
1623 */
1624
1625 PetscFunctionReturn(0);
1626}
1627
1628/**
1629 * @brief Orchestrates the calculation of all grid metrics.
1630 *
1631 * This function iterates through every UserCtx in the multigrid and multi-block
1632 * hierarchy. For each context, it calls a series of modern, modular helper
1633 * functions to compute the face metrics (Csi, Eta, Zet), the cell-centered
1634 * inverse Jacobian (Aj), and to validate the grid's orientation.
1635 *
1636 * @param simCtx The master SimCtx, containing the configured UserCtx hierarchy.
1637 * @return PetscErrorCode
1638 */
1639PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
1640{
1641 PetscErrorCode ierr;
1642 UserMG *usermg = &simCtx->usermg;
1643 MGCtx *mgctx = usermg->mgctx;
1644 PetscInt nblk = simCtx->block_number;
1645
1646 PetscFunctionBeginUser;
1647 LOG_ALLOW(GLOBAL, LOG_INFO, "Calculating grid metrics for all levels and blocks...\n");
1648
1649 // Loop through all levels and all blocks
1650 for (PetscInt level = usermg->mglevels -1 ; level >=0; level--) {
1651 for (PetscInt bi = 0; bi < nblk; bi++) {
1652 UserCtx *user = &mgctx[level].user[bi];
1653 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Calculating metrics for level %d, block %d\n", simCtx->rank, level, bi);
1654
1655 // Call the modern, modular helper functions for each UserCtx.
1656 // These functions are self-contained and operate on the data within the provided context.
1657 ierr = ComputeFaceMetrics(user); CHKERRQ(ierr);
1658 ierr = ComputeCellCenteredJacobianInverse(user); CHKERRQ(ierr);
1659 ierr = CheckAndFixGridOrientation(user); CHKERRQ(ierr);
1660 ierr = ComputeCellCentersAndSpacing(user); CHKERRQ(ierr);
1661 ierr = ComputeIFaceMetrics(user); CHKERRQ(ierr);
1662 ierr = ComputeJFaceMetrics(user); CHKERRQ(ierr);
1663 ierr = ComputeKFaceMetrics(user); CHKERRQ(ierr);
1664
1665 // Diagnostics
1666 ierr = ComputeMetricNorms(user);
1667 if (level == usermg->mglevels - 1) {
1668 ierr = ComputeMetricsDivergence(user); CHKERRQ(ierr);
1669 }
1670 }
1671 }
1672
1673 LOG_ALLOW(GLOBAL, LOG_INFO, "Grid metrics calculation complete.\n");
1674 PetscFunctionReturn(0);
1675}
1676/* ------------------------------------------------------------------------- */
1677/* End of Metric.c */
1678
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Orchestrates the calculation of all grid metrics.
Definition Metric.c:1639
MetricJacobian
PetscErrorCode MetricJacobian(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, PetscReal J[3][3], PetscReal *detJ)
Compute Jacobian matrix and its determinant at (xi,eta,zta).
Definition Metric.c:93
ComputeMetricNorms
PetscErrorCode ComputeMetricNorms(UserCtx *user)
Computes the max-min values of the grid metrics.
Definition Metric.c:1506
TrilinearBlend
static void TrilinearBlend(const Cmpnts V[8], PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Map logical coordinates to physical space using trilinear basis.
Definition Metric.c:48
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition Metric.c:1387
CheckAndFixGridOrientation
PetscErrorCode CheckAndFixGridOrientation(UserCtx *user)
Ensure a right-handed metric basis (Csi, Eta, Zet) and a positive Jacobian (Aj) over the whole domain...
Definition Metric.c:195
ComputeCellCentersAndSpacing
PetscErrorCode ComputeCellCentersAndSpacing(UserCtx *user)
Computes the physical location of cell centers and the spacing between them.
Definition Metric.c:644
MetricGetCellVertices
PetscErrorCode MetricGetCellVertices(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, Cmpnts V[8])
Extract the eight vertex coordinates of the hexahedral cell (i,j,k).
Definition Metric.c:24
ComputeJFaceMetrics
PetscErrorCode ComputeJFaceMetrics(UserCtx *user)
Computes metrics centered on constant-j faces (j-faces).
Definition Metric.c:942
ComputeMetricsDivergence
PetscErrorCode ComputeMetricsDivergence(UserCtx *user)
Computes the divergence of the grid metrics and identifies the maximum value.
Definition Metric.c:1416
ComputeFaceMetrics
PetscErrorCode ComputeFaceMetrics(UserCtx *user)
Computes the primary face metric components (Csi, Eta, Zet), including boundary extrapolation,...
Definition Metric.c:286
ComputeCellCenteredJacobianInverse
PetscErrorCode ComputeCellCenteredJacobianInverse(UserCtx *user)
Calculates the cell-centered inverse Jacobian determinant (1/J), including boundary extrapolation,...
Definition Metric.c:491
ComputeKFaceMetrics
PetscErrorCode ComputeKFaceMetrics(UserCtx *user)
Computes metrics centered on constant-k faces (k-faces).
Definition Metric.c:1145
MetricsDivergence
PetscErrorCode MetricsDivergence(UserCtx *user)
Performs a diagnostic check on the divergence of the face area metric vectors.
Definition Metric.c:1337
MetricLogicalToPhysical
PetscErrorCode MetricLogicalToPhysical(UserCtx *user, const Cmpnts ***X, PetscInt i, PetscInt j, PetscInt k, PetscReal xi, PetscReal eta, PetscReal zta, Cmpnts *Xp)
Public wrapper: map (cell index, ξ,η,ζ) to (x,y,z).
Definition Metric.c:68
ComputeIFaceMetrics
PetscErrorCode ComputeIFaceMetrics(UserCtx *user)
Computes metrics centered on constant-i faces (i-faces).
Definition Metric.c:739
MetricVelocityContravariant
PetscErrorCode MetricVelocityContravariant(const PetscReal J[3][3], PetscReal detJ, const PetscReal u[3], PetscReal uc[3])
Convert physical velocity (u,v,w) to contravariant components (u^xi, u^eta, u^zta).
Definition Metric.c:141
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:319
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:274
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:207
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:779
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:673
UserCtx::JCsi
Vec JCsi
Definition variables.h:676
UserCtx::KAj
Vec KAj
Definition variables.h:677
MGCtx::user
UserCtx * user
Definition variables.h:418
UserCtx::JEta
Vec JEta
Definition variables.h:676
UserCtx::Zet
Vec Zet
Definition variables.h:673
SimCtx::rank
PetscMPIInt rank
Definition variables.h:516
SimCtx::block_number
PetscInt block_number
Definition variables.h:562
UserCtx::lIEta
Vec lIEta
Definition variables.h:675
UserCtx::lIZet
Vec lIZet
Definition variables.h:675
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:633
UserCtx::IZet
Vec IZet
Definition variables.h:675
UserCtx::Centz
Vec Centz
Definition variables.h:674
UserCtx::IEta
Vec IEta
Definition variables.h:675
UserCtx::lZet
Vec lZet
Definition variables.h:673
SimCtx::usermg
UserMG usermg
Definition variables.h:599
UserCtx::Csi
Vec Csi
Definition variables.h:673
UserCtx::lIAj
Vec lIAj
Definition variables.h:675
UserCtx::_this
PetscInt _this
Definition variables.h:643
UserCtx::lKEta
Vec lKEta
Definition variables.h:677
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:676
Cmpnts::x
PetscScalar x
Definition variables.h:100
UserCtx::JZet
Vec JZet
Definition variables.h:676
UserCtx::Centx
Vec Centx
Definition variables.h:674
UserCtx::lKZet
Vec lKZet
Definition variables.h:677
UserCtx::Eta
Vec Eta
Definition variables.h:673
UserCtx::lJEta
Vec lJEta
Definition variables.h:676
UserCtx::lCsi
Vec lCsi
Definition variables.h:673
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:689
UserCtx::ICsi
Vec ICsi
Definition variables.h:675
Cmpnts::z
PetscScalar z
Definition variables.h:100
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:677
UserMG::mglevels
PetscInt mglevels
Definition variables.h:425
UserCtx::lJZet
Vec lJZet
Definition variables.h:676
UserCtx::IAj
Vec IAj
Definition variables.h:675
UserCtx::JAj
Vec JAj
Definition variables.h:676
UserCtx::KEta
Vec KEta
Definition variables.h:677
UserCtx::lAj
Vec lAj
Definition variables.h:673
UserCtx::lICsi
Vec lICsi
Definition variables.h:675
UserCtx::GridOrientation
PetscInt GridOrientation
Definition variables.h:643
UserCtx::info
DMDALocalInfo info
Definition variables.h:637
Cmpnts::y
PetscScalar y
Definition variables.h:100
UserCtx::lEta
Vec lEta
Definition variables.h:673
UserCtx::KZet
Vec KZet
Definition variables.h:677
UserCtx::Cent
Vec Cent
Definition variables.h:673
UserCtx::KCsi
Vec KCsi
Definition variables.h:677
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:428
UserCtx::Centy
Vec Centy
Definition variables.h:674
UserCtx::lJAj
Vec lJAj
Definition variables.h:676
UserCtx::lKAj
Vec lKAj
Definition variables.h:677
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:99
MGCtx
Context for Multigrid operations.
Definition variables.h:417
SimCtx
The master context for the entire simulation.
Definition variables.h:513
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:630
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:424