Loading...
Searching...
No Matches
#include "rhs.h"
Include dependency graph for rhs.c:
Go to the source code of this file.
Macros | |
| #define | __FUNCT__ "Convection" |
| #define | __FUNCT__ "Viscous" |
| #define | __FUNCT__ "ComputeBodyForces" |
| #define | __FUNCT__ "ComputeRHS" |
| #define | __FUNCT__ "ComputeEulerianDiffusivity" |
| #define | __FUNCT__ "ComputeEulerianDiffusivityGradient" |
Functions | |
| PetscErrorCode | Convection (UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv) |
| PetscErrorCode | Viscous (UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc) |
| PetscErrorCode | ComputeBodyForces (UserCtx *user, Vec Rct) |
| General dispatcher for applying all active body forces (momentum sources). | |
| PetscErrorCode | ComputeRHS (UserCtx *user, Vec Rhs) |
| Computes the Right-Hand-Side (RHS) of the momentum equations. | |
| PetscErrorCode | ComputeEulerianDiffusivity (UserCtx *user) |
| Computes the effective diffusivity scalar field ($\Gamma_{eff}$) on the Eulerian grid. | |
| PetscErrorCode | ComputeEulerianDiffusivityGradient (UserCtx *user) |
| Computes the gradient of the scalar Diffusivity field (Drift Vector). | |
Macro Definition Documentation
◆ __FUNCT__ [1/6]
◆ __FUNCT__ [2/6]
◆ __FUNCT__ [3/6]
◆ __FUNCT__ [4/6]
◆ __FUNCT__ [5/6]
◆ __FUNCT__ [6/6]
Function Documentation
◆ Convection()
| PetscErrorCode Convection | ( | UserCtx * | user, |
| Vec | Ucont, | ||
| Vec | Ucat, | ||
| Vec | Conv | ||
| ) |
Definition at line 5 of file rhs.c.
6{
7
8 // --- CONTEXT ACQUISITION BLOCK ---
9 // Get the master simulation context from the UserCtx.
11
12 // Create local variables to mirror the legacy globals for minimal code changes.
15 // --- END CONTEXT ACQUISITION BLOCK ---
16
17 Cmpnts ***ucont, ***ucat;
19 DMDALocalInfo info;
20 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
21 PetscInt mx, my, mz; // Dimensions in three directions
22 PetscInt i, j, k;
23 Vec Fp1, Fp2, Fp3;
24 Cmpnts ***fp1, ***fp2, ***fp3;
25 Cmpnts ***conv;
26
27 PetscReal ucon, up, um;
28 PetscReal coef = 0.125, innerblank=7.;
29
30 PetscInt lxs, lxe, lys, lye, lzs, lze, gxs, gxe, gys, gye, gzs,gze;
31
32 PetscReal ***nvert,***aj;
33
35
36 DMDAGetLocalInfo(da, &info);
37 mx = info.mx; my = info.my; mz = info.mz;
38 xs = info.xs; xe = xs + info.xm;
39 ys = info.ys; ye = ys + info.ym;
40 zs = info.zs; ze = zs + info.zm;
41 gxs = info.gxs; gxe = gxs + info.gxm;
42 gys = info.gys; gye = gys + info.gym;
43 gzs = info.gzs; gze = gzs + info.gzm;
44
45 DMDAVecGetArray(fda, Ucont, &ucont);
46 DMDAVecGetArray(fda, Ucat, &ucat);
47 DMDAVecGetArray(fda, Conv, &conv);
48 DMDAVecGetArray(da, user->lAj, &aj);
49
50 VecDuplicate(Ucont, &Fp1);
51 VecDuplicate(Ucont, &Fp2);
52 VecDuplicate(Ucont, &Fp3);
53
54 DMDAVecGetArray(fda, Fp1, &fp1);
55 DMDAVecGetArray(fda, Fp2, &fp2);
56 DMDAVecGetArray(fda, Fp3, &fp3);
57
58 DMDAVecGetArray(da, user->lNvert, &nvert);
59
60
61 /* We have two different sets of node: 1. grid node, the physical points
62 where grid lines intercross; 2. storage node, where we store variables.
63 All node without explicitly specified as "grid node" refers to
64 storage node.
65
66 The integer node is defined at cell center while half node refers to
67 the actual grid node. (The reason to choose this arrangement is we need
68 ghost node, which is half node away from boundaries, to specify boundary
69 conditions. By using this storage arrangement, the actual storage need
70 is (IM+1) * (JM + 1) * (KM+1) where IM, JM, & KM refer to the number of
71 grid nodes along i, j, k directions.)
72
73 DA, the data structure used to define the storage of 3D arrays, is defined
74 as mx * my * mz. mx = IM+1, my = JM+1, mz = KM+1.
75
76 Staggered grid arrangement is used in this solver.
77 Pressure is stored at interger node (hence the cell center) and volume
78 fluxes defined on the center of each surface of a given control volume
79 is stored on the cloest upper integer node. */
80
81 /* First we calculate the flux on cell surfaces. Stored on the upper integer
82 node. For example, along i direction, the flux are stored at node 0:mx-2*/
83
84 lxs = xs; lxe = xe;
85 lys = ys; lye = ye;
86 lzs = zs; lze = ze;
87
88 if (xs==0) lxs = xs+1;
89 if (ys==0) lys = ys+1;
90 if (zs==0) lzs = zs+1;
91
92 if (xe==mx) lxe=xe-1;
93 if (ye==my) lye=ye-1;
94 if (ze==mz) lze=ze-1;
95
96 //Mohsen Sep 2012//
97/* First update the computational ghost points velocity for periodic boundary conditions
98 just for this subroutine because of Quick scheme for velocity deravatives */
99 /* for (k=zs; k<ze; k++) { */
100/* for (j=ys; j<ye; j++) { */
101/* for (i=xs; i<xe; i++) { */
102/* if (i==1 && (j==1) && (k==0 || k==1)) */
103/* 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 ); */
104/* } */
105/* } */
106/* } */
107 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
108 if (xs==0){
109 i=xs-1;
110 for (k=gzs; k<gze; k++) {
111 for (j=gys; j<gye; j++) {
115 nvert[k][j][i]=nvert[k][j][i-2];
116 }
117 }
118 }
119 if (xe==mx){
120 i=mx;
121 for (k=gzs; k<gze; k++) {
122 for (j=gys; j<gye; j++) {
126 nvert[k][j][i]=nvert[k][j][i+2];
127 }
128 }
129 }
130 }
131 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
132 if (ys==0){
133 j=ys-1;
134 for (k=gzs; k<gze; k++) {
135 for (i=gxs; i<gxe; i++) {
139 nvert[k][j][i]=nvert[k][j-2][i];
140 }
141 }
142 }
143 if (ye==my){
144 j=my;
145 for (k=gzs; k<gze; k++) {
146 for (i=gxs; i<gxe; i++) {
150 nvert[k][j][i]=nvert[k][j+2][i];
151 }
152 }
153 }
154 }
155 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
156 if (zs==0){
157 k=zs-1;
158 for (j=gys; j<gye; j++) {
159 for (i=gxs; i<gxe; i++) {
163 nvert[k][j][i]=nvert[k-2][j][i];
164 }
165 }
166 }
167 if (ze==mz){
168 k=mz;
169 for (j=gys; j<gye; j++) {
170 for (i=gxs; i<gxe; i++) {
174 nvert[k][j][i]=nvert[k+2][j][i];
175 }
176 }
177 }
178 }
179
180 VecSet(Conv, 0.0);
181
182 /* Calculating the convective terms on cell centers.
183 First calcualte the contribution from i direction
184 The flux is evaluated by QUICK scheme */
185
186 for (k=lzs; k<lze; k++){
187 for (j=lys; j<lye; j++){
188 for (i=lxs-1; i<lxe; i++){
189
190
191 ucon = ucont[k][j][i].x * 0.5;
192
193 up = ucon + fabs(ucon);
194 um = ucon - fabs(ucon);
195
196 if (i>0 && i<mx-2 &&
197 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
198 (nvert[k][j][i-1] < 0.1 || nvert[k][j][i-1]>innerblank)) { // interial nodes
199 if ((les || central)) {
203
204 } else {
205 fp1[k][j][i].x =
206 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) +
207 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);
208 fp1[k][j][i].y =
209 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) +
210 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);
211 fp1[k][j][i].z =
212 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) +
213 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);
214 }
215 }
216 else if ((les || central) && (i==0 || i==mx-2) &&
217 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
218 (nvert[k][j][i ] < 0.1 || nvert[k][j][i ]>innerblank))
219 {
223 }
224 else if (i==0 ||(nvert[k][j][i-1] > 0.1) ) {
225 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
226 fp1[k][j][i].x =
227 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) +
228 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);
229 fp1[k][j][i].y =
230 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) +
231 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);
232 fp1[k][j][i].z =
233 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) +
234 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);
235 }else{
236 fp1[k][j][i].x =
237 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) +
238 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);
239 fp1[k][j][i].y =
240 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) +
241 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);
242 fp1[k][j][i].z =
243 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) +
244 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);
245 }
246 }
247 else if (i==mx-2 ||(nvert[k][j][i+1]) > 0.1) {
248 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
249 fp1[k][j][i].x =
250 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) +
251 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);
252 fp1[k][j][i].y =
253 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) +
254 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);
255 fp1[k][j][i].z =
256 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) +
257 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);
258 }else{
259 fp1[k][j][i].x =
260 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) +
261 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);
262 fp1[k][j][i].y =
263 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) +
264 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);
265 fp1[k][j][i].z =
266 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) +
267 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);
268 }
269 }
270 }
271 }
272 }
273
274 /* j direction */
275 for (k=lzs; k<lze; k++) {
276 for(j=lys-1; j<lye; j++) {
277 for(i=lxs; i<lxe; i++) {
278 ucon = ucont[k][j][i].y * 0.5;
279
280 up = ucon + fabs(ucon);
281 um = ucon - fabs(ucon);
282
283 if (j>0 && j<my-2 &&
284 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i] > innerblank) &&
285 (nvert[k][j-1][i] < 0.1 || nvert[k][j-1][i] > innerblank)) {
286 if ((les || central)) {
290
291 } else {
292 fp2[k][j][i].x =
293 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) +
294 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);
295 fp2[k][j][i].y =
296 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) +
297 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);
298 fp2[k][j][i].z =
299 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) +
300 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);
301 }
302 }
303 else if ((les || central) && (j==0 || i==my-2) &&
304 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i]>innerblank) &&
305 (nvert[k][j ][i] < 0.1 || nvert[k][j ][i]>innerblank))
306 {
310 }
311 else if (j==0 || (nvert[k][j-1][i]) > 0.1) {
312 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 //
313 fp2[k][j][i].x =
314 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) +
315 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);
316 fp2[k][j][i].y =
317 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) +
318 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);
319 fp2[k][j][i].z =
320 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) +
321 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);
322 }else{
323 fp2[k][j][i].x =
324 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) +
325 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);
326 fp2[k][j][i].y =
327 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) +
328 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);
329 fp2[k][j][i].z =
330 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) +
331 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);
332 }
333 }
334 else if (j==my-2 ||(nvert[k][j+1][i]) > 0.1) {
335 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//
336 fp2[k][j][i].x =
337 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) +
338 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);
339 fp2[k][j][i].y =
340 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) +
341 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);
342 fp2[k][j][i].z =
343 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) +
344 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);
345 }else{
346 fp2[k][j][i].x =
347 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) +
348 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);
349 fp2[k][j][i].y =
350 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) +
351 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);
352 fp2[k][j][i].z =
353 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) +
354 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);
355 }
356 }
357 }
358 }
359 }
360
361
362 /* k direction */
363 for (k=lzs-1; k<lze; k++) {
364 for(j=lys; j<lye; j++) {
365 for(i=lxs; i<lxe; i++) {
366 ucon = ucont[k][j][i].z * 0.5;
367
368 up = ucon + fabs(ucon);
369 um = ucon - fabs(ucon);
370
371 if (k>0 && k<mz-2 &&
372 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i] > innerblank) &&
373 (nvert[k-1][j][i] < 0.1 || nvert[k-1][j][i] > innerblank)) {
374 if ((les || central)) {
378
379 } else {
380 fp3[k][j][i].x =
381 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) +
382 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);
383 fp3[k][j][i].y =
384 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) +
385 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);
386 fp3[k][j][i].z =
387 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) +
388 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);
389 }
390 }
391 else if ((les || central) && (k==0 || k==mz-2) &&
392 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i]>innerblank) &&
393 (nvert[k ][j][i] < 0.1 || nvert[k ][j][i]>innerblank))
394 {
398 }
399 else if (k<mz-2 && (k==0 ||(nvert[k-1][j][i]) > 0.1)) {
400 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//
401 fp3[k][j][i].x =
402 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) +
403 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);
404 fp3[k][j][i].y =
405 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) +
406 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);
407 fp3[k][j][i].z =
408 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) +
409 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);
410 }else{
411 fp3[k][j][i].x =
412 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) +
413 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);
414 fp3[k][j][i].y =
415 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) +
416 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);
417 fp3[k][j][i].z =
418 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) +
419 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);
420 }
421 }
422 else if (k>0 && (k==mz-2 ||(nvert[k+1][j][i]) > 0.1)) {
423 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//
424 fp3[k][j][i].x =
425 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) +
426 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);
427 fp3[k][j][i].y =
428 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) +
429 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);
430 fp3[k][j][i].z =
431 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) +
432 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);
433 }else{
434 fp3[k][j][i].x =
435 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) +
436 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);
437 fp3[k][j][i].y =
438 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) +
439 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);
440 fp3[k][j][i].z =
441 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) +
442 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);
443 }
444 }
445 }
446 }
447 }
448
449 /* Calculate the convective terms under cartesian coordinates */
450
451 for (k=lzs; k<lze; k++) {
452 for (j=lys; j<lye; j++) {
453 for (i=lxs; i<lxe; i++) {
454 conv[k][j][i].x =
455 fp1[k][j][i].x - fp1[k][j][i-1].x +
456 fp2[k][j][i].x - fp2[k][j-1][i].x +
457 fp3[k][j][i].x - fp3[k-1][j][i].x;
458
459 conv[k][j][i].y =
460 fp1[k][j][i].y - fp1[k][j][i-1].y +
461 fp2[k][j][i].y - fp2[k][j-1][i].y +
462 fp3[k][j][i].y - fp3[k-1][j][i].y;
463
464 conv[k][j][i].z =
465 fp1[k][j][i].z - fp1[k][j][i-1].z +
466 fp2[k][j][i].z - fp2[k][j-1][i].z +
467 fp3[k][j][i].z - fp3[k-1][j][i].z;
468 }
469 }
470 }
471 /* for (k=zs; k<ze; k++) { */
472/* for (j=ys; j<ye; j++) { */
473/* for (i=xs; i<xe; i++) { */
474/* if (i==1 && (j==1) && (k==1 || k==21 || k==22|| k==200)) */
475/* 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); */
476/* } */
477/* } */
478/* } */
479
480 DMDAVecRestoreArray(fda, Ucont, &ucont);
481 DMDAVecRestoreArray(fda, Ucat, &ucat);
482 DMDAVecRestoreArray(fda, Conv, &conv);
483 DMDAVecRestoreArray(da, user->lAj, &aj);
484
485 DMDAVecRestoreArray(fda, Fp1, &fp1);
486 DMDAVecRestoreArray(fda, Fp2, &fp2);
487 DMDAVecRestoreArray(fda, Fp3, &fp3);
488 DMDAVecRestoreArray(da, user->lNvert, &nvert);
489
490 VecDestroy(&Fp1);
491 VecDestroy(&Fp2);
492 VecDestroy(&Fp3);
493
494
496
497 PROFILE_FUNCTION_END;
498 return (0);
499}
#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
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:737
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:447
Here is the caller graph for this function:
◆ Viscous()
| PetscErrorCode Viscous | ( | UserCtx * | user, |
| Vec | Ucont, | ||
| Vec | Ucat, | ||
| Vec | Visc | ||
| ) |
Definition at line 504 of file rhs.c.
505{
506
508
509 Cmpnts ***ucont, ***ucat;
510
511 Cmpnts ***csi, ***eta, ***zet;
512 Cmpnts ***icsi, ***ieta, ***izet;
513 Cmpnts ***jcsi, ***jeta, ***jzet;
514 Cmpnts ***kcsi, ***keta, ***kzet;
515
516 PetscReal ***nvert;
517
519 DMDALocalInfo info;
520 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
521 PetscInt mx, my, mz; // Dimensions in three directions
522 PetscInt i, j, k;
523 Vec Fp1, Fp2, Fp3;
524 Cmpnts ***fp1, ***fp2, ***fp3;
525 Cmpnts ***visc;
526 PetscReal ***aj, ***iaj, ***jaj, ***kaj;
527
528 PetscInt lxs, lxe, lys, lye, lzs, lze;
529
530 PetscReal ajc;
531
532 PetscReal dudc, dude, dudz, dvdc, dvde, dvdz, dwdc, dwde, dwdz;
533 PetscReal csi0, csi1, csi2, eta0, eta1, eta2, zet0, zet1, zet2;
534 PetscReal g11, g21, g31;
535 PetscReal r11, r21, r31, r12, r22, r32, r13, r23, r33;
536
537 PetscScalar solid,innerblank;
538
539 // --- CONTEXT ACQUISITION BLOCK ---
540 // Get the master simulation context from the UserCtx.
542
543 // Create local variables to mirror the legacy globals for minimal code changes.
549 solid = 0.5;
550 innerblank = 7.;
551
553
554 DMDAVecGetArray(fda, Ucont, &ucont);
555 DMDAVecGetArray(fda, Ucat, &ucat);
556 DMDAVecGetArray(fda, Visc, &visc);
557
558 DMDAVecGetArray(fda, Csi, &csi);
559 DMDAVecGetArray(fda, Eta, &eta);
560 DMDAVecGetArray(fda, Zet, &zet);
561
562 DMDAVecGetArray(fda, user->lICsi, &icsi);
563 DMDAVecGetArray(fda, user->lIEta, &ieta);
564 DMDAVecGetArray(fda, user->lIZet, &izet);
565
566 DMDAVecGetArray(fda, user->lJCsi, &jcsi);
567 DMDAVecGetArray(fda, user->lJEta, &jeta);
568 DMDAVecGetArray(fda, user->lJZet, &jzet);
569
570 DMDAVecGetArray(fda, user->lKCsi, &kcsi);
571 DMDAVecGetArray(fda, user->lKEta, &keta);
572 DMDAVecGetArray(fda, user->lKZet, &kzet);
573
574 DMDAVecGetArray(da, user->lNvert, &nvert);
575
576 VecDuplicate(Ucont, &Fp1);
577 VecDuplicate(Ucont, &Fp2);
578 VecDuplicate(Ucont, &Fp3);
579
580 DMDAVecGetArray(fda, Fp1, &fp1);
581 DMDAVecGetArray(fda, Fp2, &fp2);
582 DMDAVecGetArray(fda, Fp3, &fp3);
583
584 DMDAVecGetArray(da, user->lAj, &aj);
585
586 DMDAGetLocalInfo(da, &info);
587
588 mx = info.mx; my = info.my; mz = info.mz;
589 xs = info.xs; xe = xs + info.xm;
590 ys = info.ys; ye = ys + info.ym;
591 zs = info.zs; ze = zs + info.zm;
592
593 /* First we calculate the flux on cell surfaces. Stored on the upper integer
594 node. For example, along i direction, the flux are stored at node 0:mx-2*/
595 lxs = xs; lxe = xe;
596 lys = ys; lye = ye;
597 lzs = zs; lze = ze;
598
599 if (xs==0) lxs = xs+1;
600 if (ys==0) lys = ys+1;
601 if (zs==0) lzs = zs+1;
602
603
604 if (xe==mx) lxe=xe-1;
605 if (ye==my) lye=ye-1;
606 if (ze==mz) lze=ze-1;
607
608 VecSet(Visc,0.0);
609
610 PetscReal ***lnu_t;
611
612 if(les) {
613 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
614 } else if (rans) {
615
616 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
617 }
618
619 /* The visc flux on each surface center is stored at previous integer node */
620
621 DMDAVecGetArray(da, user->lIAj, &iaj);
622 /* for (k=zs; k<ze; k++) { */
623/* for (j=ys; j<ye; j++) { */
624/* for (i=xs; i<xe; i++) { */
625/* if (i==1 && (j==0 ||j==1 || j==2) && (k==21 || k==22|| k==20)) */
626/* 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 ); */
627/* } */
628/* } */
629/* } */
630 // i direction
631 for (k=lzs; k<lze; k++) {
632 for (j=lys; j<lye; j++) {
633 for (i=lxs-1; i<lxe; i++) {
634
638
639 if ((nvert[k][j+1][i ]> solid && nvert[k][j+1][i ]<innerblank) ||
640 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
647 }
648 else if ((nvert[k][j-1][i ]> solid && nvert[k][j-1][i ]<innerblank) ||
649 (nvert[k][j-1][i+1]> solid && nvert[k][j-1][i+1]<innerblank)) {
656 }
657 else {
664 }
665
666 if ((nvert[k+1][j][i ]> solid && nvert[k+1][j][i ]<innerblank)||
667 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
674 }
675 else if ((nvert[k-1][j][i ]> solid && nvert[k-1][j][i ]<innerblank) ||
676 (nvert[k-1][j][i+1]> solid && nvert[k-1][j][i+1]<innerblank)) {
677
684 }
685 else {
692 }
693
694 csi0 = icsi[k][j][i].x;
695 csi1 = icsi[k][j][i].y;
696 csi2 = icsi[k][j][i].z;
697
698 eta0 = ieta[k][j][i].x;
699 eta1 = ieta[k][j][i].y;
700 eta2 = ieta[k][j][i].z;
701
702 zet0 = izet[k][j][i].x;
703 zet1 = izet[k][j][i].y;
704 zet2 = izet[k][j][i].z;
705
706 g11 = csi0 * csi0 + csi1 * csi1 + csi2 * csi2;
707 g21 = eta0 * csi0 + eta1 * csi1 + eta2 * csi2;
708 g31 = zet0 * csi0 + zet1 * csi1 + zet2 * csi2;
709
710 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
711 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
712 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
713
714 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
715 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
716 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
717
718 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
719 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
720 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
721
722 ajc = iaj[k][j][i];
723
725
726 if( les || (rans && ti>0) ) {
727 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j][i+1]) ), 2.0) * Sabs;
728 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j][i+1]);
729 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;
730 fp1[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * csi0 + r21 * csi1 + r31 * csi2) * ajc * (nu_t);
731 fp1[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * csi0 + r22 * csi1 + r32 * csi2) * ajc * (nu_t);
732 fp1[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * csi0 + r23 * csi1 + r33 * csi2) * ajc * (nu_t);
733 }
734 else {
735 fp1[k][j][i].x = 0;
736 fp1[k][j][i].y = 0;
737 fp1[k][j][i].z = 0;
738 }
739
740 fp1[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * csi0 + r21 * csi1 + r31 * csi2 ) * ajc * (nu);
741 fp1[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * csi0 + r22 * csi1 + r32 * csi2 ) * ajc * (nu);
742 fp1[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * csi0 + r23 * csi1 + r33 * csi2 ) * ajc * (nu);
743
744
745 if(clark) {
746 double dc, de, dz;
747 ComputeCellCharacteristicLengthScale (ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
748 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
749
750 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
751 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
752 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
753 double t21 = t12;
754 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
755 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
756 double t31 = t13;
757 double t32 = t23;
758 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
759
760 fp1[k][j][i].x -= ( t11 * csi0 + t12 * csi1 + t13 * csi2 ) / 12.;
761 fp1[k][j][i].y -= ( t21 * csi0 + t22 * csi1 + t23 * csi2 ) / 12.;
762 fp1[k][j][i].z -= ( t31 * csi0 + t32 * csi1 + t33 * csi2 ) / 12.;
763 }
764
765 }
766 }
767 }
768 DMDAVecRestoreArray(da, user->lIAj, &iaj);
769
770
771 // j direction
772 DMDAVecGetArray(da, user->lJAj, &jaj);
773 for (k=lzs; k<lze; k++) {
774 for (j=lys-1; j<lye; j++) {
775 for (i=lxs; i<lxe; i++) {
776
777 if ((nvert[k][j ][i+1]> solid && nvert[k][j ][i+1]<innerblank)||
778 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
785 }
786 else if ((nvert[k][j ][i-1]> solid && nvert[k][j ][i-1]<innerblank) ||
787 (nvert[k][j+1][i-1]> solid && nvert[k][j+1][i-1]<innerblank)) {
794 }
795 else {
802 }
803
807
808 if ((nvert[k+1][j ][i]> solid && nvert[k+1][j ][i]<innerblank)||
809 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
816 }
817 else if ((nvert[k-1][j ][i]> solid && nvert[k-1][j ][i]<innerblank)||
818 (nvert[k-1][j+1][i]> solid && nvert[k-1][j+1][i]<innerblank)) {
825 }
826 else {
833 }
834
835 csi0 = jcsi[k][j][i].x;
836 csi1 = jcsi[k][j][i].y;
837 csi2 = jcsi[k][j][i].z;
838
839 eta0 = jeta[k][j][i].x;
840 eta1 = jeta[k][j][i].y;
841 eta2 = jeta[k][j][i].z;
842
843 zet0 = jzet[k][j][i].x;
844 zet1 = jzet[k][j][i].y;
845 zet2 = jzet[k][j][i].z;
846
847
848 g11 = csi0 * eta0 + csi1 * eta1 + csi2 * eta2;
849 g21 = eta0 * eta0 + eta1 * eta1 + eta2 * eta2;
850 g31 = zet0 * eta0 + zet1 * eta1 + zet2 * eta2;
851
852 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
853 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
854 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
855
856 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
857 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
858 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
859
860 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
861 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
862 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
863
864 // 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);
865 // 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);
866 // 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);
867 // 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);
868
869
870
871 ajc = jaj[k][j][i];
872
874
875 if( les || (rans && ti>0) ) {
876 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j+1][i]) ), 2.0) * Sabs;
877 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j+1][i]);
878 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;
879
880 fp2[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * eta0 + r21 * eta1 + r31 * eta2) * ajc * (nu_t);
881 fp2[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * eta0 + r22 * eta1 + r32 * eta2) * ajc * (nu_t);
882 fp2[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * eta0 + r23 * eta1 + r33 * eta2) * ajc * (nu_t);
883 }
884 else {
885 fp2[k][j][i].x = 0;
886 fp2[k][j][i].y = 0;
887 fp2[k][j][i].z = 0;
888 }
889
890 fp2[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * eta0 + r21 * eta1 + r31 * eta2 ) * ajc * (nu);
891 fp2[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * eta0 + r22 * eta1 + r32 * eta2 ) * ajc * (nu);
892 fp2[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * eta0 + r23 * eta1 + r33 * eta2 ) * ajc * (nu);
893
894 if(clark) {
895 double dc, de, dz;
896 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
897 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
898
899 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
900 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
901 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
902 double t21 = t12;
903 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
904 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
905 double t31 = t13;
906 double t32 = t23;
907 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
908
909 fp2[k][j][i].x -= ( t11 * eta0 + t12 * eta1 + t13 * eta2 ) / 12.;
910 fp2[k][j][i].y -= ( t21 * eta0 + t22 * eta1 + t23 * eta2 ) / 12.;
911 fp2[k][j][i].z -= ( t31 * eta0 + t32 * eta1 + t33 * eta2 ) / 12.;
912 }
913 }
914 }
915 }
916
917 DMDAVecRestoreArray(da, user->lJAj, &jaj);
918 // k direction
919
920 DMDAVecGetArray(da, user->lKAj, &kaj);
921 for (k=lzs-1; k<lze; k++) {
922 for (j=lys; j<lye; j++) {
923 for (i=lxs; i<lxe; i++) {
924 if ((nvert[k ][j][i+1]> solid && nvert[k ][j][i+1]<innerblank)||
925 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
932 }
933 else if ((nvert[k ][j][i-1]> solid && nvert[k ][j][i-1]<innerblank) ||
934 (nvert[k+1][j][i-1]> solid && nvert[k+1][j][i-1]<innerblank)) {
941 }
942 else {
949 }
950
951 if ((nvert[k ][j+1][i]> solid && nvert[k ][j+1][i]<innerblank)||
952 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
959 }
960 else if ((nvert[k ][j-1][i]> solid && nvert[k ][j-1][i]<innerblank) ||
961 (nvert[k+1][j-1][i]> solid && nvert[k+1][j-1][i]<innerblank)){
968 }
969 else {
976 }
977
981
982
983 csi0 = kcsi[k][j][i].x;
984 csi1 = kcsi[k][j][i].y;
985 csi2 = kcsi[k][j][i].z;
986
987 eta0 = keta[k][j][i].x;
988 eta1 = keta[k][j][i].y;
989 eta2 = keta[k][j][i].z;
990
991 zet0 = kzet[k][j][i].x;
992 zet1 = kzet[k][j][i].y;
993 zet2 = kzet[k][j][i].z;
994
995
996 g11 = csi0 * zet0 + csi1 * zet1 + csi2 * zet2;
997 g21 = eta0 * zet0 + eta1 * zet1 + eta2 * zet2;
998 g31 = zet0 * zet0 + zet1 * zet1 + zet2 * zet2;
999
1000 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
1001 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
1002 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
1003
1004 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
1005 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
1006 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
1007
1008 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
1009 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
1010 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
1011
1012 ajc = kaj[k][j][i];
1013
1015
1016 if( les || (rans && ti>0) ) {
1017 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k+1][j][i]) ), 2.0) * Sabs;
1018 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k+1][j][i]);
1019 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;
1020
1021 fp3[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu_t);
1022 fp3[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu_t);
1023 fp3[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu_t);
1024 }
1025 else {
1026 fp3[k][j][i].x = 0;
1027 fp3[k][j][i].y = 0;
1028 fp3[k][j][i].z = 0;
1029 }
1030 fp3[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu);//
1031 fp3[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu);//
1032 fp3[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu);//
1033
1034 if(clark) {
1035 double dc, de, dz;
1036 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
1037 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
1038
1039 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
1040 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
1041 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
1042 double t21 = t12;
1043 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
1044 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
1045 double t31 = t13;
1046 double t32 = t23;
1047 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
1048
1049 fp3[k][j][i].x -= ( t11 * zet0 + t12 * zet1 + t13 * zet2 ) / 12.;
1050 fp3[k][j][i].y -= ( t21 * zet0 + t22 * zet1 + t23 * zet2 ) / 12.;
1051 fp3[k][j][i].z -= ( t31 * zet0 + t32 * zet1 + t33 * zet2 ) / 12.;
1052 }
1053 }
1054 }
1055 }
1056
1057 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1058
1059 for (k=lzs; k<lze; k++) {
1060 for (j=lys; j<lye; j++) {
1061 for (i=lxs; i<lxe; i++) {
1062 visc[k][j][i].x =
1066
1067 visc[k][j][i].y =
1071
1072 visc[k][j][i].z =
1076
1077 }
1078 }
1079 }
1080/* for (k=zs; k<ze; k++) { */
1081/* for (j=ys; j<ye; j++) { */
1082/* for (i=xs; i<xe; i++) { */
1083/* 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); */
1084/* 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); */
1085/* 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); */
1086/* 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); */
1087/* 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); */
1088/* 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); */
1089
1090/* } */
1091/* } */
1092/* } */
1093 DMDAVecRestoreArray(fda, Ucont, &ucont);
1094 DMDAVecRestoreArray(fda, Ucat, &ucat);
1095 DMDAVecRestoreArray(fda, Visc, &visc);
1096
1097 DMDAVecRestoreArray(fda, Csi, &csi);
1098 DMDAVecRestoreArray(fda, Eta, &eta);
1099 DMDAVecRestoreArray(fda, Zet, &zet);
1100
1101 DMDAVecRestoreArray(fda, Fp1, &fp1);
1102 DMDAVecRestoreArray(fda, Fp2, &fp2);
1103 DMDAVecRestoreArray(fda, Fp3, &fp3);
1104
1105 DMDAVecRestoreArray(da, user->lAj, &aj);
1106
1107 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1108 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1109 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1110
1111 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1112 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1113 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1114
1115 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1116 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1117 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1118
1119 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1120
1121 if(les) {
1122 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1123 } else if (rans) {
1124
1125 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1126 }
1127
1128
1129 VecDestroy(&Fp1);
1130 VecDestroy(&Fp2);
1131 VecDestroy(&Fp3);
1132
1133
1135
1136 PROFILE_FUNCTION_END;
1137
1138 return(0);
1139}
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:313
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:423
Here is the call graph for this function:
Here is the caller graph for this function:
◆ ComputeBodyForces()
| PetscErrorCode ComputeBodyForces | ( | UserCtx * | user, |
| Vec | Rct | ||
| ) |
General dispatcher for applying all active body forces (momentum sources).
This function serves as a central hub for adding momentum source terms to the contravariant right-hand-side (Rct) of the momentum equations. It is called once per RHS evaluation (e.g., once per Runge-Kutta stage).
It introspects the simulation configuration to determine which, if any, body forces are active and calls their specific implementation functions.
- Parameters
-
user The UserCtx containing the simulation state for a single block. Rhs The PETSc Vec for the contravariant RHS, which will be modified in-place.
- Returns
- PetscErrorCode 0 on success.
Definition at line 1157 of file rhs.c.
1158{
1159 PetscErrorCode ierr;
1160 PetscFunctionBeginUser;
1161
1162 // --- 1. Apply momentum source for driven channel/pipe flows ---
1163 // This function will internally check if a driven flow BC is active.
1164 ierr = ComputeDrivenChannelFlowSource(user, Rct); CHKERRQ(ierr);
1165
1166 // --- 2. (Future Extension) Apply gravitational force ---
1167 // if (user->simCtx->gravityEnabled) {
1168 // ierr = ApplyGravitationalForce(user, Rhs); CHKERRQ(ierr);
1169 // }
1170
1171 PetscFunctionReturn(0);
1172}
PetscErrorCode ComputeDrivenChannelFlowSource(UserCtx *user, Vec Rct)
Applies a momentum source term to drive flow in a periodic channel or pipe.
Definition BodyForces.c:29
Here is the call graph for this function:
Here is the caller graph for this function:
◆ ComputeRHS()
| PetscErrorCode ComputeRHS | ( | UserCtx * | user, |
| Vec | Rhs | ||
| ) |
Computes the Right-Hand-Side (RHS) of the momentum equations.
Computes the Right-Hand Side (RHS) of the momentum equations.
This function orchestrates the calculation of the spatial discretization of the convective and diffusive terms. It calls specialized helper functions (Convection, Viscous) to compute these terms and then combines them with the pressure gradient to form the final RHS vector.
This is a minimally-edited version of the legacy kernel. It operates on a single UserCtx (one grid block) and retrieves all global parameters (Re, rans, les, etc.) from the master SimCtx.
- Parameters
-
user The UserCtx for the specific block being computed. Rhs The PETSc Vec where the calculated RHS will be stored.
- Returns
- PetscErrorCode 0 on success.
Definition at line 1192 of file rhs.c.
1193{
1194 PetscErrorCode ierr;
1197 DMDALocalInfo info = user->info;
1198 PetscInt i,j,k;
1199 PetscReal dpdc, dpde,dpdz;
1200 // --- Local Grid Indices and Parameters ---
1201 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1202 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1203 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1204 PetscInt lxs = (xs==0) ? xs+1 : xs;
1205 PetscInt lys = (ys==0) ? ys+1 : ys;
1206 PetscInt lzs = (zs==0) ? zs+1 : zs;
1207 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1208 PetscInt lye = (ye==my) ? ye-1 : ye;
1209 PetscInt lze = (ze==mz) ? ze-1 : ze;
1210 PetscInt mz_end = (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == CHARACTERISTIC_BC) ? mz - 2 : mz - 3;
1211
1212 // --- Array Pointers ---
1213 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1214 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert, ***nvert_o;
1215 Cmpnts ***rhs, ***rc, ***rct,***ucont;
1216
1217 // --- Temporary Vectors ---
1218 Vec Conv, Visc, Rc, Rct;
1219
1220 PetscFunctionBeginUser;
1221 PROFILE_FUNCTION_BEGIN;
1224
1225 // --- Get all necessary array pointers ---
1226 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1227 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1228 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1229 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1230 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1231 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1232 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1233 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1234 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1235 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1236 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1237 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1238 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1239 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1240 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1241 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1242 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1243 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1244 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1245
1246 // --- Create temporary work vectors ---
1247 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1248 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1249 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1250 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1251
1252 // ========================================================================
1253 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1254 // ========================================================================
1255
1256 // 1. Obtain Cartesian velocity from Contravariant velocity
1257 ierr = Contra2Cart(user); CHKERRQ(ierr);
1259
1260 // 2. Compute Convective term
1263 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1264 } else {
1266 }
1267
1268 // 3. Compute Viscous term
1270 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1271 } else {
1274 }
1275
1276 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1277 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1278
1279 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1281 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1282 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1283
1284 for (k = lzs; k < lze; k++) {
1285 for (j = lys; j < lye; j++) {
1286 for (i = lxs; i < lxe; i++) {
1287 rct[k][j][i].x = aj[k][j][i] *
1289 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1291 rct[k][j][i].y = aj[k][j][i] *
1293 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1295 rct[k][j][i].z = aj[k][j][i] *
1297 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1299 }
1300 }
1301 }
1302 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1303 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1304
1305 PetscBarrier(NULL);
1306
1307 DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct);
1308 DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct);
1309
1310 // Compute and Add Body Force term if applicable.
1311 ierr = ComputeBodyForces(user,Rct);
1312
1313 DMDAVecGetArray(fda, Rct, &rct);
1314
1315
1316 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
1317 if (xs==0){
1318 i=xs;
1319 for (k=lzs; k<lze; k++) {
1320 for (j=lys; j<lye; j++) {
1322 }
1323 }
1324 }
1325
1326 if (xe==mx){
1327 i=mx-1;
1328 for (k=lzs; k<lze; k++) {
1329 for (j=lys; j<lye; j++) {
1331 }
1332 }
1333 }
1334 }
1335
1336 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
1337 if (ys==0){
1338 j=ys;
1339 for (k=lzs; k<lze; k++) {
1340 for (i=lxs; i<lxe; i++) {
1342 }
1343 }
1344 }
1345
1346 if (ye==my){
1347 j=my-1;
1348 for (k=lzs; k<lze; k++) {
1349 for (i=lxs; i<lxe; i++) {
1351 }
1352 }
1353 }
1354 }
1355 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
1356 if (zs==0){
1357 k=zs;
1358 for (j=lys; j<lye; j++) {
1359 for (i=lxs; i<lxe; i++) {
1361 }
1362 }
1363 }
1364 if (ze==mz){
1365 k=mz-1;
1366 for (j=lys; j<lye; j++) {
1367 for (i=lxs; i<lxe; i++) {
1369 }
1370 }
1371 }
1372 }
1373
1374 // 6. Add Pressure Gradient Term and Finalize RHS
1375 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1376 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1378
1379 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1380
1381 ierr = DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1382 ierr = DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1383
1384 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1385
1386 for (k = lzs; k < lze; k++) {
1387 for (j = lys; j < lye; j++) {
1388 for (i = lxs; i < lxe; i++) {
1389 PetscReal dpdc, dpde, dpdz;
1390 dpdc = p[k][j][i+1] - p[k][j][i];
1391
1392 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) {
1393 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))) {
1394 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1395 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1396 }
1397 }
1398 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) {
1399 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1400 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1401 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1402 }
1403 }
1404 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) {
1405 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1406 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1407 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1408 }
1409 }
1410 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) {
1411 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1412 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1413 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1414 }
1415 }
1416 else {
1417 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1418 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1419 }
1420
1421 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) {
1422 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))) {
1423 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1424 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1425 }
1426 }
1427 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) {
1428 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1429 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1430 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1431 }
1432 }
1433 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) {
1434 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1435 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1436 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1437 }
1438 }
1439 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) {
1440 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1441 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1442 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1443 }
1444 }
1445 else {
1446 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1447 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1448 }
1449
1451
1452
1453 rhs[k][j][i].x -=
1457 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1458 ieta[k][j][i].y * icsi[k][j][i].y +
1459 ieta[k][j][i].z * icsi[k][j][i].z)+
1463
1464 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) {
1465 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))) {
1466 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1467 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1468 }
1469 }
1470 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) {
1471 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1472 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1473 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1474 }
1475 }
1476 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) {
1477 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1478 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1479 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1480 }
1481 }
1482 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) {
1483 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1484 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1485 p[k][j][i] - p[k][j+1][i]) * 0.5;
1486 }
1487 }
1488 else {
1489 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1490 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1491 }
1492
1493 dpde = p[k][j+1][i] - p[k][j][i];
1494
1495 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) {
1496 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))) {
1497 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1498 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1499 }
1500 }
1501 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) {
1502 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1503 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1504 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1505 }
1506 }
1507 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) {
1508 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1509 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1510 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1511 }
1512 }
1513 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) {
1514 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1515 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1516 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1517 }
1518 }
1519 else {
1520 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1521 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1522 }
1523
1525
1526
1527 rhs[k][j][i].y -=
1531 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1532 jeta[k][j][i].y * jeta[k][j][i].y +
1533 jeta[k][j][i].z * jeta[k][j][i].z) +
1537
1538 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) {
1539 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))) {
1540 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1541 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1542 }
1543 }
1544 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) {
1545 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1546 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1547 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1548 }
1549 }
1550 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) {
1551 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1552 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1553 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1554 }
1555 }
1556 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) {
1557 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1558 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1559 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1560 }
1561 }
1562 else {
1563 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1564 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1565 }
1566
1567 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) {
1568 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))) {
1569 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1570 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1571 }
1572 }
1573 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) {
1574 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1575 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1576 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1577 }
1578 }
1579 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) {
1580 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1581 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1582 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1583 }
1584 }
1585 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) {
1586 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1587 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1588 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1589 }
1590 }
1591 else {
1592 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1593 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1594 }
1595
1596 dpdz = (p[k+1][j][i] - p[k][j][i]);
1597
1599
1600 rhs[k][j][i].z -=
1604 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1605 keta[k][j][i].y * kzet[k][j][i].y +
1606 keta[k][j][i].z * kzet[k][j][i].z) +
1610
1611 }
1612 }
1613 }
1614
1615
1616 //Mohsen March 2012//
1617
1618 // rhs.x at boundaries for periodic bc at i direction//
1620 for (k=lzs; k<lze; k++) {
1621 for (j=lys; j<lye; j++) {
1622 i=xs;
1623
1624 dpdc = p[k][j][i+1] - p[k][j][i];
1625
1626 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) {
1627 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))) {
1628 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1629 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1630 }
1631 }
1632 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) {
1633 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1634 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1635 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1636 }
1637 }
1638 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) {
1639 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1640 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1641 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1642 }
1643 }
1644 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) {
1645 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1646 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1647 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1648 }
1649 }
1650 else {
1651 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1652 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1653 }
1654
1655 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) {
1656 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))) {
1657 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1658 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1659 }
1660 }
1661 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) {
1662 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1663 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1664 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1665 }
1666 }
1667 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) {
1668 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1669 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1670 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1671 }
1672 }
1673 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) {
1674 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1675 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1676 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1677 }
1678 }
1679 else {
1680 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1681 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1682 }
1683
1685 rhs[k][j][i].x -=
1689 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1690 ieta[k][j][i].y * icsi[k][j][i].y +
1691 ieta[k][j][i].z * icsi[k][j][i].z)+
1695 }
1696 }
1697 }
1698
1699// rhs.y at boundaries for periodic bc at j direction//
1701 for (k=lzs; k<lze; k++) {
1702 for (i=lxs; i<lxe; i++) {
1703
1704 j=ys;
1705
1706 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) {
1707 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))) {
1708 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1709 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1710 }
1711 }
1712 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) {
1713 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1714 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1715 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1716 }
1717 }
1718 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) {
1719 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1720 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1721 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1722 }
1723 }
1724 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) {
1725 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1726 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1727 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1728 }
1729 }
1730 else {
1731 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1732 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1733 }
1734
1735 dpde = p[k][j+1][i] - p[k][j][i];
1736
1737 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) {
1738 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))) {
1739 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1740 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1741 }
1742 }
1743 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) {
1744 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1745 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1746 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1747 }
1748 }
1749 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) {
1750 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1751 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1752 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1753 }
1754 }
1755 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) {
1756 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1757 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1758 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1759 }
1760 }
1761 else {
1762 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1763 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1764 }
1765
1767
1768 rhs[k][j][i].y -=
1772 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1773 jeta[k][j][i].y * jeta[k][j][i].y +
1774 jeta[k][j][i].z * jeta[k][j][i].z)+
1778
1779 }
1780 }
1781 }
1782
1783 // rhs.z at boundaries for periodic bc at k direction//
1785 for (j=lys; j<lye; j++) {
1786 for (i=lxs; i<lxe; i++) {
1787
1788 k=zs;
1789
1790 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) {
1791 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))) {
1792 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1793 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1794 }
1795 }
1796 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) {
1797 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1798 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1799 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1800 }
1801 }
1802 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) {
1803 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1804 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1805 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1806 }
1807 }
1808 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) {
1809 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1810 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1811 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1812 }
1813 }
1814 else {
1815 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1816 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1817 }
1818
1819 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) {
1820 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))) {
1821 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1822 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1823 }
1824 }
1825 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) {
1826 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1827 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1828 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1829 }
1830 }
1831 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) {
1832 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1833 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1834 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1835 }
1836 }
1837 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) {
1838 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1839 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1840 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1841 }
1842 }
1843 else {
1844 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1845 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1846 }
1847
1848 dpdz = (p[k+1][j][i] - p[k][j][i]);
1849
1851
1852 rhs[k][j][i].z -=
1856 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1857 keta[k][j][i].y * kzet[k][j][i].y +
1858 keta[k][j][i].z * kzet[k][j][i].z)+
1862
1863 }
1864 }
1865 }
1866
1867 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1868
1870 PetscInt TwoD = simCtx->TwoD;
1871
1873
1874 // 7. Final clean-up for immersed boundaries and 2D cases
1875 for (k=lzs; k<lze; k++) {
1876 for (j=lys; j<lye; j++) {
1877 for (i=lxs; i<lxe; i++) {
1878 if (TwoD==1)
1879 rhs[k][j][i].x =0.;
1880 else if (TwoD==2)
1881 rhs[k][j][i].y =0.;
1882 else if (TwoD==3)
1883 rhs[k][j][i].z =0.;
1884
1885 if (nvert[k][j][i]>0.1) {
1886 rhs[k][j][i].x = 0;
1887 rhs[k][j][i].y = 0;
1888 rhs[k][j][i].z = 0;
1889 }
1890 if (nvert[k][j][i+1]>0.1) {
1891 rhs[k][j][i].x=0;
1892 }
1893 if (nvert[k][j+1][i]>0.1) {
1894 rhs[k][j][i].y=0;
1895 }
1896 if (nvert[k+1][j][i]>0.1) {
1897 rhs[k][j][i].z=0;
1898 }
1899 }
1900 }
1901 }
1903
1904 // ========================================================================
1905
1906 // --- Restore all PETSc array pointers ---
1907 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1908 DMDAVecRestoreArray(fda, Rhs, &rhs);
1910
1911 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1912 DMDAVecRestoreArray(fda, user->lEta, &eta);
1913 DMDAVecRestoreArray(fda, user->lZet, &zet);
1914 DMDAVecRestoreArray(da, user->lAj, &aj);
1916
1917 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1918 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1919 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1920 DMDAVecRestoreArray(da, user->lIAj, &iaj);
1922
1923 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1924 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1925 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1926 DMDAVecRestoreArray(da, user->lJAj, &jaj);
1928
1929 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1930 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1931 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1932 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1934
1935 DMDAVecRestoreArray(da, user->lP, &p);
1937
1938 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1940
1942
1943 // --- Destroy temporary work vectors ---
1944 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1945 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1946 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1947 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1949
1952
1953 PROFILE_FUNCTION_END;
1954 PetscFunctionReturn(0);
1955}
PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
General dispatcher for applying all active body forces (momentum sources).
Definition rhs.c:1157
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2177
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1157
Here is the call graph for this function:
Here is the caller graph for this function:
◆ ComputeEulerianDiffusivity()
| PetscErrorCode ComputeEulerianDiffusivity | ( | UserCtx * | user | ) |
Computes the effective diffusivity scalar field ($\Gamma_{eff}$) on the Eulerian grid.
This function calculates the total diffusivity used to drive the stochastic motion of particles (Scalar FDF). It combines molecular diffusion and turbulent diffusion.
Formula:
![\[
\Gamma_{eff} = \underbrace{\frac{\nu}{Sc}}_{\text{Molecular}} + \underbrace{\frac{\nu_t}{Sc_t}}_{\text{Turbulent}}
\]](form_0_dark.png)
Where:
- Note
- If turbulence models are disabled,
- This function updates the local ghost values of lDiffusivity at the end to ensure gradients can be computed correctly at subdomain boundaries.
- Parameters
-
[in,out] user Pointer to the user context containing grid data and simulation parameters.
- Returns
- PetscErrorCode 0 on success.
Definition at line 1985 of file rhs.c.
1986{
1987 PetscErrorCode ierr;
1988 DM da = user->da;
1989 PetscInt i, j, k, xs, ys, zs, xm, ym, zm, xe, ye, ze;
1990 PetscInt lxs, lys, lzs, lxe, lye, lze;
1991
1992 // Pointers for 3D grid access
1993 PetscReal ***diff_arr; // Output: Diffusivity field
1994 PetscReal ***nut_arr; // Input: Eddy Viscosity field (optional)
1995
1996 // Physics parameters
1997 PetscReal nu_molecular, gamma_molecular;
1998 PetscReal nu_turbulent, gamma_turbulent;
1999 PetscReal Sc, Sct;
2000 PetscBool use_turbulence_model;
2001
2002 PetscFunctionBeginUser;
2003 PROFILE_FUNCTION_BEGIN;
2004
2005 // ------------------------------------------------------------------------
2006 // 1. Parameter Setup & Safety Checks
2007 // ------------------------------------------------------------------------
2008
2009 // Determine Molecular Viscosity (nu = 1/Re)
2010 // Guard against division by zero if Re is not set or infinite (inviscid)
2013 } else {
2014 nu_molecular = 0.0;
2015 }
2016
2017 // Set Schmidt Numbers (Default to 1.0 if not provided to prevent NaN)
2019 Sct = (user->simCtx->Turbulent_schmidt_number > 1.0e-6) ? user->simCtx->Turbulent_schmidt_number : 0.7;
2020
2021 // Pre-calculate molecular component
2022 gamma_molecular = nu_molecular / Sc;
2023
2024 // Check if a turbulence model is active (LES or RANS)
2026
2027 // ------------------------------------------------------------------------
2028 // 2. Data Access
2029 // ------------------------------------------------------------------------
2030
2031 // Get local grid boundaries
2032 DMDALocalInfo info;
2033 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
2034
2035 xs = info.xs; ys = info.ys; zs = info.zs;
2036 xm = info.xm; ym = info.ym; zm = info.zm;
2037 xe = xs + xm; ye = ys + ym; ze = zs + zm;
2038
2039 lxs = (xs == 0)? xs + 1 : xs;
2040 lys = (ys == 0)? ys + 1 : ys;
2041 lzs = (zs == 0)? zs + 1 : zs;
2042 lxe = (xe == info.mx)? xe - 1 : xe;
2043 lye = (ye == info.my)? ye - 1 : ye;
2044 lze = (ze == info.mz)? ze - 1 : ze;
2045
2046 // Get write access to the output Diffusivity array
2047 ierr = DMDAVecGetArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
2048
2049 // Get read access to Eddy Viscosity only if turbulence is active
2050 if (use_turbulence_model) {
2051 ierr = DMDAVecGetArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
2052 }
2053
2054 // ------------------------------------------------------------------------
2055 // 3. Calculation Loop
2056 // ------------------------------------------------------------------------
2057
2058 for (k = lzs; k < lze; k++) {
2059 for (j = lys; j < lye; j++) {
2060 for (i = lxs; i < lxe; i++) {
2061
2062 gamma_turbulent = 0.0;
2063
2064 if (use_turbulence_model) {
2065 // Fetch local eddy viscosity
2066 nu_turbulent = nut_arr[k][j][i];
2067
2068 // NUMERICAL SAFETY:
2069 // Some turbulence models (dynamic SGS) can locally produce
2070 // slightly negative viscosity. We clamp this to 0 to prevent
2071 // negative diffusivity, which crashes the Langevin sqrt().
2072 if (nu_turbulent < 0.0) {
2073 nu_turbulent = 0.0;
2074 }
2075
2076 gamma_turbulent = nu_turbulent / Sct;
2077 }
2078
2079 // Sum components
2080 diff_arr[k][j][i] = gamma_molecular + gamma_turbulent;
2081 }
2082 }
2083 }
2084
2085 // ------------------------------------------------------------------------
2086 // 4. Cleanup & Synchronization
2087 // ------------------------------------------------------------------------
2088
2089 // Restore arrays
2090 if (use_turbulence_model) {
2091 ierr = DMDAVecRestoreArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
2092 }
2093 ierr = DMDAVecRestoreArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
2094
2095 // Update Ghost Points
2096 // This is required because downstream operations (Drift Gradient Calculation
2097 // and Particle Interpolation) will need access to the halo regions of this field.
2099 PROFILE_FUNCTION_END;
2100 PetscFunctionReturn(0);
2101}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ ComputeEulerianDiffusivityGradient()
| PetscErrorCode ComputeEulerianDiffusivityGradient | ( | UserCtx * | user | ) |
Computes the gradient of the scalar Diffusivity field (Drift Vector).
Logic:
- Iterates over the PHYSICAL cell centers [1 ... M-2].
- Uses 2nd-Order Central Difference for interior cells.
- Uses 2nd-Order One-Sided Differences at physical boundaries (Forward at start, Backward at end) unless the boundary is Periodic.
- Transforms results to Cartesian coordinates using the Chain Rule.
Definition at line 2115 of file rhs.c.
2116{
2117 PetscErrorCode ierr;
2119 DMDALocalInfo info = user->info;
2120
2121 // 1. Determine Global Dimensions
2122 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2123
2124 // 2. Determine Local Loop Bounds (Skip Ghosts/Unused Indices)
2125 // Grid uses indices 1 to M-2 for physical cells.
2126 // Index 0 and M-1 are ghost/boundary holders.
2127
2128 // Start: If we own the global start (0), skip it and start at 1.
2129 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2130 // End: If we own the global end (mx), stop before it (mx-1), so loop covers mx-2.
2131 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
2132
2133 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2134 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
2135
2136 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2137 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
2138
2139 PetscInt i, j, k;
2140
2141 // Pointers
2142 PetscReal ***diff; // Input (Scalar)
2144 Cmpnts ***csi, ***eta, ***zet;
2145 PetscReal ***aj;
2146
2147 // Boundary Flags (Check if Periodic)
2151
2152 PetscFunctionBeginUser;
2153 PROFILE_FUNCTION_BEGIN;
2154
2155 // 3. Update Ghosts for Diffusivity
2156 // Required so that Central Differences at i=2 can read i=1,
2157 // and Central Differences at Periodic Boundaries work correctly.
2159
2160 // 4. Get Arrays (Read Only)
2161 ierr = DMDAVecGetArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2162 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2163 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2164 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2165 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2166
2168
2169 // 5. Get Output Array (Read/Write)
2170 ierr = DMDAVecGetArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2171
2173
2174 // 6. Loop over Physical Domain
2175 for (k = lzs; k < lze; k++) {
2176 for (j = lys; j < lye; j++) {
2177 for (i = lxs; i < lxe; i++) {
2178
2179 PetscReal dGdCsi, dGdEta, dGdZet;
2180
2181 // ---------------------------------------------------------
2182 // I-Direction (Csi)
2183 // ---------------------------------------------------------
2184 if (!p_x && i == 1) {
2185 // Physical Start: 2nd Order Forward Difference
2186 // Stencil: [-3, 4, -1] / 2
2187 dGdCsi = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j][i+1] - diff[k][j][i+2]) * 0.5;
2188 }
2189 else if (!p_x && i == mx - 2) {
2190 // Physical End: 2nd Order Backward Difference
2191 // Stencil: [1, -4, 3] / 2
2192 dGdCsi = (3.0 * diff[k][j][i] - 4.0 * diff[k][j][i-1] + diff[k][j][i-2]) * 0.5;
2193 }
2194 else {
2195 // Interior / Periodic: Central Difference
2196 // Stencil: [-1, 0, 1] / 2
2197 dGdCsi = (diff[k][j][i+1] - diff[k][j][i-1]) * 0.5;
2198 }
2199
2200 // ---------------------------------------------------------
2201 // J-Direction (Eta)
2202 // ---------------------------------------------------------
2203 if (!p_y && j == 1) {
2204 // Physical Start: Forward
2205 dGdEta = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j+1][i] - diff[k][j+2][i]) * 0.5;
2206 }
2207 else if (!p_y && j == my - 2) {
2208 // Physical End: Backward
2209 dGdEta = (3.0 * diff[k][j][i] - 4.0 * diff[k][j-1][i] + diff[k][j-2][i]) * 0.5;
2210 }
2211 else {
2212 // Interior: Central
2213 dGdEta = (diff[k][j+1][i] - diff[k][j-1][i]) * 0.5;
2214 }
2215
2216 // ---------------------------------------------------------
2217 // K-Direction (Zet)
2218 // ---------------------------------------------------------
2219 if (!p_z && k == 1) {
2220 // Physical Start: Forward
2221 dGdZet = (-3.0 * diff[k][j][i] + 4.0 * diff[k+1][j][i] - diff[k+2][j][i]) * 0.5;
2222 }
2223 else if (!p_z && k == mz - 2) {
2224 // Physical End: Backward
2225 dGdZet = (3.0 * diff[k][j][i] - 4.0 * diff[k-1][j][i] + diff[k-2][j][i]) * 0.5;
2226 }
2227 else {
2228 // Interior: Central
2229 dGdZet = (diff[k+1][j][i] - diff[k-1][j][i]) * 0.5;
2230 }
2231
2232 // ---------------------------------------------------------
2233 // Transform to Physical Space (Cartesian Gradient)
2234 // ---------------------------------------------------------
2235 TransformScalarDerivativesToPhysical(aj[k][j][i],
2236 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2237 dGdCsi, dGdEta, dGdZet,
2238 &grad_diff[k][j][i]);
2239 }
2240 }
2241 }
2242
2243 // 7. Restore Arrays
2244 ierr = DMDAVecRestoreArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2245 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2246 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2247 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2248 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2249 ierr = DMDAVecRestoreArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2250
2251 // 8. Update Ghosts for the Result
2252 // Important: Particles near subdomain boundaries will need to interpolate
2253 // this gradient vector, so the ghosts of the vector field must be filled.
2255
2256 PROFILE_FUNCTION_END;
2257 PetscFunctionReturn(0);
2258}
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 using the chain rule.
Definition setup.c:2779
Here is the call graph for this function:
Here is the caller graph for this function:
