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