PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
poisson.c
Go to the documentation of this file.
1// In src/poisson.c (or wherever PoissonSolver_MG is)
2#include "poisson.h" // The new header for this file
3#include "logging.h"
4
5#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user) \
6 if ((user->isc)) { \
7 ic = i; \
8 ia = 0; \
9 } \
10 else { \
11 ic = (i+1) / 2; \
12 ia = (i - 2 * (ic)) == 0 ? 1 : -1; \
13 if (i==1 || i==mx-2) ia = 0; \
14 }\
15 if ((user->jsc)) { \
16 jc = j; \
17 ja = 0; \
18 } \
19 else { \
20 jc = (j+1) / 2; \
21 ja = (j - 2 * (jc)) == 0 ? 1 : -1; \
22 if (j==1 || j==my-2) ja = 0; \
23 } \
24 if ((user->ksc)) { \
25 kc = k; \
26 ka = 0; \
27 } \
28 else { \
29 kc = (k+1) / 2; \
30 ka = (k - 2 * (kc)) == 0 ? 1 : -1; \
31 if (k==1 || k==mz-2) ka = 0; \
32 } \
33 if (ka==-1 && nvert_c[kc-1][jc][ic] > 0.1) ka = 0; \
34 else if (ka==1 && nvert_c[kc+1][jc][ic] > 0.1) ka = 0; \
35 if (ja==-1 && nvert_c[kc][jc-1][ic] > 0.1) ja = 0; \
36 else if (ja==1 && nvert_c[kc][jc+1][ic] > 0.1) ja = 0; \
37 if (ia==-1 && nvert_c[kc][jc][ic-1] > 0.1) ia = 0; \
38 else if (ia==1 && nvert_c[kc][jc][ic+1] > 0.1) ia = 0;
39
40static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
41
42{
43 PetscInt nidx;
44 DMDALocalInfo info = user->info;
45
46 PetscInt mx = info.mx, my = info.my;
47
48 AO ao;
49 DMDAGetAO(user->da, &ao);
50 nidx=i+j*mx+k*mx*my;
51
52 AOApplicationToPetsc(ao,1,&nidx);
53
54 return (nidx);
55}
56
57#undef __FUNCT__
58#define __FUNCT__ "GridRestriction"
59
60static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k,
61 PetscInt *ih, PetscInt *jh, PetscInt *kh,
62 UserCtx *user)
63{
64 PetscFunctionBeginUser;
65 PROFILE_FUNCTION_BEGIN;
66 if ((user->isc)) {
67 *ih = i;
68 }
69 else {
70 *ih = 2 * i;
71 }
72
73 if ((user->jsc)) {
74 *jh = j;
75 }
76 else {
77 *jh = 2 * j;
78 }
79
80 if ((user->ksc)) {
81 *kh = k;
82 }
83 else {
84 *kh = 2 * k;
85 }
86
87 PROFILE_FUNCTION_END;
88 PetscFunctionReturn(0);
89}
90
91#include "solvers.h" // Or your main project header
92
93#undef __FUNCT__
94#define __FUNCT__ "CorrectChannelFluxProfile"
95/**
96 * @brief Enforces a constant volumetric flux profile along the entire length of a driven periodic channel.
97 *
98 * This function is a "hard" corrector, called at the end of the projection step.
99 * The projection ensures the velocity field is divergence-free (3D continuity), but this
100 * function enforces a stricter 1D continuity condition (`Flux(plane) = constant`)
101 * required for physically realistic, fully-developed periodic channel/pipe flow.
102 *
103 * The process is as follows:
104 * 1. Introspects the boundary condition handlers to detect if a `DRIVEN_` flow is active
105 * and in which direction ('X', 'Y', or 'Z'). If none is found, it exits.
106 * 2. Measures the current volumetric flux through *every single cross-sectional plane*
107 * in the driven direction.
108 * 3. For each plane, it calculates the velocity correction required to make its flux
109 * match the global `targetVolumetricFlux` (which was set by the controller).
110 * 4. It applies this spatially-uniform (but plane-dependent) velocity correction directly
111 * to the `ucont` field, ensuring `Flux(plane) = TargetFlux` for all planes.
112 *
113 * @param user The UserCtx containing the simulation state for a single block.
114 * @return PetscErrorCode 0 on success.
115 */
116PetscErrorCode CorrectChannelFluxProfile(UserCtx *user)
117{
118 PetscErrorCode ierr;
119 SimCtx *simCtx = user->simCtx;
120
121 PetscFunctionBeginUser;
122
123 // --- Step 1: Discover if and where a driven flow is active ---
124 char drivenDirection = ' ';
125 for (int i = 0; i < 6; i++) {
126 BCHandlerType handler_type = user->boundary_faces[i].handler_type;
127 if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
128 handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
129 {
130 switch (user->boundary_faces[i].face_id) {
131 case BC_FACE_NEG_X: case BC_FACE_POS_X: drivenDirection = 'X'; break;
132 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: drivenDirection = 'Y'; break;
133 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: drivenDirection = 'Z'; break;
134 }
135 break;
136 }
137 }
138
139 // --- Step 2: Early exit if no driven flow is configured ---
140 if (drivenDirection == ' ') {
141 PetscFunctionReturn(0);
142 }
143
144 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Starting channel flux profile correction in '%c' direction...\n",
145 simCtx->rank, user->_this, drivenDirection);
146
147 // --- Step 3: Setup and Get PETSc Array Pointers ---
148 DMDALocalInfo info = user->info;
149 PetscInt i, j, k;
150 PetscInt mx = info.mx, my = info.my, mz = info.mz;
151 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
152 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
153 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
154 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
155 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
156 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
157
158 Cmpnts ***ucont, ***csi, ***eta, ***zet;
159 PetscReal ***nvert;
160 ierr = DMDAVecGetArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
161 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
162 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
163 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
164 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
165
166 // --- Step 4: Allocate Memory for Profile Arrays based on direction ---
167 PetscInt n_planes = 0;
168 switch (drivenDirection) {
169 case 'X': n_planes = mx - 1; break;
170 case 'Y': n_planes = my - 1; break;
171 case 'Z': n_planes = mz - 1; break;
172 }
173
174 PetscReal *localFluxProfile, *globalFluxProfile, *correctionProfile;
175 ierr = PetscMalloc1(n_planes, &localFluxProfile); CHKERRQ(ierr);
176 ierr = PetscMalloc1(n_planes, &globalFluxProfile); CHKERRQ(ierr);
177 ierr = PetscMalloc1(n_planes, &correctionProfile); CHKERRQ(ierr);
178 ierr = PetscMemzero(localFluxProfile, n_planes * sizeof(PetscReal)); CHKERRQ(ierr);
179
180 // --- Step 5: Calculate Total Cross-Sectional Area and Measure Flux Profile ---
181 PetscReal localArea = 0.0, globalArea = 0.0;
182
183 switch (drivenDirection) {
184 case 'X':
185 if (info.xs == 0) { // Area is calculated by rank(s) on the negative face
186 i = 0;
187 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
188 if (nvert[k][j][i + 1] < 0.1)
189 localArea += sqrt(csi[k][j][i].x*csi[k][j][i].x + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].z*csi[k][j][i].z);
190 }
191 }
192 for (i = info.xs; i < lxe; i++) {
193 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
194 if (nvert[k][j][i + 1] < 0.1) localFluxProfile[i] += ucont[k][j][i].x;
195 }
196 }
197 break;
198 case 'Y':
199 if (info.ys == 0) {
200 j = 0;
201 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
202 if (nvert[k][j + 1][i] < 0.1)
203 localArea += sqrt(eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z);
204 }
205 }
206 for (j = info.ys; j < lye; j++) {
207 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
208 if (nvert[k][j + 1][i] < 0.1) localFluxProfile[j] += ucont[k][j][i].y;
209 }
210 }
211 break;
212 case 'Z':
213 if (info.zs == 0) {
214 k = 0;
215 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
216 if (nvert[k + 1][j][i] < 0.1)
217 localArea += sqrt(zet[k][j][i].x*zet[k][j][i].x + zet[k][j][i].y*zet[k][j][i].y + zet[k][j][i].z*zet[k][j][i].z);
218 }
219 }
220 for (k = info.zs; k < lze; k++) {
221 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
222 if (nvert[k + 1][j][i] < 0.1) localFluxProfile[k] += ucont[k][j][i].z;
223 }
224 }
225 break;
226 }
227
228 ierr = MPI_Allreduce(&localArea, &globalArea, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
229 ierr = MPI_Allreduce(localFluxProfile, globalFluxProfile, n_planes, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
230
231 // --- Step 6: Calculate Correction Profile ---
232 PetscReal targetFlux = simCtx->targetVolumetricFlux;
233 if (globalArea > 1.0e-12) {
234 for (i = 0; i < n_planes; i++) {
235 correctionProfile[i] = (targetFlux - globalFluxProfile[i]) / globalArea;
236 }
237 } else {
238 ierr = PetscMemzero(correctionProfile, n_planes * sizeof(PetscReal)); CHKERRQ(ierr);
239 }
240
241 LOG_ALLOW(GLOBAL, LOG_INFO, "Channel Flux Profile Corrector Update (Dir %c):\n", drivenDirection);
242 LOG_ALLOW(GLOBAL, LOG_INFO, " - Target Flux for all planes: %.6e\n", targetFlux);
243 LOG_ALLOW(GLOBAL, LOG_INFO, " - Measured Flux at plane 0: %.6e (Correction Velocity: %.6e)\n", globalFluxProfile[0], correctionProfile[0]);
244 LOG_ALLOW(GLOBAL, LOG_INFO, " - Measured Flux at plane %d: %.6e (Correction Velocity: %.6e)\n", (n_planes-1)/2, globalFluxProfile[(n_planes-1)/2], correctionProfile[(n_planes-1)/2]);
245
246 /* TURNED OFF IN LEGACY
247 // --- Step 7: Apply Correction to Velocity Profile ---
248 switch (drivenDirection) {
249 case 'X':
250 for (i = info.xs; i < info.xs + info.xm - 1; i++) {
251 if (PetscAbs(correctionProfile[i]) > 1e-12) {
252 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
253 if (nvert[k][j][i] < 0.1) {
254 PetscReal faceArea = sqrt(csi[k][j][i].x*csi[k][j][i].x + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].z*csi[k][j][i].z);
255 ucont[k][j][i].x += correctionProfile[i] * faceArea;
256 }
257 }
258 }
259 }
260 break;
261 case 'Y':
262 for (j = info.ys; j < info.ys + info.ym - 1; j++) {
263 if (PetscAbs(correctionProfile[j]) > 1e-12) {
264 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
265 if (nvert[k][j][i] < 0.1) {
266 PetscReal faceArea = sqrt(eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z);
267 ucont[k][j][i].y += correctionProfile[j] * faceArea;
268 }
269 }
270 }
271 }
272 break;
273 case 'Z':
274 for (k = info.zs; k < info.zs + info.zm - 1; k++) {
275 if (PetscAbs(correctionProfile[k]) > 1e-12) {
276 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
277 if (nvert[k][j][i] < 0.1) {
278 PetscReal faceArea = sqrt(zet[k][j][i].x*zet[k][j][i].x + zet[k][j][i].y*zet[k][j][i].y + zet[k][j][i].z*zet[k][j][i].z);
279 ucont[k][j][i].z += correctionProfile[k] * faceArea;
280 }
281 }
282 }
283 }
284 break;
285 }
286 */
287
288 // --- Step 8: Cleanup and Restore ---
289 ierr = PetscFree(localFluxProfile); CHKERRQ(ierr);
290 ierr = PetscFree(globalFluxProfile); CHKERRQ(ierr);
291 ierr = PetscFree(correctionProfile); CHKERRQ(ierr);
292
293 ierr = DMDAVecRestoreArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
294 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
295 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
296 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
297 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
298
299 //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Channel flux profile correction complete.\n",
300 // simCtx->rank, user->_this);
301
302 PetscFunctionReturn(0);
303}
304
305//Cell Center
306#define CP 0
307// Face Centers
308#define EP 1
309#define WP 2
310#define NP 3
311#define SP 4
312#define TP 5
313#define BP 6
314
315// Edge Centers
316#define NE 7
317#define SE 8
318#define NW 9
319#define SW 10
320#define TN 11
321#define BN 12
322#define TS 13
323#define BS 14
324#define TE 15
325#define BE 16
326#define TW 17
327#define BW 18
328
329#undef __FUNCT__
330#define __FUNCT__ "Projection"
331/**
332 * @brief Performs the projection step to enforce an incompressible velocity field.
333 *
334 * @details
335 * This function executes the final "correction" stage of a fractional-step (projection)
336 * method for solving the incompressible Navier-Stokes equations. After solving the
337 * Pressure Poisson Equation to get a pressure correction `Phi`, this function uses
338 * `Phi` to correct the intermediate velocity field, making it divergence-free.
339 *
340 * The main steps are:
341 * 1. **Calculate Pressure Gradient**: At each velocity point (on the cell faces), it
342 * computes the gradient of the pressure correction field (`∇Φ`). This is done using
343 * finite differences on a generalized curvilinear grid.
344 *
345 * 2. **Correct Velocity Field**: It updates the contravariant velocity components `Ucont`
346 * by subtracting the scaled pressure gradient. The update rule is:
347 * `U_new = U_intermediate - Δt * G * ∇Φ`, where `G` is the metric tensor `g_ij`
348 * that transforms the gradient into contravariant components.
349 *
350 * 3. **Boundary-Aware Gradients**: The gradient calculation is highly sensitive to
351 * boundaries. The code uses complex, dynamic stencils that switch from centered
352 * differences to one-sided differences if a standard stencil would use a point
353 * inside an immersed solid (`nvert > 0.1`). This preserves accuracy at the
354 * fluid-solid interface.
355 *
356 * 4. **Periodic Boundary Updates**: Includes dedicated loops to correctly calculate
357 * the pressure gradient at the domain edges for periodic boundary conditions.
358 *
359 * 5. **Finalization**: After correcting `Ucont`, it calls helper functions to convert
360 * the velocity back to Cartesian coordinates (`Contra2Cart`) and update all
361 * ghost cell values.
362 *
363 * @param[in, out] user Pointer to the UserCtx struct, containing simulation context,
364 * grid metrics, and the velocity/pressure fields.
365 *
366 * @return PetscErrorCode 0 on success, or a non-zero error code on failure.
367 */
368PetscErrorCode Projection(UserCtx *user)
369{
370 PetscFunctionBeginUser;
371 PROFILE_FUNCTION_BEGIN;
372 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering Projection step to correct velocity field.\n");
373
374 //================================================================================
375 // Section 1: Initialization and Data Acquisition
376 //================================================================================
377
378 // --- Get simulation and grid context ---
379 SimCtx *simCtx = user->simCtx;
380 DM da = user->da, fda = user->fda;
381 DMDALocalInfo info = user->info;
382
383 // --- Grid dimensions ---
384 PetscInt mx = info.mx, my = info.my, mz = info.mz;
385 PetscInt xs = info.xs, xe = info.xs + info.xm;
386 PetscInt ys = info.ys, ye = info.ys + info.ym;
387 PetscInt zs = info.zs, ze = info.zs + info.zm;
388
389 // --- Loop bounds (excluding outer ghost layers) ---
390 PetscInt lxs = (xs == 0) ? xs + 1 : xs;
391 PetscInt lxe = (xe == mx) ? xe - 1 : xe;
392 PetscInt lys = (ys == 0) ? ys + 1 : ys;
393 PetscInt lye = (ye == my) ? ye - 1 : ye;
394 PetscInt lzs = (zs == 0) ? zs + 1 : zs;
395 PetscInt lze = (ze == mz) ? ze - 1 : ze;
396
397 // --- Get direct pointer access to grid metric and field data ---
398 Cmpnts ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
399 PetscReal ***iaj, ***jaj, ***kaj, ***p, ***nvert;
400 Cmpnts ***ucont;
401 DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
402 DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
403 DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
404 DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
405 DMDAVecGetArray(da, user->lNvert, &nvert);
406 DMDAVecGetArray(da, user->lPhi, &p); // Note: using lPhi, which is the pressure correction
407 //DMDAVecGetArray(da,user->lP,&p);
408 DMDAVecGetArray(fda, user->Ucont, &ucont);
409
410 // --- Constants for clarity ---
411 const PetscReal IBM_FLUID_THRESHOLD = 0.1;
412 const PetscReal scale = simCtx->dt * 1.0 / COEF_TIME_ACCURACY; // simCtx->st replaced by 1.0.
413
414 LOG_ALLOW(GLOBAL,LOG_DEBUG," Starting velocity correction: Scale = %le .\n",scale);
415
416 //================================================================================
417 // Section 2: Correct Velocity Components
418 //================================================================================
419 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating pressure gradients and correcting velocity components.\n");
420
421 // --- Main loop over interior domain points ---
422 for (PetscInt k = lzs; k < lze; k++) {
423 for (PetscInt j = lys; j < lye; j++) {
424 for (PetscInt i = lxs; i < lxe; i++) {
425
426 // --- Correct U_contravariant (x-component of velocity) ---
427 PetscInt i_end = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) ? mx - 1 : mx - 2;
428 if (i < i_end) {
429
430 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j][i + 1] > IBM_FLUID_THRESHOLD)) {
431 // Compute pressure derivatives (dp/d_csi, dp/d_eta, dp/d_zet) at the i-face
432
433 PetscReal dpdc = p[k][j][i + 1] - p[k][j][i];
434 PetscReal dpde = 0.0, dpdz = 0.0;
435
436 // Boundary-aware stencil for dp/d_eta
437 if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
438 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
439 dpde = (p[k][j][i] + p[k][j][i+1] -
440 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
441 }
442 }
443
444 else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
445 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) { dpde = (p[k][j][i] + p[k][j][i+1] - p[k][j-1][i] - p[k][j-1][i+1]) * 0.5; }
446 }
447
448 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
449 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
450 }
451
452 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
453 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
454 }
455
456 else { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j-1][i] - p[k][j-1][i+1]) * 0.25; }
457
458 // Boundary-aware stencil for dp/d_zet
459 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
460 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) { dpdz = (p[k][j][i] + p[k][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.5; }
461 }
462
463 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
464 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) { dpdz = (p[k][j][i] + p[k][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.5; }
465 }
466
467 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
468 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
469 }
470
471 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
472 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
473 }
474
475 else { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.25; }
476
477 // Apply the correction: U_new = U_old - dt * (g11*dpdc + g12*dpde + g13*dpdz)
478
479
480
481 PetscReal grad_p_x = (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x + icsi[k][j][i].y * icsi[k][j][i].y
482 + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
483 dpde * (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y
484 + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
485 dpdz * (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y
486 + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i]);
487
488 PetscReal correction = grad_p_x*scale;
489 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Csi Direction: %le.\n",correction);
490 ucont[k][j][i].x -= correction;
491
492 }
493 }
494
495 // --- Correct V_contravariant (y-component of velocity) ---
496 PetscInt j_end = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) ? my - 1 : my - 2;
497 if (j < j_end) {
498 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j + 1][i] > IBM_FLUID_THRESHOLD)) {
499 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
500 dpde = p[k][j + 1][i] - p[k][j][i];
501
502 // Boundary-aware stencil for dp/d_csi
503 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
504 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) { dpdc = (p[k][j][i] + p[k][j+1][i] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.5; }
505 } else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
506 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) { dpdc = (p[k][j][i] + p[k][j+1][i] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.5; }
507 } else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
508 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
509 } else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
510 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
511 } else { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.25; }
512
513 // Boundary-aware stencil for dp/d_zet
514 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
515 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) { dpdz = (p[k][j][i] + p[k][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.5; }
516 } else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
517 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) { dpdz = (p[k][j][i] + p[k][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.5; }
518 } else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
519 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
520 } else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
521 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
522 } else { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.25; }
523
524 PetscReal grad_p_y = (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x + jcsi[k][j][i].y * jeta[k][j][i].y + jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
525 dpde * (jeta[k][j][i].x * jeta[k][j][i].x + jeta[k][j][i].y * jeta[k][j][i].y + jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
526 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x + jzet[k][j][i].y * jeta[k][j][i].y + jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i]);
527
528 PetscReal correction = grad_p_y*scale;
529 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Eta Direction: %le.\n",correction);
530 ucont[k][j][i].y -= correction;
531 }
532 }
533
534 // --- Correct W_contravariant (z-component of velocity) ---
535 PetscInt k_end = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) ? mz - 1 : mz - 2;
536 if (k < k_end) {
537 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k + 1][j][i] > IBM_FLUID_THRESHOLD)) {
538 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
539 dpdz = p[k + 1][j][i] - p[k][j][i];
540
541 // Boundary-aware stencil for dp/d_csi
542 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
543 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) { dpdc = (p[k][j][i] + p[k+1][j][i] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.5; }
544 } else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
545 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) { dpdc = (p[k][j][i] + p[k+1][j][i] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.5; }
546 } else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
547 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
548 } else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
549 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
550 } else { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.25; }
551
552 // Boundary-aware stencil for dp/d_eta
553 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
554 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) { dpde = (p[k][j][i] + p[k+1][j][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.5; }
555 } else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
556 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) { dpde = (p[k][j][i] + p[k+1][j][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.5; }
557 } else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
558 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
559 } else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
560 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
561 } else { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.25; }
562
563 PetscReal grad_p_z = (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x + kcsi[k][j][i].y * kzet[k][j][i].y + kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
564 dpde * (keta[k][j][i].x * kzet[k][j][i].x + keta[k][j][i].y * kzet[k][j][i].y + keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
565 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x + kzet[k][j][i].y * kzet[k][j][i].y + kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i]);
566
567 // ========================= DEBUG PRINT =========================
568 LOG_LOOP_ALLOW_EXACT(GLOBAL, LOG_DEBUG, k, 5,
569 "[k=%d, j=%d, i=%d] ---- Neighbor Pressures ----\n"
570 " Central Z-Neighbors: p[k+1][j][i] = %g | p[k][j][i] = %g\n"
571 " Eta-Stencil (Y-dir): p[k][j-1][i] = %g, p[k+1][j-1][i] = %g | p[k][j+1][i] = %g, p[k+1][j+1][i] = %g\n"
572 " Csi-Stencil (X-dir): p[k][j][i-1] = %g, p[k+1][j][i-1] = %g | p[k][j][i+1] = %g, p[k+1][j][i+1] = %g\n",
573 k, j, i,
574 p[k + 1][j][i], p[k][j][i],
575 p[k][j - 1][i], p[k + 1][j - 1][i], p[k][j + 1][i], p[k + 1][j + 1][i],
576 p[k][j][i - 1], p[k + 1][j][i - 1], p[k][j][i + 1], p[k + 1][j][i + 1]);
577 // ======================= END DEBUG PRINT =======================
578
579 LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," dpdc: %le | dpde: %le | dpdz: %le.\n",dpdc,dpde,dpdz);
580 PetscReal correction = grad_p_z*scale;
581 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Zet Direction: %le.\n",correction);
582 ucont[k][j][i].z -= correction;
583 }
584 }
585 }
586 }
587 }
588
589 // --- Explicit correction for periodic boundaries (if necessary) ---
590 // The main loop handles the interior, but this handles the first physical layer at periodic boundaries.
591 // Note: This logic is largely duplicated from the main loop and could be merged, but is preserved for fidelity.
592 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
593 for (PetscInt k=lzs; k<lze; k++) {
594 for (PetscInt j=lys; j<lye; j++) {
595 PetscInt i=xs;
596
597 PetscReal dpdc = p[k][j][i+1] - p[k][j][i];
598
599 PetscReal dpde = 0.;
600 PetscReal dpdz = 0.;
601
602 if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
603 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
604 dpde = (p[k][j ][i] + p[k][j ][i+1] -
605 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
606 }
607 }
608 else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
609 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
610 dpde = (p[k][j ][i] + p[k][j ][i+1] -
611 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
612 }
613 }
614 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
615 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
616 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
617 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
618 }
619 }
620 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
621 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
622 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
623 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
624 }
625 }
626 else {
627 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
628 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
629 }
630
631 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
632 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
633 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
634 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
635 }
636 }
637 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
638 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
639 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
640 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
641 }
642 }
643 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
644 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
645 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
646 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
647 }
648 }
649 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
650 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
651 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
652 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
653 }
654 }
655 else {
656 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
657 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
658 }
659
660
661
662 if (!(nvert[k][j][i] + nvert[k][j][i+1])) {
663 ucont[k][j][i].x -=
664 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
665 icsi[k][j][i].y * icsi[k][j][i].y +
666 icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
667 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
668 ieta[k][j][i].y * icsi[k][j][i].y +
669 ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
670 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
671 izet[k][j][i].y * icsi[k][j][i].y +
672 izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i])
673 * scale;
674
675 }
676 }
677 }
678 }
679 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
680
681 for (PetscInt k=lzs; k<lze; k++) {
682 for (PetscInt i=lxs; i<lxe; i++) {
683 PetscInt j=ys;
684
685 PetscReal dpdc = 0.;
686 PetscReal dpdz = 0.;
687 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
688 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
689 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
690 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
691 }
692 }
693 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
694 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
695 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
696 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
697 }
698 }
699 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
700 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
701 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
702 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
703 }
704 }
705 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
706 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
707 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
708 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
709 }
710 }
711 else {
712 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
713 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
714 }
715
716 PetscReal dpde = p[k][j+1][i] - p[k][j][i];
717
718 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
719 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
720 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
721 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
722 }
723 }
724 else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
725 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
726 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
727 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
728 }
729 }
730 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
731 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
732 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
733 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
734 }
735 }
736 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
737 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
738 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
739 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
740 }
741 }
742 else {
743 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
744 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
745 }
746
747 if (!(nvert[k][j][i] + nvert[k][j+1][i])) {
748 ucont[k][j][i].y -=
749 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
750 jcsi[k][j][i].y * jeta[k][j][i].y +
751 jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
752 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
753 jeta[k][j][i].y * jeta[k][j][i].y +
754 jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
755 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
756 jzet[k][j][i].y * jeta[k][j][i].y +
757 jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i])
758 * scale;
759 }
760 }
761 }
762 }
763
764 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
765 for (PetscInt j=lys; j<lye; j++) {
766 for (PetscInt i=lxs; i<lxe; i++) {
767
768 PetscInt k=zs;
769 PetscReal dpdc = 0.;
770 PetscReal dpde = 0.;
771
772 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
773 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
774 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
775 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
776 }
777 }
778 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
779 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
780 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
781 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
782 }
783 }
784 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
785 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
786 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
787 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
788 }
789 }
790 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
791 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
792 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
793 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
794 }
795 }
796 else {
797 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
798 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
799 }
800
801 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
802 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
803 dpde = (p[k][j ][i] + p[k+1][j ][i] -
804 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
805 }
806 }
807 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
808 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
809 dpde = (p[k][j ][i] + p[k+1][j ][i] -
810 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
811 }
812 }
813 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
814 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
815 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
816 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
817 }
818 }
819 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
820 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
821 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
822 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
823 }
824 }
825 else {
826 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
827 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
828 }
829
830 PetscReal dpdz = p[k+1][j][i] - p[k][j][i];
831
832 if (!(nvert[k][j][i] + nvert[k+1][j][i])) {
833
834 ucont[k][j][i].z -=
835 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
836 kcsi[k][j][i].y * kzet[k][j][i].y +
837 kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
838 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
839 keta[k][j][i].y * kzet[k][j][i].y +
840 keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
841 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
842 kzet[k][j][i].y * kzet[k][j][i].y +
843 kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i])
844 * scale;
845
846 }
847 }
848 }
849 }
850
851 // Corrects Flux Profile for Driven Flows if applicable.
852 CorrectChannelFluxProfile(user);
853
854 //================================================================================
855 // Section 3: Finalization and Cleanup
856 //================================================================================
857
858 // --- Restore access to all PETSc vector arrays ---
859 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
860 // DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
861 //DMDAVecRestoreArray(da, user->lAj, &aj);
862 DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
863 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
864 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
865 DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
866 DMDAVecRestoreArray(da, user->lP, &p);
867 DMDAVecRestoreArray(da, user->lNvert, &nvert);
868
869 // --- Update ghost cells for the newly corrected velocity field ---
870 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating ghost cells for corrected velocity.\n");
871 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
872 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
873
874 // --- Convert velocity to Cartesian and update ghost nodes ---
875 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Converting velocity to Cartesian and finalizing ghost nodes.\n");
876 Contra2Cart(user);
877 UpdateLocalGhosts(user,"Ucat");
878 RefreshBoundaryGhostCells(user);
879 //GhostNodeVelocity(user);
880
881 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting Projection step.\n");
882 PROFILE_FUNCTION_END;
883 PetscFunctionReturn(0);
884}
885
886#undef __FUNCT__
887#define __FUNCT__ "UpdatePressure"
888/**
889 * @brief Updates the pressure field and handles periodic boundary conditions.
890 *
891 * @details
892 * This function is a core part of the projection method in a CFD solver. It is
893 * called after the Pressure Poisson Equation has been solved to obtain a pressure
894 * correction field, `Phi`.
895 *
896 * The function performs two main operations:
897 *
898 * 1. **Core Pressure Update**: It updates the main pressure field `P` by adding the
899 * pressure correction `Phi` at every grid point in the fluid domain. The
900 * operation is `P_new = P_old + Phi`.
901 *
902 * 2. **Periodic Boundary Synchronization**: If any of the domain boundaries are
903 * periodic (`bctype == 7`), this function manually updates the pressure values
904 * in the ghost cells. It copies values from the physical cells on one side of
905 * the domain to the corresponding ghost cells on the opposite side. This ensures
906 * that subsequent calculations involving pressure gradients are correct across
907 * periodic boundaries. The refactored version consolidates the original code's
908 * redundant updates into a single, efficient pass.
909 *
910 * @param[in, out] user Pointer to the UserCtx struct, containing simulation context and
911 * the pressure vectors `P` and `Phi`.
912 *
913 * @return PetscErrorCode 0 on success, or a non-zero error code on failure.
914 */
915PetscErrorCode UpdatePressure(UserCtx *user)
916{
917 PetscFunctionBeginUser;
918 PROFILE_FUNCTION_BEGIN;
919 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering UpdatePressure.\n");
920
921 //================================================================================
922 // Section 1: Initialization and Data Acquisition
923 //================================================================================
924 DM da = user->da;
925 DMDALocalInfo info = user->info;
926
927 // Local grid extents for the main update loop
928 PetscInt xs = info.xs, xe = info.xs + info.xm;
929 PetscInt ys = info.ys, ye = info.ys + info.ym;
930 PetscInt zs = info.zs, ze = info.zs + info.zm;
931 PetscInt mx = info.mx, my = info.my, mz = info.mz;
932
933 PetscInt i,j,k;
934
935 // --- Get direct pointer access to PETSc vector data for performance ---
936 PetscReal ***p, ***phi;
937 DMDAVecGetArray(da, user->P, &p);
938 DMDAVecGetArray(da, user->Phi, &phi);
939
940 //================================================================================
941 // Section 2: Core Pressure Update
942 //================================================================================
943 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Performing core pressure update (P_new = P_old + Phi).\n");
944 for (PetscInt k = zs; k < ze; k++) {
945 for (PetscInt j = ys; j < ye; j++) {
946 for (PetscInt i = xs; i < xe; i++) {
947 // This is the fundamental pressure update in a projection method.
948 p[k][j][i] += phi[k][j][i];
949 }
950 }
951 }
952
953 // Restore arrays now that the core computation is done.
954 DMDAVecRestoreArray(da, user->Phi, &phi);
955 DMDAVecRestoreArray(da, user->P, &p);
956
957
958 //================================================================================
959 // Section 3: Handle Periodic Boundary Condition Synchronization
960 //================================================================================
961 // This block is executed only if at least one boundary is periodic.
962 // The original code contained many redundant Get/Restore and update calls.
963 // This refactored version performs the same effective logic but in a single,
964 // efficient pass, which is numerically equivalent and much cleaner.
965 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC ||
966 user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC ||
967 user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC)
968 {
969 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Synchronizing ghost cells for periodic boundaries.\n");
970
971 // First, update the local vectors (including ghost regions) with the new global data.
972 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
973 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
974 DMGlobalToLocalBegin(da, user->Phi, INSERT_VALUES, user->lPhi);
975 DMGlobalToLocalEnd(da, user->Phi, INSERT_VALUES, user->lPhi);
976
977 // Get pointers to all necessary local and global arrays ONCE.
978 PetscReal ***lp, ***lphi;
979 DMDAVecGetArray(da, user->lP, &lp);
980 DMDAVecGetArray(da, user->lPhi, &lphi);
981 DMDAVecGetArray(da, user->P, &p);
982 DMDAVecGetArray(da, user->Phi, &phi);
983
984 // --- X-Direction Periodic Update ---
985 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
986 // Update left boundary physical cells from right boundary ghost cells
987 if (xs == 0) {
988 PetscInt i = 0;
989 for (k = zs; k < ze; k++) {
990 for (j = ys; j < ye; j++) {
991 // Note: Accessing lp[...][i-2] reads from the ghost cell region.
992 p[k][j][i] = lp[k][j][i - 2];
993 phi[k][j][i] = lphi[k][j][i - 2];
994 }
995 }
996 }
997 // Update right boundary physical cells from left boundary ghost cells
998 if (xe == mx) {
999 PetscInt i = mx - 1;
1000 for (k = zs; k < ze; k++) {
1001 for (j = ys; j < ye; j++) {
1002 p[k][j][i] = lp[k][j][i + 2];
1003 phi[k][j][i] = lphi[k][j][i + 2];
1004 }
1005 }
1006 }
1007 }
1008
1009 // --- Y-Direction Periodic Update ---
1010 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
1011 // Update bottom boundary
1012 if (ys == 0) {
1013 PetscInt j = 0;
1014 for (k = zs; k < ze; k++) {
1015 for (i = xs; i < xe; i++) {
1016 p[k][j][i] = lp[k][j - 2][i];
1017 phi[k][j][i] = lphi[k][j - 2][i];
1018 }
1019 }
1020 }
1021 // Update top boundary
1022 if (ye == my) {
1023 PetscInt j = my - 1;
1024 for (k = zs; k < ze; k++) {
1025 for (i = xs; i < xe; i++) {
1026 p[k][j][i] = lp[k][j + 2][i];
1027 phi[k][j][i] = lphi[k][j + 2][i];
1028 }
1029 }
1030 }
1031 }
1032
1033 // --- Z-Direction Periodic Update ---
1034 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
1035 // Update front boundary
1036 if (zs == 0) {
1037 PetscInt k = 0;
1038 for (j = ys; j < ye; j++) {
1039 for (i = xs; i < xe; i++) {
1040 p[k][j][i] = lp[k - 2][j][i];
1041 phi[k][j][i] = lphi[k - 2][j][i];
1042 }
1043 }
1044 }
1045 // Update back boundary
1046 if (ze == mz) {
1047 PetscInt k = mz - 1;
1048 for (j = ys; j < ye; j++) {
1049 for (i = xs; i < xe; i++) {
1050 p[k][j][i] = lp[k + 2][j][i];
1051 phi[k][j][i] = lphi[k + 2][j][i];
1052 }
1053 }
1054 }
1055 }
1056
1057 // Restore all arrays ONCE.
1058 DMDAVecRestoreArray(da, user->lP, &lp);
1059 DMDAVecRestoreArray(da, user->lPhi, &lphi);
1060 DMDAVecRestoreArray(da, user->P, &p);
1061 DMDAVecRestoreArray(da, user->Phi, &phi);
1062
1063 // After manually updating the physical boundary cells, we must call
1064 // DMGlobalToLocal again to ensure all processes have the updated ghost
1065 // values for the *next* function that needs them.
1066 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
1067 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
1068 DMGlobalToLocalBegin(da, user->Phi, INSERT_VALUES, user->lPhi);
1069 DMGlobalToLocalEnd(da, user->Phi, INSERT_VALUES, user->lPhi);
1070 }
1071
1072 //================================================================================
1073 // Section 4: Final Cleanup (pointers already restored)
1074 //================================================================================
1075
1076 UpdateLocalGhosts(user,"P");
1077 UpdateLocalGhosts(user,"Phi");
1078
1079 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting UpdatePressure.\n");
1080 PROFILE_FUNCTION_END;
1081 PetscFunctionReturn(0);
1082}
1083
1084
1085
1086static PetscErrorCode mymatmultadd(Mat mat, Vec v1, Vec v2)
1087{
1088 Vec vt;
1089 VecDuplicate(v2, &vt);
1090 MatMult(mat, v1, vt);
1091 VecAYPX(v2, 1., vt);
1092 VecDestroy(&vt);
1093 return(0);
1094}
1095
1096
1097#undef __FUNCT__
1098#define __FUNCT__ "PoissonNullSpaceFunction"
1099PetscErrorCode PoissonNullSpaceFunction(MatNullSpace nullsp,Vec X, void *ctx)
1100{
1101 UserCtx *user = (UserCtx*)ctx;
1102
1103 DM da = user->da;
1104
1105 DMDALocalInfo info = user->info;
1106 PetscInt xs = info.xs, xe = info.xs + info.xm;
1107 PetscInt ys = info.ys, ye = info.ys + info.ym;
1108 PetscInt zs = info.zs, ze = info.zs + info.zm;
1109 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1110 PetscInt lxs, lxe, lys, lye, lzs, lze;
1111
1112 PetscReal ***x, ***nvert;
1113 PetscInt i, j, k;
1114
1115/* /\* First remove a constant from the Vec field X *\/ */
1116
1117
1118 /* Then apply boundary conditions */
1119 DMDAVecGetArray(da, X, &x);
1120 DMDAVecGetArray(da, user->lNvert, &nvert);
1121
1122 lxs = xs; lxe = xe;
1123 lys = ys; lye = ye;
1124 lzs = zs; lze = ze;
1125
1126 if (xs==0) lxs = xs+1;
1127 if (ys==0) lys = ys+1;
1128 if (zs==0) lzs = zs+1;
1129
1130 if (xe==mx) lxe = xe-1;
1131 if (ye==my) lye = ye-1;
1132 if (ze==mz) lze = ze-1;
1133
1134 PetscReal lsum, sum;
1135 PetscReal lnum, num;
1136
1137 if (user->multinullspace) PetscPrintf(PETSC_COMM_WORLD, "MultiNullSpace!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1138 if (!user->multinullspace) {
1139 lsum = 0;
1140 lnum = 0;
1141 for (k=lzs; k<lze; k++) {
1142 for (j=lys; j<lye; j++) {
1143 for (i=lxs; i<lxe; i++) {
1144 if (nvert[k][j][i] < 0.1) {
1145 lsum += x[k][j][i];
1146 lnum ++;
1147 }
1148 }
1149 }
1150 }
1151
1152 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1153 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1154 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1155/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1156 sum = sum / (-1.0 * num);
1157
1158 for (k=lzs; k<lze; k++) {
1159 for (j=lys; j<lye; j++) {
1160 for (i=lxs; i<lxe; i++) {
1161 if (nvert[k][j][i] < 0.1) {
1162 x[k][j][i] +=sum;
1163 }
1164 }
1165 }
1166 }
1167 }
1168 else {
1169 lsum = 0;
1170 lnum = 0;
1171 for (j=lys; j<lye; j++) {
1172 for (i=lxs; i<lxe; i++) {
1173 for (k=lzs; k<lze; k++) {
1174 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1175 lsum += x[k][j][i];
1176 lnum ++;
1177 }
1178 }
1179 }
1180 }
1181 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1182 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1183 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1184/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1185 sum /= -num;
1186 for (j=lys; j<lye; j++) {
1187 for (i=lxs; i<lxe; i++) {
1188 for (k=lzs; k<lze; k++) {
1189 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1190 x[k][j][i] += sum;
1191 }
1192 }
1193 }
1194 }
1195
1196 lsum = 0;
1197 lnum = 0;
1198 for (j=lys; j<lye; j++) {
1199 for (i=lxs; i<lxe; i++) {
1200 for (k=lzs; k<lze; k++) {
1201 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1202 lsum += x[k][j][i];
1203 lnum ++;
1204 }
1205 }
1206 }
1207 }
1208 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1209 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1210 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1211/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1212 sum /= -num;
1213 for (j=lys; j<lye; j++) {
1214 for (i=lxs; i<lxe; i++) {
1215 for (k=lzs; k<lze; k++) {
1216 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1217 x[k][j][i] += sum;
1218 }
1219 }
1220 }
1221 }
1222
1223 } //if multinullspace
1224 if (zs == 0) {
1225 k = 0;
1226 for (j=ys; j<ye; j++) {
1227 for (i=xs; i<xe; i++) {
1228 x[k][j][i] = 0.;
1229 }
1230 }
1231 }
1232
1233 if (ze == mz) {
1234 k = mz-1;
1235 for (j=ys; j<ye; j++) {
1236 for (i=xs; i<xe; i++) {
1237 x[k][j][i] = 0.;
1238 }
1239 }
1240 }
1241
1242 if (ys == 0) {
1243 j = 0;
1244 for (k=zs; k<ze; k++) {
1245 for (i=xs; i<xe; i++) {
1246 x[k][j][i] = 0.;
1247 }
1248 }
1249 }
1250
1251 if (ye == my) {
1252 j = my-1;
1253 for (k=zs; k<ze; k++) {
1254 for (i=xs; i<xe; i++) {
1255 x[k][j][i] = 0.;
1256 }
1257 }
1258 }
1259
1260 if (xs == 0) {
1261 i = 0;
1262 for (k=zs; k<ze; k++) {
1263 for (j=ys; j<ye; j++) {
1264 x[k][j][i] = 0.;
1265 }
1266 }
1267 }
1268
1269 if (xe == mx) {
1270 i = mx-1;
1271 for (k=zs; k<ze; k++) {
1272 for (j=ys; j<ye; j++) {
1273 x[k][j][i] = 0.;
1274 }
1275 }
1276 }
1277
1278 for (k=zs; k<ze; k++) {
1279 for (j=ys; j<ye; j++) {
1280 for (i=xs; i<xe; i++) {
1281 if (nvert[k][j][i] > 0.1)
1282 x[k][j][i] = 0.;
1283 }
1284 }
1285 }
1286 DMDAVecRestoreArray(da, X, &x);
1287 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1288
1289 return 0;
1290}
1291
1292PetscErrorCode MyInterpolation(Mat A, Vec X, Vec F)
1293{
1294 UserCtx *user;
1295
1296 MatShellGetContext(A, (void**)&user);
1297
1298
1299
1300 DM da = user->da;
1301
1302 DM da_c = *user->da_c;
1303
1304 DMDALocalInfo info = user->info;
1305 PetscInt xs = info.xs, xe = info.xs + info.xm;
1306 PetscInt ys = info.ys, ye = info.ys + info.ym;
1307 PetscInt zs = info.zs, ze = info.zs + info.zm;
1308 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1309 PetscInt lxs, lxe, lys, lye, lzs, lze;
1310
1311 PetscReal ***f, ***x, ***nvert, ***nvert_c;
1312 PetscInt i, j, k, ic, jc, kc, ia, ja, ka;
1313
1314 lxs = xs; lxe = xe;
1315 lys = ys; lye = ye;
1316 lzs = zs; lze = ze;
1317
1318 if (xs==0) lxs = xs+1;
1319 if (ys==0) lys = ys+1;
1320 if (zs==0) lzs = zs+1;
1321
1322 if (xe==mx) lxe = xe-1;
1323 if (ye==my) lye = ye-1;
1324 if (ze==mz) lze = ze-1;
1325
1326
1327 DMDAVecGetArray(da, F, &f);
1328
1329
1330 Vec lX;
1331 DMCreateLocalVector(da_c, &lX);
1332
1333 DMGlobalToLocalBegin(da_c, X, INSERT_VALUES, lX);
1334 DMGlobalToLocalEnd(da_c, X, INSERT_VALUES, lX);
1335 DMDAVecGetArray(da_c, lX, &x);
1336
1337 DMDAVecGetArray(da, user->lNvert, &nvert);
1338 DMDAVecGetArray(da_c, *(user->lNvert_c), &nvert_c);
1339 for (k=lzs; k<lze; k++) {
1340 for (j=lys; j<lye; j++) {
1341 for (i=lxs; i<lxe; i++) {
1342
1343 GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user);
1344
1345 f[k][j][i] = (x[kc ][jc ][ic ] * 9 +
1346 x[kc ][jc+ja][ic ] * 3 +
1347 x[kc ][jc ][ic+ia] * 3 +
1348 x[kc ][jc+ja][ic+ia]) * 3./64. +
1349 (x[kc+ka][jc ][ic ] * 9 +
1350 x[kc+ka][jc+ja][ic ] * 3 +
1351 x[kc+ka][jc ][ic+ia] * 3 +
1352 x[kc+ka][jc+ja][ic+ia]) /64.;
1353 }
1354 }
1355 }
1356
1357 for (k=zs; k<ze; k++) {
1358 for (j=ys; j<ye; j++) {
1359 for (i=xs; i<xe; i++) {
1360
1361 if (i==0) {
1362 f[k][j][i] = 0.;//-f[k][j][i+1];
1363 }
1364 else if (i==mx-1) {
1365 f[k][j][i] = 0.;//-f[k][j][i-1];
1366 }
1367 else if (j==0) {
1368 f[k][j][i] = 0.;//-f[k][j+1][i];
1369 }
1370 else if (j==my-1) {
1371 f[k][j][i] = 0.;//-f[k][j-1][i];
1372 }
1373 else if (k==0) {
1374 f[k][j][i] = 0.;//-f[k+1][j][i];
1375 }
1376 else if (k==mz-1) {
1377 f[k][j][i] = 0.;//-f[k-1][j][i];
1378 }
1379 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
1380
1381 }
1382 }
1383 }
1384
1385 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1386 DMDAVecRestoreArray(da_c, *(user->lNvert_c), &nvert_c);
1387
1388 DMDAVecRestoreArray(da_c, lX, &x);
1389
1390 VecDestroy(&lX);
1391 DMDAVecRestoreArray(da, F, &f);
1392
1393
1394
1395 return 0;
1396
1397}
1398
1399static PetscErrorCode RestrictResidual_SolidAware(Mat A, Vec X, Vec F)
1400{
1401 UserCtx *user;
1402 MatShellGetContext(A, (void**)&user);
1403
1404 DM da = user->da;
1405 DM da_f = *user->da_f;
1406
1407 DMDALocalInfo info;
1408 DMDAGetLocalInfo(da, &info);
1409 PetscInt xs = info.xs, xe = info.xs + info.xm;
1410 PetscInt ys = info.ys, ye = info.ys + info.ym;
1411 PetscInt zs = info.zs, ze = info.zs + info.zm;
1412 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1413
1414 PetscReal ***f, ***x, ***nvert;
1415 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
1416
1417 DMDAVecGetArray(da, F, &f);
1418
1419 Vec lX;
1420 DMCreateLocalVector(da_f, &lX);
1421 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
1422 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
1423 DMDAVecGetArray(da_f, lX, &x);
1424
1425 DMDAVecGetArray(da, user->lNvert, &nvert);
1426
1427 PetscReal ***nvert_f;
1428 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
1429
1430 if ((user->isc)) ia = 0;
1431 else ia = 1;
1432
1433 if ((user->jsc)) ja = 0;
1434 else ja = 1;
1435
1436 if ((user->ksc)) ka = 0;
1437 else ka = 1;
1438
1439 for (k=zs; k<ze; k++) {
1440 for (j=ys; j<ye; j++) {
1441 for (i=xs; i<xe; i++) {
1442 // --- CORRECTED LOGIC ---
1443 // First, check if the current point is a boundary point.
1444 // If it is, it does not contribute to the coarse grid residual.
1445 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1 || nvert[k][j][i] > 0.1) {
1446 f[k][j][i] = 0.0;
1447 }
1448 // Only if it's a true interior fluid point, perform the restriction.
1449 else {
1450 GridRestriction(i, j, k, &ih, &jh, &kh, user);
1451 f[k][j][i] = 0.125 *
1452 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
1453 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
1454 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
1455 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
1456 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
1457 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
1458 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
1459 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
1460 }
1461 }
1462 }
1463 }
1464
1465 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
1466 DMDAVecRestoreArray(da_f, lX, &x);
1467 VecDestroy(&lX);
1468 DMDAVecRestoreArray(da, F, &f);
1469 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1470
1471 return 0;
1472}
1473
1474PetscErrorCode MyRestriction(Mat A, Vec X, Vec F)
1475{
1476 UserCtx *user;
1477
1478 MatShellGetContext(A, (void**)&user);
1479
1480
1481 DM da = user->da;
1482
1483 DM da_f = *user->da_f;
1484
1485 DMDALocalInfo info;
1486 DMDAGetLocalInfo(da, &info);
1487 PetscInt xs = info.xs, xe = info.xs + info.xm;
1488 PetscInt ys = info.ys, ye = info.ys + info.ym;
1489 PetscInt zs = info.zs, ze = info.zs + info.zm;
1490 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1491 // PetscInt lxs, lxe, lys, lye, lzs, lze;
1492
1493 PetscReal ***f, ***x, ***nvert;
1494 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
1495
1496 DMDAVecGetArray(da, F, &f);
1497
1498 Vec lX;
1499
1500 DMCreateLocalVector(da_f, &lX);
1501 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
1502 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
1503 DMDAVecGetArray(da_f, lX, &x);
1504
1505 DMDAVecGetArray(da, user->lNvert, &nvert);
1506
1507 PetscReal ***nvert_f;
1508 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
1509
1510 if ((user->isc)) ia = 0;
1511 else ia = 1;
1512
1513 if ((user->jsc)) ja = 0;
1514 else ja = 1;
1515
1516 if ((user->ksc)) ka = 0;
1517 else ka = 1;
1518
1519 for (k=zs; k<ze; k++) {
1520 for (j=ys; j<ye; j++) {
1521 for (i=xs; i<xe; i++) {
1522 if (k==0) {
1523 f[k][j][i] = 0.;
1524 }
1525 else if (k==mz-1) {
1526 f[k][j][i] = 0.;
1527 }
1528 else if (j==0) {
1529 f[k][j][i] = 0.;
1530 }
1531 else if (j==my-1) {
1532 f[k][j][i] = 0.;
1533 }
1534 else if (i==0) {
1535 f[k][j][i] = 0.;
1536 }
1537 else if (i==mx-1) {
1538 f[k][j][i] = 0.;
1539 }
1540 else {
1541 GridRestriction(i, j, k, &ih, &jh, &kh, user);
1542 f[k][j][i] = 0.125 *
1543 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
1544 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
1545 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
1546 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
1547 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
1548 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
1549 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
1550 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
1551
1552
1553
1554 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
1555 }
1556 }
1557 }
1558 }
1559
1560
1561 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
1562
1563 DMDAVecRestoreArray(da_f, lX, &x);
1564 VecDestroy(&lX);
1565
1566 DMDAVecRestoreArray(da, F, &f);
1567 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1568
1569
1570 return 0;
1571}
1572
1573#undef __FUNCT__
1574#define __FUNCT__ "PoissonLHSNew"
1575
1576/**
1577 * @brief Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.
1578 *
1579 * @details
1580 * This function constructs the sparse matrix `A` (the LHS) for the linear system `Ax = B`,
1581 * which is the Pressure Poisson Equation (PPE). The matrix `A` represents a discrete
1582 * version of the negative Laplacian operator (-∇²), tailored for a general curvilinear,
1583 * staggered grid.
1584 *
1585 * The assembly process is highly complex and follows these main steps:
1586 *
1587 * 1. **Matrix Initialization**: On the first call, it allocates a sparse PETSc matrix `A`
1588 * pre-allocating space for a 19-point stencil per row. On subsequent calls, it
1589 * simply zeroes out the existing matrix.
1590 *
1591 * 2. **Metric Tensor Calculation**: It first computes the 9 components of the metric
1592 * tensor (`g11`, `g12`, ..., `g33`) at the cell faces. These `g_ij` coefficients
1593 * are essential for defining the Laplacian operator in generalized curvilinear
1594 * coordinates and account for grid stretching and non-orthogonality.
1595 *
1596 * 3. **Stencil Assembly Loop**: The function iterates through every grid point `(i,j,k)`.
1597 * For each point, it determines the matrix row entries based on its status:
1598 * - **Boundary/Solid Point**: If the point is on a domain boundary or inside an
1599 * immersed solid (`nvert > 0.1`), it sets an identity row (`A(row,row) = 1`),
1600 * effectively removing it from the linear solve.
1601 * - **Fluid Point**: For a fluid point, it computes the 19 coefficients of the
1602 * finite volume stencil. This involves summing contributions from each of the
1603 * six faces of the control volume around the point.
1604 *
1605 * 4. **Boundary-Aware Stencils**: The stencil calculation is critically dependent on the
1606 * state of neighboring cells. The code contains intricate logic to check if neighbors
1607 * are fluid or solid. If a standard centered-difference stencil would cross into a
1608 * solid, the scheme is automatically adapted to a one-sided difference to maintain
1609 * accuracy at the fluid-solid interface.
1610 *
1611 * 5. **Matrix Finalization**: After all rows corresponding to the local processor's
1612 * grid points have been set, `MatAssemblyBegin/End` is called to finalize the
1613 * global matrix, making it ready for use by a linear solver.
1614 *
1615 * @param[in, out] user Pointer to the UserCtx struct, containing all simulation context,
1616 * grid data, and the matrix `A`.
1617 *
1618 * @return PetscErrorCode 0 on success, or a non-zero error code on failure.
1619 */
1620PetscErrorCode PoissonLHSNew(UserCtx *user)
1621{
1622 PetscFunctionBeginUser;
1623 PROFILE_FUNCTION_BEGIN;
1624 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonLHSNew to assemble Laplacian matrix.\n");
1625 PetscErrorCode ierr;
1626 //================================================================================
1627 // Section 1: Initialization and Data Acquisition
1628 //================================================================================
1629
1630
1631 // --- Get simulation and grid context ---
1632 SimCtx *simCtx = user->simCtx;
1633 DM da = user->da, fda = user->fda;
1634 DMDALocalInfo info = user->info;
1635 PetscInt IM = user->IM, JM = user->JM, KM = user->KM;
1636 PetscInt i,j,k;
1637
1638 // --- Grid dimensions ---
1639 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1640 PetscInt xs = info.xs, xe = info.xs + info.xm;
1641 PetscInt ys = info.ys, ye = info.ys + info.ym;
1642 PetscInt zs = info.zs, ze = info.zs + info.zm;
1643 PetscInt gxs = info.gxs, gxe = gxs + info.gxm;
1644 PetscInt gys = info.gys, gye = gys + info.gym;
1645 PetscInt gzs = info.gzs, gze = gzs + info.gzm;
1646
1647 // --- Define constants for clarity ---
1648 const PetscReal IBM_FLUID_THRESHOLD = 0.1;
1649
1650 // --- Allocate the LHS matrix A on the first call ---
1651 if (!user->assignedA) {
1652 LOG_ALLOW(GLOBAL, LOG_INFO, "First call: Creating LHS matrix 'A' with 19-point stencil preallocation.\n");
1653 PetscInt N = mx * my * mz; // Total size
1654 PetscInt M; // Local size
1655 VecGetLocalSize(user->Phi, &M);
1656 // Create a sparse AIJ matrix, preallocating for 19 non-zeros per row (d=diagonal, o=off-diagonal)
1657 MatCreateAIJ(PETSC_COMM_WORLD, M, M, N, N, 19, PETSC_NULL, 19, PETSC_NULL, &(user->A));
1658 user->assignedA = PETSC_TRUE;
1659 }
1660
1661 // Zero out matrix entries from the previous solve
1662 MatZeroEntries(user->A);
1663
1664 // --- Get direct pointer access to grid metric data ---
1665 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1666 PetscReal ***aj, ***iaj, ***jaj, ***kaj, ***nvert;
1667 DMDAVecGetArray(fda, user->lCsi, &csi); DMDAVecGetArray(fda, user->lEta, &eta); DMDAVecGetArray(fda, user->lZet, &zet);
1668 DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
1669 DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
1670 DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
1671 DMDAVecGetArray(da, user->lAj, &aj); DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
1672 DMDAVecGetArray(da, user->lNvert, &nvert);
1673
1674 // --- Create temporary vectors for the metric tensor components G_ij ---
1675 Vec G11, G12, G13, G21, G22, G23, G31, G32, G33;
1676 PetscReal ***g11, ***g12, ***g13, ***g21, ***g22, ***g23, ***g31, ***g32, ***g33;
1677 VecDuplicate(user->lAj, &G11); VecDuplicate(user->lAj, &G12); VecDuplicate(user->lAj, &G13);
1678 VecDuplicate(user->lAj, &G21); VecDuplicate(user->lAj, &G22); VecDuplicate(user->lAj, &G23);
1679 VecDuplicate(user->lAj, &G31); VecDuplicate(user->lAj, &G32); VecDuplicate(user->lAj, &G33);
1680 DMDAVecGetArray(da, G11, &g11); DMDAVecGetArray(da, G12, &g12); DMDAVecGetArray(da, G13, &g13);
1681 DMDAVecGetArray(da, G21, &g21); DMDAVecGetArray(da, G22, &g22); DMDAVecGetArray(da, G23, &g23);
1682 DMDAVecGetArray(da, G31, &g31); DMDAVecGetArray(da, G32, &g32); DMDAVecGetArray(da, G33, &g33);
1683
1684 //================================================================================
1685 // Section 2: Pre-compute Metric Tensor Coefficients (g_ij)
1686 //================================================================================
1687 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pre-computing metric tensor components (g_ij).\n");
1688 for (k = gzs; k < gze; k++) {
1689 for (j = gys; j < gye; j++) {
1690 for (i = gxs; i < gxe; i++) {
1691 // These coefficients represent the dot products of the grid's contravariant base vectors,
1692 // scaled by face area. They are the core of the Laplacian operator on a curvilinear grid.
1693 if(i>-1 && j>-1 && k>-1 && i<IM+1 && j<JM+1 && k<KM+1){
1694 g11[k][j][i] = (icsi[k][j][i].x * icsi[k][j][i].x + icsi[k][j][i].y * icsi[k][j][i].y + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
1695 g12[k][j][i] = (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
1696 g13[k][j][i] = (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
1697 g21[k][j][i] = (jcsi[k][j][i].x * jeta[k][j][i].x + jcsi[k][j][i].y * jeta[k][j][i].y + jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
1698 g22[k][j][i] = (jeta[k][j][i].x * jeta[k][j][i].x + jeta[k][j][i].y * jeta[k][j][i].y + jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
1699 g23[k][j][i] = (jzet[k][j][i].x * jeta[k][j][i].x + jzet[k][j][i].y * jeta[k][j][i].y + jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
1700 g31[k][j][i] = (kcsi[k][j][i].x * kzet[k][j][i].x + kcsi[k][j][i].y * kzet[k][j][i].y + kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
1701 g32[k][j][i] = (keta[k][j][i].x * kzet[k][j][i].x + keta[k][j][i].y * kzet[k][j][i].y + keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
1702 g33[k][j][i] = (kzet[k][j][i].x * kzet[k][j][i].x + kzet[k][j][i].y * kzet[k][j][i].y + kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
1703 }
1704 }
1705 }
1706 }
1707
1708 //================================================================================
1709 // Section 3: Assemble the LHS Matrix A
1710 //================================================================================
1711 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling the LHS matrix A using a 19-point stencil.\n");
1712
1713 // --- Define domain boundaries for stencil logic, accounting for periodic BCs ---
1714 PetscInt x_str, x_end, y_str, y_end, z_str, z_end;
1715 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) { x_end = mx - 1; x_str = 0; }
1716 else { x_end = mx - 2; x_str = 1; }
1717 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) { y_end = my - 1; y_str = 0; }
1718 else { y_end = my - 2; y_str = 1; }
1719 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) { z_end = mz - 1; z_str = 0; }
1720 else { z_end = mz - 2; z_str = 1; }
1721
1722 // --- Main assembly loop over all local grid points ---
1723 for (k = zs; k < ze; k++) {
1724 for (j = ys; j < ye; j++) {
1725 for (i = xs; i < xe; i++) {
1726 PetscScalar vol[19]; // Holds the 19 stencil coefficient values for the current row
1727 PetscInt idx[19]; // Holds the 19 global column indices for the current row
1728 PetscInt row = Gidx(i, j, k, user); // Global index for the current row
1729
1730 // --- Handle Domain Boundary and Immersed Solid Points ---
1731 // For these points, we don't solve the Poisson equation. We set an identity
1732 // row (A_ii = 1) to effectively fix the pressure value (usually to 0).
1733 if (i == 0 || i == mx - 1 || j == 0 || j == my - 1 || k == 0 || k == mz - 1 || nvert[k][j][i] > IBM_FLUID_THRESHOLD) {
1734 vol[CP] = 1.0;
1735 idx[CP] = row;
1736 MatSetValues(user->A, 1, &row, 1, &idx[CP], &vol[CP], INSERT_VALUES);
1737 }
1738 // --- Handle Fluid Points ---
1739 else {
1740 for (PetscInt m = 0; m < 19; m++) {
1741 vol[m] = 0.0;
1742 }
1743
1744 /************************************************************************
1745 * EAST FACE CONTRIBUTION (between i and i+1)
1746 ************************************************************************/
1747 if (nvert[k][j][i + 1] < IBM_FLUID_THRESHOLD && i != x_end) { // East neighbor is fluid
1748 // Primary derivative term: d/d_csi (g11 * dP/d_csi)
1749 vol[CP] -= g11[k][j][i];
1750 vol[EP] += g11[k][j][i];
1751
1752 // Cross-derivative term: d/d_csi (g12 * dP/d_eta).
1753 // This requires an average of dP/d_eta. If a neighbor is solid, the stencil
1754 // dynamically switches to a one-sided difference to avoid using solid points.
1755 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
1756 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1757 vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
1758 vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
1759 }
1760 }
1761 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
1762 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1763 vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
1764 vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
1765 }
1766 }
1767 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1768 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1769 vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
1770 vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
1771 }
1772 }
1773 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1774 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1775 vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
1776 vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
1777 }
1778 }
1779 else { // Centered difference
1780 vol[NP] += g12[k][j][i] * 0.25; vol[NE] += g12[k][j][i] * 0.25;
1781 vol[SP] -= g12[k][j][i] * 0.25; vol[SE] -= g12[k][j][i] * 0.25;
1782 }
1783
1784 // Cross-derivative term: d/d_csi (g13 * dP/d_zet)
1785 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1786 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1787 vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
1788 vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
1789 }
1790 }
1791 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1792 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1793 vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
1794 vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
1795 }
1796 }
1797 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1798 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1799 vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
1800 vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
1801 }
1802 }
1803 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1804 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1805 vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
1806 vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
1807 }
1808 }
1809 else { // Centered difference
1810 vol[TP] += g13[k][j][i] * 0.25; vol[TE] += g13[k][j][i] * 0.25;
1811 vol[BP] -= g13[k][j][i] * 0.25; vol[BE] -= g13[k][j][i] * 0.25;
1812 }
1813 }
1814
1815 /************************************************************************
1816 * WEST FACE CONTRIBUTION (between i-1 and i)
1817 ************************************************************************/
1818 if (nvert[k][j][i-1] < IBM_FLUID_THRESHOLD && i != x_str) {
1819 vol[CP] -= g11[k][j][i-1];
1820 vol[WP] += g11[k][j][i-1];
1821
1822 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
1823 if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1824 vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
1825 vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
1826 }
1827 }
1828 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
1829 if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1) {
1830 vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
1831 vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
1832 }
1833 }
1834 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
1835 if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
1836 vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
1837 vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
1838 }
1839 }
1840 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
1841 if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
1842 vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
1843 vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
1844 }
1845 }
1846 else {
1847 vol[NP] -= g12[k][j][i-1] * 0.25; vol[NW] -= g12[k][j][i-1] * 0.25;
1848 vol[SP] += g12[k][j][i-1] * 0.25; vol[SW] += g12[k][j][i-1] * 0.25;
1849 }
1850
1851 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
1852 if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1853 vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
1854 vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
1855 }
1856 }
1857 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
1858 if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1) {
1859 vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
1860 vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
1861 }
1862 }
1863 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
1864 if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
1865 vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
1866 vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
1867 }
1868 }
1869 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
1870 if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
1871 vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
1872 vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
1873 }
1874 }
1875 else {
1876 vol[TP] -= g13[k][j][i-1] * 0.25; vol[TW] -= g13[k][j][i-1] * 0.25;
1877 vol[BP] += g13[k][j][i-1] * 0.25; vol[BW] += g13[k][j][i-1] * 0.25;
1878 }
1879 }
1880
1881 /************************************************************************
1882 * NORTH FACE CONTRIBUTION (between j and j+1)
1883 ************************************************************************/
1884 if (nvert[k][j+1][i] < IBM_FLUID_THRESHOLD && j != y_end) {
1885 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1886 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1887 vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
1888 vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
1889 }
1890 }
1891 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1892 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1893 vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
1894 vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
1895 }
1896 }
1897 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1898 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1899 vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
1900 vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
1901 }
1902 }
1903 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1904 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1905 vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
1906 vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
1907 }
1908 }
1909 else {
1910 vol[EP] += g21[k][j][i] * 0.25; vol[NE] += g21[k][j][i] * 0.25;
1911 vol[WP] -= g21[k][j][i] * 0.25; vol[NW] -= g21[k][j][i] * 0.25;
1912 }
1913
1914 vol[CP] -= g22[k][j][i];
1915 vol[NP] += g22[k][j][i];
1916
1917 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1918 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1919 vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
1920 vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
1921 }
1922 }
1923 else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1924 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1925 vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
1926 vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
1927 }
1928 }
1929 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1930 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1931 vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
1932 vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
1933 }
1934 }
1935 else if ((k == 1 || k==mz-2 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1936 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1937 vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
1938 vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
1939 }
1940 }
1941 else {
1942 vol[TP] += g23[k][j][i] * 0.25; vol[TN] += g23[k][j][i] * 0.25;
1943 vol[BP] -= g23[k][j][i] * 0.25; vol[BN] -= g23[k][j][i] * 0.25;
1944 }
1945 }
1946
1947 /************************************************************************
1948 * SOUTH FACE CONTRIBUTION (between j-1 and j)
1949 ************************************************************************/
1950 if (nvert[k][j-1][i] < IBM_FLUID_THRESHOLD && j != y_str) {
1951 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
1952 if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1953 vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
1954 vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
1955 }
1956 }
1957 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
1958 if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1) {
1959 vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
1960 vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
1961 }
1962 }
1963 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
1964 if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
1965 vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
1966 vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
1967 }
1968 }
1969 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
1970 if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
1971 vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
1972 vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
1973 }
1974 }
1975 else {
1976 vol[EP] -= g21[k][j-1][i] * 0.25; vol[SE] -= g21[k][j-1][i] * 0.25;
1977 vol[WP] += g21[k][j-1][i] * 0.25; vol[SW] += g21[k][j-1][i] * 0.25;
1978 }
1979
1980 vol[CP] -= g22[k][j-1][i];
1981 vol[SP] += g22[k][j-1][i];
1982
1983 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
1984 if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1985 vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
1986 vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
1987 }
1988 }
1989 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
1990 if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 ) {
1991 vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
1992 vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
1993 }
1994 }
1995 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
1996 if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
1997 vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
1998 vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
1999 }
2000 }
2001 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
2002 if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
2003 vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
2004 vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
2005 }
2006 }
2007 else {
2008 vol[TP] -= g23[k][j-1][i] * 0.25; vol[TS] -= g23[k][j-1][i] * 0.25;
2009 vol[BP] += g23[k][j-1][i] * 0.25; vol[BS] += g23[k][j-1][i] * 0.25;
2010 }
2011 }
2012
2013 /************************************************************************
2014 * TOP FACE CONTRIBUTION (between k and k+1)
2015 ************************************************************************/
2016 if (nvert[k+1][j][i] < IBM_FLUID_THRESHOLD && k != z_end) {
2017 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
2018 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
2019 vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
2020 vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
2021 }
2022 }
2023 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
2024 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
2025 vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
2026 vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
2027 }
2028 }
2029 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
2030 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
2031 vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
2032 vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
2033 }
2034 }
2035 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
2036 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
2037 vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
2038 vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
2039 }
2040 }
2041 else {
2042 vol[EP] += g31[k][j][i] * 0.25; vol[TE] += g31[k][j][i] * 0.25;
2043 vol[WP] -= g31[k][j][i] * 0.25; vol[TW] -= g31[k][j][i] * 0.25;
2044 }
2045
2046 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
2047 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
2048 vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
2049 vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
2050 }
2051 }
2052 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
2053 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
2054 vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
2055 vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
2056 }
2057 }
2058 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
2059 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
2060 vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
2061 vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
2062 }
2063 }
2064 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
2065 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
2066 vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
2067 vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
2068 }
2069 }
2070 else {
2071 vol[NP] += g32[k][j][i] * 0.25; vol[TN] += g32[k][j][i] * 0.25;
2072 vol[SP] -= g32[k][j][i] * 0.25; vol[TS] -= g32[k][j][i] * 0.25;
2073 }
2074
2075 vol[CP] -= g33[k][j][i];
2076 vol[TP] += g33[k][j][i];
2077 }
2078
2079 /************************************************************************
2080 * BOTTOM FACE CONTRIBUTION (between k-1 and k)
2081 ************************************************************************/
2082 if (nvert[k-1][j][i] < IBM_FLUID_THRESHOLD && k != z_str) {
2083 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
2084 if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
2085 vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
2086 vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
2087 }
2088 }
2089 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
2090 if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1) {
2091 vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
2092 vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
2093 }
2094 }
2095 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
2096 if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
2097 vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
2098 vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
2099 }
2100 }
2101 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
2102 if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
2103 vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
2104 vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
2105 }
2106 }
2107 else {
2108 vol[EP] -= g31[k-1][j][i] * 0.25; vol[BE] -= g31[k-1][j][i] * 0.25;
2109 vol[WP] += g31[k-1][j][i] * 0.25; vol[BW] += g31[k-1][j][i] * 0.25;
2110 }
2111
2112 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
2113 if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
2114 vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
2115 vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
2116 }
2117 }
2118 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
2119 if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1) {
2120 vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
2121 vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
2122 }
2123 }
2124 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
2125 if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
2126 vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
2127 vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
2128 }
2129 }
2130 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
2131 if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
2132 vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
2133 vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
2134 }
2135 }
2136 else {
2137 vol[NP] -= g32[k-1][j][i] * 0.25; vol[BN] -= g32[k-1][j][i] * 0.25;
2138 vol[SP] += g32[k-1][j][i] * 0.25; vol[BS] += g32[k-1][j][i] * 0.25;
2139 }
2140
2141 vol[CP] -= g33[k-1][j][i];
2142 vol[BP] += g33[k-1][j][i];
2143 }
2144
2145 // --- Final scaling and insertion into the matrix ---
2146
2147 // Scale all stencil coefficients by the negative cell volume (-aj).
2148 for (PetscInt m = 0; m < 19; m++) {
2149 vol[m] *= -aj[k][j][i];
2150 }
2151
2152 // Set the global column indices for the 19 stencil points, handling periodic BCs.
2153 idx[CP] = Gidx(i, j, k, user);
2154 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[EP] = Gidx(1, j, k, user); else idx[EP] = Gidx(i+1, j, k, user);
2155 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[WP] = Gidx(mx-2, j, k, user); else idx[WP] = Gidx(i-1, j, k, user);
2156 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[NP] = Gidx(i, 1, k, user); else idx[NP] = Gidx(i, j+1, k, user);
2157 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[SP] = Gidx(i, my-2, k, user); else idx[SP] = Gidx(i, j-1, k, user);
2158 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TP] = Gidx(i, j, 1, user); else idx[TP] = Gidx(i, j, k+1, user);
2159 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BP] = Gidx(i, j, mz-2, user); else idx[BP] = Gidx(i, j, k-1, user);
2160 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==mx-2 && j==my-2) idx[NE] = Gidx(1, 1, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[NE] = Gidx(1, j+1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[NE] = Gidx(i+1, 1, k, user); else idx[NE] = Gidx(i+1, j+1, k, user);
2161 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==mx-2 && j==1) idx[SE] = Gidx(1, my-2, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[SE] = Gidx(1, j-1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[SE] = Gidx(i+1, my-2, k, user); else idx[SE] = Gidx(i+1, j-1, k, user);
2162 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==1 && j==my-2) idx[NW] = Gidx(mx-2, 1, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[NW] = Gidx(mx-2, j+1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[NW] = Gidx(i-1, 1, k, user); else idx[NW] = Gidx(i-1, j+1, k, user);
2163 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==1 && j==1) idx[SW] = Gidx(mx-2, my-2, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[SW] = Gidx(mx-2, j-1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[SW] = Gidx(i-1, my-2, k, user); else idx[SW] = Gidx(i-1, j-1, k, user);
2164 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==my-2 && k==mz-2) idx[TN] = Gidx(i, 1, 1, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[TN] = Gidx(i, 1, k+1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TN] = Gidx(i, j+1, 1, user); else idx[TN] = Gidx(i, j+1, k+1, user);
2165 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==my-2 && k==1) idx[BN] = Gidx(i, 1, mz-2, user); else if(user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[BN] = Gidx(i, 1, k-1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BN] = Gidx(i, j+1, mz-2, user); else idx[BN] = Gidx(i, j+1, k-1, user);
2166 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==1 && k==mz-2) idx[TS] = Gidx(i, my-2, 1, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[TS] = Gidx(i, my-2, k+1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TS] = Gidx(i, j-1, 1, user); else idx[TS] = Gidx(i, j-1, k+1, user);
2167 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==1 && k==1) idx[BS] = Gidx(i, my-2, mz-2, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[BS] = Gidx(i, my-2, k-1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BS] = Gidx(i, j-1, mz-2, user); else idx[BS] = Gidx(i, j-1, k-1, user);
2168 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==mx-2 && k==mz-2) idx[TE] = Gidx(1, j, 1, user); else if(user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[TE] = Gidx(1, j, k+1, user); else if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TE] = Gidx(i+1, j, 1, user); else idx[TE] = Gidx(i+1, j, k+1, user);
2169 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==mx-2 && k==1) idx[BE] = Gidx(1, j, mz-2, user); else if(user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[BE] = Gidx(1, j, k-1, user); else if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BE] = Gidx(i+1, j, mz-2, user); else idx[BE] = Gidx(i+1, j, k-1, user);
2170 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==1 && k==mz-2) idx[TW] = Gidx(mx-2, j, 1, user); else if(user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[TW] = Gidx(mx-2, j, k+1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TW] = Gidx(i-1, j, 1, user); else idx[TW] = Gidx(i-1, j, k+1, user);
2171 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==1 && k==1) idx[BW] = Gidx(mx-2, j, mz-2, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[BW] = Gidx(mx-2, j, k-1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BW] = Gidx(i-1, j, mz-2, user); else idx[BW] = Gidx(i-1, j, k-1, user);
2172
2173 // Insert the computed row into the matrix A.
2174 MatSetValues(user->A, 1, &row, 19, idx, vol, INSERT_VALUES);
2175 }
2176 }
2177 }
2178 }
2179
2180 //================================================================================
2181 // Section 4: Finalize Matrix and Cleanup
2182 //================================================================================
2183
2184 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing matrix assembly.\n");
2185 MatAssemblyBegin(user->A, MAT_FINAL_ASSEMBLY);
2186 MatAssemblyEnd(user->A, MAT_FINAL_ASSEMBLY);
2187
2188 PetscReal max_A;
2189
2190 ierr = MatNorm(user->A,NORM_INFINITY,&max_A);CHKERRQ(ierr);
2191
2192 LOG_ALLOW(GLOBAL,LOG_DEBUG," Max value in A matrix for level %d = %le.\n",user->thislevel,max_A);
2193
2194 // if (get_log_level() >= LOG_DEBUG) {
2195 // ierr = MatView(user->A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
2196 // }
2197
2198 // --- Restore access to all PETSc vectors and destroy temporary ones ---
2199 DMDAVecRestoreArray(da, G11, &g11); DMDAVecRestoreArray(da, G12, &g12); DMDAVecRestoreArray(da, G13, &g13);
2200 DMDAVecRestoreArray(da, G21, &g21); DMDAVecRestoreArray(da, G22, &g22); DMDAVecRestoreArray(da, G23, &g23);
2201 DMDAVecRestoreArray(da, G31, &g31); DMDAVecRestoreArray(da, G32, &g32); DMDAVecRestoreArray(da, G33, &g33);
2202
2203 VecDestroy(&G11); VecDestroy(&G12); VecDestroy(&G13);
2204 VecDestroy(&G21); VecDestroy(&G22); VecDestroy(&G23);
2205 VecDestroy(&G31); VecDestroy(&G32); VecDestroy(&G33);
2206
2207 DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
2208 DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
2209 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
2210 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
2211 DMDAVecRestoreArray(da, user->lAj, &aj); DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
2212 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2213
2214 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting PoissonLHSNew.\n");
2215 PROFILE_FUNCTION_END;
2216 PetscFunctionReturn(0);
2217}
2218
2219
2220#undef __FUNCT__
2221#define __FUNCT__ "PoissonRHS"
2222
2223PetscErrorCode PoissonRHS(UserCtx *user, Vec B)
2224{
2225 DMDALocalInfo info = user->info;
2226 PetscInt xs = info.xs, xe = info.xs + info.xm;
2227 PetscInt ys = info.ys, ye = info.ys + info.ym;
2228 PetscInt zs = info.zs, ze = info.zs + info.zm;
2229 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2230
2231 PetscInt i, j, k;
2232 PetscReal ***nvert, ***aj, ***rb, dt = user->simCtx->dt;
2233 struct Components{
2234 PetscReal x;
2235 PetscReal y;
2236 PetscReal z;
2237 } *** ucont;
2238
2239 PROFILE_FUNCTION_BEGIN;
2240
2241 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonRHS to compute pressure equation RHS.\n");
2242
2243 DMDAVecGetArray(user->da, B, &rb);
2244 DMDAVecGetArray(user->fda, user->lUcont, &ucont);
2245 DMDAVecGetArray(user->da, user->lNvert, &nvert);
2246 DMDAVecGetArray(user->da, user->lAj, &aj);
2247
2248
2249 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing RHS values for each cell.\n");
2250
2251 for (k=zs; k<ze; k++) {
2252 for (j=ys; j<ye; j++) {
2253 for (i=xs; i<xe; i++) {
2254
2255 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1) {
2256 rb[k][j][i] = 0.;
2257 }
2258 else if (nvert[k][j][i] > 0.1) {
2259 rb[k][j][i] = 0;
2260 }
2261 else {
2262 rb[k][j][i] = -(ucont[k][j][i].x - ucont[k][j][i-1].x +
2263 ucont[k][j][i].y - ucont[k][j-1][i].y +
2264 ucont[k][j][i].z - ucont[k-1][j][i].z) / dt
2265 * aj[k][j][i] / 1.0 * COEF_TIME_ACCURACY; // user->simCtx->st replaced by 1.0.
2266
2267 }
2268 }
2269 }
2270 }
2271
2272
2273 // --- Check the solvability condition for the Poisson equation ---
2274 // The global sum of the RHS (proportional to the total divergence) must be zero.
2275 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying solvability condition (sum of RHS terms).\n");
2276 PetscReal lsum=0., sum=0.;
2277
2278 for (k=zs; k<ze; k++) {
2279 for (j=ys; j<ye; j++) {
2280 for (i=xs; i<xe; i++) {
2281
2282 lsum += rb[k][j][i] / aj[k][j][i]* dt/COEF_TIME_ACCURACY;
2283
2284 }
2285 }
2286 }
2287
2288 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2289
2290 LOG_ALLOW(GLOBAL, LOG_INFO, "Global Sum of RHS (Divergence Check): %le\n", sum);
2291
2292 user->simCtx->summationRHS = sum;
2293
2294 DMDAVecRestoreArray(user->fda, user->lUcont, &ucont);
2295 DMDAVecRestoreArray(user->da, user->lNvert, &nvert);
2296 DMDAVecRestoreArray(user->da, user->lAj, &aj);
2297 DMDAVecRestoreArray(user->da, B, &rb);
2298
2299
2300 PROFILE_FUNCTION_END;
2301 return 0;
2302}
2303
2304PetscErrorCode VolumeFlux_rev(UserCtx *user, PetscReal *ibm_Flux,
2305 PetscReal *ibm_Area, PetscInt flg)
2306{
2307
2308 // --- CONTEXT ACQUISITION BLOCK ---
2309 // Get the master simulation context from the UserCtx.
2310 SimCtx *simCtx = user->simCtx;
2311
2312 // Create a local variable to mirror the legacy global for minimal code changes.
2313 const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
2314 // --- END CONTEXT ACQUISITION BLOCK ---
2315
2316 DM da = user->da, fda = user->fda;
2317
2318 DMDALocalInfo info = user->info;
2319
2320 PetscInt xs = info.xs, xe = info.xs + info.xm;
2321 PetscInt ys = info.ys, ye = info.ys + info.ym;
2322 PetscInt zs = info.zs, ze = info.zs + info.zm;
2323 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2324
2325 PetscInt i, j, k;
2326 PetscInt lxs, lys, lzs, lxe, lye, lze;
2327
2328 lxs = xs; lxe = xe;
2329 lys = ys; lye = ye;
2330 lzs = zs; lze = ze;
2331
2332 if (xs==0) lxs = xs+1;
2333 if (ys==0) lys = ys+1;
2334 if (zs==0) lzs = zs+1;
2335
2336 if (xe==mx) lxe = xe-1;
2337 if (ye==my) lye = ye-1;
2338 if (ze==mz) lze = ze-1;
2339
2340 PetscReal ***nvert, ibmval=1.5;
2341 Cmpnts ***ucor, ***csi, ***eta, ***zet;
2342 DMDAVecGetArray(fda, user->Ucont, &ucor);
2343 DMDAVecGetArray(fda, user->lCsi, &csi);
2344 DMDAVecGetArray(fda, user->lEta, &eta);
2345 DMDAVecGetArray(fda, user->lZet, &zet);
2346 DMDAVecGetArray(da, user->lNvert, &nvert);
2347
2348 PetscReal libm_Flux, libm_area;
2349 libm_Flux = 0;
2350 libm_area = 0;
2351 for (k=lzs; k<lze; k++) {
2352 for (j=lys; j<lye; j++) {
2353 for (i=lxs; i<lxe; i++) {
2354 if (nvert[k][j][i] < 0.1) {
2355 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
2356 libm_Flux += ucor[k][j][i].x;
2357 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2358 csi[k][j][i].y * csi[k][j][i].y +
2359 csi[k][j][i].z * csi[k][j][i].z);
2360
2361 }
2362 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
2363 libm_Flux += ucor[k][j][i].y;
2364 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2365 eta[k][j][i].y * eta[k][j][i].y +
2366 eta[k][j][i].z * eta[k][j][i].z);
2367 }
2368 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
2369 libm_Flux += ucor[k][j][i].z;
2370 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2371 zet[k][j][i].y * zet[k][j][i].y +
2372 zet[k][j][i].z * zet[k][j][i].z);
2373 }
2374 }
2375
2376 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
2377 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
2378 libm_Flux -= ucor[k][j][i].x;
2379 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2380 csi[k][j][i].y * csi[k][j][i].y +
2381 csi[k][j][i].z * csi[k][j][i].z);
2382
2383 }
2384 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
2385 libm_Flux -= ucor[k][j][i].y;
2386 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2387 eta[k][j][i].y * eta[k][j][i].y +
2388 eta[k][j][i].z * eta[k][j][i].z);
2389 }
2390 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
2391 libm_Flux -= ucor[k][j][i].z;
2392 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2393 zet[k][j][i].y * zet[k][j][i].y +
2394 zet[k][j][i].z * zet[k][j][i].z);
2395 }
2396 }
2397
2398 }
2399 }
2400 }
2401
2402 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2403 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2404
2405 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
2406/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
2407 PetscPrintf(PETSC_COMM_WORLD, "IBMFlux REV %le %le\n", *ibm_Flux, *ibm_Area);
2408
2409 PetscReal correction;
2410
2411 if (*ibm_Area > 1.e-15) {
2412 if (flg)
2413 correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
2414 else
2415 correction = *ibm_Flux / *ibm_Area;
2416 }
2417 else {
2418 correction = 0;
2419 }
2420
2421 for (k=lzs; k<lze; k++) {
2422 for (j=lys; j<lye; j++) {
2423 for (i=lxs; i<lxe; i++) {
2424 if (nvert[k][j][i] < 0.1) {
2425 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
2426 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
2427 csi[k][j][i].y * csi[k][j][i].y +
2428 csi[k][j][i].z * csi[k][j][i].z) *
2429 correction;
2430
2431 }
2432 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
2433 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
2434 eta[k][j][i].y * eta[k][j][i].y +
2435 eta[k][j][i].z * eta[k][j][i].z) *
2436 correction;
2437 }
2438 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
2439 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
2440 zet[k][j][i].y * zet[k][j][i].y +
2441 zet[k][j][i].z * zet[k][j][i].z) *
2442 correction;
2443 }
2444 }
2445
2446 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
2447 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
2448 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2449 csi[k][j][i].y * csi[k][j][i].y +
2450 csi[k][j][i].z * csi[k][j][i].z) *
2451 correction;
2452
2453 }
2454 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
2455 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2456 eta[k][j][i].y * eta[k][j][i].y +
2457 eta[k][j][i].z * eta[k][j][i].z) *
2458 correction;
2459 }
2460 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
2461 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2462 zet[k][j][i].y * zet[k][j][i].y +
2463 zet[k][j][i].z * zet[k][j][i].z) *
2464 correction;
2465 }
2466 }
2467
2468 }
2469 }
2470 }
2471
2472
2473
2474 libm_Flux = 0;
2475 libm_area = 0;
2476 for (k=lzs; k<lze; k++) {
2477 for (j=lys; j<lye; j++) {
2478 for (i=lxs; i<lxe; i++) {
2479 if (nvert[k][j][i] < 0.1) {
2480 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
2481 libm_Flux += ucor[k][j][i].x;
2482 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2483 csi[k][j][i].y * csi[k][j][i].y +
2484 csi[k][j][i].z * csi[k][j][i].z);
2485
2486 }
2487 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
2488 libm_Flux += ucor[k][j][i].y;
2489 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2490 eta[k][j][i].y * eta[k][j][i].y +
2491 eta[k][j][i].z * eta[k][j][i].z);
2492 }
2493 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
2494 libm_Flux += ucor[k][j][i].z;
2495 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2496 zet[k][j][i].y * zet[k][j][i].y +
2497 zet[k][j][i].z * zet[k][j][i].z);
2498 }
2499 }
2500
2501 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
2502 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
2503 libm_Flux -= ucor[k][j][i].x;
2504 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2505 csi[k][j][i].y * csi[k][j][i].y +
2506 csi[k][j][i].z * csi[k][j][i].z);
2507
2508 }
2509 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
2510 libm_Flux -= ucor[k][j][i].y;
2511 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2512 eta[k][j][i].y * eta[k][j][i].y +
2513 eta[k][j][i].z * eta[k][j][i].z);
2514 }
2515 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
2516 libm_Flux -= ucor[k][j][i].z;
2517 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2518 zet[k][j][i].y * zet[k][j][i].y +
2519 zet[k][j][i].z * zet[k][j][i].z);
2520 }
2521 }
2522
2523 }
2524 }
2525 }
2526
2527 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2528 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2529
2530 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
2531/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
2532 PetscPrintf(PETSC_COMM_WORLD, "IBMFlux22 REV %le %le\n", *ibm_Flux, *ibm_Area);
2533
2534 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2535 DMDAVecRestoreArray(fda, user->lCsi, &csi);
2536 DMDAVecRestoreArray(fda, user->lEta, &eta);
2537 DMDAVecRestoreArray(fda, user->lZet, &zet);
2538 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
2539
2540 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2541 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
2542 return 0;
2543}
2544
2545
2546PetscErrorCode VolumeFlux(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
2547{
2548 // --- CONTEXT ACQUISITION BLOCK ---
2549 // Get the master simulation context from the UserCtx.
2550 SimCtx *simCtx = user->simCtx;
2551
2552 // Create local variables to mirror the legacy globals for minimal code changes.
2553 const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
2554 const PetscInt channelz = simCtx->channelz;
2555 const PetscInt fish = simCtx->fish;
2556 // --- END CONTEXT ACQUISITION BLOCK ---
2557
2558 DM da = user->da, fda = user->fda;
2559
2560 DMDALocalInfo info = user->info;
2561
2562 PetscInt xs = info.xs, xe = info.xs + info.xm;
2563 PetscInt ys = info.ys, ye = info.ys + info.ym;
2564 PetscInt zs = info.zs, ze = info.zs + info.zm;
2565 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2566
2567 PetscInt i, j, k,ibi;
2568 PetscInt lxs, lys, lzs, lxe, lye, lze;
2569
2570 PetscInt gxs, gxe, gys, gye, gzs, gze;
2571
2572 gxs = info.gxs; gxe = gxs + info.gxm;
2573 gys = info.gys; gye = gys + info.gym;
2574 gzs = info.gzs; gze = gzs + info.gzm;
2575
2576 lxs = xs; lxe = xe;
2577 lys = ys; lye = ye;
2578 lzs = zs; lze = ze;
2579
2580 if (xs==0) lxs = xs+1;
2581 if (ys==0) lys = ys+1;
2582 if (zs==0) lzs = zs+1;
2583
2584 if (xe==mx) lxe = xe-1;
2585 if (ye==my) lye = ye-1;
2586 if (ze==mz) lze = ze-1;
2587
2588 PetscReal epsilon=1.e-8;
2589 PetscReal ***nvert, ibmval=1.9999;
2590
2591 struct Components {
2592 PetscReal x;
2593 PetscReal y;
2594 PetscReal z;
2595 }***ucor, ***lucor,***csi, ***eta, ***zet;
2596
2597
2598 PetscInt xend=mx-2 ,yend=my-2,zend=mz-2;
2599
2600 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) xend=mx-1;
2601 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) yend=my-1;
2602 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) zend=mz-1;
2603
2604 DMDAVecGetArray(fda, user->Ucont, &ucor);
2605 DMDAVecGetArray(fda, user->lCsi, &csi);
2606 DMDAVecGetArray(fda, user->lEta, &eta);
2607 DMDAVecGetArray(fda, user->lZet, &zet);
2608 DMDAVecGetArray(da, user->lNvert, &nvert);
2609
2610 PetscReal libm_Flux, libm_area, libm_Flux_abs=0., ibm_Flux_abs;
2611 libm_Flux = 0;
2612 libm_area = 0;
2613
2614 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering VolumeFlux to enforce no-penetration condition.\n");
2615
2616 //Mohsen March 2017
2617 PetscReal *lIB_Flux, *lIB_area,*IB_Flux,*IB_Area;
2618 if (NumberOfBodies > 1) {
2619
2620 lIB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
2621 lIB_area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
2622 IB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
2623 IB_Area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
2624
2625
2626 for (ibi=0; ibi<NumberOfBodies; ibi++) {
2627 lIB_Flux[ibi]=0.0;
2628 lIB_area[ibi]=0.0;
2629 IB_Flux[ibi]=0.0;
2630 IB_Area[ibi]=0.0;
2631 }
2632 }
2633
2634
2635 //================================================================================
2636 // PASS 1: Calculate Uncorrected Flux and Area
2637 // This pass measures the total fluid "leakage" across the immersed boundary.
2638 //================================================================================
2639 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 1: Measuring uncorrected flux and area.\n");
2640
2641 for (k=lzs; k<lze; k++) {
2642 for (j=lys; j<lye; j++) {
2643 for (i=lxs; i<lxe; i++) {
2644 if (nvert[k][j][i] < 0.1) {
2645 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
2646
2647 if (fabs(ucor[k][j][i].x)>epsilon) {
2648 libm_Flux += ucor[k][j][i].x;
2649 if (flg==3)
2650 libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
2651 csi[k][j][i].y * csi[k][j][i].y +
2652 csi[k][j][i].z * csi[k][j][i].z);
2653 else
2654 libm_Flux_abs += fabs(ucor[k][j][i].x);
2655
2656 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2657 csi[k][j][i].y * csi[k][j][i].y +
2658 csi[k][j][i].z * csi[k][j][i].z);
2659
2660 if (NumberOfBodies > 1) {
2661
2662 ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
2663 lIB_Flux[ibi] += ucor[k][j][i].x;
2664 lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2665 csi[k][j][i].y * csi[k][j][i].y +
2666 csi[k][j][i].z * csi[k][j][i].z);
2667 }
2668 } else
2669 ucor[k][j][i].x=0.;
2670
2671 }
2672 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
2673
2674 if (fabs(ucor[k][j][i].y)>epsilon) {
2675 libm_Flux += ucor[k][j][i].y;
2676 if (flg==3)
2677 libm_Flux_abs += fabs(ucor[k][j][i].y)/sqrt(eta[k][j][i].x * eta[k][j][i].x +
2678 eta[k][j][i].y * eta[k][j][i].y +
2679 eta[k][j][i].z * eta[k][j][i].z);
2680 else
2681 libm_Flux_abs += fabs(ucor[k][j][i].y);
2682 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2683 eta[k][j][i].y * eta[k][j][i].y +
2684 eta[k][j][i].z * eta[k][j][i].z);
2685 if (NumberOfBodies > 1) {
2686
2687 ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
2688
2689 lIB_Flux[ibi] += ucor[k][j][i].y;
2690 lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2691 eta[k][j][i].y * eta[k][j][i].y +
2692 eta[k][j][i].z * eta[k][j][i].z);
2693 }
2694 } else
2695 ucor[k][j][i].y=0.;
2696 }
2697 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
2698
2699 if (fabs(ucor[k][j][i].z)>epsilon) {
2700 libm_Flux += ucor[k][j][i].z;
2701 if (flg==3)
2702 libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
2703 zet[k][j][i].y * zet[k][j][i].y +
2704 zet[k][j][i].z * zet[k][j][i].z);
2705 else
2706 libm_Flux_abs += fabs(ucor[k][j][i].z);
2707 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2708 zet[k][j][i].y * zet[k][j][i].y +
2709 zet[k][j][i].z * zet[k][j][i].z);
2710
2711 if (NumberOfBodies > 1) {
2712
2713 ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
2714 lIB_Flux[ibi] += ucor[k][j][i].z;
2715 lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2716 zet[k][j][i].y * zet[k][j][i].y +
2717 zet[k][j][i].z * zet[k][j][i].z);
2718 }
2719 }else
2720 ucor[k][j][i].z=0.;
2721 }
2722 }
2723
2724 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
2725
2726 if (nvert[k][j][i+1] < 0.1 && i < xend) {
2727 if (fabs(ucor[k][j][i].x)>epsilon) {
2728 libm_Flux -= ucor[k][j][i].x;
2729 if (flg==3)
2730 libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
2731 csi[k][j][i].y * csi[k][j][i].y +
2732 csi[k][j][i].z * csi[k][j][i].z);
2733 else
2734 libm_Flux_abs += fabs(ucor[k][j][i].x);
2735 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2736 csi[k][j][i].y * csi[k][j][i].y +
2737 csi[k][j][i].z * csi[k][j][i].z);
2738 if (NumberOfBodies > 1) {
2739
2740 ibi=(int)((nvert[k][j][i]-1.0)*1001);
2741 lIB_Flux[ibi] -= ucor[k][j][i].x;
2742 lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2743 csi[k][j][i].y * csi[k][j][i].y +
2744 csi[k][j][i].z * csi[k][j][i].z);
2745 }
2746
2747 }else
2748 ucor[k][j][i].x=0.;
2749 }
2750 if (nvert[k][j+1][i] < 0.1 && j < yend) {
2751 if (fabs(ucor[k][j][i].y)>epsilon) {
2752 libm_Flux -= ucor[k][j][i].y;
2753 if (flg==3)
2754 libm_Flux_abs += fabs(ucor[k][j][i].y)/ sqrt(eta[k][j][i].x * eta[k][j][i].x +
2755 eta[k][j][i].y * eta[k][j][i].y +
2756 eta[k][j][i].z * eta[k][j][i].z);
2757 else
2758 libm_Flux_abs += fabs(ucor[k][j][i].y);
2759 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2760 eta[k][j][i].y * eta[k][j][i].y +
2761 eta[k][j][i].z * eta[k][j][i].z);
2762 if (NumberOfBodies > 1) {
2763
2764 ibi=(int)((nvert[k][j][i]-1.0)*1001);
2765 lIB_Flux[ibi] -= ucor[k][j][i].y;
2766 lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2767 eta[k][j][i].y * eta[k][j][i].y +
2768 eta[k][j][i].z * eta[k][j][i].z);
2769 }
2770 }else
2771 ucor[k][j][i].y=0.;
2772 }
2773 if (nvert[k+1][j][i] < 0.1 && k < zend) {
2774 if (fabs(ucor[k][j][i].z)>epsilon) {
2775 libm_Flux -= ucor[k][j][i].z;
2776 if (flg==3)
2777 libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
2778 zet[k][j][i].y * zet[k][j][i].y +
2779 zet[k][j][i].z * zet[k][j][i].z);
2780 else
2781 libm_Flux_abs += fabs(ucor[k][j][i].z);
2782 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2783 zet[k][j][i].y * zet[k][j][i].y +
2784 zet[k][j][i].z * zet[k][j][i].z);
2785 if (NumberOfBodies > 1) {
2786
2787 ibi=(int)((nvert[k][j][i]-1.0)*1001);
2788 lIB_Flux[ibi] -= ucor[k][j][i].z;
2789 lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2790 zet[k][j][i].y * zet[k][j][i].y +
2791 zet[k][j][i].z * zet[k][j][i].z);
2792 }
2793 }else
2794 ucor[k][j][i].z=0.;
2795 }
2796 }
2797
2798 }
2799 }
2800 }
2801
2802 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2803 MPI_Allreduce(&libm_Flux_abs, &ibm_Flux_abs,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2804 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2805
2806 if (NumberOfBodies > 1) {
2807 MPI_Allreduce(lIB_Flux,IB_Flux,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
2808 MPI_Allreduce(lIB_area,IB_Area,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
2809 }
2810
2811 PetscReal correction;
2812
2813 PetscReal *Correction;
2814 if (NumberOfBodies > 1) {
2815 Correction=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
2816 for (ibi=0; ibi<NumberOfBodies; ibi++) Correction[ibi]=0.0;
2817 }
2818
2819 if (*ibm_Area > 1.e-15) {
2820 if (flg>1)
2821 correction = (*ibm_Flux + user->FluxIntpSum)/ ibm_Flux_abs;
2822 else if (flg)
2823 correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
2824 else
2825 correction = *ibm_Flux / *ibm_Area;
2826 if (NumberOfBodies > 1)
2827 for (ibi=0; ibi<NumberOfBodies; ibi++) if (IB_Area[ibi]>1.e-15) Correction[ibi] = IB_Flux[ibi] / IB_Area[ibi];
2828 }
2829 else {
2830 correction = 0;
2831 }
2832 // --- Log the uncorrected results and calculated correction ---
2833 LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Uncorrected Flux: %g, Area: %g, Correction: %g\n", *ibm_Flux, *ibm_Area, correction);
2834 if (NumberOfBodies>1){
2835 for (ibi=0; ibi<NumberOfBodies; ibi++) LOG_ALLOW(GLOBAL, LOG_INFO, " [Body %d] Uncorrected Flux: %g, Area: %g, Correction: %g\n", ibi, IB_Flux[ibi], IB_Area[ibi], Correction[ibi]);
2836 }
2837
2838 //================================================================================
2839 // PASS 2: Apply Correction to Velocity Field
2840 // This pass modifies the velocity at the boundary to enforce zero net flux.
2841 //================================================================================
2842 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 2: Applying velocity corrections at the boundary.\n");
2843
2844 for (k=lzs; k<lze; k++) {
2845 for (j=lys; j<lye; j++) {
2846 for (i=lxs; i<lxe; i++) {
2847 if (nvert[k][j][i] < 0.1) {
2848 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] <ibmval && i < xend) {
2849 if (fabs(ucor[k][j][i].x)>epsilon){
2850 if (flg==3)
2851 ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x)/
2852 sqrt(csi[k][j][i].x * csi[k][j][i].x +
2853 csi[k][j][i].y * csi[k][j][i].y +
2854 csi[k][j][i].z * csi[k][j][i].z);
2855 else if (flg==2)
2856 ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x);
2857 else if (NumberOfBodies > 1) {
2858 ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
2859 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
2860 csi[k][j][i].y * csi[k][j][i].y +
2861 csi[k][j][i].z * csi[k][j][i].z) *
2862 Correction[ibi];
2863 }
2864 else
2865 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
2866 csi[k][j][i].y * csi[k][j][i].y +
2867 csi[k][j][i].z * csi[k][j][i].z) *
2868 correction;
2869 }
2870 }
2871 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
2872 if (fabs(ucor[k][j][i].y)>epsilon) {
2873 if (flg==3)
2874 ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y)/
2875 sqrt(eta[k][j][i].x * eta[k][j][i].x +
2876 eta[k][j][i].y * eta[k][j][i].y +
2877 eta[k][j][i].z * eta[k][j][i].z);
2878 else if (flg==2)
2879 ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y);
2880 else if (NumberOfBodies > 1) {
2881 ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
2882 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
2883 eta[k][j][i].y * eta[k][j][i].y +
2884 eta[k][j][i].z * eta[k][j][i].z) *
2885 Correction[ibi];
2886 }
2887 else
2888 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
2889 eta[k][j][i].y * eta[k][j][i].y +
2890 eta[k][j][i].z * eta[k][j][i].z) *
2891 correction;
2892 }
2893 }
2894 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
2895 if (fabs(ucor[k][j][i].z)>epsilon) {
2896 if (flg==3)
2897 ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z)/
2898 sqrt(zet[k][j][i].x * zet[k][j][i].x +
2899 zet[k][j][i].y * zet[k][j][i].y +
2900 zet[k][j][i].z * zet[k][j][i].z);
2901 else if (flg==2)
2902 ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z);
2903 else if (NumberOfBodies > 1) {
2904 ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
2905 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
2906 zet[k][j][i].y * zet[k][j][i].y +
2907 zet[k][j][i].z * zet[k][j][i].z) *
2908 Correction[ibi];
2909 }
2910 else
2911 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
2912 zet[k][j][i].y * zet[k][j][i].y +
2913 zet[k][j][i].z * zet[k][j][i].z) *
2914 correction;
2915 }
2916 }
2917 }
2918
2919 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
2920 if (nvert[k][j][i+1] < 0.1 && i < xend) {
2921 if (fabs(ucor[k][j][i].x)>epsilon) {
2922 if (flg==3)
2923 ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x)/
2924 sqrt(csi[k][j][i].x * csi[k][j][i].x +
2925 csi[k][j][i].y * csi[k][j][i].y +
2926 csi[k][j][i].z * csi[k][j][i].z);
2927 else if (flg==2)
2928 ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x);
2929 else if (NumberOfBodies > 1) {
2930 ibi=(int)((nvert[k][j][i]-1.0)*1001);
2931 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2932 csi[k][j][i].y * csi[k][j][i].y +
2933 csi[k][j][i].z * csi[k][j][i].z) *
2934 Correction[ibi];
2935 }
2936 else
2937 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2938 csi[k][j][i].y * csi[k][j][i].y +
2939 csi[k][j][i].z * csi[k][j][i].z) *
2940 correction;
2941 }
2942 }
2943 if (nvert[k][j+1][i] < 0.1 && j < yend) {
2944 if (fabs(ucor[k][j][i].y)>epsilon) {
2945 if (flg==3)
2946 ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y)/
2947 sqrt(eta[k][j][i].x * eta[k][j][i].x +
2948 eta[k][j][i].y * eta[k][j][i].y +
2949 eta[k][j][i].z * eta[k][j][i].z);
2950 else if (flg==2)
2951 ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y);
2952 else if (NumberOfBodies > 1) {
2953 ibi=(int)((nvert[k][j][i]-1.0)*1001);
2954 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2955 eta[k][j][i].y * eta[k][j][i].y +
2956 eta[k][j][i].z * eta[k][j][i].z) *
2957 Correction[ibi];
2958 }
2959 else
2960 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2961 eta[k][j][i].y * eta[k][j][i].y +
2962 eta[k][j][i].z * eta[k][j][i].z) *
2963 correction;
2964 }
2965 }
2966 if (nvert[k+1][j][i] < 0.1 && k < zend) {
2967 if (fabs(ucor[k][j][i].z)>epsilon) {
2968 if (flg==3)
2969 ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z)/
2970 sqrt(zet[k][j][i].x * zet[k][j][i].x +
2971 zet[k][j][i].y * zet[k][j][i].y +
2972 zet[k][j][i].z * zet[k][j][i].z);
2973 else if (flg==2)
2974 ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z);
2975 else if (NumberOfBodies > 1) {
2976 ibi=(int)((nvert[k][j][i]-1.0)*1001);
2977 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2978 zet[k][j][i].y * zet[k][j][i].y +
2979 zet[k][j][i].z * zet[k][j][i].z) *
2980 Correction[ibi];
2981 }
2982 else
2983 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2984 zet[k][j][i].y * zet[k][j][i].y +
2985 zet[k][j][i].z * zet[k][j][i].z) *
2986 correction;
2987 }
2988 }
2989 }
2990
2991 }
2992 }
2993 }
2994
2995 //================================================================================
2996 // PASS 3: Verification
2997 // This optional pass recalculates the flux to confirm the correction was successful.
2998 //================================================================================
2999 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 3: Verifying corrected flux.\n");
3000
3001 libm_Flux = 0;
3002 libm_area = 0;
3003 for (k=lzs; k<lze; k++) {
3004 for (j=lys; j<lye; j++) {
3005 for (i=lxs; i<lxe; i++) {
3006 if (nvert[k][j][i] < 0.1) {
3007 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
3008 libm_Flux += ucor[k][j][i].x;
3009 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3010 csi[k][j][i].y * csi[k][j][i].y +
3011 csi[k][j][i].z * csi[k][j][i].z);
3012
3013 }
3014 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3015 libm_Flux += ucor[k][j][i].y;
3016 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3017 eta[k][j][i].y * eta[k][j][i].y +
3018 eta[k][j][i].z * eta[k][j][i].z);
3019 }
3020 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3021 libm_Flux += ucor[k][j][i].z;
3022 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3023 zet[k][j][i].y * zet[k][j][i].y +
3024 zet[k][j][i].z * zet[k][j][i].z);
3025 }
3026 }
3027
3028 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3029 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3030 libm_Flux -= ucor[k][j][i].x;
3031 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3032 csi[k][j][i].y * csi[k][j][i].y +
3033 csi[k][j][i].z * csi[k][j][i].z);
3034
3035 }
3036 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3037 libm_Flux -= ucor[k][j][i].y;
3038 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3039 eta[k][j][i].y * eta[k][j][i].y +
3040 eta[k][j][i].z * eta[k][j][i].z);
3041 }
3042 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3043 libm_Flux -= ucor[k][j][i].z;
3044 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3045 zet[k][j][i].y * zet[k][j][i].y +
3046 zet[k][j][i].z * zet[k][j][i].z);
3047 }
3048 }
3049
3050 }
3051 }
3052 }
3053
3054 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3055 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3056
3057 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
3058/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
3059 LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Corrected (Verified) Flux: %g, Area: %g\n", *ibm_Flux, *ibm_Area);
3060
3061
3062 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
3063 if (xe==mx){
3064 i=mx-2;
3065 for (k=lzs; k<lze; k++) {
3066 for (j=lys; j<lye; j++) {
3067 // if(j>0 && k>0 && j<user->JM && k<user->KM){
3068 if ((nvert[k][j][i]>ibmval && nvert[k][j][i+1]<0.1) || (nvert[k][j][i]<0.1 && nvert[k][j][i+1]>ibmval)) ucor[k][j][i].x=0.0;
3069
3070 // }
3071 }
3072 }
3073 }
3074 }
3075
3076 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
3077 if (ye==my){
3078 j=my-2;
3079 for (k=lzs; k<lze; k++) {
3080 for (i=lxs; i<lxe; i++) {
3081 // if(i>0 && k>0 && i<user->IM && k<user->KM){
3082 if ((nvert[k][j][i]>ibmval && nvert[k][j+1][i]<0.1) || (nvert[k][j][i]<0.1 && nvert[k][j+1][i]>ibmval)) ucor[k][j][i].y=0.0;
3083 // }
3084 }
3085 }
3086 }
3087 }
3088
3089 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
3090 if (ze==mz){
3091 k=mz-2;
3092 for (j=lys; j<lye; j++) {
3093 for (i=lxs; i<lxe; i++) {
3094 // if(i>0 && j>0 && i<user->IM && j<user->JM){
3095 if ((nvert[k][j][i]>ibmval && nvert[k+1][j][i]<0.1) || (nvert[k][j][i]<0.1 && nvert[k+1][j][i]>ibmval)) ucor[k][j][i].z=0.0;
3096 // }
3097 }
3098 }
3099 }
3100 }
3101
3102
3103 DMDAVecRestoreArray(da, user->lNvert, &nvert);
3104 DMDAVecRestoreArray(fda, user->lCsi, &csi);
3105 DMDAVecRestoreArray(fda, user->lEta, &eta);
3106 DMDAVecRestoreArray(fda, user->lZet, &zet);
3107 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3108
3109 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3110 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3111
3112 /* periodci boundary condition update corrected flux */
3113 //Mohsen Dec 2015
3114 DMDAVecGetArray(fda, user->lUcont, &lucor);
3115 DMDAVecGetArray(fda, user->Ucont, &ucor);
3116
3117 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
3118 if (xs==0){
3119 i=xs;
3120 for (k=zs; k<ze; k++) {
3121 for (j=ys; j<ye; j++) {
3122 if(j>0 && k>0 && j<user->JM && k<user->KM){
3123 ucor[k][j][i].x=lucor[k][j][i-2].x;
3124 }
3125 }
3126 }
3127 }
3128 }
3129 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
3130 if (ys==0){
3131 j=ys;
3132 for (k=zs; k<ze; k++) {
3133 for (i=xs; i<xe; i++) {
3134 if(i>0 && k>0 && i<user->IM && k<user->KM){
3135 ucor[k][j][i].y=lucor[k][j-2][i].y;
3136 }
3137 }
3138 }
3139 }
3140 }
3141 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
3142 if (zs==0){
3143 k=zs;
3144 for (j=ys; j<ye; j++) {
3145 for (i=xs; i<xe; i++) {
3146 if(i>0 && j>0 && i<user->IM && j<user->JM){
3147 ucor[k][j][i].z=lucor[k-2][j][i].z;
3148 }
3149 }
3150 }
3151 }
3152 }
3153 DMDAVecRestoreArray(fda, user->lUcont, &lucor);
3154 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3155
3156 /* DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); */
3157/* DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); */
3158 DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont);
3159 DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont);
3160
3161 if (NumberOfBodies > 1) {
3162 free(lIB_Flux);
3163 free(lIB_area);
3164 free(IB_Flux);
3165 free(IB_Area);
3166 free(Correction);
3167 }
3168
3169 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting VolumeFlux.\n");
3170
3171 return 0;
3172}
3173
3174PetscErrorCode FullyBlocked(UserCtx *user)
3175{
3176 DM da = user->da;
3177 Vec nNvert;
3178 DMDALocalInfo info = user->info;
3179
3180 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3181
3182 PetscInt i, j, k;
3183
3184 PetscInt *KSKE = user->KSKE;
3185 PetscReal ***nvert;
3186 PetscBool *Blocked;
3187
3188 DMDACreateNaturalVector(da, &nNvert);
3189 DMDAGlobalToNaturalBegin(da, user->Nvert, INSERT_VALUES, nNvert);
3190 DMDAGlobalToNaturalEnd(da, user->Nvert, INSERT_VALUES, nNvert);
3191
3192 VecScatter ctx;
3193 Vec Zvert;
3194 VecScatterCreateToZero(nNvert, &ctx, &Zvert);
3195
3196 VecScatterBegin(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3197 VecScatterEnd(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3198
3199 VecScatterDestroy(&ctx);
3200 VecDestroy(&nNvert);
3201
3202 PetscInt rank;
3203 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
3204
3205 if (!rank) {
3206
3207 VecGetArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3208 PetscMalloc(mx*my*sizeof(PetscBool), &Blocked);
3209 for (j=1; j<my-1; j++) {
3210 for (i=1; i<mx-1; i++) {
3211 Blocked[j*mx+i] = PETSC_FALSE;
3212 for (k=0; k<mz; k++) {
3213 if (nvert[k][j][i] > 0.1) {
3214 if (!Blocked[j*mx+i]) {
3215 KSKE[2*(j*mx+i)] = k;
3216 Blocked[j*mx+i] = PETSC_TRUE;
3217 }
3218 else {
3219 KSKE[2*(j*mx+i)] = PetscMin(KSKE[2*(j*mx+i)], k);
3220 }
3221 }
3222 }
3223 }
3224 }
3225
3226
3227 user->multinullspace = PETSC_TRUE;
3228 for (j=1; j<my-1; j++) {
3229 for (i=1; i<mx-1; i++) {
3230 if (!Blocked[j*mx+i]) {
3231 user->multinullspace = PETSC_FALSE;
3232 break;
3233 }
3234 }
3235 }
3236 PetscFree(Blocked);
3237 VecRestoreArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3238 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3239 if (user->multinullspace) {
3240 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3241
3242 }
3243 }
3244 else {
3245 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3246 if (user->multinullspace) {
3247 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3248 }
3249 }
3250
3251
3252
3253 VecDestroy(&Zvert);
3254 return 0;
3255}
3256
3257PetscErrorCode MyNvertRestriction(UserCtx *user_h, UserCtx *user_c)
3258{
3259 // DA da = user_c->da, fda = user_c->fda;
3260
3261
3262
3263 DMDALocalInfo info = user_c->info;
3264 PetscInt xs = info.xs, xe = info.xs + info.xm;
3265 PetscInt ys = info.ys, ye = info.ys + info.ym;
3266 PetscInt zs = info.zs, ze = info.zs + info.zm;
3267 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3268
3269 PetscInt i,j,k;
3270 PetscInt ih, jh, kh, ia, ja, ka;
3271 PetscInt lxs, lxe, lys, lye, lzs, lze;
3272
3273 PetscReal ***nvert, ***nvert_h;
3274
3275 DMDAVecGetArray(user_h->da, user_h->lNvert, &nvert_h);
3276 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert);
3277
3278 lxs = xs; lxe = xe;
3279 lys = ys; lye = ye;
3280 lzs = zs; lze = ze;
3281
3282 if (xs==0) lxs = xs+1;
3283 if (ys==0) lys = ys+1;
3284 if (zs==0) lzs = zs+1;
3285
3286 if (xe==mx) lxe = xe-1;
3287 if (ye==my) lye = ye-1;
3288 if (ze==mz) lze = ze-1;
3289
3290 if ((user_c->isc)) ia = 0;
3291 else ia = 1;
3292
3293 if ((user_c->jsc)) ja = 0;
3294 else ja = 1;
3295
3296 if ((user_c->ksc)) ka = 0;
3297 else ka = 1;
3298
3299 VecSet(user_c->Nvert, 0.);
3300 if (user_c->thislevel > 0) {
3301 for (k=lzs; k<lze; k++) {
3302 for (j=lys; j<lye; j++) {
3303 for (i=lxs; i<lxe; i++) {
3304 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3305 if (nvert_h[kh ][jh ][ih ] *
3306 nvert_h[kh ][jh ][ih-ia] *
3307 nvert_h[kh ][jh-ja][ih ] *
3308 nvert_h[kh-ka][jh ][ih ] *
3309 nvert_h[kh ][jh-ja][ih-ia] *
3310 nvert_h[kh-ka][jh ][ih-ia] *
3311 nvert_h[kh-ka][jh-ja][ih ] *
3312 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3313 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3314 }
3315 }
3316 }
3317 }
3318 }
3319 else {
3320 for (k=lzs; k<lze; k++) {
3321 for (j=lys; j<lye; j++) {
3322 for (i=lxs; i<lxe; i++) {
3323 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3324 if (nvert_h[kh ][jh ][ih ] *
3325 nvert_h[kh ][jh ][ih-ia] *
3326 nvert_h[kh ][jh-ja][ih ] *
3327 nvert_h[kh-ka][jh ][ih ] *
3328 nvert_h[kh ][jh-ja][ih-ia] *
3329 nvert_h[kh-ka][jh ][ih-ia] *
3330 nvert_h[kh-ka][jh-ja][ih ] *
3331 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3332 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3333 }
3334 }
3335 }
3336 }
3337 }
3338 DMDAVecRestoreArray(user_h->da, user_h->lNvert, &nvert_h);
3339 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert);
3340
3341 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3342 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3343 //Mohsen Dec 2015
3344 DMDAVecGetArray(user_c->da, user_c->lNvert, &nvert);
3345 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert_h);
3346
3347 for (k=lzs; k<lze; k++) {
3348 for (j=lys; j<lye; j++) {
3349 for (i=lxs; i<lxe; i++) {
3350 if (nvert_h[k][j][i] < 0.1) {
3351 if (nvert[k][j][i+1] + nvert[k][j][i-1] > 1.1 &&
3352 nvert[k][j+1][i] + nvert[k][j-1][i] > 1.1 &&
3353 nvert[k+1][j][i] + nvert[k-1][j][i] > 1.1) {
3354 nvert_h[k][j][i] = 1.;
3355 }
3356 }
3357 }
3358 }
3359 }
3360
3361 DMDAVecRestoreArray(user_c->da, user_c->lNvert, &nvert);
3362 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert_h);
3363 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3364 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3365 /* DMLocalToGlobalBegin(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3366/* DMLocalToGlobalEnd(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3367 return 0;
3368}
3369
3370#undef __FUNCT__
3371#define __FUNCT__ "PoissonSolver_MG"
3372PetscErrorCode PoissonSolver_MG(UserMG *usermg)
3373{
3374 // --- CONTEXT ACQUISITION BLOCK ---
3375 // Get the master simulation context from the first block's UserCtx on the finest level.
3376 // This provides access to all former global variables.
3377 SimCtx *simCtx = usermg->mgctx[0].user[0].simCtx;
3378
3379 // Create local variables to mirror the legacy globals for minimal code changes.
3380 const PetscInt block_number = simCtx->block_number;
3381 const PetscInt immersed = simCtx->immersed;
3382 const PetscInt MHV = simCtx->MHV;
3383 const PetscInt LV = simCtx->LV;
3384 PetscMPIInt rank = simCtx->rank;
3385 // --- END CONTEXT ACQUISITION BLOCK ---
3386
3387 PetscErrorCode ierr;
3388 PetscInt l, bi;
3389 MGCtx *mgctx = usermg->mgctx;
3390 KSP mgksp, subksp;
3391 PC mgpc, subpc;
3392 UserCtx *user;
3393
3394 PetscFunctionBeginUser; // Moved to after variable declarations
3395 PROFILE_FUNCTION_BEGIN;
3396 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting Multigrid Poisson Solve...\n");
3397
3398 for (bi = 0; bi < block_number; bi++) {
3399
3400 // ====================================================================
3401 // SECTION: Immersed Boundary Specific Setup (Conditional)
3402 // ====================================================================
3403 if (immersed) {
3404 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Performing IBM pre-solve setup (Nvert restriction, etc.).\n", bi);
3405 for (l = usermg->mglevels - 1; l > 0; l--) {
3406 mgctx[l].user[bi].multinullspace = PETSC_FALSE;
3407 MyNvertRestriction(&mgctx[l].user[bi], &mgctx[l-1].user[bi]);
3408 }
3409 // Coarsest level check for disconnected domains due to IBM
3410 l = 0;
3411 user = mgctx[l].user;
3412 ierr = PetscMalloc1(user[bi].info.mx * user[bi].info.my * 2, &user[bi].KSKE); CHKERRQ(ierr);
3413 FullyBlocked(&user[bi]);
3414 }
3415
3416
3417 l = usermg->mglevels - 1;
3418 user = mgctx[l].user;
3419
3420 // We are solving the linear system AX=B where A = Laplacian Operator Matrix; X = Unknown Phi (Pressure Correction) and B = RHS (Flux Divergence based)
3421
3422 // --- 1. Compute RHS of the Poisson Equation ---
3423 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Computing Poisson RHS...\n", bi);
3424 ierr = VecDuplicate(user[bi].P, &user[bi].B); CHKERRQ(ierr);
3425
3426 PetscReal ibm_Flux, ibm_Area;
3427 PetscInt flg = immersed - 1;
3428
3429 // Calculate volume flux source terms (often from IBM)
3430 VolumeFlux(&user[bi], &ibm_Flux, &ibm_Area, flg);
3431 if (MHV || LV) {
3432 flg = ((MHV > 1 || LV) && bi == 0) ? 1 : 0;
3433 VolumeFlux_rev(&user[bi], &ibm_Flux, &ibm_Area, flg);
3434 }
3435 // Calculate the main flux divergence term B.
3436 PoissonRHS(&user[bi], user[bi].B);
3437
3438 // --- 2. Assemble LHS Matrix (Laplacian) on all MG levels ---
3439 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Assembling Poisson LHS on all levels...\n", bi);
3440 for (l = usermg->mglevels - 1; l >= 0; l--) {
3441 user = mgctx[l].user;
3442 LOG_ALLOW(GLOBAL,LOG_DEBUG," Calculating LHS for Level %d.\n",l);
3443 PoissonLHSNew(&user[bi]);
3444 }
3445
3446 // --- 3. Setup PETSc KSP and PCMG (Multigrid Preconditioner) ---
3447 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Configuring KSP and PCMG...\n", bi);
3448
3449 ierr = KSPCreate(PETSC_COMM_WORLD, &mgksp); CHKERRQ(ierr);
3450 ierr = KSPAppendOptionsPrefix(mgksp, "ps_"); CHKERRQ(ierr);
3451
3452 // =======================================================================
3453 DualMonitorCtx *monctx;
3454 char filen[128];
3455 PetscBool has_monitor_option;
3456
3457 // 1. Allocate the context and set it up.
3458 ierr = PetscNew(&monctx); CHKERRQ(ierr);
3459
3460 monctx->step = simCtx->step;
3461 monctx->block_id = bi;
3462 monctx->file_handle = NULL;
3463
3464 // Only rank 0 handles the file.
3465 if (!rank) {
3466 sprintf(filen, "%s/Poisson_Solver_Convergence_History_Block_%d.log", simCtx->log_dir,bi);
3467 // On the very first step of the entire simulation, TRUNCATE the file.
3468 if (simCtx->step == simCtx->StartStep + 1) {
3469 monctx->file_handle = fopen(filen, "w");
3470 } else { // For all subsequent steps, APPEND to the file.
3471 monctx->file_handle = fopen(filen, "a");
3472 }
3473
3474 if (monctx->file_handle) {
3475 PetscFPrintf(PETSC_COMM_SELF, monctx->file_handle, "--- Convergence for Timestep %d, Block %d ---\n", (int)simCtx->step, bi);
3476 } else {
3477 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open KSP monitor log file: %s", filen);
3478 }
3479 }
3480
3481 ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &has_monitor_option); CHKERRQ(ierr);
3482 monctx->log_to_console = has_monitor_option;
3483
3484 ierr = KSPMonitorSet(mgksp, DualKSPMonitor, monctx, DualMonitorDestroy); CHKERRQ(ierr);
3485 // =======================================================================
3486
3487 ierr = KSPGetPC(mgksp, &mgpc); CHKERRQ(ierr);
3488 ierr = PCSetType(mgpc, PCMG); CHKERRQ(ierr);
3489
3490 ierr = PCMGSetLevels(mgpc, usermg->mglevels, PETSC_NULL); CHKERRQ(ierr);
3491 ierr = PCMGSetCycleType(mgpc, PC_MG_CYCLE_V); CHKERRQ(ierr);
3492 ierr = PCMGSetType(mgpc, PC_MG_MULTIPLICATIVE); CHKERRQ(ierr);
3493
3494 // --- 4. Define Restriction and Interpolation Operators for MG ---
3495 for (l = usermg->mglevels - 1; l > 0; l--) {
3496
3497 // Get stable pointers directly from the main mgctx array.
3498 // These pointers point to memory that will persist.
3499 UserCtx *fine_user_ctx = &mgctx[l].user[bi];
3500 UserCtx *coarse_user_ctx = &mgctx[l-1].user[bi];
3501
3502 // --- Configure the context pointers ---
3503 // The coarse UserCtx needs to know about the fine grid for restriction.
3504 coarse_user_ctx->da_f = &(fine_user_ctx->da);
3505 coarse_user_ctx->user_f = fine_user_ctx;
3506
3507 // The fine UserCtx needs to know about the coarse grid for interpolation.
3508 fine_user_ctx->da_c = &(coarse_user_ctx->da);
3509 fine_user_ctx->user_c = coarse_user_ctx;
3510 fine_user_ctx->lNvert_c = &(coarse_user_ctx->lNvert);
3511
3512 // --- Get matrix dimensions ---
3513 PetscInt m_c = (coarse_user_ctx->info.xm * coarse_user_ctx->info.ym * coarse_user_ctx->info.zm);
3514 PetscInt m_f = (fine_user_ctx->info.xm * fine_user_ctx->info.ym * fine_user_ctx->info.zm);
3515 PetscInt M_c = (coarse_user_ctx->info.mx * coarse_user_ctx->info.my * coarse_user_ctx->info.mz);
3516 PetscInt M_f = (fine_user_ctx->info.mx * fine_user_ctx->info.my * fine_user_ctx->info.mz);
3517
3518 LOG_ALLOW(GLOBAL,LOG_DEBUG,"level = %d; m_c = %d; m_f = %d; M_c = %d; M_f = %d.\n",l,m_c,m_f,M_c,M_f);
3519 // --- Create the MatShell objects ---
3520 // Pass the STABLE pointer coarse_user_ctx as the context for restriction.
3521 ierr = MatCreateShell(PETSC_COMM_WORLD, m_c, m_f, M_c, M_f, (void*)coarse_user_ctx, &fine_user_ctx->MR); CHKERRQ(ierr);
3522
3523 // Pass the STABLE pointer fine_user_ctx as the context for interpolation.
3524 ierr = MatCreateShell(PETSC_COMM_WORLD, m_f, m_c, M_f, M_c, (void*)fine_user_ctx, &fine_user_ctx->MP); CHKERRQ(ierr);
3525
3526 // --- Set the operations for the MatShells ---
3527 ierr = MatShellSetOperation(fine_user_ctx->MR, MATOP_MULT, (void(*)(void))RestrictResidual_SolidAware); CHKERRQ(ierr);
3528 ierr = MatShellSetOperation(fine_user_ctx->MP, MATOP_MULT, (void(*)(void))MyInterpolation); CHKERRQ(ierr);
3529
3530 // --- Register the operators with PCMG ---
3531 ierr = PCMGSetRestriction(mgpc, l, fine_user_ctx->MR); CHKERRQ(ierr);
3532 ierr = PCMGSetInterpolation(mgpc, l, fine_user_ctx->MP); CHKERRQ(ierr);
3533
3534 }
3535
3536 // --- 5. Configure Solvers on Each MG Level ---
3537 for (l = usermg->mglevels - 1; l >= 0; l--) {
3538 user = mgctx[l].user;
3539 if (l > 0) { // Smoother for fine levels
3540 ierr = PCMGGetSmoother(mgpc, l, &subksp); CHKERRQ(ierr);
3541 } else { // Direct or iterative solver for the coarsest level
3542 ierr = PCMGGetCoarseSolve(mgpc, &subksp); CHKERRQ(ierr);
3543 ierr = KSPSetTolerances(subksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, 40); CHKERRQ(ierr);
3544 }
3545
3546 ierr = KSPSetOperators(subksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
3547 ierr = KSPSetFromOptions(subksp); CHKERRQ(ierr);
3548 ierr = KSPGetPC(subksp, &subpc); CHKERRQ(ierr);
3549 ierr = PCSetType(subpc, PCBJACOBI); CHKERRQ(ierr);
3550
3551 KSP *subsubksp;
3552 PC subsubpc;
3553 PetscInt nlocal;
3554
3555 // This logic is required for both the smoother and the coarse solve
3556 // since both use PCBJACOBI.
3557 ierr = KSPSetUp(subksp); CHKERRQ(ierr); // Set up KSP to allow access to sub-KSPs
3558 ierr = PCBJacobiGetSubKSP(subpc, &nlocal, NULL, &subsubksp); CHKERRQ(ierr);
3559
3560 for (PetscInt abi = 0; abi < nlocal; abi++) {
3561 ierr = KSPGetPC(subsubksp[abi], &subsubpc); CHKERRQ(ierr);
3562 // Add the critical shift amount
3563 ierr = PCFactorSetShiftAmount(subsubpc, 1.e-10); CHKERRQ(ierr);
3564 }
3565
3566 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &user[bi].nullsp); CHKERRQ(ierr);
3567 ierr = MatNullSpaceSetFunction(user[bi].nullsp, PoissonNullSpaceFunction, &user[bi]); CHKERRQ(ierr);
3568 ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
3569
3570 ierr = PCMGSetResidual(mgpc, l, PCMGResidualDefault, user[bi].A); CHKERRQ(ierr);
3571 ierr = KSPSetUp(subksp); CHKERRQ(ierr);
3572
3573 if (l < usermg->mglevels - 1) {
3574 ierr = MatCreateVecs(user[bi].A, &user[bi].R, PETSC_NULL); CHKERRQ(ierr);
3575 ierr = PCMGSetRhs(mgpc, l, user[bi].R); CHKERRQ(ierr);
3576 }
3577 }
3578
3579 // --- 6. Set Final KSP Operators and Solve ---
3580 l = usermg->mglevels - 1;
3581 user = mgctx[l].user;
3582
3583 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Setting KSP operators and solving...\n", bi);
3584 ierr = KSPSetOperators(mgksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
3585 ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
3586 ierr = KSPSetFromOptions(mgksp); CHKERRQ(ierr);
3587 ierr = KSPSetUp(mgksp); CHKERRQ(ierr);
3588 ierr = KSPSolve(mgksp, user[bi].B, user[bi].Phi); CHKERRQ(ierr);
3589
3590 // --- 7. Cleanup for this block ---
3591 for (l = usermg->mglevels - 1; l >= 0; l--) {
3592 user = mgctx[l].user;
3593 MatNullSpaceDestroy(&user[bi].nullsp);
3594 MatDestroy(&user[bi].A);
3595 user[bi].assignedA = PETSC_FALSE;
3596 if (l > 0) {
3597 MatDestroy(&user[bi].MR);
3598 MatDestroy(&user[bi].MP);
3599 } else if (l==0 && immersed) {
3600 PetscFree(user[bi].KSKE);
3601 }
3602 if (l < usermg->mglevels - 1) {
3603 VecDestroy(&user[bi].R);
3604 }
3605 }
3606
3607 KSPDestroy(&mgksp);
3608 VecDestroy(&mgctx[usermg->mglevels-1].user[bi].B);
3609
3610 } // End of loop over blocks
3611
3612 LOG_ALLOW(GLOBAL, LOG_INFO, "Multigrid Poisson Solve complete.\n");
3613 PROFILE_FUNCTION_END;
3614 PetscFunctionReturn(0);
3615}
RefreshBoundaryGhostCells
PetscErrorCode RefreshBoundaryGhostCells(UserCtx *user)
(Public) Orchestrates the "light" refresh of all boundary ghost cells after the projection step.
Definition Boundaries.c:3106
logging.h
Logging utilities and macros for PETSc-based applications.
DualMonitorDestroy
PetscErrorCode DualMonitorDestroy(void **ctx)
Destroys the DualMonitorCtx.
Definition logging.c:787
DualMonitorCtx::log_to_console
PetscBool log_to_console
Definition logging.h:58
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:746
DualMonitorCtx::step
PetscInt step
Definition logging.h:60
LOG_LOOP_ALLOW_EXACT
#define LOG_LOOP_ALLOW_EXACT(scope, level, var, val, fmt,...)
Logs a custom message if a variable equals a specific value.
Definition logging.h:348
DualKSPMonitor
PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
A custom KSP monitor that logs to a file and optionally to the console.
Definition logging.c:819
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:737
DualMonitorCtx::file_handle
FILE * file_handle
Definition logging.h:57
DualMonitorCtx::block_id
PetscInt block_id
Definition logging.h:61
DualMonitorCtx
Context for a dual-purpose KSP monitor.
Definition logging.h:56
TW
#define TW
Definition poisson.c:326
MyNvertRestriction
PetscErrorCode MyNvertRestriction(UserCtx *user_h, UserCtx *user_c)
Definition poisson.c:3257
SE
#define SE
Definition poisson.c:317
BN
#define BN
Definition poisson.c:321
PoissonNullSpaceFunction
PetscErrorCode PoissonNullSpaceFunction(MatNullSpace nullsp, Vec X, void *ctx)
The callback function for PETSc's MatNullSpace object.
Definition poisson.c:1099
PoissonLHSNew
PetscErrorCode PoissonLHSNew(UserCtx *user)
Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.
Definition poisson.c:1620
WP
#define WP
Definition poisson.c:309
Projection
PetscErrorCode Projection(UserCtx *user)
Performs the projection step to enforce an incompressible velocity field.
Definition poisson.c:368
VolumeFlux_rev
PetscErrorCode VolumeFlux_rev(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
A specialized version of VolumeFlux, likely for reversed normals.
Definition poisson.c:2304
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition poisson.c:40
mymatmultadd
static PetscErrorCode mymatmultadd(Mat mat, Vec v1, Vec v2)
Definition poisson.c:1086
SW
#define SW
Definition poisson.c:319
BS
#define BS
Definition poisson.c:323
NE
#define NE
Definition poisson.c:316
MyRestriction
PetscErrorCode MyRestriction(Mat A, Vec X, Vec F)
The callback function for the multigrid restriction operator (MatShell).
Definition poisson.c:1474
CP
#define CP
Definition poisson.c:306
BE
#define BE
Definition poisson.c:325
UpdatePressure
PetscErrorCode UpdatePressure(UserCtx *user)
Updates the pressure field and handles periodic boundary conditions.
Definition poisson.c:915
BP
#define BP
Definition poisson.c:313
GridInterpolation
#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user)
Definition poisson.c:5
BW
#define BW
Definition poisson.c:327
RestrictResidual_SolidAware
static PetscErrorCode RestrictResidual_SolidAware(Mat A, Vec X, Vec F)
Definition poisson.c:1399
VolumeFlux
PetscErrorCode VolumeFlux(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
Calculates the net flux across the immersed boundary surface.
Definition poisson.c:2546
TE
#define TE
Definition poisson.c:324
PoissonRHS
PetscErrorCode PoissonRHS(UserCtx *user, Vec B)
Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermedia...
Definition poisson.c:2223
CorrectChannelFluxProfile
PetscErrorCode CorrectChannelFluxProfile(UserCtx *user)
Enforces a constant volumetric flux profile along the entire length of a driven periodic channel.
Definition poisson.c:116
MyInterpolation
PetscErrorCode MyInterpolation(Mat A, Vec X, Vec F)
The callback function for the multigrid interpolation operator (MatShell).
Definition poisson.c:1292
TS
#define TS
Definition poisson.c:322
FullyBlocked
PetscErrorCode FullyBlocked(UserCtx *user)
Definition poisson.c:3174
NP
#define NP
Definition poisson.c:310
PoissonSolver_MG
PetscErrorCode PoissonSolver_MG(UserMG *usermg)
Solves the pressure-Poisson equation using a geometric multigrid method.
Definition poisson.c:3372
EP
#define EP
Definition poisson.c:308
TN
#define TN
Definition poisson.c:320
SP
#define SP
Definition poisson.c:311
TP
#define TP
Definition poisson.c:312
NW
#define NW
Definition poisson.c:318
GridRestriction
static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k, PetscInt *ih, PetscInt *jh, PetscInt *kh, UserCtx *user)
Definition poisson.c:60
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2177
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1157
SimCtx::MHV
PetscInt MHV
Definition variables.h:620
UserCtx::isc
PetscInt isc
Definition variables.h:741
PERIODIC
@ PERIODIC
Definition variables.h:219
MGCtx::user
UserCtx * user
Definition variables.h:474
SimCtx::rank
PetscMPIInt rank
Definition variables.h:588
SimCtx::fish
PetscInt fish
Definition variables.h:620
SimCtx::LV
PetscInt LV
Definition variables.h:620
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:746
SimCtx::block_number
PetscInt block_number
Definition variables.h:649
UserCtx::lIEta
Vec lIEta
Definition variables.h:778
UserCtx::KSKE
PetscInt * KSKE
Definition variables.h:768
SimCtx::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition variables.h:661
UserCtx::lIZet
Vec lIZet
Definition variables.h:778
UserCtx::user_f
UserCtx * user_f
Definition variables.h:793
SimCtx::channelz
PetscInt channelz
Definition variables.h:621
UserCtx::lNvert
Vec lNvert
Definition variables.h:755
UserCtx::Phi
Vec Phi
Definition variables.h:755
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:731
UserCtx::ksc
PetscInt ksc
Definition variables.h:741
UserCtx::KM
PetscInt KM
Definition variables.h:737
UserCtx::lZet
Vec lZet
Definition variables.h:776
UserCtx::assignedA
PetscBool assignedA
Definition variables.h:772
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:230
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:246
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:245
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:296
UserCtx::lIAj
Vec lIAj
Definition variables.h:778
UserCtx::_this
PetscInt _this
Definition variables.h:741
UserCtx::lKEta
Vec lKEta
Definition variables.h:780
SimCtx::dt
PetscReal dt
Definition variables.h:599
UserCtx::jsc
PetscInt jsc
Definition variables.h:741
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:779
UserCtx::Ucont
Vec Ucont
Definition variables.h:755
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:595
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::lPhi
Vec lPhi
Definition variables.h:755
UserCtx::da_c
DM * da_c
Definition variables.h:794
UserCtx::user_c
UserCtx * user_c
Definition variables.h:793
SimCtx::NumberOfBodies
PetscInt NumberOfBodies
Definition variables.h:640
UserCtx::lKZet
Vec lKZet
Definition variables.h:780
SimCtx::log_dir
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:609
UserCtx::da_f
DM * da_f
Definition variables.h:794
UserCtx::lJEta
Vec lJEta
Definition variables.h:779
UserCtx::lCsi
Vec lCsi
Definition variables.h:776
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:792
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:780
UserCtx::JM
PetscInt JM
Definition variables.h:737
UserMG::mglevels
PetscInt mglevels
Definition variables.h:481
UserCtx::lJZet
Vec lJZet
Definition variables.h:779
UserCtx::lUcont
Vec lUcont
Definition variables.h:755
SimCtx::step
PetscInt step
Definition variables.h:593
UserCtx::lAj
Vec lAj
Definition variables.h:776
UserCtx::lICsi
Vec lICsi
Definition variables.h:778
UserCtx::info
DMDALocalInfo info
Definition variables.h:735
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::IM
PetscInt IM
Definition variables.h:737
UserCtx::FluxIntpSum
PetscReal FluxIntpSum
Definition variables.h:752
UserCtx::lEta
Vec lEta
Definition variables.h:776
UserCtx::Nvert
Vec Nvert
Definition variables.h:755
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:484
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:295
UserCtx::multinullspace
PetscBool multinullspace
Definition variables.h:769
UserCtx::lJAj
Vec lJAj
Definition variables.h:779
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:704
UserCtx::lKAj
Vec lKAj
Definition variables.h:780
SimCtx::immersed
PetscInt immersed
Definition variables.h:614
COEF_TIME_ACCURACY
#define COEF_TIME_ACCURACY
Coefficient controlling the temporal accuracy scheme (e.g., 1.5 for 2nd Order Backward Difference).
Definition variables.h:57
UserCtx::lNvert_c
Vec * lNvert_c
Definition variables.h:795
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:204
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:206
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:205
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:206
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:204
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:205
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
MGCtx
Context for Multigrid operations.
Definition variables.h:473
SimCtx
The master context for the entire simulation.
Definition variables.h:585
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:728
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:480