PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
 
Loading...
Searching...
No Matches
Macros | Functions
rhs.c File Reference
#include "rhs.h"
#include "verification_sources.h"
#include "les.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)
 Implementation of Convection().
 
PetscErrorCode Viscous (UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
 Implementation of Viscous().
 
PetscErrorCode ComputeBodyForces (UserCtx *user, Vec Rct)
 Internal helper implementation: ComputeBodyForces().
 
PetscErrorCode ComputeRHS (UserCtx *user, Vec Rhs)
 Internal helper implementation: ComputeRHS().
 
PetscErrorCode ComputeEulerianDiffusivity (UserCtx *user)
 Implementation of ComputeEulerianDiffusivity().
 
PetscErrorCode ComputeEulerianDiffusivityGradient (UserCtx *user)
 Internal helper implementation: ComputeEulerianDiffusivityGradient().
 

Macro Definition Documentation

◆ __FUNCT__ [1/6]

#define __FUNCT__   "Convection"

Definition at line 6 of file rhs.c.

◆ __FUNCT__ [2/6]

#define __FUNCT__   "Viscous"

Definition at line 6 of file rhs.c.

◆ __FUNCT__ [3/6]

#define __FUNCT__   "ComputeBodyForces"

Definition at line 6 of file rhs.c.

◆ __FUNCT__ [4/6]

#define __FUNCT__   "ComputeRHS"

Definition at line 6 of file rhs.c.

◆ __FUNCT__ [5/6]

#define __FUNCT__   "ComputeEulerianDiffusivity"

Definition at line 6 of file rhs.c.

◆ __FUNCT__ [6/6]

#define __FUNCT__   "ComputeEulerianDiffusivityGradient"

Definition at line 6 of file rhs.c.

Function Documentation

◆ Convection()

PetscErrorCode Convection ( UserCtx *  user,
Vec  Ucont,
Vec  Ucat,
Vec  Conv 
)

Implementation of Convection().

Computes the convective contribution to the contravariant momentum RHS.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/rhs.h.

See also
Convection()

Definition at line 14 of file rhs.c.

15{
16 PetscErrorCode ierr;
17
18 // --- CONTEXT ACQUISITION BLOCK ---
19 // Get the master simulation context from the UserCtx.
20 SimCtx *simCtx = user->simCtx;
21
22 // Create local variables to mirror the legacy globals for minimal code changes.
23 const LESModelType les = simCtx->les;
24 const PetscInt central = simCtx->central; // Get this from SimCtx now
25 // --- END CONTEXT ACQUISITION BLOCK ---
26
27 Cmpnts ***ucont, ***ucat;
28 DM da = user->da, fda = user->fda;
29 DMDALocalInfo info;
30 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
31 PetscInt mx, my, mz; // Dimensions in three directions
32 PetscInt i, j, k;
33 Vec Fp1, Fp2, Fp3;
34 Cmpnts ***fp1, ***fp2, ***fp3;
35 Cmpnts ***conv;
36
37 PetscReal ucon, up, um;
38 PetscReal coef = 0.125, innerblank=7.;
39
40 PetscInt lxs, lxe, lys, lye, lzs, lze;
41
42 PetscReal ***nvert,***aj;
43
44 PROFILE_FUNCTION_BEGIN;
45
46 DMDAGetLocalInfo(da, &info);
47 mx = info.mx; my = info.my; mz = info.mz;
48 xs = info.xs; xe = xs + info.xm;
49 ys = info.ys; ye = ys + info.ym;
50 zs = info.zs; ze = zs + info.zm;
51 ierr = PreparePeriodicQuickStencilFields(user, Ucat, user->lNvert); CHKERRQ(ierr);
52
53 DMDAVecGetArray(fda, Ucont, &ucont);
54 DMDAVecGetArray(fda, Ucat, &ucat);
55 DMDAVecGetArray(fda, Conv, &conv);
56 DMDAVecGetArray(da, user->lAj, &aj);
57
58 VecDuplicate(Ucont, &Fp1);
59 VecDuplicate(Ucont, &Fp2);
60 VecDuplicate(Ucont, &Fp3);
61
62 DMDAVecGetArray(fda, Fp1, &fp1);
63 DMDAVecGetArray(fda, Fp2, &fp2);
64 DMDAVecGetArray(fda, Fp3, &fp3);
65
66 DMDAVecGetArray(da, user->lNvert, &nvert);
67
68
69 /* We have two different sets of node: 1. grid node, the physical points
70 where grid lines intercross; 2. storage node, where we store variables.
71 All node without explicitly specified as "grid node" refers to
72 storage node.
73
74 The integer node is defined at cell center while half node refers to
75 the actual grid node. (The reason to choose this arrangement is we need
76 ghost node, which is half node away from boundaries, to specify boundary
77 conditions. By using this storage arrangement, the actual storage need
78 is (IM+1) * (JM + 1) * (KM+1) where IM, JM, & KM refer to the number of
79 grid nodes along i, j, k directions.)
80
81 DA, the data structure used to define the storage of 3D arrays, is defined
82 as mx * my * mz. mx = IM+1, my = JM+1, mz = KM+1.
83
84 Staggered grid arrangement is used in this solver.
85 Pressure is stored at interger node (hence the cell center) and volume
86 fluxes defined on the center of each surface of a given control volume
87 is stored on the cloest upper integer node. */
88
89 /* First we calculate the flux on cell surfaces. Stored on the upper integer
90 node. For example, along i direction, the flux are stored at node 0:mx-2*/
91
92 lxs = xs; lxe = xe;
93 lys = ys; lye = ye;
94 lzs = zs; lze = ze;
95
96 if (xs==0) lxs = xs+1;
97 if (ys==0) lys = ys+1;
98 if (zs==0) lzs = zs+1;
99
100 if (xe==mx) lxe=xe-1;
101 if (ye==my) lye=ye-1;
102 if (ze==mz) lze=ze-1;
103
104 VecSet(Conv, 0.0);
105
106 /* Calculating the convective terms on cell centers.
107 First calcualte the contribution from i direction
108 The flux is evaluated by QUICK scheme */
109
110 for (k=lzs; k<lze; k++){
111 for (j=lys; j<lye; j++){
112 for (i=lxs-1; i<lxe; i++){
113
114
115 ucon = ucont[k][j][i].x * 0.5;
116
117 up = ucon + fabs(ucon);
118 um = ucon - fabs(ucon);
119
120 if (i>0 && i<mx-2 &&
121 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
122 (nvert[k][j][i-1] < 0.1 || nvert[k][j][i-1]>innerblank)) { // interial nodes
123 if ((les || central)) {
124 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
125 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
126 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
127
128 } else {
129 fp1[k][j][i].x =
130 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) +
131 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);
132 fp1[k][j][i].y =
133 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) +
134 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);
135 fp1[k][j][i].z =
136 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) +
137 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);
138 }
139 }
140 else if ((les || central) && (i==0 || i==mx-2) &&
141 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
142 (nvert[k][j][i ] < 0.1 || nvert[k][j][i ]>innerblank))
143 {
144 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
145 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
146 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
147 }
148 else if (i==0 ||(nvert[k][j][i-1] > 0.1) ) {
149 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
150 fp1[k][j][i].x =
151 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) +
152 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);
153 fp1[k][j][i].y =
154 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) +
155 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);
156 fp1[k][j][i].z =
157 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) +
158 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);
159 }else{
160 fp1[k][j][i].x =
161 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) +
162 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);
163 fp1[k][j][i].y =
164 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) +
165 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);
166 fp1[k][j][i].z =
167 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) +
168 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);
169 }
170 }
171 else if (i==mx-2 ||(nvert[k][j][i+1]) > 0.1) {
172 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
173 fp1[k][j][i].x =
174 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) +
175 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);
176 fp1[k][j][i].y =
177 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) +
178 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);
179 fp1[k][j][i].z =
180 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) +
181 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);
182 }else{
183 fp1[k][j][i].x =
184 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) +
185 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);
186 fp1[k][j][i].y =
187 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) +
188 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);
189 fp1[k][j][i].z =
190 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) +
191 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);
192 }
193 }
194 }
195 }
196 }
197
198 /* j direction */
199 for (k=lzs; k<lze; k++) {
200 for(j=lys-1; j<lye; j++) {
201 for(i=lxs; i<lxe; i++) {
202 ucon = ucont[k][j][i].y * 0.5;
203
204 up = ucon + fabs(ucon);
205 um = ucon - fabs(ucon);
206
207 if (j>0 && j<my-2 &&
208 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i] > innerblank) &&
209 (nvert[k][j-1][i] < 0.1 || nvert[k][j-1][i] > innerblank)) {
210 if ((les || central)) {
211 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
212 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
213 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
214
215 } else {
216 fp2[k][j][i].x =
217 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) +
218 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);
219 fp2[k][j][i].y =
220 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) +
221 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);
222 fp2[k][j][i].z =
223 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) +
224 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);
225 }
226 }
227 else if ((les || central) && (j==0 || j==my-2) &&
228 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i]>innerblank) &&
229 (nvert[k][j ][i] < 0.1 || nvert[k][j ][i]>innerblank))
230 {
231 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
232 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
233 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
234 }
235 else if (j==0 || (nvert[k][j-1][i]) > 0.1) {
236 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 //
237 fp2[k][j][i].x =
238 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) +
239 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);
240 fp2[k][j][i].y =
241 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) +
242 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);
243 fp2[k][j][i].z =
244 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) +
245 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);
246 }else{
247 fp2[k][j][i].x =
248 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) +
249 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);
250 fp2[k][j][i].y =
251 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) +
252 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);
253 fp2[k][j][i].z =
254 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) +
255 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);
256 }
257 }
258 else if (j==my-2 ||(nvert[k][j+1][i]) > 0.1) {
259 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//
260 fp2[k][j][i].x =
261 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) +
262 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);
263 fp2[k][j][i].y =
264 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) +
265 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);
266 fp2[k][j][i].z =
267 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) +
268 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);
269 }else{
270 fp2[k][j][i].x =
271 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) +
272 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);
273 fp2[k][j][i].y =
274 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) +
275 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);
276 fp2[k][j][i].z =
277 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) +
278 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);
279 }
280 }
281 }
282 }
283 }
284
285
286 /* k direction */
287 for (k=lzs-1; k<lze; k++) {
288 for(j=lys; j<lye; j++) {
289 for(i=lxs; i<lxe; i++) {
290 ucon = ucont[k][j][i].z * 0.5;
291
292 up = ucon + fabs(ucon);
293 um = ucon - fabs(ucon);
294
295 if (k>0 && k<mz-2 &&
296 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i] > innerblank) &&
297 (nvert[k-1][j][i] < 0.1 || nvert[k-1][j][i] > innerblank)) {
298 if ((les || central)) {
299 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
300 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
301 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
302
303 } else {
304 fp3[k][j][i].x =
305 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) +
306 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);
307 fp3[k][j][i].y =
308 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) +
309 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);
310 fp3[k][j][i].z =
311 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) +
312 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);
313 }
314 }
315 else if ((les || central) && (k==0 || k==mz-2) &&
316 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i]>innerblank) &&
317 (nvert[k ][j][i] < 0.1 || nvert[k ][j][i]>innerblank))
318 {
319 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
320 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
321 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
322 }
323 else if (k<mz-2 && (k==0 ||(nvert[k-1][j][i]) > 0.1)) {
324 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//
325 fp3[k][j][i].x =
326 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) +
327 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);
328 fp3[k][j][i].y =
329 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) +
330 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);
331 fp3[k][j][i].z =
332 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) +
333 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);
334 }else{
335 fp3[k][j][i].x =
336 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) +
337 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);
338 fp3[k][j][i].y =
339 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) +
340 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);
341 fp3[k][j][i].z =
342 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) +
343 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);
344 }
345 }
346 else if (k>0 && (k==mz-2 ||(nvert[k+1][j][i]) > 0.1)) {
347 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//
348 fp3[k][j][i].x =
349 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) +
350 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);
351 fp3[k][j][i].y =
352 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) +
353 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);
354 fp3[k][j][i].z =
355 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) +
356 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);
357 }else{
358 fp3[k][j][i].x =
359 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) +
360 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);
361 fp3[k][j][i].y =
362 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) +
363 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);
364 fp3[k][j][i].z =
365 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) +
366 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);
367 }
368 }
369 }
370 }
371 }
372
373 /* Calculate the convective terms under cartesian coordinates */
374
375 for (k=lzs; k<lze; k++) {
376 for (j=lys; j<lye; j++) {
377 for (i=lxs; i<lxe; i++) {
378 conv[k][j][i].x =
379 fp1[k][j][i].x - fp1[k][j][i-1].x +
380 fp2[k][j][i].x - fp2[k][j-1][i].x +
381 fp3[k][j][i].x - fp3[k-1][j][i].x;
382
383 conv[k][j][i].y =
384 fp1[k][j][i].y - fp1[k][j][i-1].y +
385 fp2[k][j][i].y - fp2[k][j-1][i].y +
386 fp3[k][j][i].y - fp3[k-1][j][i].y;
387
388 conv[k][j][i].z =
389 fp1[k][j][i].z - fp1[k][j][i-1].z +
390 fp2[k][j][i].z - fp2[k][j-1][i].z +
391 fp3[k][j][i].z - fp3[k-1][j][i].z;
392 }
393 }
394 }
395 /* for (k=zs; k<ze; k++) { */
396/* for (j=ys; j<ye; j++) { */
397/* for (i=xs; i<xe; i++) { */
398/* if (i==1 && (j==1) && (k==1 || k==21 || k==22|| k==200)) */
399/* 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); */
400/* } */
401/* } */
402/* } */
403
404 DMDAVecRestoreArray(fda, Ucont, &ucont);
405 DMDAVecRestoreArray(fda, Ucat, &ucat);
406 DMDAVecRestoreArray(fda, Conv, &conv);
407 DMDAVecRestoreArray(da, user->lAj, &aj);
408
409 DMDAVecRestoreArray(fda, Fp1, &fp1);
410 DMDAVecRestoreArray(fda, Fp2, &fp2);
411 DMDAVecRestoreArray(fda, Fp3, &fp3);
412 DMDAVecRestoreArray(da, user->lNvert, &nvert);
413
414 VecDestroy(&Fp1);
415 VecDestroy(&Fp2);
416 VecDestroy(&Fp3);
417
418
419 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Convective term calculated .\n");
420
421 PROFILE_FUNCTION_END;
422 return (0);
423}
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:2060
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:894
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:885
LESModelType
LESModelType
Identifies the subgrid-scale closure evaluated during a timestep.
Definition variables.h:548
PERIODIC
@ PERIODIC
Definition variables.h:318
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:1099
UserCtx::lNvert
Vec lNvert
Definition variables.h:1113
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:1077
Cmpnts::x
PetscScalar x
Definition variables.h:122
Cmpnts::z
PetscScalar z
Definition variables.h:122
SimCtx::central
PetscInt central
Definition variables.h:904
Cmpnts::y
PetscScalar y
Definition variables.h:122
SimCtx::les
PetscInt les
Active LES closure; an LESModelType value.
Definition variables.h:984
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:392
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:288
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:290
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:289
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:121
SimCtx
The master context for the entire simulation.
Definition variables.h:859
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Viscous()

PetscErrorCode Viscous ( UserCtx *  user,
Vec  Ucont,
Vec  Ucat,
Vec  Visc 
)

Implementation of Viscous().

Computes the viscous contribution to the contravariant momentum RHS.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/rhs.h.

See also
Viscous()

Definition at line 435 of file rhs.c.

436{
437
438 Vec Csi = user->lCsi, Eta = user->lEta, Zet = user->lZet;
439
440 Cmpnts ***ucont, ***ucat;
441
442 Cmpnts ***csi, ***eta, ***zet;
443 Cmpnts ***icsi, ***ieta, ***izet;
444 Cmpnts ***jcsi, ***jeta, ***jzet;
445 Cmpnts ***kcsi, ***keta, ***kzet;
446
447 PetscReal ***nvert;
448
449 DM da = user->da, fda = user->fda;
450 DMDALocalInfo info;
451 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
452 PetscInt mx, my, mz; // Dimensions in three directions
453 PetscInt i, j, k;
454 Vec Fp1, Fp2, Fp3;
455 Cmpnts ***fp1, ***fp2, ***fp3;
456 Cmpnts ***visc;
457 PetscReal ***aj, ***iaj, ***jaj, ***kaj;
458
459 PetscInt lxs, lxe, lys, lye, lzs, lze;
460
461 PetscReal ajc;
462
463 PetscReal dudc, dude, dudz, dvdc, dvde, dvdz, dwdc, dwde, dwdz;
464 PetscReal csi0, csi1, csi2, eta0, eta1, eta2, zet0, zet1, zet2;
465 PetscReal g11, g21, g31;
466 PetscReal r11, r21, r31, r12, r22, r32, r13, r23, r33;
467
468 PetscScalar solid,innerblank;
469
470 // --- CONTEXT ACQUISITION BLOCK ---
471 // Get the master simulation context from the UserCtx.
472 SimCtx *simCtx = user->simCtx;
473
474 // Create local variables to mirror the legacy globals for minimal code changes.
475 const LESModelType les = simCtx->les;
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->les_gradient_model;
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 PetscReal ***lnu_wall = NULL;
543
544 if(les) {
545 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
546 }
547 /* A wall model replaces the subgrid viscosity at its own face rather than adding to
548 it, so the flux there carries the stress the model computed instead of the
549 molecular fraction of it. */
550 if (user->simCtx->wallfunction) {
551 DMDAVecGetArray(da, user->lNu_Wall, &lnu_wall);
552 }
553
554 /* The visc flux on each surface center is stored at previous integer node */
555
556 DMDAVecGetArray(da, user->lIAj, &iaj);
557 /* for (k=zs; k<ze; k++) { */
558/* for (j=ys; j<ye; j++) { */
559/* for (i=xs; i<xe; i++) { */
560/* if (i==1 && (j==0 ||j==1 || j==2) && (k==21 || k==22|| k==20)) */
561/* 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 ); */
562/* } */
563/* } */
564/* } */
565 // i direction
566 for (k=lzs; k<lze; k++) {
567 for (j=lys; j<lye; j++) {
568 for (i=lxs-1; i<lxe; i++) {
569
570 dudc = ucat[k][j][i+1].x - ucat[k][j][i].x;
571 dvdc = ucat[k][j][i+1].y - ucat[k][j][i].y;
572 dwdc = ucat[k][j][i+1].z - ucat[k][j][i].z;
573
574 if ((nvert[k][j+1][i ]> solid && nvert[k][j+1][i ]<innerblank) ||
575 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
576 dude = (ucat[k][j ][i+1].x + ucat[k][j ][i].x -
577 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.5;
578 dvde = (ucat[k][j ][i+1].y + ucat[k][j ][i].y -
579 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.5;
580 dwde = (ucat[k][j ][i+1].z + ucat[k][j ][i].z -
581 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.5;
582 }
583 else if ((nvert[k][j-1][i ]> solid && nvert[k][j-1][i ]<innerblank) ||
584 (nvert[k][j-1][i+1]> solid && nvert[k][j-1][i+1]<innerblank)) {
585 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
586 ucat[k][j ][i+1].x - ucat[k][j ][i].x) * 0.5;
587 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
588 ucat[k][j ][i+1].y - ucat[k][j ][i].y) * 0.5;
589 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
590 ucat[k][j ][i+1].z - ucat[k][j ][i].z) * 0.5;
591 }
592 else {
593 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
594 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.25;
595 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
596 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.25;
597 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
598 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.25;
599 }
600
601 if ((nvert[k+1][j][i ]> solid && nvert[k+1][j][i ]<innerblank)||
602 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
603 dudz = (ucat[k ][j][i+1].x + ucat[k ][j][i].x -
604 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.5;
605 dvdz = (ucat[k ][j][i+1].y + ucat[k ][j][i].y -
606 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.5;
607 dwdz = (ucat[k ][j][i+1].z + ucat[k ][j][i].z -
608 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.5;
609 }
610 else if ((nvert[k-1][j][i ]> solid && nvert[k-1][j][i ]<innerblank) ||
611 (nvert[k-1][j][i+1]> solid && nvert[k-1][j][i+1]<innerblank)) {
612
613 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
614 ucat[k ][j][i+1].x - ucat[k ][j][i].x) * 0.5;
615 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
616 ucat[k ][j][i+1].y - ucat[k ][j][i].y) * 0.5;
617 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
618 ucat[k ][j][i+1].z - ucat[k ][j][i].z) * 0.5;
619 }
620 else {
621 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
622 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.25;
623 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
624 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.25;
625 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
626 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.25;
627 }
628
629 csi0 = icsi[k][j][i].x;
630 csi1 = icsi[k][j][i].y;
631 csi2 = icsi[k][j][i].z;
632
633 eta0 = ieta[k][j][i].x;
634 eta1 = ieta[k][j][i].y;
635 eta2 = ieta[k][j][i].z;
636
637 zet0 = izet[k][j][i].x;
638 zet1 = izet[k][j][i].y;
639 zet2 = izet[k][j][i].z;
640
641 g11 = csi0 * csi0 + csi1 * csi1 + csi2 * csi2;
642 g21 = eta0 * csi0 + eta1 * csi1 + eta2 * csi2;
643 g31 = zet0 * csi0 + zet1 * csi1 + zet2 * csi2;
644
645 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
646 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
647 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
648
649 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
650 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
651 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
652
653 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
654 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
655 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
656
657 ajc = iaj[k][j][i];
658
659 double nu = 1./ren, nu_t=0;
660
661 if( les ) {
662 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j][i+1]) ), 2.0) * Sabs;
663 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j][i+1]);
664 /* Zero is right for a wall-resolved run, where the subgrid stress vanishes at
665 the wall. With a wall model the same face is where the modelled stress has to
666 enter, so it takes the model's effective viscosity instead. */
667 if ( user->boundary_faces[BC_FACE_NEG_X].mathematical_type == WALL && i==0 ) nu_t = lnu_wall ? lnu_wall[k][j][1] : 0.0;
668 if ( user->boundary_faces[BC_FACE_POS_X].mathematical_type == WALL && i==mx-2 ) nu_t = lnu_wall ? lnu_wall[k][j][mx-2] : 0.0;
669 fp1[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * csi0 + r21 * csi1 + r31 * csi2) * ajc * (nu_t);
670 fp1[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * csi0 + r22 * csi1 + r32 * csi2) * ajc * (nu_t);
671 fp1[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * csi0 + r23 * csi1 + r33 * csi2) * ajc * (nu_t);
672 }
673 else {
674 fp1[k][j][i].x = 0;
675 fp1[k][j][i].y = 0;
676 fp1[k][j][i].z = 0;
677 }
678
679 fp1[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * csi0 + r21 * csi1 + r31 * csi2 ) * ajc * (nu);
680 fp1[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * csi0 + r22 * csi1 + r32 * csi2 ) * ajc * (nu);
681 fp1[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * csi0 + r23 * csi1 + r33 * csi2 ) * ajc * (nu);
682
683
684 if(clark) {
685 /* dudc, dude, dudz are differences between neighbouring cells along each grid
686 direction, so each already equals that direction's spacing times the physical
687 derivative. Their outer products are therefore Delta_k^2 du_i/dx_k du_j/dx_k
688 as they stand. Multiplying by a squared cell extent as well - which this term
689 did, carried over from the legacy - made the stress scale as Delta^4 and left
690 it smaller than intended by Delta^2, 1e-4 to 1e-6 on a production mesh; the
691 extents were also Cartesian, so they did not belong to the directions they
692 were paired with once the grid turned. */
693
694 /* The gradient-model tensor is the sum over the three computational
695 directions of the per-cell velocity difference's outer product with itself;
696 the cell's own extent in each direction is already inside that difference.
697 It is symmetric by construction, so it is carried as one. */
698 const Cmpnts grad_csi = { dudc, dvdc, dwdc };
699 const Cmpnts grad_eta = { dude, dvde, dwde };
700 const Cmpnts grad_zet = { dudz, dvdz, dwdz };
701 const SymTensor gradient_tensor =
702 SymTensorCombine(1.0, SymTensorSelfOuter(grad_csi), 1.0,
703 SymTensorCombine(1.0, SymTensorSelfOuter(grad_eta),
704 1.0, SymTensorSelfOuter(grad_zet)));
705
706 {
707 const Cmpnts face_normal = { csi0, csi1, csi2 };
708 const Cmpnts gradient_flux = SymTensorTimesVector(gradient_tensor, face_normal);
709
710 fp1[k][j][i].x -= gradient_flux.x / 12.;
711 fp1[k][j][i].y -= gradient_flux.y / 12.;
712 fp1[k][j][i].z -= gradient_flux.z / 12.;
713 }
714 }
715
716 }
717 }
718 }
719 DMDAVecRestoreArray(da, user->lIAj, &iaj);
720
721
722 // j direction
723 DMDAVecGetArray(da, user->lJAj, &jaj);
724 for (k=lzs; k<lze; k++) {
725 for (j=lys-1; j<lye; j++) {
726 for (i=lxs; i<lxe; i++) {
727
728 if ((nvert[k][j ][i+1]> solid && nvert[k][j ][i+1]<innerblank)||
729 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
730 dudc = (ucat[k][j+1][i ].x + ucat[k][j][i ].x -
731 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.5;
732 dvdc = (ucat[k][j+1][i ].y + ucat[k][j][i ].y -
733 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.5;
734 dwdc = (ucat[k][j+1][i ].z + ucat[k][j][i ].z -
735 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.5;
736 }
737 else if ((nvert[k][j ][i-1]> solid && nvert[k][j ][i-1]<innerblank) ||
738 (nvert[k][j+1][i-1]> solid && nvert[k][j+1][i-1]<innerblank)) {
739 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
740 ucat[k][j+1][i ].x - ucat[k][j][i ].x) * 0.5;
741 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
742 ucat[k][j+1][i ].y - ucat[k][j][i ].y) * 0.5;
743 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
744 ucat[k][j+1][i ].z - ucat[k][j][i ].z) * 0.5;
745 }
746 else {
747 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
748 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.25;
749 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
750 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.25;
751 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
752 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.25;
753 }
754
755 dude = ucat[k][j+1][i].x - ucat[k][j][i].x;
756 dvde = ucat[k][j+1][i].y - ucat[k][j][i].y;
757 dwde = ucat[k][j+1][i].z - ucat[k][j][i].z;
758
759 if ((nvert[k+1][j ][i]> solid && nvert[k+1][j ][i]<innerblank)||
760 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
761 dudz = (ucat[k ][j+1][i].x + ucat[k ][j][i].x -
762 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.5;
763 dvdz = (ucat[k ][j+1][i].y + ucat[k ][j][i].y -
764 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.5;
765 dwdz = (ucat[k ][j+1][i].z + ucat[k ][j][i].z -
766 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.5;
767 }
768 else if ((nvert[k-1][j ][i]> solid && nvert[k-1][j ][i]<innerblank)||
769 (nvert[k-1][j+1][i]> solid && nvert[k-1][j+1][i]<innerblank)) {
770 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
771 ucat[k ][j+1][i].x - ucat[k ][j][i].x) * 0.5;
772 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
773 ucat[k ][j+1][i].y - ucat[k ][j][i].y) * 0.5;
774 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
775 ucat[k ][j+1][i].z - ucat[k ][j][i].z) * 0.5;
776 }
777 else {
778 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
779 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.25;
780 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
781 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.25;
782 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
783 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.25;
784 }
785
786 csi0 = jcsi[k][j][i].x;
787 csi1 = jcsi[k][j][i].y;
788 csi2 = jcsi[k][j][i].z;
789
790 eta0 = jeta[k][j][i].x;
791 eta1 = jeta[k][j][i].y;
792 eta2 = jeta[k][j][i].z;
793
794 zet0 = jzet[k][j][i].x;
795 zet1 = jzet[k][j][i].y;
796 zet2 = jzet[k][j][i].z;
797
798
799 g11 = csi0 * eta0 + csi1 * eta1 + csi2 * eta2;
800 g21 = eta0 * eta0 + eta1 * eta1 + eta2 * eta2;
801 g31 = zet0 * eta0 + zet1 * eta1 + zet2 * eta2;
802
803 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
804 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
805 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
806
807 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
808 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
809 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
810
811 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
812 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
813 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
814
815 // 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);
816 // 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);
817 // 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);
818 // 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);
819
820
821
822 ajc = jaj[k][j][i];
823
824 double nu = 1./ren, nu_t = 0;
825
826 if( les ) {
827 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j+1][i]) ), 2.0) * Sabs;
828 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j+1][i]);
829 /* Zero is right for a wall-resolved run, where the subgrid stress vanishes at
830 the wall. With a wall model the same face is where the modelled stress has to
831 enter, so it takes the model's effective viscosity instead. */
832 if ( user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == WALL && j==0 ) nu_t = lnu_wall ? lnu_wall[k][1][i] : 0.0;
833 if ( user->boundary_faces[BC_FACE_POS_Y].mathematical_type == WALL && j==my-2 ) nu_t = lnu_wall ? lnu_wall[k][my-2][i] : 0.0;
834
835 fp2[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * eta0 + r21 * eta1 + r31 * eta2) * ajc * (nu_t);
836 fp2[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * eta0 + r22 * eta1 + r32 * eta2) * ajc * (nu_t);
837 fp2[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * eta0 + r23 * eta1 + r33 * eta2) * ajc * (nu_t);
838 }
839 else {
840 fp2[k][j][i].x = 0;
841 fp2[k][j][i].y = 0;
842 fp2[k][j][i].z = 0;
843 }
844
845 fp2[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * eta0 + r21 * eta1 + r31 * eta2 ) * ajc * (nu);
846 fp2[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * eta0 + r22 * eta1 + r32 * eta2 ) * ajc * (nu);
847 fp2[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * eta0 + r23 * eta1 + r33 * eta2 ) * ajc * (nu);
848
849 if(clark) {
850 /* dudc, dude, dudz are differences between neighbouring cells along each grid
851 direction, so each already equals that direction's spacing times the physical
852 derivative. Their outer products are therefore Delta_k^2 du_i/dx_k du_j/dx_k
853 as they stand. Multiplying by a squared cell extent as well - which this term
854 did, carried over from the legacy - made the stress scale as Delta^4 and left
855 it smaller than intended by Delta^2, 1e-4 to 1e-6 on a production mesh; the
856 extents were also Cartesian, so they did not belong to the directions they
857 were paired with once the grid turned. */
858
859 /* The gradient-model tensor is the sum over the three computational
860 directions of the per-cell velocity difference's outer product with itself;
861 the cell's own extent in each direction is already inside that difference.
862 It is symmetric by construction, so it is carried as one. */
863 const Cmpnts grad_csi = { dudc, dvdc, dwdc };
864 const Cmpnts grad_eta = { dude, dvde, dwde };
865 const Cmpnts grad_zet = { dudz, dvdz, dwdz };
866 const SymTensor gradient_tensor =
867 SymTensorCombine(1.0, SymTensorSelfOuter(grad_csi), 1.0,
868 SymTensorCombine(1.0, SymTensorSelfOuter(grad_eta),
869 1.0, SymTensorSelfOuter(grad_zet)));
870
871 {
872 const Cmpnts face_normal = { eta0, eta1, eta2 };
873 const Cmpnts gradient_flux = SymTensorTimesVector(gradient_tensor, face_normal);
874
875 fp2[k][j][i].x -= gradient_flux.x / 12.;
876 fp2[k][j][i].y -= gradient_flux.y / 12.;
877 fp2[k][j][i].z -= gradient_flux.z / 12.;
878 }
879 }
880 }
881 }
882 }
883
884 DMDAVecRestoreArray(da, user->lJAj, &jaj);
885 // k direction
886
887 DMDAVecGetArray(da, user->lKAj, &kaj);
888 for (k=lzs-1; k<lze; k++) {
889 for (j=lys; j<lye; j++) {
890 for (i=lxs; i<lxe; i++) {
891 if ((nvert[k ][j][i+1]> solid && nvert[k ][j][i+1]<innerblank)||
892 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
893 dudc = (ucat[k+1][j][i ].x + ucat[k][j][i ].x -
894 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.5;
895 dvdc = (ucat[k+1][j][i ].y + ucat[k][j][i ].y -
896 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.5;
897 dwdc = (ucat[k+1][j][i ].z + ucat[k][j][i ].z -
898 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.5;
899 }
900 else if ((nvert[k ][j][i-1]> solid && nvert[k ][j][i-1]<innerblank) ||
901 (nvert[k+1][j][i-1]> solid && nvert[k+1][j][i-1]<innerblank)) {
902 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
903 ucat[k+1][j][i ].x - ucat[k][j][i ].x) * 0.5;
904 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
905 ucat[k+1][j][i ].y - ucat[k][j][i ].y) * 0.5;
906 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
907 ucat[k+1][j][i ].z - ucat[k][j][i ].z) * 0.5;
908 }
909 else {
910 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
911 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.25;
912 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
913 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.25;
914 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
915 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.25;
916 }
917
918 if ((nvert[k ][j+1][i]> solid && nvert[k ][j+1][i]<innerblank)||
919 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
920 dude = (ucat[k+1][j ][i].x + ucat[k][j ][i].x -
921 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.5;
922 dvde = (ucat[k+1][j ][i].y + ucat[k][j ][i].y -
923 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.5;
924 dwde = (ucat[k+1][j ][i].z + ucat[k][j ][i].z -
925 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.5;
926 }
927 else if ((nvert[k ][j-1][i]> solid && nvert[k ][j-1][i]<innerblank) ||
928 (nvert[k+1][j-1][i]> solid && nvert[k+1][j-1][i]<innerblank)){
929 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
930 ucat[k+1][j ][i].x - ucat[k][j ][i].x) * 0.5;
931 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
932 ucat[k+1][j ][i].y - ucat[k][j ][i].y) * 0.5;
933 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
934 ucat[k+1][j ][i].z - ucat[k][j ][i].z) * 0.5;
935 }
936 else {
937 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
938 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.25;
939 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
940 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.25;
941 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
942 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.25;
943 }
944
945 dudz = ucat[k+1][j][i].x - ucat[k][j][i].x;
946 dvdz = ucat[k+1][j][i].y - ucat[k][j][i].y;
947 dwdz = ucat[k+1][j][i].z - ucat[k][j][i].z;
948
949
950 csi0 = kcsi[k][j][i].x;
951 csi1 = kcsi[k][j][i].y;
952 csi2 = kcsi[k][j][i].z;
953
954 eta0 = keta[k][j][i].x;
955 eta1 = keta[k][j][i].y;
956 eta2 = keta[k][j][i].z;
957
958 zet0 = kzet[k][j][i].x;
959 zet1 = kzet[k][j][i].y;
960 zet2 = kzet[k][j][i].z;
961
962
963 g11 = csi0 * zet0 + csi1 * zet1 + csi2 * zet2;
964 g21 = eta0 * zet0 + eta1 * zet1 + eta2 * zet2;
965 g31 = zet0 * zet0 + zet1 * zet1 + zet2 * zet2;
966
967 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
968 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
969 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
970
971 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
972 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
973 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
974
975 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
976 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
977 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
978
979 ajc = kaj[k][j][i];
980
981 double nu = 1./ren, nu_t =0;
982
983 if( les ) {
984 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k+1][j][i]) ), 2.0) * Sabs;
985 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k+1][j][i]);
986 /* Zero is right for a wall-resolved run, where the subgrid stress vanishes at
987 the wall. With a wall model the same face is where the modelled stress has to
988 enter, so it takes the model's effective viscosity instead. */
989 if ( user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == WALL && k==0 ) nu_t = lnu_wall ? lnu_wall[1][j][i] : 0.0;
990 if ( user->boundary_faces[BC_FACE_POS_Z].mathematical_type == WALL && k==mz-2 ) nu_t = lnu_wall ? lnu_wall[mz-2][j][i] : 0.0;
991
992 fp3[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu_t);
993 fp3[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu_t);
994 fp3[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu_t);
995 }
996 else {
997 fp3[k][j][i].x = 0;
998 fp3[k][j][i].y = 0;
999 fp3[k][j][i].z = 0;
1000 }
1001 fp3[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu);//
1002 fp3[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu);//
1003 fp3[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu);//
1004
1005 if(clark) {
1006 /* dudc, dude, dudz are differences between neighbouring cells along each grid
1007 direction, so each already equals that direction's spacing times the physical
1008 derivative. Their outer products are therefore Delta_k^2 du_i/dx_k du_j/dx_k
1009 as they stand. Multiplying by a squared cell extent as well - which this term
1010 did, carried over from the legacy - made the stress scale as Delta^4 and left
1011 it smaller than intended by Delta^2, 1e-4 to 1e-6 on a production mesh; the
1012 extents were also Cartesian, so they did not belong to the directions they
1013 were paired with once the grid turned. */
1014
1015 /* The gradient-model tensor is the sum over the three computational
1016 directions of the per-cell velocity difference's outer product with itself;
1017 the cell's own extent in each direction is already inside that difference.
1018 It is symmetric by construction, so it is carried as one. */
1019 const Cmpnts grad_csi = { dudc, dvdc, dwdc };
1020 const Cmpnts grad_eta = { dude, dvde, dwde };
1021 const Cmpnts grad_zet = { dudz, dvdz, dwdz };
1022 const SymTensor gradient_tensor =
1023 SymTensorCombine(1.0, SymTensorSelfOuter(grad_csi), 1.0,
1024 SymTensorCombine(1.0, SymTensorSelfOuter(grad_eta),
1025 1.0, SymTensorSelfOuter(grad_zet)));
1026
1027 {
1028 const Cmpnts face_normal = { zet0, zet1, zet2 };
1029 const Cmpnts gradient_flux = SymTensorTimesVector(gradient_tensor, face_normal);
1030
1031 fp3[k][j][i].x -= gradient_flux.x / 12.;
1032 fp3[k][j][i].y -= gradient_flux.y / 12.;
1033 fp3[k][j][i].z -= gradient_flux.z / 12.;
1034 }
1035 }
1036 }
1037 }
1038 }
1039
1040 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1041
1042 for (k=lzs; k<lze; k++) {
1043 for (j=lys; j<lye; j++) {
1044 for (i=lxs; i<lxe; i++) {
1045 visc[k][j][i].x =
1046 (fp1[k][j][i].x - fp1[k][j][i-1].x +
1047 fp2[k][j][i].x - fp2[k][j-1][i].x +
1048 fp3[k][j][i].x - fp3[k-1][j][i].x);
1049
1050 visc[k][j][i].y =
1051 (fp1[k][j][i].y - fp1[k][j][i-1].y +
1052 fp2[k][j][i].y - fp2[k][j-1][i].y +
1053 fp3[k][j][i].y - fp3[k-1][j][i].y);
1054
1055 visc[k][j][i].z =
1056 (fp1[k][j][i].z - fp1[k][j][i-1].z +
1057 fp2[k][j][i].z - fp2[k][j-1][i].z +
1058 fp3[k][j][i].z - fp3[k-1][j][i].z);
1059
1060 }
1061 }
1062 }
1063/* for (k=zs; k<ze; k++) { */
1064/* for (j=ys; j<ye; j++) { */
1065/* for (i=xs; i<xe; i++) { */
1066/* 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); */
1067/* 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); */
1068/* 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); */
1069/* 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); */
1070/* 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); */
1071/* 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); */
1072
1073/* } */
1074/* } */
1075/* } */
1076 DMDAVecRestoreArray(fda, Ucont, &ucont);
1077 DMDAVecRestoreArray(fda, Ucat, &ucat);
1078 DMDAVecRestoreArray(fda, Visc, &visc);
1079
1080 DMDAVecRestoreArray(fda, Csi, &csi);
1081 DMDAVecRestoreArray(fda, Eta, &eta);
1082 DMDAVecRestoreArray(fda, Zet, &zet);
1083
1084 DMDAVecRestoreArray(fda, Fp1, &fp1);
1085 DMDAVecRestoreArray(fda, Fp2, &fp2);
1086 DMDAVecRestoreArray(fda, Fp3, &fp3);
1087
1088 DMDAVecRestoreArray(da, user->lAj, &aj);
1089
1090 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1091 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1092 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1093
1094 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1095 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1096 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1097
1098 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1099 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1100 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1101
1102 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1103
1104 if(les) {
1105 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1106 }
1107 /* Opened independently of the turbulence model, so closed independently of it. */
1108 if (lnu_wall) DMDAVecRestoreArray(da, user->lNu_Wall, &lnu_wall);
1109
1110
1111 VecDestroy(&Fp1);
1112 VecDestroy(&Fp2);
1113 VecDestroy(&Fp3);
1114
1115
1116 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Viscous terms calculated .\n");
1117
1118 PROFILE_FUNCTION_END;
1119
1120 return(0);
1121}
SymTensorTimesVector
Cmpnts SymTensorTimesVector(SymTensor t, Cmpnts v)
Applies a symmetric tensor to a vector, returning t_ij v_j.
Definition les.c:232
SymTensorCombine
SymTensor SymTensorCombine(PetscReal a, SymTensor x, PetscReal b, SymTensor y)
Forms the linear combination a*x + b*y.
Definition les.c:165
SymTensorSelfOuter
SymTensor SymTensorSelfOuter(Cmpnts v)
Forms a symmetric tensor from a vector's outer product with itself.
Definition les.c:144
WALL
@ WALL
Definition variables.h:312
UserCtx::lNu_Wall
Vec lNu_Wall
Definition variables.h:1110
UserCtx::lIEta
Vec lIEta
Definition variables.h:1151
UserCtx::lIZet
Vec lIZet
Definition variables.h:1151
SimCtx::ren
PetscReal ren
Definition variables.h:906
UserCtx::lKEta
Vec lKEta
Definition variables.h:1153
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:1152
UserCtx::lKZet
Vec lKZet
Definition variables.h:1153
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:1156
UserCtx::lJEta
Vec lJEta
Definition variables.h:1152
SimCtx::les_gradient_model
PetscInt les_gradient_model
Add the Clark gradient (tensor-diffusivity) term to the viscous flux.
Definition variables.h:986
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:1153
SimCtx::wallfunction
PetscInt wallfunction
Enable wall functions on WALL faces.
Definition variables.h:985
UserCtx::lJZet
Vec lJZet
Definition variables.h:1152
SimCtx::step
PetscInt step
Definition variables.h:867
UserCtx::lICsi
Vec lICsi
Definition variables.h:1151
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:290
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:289
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:288
SymTensor
A symmetric second-order tensor stored by its six independent components.
Definition variables.h:138
nu_t
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:303
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeBodyForces()

PetscErrorCode ComputeBodyForces ( UserCtx *  user,
Vec  Rct 
)

Internal helper implementation: ComputeBodyForces().

General dispatcher for applying all active body forces (momentum sources).

Local to this translation unit.

Definition at line 1129 of file rhs.c.

1130{
1131 PetscErrorCode ierr;
1132 PetscFunctionBeginUser;
1133
1134 // --- 1. Apply momentum source for driven channel/pipe flows ---
1135 // This function will internally check if a driven flow BC is active.
1136 ierr = ComputeDrivenChannelFlowSource(user, Rct); CHKERRQ(ierr);
1137
1138 // --- 2. (Future Extension) Apply gravitational force ---
1139 // if (user->simCtx->gravityEnabled) {
1140 // ierr = ApplyGravitationalForce(user, Rhs); CHKERRQ(ierr);
1141 // }
1142 //
1143 // Adding a body force here: see the contract in include/BodyForces.h.
1144 // In particular, this function runs once per RESIDUAL EVALUATION, not once
1145 // per timestep, so any force carrying state across calls (a filter, ramp,
1146 // moving average, or integral term) must gate its update on simCtx->step.
1147
1148 PetscFunctionReturn(0);
1149}
ComputeDrivenChannelFlowSource
PetscErrorCode ComputeDrivenChannelFlowSource(UserCtx *user, Vec Rct)
Applies a momentum source term to drive flow in a periodic channel or pipe.
Definition BodyForces.c:37
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeRHS()

PetscErrorCode ComputeRHS ( UserCtx *  user,
Vec  Rhs 
)

Internal helper implementation: ComputeRHS().

Computes the Right-Hand Side (RHS) of the momentum equations.

Local to this translation unit.

Definition at line 1157 of file rhs.c.

1158{
1159 PetscErrorCode ierr;
1160 SimCtx *simCtx = user->simCtx;
1161 DM da = user->da, fda = user->fda;
1162 DMDALocalInfo info = user->info;
1163 PetscInt i,j,k;
1164 // --- Local Grid Indices and Parameters ---
1165 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1166 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1167 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1168 PetscInt lxs = (xs==0) ? xs+1 : xs;
1169 PetscInt lys = (ys==0) ? ys+1 : ys;
1170 PetscInt lzs = (zs==0) ? zs+1 : zs;
1171 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1172 PetscInt lye = (ye==my) ? ye-1 : ye;
1173 PetscInt lze = (ze==mz) ? ze-1 : ze;
1174
1175 // --- Array Pointers ---
1176 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1177 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert;
1178 Cmpnts ***rhs, ***rc, ***rct;
1179
1180 // --- Temporary Vectors ---
1181 Vec Conv, Visc, Rc, Rct;
1182
1183 PetscFunctionBeginUser;
1184 PROFILE_FUNCTION_BEGIN;
1185 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1186 simCtx->rank, user->_this);
1187
1188 // --- Get all necessary array pointers ---
1189 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1190 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1191 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1192 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1193 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1194 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1195 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1196 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1197 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1198 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1199 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1200 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1201 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1202 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1203 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1204 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1205 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1206 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1207 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1208
1209 // --- Create temporary work vectors ---
1210 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1211 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1212 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1213 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1214
1215 // ========================================================================
1216 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1217 // ========================================================================
1218
1219 // 1. Obtain Cartesian velocity from Contravariant velocity
1220 ierr = Contra2Cart(user); CHKERRQ(ierr);
1221 {
1222 const FieldId cell_fields[] = {FIELD_ID_UCAT};
1223 ierr = SynchronizePeriodicCellFields(user, 1, cell_fields); CHKERRQ(ierr);
1224 }
1225 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
1226
1227 // 2. Compute Convective term
1228 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1229 if (simCtx->moveframe || simCtx->rotateframe) {
1230 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1231 } else {
1232 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1233 }
1234
1235 // 3. Compute Viscous term
1236 if (simCtx->invicid) {
1237 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1238 } else {
1239 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1240 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1241 }
1242
1243 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1244 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1245
1246 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1247 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1248 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1249 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1250
1251 for (k = lzs; k < lze; k++) {
1252 for (j = lys; j < lye; j++) {
1253 for (i = lxs; i < lxe; i++) {
1254 rct[k][j][i].x = aj[k][j][i] *
1255 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1256 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1257 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1258 rct[k][j][i].y = aj[k][j][i] *
1259 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1260 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1261 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1262 rct[k][j][i].z = aj[k][j][i] *
1263 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1264 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1265 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1266 }
1267 }
1268 }
1269 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1270 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1271
1272 PetscBarrier(NULL);
1273
1274 // Compute and Add Body Force term if applicable.
1275 ierr = ComputeBodyForces(user,Rct); CHKERRQ(ierr);
1276 ierr = SynchronizePeriodicLocalStaggeredField(user, Rct); CHKERRQ(ierr);
1277
1278 // 6. Add Pressure Gradient Term and Finalize RHS
1279 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1280 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1281 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1282
1283 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1284
1285 for (k = lzs; k < lze; k++) {
1286 for (j = lys; j < lye; j++) {
1287 for (i = lxs; i < lxe; i++) {
1288 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1289 dpdc = p[k][j][i+1] - p[k][j][i];
1290
1291 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) {
1292 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))) {
1293 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1294 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1295 }
1296 }
1297 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) {
1298 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1299 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1300 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1301 }
1302 }
1303 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) {
1304 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1305 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1306 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1307 }
1308 }
1309 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) {
1310 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1311 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1312 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1313 }
1314 }
1315 else {
1316 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1317 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1318 }
1319
1320 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) {
1321 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))) {
1322 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1323 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1324 }
1325 }
1326 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) {
1327 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1328 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1329 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1330 }
1331 }
1332 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) {
1333 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1334 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1335 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1336 }
1337 }
1338 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) {
1339 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1340 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1341 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1342 }
1343 }
1344 else {
1345 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1346 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1347 }
1348
1349 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1350
1351
1352 rhs[k][j][i].x -=
1353 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1354 icsi[k][j][i].y * icsi[k][j][i].y +
1355 icsi[k][j][i].z * icsi[k][j][i].z)+
1356 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1357 ieta[k][j][i].y * icsi[k][j][i].y +
1358 ieta[k][j][i].z * icsi[k][j][i].z)+
1359 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1360 izet[k][j][i].y * icsi[k][j][i].y +
1361 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1362
1363 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) {
1364 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))) {
1365 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1366 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1367 }
1368 }
1369 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) {
1370 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1371 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1372 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1373 }
1374 }
1375 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) {
1376 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1377 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1378 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1379 }
1380 }
1381 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) {
1382 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1383 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1384 p[k][j][i] - p[k][j+1][i]) * 0.5;
1385 }
1386 }
1387 else {
1388 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1389 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1390 }
1391
1392 dpde = p[k][j+1][i] - p[k][j][i];
1393
1394 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) {
1395 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))) {
1396 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1397 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1398 }
1399 }
1400 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) {
1401 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1402 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1403 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1404 }
1405 }
1406 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) {
1407 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1408 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1409 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1410 }
1411 }
1412 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) {
1413 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1414 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1415 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1416 }
1417 }
1418 else {
1419 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1420 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1421 }
1422
1423 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1424
1425
1426 rhs[k][j][i].y -=
1427 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1428 jcsi[k][j][i].y * jeta[k][j][i].y +
1429 jcsi[k][j][i].z * jeta[k][j][i].z) +
1430 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1431 jeta[k][j][i].y * jeta[k][j][i].y +
1432 jeta[k][j][i].z * jeta[k][j][i].z) +
1433 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1434 jzet[k][j][i].y * jeta[k][j][i].y +
1435 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1436
1437 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) {
1438 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))) {
1439 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1440 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1441 }
1442 }
1443 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) {
1444 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1445 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1446 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1447 }
1448 }
1449 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) {
1450 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1451 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1452 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1453 }
1454 }
1455 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) {
1456 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1457 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1458 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1459 }
1460 }
1461 else {
1462 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1463 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1464 }
1465
1466 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) {
1467 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))) {
1468 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1469 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1470 }
1471 }
1472 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) {
1473 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1474 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1475 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1476 }
1477 }
1478 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) {
1479 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1480 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1481 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1482 }
1483 }
1484 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) {
1485 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1486 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1487 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1488 }
1489 }
1490 else {
1491 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1492 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1493 }
1494
1495 dpdz = (p[k+1][j][i] - p[k][j][i]);
1496
1497 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1498
1499 rhs[k][j][i].z -=
1500 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1501 kcsi[k][j][i].y * kzet[k][j][i].y +
1502 kcsi[k][j][i].z * kzet[k][j][i].z) +
1503 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1504 keta[k][j][i].y * kzet[k][j][i].y +
1505 keta[k][j][i].z * kzet[k][j][i].z) +
1506 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1507 kzet[k][j][i].y * kzet[k][j][i].y +
1508 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1509
1510 }
1511 }
1512 }
1513
1514
1515 //Mohsen March 2012//
1516
1517 // rhs.x at boundaries for periodic bc at i direction//
1518 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs==0){
1519 for (k=lzs; k<lze; k++) {
1520 for (j=lys; j<lye; j++) {
1521 i=xs;
1522 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1523
1524 dpdc = p[k][j][i+1] - p[k][j][i];
1525
1526 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) {
1527 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))) {
1528 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1529 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1530 }
1531 }
1532 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) {
1533 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1534 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1535 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1536 }
1537 }
1538 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) {
1539 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1540 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1541 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1542 }
1543 }
1544 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) {
1545 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1546 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1547 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1548 }
1549 }
1550 else {
1551 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1552 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1553 }
1554
1555 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) {
1556 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))) {
1557 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1558 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1559 }
1560 }
1561 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) {
1562 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1563 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1564 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1565 }
1566 }
1567 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) {
1568 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1569 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1570 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1571 }
1572 }
1573 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) {
1574 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1575 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1576 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1577 }
1578 }
1579 else {
1580 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1581 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1582 }
1583
1584 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1585 rhs[k][j][i].x -=
1586 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1587 icsi[k][j][i].y * icsi[k][j][i].y +
1588 icsi[k][j][i].z * icsi[k][j][i].z)+
1589 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1590 ieta[k][j][i].y * icsi[k][j][i].y +
1591 ieta[k][j][i].z * icsi[k][j][i].z)+
1592 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1593 izet[k][j][i].y * icsi[k][j][i].y +
1594 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1595 }
1596 }
1597 }
1598
1599// rhs.y at boundaries for periodic bc at j direction//
1600 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys==0){
1601 for (k=lzs; k<lze; k++) {
1602 for (i=lxs; i<lxe; i++) {
1603
1604 j=ys;
1605 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1606
1607 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) {
1608 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))) {
1609 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1610 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1611 }
1612 }
1613 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) {
1614 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1615 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1616 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1617 }
1618 }
1619 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) {
1620 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1621 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1622 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1623 }
1624 }
1625 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) {
1626 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1627 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1628 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1629 }
1630 }
1631 else {
1632 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1633 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1634 }
1635
1636 dpde = p[k][j+1][i] - p[k][j][i];
1637
1638 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) {
1639 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))) {
1640 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1641 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1642 }
1643 }
1644 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) {
1645 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1646 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1647 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1648 }
1649 }
1650 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) {
1651 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1652 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1653 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1654 }
1655 }
1656 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) {
1657 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1658 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1659 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1660 }
1661 }
1662 else {
1663 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1664 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1665 }
1666
1667 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1668
1669 rhs[k][j][i].y -=
1670 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1671 jcsi[k][j][i].y * jeta[k][j][i].y +
1672 jcsi[k][j][i].z * jeta[k][j][i].z)+
1673 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1674 jeta[k][j][i].y * jeta[k][j][i].y +
1675 jeta[k][j][i].z * jeta[k][j][i].z)+
1676 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1677 jzet[k][j][i].y * jeta[k][j][i].y +
1678 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1679
1680 }
1681 }
1682 }
1683
1684 // rhs.z at boundaries for periodic bc at k direction//
1685 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC&& zs==0){
1686 for (j=lys; j<lye; j++) {
1687 for (i=lxs; i<lxe; i++) {
1688
1689 k=zs;
1690 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1691
1692 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) {
1693 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))) {
1694 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1695 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1696 }
1697 }
1698 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) {
1699 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1700 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1701 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1702 }
1703 }
1704 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) {
1705 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1706 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1707 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1708 }
1709 }
1710 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) {
1711 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1712 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1713 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1714 }
1715 }
1716 else {
1717 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1718 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1719 }
1720
1721 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) {
1722 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))) {
1723 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1724 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1725 }
1726 }
1727 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) {
1728 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1729 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1730 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1731 }
1732 }
1733 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) {
1734 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1735 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1736 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1737 }
1738 }
1739 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) {
1740 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1741 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1742 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1743 }
1744 }
1745 else {
1746 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1747 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1748 }
1749
1750 dpdz = (p[k+1][j][i] - p[k][j][i]);
1751
1752 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1753
1754 rhs[k][j][i].z -=
1755 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1756 kcsi[k][j][i].y * kzet[k][j][i].y +
1757 kcsi[k][j][i].z * kzet[k][j][i].z)+
1758 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1759 keta[k][j][i].y * kzet[k][j][i].y +
1760 keta[k][j][i].z * kzet[k][j][i].z)+
1761 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1762 kzet[k][j][i].y * kzet[k][j][i].y +
1763 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1764
1765 }
1766 }
1767 }
1768
1769 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1770
1771 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1772 PetscInt TwoD = simCtx->TwoD;
1773
1774 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1775
1776 // 7. Final clean-up for immersed boundaries and 2D cases
1777 for (k=lzs; k<lze; k++) {
1778 for (j=lys; j<lye; j++) {
1779 for (i=lxs; i<lxe; i++) {
1780 if (TwoD==1)
1781 rhs[k][j][i].x =0.;
1782 else if (TwoD==2)
1783 rhs[k][j][i].y =0.;
1784 else if (TwoD==3)
1785 rhs[k][j][i].z =0.;
1786
1787 if (nvert[k][j][i]>0.1) {
1788 rhs[k][j][i].x = 0;
1789 rhs[k][j][i].y = 0;
1790 rhs[k][j][i].z = 0;
1791 }
1792 if (nvert[k][j][i+1]>0.1) {
1793 rhs[k][j][i].x=0;
1794 }
1795 if (nvert[k][j+1][i]>0.1) {
1796 rhs[k][j][i].y=0;
1797 }
1798 if (nvert[k+1][j][i]>0.1) {
1799 rhs[k][j][i].z=0;
1800 }
1801 }
1802 }
1803 }
1804 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1805
1806 // ========================================================================
1807
1808 // --- Restore all PETSc array pointers ---
1809 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1810 ierr = DMDAVecRestoreArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1811 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1812
1813 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1814 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1815 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1816 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1817 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1818
1819 ierr = DMDAVecRestoreArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1820 ierr = DMDAVecRestoreArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1821 ierr = DMDAVecRestoreArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1822 ierr = DMDAVecRestoreArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1823 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1824
1825 ierr = DMDAVecRestoreArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1826 ierr = DMDAVecRestoreArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1827 ierr = DMDAVecRestoreArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1828 ierr = DMDAVecRestoreArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1829 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1830
1831 ierr = DMDAVecRestoreArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1832 ierr = DMDAVecRestoreArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1833 ierr = DMDAVecRestoreArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1834 ierr = DMDAVecRestoreArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1835 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1836
1837 ierr = DMDAVecRestoreArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1838 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1839
1840 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1841 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1842
1843 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1844
1845 // --- Destroy temporary work vectors ---
1846 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1847 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1848 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1849 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1850 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1851
1852 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1853 simCtx->rank, user->_this);
1854
1855 PROFILE_FUNCTION_END;
1856 PetscFunctionReturn(0);
1857}
SynchronizePeriodicLocalStaggeredField
PetscErrorCode SynchronizePeriodicLocalStaggeredField(UserCtx *user, Vec local_field)
Synchronizes one local-only component-staggered periodic work field.
Definition Boundaries.c:2126
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:1647
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
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
ComputeBodyForces
PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
Internal helper implementation: ComputeBodyForces().
Definition rhs.c:1129
Viscous
PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
Implementation of Viscous().
Definition rhs.c:435
Convection
PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
Implementation of Convection().
Definition rhs.c:14
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:3300
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:2489
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:892
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:892
SimCtx::rank
PetscMPIInt rank
Definition variables.h:862
UserCtx::_this
PetscInt _this
Definition variables.h:1092
SimCtx::invicid
PetscInt invicid
Definition variables.h:892
UserCtx::lUcont
Vec lUcont
Definition variables.h:1113
UserCtx::info
DMDALocalInfo info
Definition variables.h:1086
UserCtx::lUcat
Vec lUcat
Definition variables.h:1113
SimCtx::rotateframe
PetscInt rotateframe
moveframe/rotateframe are refused at setup.
Definition variables.h:892
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeEulerianDiffusivity()

PetscErrorCode ComputeEulerianDiffusivity ( UserCtx *  user)

Implementation of ComputeEulerianDiffusivity().

Computes the effective diffusivity scalar field (Gamma_eff) on the Eulerian grid.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/rhs.h.

See also
ComputeEulerianDiffusivity()

Definition at line 1867 of file rhs.c.

1868{
1869 PetscErrorCode ierr;
1870 DM da = user->da;
1871 PetscInt i, j, k, xs, ys, zs, xm, ym, zm, xe, ye, ze;
1872 PetscInt lxs, lys, lzs, lxe, lye, lze;
1873
1874 // Pointers for 3D grid access
1875 PetscReal ***diff_arr; // Output: Diffusivity field
1876 PetscReal ***nut_arr; // Input: Eddy Viscosity field (optional)
1877
1878 // Physics parameters
1879 PetscReal nu_molecular, gamma_molecular;
1880 PetscReal nu_turbulent, gamma_turbulent;
1881 PetscReal Sc, Sct;
1882 PetscBool use_turbulence_model;
1883
1884 PetscFunctionBeginUser;
1885 PROFILE_FUNCTION_BEGIN;
1886
1887 if (VerificationDiffusivityOverrideActive(user->simCtx)) {
1888 ierr = ApplyVerificationDiffusivityOverride(user); CHKERRQ(ierr);
1889 PROFILE_FUNCTION_END;
1890 PetscFunctionReturn(0);
1891 }
1892
1893 // ------------------------------------------------------------------------
1894 // 1. Parameter Setup & Safety Checks
1895 // ------------------------------------------------------------------------
1896
1897 // Determine Molecular Viscosity (nu = 1/Re)
1898 // Guard against division by zero if Re is not set or infinite (inviscid)
1899 if (user->simCtx->ren > 1.0e-12) {
1900 nu_molecular = 1.0 / user->simCtx->ren;
1901 } else {
1902 nu_molecular = 0.0;
1903 }
1904
1905 // Set Schmidt Numbers (Default to 1.0 if not provided to prevent NaN)
1906 Sc = (user->simCtx->schmidt_number > 1.0e-6) ? user->simCtx->schmidt_number : 1.0;
1907 Sct = (user->simCtx->Turbulent_schmidt_number > 1.0e-6) ? user->simCtx->Turbulent_schmidt_number : 0.7;
1908
1909 // Pre-calculate molecular component
1910 gamma_molecular = nu_molecular / Sc;
1911
1912 // Check if a turbulence model is active
1913 use_turbulence_model = user->simCtx->les ? PETSC_TRUE : PETSC_FALSE;
1914
1915 // ------------------------------------------------------------------------
1916 // 2. Data Access
1917 // ------------------------------------------------------------------------
1918
1919 // Get local grid boundaries
1920 DMDALocalInfo info;
1921 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
1922
1923 xs = info.xs; ys = info.ys; zs = info.zs;
1924 xm = info.xm; ym = info.ym; zm = info.zm;
1925 xe = xs + xm; ye = ys + ym; ze = zs + zm;
1926
1927 lxs = (xs == 0)? xs + 1 : xs;
1928 lys = (ys == 0)? ys + 1 : ys;
1929 lzs = (zs == 0)? zs + 1 : zs;
1930 lxe = (xe == info.mx)? xe - 1 : xe;
1931 lye = (ye == info.my)? ye - 1 : ye;
1932 lze = (ze == info.mz)? ze - 1 : ze;
1933
1934 // Get write access to the output Diffusivity array
1935 ierr = DMDAVecGetArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
1936
1937 // Get read access to Eddy Viscosity only if turbulence is active
1938 if (use_turbulence_model) {
1939 ierr = DMDAVecGetArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
1940 }
1941
1942 // ------------------------------------------------------------------------
1943 // 3. Calculation Loop
1944 // ------------------------------------------------------------------------
1945
1946 for (k = lzs; k < lze; k++) {
1947 for (j = lys; j < lye; j++) {
1948 for (i = lxs; i < lxe; i++) {
1949
1950 gamma_turbulent = 0.0;
1951
1952 if (use_turbulence_model) {
1953 // Fetch local eddy viscosity
1954 nu_turbulent = nut_arr[k][j][i];
1955
1956 // NUMERICAL SAFETY:
1957 // Some turbulence models (dynamic SGS) can locally produce
1958 // slightly negative viscosity. We clamp this to 0 to prevent
1959 // negative diffusivity, which crashes the Langevin sqrt().
1960 if (nu_turbulent < 0.0) {
1961 nu_turbulent = 0.0;
1962 }
1963
1964 gamma_turbulent = nu_turbulent / Sct;
1965 }
1966
1967 // Sum components
1968 diff_arr[k][j][i] = gamma_molecular + gamma_turbulent;
1969 }
1970 }
1971 }
1972
1973 // ------------------------------------------------------------------------
1974 // 4. Cleanup & Synchronization
1975 // ------------------------------------------------------------------------
1976
1977 // Restore arrays
1978 if (use_turbulence_model) {
1979 ierr = DMDAVecRestoreArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
1980 }
1981 ierr = DMDAVecRestoreArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
1982
1983 // Update Ghost Points
1984 // This is required because downstream operations (Drift Gradient Calculation
1985 // and Particle Interpolation) will need access to the halo regions of this field.
1986 const FieldId periodic_fields[] = {FIELD_ID_DIFFUSIVITY};
1987 ierr = SynchronizePeriodicCellFields(user, 1, periodic_fields); CHKERRQ(ierr);
1988 ierr = UpdateLocalGhosts(user, FIELD_ID_DIFFUSIVITY); CHKERRQ(ierr);
1989 PROFILE_FUNCTION_END;
1990 PetscFunctionReturn(0);
1991}
FIELD_ID_DIFFUSIVITY
@ FIELD_ID_DIFFUSIVITY
Definition field_catalog.h:35
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:948
SimCtx::Turbulent_schmidt_number
PetscReal Turbulent_schmidt_number
Definition variables.h:948
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:1116
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeEulerianDiffusivityGradient()

PetscErrorCode ComputeEulerianDiffusivityGradient ( UserCtx *  user)

Internal helper implementation: ComputeEulerianDiffusivityGradient().

Computes the Eulerian gradient of the effective diffusivity field.

Local to this translation unit.

Definition at line 1999 of file rhs.c.

2000{
2001 PetscErrorCode ierr;
2002 DM da = user->da, fda = user->fda;
2003 DMDALocalInfo info = user->info;
2004
2005 // 1. Determine Global Dimensions
2006 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2007
2008 // 2. Determine Local Loop Bounds (Skip Ghosts/Unused Indices)
2009 // Grid uses indices 1 to M-2 for physical cells.
2010 // Index 0 and M-1 are ghost/boundary holders.
2011
2012 // Start: If we own the global start (0), skip it and start at 1.
2013 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2014 // End: If we own the global end (mx), stop before it (mx-1), so loop covers mx-2.
2015 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
2016
2017 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2018 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
2019
2020 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2021 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
2022
2023 PetscInt i, j, k;
2024
2025 // Pointers
2026 PetscReal ***diff; // Input (Scalar)
2027 Cmpnts ***grad_diff; // Output (Vector)
2028 Cmpnts ***csi, ***eta, ***zet;
2029 PetscReal ***aj;
2030
2031 // Boundary Flags (Check if Periodic)
2032 PetscBool p_x = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
2033 PetscBool p_y = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
2034 PetscBool p_z = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
2035
2036 PetscFunctionBeginUser;
2037 PROFILE_FUNCTION_BEGIN;
2038
2039 // 3. Update Ghosts for Diffusivity
2040 // Required so that Central Differences at i=2 can read i=1,
2041 // and Central Differences at Periodic Boundaries work correctly.
2042 ierr = UpdateLocalGhosts(user, FIELD_ID_DIFFUSIVITY); CHKERRQ(ierr);
2043
2044 // 4. Get Arrays (Read Only)
2045 ierr = DMDAVecGetArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2046 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2047 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2048 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2049 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2050
2051 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity and Metrics arrays accessed successfully! .\n");
2052
2053 // 5. Get Output Array (Read/Write)
2054 ierr = DMDAVecGetArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2055
2056 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity Gradient array accessed successfully! .\n");
2057
2058 // 6. Loop over Physical Domain
2059 for (k = lzs; k < lze; k++) {
2060 for (j = lys; j < lye; j++) {
2061 for (i = lxs; i < lxe; i++) {
2062
2063 PetscReal dGdCsi, dGdEta, dGdZet;
2064
2065 // ---------------------------------------------------------
2066 // I-Direction (Csi)
2067 // ---------------------------------------------------------
2068 if (!p_x && i == 1) {
2069 // Physical Start: 2nd Order Forward Difference
2070 // Stencil: [-3, 4, -1] / 2
2071 dGdCsi = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j][i+1] - diff[k][j][i+2]) * 0.5;
2072 }
2073 else if (!p_x && i == mx - 2) {
2074 // Physical End: 2nd Order Backward Difference
2075 // Stencil: [1, -4, 3] / 2
2076 dGdCsi = (3.0 * diff[k][j][i] - 4.0 * diff[k][j][i-1] + diff[k][j][i-2]) * 0.5;
2077 }
2078 else {
2079 // Interior / Periodic: Central Difference
2080 // Stencil: [-1, 0, 1] / 2
2081 dGdCsi = (diff[k][j][i+1] - diff[k][j][i-1]) * 0.5;
2082 }
2083
2084 // ---------------------------------------------------------
2085 // J-Direction (Eta)
2086 // ---------------------------------------------------------
2087 if (!p_y && j == 1) {
2088 // Physical Start: Forward
2089 dGdEta = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j+1][i] - diff[k][j+2][i]) * 0.5;
2090 }
2091 else if (!p_y && j == my - 2) {
2092 // Physical End: Backward
2093 dGdEta = (3.0 * diff[k][j][i] - 4.0 * diff[k][j-1][i] + diff[k][j-2][i]) * 0.5;
2094 }
2095 else {
2096 // Interior: Central
2097 dGdEta = (diff[k][j+1][i] - diff[k][j-1][i]) * 0.5;
2098 }
2099
2100 // ---------------------------------------------------------
2101 // K-Direction (Zet)
2102 // ---------------------------------------------------------
2103 if (!p_z && k == 1) {
2104 // Physical Start: Forward
2105 dGdZet = (-3.0 * diff[k][j][i] + 4.0 * diff[k+1][j][i] - diff[k+2][j][i]) * 0.5;
2106 }
2107 else if (!p_z && k == mz - 2) {
2108 // Physical End: Backward
2109 dGdZet = (3.0 * diff[k][j][i] - 4.0 * diff[k-1][j][i] + diff[k-2][j][i]) * 0.5;
2110 }
2111 else {
2112 // Interior: Central
2113 dGdZet = (diff[k+1][j][i] - diff[k-1][j][i]) * 0.5;
2114 }
2115
2116 // ---------------------------------------------------------
2117 // Transform to Physical Space (Cartesian Gradient)
2118 // ---------------------------------------------------------
2119 TransformScalarDerivativesToPhysical(aj[k][j][i],
2120 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2121 dGdCsi, dGdEta, dGdZet,
2122 &grad_diff[k][j][i]);
2123 }
2124 }
2125 }
2126
2127 // 7. Restore Arrays
2128 ierr = DMDAVecRestoreArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2129 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2130 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2131 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2132 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2133 ierr = DMDAVecRestoreArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2134
2135 // 8. Update Ghosts for the Result
2136 // Important: Particles near subdomain boundaries will need to interpolate
2137 // this gradient vector, so the ghosts of the vector field must be filled.
2138 ierr = UpdateLocalGhosts(user, FIELD_ID_DIFFUSIVITY_GRADIENT); CHKERRQ(ierr);
2139
2140 PROFILE_FUNCTION_END;
2141 PetscFunctionReturn(0);
2142}
FIELD_ID_DIFFUSIVITY_GRADIENT
@ FIELD_ID_DIFFUSIVITY_GRADIENT
Definition field_catalog.h:36
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:3983
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:1117
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:1116
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8