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 // Adding a body force here: see the contract in include/BodyForces.h.
1092 // In particular, this function runs once per RESIDUAL EVALUATION, not once
1093 // per timestep, so any force carrying state across calls (a filter, ramp,
1094 // moving average, or integral term) must gate its update on simCtx->step.
1095
1096 PetscFunctionReturn(0);
1097}
1098
1099#undef __FUNCT__
1100#define __FUNCT__ "ComputeRHS"
1101/**
1102 * @brief Internal helper implementation: `ComputeRHS()`.
1103 * @details Local to this translation unit.
1104 */
1105PetscErrorCode ComputeRHS(UserCtx *user, Vec Rhs)
1106{
1107 PetscErrorCode ierr;
1108 SimCtx *simCtx = user->simCtx;
1109 DM da = user->da, fda = user->fda;
1110 DMDALocalInfo info = user->info;
1111 PetscInt i,j,k;
1112 // --- Local Grid Indices and Parameters ---
1113 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1114 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1115 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1116 PetscInt lxs = (xs==0) ? xs+1 : xs;
1117 PetscInt lys = (ys==0) ? ys+1 : ys;
1118 PetscInt lzs = (zs==0) ? zs+1 : zs;
1119 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1120 PetscInt lye = (ye==my) ? ye-1 : ye;
1121 PetscInt lze = (ze==mz) ? ze-1 : ze;
1122
1123 // --- Array Pointers ---
1124 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1125 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert;
1126 Cmpnts ***rhs, ***rc, ***rct;
1127
1128 // --- Temporary Vectors ---
1129 Vec Conv, Visc, Rc, Rct;
1130
1131 PetscFunctionBeginUser;
1132 PROFILE_FUNCTION_BEGIN;
1133 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1134 simCtx->rank, user->_this);
1135
1136 // --- Get all necessary array pointers ---
1137 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1138 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1139 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1140 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1141 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1142 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1143 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1144 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1145 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1146 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1147 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1148 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1149 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1150 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1151 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1152 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1153 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1154 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1155 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1156
1157 // --- Create temporary work vectors ---
1158 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1159 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1160 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1161 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1162
1163 // ========================================================================
1164 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1165 // ========================================================================
1166
1167 // 1. Obtain Cartesian velocity from Contravariant velocity
1168 ierr = Contra2Cart(user); CHKERRQ(ierr);
1169 {
1170 const FieldId cell_fields[] = {FIELD_ID_UCAT};
1171 ierr = SynchronizePeriodicCellFields(user, 1, cell_fields); CHKERRQ(ierr);
1172 }
1173 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
1174
1175 // 2. Compute Convective term
1176 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1177 if (simCtx->moveframe || simCtx->rotateframe) {
1178 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1179 } else {
1180 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1181 }
1182
1183 // 3. Compute Viscous term
1184 if (simCtx->invicid) {
1185 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1186 } else {
1187 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1188 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1189 }
1190
1191 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1192 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1193
1194 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1195 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1196 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1197 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1198
1199 for (k = lzs; k < lze; k++) {
1200 for (j = lys; j < lye; j++) {
1201 for (i = lxs; i < lxe; i++) {
1202 rct[k][j][i].x = aj[k][j][i] *
1203 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1204 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1205 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1206 rct[k][j][i].y = aj[k][j][i] *
1207 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1208 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1209 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1210 rct[k][j][i].z = aj[k][j][i] *
1211 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1212 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1213 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1214 }
1215 }
1216 }
1217 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1218 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1219
1220 PetscBarrier(NULL);
1221
1222 // Compute and Add Body Force term if applicable.
1223 ierr = ComputeBodyForces(user,Rct); CHKERRQ(ierr);
1224 ierr = SynchronizePeriodicLocalStaggeredField(user, Rct); CHKERRQ(ierr);
1225
1226 // 6. Add Pressure Gradient Term and Finalize RHS
1227 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1228 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1229 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1230
1231 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1232
1233 for (k = lzs; k < lze; k++) {
1234 for (j = lys; j < lye; j++) {
1235 for (i = lxs; i < lxe; i++) {
1236 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1237 dpdc = p[k][j][i+1] - p[k][j][i];
1238
1239 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) {
1240 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))) {
1241 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1242 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1243 }
1244 }
1245 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) {
1246 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1247 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1248 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1249 }
1250 }
1251 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) {
1252 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1253 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1254 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1255 }
1256 }
1257 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) {
1258 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1259 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1260 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1261 }
1262 }
1263 else {
1264 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1265 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1266 }
1267
1268 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) {
1269 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))) {
1270 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1271 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1272 }
1273 }
1274 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) {
1275 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1276 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1277 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1278 }
1279 }
1280 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) {
1281 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1282 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1283 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1284 }
1285 }
1286 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) {
1287 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1288 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1289 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1290 }
1291 }
1292 else {
1293 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1294 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1295 }
1296
1297 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1298
1299
1300 rhs[k][j][i].x -=
1301 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1302 icsi[k][j][i].y * icsi[k][j][i].y +
1303 icsi[k][j][i].z * icsi[k][j][i].z)+
1304 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1305 ieta[k][j][i].y * icsi[k][j][i].y +
1306 ieta[k][j][i].z * icsi[k][j][i].z)+
1307 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1308 izet[k][j][i].y * icsi[k][j][i].y +
1309 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1310
1311 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) {
1312 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))) {
1313 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1314 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1315 }
1316 }
1317 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) {
1318 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1319 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1320 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1321 }
1322 }
1323 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) {
1324 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1325 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1326 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1327 }
1328 }
1329 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) {
1330 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1331 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1332 p[k][j][i] - p[k][j+1][i]) * 0.5;
1333 }
1334 }
1335 else {
1336 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1337 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1338 }
1339
1340 dpde = p[k][j+1][i] - p[k][j][i];
1341
1342 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) {
1343 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))) {
1344 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1345 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1346 }
1347 }
1348 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) {
1349 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1350 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1351 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1352 }
1353 }
1354 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) {
1355 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1356 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1357 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1358 }
1359 }
1360 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) {
1361 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1362 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1363 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1364 }
1365 }
1366 else {
1367 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1368 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1369 }
1370
1371 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1372
1373
1374 rhs[k][j][i].y -=
1375 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1376 jcsi[k][j][i].y * jeta[k][j][i].y +
1377 jcsi[k][j][i].z * jeta[k][j][i].z) +
1378 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1379 jeta[k][j][i].y * jeta[k][j][i].y +
1380 jeta[k][j][i].z * jeta[k][j][i].z) +
1381 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1382 jzet[k][j][i].y * jeta[k][j][i].y +
1383 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1384
1385 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) {
1386 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))) {
1387 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1388 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1389 }
1390 }
1391 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) {
1392 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1393 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1394 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1395 }
1396 }
1397 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) {
1398 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1399 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1400 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1401 }
1402 }
1403 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) {
1404 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1405 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1406 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1407 }
1408 }
1409 else {
1410 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1411 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1412 }
1413
1414 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) {
1415 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))) {
1416 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1417 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1418 }
1419 }
1420 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) {
1421 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1422 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1423 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1424 }
1425 }
1426 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) {
1427 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1428 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1429 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1430 }
1431 }
1432 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) {
1433 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1434 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1435 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1436 }
1437 }
1438 else {
1439 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1440 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1441 }
1442
1443 dpdz = (p[k+1][j][i] - p[k][j][i]);
1444
1445 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1446
1447 rhs[k][j][i].z -=
1448 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1449 kcsi[k][j][i].y * kzet[k][j][i].y +
1450 kcsi[k][j][i].z * kzet[k][j][i].z) +
1451 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1452 keta[k][j][i].y * kzet[k][j][i].y +
1453 keta[k][j][i].z * kzet[k][j][i].z) +
1454 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1455 kzet[k][j][i].y * kzet[k][j][i].y +
1456 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1457
1458 }
1459 }
1460 }
1461
1462
1463 //Mohsen March 2012//
1464
1465 // rhs.x at boundaries for periodic bc at i direction//
1466 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs==0){
1467 for (k=lzs; k<lze; k++) {
1468 for (j=lys; j<lye; j++) {
1469 i=xs;
1470 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1471
1472 dpdc = p[k][j][i+1] - p[k][j][i];
1473
1474 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) {
1475 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))) {
1476 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1477 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1478 }
1479 }
1480 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) {
1481 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1482 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1483 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1484 }
1485 }
1486 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) {
1487 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1488 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1489 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1490 }
1491 }
1492 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) {
1493 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1494 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1495 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1496 }
1497 }
1498 else {
1499 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1500 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1501 }
1502
1503 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) {
1504 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))) {
1505 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1506 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1507 }
1508 }
1509 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) {
1510 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1511 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1512 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1513 }
1514 }
1515 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) {
1516 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1517 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1518 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1519 }
1520 }
1521 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) {
1522 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1523 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1524 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1525 }
1526 }
1527 else {
1528 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1529 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1530 }
1531
1532 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1533 rhs[k][j][i].x -=
1534 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1535 icsi[k][j][i].y * icsi[k][j][i].y +
1536 icsi[k][j][i].z * icsi[k][j][i].z)+
1537 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1538 ieta[k][j][i].y * icsi[k][j][i].y +
1539 ieta[k][j][i].z * icsi[k][j][i].z)+
1540 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1541 izet[k][j][i].y * icsi[k][j][i].y +
1542 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1543 }
1544 }
1545 }
1546
1547// rhs.y at boundaries for periodic bc at j direction//
1548 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys==0){
1549 for (k=lzs; k<lze; k++) {
1550 for (i=lxs; i<lxe; i++) {
1551
1552 j=ys;
1553 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1554
1555 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) {
1556 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))) {
1557 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1558 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1559 }
1560 }
1561 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) {
1562 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1563 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1564 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1565 }
1566 }
1567 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) {
1568 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1569 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1570 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1571 }
1572 }
1573 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) {
1574 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1575 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1576 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1577 }
1578 }
1579 else {
1580 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1581 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1582 }
1583
1584 dpde = p[k][j+1][i] - p[k][j][i];
1585
1586 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) {
1587 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))) {
1588 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1589 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1590 }
1591 }
1592 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) {
1593 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1594 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1595 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1596 }
1597 }
1598 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) {
1599 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1600 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1601 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1602 }
1603 }
1604 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) {
1605 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1606 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1607 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1608 }
1609 }
1610 else {
1611 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1612 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1613 }
1614
1615 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1616
1617 rhs[k][j][i].y -=
1618 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1619 jcsi[k][j][i].y * jeta[k][j][i].y +
1620 jcsi[k][j][i].z * jeta[k][j][i].z)+
1621 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1622 jeta[k][j][i].y * jeta[k][j][i].y +
1623 jeta[k][j][i].z * jeta[k][j][i].z)+
1624 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1625 jzet[k][j][i].y * jeta[k][j][i].y +
1626 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1627
1628 }
1629 }
1630 }
1631
1632 // rhs.z at boundaries for periodic bc at k direction//
1633 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC&& zs==0){
1634 for (j=lys; j<lye; j++) {
1635 for (i=lxs; i<lxe; i++) {
1636
1637 k=zs;
1638 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1639
1640 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) {
1641 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))) {
1642 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1643 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1644 }
1645 }
1646 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) {
1647 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1648 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1649 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1650 }
1651 }
1652 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) {
1653 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1654 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1655 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1656 }
1657 }
1658 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) {
1659 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1660 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1661 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1662 }
1663 }
1664 else {
1665 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1666 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1667 }
1668
1669 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) {
1670 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))) {
1671 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1672 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1673 }
1674 }
1675 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) {
1676 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1677 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1678 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1679 }
1680 }
1681 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) {
1682 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1683 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1684 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1685 }
1686 }
1687 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) {
1688 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1689 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1690 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1691 }
1692 }
1693 else {
1694 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1695 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1696 }
1697
1698 dpdz = (p[k+1][j][i] - p[k][j][i]);
1699
1700 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1701
1702 rhs[k][j][i].z -=
1703 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1704 kcsi[k][j][i].y * kzet[k][j][i].y +
1705 kcsi[k][j][i].z * kzet[k][j][i].z)+
1706 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1707 keta[k][j][i].y * kzet[k][j][i].y +
1708 keta[k][j][i].z * kzet[k][j][i].z)+
1709 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1710 kzet[k][j][i].y * kzet[k][j][i].y +
1711 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1712
1713 }
1714 }
1715 }
1716
1717 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1718
1719 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1720 PetscInt TwoD = simCtx->TwoD;
1721
1722 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1723
1724 // 7. Final clean-up for immersed boundaries and 2D cases
1725 for (k=lzs; k<lze; k++) {
1726 for (j=lys; j<lye; j++) {
1727 for (i=lxs; i<lxe; i++) {
1728 if (TwoD==1)
1729 rhs[k][j][i].x =0.;
1730 else if (TwoD==2)
1731 rhs[k][j][i].y =0.;
1732 else if (TwoD==3)
1733 rhs[k][j][i].z =0.;
1734
1735 if (nvert[k][j][i]>0.1) {
1736 rhs[k][j][i].x = 0;
1737 rhs[k][j][i].y = 0;
1738 rhs[k][j][i].z = 0;
1739 }
1740 if (nvert[k][j][i+1]>0.1) {
1741 rhs[k][j][i].x=0;
1742 }
1743 if (nvert[k][j+1][i]>0.1) {
1744 rhs[k][j][i].y=0;
1745 }
1746 if (nvert[k+1][j][i]>0.1) {
1747 rhs[k][j][i].z=0;
1748 }
1749 }
1750 }
1751 }
1752 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1753
1754 // ========================================================================
1755
1756 // --- Restore all PETSc array pointers ---
1757 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1758 ierr = DMDAVecRestoreArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1759 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1760
1761 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1762 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1763 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1764 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1765 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1766
1767 ierr = DMDAVecRestoreArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1768 ierr = DMDAVecRestoreArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1769 ierr = DMDAVecRestoreArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1770 ierr = DMDAVecRestoreArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1771 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1772
1773 ierr = DMDAVecRestoreArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1774 ierr = DMDAVecRestoreArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1775 ierr = DMDAVecRestoreArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1776 ierr = DMDAVecRestoreArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1777 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1778
1779 ierr = DMDAVecRestoreArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1780 ierr = DMDAVecRestoreArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1781 ierr = DMDAVecRestoreArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1782 ierr = DMDAVecRestoreArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1783 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1784
1785 ierr = DMDAVecRestoreArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1786 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1787
1788 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1789 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1790
1791 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1792
1793 // --- Destroy temporary work vectors ---
1794 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1795 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1796 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1797 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1798 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1799
1800 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1801 simCtx->rank, user->_this);
1802
1803 PROFILE_FUNCTION_END;
1804 PetscFunctionReturn(0);
1805}
1806
1807#undef __FUNCT__
1808#define __FUNCT__ "ComputeEulerianDiffusivity"
1809/**
1810 * @brief Implementation of \ref ComputeEulerianDiffusivity().
1811 * @details Full API contract (arguments, ownership, side effects) is documented with
1812 * the header declaration in `include/rhs.h`.
1813 * @see ComputeEulerianDiffusivity()
1814 */
1815PetscErrorCode ComputeEulerianDiffusivity(UserCtx *user)
1816{
1817 PetscErrorCode ierr;
1818 DM da = user->da;
1819 PetscInt i, j, k, xs, ys, zs, xm, ym, zm, xe, ye, ze;
1820 PetscInt lxs, lys, lzs, lxe, lye, lze;
1821
1822 // Pointers for 3D grid access
1823 PetscReal ***diff_arr; // Output: Diffusivity field
1824 PetscReal ***nut_arr; // Input: Eddy Viscosity field (optional)
1825
1826 // Physics parameters
1827 PetscReal nu_molecular, gamma_molecular;
1828 PetscReal nu_turbulent, gamma_turbulent;
1829 PetscReal Sc, Sct;
1830 PetscBool use_turbulence_model;
1831
1832 PetscFunctionBeginUser;
1833 PROFILE_FUNCTION_BEGIN;
1834
1835 if (VerificationDiffusivityOverrideActive(user->simCtx)) {
1836 ierr = ApplyVerificationDiffusivityOverride(user); CHKERRQ(ierr);
1837 PROFILE_FUNCTION_END;
1838 PetscFunctionReturn(0);
1839 }
1840
1841 // ------------------------------------------------------------------------
1842 // 1. Parameter Setup & Safety Checks
1843 // ------------------------------------------------------------------------
1844
1845 // Determine Molecular Viscosity (nu = 1/Re)
1846 // Guard against division by zero if Re is not set or infinite (inviscid)
1847 if (user->simCtx->ren > 1.0e-12) {
1848 nu_molecular = 1.0 / user->simCtx->ren;
1849 } else {
1850 nu_molecular = 0.0;
1851 }
1852
1853 // Set Schmidt Numbers (Default to 1.0 if not provided to prevent NaN)
1854 Sc = (user->simCtx->schmidt_number > 1.0e-6) ? user->simCtx->schmidt_number : 1.0;
1855 Sct = (user->simCtx->Turbulent_schmidt_number > 1.0e-6) ? user->simCtx->Turbulent_schmidt_number : 0.7;
1856
1857 // Pre-calculate molecular component
1858 gamma_molecular = nu_molecular / Sc;
1859
1860 // Check if a turbulence model is active (LES or RANS)
1861 use_turbulence_model = (user->simCtx->les || user->simCtx->rans) ? PETSC_TRUE : PETSC_FALSE;
1862
1863 // ------------------------------------------------------------------------
1864 // 2. Data Access
1865 // ------------------------------------------------------------------------
1866
1867 // Get local grid boundaries
1868 DMDALocalInfo info;
1869 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
1870
1871 xs = info.xs; ys = info.ys; zs = info.zs;
1872 xm = info.xm; ym = info.ym; zm = info.zm;
1873 xe = xs + xm; ye = ys + ym; ze = zs + zm;
1874
1875 lxs = (xs == 0)? xs + 1 : xs;
1876 lys = (ys == 0)? ys + 1 : ys;
1877 lzs = (zs == 0)? zs + 1 : zs;
1878 lxe = (xe == info.mx)? xe - 1 : xe;
1879 lye = (ye == info.my)? ye - 1 : ye;
1880 lze = (ze == info.mz)? ze - 1 : ze;
1881
1882 // Get write access to the output Diffusivity array
1883 ierr = DMDAVecGetArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
1884
1885 // Get read access to Eddy Viscosity only if turbulence is active
1886 if (use_turbulence_model) {
1887 ierr = DMDAVecGetArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
1888 }
1889
1890 // ------------------------------------------------------------------------
1891 // 3. Calculation Loop
1892 // ------------------------------------------------------------------------
1893
1894 for (k = lzs; k < lze; k++) {
1895 for (j = lys; j < lye; j++) {
1896 for (i = lxs; i < lxe; i++) {
1897
1898 gamma_turbulent = 0.0;
1899
1900 if (use_turbulence_model) {
1901 // Fetch local eddy viscosity
1902 nu_turbulent = nut_arr[k][j][i];
1903
1904 // NUMERICAL SAFETY:
1905 // Some turbulence models (dynamic SGS) can locally produce
1906 // slightly negative viscosity. We clamp this to 0 to prevent
1907 // negative diffusivity, which crashes the Langevin sqrt().
1908 if (nu_turbulent < 0.0) {
1909 nu_turbulent = 0.0;
1910 }
1911
1912 gamma_turbulent = nu_turbulent / Sct;
1913 }
1914
1915 // Sum components
1916 diff_arr[k][j][i] = gamma_molecular + gamma_turbulent;
1917 }
1918 }
1919 }
1920
1921 // ------------------------------------------------------------------------
1922 // 4. Cleanup & Synchronization
1923 // ------------------------------------------------------------------------
1924
1925 // Restore arrays
1926 if (use_turbulence_model) {
1927 ierr = DMDAVecRestoreArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
1928 }
1929 ierr = DMDAVecRestoreArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
1930
1931 // Update Ghost Points
1932 // This is required because downstream operations (Drift Gradient Calculation
1933 // and Particle Interpolation) will need access to the halo regions of this field.
1934 const FieldId periodic_fields[] = {FIELD_ID_DIFFUSIVITY};
1935 ierr = SynchronizePeriodicCellFields(user, 1, periodic_fields); CHKERRQ(ierr);
1936 ierr = UpdateLocalGhosts(user, FIELD_ID_DIFFUSIVITY); CHKERRQ(ierr);
1937 PROFILE_FUNCTION_END;
1938 PetscFunctionReturn(0);
1939}
1940
1941#undef __FUNCT__
1942#define __FUNCT__ "ComputeEulerianDiffusivityGradient"
1943/**
1944 * @brief Internal helper implementation: `ComputeEulerianDiffusivityGradient()`.
1945 * @details Local to this translation unit.
1946 */
1947PetscErrorCode ComputeEulerianDiffusivityGradient(UserCtx *user)
1948{
1949 PetscErrorCode ierr;
1950 DM da = user->da, fda = user->fda;
1951 DMDALocalInfo info = user->info;
1952
1953 // 1. Determine Global Dimensions
1954 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1955
1956 // 2. Determine Local Loop Bounds (Skip Ghosts/Unused Indices)
1957 // Grid uses indices 1 to M-2 for physical cells.
1958 // Index 0 and M-1 are ghost/boundary holders.
1959
1960 // Start: If we own the global start (0), skip it and start at 1.
1961 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
1962 // End: If we own the global end (mx), stop before it (mx-1), so loop covers mx-2.
1963 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
1964
1965 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
1966 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
1967
1968 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
1969 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
1970
1971 PetscInt i, j, k;
1972
1973 // Pointers
1974 PetscReal ***diff; // Input (Scalar)
1975 Cmpnts ***grad_diff; // Output (Vector)
1976 Cmpnts ***csi, ***eta, ***zet;
1977 PetscReal ***aj;
1978
1979 // Boundary Flags (Check if Periodic)
1980 PetscBool p_x = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
1981 PetscBool p_y = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
1982 PetscBool p_z = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
1983
1984 PetscFunctionBeginUser;
1985 PROFILE_FUNCTION_BEGIN;
1986
1987 // 3. Update Ghosts for Diffusivity
1988 // Required so that Central Differences at i=2 can read i=1,
1989 // and Central Differences at Periodic Boundaries work correctly.
1990 ierr = UpdateLocalGhosts(user, FIELD_ID_DIFFUSIVITY); CHKERRQ(ierr);
1991
1992 // 4. Get Arrays (Read Only)
1993 ierr = DMDAVecGetArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
1994 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1995 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1996 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1997 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1998
1999 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity and Metrics arrays accessed successfully! .\n");
2000
2001 // 5. Get Output Array (Read/Write)
2002 ierr = DMDAVecGetArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2003
2004 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity Gradient array accessed successfully! .\n");
2005
2006 // 6. Loop over Physical Domain
2007 for (k = lzs; k < lze; k++) {
2008 for (j = lys; j < lye; j++) {
2009 for (i = lxs; i < lxe; i++) {
2010
2011 PetscReal dGdCsi, dGdEta, dGdZet;
2012
2013 // ---------------------------------------------------------
2014 // I-Direction (Csi)
2015 // ---------------------------------------------------------
2016 if (!p_x && i == 1) {
2017 // Physical Start: 2nd Order Forward Difference
2018 // Stencil: [-3, 4, -1] / 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 if (!p_x && i == mx - 2) {
2022 // Physical End: 2nd Order Backward Difference
2023 // Stencil: [1, -4, 3] / 2
2024 dGdCsi = (3.0 * diff[k][j][i] - 4.0 * diff[k][j][i-1] + diff[k][j][i-2]) * 0.5;
2025 }
2026 else {
2027 // Interior / Periodic: Central Difference
2028 // Stencil: [-1, 0, 1] / 2
2029 dGdCsi = (diff[k][j][i+1] - diff[k][j][i-1]) * 0.5;
2030 }
2031
2032 // ---------------------------------------------------------
2033 // J-Direction (Eta)
2034 // ---------------------------------------------------------
2035 if (!p_y && j == 1) {
2036 // Physical Start: Forward
2037 dGdEta = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j+1][i] - diff[k][j+2][i]) * 0.5;
2038 }
2039 else if (!p_y && j == my - 2) {
2040 // Physical End: Backward
2041 dGdEta = (3.0 * diff[k][j][i] - 4.0 * diff[k][j-1][i] + diff[k][j-2][i]) * 0.5;
2042 }
2043 else {
2044 // Interior: Central
2045 dGdEta = (diff[k][j+1][i] - diff[k][j-1][i]) * 0.5;
2046 }
2047
2048 // ---------------------------------------------------------
2049 // K-Direction (Zet)
2050 // ---------------------------------------------------------
2051 if (!p_z && k == 1) {
2052 // Physical Start: Forward
2053 dGdZet = (-3.0 * diff[k][j][i] + 4.0 * diff[k+1][j][i] - diff[k+2][j][i]) * 0.5;
2054 }
2055 else if (!p_z && k == mz - 2) {
2056 // Physical End: Backward
2057 dGdZet = (3.0 * diff[k][j][i] - 4.0 * diff[k-1][j][i] + diff[k-2][j][i]) * 0.5;
2058 }
2059 else {
2060 // Interior: Central
2061 dGdZet = (diff[k+1][j][i] - diff[k-1][j][i]) * 0.5;
2062 }
2063
2064 // ---------------------------------------------------------
2065 // Transform to Physical Space (Cartesian Gradient)
2066 // ---------------------------------------------------------
2067 TransformScalarDerivativesToPhysical(aj[k][j][i],
2068 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2069 dGdCsi, dGdEta, dGdZet,
2070 &grad_diff[k][j][i]);
2071 }
2072 }
2073 }
2074
2075 // 7. Restore Arrays
2076 ierr = DMDAVecRestoreArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2077 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2078 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2079 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2080 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2081 ierr = DMDAVecRestoreArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2082
2083 // 8. Update Ghosts for the Result
2084 // Important: Particles near subdomain boundaries will need to interpolate
2085 // this gradient vector, so the ghosts of the vector field must be filled.
2086 ierr = UpdateLocalGhosts(user, FIELD_ID_DIFFUSIVITY_GRADIENT); CHKERRQ(ierr);
2087
2088 PROFILE_FUNCTION_END;
2089 PetscFunctionReturn(0);
2090}
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:2034
SynchronizePeriodicLocalStaggeredField
PetscErrorCode SynchronizePeriodicLocalStaggeredField(UserCtx *user, Vec local_field)
Synchronizes one local-only component-staggered periodic work field.
Definition Boundaries.c:2100
SynchronizePeriodicCellFields
PetscErrorCode SynchronizePeriodicCellFields(UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
Synchronizes periodic endpoint cells for a list of cell-centered fields.
Definition Boundaries.c:1621
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:282
FieldId
FieldId
Compile-time identity for a catalogued Eulerian field.
Definition field_catalog.h:23
FIELD_ID_UCAT
@ FIELD_ID_UCAT
Definition field_catalog.h:26
FIELD_ID_DIFFUSIVITY_GRADIENT
@ FIELD_ID_DIFFUSIVITY_GRADIENT
Definition field_catalog.h:34
FIELD_ID_DIFFUSIVITY
@ FIELD_ID_DIFFUSIVITY
Definition field_catalog.h:33
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:859
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:850
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:1815
ComputeEulerianDiffusivityGradient
PetscErrorCode ComputeEulerianDiffusivityGradient(UserCtx *user)
Internal helper implementation: ComputeEulerianDiffusivityGradient().
Definition rhs.c:1947
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:1105
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2649
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:3324
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1838
LESModelType
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:520
SimCtx::clark
PetscInt clark
Definition variables.h:822
PERIODIC
@ PERIODIC
Definition variables.h:292
WALL
@ WALL
Definition variables.h:286
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:727
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:727
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:787
SimCtx::rank
PetscMPIInt rank
Definition variables.h:698
SimCtx::Turbulent_schmidt_number
PetscReal Turbulent_schmidt_number
Definition variables.h:787
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:931
UserCtx::lIEta
Vec lIEta
Definition variables.h:977
UserCtx::lIZet
Vec lIZet
Definition variables.h:977
UserCtx::lNvert
Vec lNvert
Definition variables.h:939
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
SimCtx::rans
PetscInt rans
Definition variables.h:821
UserCtx::lZet
Vec lZet
Definition variables.h:974
SimCtx::ren
PetscReal ren
Definition variables.h:744
UserCtx::lIAj
Vec lIAj
Definition variables.h:977
UserCtx::_this
PetscInt _this
Definition variables.h:924
UserCtx::lKEta
Vec lKEta
Definition variables.h:979
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:943
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:978
Cmpnts::x
PetscScalar x
Definition variables.h:103
SimCtx::invicid
PetscInt invicid
Definition variables.h:727
UserCtx::lKZet
Vec lKZet
Definition variables.h:979
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:982
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:982
UserCtx::lJEta
Vec lJEta
Definition variables.h:978
UserCtx::lCsi
Vec lCsi
Definition variables.h:974
Cmpnts::z
PetscScalar z
Definition variables.h:103
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:979
SimCtx::central
PetscInt central
Definition variables.h:742
UserCtx::lJZet
Vec lJZet
Definition variables.h:978
UserCtx::lUcont
Vec lUcont
Definition variables.h:939
SimCtx::step
PetscInt step
Definition variables.h:703
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:942
UserCtx::lAj
Vec lAj
Definition variables.h:974
UserCtx::lICsi
Vec lICsi
Definition variables.h:977
UserCtx::info
DMDALocalInfo info
Definition variables.h:918
UserCtx::lUcat
Vec lUcat
Definition variables.h:939
Cmpnts::y
PetscScalar y
Definition variables.h:103
UserCtx::lEta
Vec lEta
Definition variables.h:974
SimCtx::les
PetscInt les
Definition variables.h:821
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:942
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:368
UserCtx::lJAj
Vec lJAj
Definition variables.h:978
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:727
UserCtx::lKAj
Vec lKAj
Definition variables.h:979
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:262
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:264
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:263
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:264
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:262
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:263
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102
SimCtx
The master context for the entire simulation.
Definition variables.h:695
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:906
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