PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
rhs.c
Go to the documentation of this file.
1#include "rhs.h"
2#include "verification_sources.h"
3
4#undef __FUNCT__
5#define __FUNCT__ "Convection"
6/**
7 * @brief Implementation of \ref Convection().
8 * @details Full API contract (arguments, ownership, side effects) is documented with
9 * the header declaration in `include/rhs.h`.
10 * @see Convection()
11 */
12
13PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
14{
15 PetscErrorCode ierr;
16
17 // --- CONTEXT ACQUISITION BLOCK ---
18 // Get the master simulation context from the UserCtx.
19 SimCtx *simCtx = user->simCtx;
20
21 // Create local variables to mirror the legacy globals for minimal code changes.
22 const LESModelType les = simCtx->les;
23 const PetscInt central = simCtx->central; // Get this from SimCtx now
24 // --- END CONTEXT ACQUISITION BLOCK ---
25
26 Cmpnts ***ucont, ***ucat;
27 DM da = user->da, fda = user->fda;
28 DMDALocalInfo info;
29 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
30 PetscInt mx, my, mz; // Dimensions in three directions
31 PetscInt i, j, k;
32 Vec Fp1, Fp2, Fp3;
33 Cmpnts ***fp1, ***fp2, ***fp3;
34 Cmpnts ***conv;
35
36 PetscReal ucon, up, um;
37 PetscReal coef = 0.125, innerblank=7.;
38
39 PetscInt lxs, lxe, lys, lye, lzs, lze;
40
41 PetscReal ***nvert,***aj;
42
43 PROFILE_FUNCTION_BEGIN;
44
45 DMDAGetLocalInfo(da, &info);
46 mx = info.mx; my = info.my; mz = info.mz;
47 xs = info.xs; xe = xs + info.xm;
48 ys = info.ys; ye = ys + info.ym;
49 zs = info.zs; ze = zs + info.zm;
50 ierr = PreparePeriodicQuickStencilFields(user, Ucat, user->lNvert); CHKERRQ(ierr);
51
52 DMDAVecGetArray(fda, Ucont, &ucont);
53 DMDAVecGetArray(fda, Ucat, &ucat);
54 DMDAVecGetArray(fda, Conv, &conv);
55 DMDAVecGetArray(da, user->lAj, &aj);
56
57 VecDuplicate(Ucont, &Fp1);
58 VecDuplicate(Ucont, &Fp2);
59 VecDuplicate(Ucont, &Fp3);
60
61 DMDAVecGetArray(fda, Fp1, &fp1);
62 DMDAVecGetArray(fda, Fp2, &fp2);
63 DMDAVecGetArray(fda, Fp3, &fp3);
64
65 DMDAVecGetArray(da, user->lNvert, &nvert);
66
67
68 /* We have two different sets of node: 1. grid node, the physical points
69 where grid lines intercross; 2. storage node, where we store variables.
70 All node without explicitly specified as "grid node" refers to
71 storage node.
72
73 The integer node is defined at cell center while half node refers to
74 the actual grid node. (The reason to choose this arrangement is we need
75 ghost node, which is half node away from boundaries, to specify boundary
76 conditions. By using this storage arrangement, the actual storage need
77 is (IM+1) * (JM + 1) * (KM+1) where IM, JM, & KM refer to the number of
78 grid nodes along i, j, k directions.)
79
80 DA, the data structure used to define the storage of 3D arrays, is defined
81 as mx * my * mz. mx = IM+1, my = JM+1, mz = KM+1.
82
83 Staggered grid arrangement is used in this solver.
84 Pressure is stored at interger node (hence the cell center) and volume
85 fluxes defined on the center of each surface of a given control volume
86 is stored on the cloest upper integer node. */
87
88 /* First we calculate the flux on cell surfaces. Stored on the upper integer
89 node. For example, along i direction, the flux are stored at node 0:mx-2*/
90
91 lxs = xs; lxe = xe;
92 lys = ys; lye = ye;
93 lzs = zs; lze = ze;
94
95 if (xs==0) lxs = xs+1;
96 if (ys==0) lys = ys+1;
97 if (zs==0) lzs = zs+1;
98
99 if (xe==mx) lxe=xe-1;
100 if (ye==my) lye=ye-1;
101 if (ze==mz) lze=ze-1;
102
103 VecSet(Conv, 0.0);
104
105 /* Calculating the convective terms on cell centers.
106 First calcualte the contribution from i direction
107 The flux is evaluated by QUICK scheme */
108
109 for (k=lzs; k<lze; k++){
110 for (j=lys; j<lye; j++){
111 for (i=lxs-1; i<lxe; i++){
112
113
114 ucon = ucont[k][j][i].x * 0.5;
115
116 up = ucon + fabs(ucon);
117 um = ucon - fabs(ucon);
118
119 if (i>0 && i<mx-2 &&
120 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
121 (nvert[k][j][i-1] < 0.1 || nvert[k][j][i-1]>innerblank)) { // interial nodes
122 if ((les || central)) {
123 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
124 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
125 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
126
127 } else {
128 fp1[k][j][i].x =
129 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
130 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
131 fp1[k][j][i].y =
132 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
133 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
134 fp1[k][j][i].z =
135 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
136 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
137 }
138 }
139 else if ((les || central) && (i==0 || i==mx-2) &&
140 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
141 (nvert[k][j][i ] < 0.1 || nvert[k][j][i ]>innerblank))
142 {
143 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
144 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
145 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
146 }
147 else if (i==0 ||(nvert[k][j][i-1] > 0.1) ) {
148 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==0 && (nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
149 fp1[k][j][i].x =
150 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
151 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
152 fp1[k][j][i].y =
153 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
154 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
155 fp1[k][j][i].z =
156 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
157 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
158 }else{
159 fp1[k][j][i].x =
160 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
161 up * (coef * (-ucat[k][j][i ].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
162 fp1[k][j][i].y =
163 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
164 up * (coef * (-ucat[k][j][i ].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
165 fp1[k][j][i].z =
166 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
167 up * (coef * (-ucat[k][j][i ].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
168 }
169 }
170 else if (i==mx-2 ||(nvert[k][j][i+1]) > 0.1) {
171 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2 &&(nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
172 fp1[k][j][i].x =
173 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
174 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
175 fp1[k][j][i].y =
176 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
177 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
178 fp1[k][j][i].z =
179 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
180 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
181 }else{
182 fp1[k][j][i].x =
183 um * (coef * (-ucat[k][j][i+1].x -2. * ucat[k][j][i+1].x +3. * ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
184 up * (coef * (-ucat[k][j][i-1].x -2. * ucat[k][j][i ].x +3. * ucat[k][j][i+1].x) +ucat[k][j][i ].x);
185 fp1[k][j][i].y =
186 um * (coef * (-ucat[k][j][i+1].y -2. * ucat[k][j][i+1].y +3. * ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
187 up * (coef * (-ucat[k][j][i-1].y -2. * ucat[k][j][i ].y +3. * ucat[k][j][i+1].y) +ucat[k][j][i ].y);
188 fp1[k][j][i].z =
189 um * (coef * (-ucat[k][j][i+1].z -2. * ucat[k][j][i+1].z +3. * ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
190 up * (coef * (-ucat[k][j][i-1].z -2. * ucat[k][j][i ].z +3. * ucat[k][j][i+1].z) +ucat[k][j][i ].z);
191 }
192 }
193 }
194 }
195 }
196
197 /* j direction */
198 for (k=lzs; k<lze; k++) {
199 for(j=lys-1; j<lye; j++) {
200 for(i=lxs; i<lxe; i++) {
201 ucon = ucont[k][j][i].y * 0.5;
202
203 up = ucon + fabs(ucon);
204 um = ucon - fabs(ucon);
205
206 if (j>0 && j<my-2 &&
207 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i] > innerblank) &&
208 (nvert[k][j-1][i] < 0.1 || nvert[k][j-1][i] > innerblank)) {
209 if ((les || central)) {
210 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
211 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
212 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
213
214 } else {
215 fp2[k][j][i].x =
216 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
217 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
218 fp2[k][j][i].y =
219 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
220 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
221 fp2[k][j][i].z =
222 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
223 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
224 }
225 }
226 else if ((les || central) && (j==0 || j==my-2) &&
227 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i]>innerblank) &&
228 (nvert[k][j ][i] < 0.1 || nvert[k][j ][i]>innerblank))
229 {
230 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
231 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
232 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
233 }
234 else if (j==0 || (nvert[k][j-1][i]) > 0.1) {
235 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==0 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12 //
236 fp2[k][j][i].x =
237 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
238 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
239 fp2[k][j][i].y =
240 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
241 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
242 fp2[k][j][i].z =
243 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
244 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
245 }else{
246 fp2[k][j][i].x =
247 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
248 up * (coef * (-ucat[k][j ][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
249 fp2[k][j][i].y =
250 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
251 up * (coef * (-ucat[k][j ][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
252 fp2[k][j][i].z =
253 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
254 up * (coef * (-ucat[k][j ][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
255 }
256 }
257 else if (j==my-2 ||(nvert[k][j+1][i]) > 0.1) {
258 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12//
259 fp2[k][j][i].x =
260 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
261 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
262 fp2[k][j][i].y =
263 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
264 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
265 fp2[k][j][i].z =
266 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
267 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
268 }else{
269 fp2[k][j][i].x =
270 um * (coef * (-ucat[k][j+1][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
271 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
272 fp2[k][j][i].y =
273 um * (coef * (-ucat[k][j+1][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
274 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j][i ].y);
275 fp2[k][j][i].z =
276 um * (coef * (-ucat[k][j+1][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
277 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j][i ].z);
278 }
279 }
280 }
281 }
282 }
283
284
285 /* k direction */
286 for (k=lzs-1; k<lze; k++) {
287 for(j=lys; j<lye; j++) {
288 for(i=lxs; i<lxe; i++) {
289 ucon = ucont[k][j][i].z * 0.5;
290
291 up = ucon + fabs(ucon);
292 um = ucon - fabs(ucon);
293
294 if (k>0 && k<mz-2 &&
295 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i] > innerblank) &&
296 (nvert[k-1][j][i] < 0.1 || nvert[k-1][j][i] > innerblank)) {
297 if ((les || central)) {
298 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
299 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
300 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
301
302 } else {
303 fp3[k][j][i].x =
304 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
305 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
306 fp3[k][j][i].y =
307 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
308 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
309 fp3[k][j][i].z =
310 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
311 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
312 }
313 }
314 else if ((les || central) && (k==0 || k==mz-2) &&
315 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i]>innerblank) &&
316 (nvert[k ][j][i] < 0.1 || nvert[k ][j][i]>innerblank))
317 {
318 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
319 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
320 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
321 }
322 else if (k<mz-2 && (k==0 ||(nvert[k-1][j][i]) > 0.1)) {
323 if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==0 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
324 fp3[k][j][i].x =
325 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
326 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
327 fp3[k][j][i].y =
328 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
329 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
330 fp3[k][j][i].z =
331 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
332 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
333 }else{
334 fp3[k][j][i].x =
335 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
336 up * (coef * (-ucat[k ][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k][j][i ].x);
337 fp3[k][j][i].y =
338 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
339 up * (coef * (-ucat[k ][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k][j][i ].y);
340 fp3[k][j][i].z =
341 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
342 up * (coef * (-ucat[k ][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k][j][i ].z);
343 }
344 }
345 else if (k>0 && (k==mz-2 ||(nvert[k+1][j][i]) > 0.1)) {
346 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
347 fp3[k][j][i].x =
348 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
349 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
350 fp3[k][j][i].y =
351 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
352 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
353 fp3[k][j][i].z =
354 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
355 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
356 }else{
357 fp3[k][j][i].x =
358 um * (coef * (-ucat[k+1][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
359 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k][j][i ].x);
360 fp3[k][j][i].y =
361 um * (coef * (-ucat[k+1][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
362 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k][j][i ].y);
363 fp3[k][j][i].z =
364 um * (coef * (-ucat[k+1][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
365 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k][j][i ].z);
366 }
367 }
368 }
369 }
370 }
371
372 /* Calculate the convective terms under cartesian coordinates */
373
374 for (k=lzs; k<lze; k++) {
375 for (j=lys; j<lye; j++) {
376 for (i=lxs; i<lxe; i++) {
377 conv[k][j][i].x =
378 fp1[k][j][i].x - fp1[k][j][i-1].x +
379 fp2[k][j][i].x - fp2[k][j-1][i].x +
380 fp3[k][j][i].x - fp3[k-1][j][i].x;
381
382 conv[k][j][i].y =
383 fp1[k][j][i].y - fp1[k][j][i-1].y +
384 fp2[k][j][i].y - fp2[k][j-1][i].y +
385 fp3[k][j][i].y - fp3[k-1][j][i].y;
386
387 conv[k][j][i].z =
388 fp1[k][j][i].z - fp1[k][j][i-1].z +
389 fp2[k][j][i].z - fp2[k][j-1][i].z +
390 fp3[k][j][i].z - fp3[k-1][j][i].z;
391 }
392 }
393 }
394 /* for (k=zs; k<ze; k++) { */
395/* for (j=ys; j<ye; j++) { */
396/* for (i=xs; i<xe; i++) { */
397/* if (i==1 && (j==1) && (k==1 || k==21 || k==22|| k==200)) */
398/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d conv.y is %.15le conv.z is %.15le \n",i,j,k,conv[k][j][i].y,conv[k][j][i].z); */
399/* } */
400/* } */
401/* } */
402
403 DMDAVecRestoreArray(fda, Ucont, &ucont);
404 DMDAVecRestoreArray(fda, Ucat, &ucat);
405 DMDAVecRestoreArray(fda, Conv, &conv);
406 DMDAVecRestoreArray(da, user->lAj, &aj);
407
408 DMDAVecRestoreArray(fda, Fp1, &fp1);
409 DMDAVecRestoreArray(fda, Fp2, &fp2);
410 DMDAVecRestoreArray(fda, Fp3, &fp3);
411 DMDAVecRestoreArray(da, user->lNvert, &nvert);
412
413 VecDestroy(&Fp1);
414 VecDestroy(&Fp2);
415 VecDestroy(&Fp3);
416
417
418 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Convective term calculated .\n");
419
420 PROFILE_FUNCTION_END;
421 return (0);
422}
423
424
425#undef __FUNCT__
426#define __FUNCT__ "Viscous"
427/**
428 * @brief Implementation of \ref Viscous().
429 * @details Full API contract (arguments, ownership, side effects) is documented with
430 * the header declaration in `include/rhs.h`.
431 * @see Viscous()
432 */
433
434PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
435{
436
437 Vec Csi = user->lCsi, Eta = user->lEta, Zet = user->lZet;
438
439 Cmpnts ***ucont, ***ucat;
440
441 Cmpnts ***csi, ***eta, ***zet;
442 Cmpnts ***icsi, ***ieta, ***izet;
443 Cmpnts ***jcsi, ***jeta, ***jzet;
444 Cmpnts ***kcsi, ***keta, ***kzet;
445
446 PetscReal ***nvert;
447
448 DM da = user->da, fda = user->fda;
449 DMDALocalInfo info;
450 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
451 PetscInt mx, my, mz; // Dimensions in three directions
452 PetscInt i, j, k;
453 Vec Fp1, Fp2, Fp3;
454 Cmpnts ***fp1, ***fp2, ***fp3;
455 Cmpnts ***visc;
456 PetscReal ***aj, ***iaj, ***jaj, ***kaj;
457
458 PetscInt lxs, lxe, lys, lye, lzs, lze;
459
460 PetscReal ajc;
461
462 PetscReal dudc, dude, dudz, dvdc, dvde, dvdz, dwdc, dwde, dwdz;
463 PetscReal csi0, csi1, csi2, eta0, eta1, eta2, zet0, zet1, zet2;
464 PetscReal g11, g21, g31;
465 PetscReal r11, r21, r31, r12, r22, r32, r13, r23, r33;
466
467 PetscScalar solid,innerblank;
468
469 // --- CONTEXT ACQUISITION BLOCK ---
470 // Get the master simulation context from the UserCtx.
471 SimCtx *simCtx = user->simCtx;
472
473 // Create local variables to mirror the legacy globals for minimal code changes.
474 const LESModelType les = simCtx->les;
475 const PetscInt rans = simCtx->rans;
476 const PetscInt ti = simCtx->step; // Assuming simCtx->step is the new integer time counter
477 const PetscReal ren = simCtx->ren;
478 const PetscInt clark = simCtx->clark;
479 solid = 0.5;
480 innerblank = 7.;
481
482 PROFILE_FUNCTION_BEGIN;
483
484 DMDAVecGetArray(fda, Ucont, &ucont);
485 DMDAVecGetArray(fda, Ucat, &ucat);
486 DMDAVecGetArray(fda, Visc, &visc);
487
488 DMDAVecGetArray(fda, Csi, &csi);
489 DMDAVecGetArray(fda, Eta, &eta);
490 DMDAVecGetArray(fda, Zet, &zet);
491
492 DMDAVecGetArray(fda, user->lICsi, &icsi);
493 DMDAVecGetArray(fda, user->lIEta, &ieta);
494 DMDAVecGetArray(fda, user->lIZet, &izet);
495
496 DMDAVecGetArray(fda, user->lJCsi, &jcsi);
497 DMDAVecGetArray(fda, user->lJEta, &jeta);
498 DMDAVecGetArray(fda, user->lJZet, &jzet);
499
500 DMDAVecGetArray(fda, user->lKCsi, &kcsi);
501 DMDAVecGetArray(fda, user->lKEta, &keta);
502 DMDAVecGetArray(fda, user->lKZet, &kzet);
503
504 DMDAVecGetArray(da, user->lNvert, &nvert);
505
506 VecDuplicate(Ucont, &Fp1);
507 VecDuplicate(Ucont, &Fp2);
508 VecDuplicate(Ucont, &Fp3);
509
510 DMDAVecGetArray(fda, Fp1, &fp1);
511 DMDAVecGetArray(fda, Fp2, &fp2);
512 DMDAVecGetArray(fda, Fp3, &fp3);
513
514 DMDAVecGetArray(da, user->lAj, &aj);
515
516 DMDAGetLocalInfo(da, &info);
517
518 mx = info.mx; my = info.my; mz = info.mz;
519 xs = info.xs; xe = xs + info.xm;
520 ys = info.ys; ye = ys + info.ym;
521 zs = info.zs; ze = zs + info.zm;
522
523 /* First we calculate the flux on cell surfaces. Stored on the upper integer
524 node. For example, along i direction, the flux are stored at node 0:mx-2*/
525 lxs = xs; lxe = xe;
526 lys = ys; lye = ye;
527 lzs = zs; lze = ze;
528
529 if (xs==0) lxs = xs+1;
530 if (ys==0) lys = ys+1;
531 if (zs==0) lzs = zs+1;
532
533
534 if (xe==mx) lxe=xe-1;
535 if (ye==my) lye=ye-1;
536 if (ze==mz) lze=ze-1;
537
538 VecSet(Visc,0.0);
539
540 PetscReal ***lnu_t;
541
542 if(les) {
543 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
544 } else if (rans) {
545
546 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
547 }
548
549 /* The visc flux on each surface center is stored at previous integer node */
550
551 DMDAVecGetArray(da, user->lIAj, &iaj);
552 /* for (k=zs; k<ze; k++) { */
553/* for (j=ys; j<ye; j++) { */
554/* for (i=xs; i<xe; i++) { */
555/* if (i==1 && (j==0 ||j==1 || j==2) && (k==21 || k==22|| k==20)) */
556/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d u is %.15le v is %.15le w is %.15le \n",i,j,k,ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z ); */
557/* } */
558/* } */
559/* } */
560 // i direction
561 for (k=lzs; k<lze; k++) {
562 for (j=lys; j<lye; j++) {
563 for (i=lxs-1; i<lxe; i++) {
564
565 dudc = ucat[k][j][i+1].x - ucat[k][j][i].x;
566 dvdc = ucat[k][j][i+1].y - ucat[k][j][i].y;
567 dwdc = ucat[k][j][i+1].z - ucat[k][j][i].z;
568
569 if ((nvert[k][j+1][i ]> solid && nvert[k][j+1][i ]<innerblank) ||
570 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
571 dude = (ucat[k][j ][i+1].x + ucat[k][j ][i].x -
572 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.5;
573 dvde = (ucat[k][j ][i+1].y + ucat[k][j ][i].y -
574 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.5;
575 dwde = (ucat[k][j ][i+1].z + ucat[k][j ][i].z -
576 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.5;
577 }
578 else if ((nvert[k][j-1][i ]> solid && nvert[k][j-1][i ]<innerblank) ||
579 (nvert[k][j-1][i+1]> solid && nvert[k][j-1][i+1]<innerblank)) {
580 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
581 ucat[k][j ][i+1].x - ucat[k][j ][i].x) * 0.5;
582 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
583 ucat[k][j ][i+1].y - ucat[k][j ][i].y) * 0.5;
584 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
585 ucat[k][j ][i+1].z - ucat[k][j ][i].z) * 0.5;
586 }
587 else {
588 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
589 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.25;
590 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
591 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.25;
592 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
593 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.25;
594 }
595
596 if ((nvert[k+1][j][i ]> solid && nvert[k+1][j][i ]<innerblank)||
597 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
598 dudz = (ucat[k ][j][i+1].x + ucat[k ][j][i].x -
599 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.5;
600 dvdz = (ucat[k ][j][i+1].y + ucat[k ][j][i].y -
601 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.5;
602 dwdz = (ucat[k ][j][i+1].z + ucat[k ][j][i].z -
603 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.5;
604 }
605 else if ((nvert[k-1][j][i ]> solid && nvert[k-1][j][i ]<innerblank) ||
606 (nvert[k-1][j][i+1]> solid && nvert[k-1][j][i+1]<innerblank)) {
607
608 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
609 ucat[k ][j][i+1].x - ucat[k ][j][i].x) * 0.5;
610 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
611 ucat[k ][j][i+1].y - ucat[k ][j][i].y) * 0.5;
612 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
613 ucat[k ][j][i+1].z - ucat[k ][j][i].z) * 0.5;
614 }
615 else {
616 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
617 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.25;
618 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
619 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.25;
620 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
621 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.25;
622 }
623
624 csi0 = icsi[k][j][i].x;
625 csi1 = icsi[k][j][i].y;
626 csi2 = icsi[k][j][i].z;
627
628 eta0 = ieta[k][j][i].x;
629 eta1 = ieta[k][j][i].y;
630 eta2 = ieta[k][j][i].z;
631
632 zet0 = izet[k][j][i].x;
633 zet1 = izet[k][j][i].y;
634 zet2 = izet[k][j][i].z;
635
636 g11 = csi0 * csi0 + csi1 * csi1 + csi2 * csi2;
637 g21 = eta0 * csi0 + eta1 * csi1 + eta2 * csi2;
638 g31 = zet0 * csi0 + zet1 * csi1 + zet2 * csi2;
639
640 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
641 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
642 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
643
644 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
645 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
646 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
647
648 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
649 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
650 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
651
652 ajc = iaj[k][j][i];
653
654 double nu = 1./ren, nu_t=0;
655
656 if( les || (rans && ti>0) ) {
657 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j][i+1]) ), 2.0) * Sabs;
658 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j][i+1]);
659 if ( (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == WALL && i==0) || (user->boundary_faces[BC_FACE_POS_X].mathematical_type == WALL && i==mx-2) ) nu_t=0;
660 fp1[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * csi0 + r21 * csi1 + r31 * csi2) * ajc * (nu_t);
661 fp1[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * csi0 + r22 * csi1 + r32 * csi2) * ajc * (nu_t);
662 fp1[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * csi0 + r23 * csi1 + r33 * csi2) * ajc * (nu_t);
663 }
664 else {
665 fp1[k][j][i].x = 0;
666 fp1[k][j][i].y = 0;
667 fp1[k][j][i].z = 0;
668 }
669
670 fp1[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * csi0 + r21 * csi1 + r31 * csi2 ) * ajc * (nu);
671 fp1[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * csi0 + r22 * csi1 + r32 * csi2 ) * ajc * (nu);
672 fp1[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * csi0 + r23 * csi1 + r33 * csi2 ) * ajc * (nu);
673
674
675 if(clark) {
676 double dc, de, dz;
677 ComputeCellCharacteristicLengthScale (ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
678 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
679
680 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
681 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
682 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
683 double t21 = t12;
684 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
685 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
686 double t31 = t13;
687 double t32 = t23;
688 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
689
690 fp1[k][j][i].x -= ( t11 * csi0 + t12 * csi1 + t13 * csi2 ) / 12.;
691 fp1[k][j][i].y -= ( t21 * csi0 + t22 * csi1 + t23 * csi2 ) / 12.;
692 fp1[k][j][i].z -= ( t31 * csi0 + t32 * csi1 + t33 * csi2 ) / 12.;
693 }
694
695 }
696 }
697 }
698 DMDAVecRestoreArray(da, user->lIAj, &iaj);
699
700
701 // j direction
702 DMDAVecGetArray(da, user->lJAj, &jaj);
703 for (k=lzs; k<lze; k++) {
704 for (j=lys-1; j<lye; j++) {
705 for (i=lxs; i<lxe; i++) {
706
707 if ((nvert[k][j ][i+1]> solid && nvert[k][j ][i+1]<innerblank)||
708 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
709 dudc = (ucat[k][j+1][i ].x + ucat[k][j][i ].x -
710 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.5;
711 dvdc = (ucat[k][j+1][i ].y + ucat[k][j][i ].y -
712 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.5;
713 dwdc = (ucat[k][j+1][i ].z + ucat[k][j][i ].z -
714 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.5;
715 }
716 else if ((nvert[k][j ][i-1]> solid && nvert[k][j ][i-1]<innerblank) ||
717 (nvert[k][j+1][i-1]> solid && nvert[k][j+1][i-1]<innerblank)) {
718 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
719 ucat[k][j+1][i ].x - ucat[k][j][i ].x) * 0.5;
720 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
721 ucat[k][j+1][i ].y - ucat[k][j][i ].y) * 0.5;
722 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
723 ucat[k][j+1][i ].z - ucat[k][j][i ].z) * 0.5;
724 }
725 else {
726 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
727 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.25;
728 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
729 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.25;
730 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
731 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.25;
732 }
733
734 dude = ucat[k][j+1][i].x - ucat[k][j][i].x;
735 dvde = ucat[k][j+1][i].y - ucat[k][j][i].y;
736 dwde = ucat[k][j+1][i].z - ucat[k][j][i].z;
737
738 if ((nvert[k+1][j ][i]> solid && nvert[k+1][j ][i]<innerblank)||
739 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
740 dudz = (ucat[k ][j+1][i].x + ucat[k ][j][i].x -
741 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.5;
742 dvdz = (ucat[k ][j+1][i].y + ucat[k ][j][i].y -
743 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.5;
744 dwdz = (ucat[k ][j+1][i].z + ucat[k ][j][i].z -
745 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.5;
746 }
747 else if ((nvert[k-1][j ][i]> solid && nvert[k-1][j ][i]<innerblank)||
748 (nvert[k-1][j+1][i]> solid && nvert[k-1][j+1][i]<innerblank)) {
749 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
750 ucat[k ][j+1][i].x - ucat[k ][j][i].x) * 0.5;
751 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
752 ucat[k ][j+1][i].y - ucat[k ][j][i].y) * 0.5;
753 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
754 ucat[k ][j+1][i].z - ucat[k ][j][i].z) * 0.5;
755 }
756 else {
757 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
758 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.25;
759 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
760 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.25;
761 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
762 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.25;
763 }
764
765 csi0 = jcsi[k][j][i].x;
766 csi1 = jcsi[k][j][i].y;
767 csi2 = jcsi[k][j][i].z;
768
769 eta0 = jeta[k][j][i].x;
770 eta1 = jeta[k][j][i].y;
771 eta2 = jeta[k][j][i].z;
772
773 zet0 = jzet[k][j][i].x;
774 zet1 = jzet[k][j][i].y;
775 zet2 = jzet[k][j][i].z;
776
777
778 g11 = csi0 * eta0 + csi1 * eta1 + csi2 * eta2;
779 g21 = eta0 * eta0 + eta1 * eta1 + eta2 * eta2;
780 g31 = zet0 * eta0 + zet1 * eta1 + zet2 * eta2;
781
782 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
783 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
784 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
785
786 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
787 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
788 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
789
790 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
791 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
792 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
793
794 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d dvdc is %.15le dvde is %.15le dvdz is %.15le \n",i,j,k,dvdc,dvde,dvdz);
795 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d dwdc is %.15le dwde is %.15le dwdz is %.15le \n",i,j,k,dwdc,dwde,dwdz);
796 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d jcsi is %.15le jeta is %.15le jzet is %.15le \n",i,j,k,jcsi[k][j][i].z,jeta[k][j][i].z,jzet[k][j][i].z);
797 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d r13 is %.15le r23 is %.15le r33 is %.15le \n",i,j,k,r13,r23,r33);
798
799
800
801 ajc = jaj[k][j][i];
802
803 double nu = 1./ren, nu_t = 0;
804
805 if( les || (rans && ti>0) ) {
806 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j+1][i]) ), 2.0) * Sabs;
807 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j+1][i]);
808 if ( (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == WALL && j==0) || (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == WALL && j==my-2) ) nu_t=0;
809
810 fp2[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * eta0 + r21 * eta1 + r31 * eta2) * ajc * (nu_t);
811 fp2[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * eta0 + r22 * eta1 + r32 * eta2) * ajc * (nu_t);
812 fp2[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * eta0 + r23 * eta1 + r33 * eta2) * ajc * (nu_t);
813 }
814 else {
815 fp2[k][j][i].x = 0;
816 fp2[k][j][i].y = 0;
817 fp2[k][j][i].z = 0;
818 }
819
820 fp2[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * eta0 + r21 * eta1 + r31 * eta2 ) * ajc * (nu);
821 fp2[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * eta0 + r22 * eta1 + r32 * eta2 ) * ajc * (nu);
822 fp2[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * eta0 + r23 * eta1 + r33 * eta2 ) * ajc * (nu);
823
824 if(clark) {
825 double dc, de, dz;
826 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
827 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
828
829 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
830 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
831 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
832 double t21 = t12;
833 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
834 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
835 double t31 = t13;
836 double t32 = t23;
837 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
838
839 fp2[k][j][i].x -= ( t11 * eta0 + t12 * eta1 + t13 * eta2 ) / 12.;
840 fp2[k][j][i].y -= ( t21 * eta0 + t22 * eta1 + t23 * eta2 ) / 12.;
841 fp2[k][j][i].z -= ( t31 * eta0 + t32 * eta1 + t33 * eta2 ) / 12.;
842 }
843 }
844 }
845 }
846
847 DMDAVecRestoreArray(da, user->lJAj, &jaj);
848 // k direction
849
850 DMDAVecGetArray(da, user->lKAj, &kaj);
851 for (k=lzs-1; k<lze; k++) {
852 for (j=lys; j<lye; j++) {
853 for (i=lxs; i<lxe; i++) {
854 if ((nvert[k ][j][i+1]> solid && nvert[k ][j][i+1]<innerblank)||
855 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
856 dudc = (ucat[k+1][j][i ].x + ucat[k][j][i ].x -
857 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.5;
858 dvdc = (ucat[k+1][j][i ].y + ucat[k][j][i ].y -
859 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.5;
860 dwdc = (ucat[k+1][j][i ].z + ucat[k][j][i ].z -
861 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.5;
862 }
863 else if ((nvert[k ][j][i-1]> solid && nvert[k ][j][i-1]<innerblank) ||
864 (nvert[k+1][j][i-1]> solid && nvert[k+1][j][i-1]<innerblank)) {
865 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
866 ucat[k+1][j][i ].x - ucat[k][j][i ].x) * 0.5;
867 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
868 ucat[k+1][j][i ].y - ucat[k][j][i ].y) * 0.5;
869 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
870 ucat[k+1][j][i ].z - ucat[k][j][i ].z) * 0.5;
871 }
872 else {
873 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
874 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.25;
875 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
876 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.25;
877 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
878 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.25;
879 }
880
881 if ((nvert[k ][j+1][i]> solid && nvert[k ][j+1][i]<innerblank)||
882 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
883 dude = (ucat[k+1][j ][i].x + ucat[k][j ][i].x -
884 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.5;
885 dvde = (ucat[k+1][j ][i].y + ucat[k][j ][i].y -
886 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.5;
887 dwde = (ucat[k+1][j ][i].z + ucat[k][j ][i].z -
888 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.5;
889 }
890 else if ((nvert[k ][j-1][i]> solid && nvert[k ][j-1][i]<innerblank) ||
891 (nvert[k+1][j-1][i]> solid && nvert[k+1][j-1][i]<innerblank)){
892 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
893 ucat[k+1][j ][i].x - ucat[k][j ][i].x) * 0.5;
894 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
895 ucat[k+1][j ][i].y - ucat[k][j ][i].y) * 0.5;
896 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
897 ucat[k+1][j ][i].z - ucat[k][j ][i].z) * 0.5;
898 }
899 else {
900 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
901 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.25;
902 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
903 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.25;
904 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
905 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.25;
906 }
907
908 dudz = ucat[k+1][j][i].x - ucat[k][j][i].x;
909 dvdz = ucat[k+1][j][i].y - ucat[k][j][i].y;
910 dwdz = ucat[k+1][j][i].z - ucat[k][j][i].z;
911
912
913 csi0 = kcsi[k][j][i].x;
914 csi1 = kcsi[k][j][i].y;
915 csi2 = kcsi[k][j][i].z;
916
917 eta0 = keta[k][j][i].x;
918 eta1 = keta[k][j][i].y;
919 eta2 = keta[k][j][i].z;
920
921 zet0 = kzet[k][j][i].x;
922 zet1 = kzet[k][j][i].y;
923 zet2 = kzet[k][j][i].z;
924
925
926 g11 = csi0 * zet0 + csi1 * zet1 + csi2 * zet2;
927 g21 = eta0 * zet0 + eta1 * zet1 + eta2 * zet2;
928 g31 = zet0 * zet0 + zet1 * zet1 + zet2 * zet2;
929
930 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
931 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
932 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
933
934 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
935 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
936 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
937
938 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
939 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
940 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
941
942 ajc = kaj[k][j][i];
943
944 double nu = 1./ren, nu_t =0;
945
946 if( les || (rans && ti>0) ) {
947 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k+1][j][i]) ), 2.0) * Sabs;
948 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k+1][j][i]);
949 if ( (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == WALL && k==0) || (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == WALL && k==mz-2) ) nu_t=0;
950
951 fp3[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu_t);
952 fp3[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu_t);
953 fp3[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu_t);
954 }
955 else {
956 fp3[k][j][i].x = 0;
957 fp3[k][j][i].y = 0;
958 fp3[k][j][i].z = 0;
959 }
960 fp3[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu);//
961 fp3[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu);//
962 fp3[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu);//
963
964 if(clark) {
965 double dc, de, dz;
966 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
967 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
968
969 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
970 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
971 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
972 double t21 = t12;
973 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
974 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
975 double t31 = t13;
976 double t32 = t23;
977 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
978
979 fp3[k][j][i].x -= ( t11 * zet0 + t12 * zet1 + t13 * zet2 ) / 12.;
980 fp3[k][j][i].y -= ( t21 * zet0 + t22 * zet1 + t23 * zet2 ) / 12.;
981 fp3[k][j][i].z -= ( t31 * zet0 + t32 * zet1 + t33 * zet2 ) / 12.;
982 }
983 }
984 }
985 }
986
987 DMDAVecRestoreArray(da, user->lKAj, &kaj);
988
989 for (k=lzs; k<lze; k++) {
990 for (j=lys; j<lye; j++) {
991 for (i=lxs; i<lxe; i++) {
992 visc[k][j][i].x =
993 (fp1[k][j][i].x - fp1[k][j][i-1].x +
994 fp2[k][j][i].x - fp2[k][j-1][i].x +
995 fp3[k][j][i].x - fp3[k-1][j][i].x);
996
997 visc[k][j][i].y =
998 (fp1[k][j][i].y - fp1[k][j][i-1].y +
999 fp2[k][j][i].y - fp2[k][j-1][i].y +
1000 fp3[k][j][i].y - fp3[k-1][j][i].y);
1001
1002 visc[k][j][i].z =
1003 (fp1[k][j][i].z - fp1[k][j][i-1].z +
1004 fp2[k][j][i].z - fp2[k][j-1][i].z +
1005 fp3[k][j][i].z - fp3[k-1][j][i].z);
1006
1007 }
1008 }
1009 }
1010/* for (k=zs; k<ze; k++) { */
1011/* for (j=ys; j<ye; j++) { */
1012/* for (i=xs; i<xe; i++) { */
1013/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp1.z is %.15le \n",i,j,k,fp1[k][j][i].z); */
1014/* if (i==0 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp1.z is %.15le \n",i,j,k,fp1[k][j][i].z); */
1015/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp2.z is %.15le \n",i,j,k,fp2[k][j][i].z); */
1016/* if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp2.z is %.15le \n",i,j,k,fp2[k][j][i].z); */
1017/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp3.z is %.15le \n",i,j,k,fp3[k][j][i].z); */
1018/* if (i==1 && j==1 && k==20) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp3.z is %.15le \n",i,j,k,fp3[k][j][i].z); */
1019
1020/* } */
1021/* } */
1022/* } */
1023 DMDAVecRestoreArray(fda, Ucont, &ucont);
1024 DMDAVecRestoreArray(fda, Ucat, &ucat);
1025 DMDAVecRestoreArray(fda, Visc, &visc);
1026
1027 DMDAVecRestoreArray(fda, Csi, &csi);
1028 DMDAVecRestoreArray(fda, Eta, &eta);
1029 DMDAVecRestoreArray(fda, Zet, &zet);
1030
1031 DMDAVecRestoreArray(fda, Fp1, &fp1);
1032 DMDAVecRestoreArray(fda, Fp2, &fp2);
1033 DMDAVecRestoreArray(fda, Fp3, &fp3);
1034
1035 DMDAVecRestoreArray(da, user->lAj, &aj);
1036
1037 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1038 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1039 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1040
1041 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1042 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1043 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1044
1045 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1046 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1047 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1048
1049 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1050
1051 if(les) {
1052 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1053 } else if (rans) {
1054
1055 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1056 }
1057
1058
1059 VecDestroy(&Fp1);
1060 VecDestroy(&Fp2);
1061 VecDestroy(&Fp3);
1062
1063
1064 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Viscous terms calculated .\n");
1065
1066 PROFILE_FUNCTION_END;
1067
1068 return(0);
1069}
1070
1071#undef __FUNCT__
1072#define __FUNCT__ "ComputeBodyForces"
1073/**
1074 * @brief Internal helper implementation: `ComputeBodyForces()`.
1075 * @details Local to this translation unit.
1076 */
1077PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
1078{
1079 PetscErrorCode ierr;
1080 PetscFunctionBeginUser;
1081
1082 // --- 1. Apply momentum source for driven channel/pipe flows ---
1083 // This function will internally check if a driven flow BC is active.
1084 ierr = ComputeDrivenChannelFlowSource(user, Rct); CHKERRQ(ierr);
1085
1086 // --- 2. (Future Extension) Apply gravitational force ---
1087 // if (user->simCtx->gravityEnabled) {
1088 // ierr = ApplyGravitationalForce(user, Rhs); CHKERRQ(ierr);
1089 // }
1090
1091 PetscFunctionReturn(0);
1092}
1093
1094#undef __FUNCT__
1095#define __FUNCT__ "ComputeRHS"
1096/**
1097 * @brief Internal helper implementation: `ComputeRHS()`.
1098 * @details Local to this translation unit.
1099 */
1100PetscErrorCode ComputeRHS(UserCtx *user, Vec Rhs)
1101{
1102 PetscErrorCode ierr;
1103 SimCtx *simCtx = user->simCtx;
1104 DM da = user->da, fda = user->fda;
1105 DMDALocalInfo info = user->info;
1106 PetscInt i,j,k;
1107 // --- Local Grid Indices and Parameters ---
1108 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1109 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1110 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1111 PetscInt lxs = (xs==0) ? xs+1 : xs;
1112 PetscInt lys = (ys==0) ? ys+1 : ys;
1113 PetscInt lzs = (zs==0) ? zs+1 : zs;
1114 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1115 PetscInt lye = (ye==my) ? ye-1 : ye;
1116 PetscInt lze = (ze==mz) ? ze-1 : ze;
1117
1118 // --- Array Pointers ---
1119 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1120 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert;
1121 Cmpnts ***rhs, ***rc, ***rct;
1122
1123 // --- Temporary Vectors ---
1124 Vec Conv, Visc, Rc, Rct;
1125
1126 PetscFunctionBeginUser;
1127 PROFILE_FUNCTION_BEGIN;
1128 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1129 simCtx->rank, user->_this);
1130
1131 // --- Get all necessary array pointers ---
1132 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1133 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1134 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1135 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1136 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1137 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1138 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1139 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1140 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1141 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1142 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1143 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1144 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1145 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1146 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1147 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1148 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1149 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1150 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1151
1152 // --- Create temporary work vectors ---
1153 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1154 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1155 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1156 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1157
1158 // ========================================================================
1159 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1160 // ========================================================================
1161
1162 // 1. Obtain Cartesian velocity from Contravariant velocity
1163 ierr = Contra2Cart(user); CHKERRQ(ierr);
1164 {
1165 const char *cell_fields[] = {"Ucat"};
1166 ierr = SynchronizePeriodicCellFields(user, 1, cell_fields); CHKERRQ(ierr);
1167 }
1168 ierr = UpdateLocalGhosts(user,"Ucat"); CHKERRQ(ierr);
1169
1170 // 2. Compute Convective term
1171 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1172 if (simCtx->moveframe || simCtx->rotateframe) {
1173 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1174 } else {
1175 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1176 }
1177
1178 // 3. Compute Viscous term
1179 if (simCtx->invicid) {
1180 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1181 } else {
1182 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1183 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1184 }
1185
1186 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1187 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1188
1189 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1190 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1191 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1192 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1193
1194 for (k = lzs; k < lze; k++) {
1195 for (j = lys; j < lye; j++) {
1196 for (i = lxs; i < lxe; i++) {
1197 rct[k][j][i].x = aj[k][j][i] *
1198 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1199 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1200 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1201 rct[k][j][i].y = aj[k][j][i] *
1202 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1203 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1204 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1205 rct[k][j][i].z = aj[k][j][i] *
1206 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1207 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1208 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1209 }
1210 }
1211 }
1212 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1213 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1214
1215 PetscBarrier(NULL);
1216
1217 // Compute and Add Body Force term if applicable.
1218 ierr = ComputeBodyForces(user,Rct); CHKERRQ(ierr);
1219 ierr = SynchronizePeriodicLocalStaggeredField(user, Rct); CHKERRQ(ierr);
1220
1221 // 6. Add Pressure Gradient Term and Finalize RHS
1222 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1223 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1224 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1225
1226 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1227
1228 for (k = lzs; k < lze; k++) {
1229 for (j = lys; j < lye; j++) {
1230 for (i = lxs; i < lxe; i++) {
1231 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1232 dpdc = p[k][j][i+1] - p[k][j][i];
1233
1234 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) {
1235 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))) {
1236 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1237 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1238 }
1239 }
1240 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) {
1241 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1242 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1243 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1244 }
1245 }
1246 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) {
1247 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1248 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1249 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1250 }
1251 }
1252 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) {
1253 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1254 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1255 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1256 }
1257 }
1258 else {
1259 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1260 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1261 }
1262
1263 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) {
1264 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))) {
1265 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1266 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1267 }
1268 }
1269 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) {
1270 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1271 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1272 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1273 }
1274 }
1275 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) {
1276 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1277 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1278 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1279 }
1280 }
1281 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) {
1282 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1283 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1284 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1285 }
1286 }
1287 else {
1288 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1289 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1290 }
1291
1292 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1293
1294
1295 rhs[k][j][i].x -=
1296 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1297 icsi[k][j][i].y * icsi[k][j][i].y +
1298 icsi[k][j][i].z * icsi[k][j][i].z)+
1299 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1300 ieta[k][j][i].y * icsi[k][j][i].y +
1301 ieta[k][j][i].z * icsi[k][j][i].z)+
1302 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1303 izet[k][j][i].y * icsi[k][j][i].y +
1304 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1305
1306 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) {
1307 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))) {
1308 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1309 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1310 }
1311 }
1312 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) {
1313 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1314 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1315 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1316 }
1317 }
1318 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) {
1319 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1320 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1321 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1322 }
1323 }
1324 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) {
1325 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1326 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1327 p[k][j][i] - p[k][j+1][i]) * 0.5;
1328 }
1329 }
1330 else {
1331 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1332 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1333 }
1334
1335 dpde = p[k][j+1][i] - p[k][j][i];
1336
1337 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) {
1338 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))) {
1339 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1340 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1341 }
1342 }
1343 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) {
1344 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1345 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1346 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1347 }
1348 }
1349 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) {
1350 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1351 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1352 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1353 }
1354 }
1355 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) {
1356 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1357 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1358 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1359 }
1360 }
1361 else {
1362 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1363 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1364 }
1365
1366 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1367
1368
1369 rhs[k][j][i].y -=
1370 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1371 jcsi[k][j][i].y * jeta[k][j][i].y +
1372 jcsi[k][j][i].z * jeta[k][j][i].z) +
1373 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1374 jeta[k][j][i].y * jeta[k][j][i].y +
1375 jeta[k][j][i].z * jeta[k][j][i].z) +
1376 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1377 jzet[k][j][i].y * jeta[k][j][i].y +
1378 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1379
1380 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) {
1381 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))) {
1382 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1383 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1384 }
1385 }
1386 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) {
1387 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1388 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1389 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1390 }
1391 }
1392 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) {
1393 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1394 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1395 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1396 }
1397 }
1398 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) {
1399 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1400 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1401 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1402 }
1403 }
1404 else {
1405 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1406 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1407 }
1408
1409 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) {
1410 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))) {
1411 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1412 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1413 }
1414 }
1415 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) {
1416 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1417 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1418 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1419 }
1420 }
1421 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) {
1422 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1423 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1424 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1425 }
1426 }
1427 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) {
1428 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1429 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1430 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1431 }
1432 }
1433 else {
1434 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1435 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1436 }
1437
1438 dpdz = (p[k+1][j][i] - p[k][j][i]);
1439
1440 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1441
1442 rhs[k][j][i].z -=
1443 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1444 kcsi[k][j][i].y * kzet[k][j][i].y +
1445 kcsi[k][j][i].z * kzet[k][j][i].z) +
1446 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1447 keta[k][j][i].y * kzet[k][j][i].y +
1448 keta[k][j][i].z * kzet[k][j][i].z) +
1449 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1450 kzet[k][j][i].y * kzet[k][j][i].y +
1451 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1452
1453 }
1454 }
1455 }
1456
1457
1458 //Mohsen March 2012//
1459
1460 // rhs.x at boundaries for periodic bc at i direction//
1461 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs==0){
1462 for (k=lzs; k<lze; k++) {
1463 for (j=lys; j<lye; j++) {
1464 i=xs;
1465 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1466
1467 dpdc = p[k][j][i+1] - p[k][j][i];
1468
1469 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) {
1470 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))) {
1471 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1472 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1473 }
1474 }
1475 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) {
1476 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1477 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1478 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1479 }
1480 }
1481 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) {
1482 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1483 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1484 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1485 }
1486 }
1487 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) {
1488 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1489 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1490 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1491 }
1492 }
1493 else {
1494 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1495 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1496 }
1497
1498 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) {
1499 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))) {
1500 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1501 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1502 }
1503 }
1504 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) {
1505 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1506 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1507 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1508 }
1509 }
1510 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) {
1511 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1512 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1513 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1514 }
1515 }
1516 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) {
1517 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1518 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1519 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1520 }
1521 }
1522 else {
1523 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1524 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1525 }
1526
1527 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1528 rhs[k][j][i].x -=
1529 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1530 icsi[k][j][i].y * icsi[k][j][i].y +
1531 icsi[k][j][i].z * icsi[k][j][i].z)+
1532 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1533 ieta[k][j][i].y * icsi[k][j][i].y +
1534 ieta[k][j][i].z * icsi[k][j][i].z)+
1535 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1536 izet[k][j][i].y * icsi[k][j][i].y +
1537 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1538 }
1539 }
1540 }
1541
1542// rhs.y at boundaries for periodic bc at j direction//
1543 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys==0){
1544 for (k=lzs; k<lze; k++) {
1545 for (i=lxs; i<lxe; i++) {
1546
1547 j=ys;
1548 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1549
1550 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) {
1551 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))) {
1552 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1553 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1554 }
1555 }
1556 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) {
1557 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1558 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1559 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1560 }
1561 }
1562 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) {
1563 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1564 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1565 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1566 }
1567 }
1568 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) {
1569 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1570 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1571 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1572 }
1573 }
1574 else {
1575 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1576 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1577 }
1578
1579 dpde = p[k][j+1][i] - p[k][j][i];
1580
1581 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) {
1582 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))) {
1583 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1584 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1585 }
1586 }
1587 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) {
1588 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1589 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1590 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1591 }
1592 }
1593 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) {
1594 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1595 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1596 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1597 }
1598 }
1599 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) {
1600 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1601 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1602 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1603 }
1604 }
1605 else {
1606 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1607 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1608 }
1609
1610 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1611
1612 rhs[k][j][i].y -=
1613 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1614 jcsi[k][j][i].y * jeta[k][j][i].y +
1615 jcsi[k][j][i].z * jeta[k][j][i].z)+
1616 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1617 jeta[k][j][i].y * jeta[k][j][i].y +
1618 jeta[k][j][i].z * jeta[k][j][i].z)+
1619 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1620 jzet[k][j][i].y * jeta[k][j][i].y +
1621 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1622
1623 }
1624 }
1625 }
1626
1627 // rhs.z at boundaries for periodic bc at k direction//
1628 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC&& zs==0){
1629 for (j=lys; j<lye; j++) {
1630 for (i=lxs; i<lxe; i++) {
1631
1632 k=zs;
1633 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1634
1635 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) {
1636 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))) {
1637 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1638 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1639 }
1640 }
1641 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) {
1642 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1643 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1644 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1645 }
1646 }
1647 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) {
1648 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1649 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1650 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1651 }
1652 }
1653 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) {
1654 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1655 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1656 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1657 }
1658 }
1659 else {
1660 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1661 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1662 }
1663
1664 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) {
1665 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))) {
1666 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1667 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1668 }
1669 }
1670 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) {
1671 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1672 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1673 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1674 }
1675 }
1676 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) {
1677 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1678 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1679 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1680 }
1681 }
1682 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) {
1683 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1684 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1685 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1686 }
1687 }
1688 else {
1689 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1690 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1691 }
1692
1693 dpdz = (p[k+1][j][i] - p[k][j][i]);
1694
1695 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1696
1697 rhs[k][j][i].z -=
1698 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1699 kcsi[k][j][i].y * kzet[k][j][i].y +
1700 kcsi[k][j][i].z * kzet[k][j][i].z)+
1701 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1702 keta[k][j][i].y * kzet[k][j][i].y +
1703 keta[k][j][i].z * kzet[k][j][i].z)+
1704 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1705 kzet[k][j][i].y * kzet[k][j][i].y +
1706 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1707
1708 }
1709 }
1710 }
1711
1712 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1713
1714 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1715 PetscInt TwoD = simCtx->TwoD;
1716
1717 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1718
1719 // 7. Final clean-up for immersed boundaries and 2D cases
1720 for (k=lzs; k<lze; k++) {
1721 for (j=lys; j<lye; j++) {
1722 for (i=lxs; i<lxe; i++) {
1723 if (TwoD==1)
1724 rhs[k][j][i].x =0.;
1725 else if (TwoD==2)
1726 rhs[k][j][i].y =0.;
1727 else if (TwoD==3)
1728 rhs[k][j][i].z =0.;
1729
1730 if (nvert[k][j][i]>0.1) {
1731 rhs[k][j][i].x = 0;
1732 rhs[k][j][i].y = 0;
1733 rhs[k][j][i].z = 0;
1734 }
1735 if (nvert[k][j][i+1]>0.1) {
1736 rhs[k][j][i].x=0;
1737 }
1738 if (nvert[k][j+1][i]>0.1) {
1739 rhs[k][j][i].y=0;
1740 }
1741 if (nvert[k+1][j][i]>0.1) {
1742 rhs[k][j][i].z=0;
1743 }
1744 }
1745 }
1746 }
1747 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1748
1749 // ========================================================================
1750
1751 // --- Restore all PETSc array pointers ---
1752 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1753 ierr = DMDAVecRestoreArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1754 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1755
1756 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1757 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1758 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1759 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1760 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1761
1762 ierr = DMDAVecRestoreArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1763 ierr = DMDAVecRestoreArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1764 ierr = DMDAVecRestoreArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1765 ierr = DMDAVecRestoreArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1766 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1767
1768 ierr = DMDAVecRestoreArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1769 ierr = DMDAVecRestoreArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1770 ierr = DMDAVecRestoreArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1771 ierr = DMDAVecRestoreArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1772 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1773
1774 ierr = DMDAVecRestoreArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1775 ierr = DMDAVecRestoreArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1776 ierr = DMDAVecRestoreArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1777 ierr = DMDAVecRestoreArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1778 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1779
1780 ierr = DMDAVecRestoreArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1781 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1782
1783 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1784 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1785
1786 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1787
1788 // --- Destroy temporary work vectors ---
1789 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1790 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1791 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1792 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1793 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1794
1795 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1796 simCtx->rank, user->_this);
1797
1798 PROFILE_FUNCTION_END;
1799 PetscFunctionReturn(0);
1800}
1801
1802#undef __FUNCT__
1803#define __FUNCT__ "ComputeEulerianDiffusivity"
1804/**
1805 * @brief Implementation of \ref ComputeEulerianDiffusivity().
1806 * @details Full API contract (arguments, ownership, side effects) is documented with
1807 * the header declaration in `include/rhs.h`.
1808 * @see ComputeEulerianDiffusivity()
1809 */
1810PetscErrorCode ComputeEulerianDiffusivity(UserCtx *user)
1811{
1812 PetscErrorCode ierr;
1813 DM da = user->da;
1814 PetscInt i, j, k, xs, ys, zs, xm, ym, zm, xe, ye, ze;
1815 PetscInt lxs, lys, lzs, lxe, lye, lze;
1816
1817 // Pointers for 3D grid access
1818 PetscReal ***diff_arr; // Output: Diffusivity field
1819 PetscReal ***nut_arr; // Input: Eddy Viscosity field (optional)
1820
1821 // Physics parameters
1822 PetscReal nu_molecular, gamma_molecular;
1823 PetscReal nu_turbulent, gamma_turbulent;
1824 PetscReal Sc, Sct;
1825 PetscBool use_turbulence_model;
1826
1827 PetscFunctionBeginUser;
1828 PROFILE_FUNCTION_BEGIN;
1829
1830 if (VerificationDiffusivityOverrideActive(user->simCtx)) {
1831 ierr = ApplyVerificationDiffusivityOverride(user); CHKERRQ(ierr);
1832 PROFILE_FUNCTION_END;
1833 PetscFunctionReturn(0);
1834 }
1835
1836 // ------------------------------------------------------------------------
1837 // 1. Parameter Setup & Safety Checks
1838 // ------------------------------------------------------------------------
1839
1840 // Determine Molecular Viscosity (nu = 1/Re)
1841 // Guard against division by zero if Re is not set or infinite (inviscid)
1842 if (user->simCtx->ren > 1.0e-12) {
1843 nu_molecular = 1.0 / user->simCtx->ren;
1844 } else {
1845 nu_molecular = 0.0;
1846 }
1847
1848 // Set Schmidt Numbers (Default to 1.0 if not provided to prevent NaN)
1849 Sc = (user->simCtx->schmidt_number > 1.0e-6) ? user->simCtx->schmidt_number : 1.0;
1850 Sct = (user->simCtx->Turbulent_schmidt_number > 1.0e-6) ? user->simCtx->Turbulent_schmidt_number : 0.7;
1851
1852 // Pre-calculate molecular component
1853 gamma_molecular = nu_molecular / Sc;
1854
1855 // Check if a turbulence model is active (LES or RANS)
1856 use_turbulence_model = (user->simCtx->les || user->simCtx->rans) ? PETSC_TRUE : PETSC_FALSE;
1857
1858 // ------------------------------------------------------------------------
1859 // 2. Data Access
1860 // ------------------------------------------------------------------------
1861
1862 // Get local grid boundaries
1863 DMDALocalInfo info;
1864 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
1865
1866 xs = info.xs; ys = info.ys; zs = info.zs;
1867 xm = info.xm; ym = info.ym; zm = info.zm;
1868 xe = xs + xm; ye = ys + ym; ze = zs + zm;
1869
1870 lxs = (xs == 0)? xs + 1 : xs;
1871 lys = (ys == 0)? ys + 1 : ys;
1872 lzs = (zs == 0)? zs + 1 : zs;
1873 lxe = (xe == info.mx)? xe - 1 : xe;
1874 lye = (ye == info.my)? ye - 1 : ye;
1875 lze = (ze == info.mz)? ze - 1 : ze;
1876
1877 // Get write access to the output Diffusivity array
1878 ierr = DMDAVecGetArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
1879
1880 // Get read access to Eddy Viscosity only if turbulence is active
1881 if (use_turbulence_model) {
1882 ierr = DMDAVecGetArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
1883 }
1884
1885 // ------------------------------------------------------------------------
1886 // 3. Calculation Loop
1887 // ------------------------------------------------------------------------
1888
1889 for (k = lzs; k < lze; k++) {
1890 for (j = lys; j < lye; j++) {
1891 for (i = lxs; i < lxe; i++) {
1892
1893 gamma_turbulent = 0.0;
1894
1895 if (use_turbulence_model) {
1896 // Fetch local eddy viscosity
1897 nu_turbulent = nut_arr[k][j][i];
1898
1899 // NUMERICAL SAFETY:
1900 // Some turbulence models (dynamic SGS) can locally produce
1901 // slightly negative viscosity. We clamp this to 0 to prevent
1902 // negative diffusivity, which crashes the Langevin sqrt().
1903 if (nu_turbulent < 0.0) {
1904 nu_turbulent = 0.0;
1905 }
1906
1907 gamma_turbulent = nu_turbulent / Sct;
1908 }
1909
1910 // Sum components
1911 diff_arr[k][j][i] = gamma_molecular + gamma_turbulent;
1912 }
1913 }
1914 }
1915
1916 // ------------------------------------------------------------------------
1917 // 4. Cleanup & Synchronization
1918 // ------------------------------------------------------------------------
1919
1920 // Restore arrays
1921 if (use_turbulence_model) {
1922 ierr = DMDAVecRestoreArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
1923 }
1924 ierr = DMDAVecRestoreArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
1925
1926 // Update Ghost Points
1927 // This is required because downstream operations (Drift Gradient Calculation
1928 // and Particle Interpolation) will need access to the halo regions of this field.
1929 const char *periodic_fields[] = {"Diffusivity"};
1930 ierr = SynchronizePeriodicCellFields(user, 1, periodic_fields); CHKERRQ(ierr);
1931 ierr = UpdateLocalGhosts(user,"Diffusivity"); CHKERRQ(ierr);
1932 PROFILE_FUNCTION_END;
1933 PetscFunctionReturn(0);
1934}
1935
1936#undef __FUNCT__
1937#define __FUNCT__ "ComputeEulerianDiffusivityGradient"
1938/**
1939 * @brief Internal helper implementation: `ComputeEulerianDiffusivityGradient()`.
1940 * @details Local to this translation unit.
1941 */
1942PetscErrorCode ComputeEulerianDiffusivityGradient(UserCtx *user)
1943{
1944 PetscErrorCode ierr;
1945 DM da = user->da, fda = user->fda;
1946 DMDALocalInfo info = user->info;
1947
1948 // 1. Determine Global Dimensions
1949 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1950
1951 // 2. Determine Local Loop Bounds (Skip Ghosts/Unused Indices)
1952 // Grid uses indices 1 to M-2 for physical cells.
1953 // Index 0 and M-1 are ghost/boundary holders.
1954
1955 // Start: If we own the global start (0), skip it and start at 1.
1956 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
1957 // End: If we own the global end (mx), stop before it (mx-1), so loop covers mx-2.
1958 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
1959
1960 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
1961 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
1962
1963 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
1964 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
1965
1966 PetscInt i, j, k;
1967
1968 // Pointers
1969 PetscReal ***diff; // Input (Scalar)
1970 Cmpnts ***grad_diff; // Output (Vector)
1971 Cmpnts ***csi, ***eta, ***zet;
1972 PetscReal ***aj;
1973
1974 // Boundary Flags (Check if Periodic)
1975 PetscBool p_x = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
1976 PetscBool p_y = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
1977 PetscBool p_z = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
1978
1979 PetscFunctionBeginUser;
1980 PROFILE_FUNCTION_BEGIN;
1981
1982 // 3. Update Ghosts for Diffusivity
1983 // Required so that Central Differences at i=2 can read i=1,
1984 // and Central Differences at Periodic Boundaries work correctly.
1985 ierr = UpdateLocalGhosts(user, "Diffusivity"); CHKERRQ(ierr);
1986
1987 // 4. Get Arrays (Read Only)
1988 ierr = DMDAVecGetArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
1989 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1990 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1991 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1992 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1993
1994 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity and Metrics arrays accessed successfully! .\n");
1995
1996 // 5. Get Output Array (Read/Write)
1997 ierr = DMDAVecGetArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
1998
1999 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity Gradient array accessed successfully! .\n");
2000
2001 // 6. Loop over Physical Domain
2002 for (k = lzs; k < lze; k++) {
2003 for (j = lys; j < lye; j++) {
2004 for (i = lxs; i < lxe; i++) {
2005
2006 PetscReal dGdCsi, dGdEta, dGdZet;
2007
2008 // ---------------------------------------------------------
2009 // I-Direction (Csi)
2010 // ---------------------------------------------------------
2011 if (!p_x && i == 1) {
2012 // Physical Start: 2nd Order Forward Difference
2013 // Stencil: [-3, 4, -1] / 2
2014 dGdCsi = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j][i+1] - diff[k][j][i+2]) * 0.5;
2015 }
2016 else if (!p_x && i == mx - 2) {
2017 // Physical End: 2nd Order Backward Difference
2018 // Stencil: [1, -4, 3] / 2
2019 dGdCsi = (3.0 * diff[k][j][i] - 4.0 * diff[k][j][i-1] + diff[k][j][i-2]) * 0.5;
2020 }
2021 else {
2022 // Interior / Periodic: Central Difference
2023 // Stencil: [-1, 0, 1] / 2
2024 dGdCsi = (diff[k][j][i+1] - diff[k][j][i-1]) * 0.5;
2025 }
2026
2027 // ---------------------------------------------------------
2028 // J-Direction (Eta)
2029 // ---------------------------------------------------------
2030 if (!p_y && j == 1) {
2031 // Physical Start: Forward
2032 dGdEta = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j+1][i] - diff[k][j+2][i]) * 0.5;
2033 }
2034 else if (!p_y && j == my - 2) {
2035 // Physical End: Backward
2036 dGdEta = (3.0 * diff[k][j][i] - 4.0 * diff[k][j-1][i] + diff[k][j-2][i]) * 0.5;
2037 }
2038 else {
2039 // Interior: Central
2040 dGdEta = (diff[k][j+1][i] - diff[k][j-1][i]) * 0.5;
2041 }
2042
2043 // ---------------------------------------------------------
2044 // K-Direction (Zet)
2045 // ---------------------------------------------------------
2046 if (!p_z && k == 1) {
2047 // Physical Start: Forward
2048 dGdZet = (-3.0 * diff[k][j][i] + 4.0 * diff[k+1][j][i] - diff[k+2][j][i]) * 0.5;
2049 }
2050 else if (!p_z && k == mz - 2) {
2051 // Physical End: Backward
2052 dGdZet = (3.0 * diff[k][j][i] - 4.0 * diff[k-1][j][i] + diff[k-2][j][i]) * 0.5;
2053 }
2054 else {
2055 // Interior: Central
2056 dGdZet = (diff[k+1][j][i] - diff[k-1][j][i]) * 0.5;
2057 }
2058
2059 // ---------------------------------------------------------
2060 // Transform to Physical Space (Cartesian Gradient)
2061 // ---------------------------------------------------------
2062 TransformScalarDerivativesToPhysical(aj[k][j][i],
2063 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2064 dGdCsi, dGdEta, dGdZet,
2065 &grad_diff[k][j][i]);
2066 }
2067 }
2068 }
2069
2070 // 7. Restore Arrays
2071 ierr = DMDAVecRestoreArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2072 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2073 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2074 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2075 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2076 ierr = DMDAVecRestoreArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2077
2078 // 8. Update Ghosts for the Result
2079 // Important: Particles near subdomain boundaries will need to interpolate
2080 // this gradient vector, so the ghosts of the vector field must be filled.
2081 ierr = UpdateLocalGhosts(user, "DiffusivityGradient"); CHKERRQ(ierr);
2082
2083 PROFILE_FUNCTION_END;
2084 PetscFunctionReturn(0);
2085}
ComputeDrivenChannelFlowSource
PetscErrorCode ComputeDrivenChannelFlowSource(UserCtx *user, Vec Rct)
Applies a momentum source term to drive flow in a periodic channel or pipe.
Definition BodyForces.c:14
PreparePeriodicQuickStencilFields
PetscErrorCode PreparePeriodicQuickStencilFields(UserCtx *user, Vec local_vector_field, Vec local_scalar_field)
Repairs the outer adjacent periodic ghosts used by QUICK cell stencils.
Definition Boundaries.c:2093
SynchronizePeriodicLocalStaggeredField
PetscErrorCode SynchronizePeriodicLocalStaggeredField(UserCtx *user, Vec local_field)
Synchronizes one local-only component-staggered periodic work field.
Definition Boundaries.c:2159
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
ComputeCellCharacteristicLengthScale
PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.
Definition Metric.c:283
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:827
LOG_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
ComputeBodyForces
PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
Internal helper implementation: ComputeBodyForces().
Definition rhs.c:1077
Viscous
PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
Implementation of Viscous().
Definition rhs.c:434
ComputeEulerianDiffusivity
PetscErrorCode ComputeEulerianDiffusivity(UserCtx *user)
Implementation of ComputeEulerianDiffusivity().
Definition rhs.c:1810
ComputeEulerianDiffusivityGradient
PetscErrorCode ComputeEulerianDiffusivityGradient(UserCtx *user)
Internal helper implementation: ComputeEulerianDiffusivityGradient().
Definition rhs.c:1942
Convection
PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
Implementation of Convection().
Definition rhs.c:13
ComputeRHS
PetscErrorCode ComputeRHS(UserCtx *user, Vec Rhs)
Internal helper implementation: ComputeRHS().
Definition rhs.c:1100
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2746
TransformScalarDerivativesToPhysical
void TransformScalarDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
Transforms scalar derivatives from computational space to physical space.
Definition setup.c:3425
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1755
LESModelType
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:518
SimCtx::clark
PetscInt clark
Definition variables.h:790
PERIODIC
@ PERIODIC
Definition variables.h:290
WALL
@ WALL
Definition variables.h:284
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:715
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:715
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:765
SimCtx::rank
PetscMPIInt rank
Definition variables.h:687
SimCtx::Turbulent_schmidt_number
PetscReal Turbulent_schmidt_number
Definition variables.h:765
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:896
UserCtx::lIEta
Vec lIEta
Definition variables.h:930
UserCtx::lIZet
Vec lIZet
Definition variables.h:930
UserCtx::lNvert
Vec lNvert
Definition variables.h:904
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:879
SimCtx::rans
PetscInt rans
Definition variables.h:789
UserCtx::lZet
Vec lZet
Definition variables.h:927
SimCtx::ren
PetscReal ren
Definition variables.h:732
UserCtx::lIAj
Vec lIAj
Definition variables.h:930
UserCtx::_this
PetscInt _this
Definition variables.h:889
UserCtx::lKEta
Vec lKEta
Definition variables.h:932
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:908
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:931
Cmpnts::x
PetscScalar x
Definition variables.h:101
SimCtx::invicid
PetscInt invicid
Definition variables.h:715
UserCtx::lKZet
Vec lKZet
Definition variables.h:932
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:935
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:935
UserCtx::lJEta
Vec lJEta
Definition variables.h:931
UserCtx::lCsi
Vec lCsi
Definition variables.h:927
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:932
SimCtx::central
PetscInt central
Definition variables.h:730
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::Diffusivity
Vec Diffusivity
Definition variables.h:907
UserCtx::lAj
Vec lAj
Definition variables.h:927
UserCtx::lICsi
Vec lICsi
Definition variables.h:930
UserCtx::info
DMDALocalInfo info
Definition variables.h:883
UserCtx::lUcat
Vec lUcat
Definition variables.h:904
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:927
SimCtx::les
PetscInt les
Definition variables.h:789
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:907
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:366
UserCtx::lJAj
Vec lJAj
Definition variables.h:931
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:715
UserCtx::lKAj
Vec lKAj
Definition variables.h:932
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:260
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:262
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:261
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:262
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:260
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:261
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
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
verification_sources.h
ApplyVerificationDiffusivityOverride
PetscErrorCode ApplyVerificationDiffusivityOverride(UserCtx *user)
Populates the Eulerian diffusivity field from a verification-only source override.
Definition verification_sources.c:31
VerificationDiffusivityOverrideActive
PetscBool VerificationDiffusivityOverrideActive(const SimCtx *simCtx)
Reports whether a verification-only diffusivity override is active.
Definition verification_sources.c:11
nu_t
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:303
Generated by doxygen 1.9.8