PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
les.h File Reference
#include "variables.h"
#include "logging.h"
#include "Metric.h"
#include "setup.h"
#include "Filter.h"
Include dependency graph for les.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

PetscErrorCode ComputeSmagorinskyConstant (UserCtx *user)
 Computes the dynamic Smagorinsky constant (Cs) for the LES model.
 
PetscErrorCode ComputeEddyViscosityLES (UserCtx *user)
 Computes the turbulent eddy viscosity (Nu_t) for the LES model.
 

Function Documentation

◆ ComputeSmagorinskyConstant()

PetscErrorCode ComputeSmagorinskyConstant ( UserCtx *  user)

Computes the dynamic Smagorinsky constant (Cs) for the LES model.

This function implements the dynamic procedure for calculating the Smagorinsky constant at each grid point. It involves test-filtering the velocity field and strain rate tensor to determine the constant locally. The calculation can be computationally expensive and is typically run less frequently than every time step, controlled by simCtx->dynamic_freq.

Parameters
userThe user context for a single computational block.
Returns
PetscErrorCode 0 on success.

Definition at line 23 of file les.c.

24{
25 PetscErrorCode ierr;
26 SimCtx *simCtx = user->simCtx; // Get the global context for simulation parameters
27
28 PetscFunctionBeginUser;
29 PROFILE_FUNCTION_BEGIN;
30
31 // --- Initial Condition / Simple Model Handling ---
32
33 // For the first couple of steps of a simulation from t=0, the flow is not developed
34 // enough for the dynamic procedure to be stable. Set Cs to zero.
35 if (simCtx->step < 2 && simCtx->StartStep == 0) {
36 ierr = VecSet(user->CS, 0.0); CHKERRQ(ierr);
37 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Setting Smagorinsky coefficient Cs=0.0 for initial steps (step=%d)\n", simCtx->step);
38 PROFILE_FUNCTION_END;
39 PetscFunctionReturn(0);
40 }
41
42 // If the user requests the non-dynamic, constant-coefficient Smagorinsky model (les=1),
43 // set a constant value and exit.
44 if (simCtx->les == 1) {
45 LOG_ALLOW(GLOBAL,LOG_INFO,"Using constant-coefficient Smagorinsky model with Cs=%.4f \n", simCtx->Const_CS);
46 ierr = VecSet(user->CS, simCtx->Const_CS); CHKERRQ(ierr); // A typical constant value
47 PROFILE_FUNCTION_END;
48 PetscFunctionReturn(0);
49 }
50
51 // --- Setup for Dynamic Procedure ---
52
53 DM da = user->da, fda = user->fda;
54 DMDALocalInfo info;
55 PetscInt i, j, k, p, q, r;
56 PetscInt xs, xe, ys, ye, zs, ze; // Local grid corner indices
57 PetscInt mx, my, mz; // Global grid dimensions
58 PetscInt lxs, lxe, lys, lye, lzs, lze; // Local interior loop bounds
59
60 // Temporary PETSc vectors required for the calculation. These are allocated
61 // here and destroyed at the end to keep the UserCtx clean.
62 Vec lSx, lSy, lSz, lS_abs;
63 Vec lLM, lMM; // Leonard & cross-term stress tensors
64
65 // Get grid dimensions and local bounds
66 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
67 mx = info.mx; my = info.my; mz = info.mz;
68 xs = info.xs; xe = xs + info.xm;
69 ys = info.ys; ye = ys + info.ym;
70 zs = info.zs; ze = zs + info.zm;
71
72 // Define loop bounds to exclude physical boundaries where stencils are invalid.
73 lxs = xs; lxe = xe; lys = ys; lye = ye; lzs = zs; lze = ze;
74 if (xs==0) lxs = xs+1;
75 if (ys==0) lys = ys+1;
76 if (zs==0) lzs = zs+1;
77 if (xe==mx) lxe = xe-1;
78 if (ye==my) lye = ye-1;
79 if (ze==mz) lze = ze-1;
80
81 // Allocate temporary vectors
82 ierr = VecDuplicate(user->lUcont, &lSx); CHKERRQ(ierr);
83 ierr = VecDuplicate(user->lUcont, &lSy); CHKERRQ(ierr);
84 ierr = VecDuplicate(user->lUcont, &lSz); CHKERRQ(ierr);
85 ierr = VecDuplicate(user->lP, &lS_abs); CHKERRQ(ierr);
86 ierr = VecDuplicate(user->lP, &lLM); CHKERRQ(ierr);
87 ierr = VecDuplicate(user->lP, &lMM); CHKERRQ(ierr);
88 ierr = VecSet(lLM, 0.0); CHKERRQ(ierr);
89 ierr = VecSet(lMM, 0.0); CHKERRQ(ierr);
90
91 // Get read/write access to PETSc vector data as multidimensional arrays
92 Cmpnts ***ucat, ***csi, ***eta, ***zet;
93 PetscReal ***nvert, ***Cs_arr, ***aj, ***S_abs_arr, ***LM_arr, ***MM_arr;
94 Cmpnts ***Sx_arr, ***Sy_arr, ***Sz_arr;
95
96 ierr = DMDAVecGetArray(fda, user->lUcat, &ucat); CHKERRQ(ierr);
97 ierr = DMDAVecGetArrayRead(fda, user->lCsi, (Cmpnts***)&csi); CHKERRQ(ierr);
98 ierr = DMDAVecGetArrayRead(fda, user->lEta, (Cmpnts***)&eta); CHKERRQ(ierr);
99 ierr = DMDAVecGetArrayRead(fda, user->lZet, (Cmpnts***)&zet); CHKERRQ(ierr);
100 ierr = DMDAVecGetArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
101 ierr = DMDAVecGetArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
102 ierr = DMDAVecGetArray(da, user->CS, &Cs_arr); CHKERRQ(ierr);
103
104 ierr = DMDAVecGetArray(fda, lSx, &Sx_arr); CHKERRQ(ierr);
105 ierr = DMDAVecGetArray(fda, lSy, &Sy_arr); CHKERRQ(ierr);
106 ierr = DMDAVecGetArray(fda, lSz, &Sz_arr); CHKERRQ(ierr);
107 ierr = DMDAVecGetArray(da, lS_abs, &S_abs_arr); CHKERRQ(ierr);
108 ierr = DMDAVecGetArray(da, lLM, &LM_arr); CHKERRQ(ierr);
109 ierr = DMDAVecGetArray(da, lMM, &MM_arr); CHKERRQ(ierr);
110
111 // --- 1. Compute and store the strain rate tensor |S| and velocity gradients at all points ---
112 for (k=lzs; k<lze; k++)
113 for (j=lys; j<lye; j++)
114 for (i=lxs; i<lxe; i++) {
115 if( nvert[k][j][i] > 1.1) continue; // Skip points inside solid bodies
116
117 Cmpnts dudx, dvdx, dwdx;
118 ierr = ComputeVectorFieldDerivatives(user, i, j, k, ucat, &dudx, &dvdx, &dwdx); CHKERRQ(ierr);
119
120 double Sxx = dudx.x;
121 double Sxy = 0.5 * (dudx.y + dvdx.x);
122 double Sxz = 0.5 * (dudx.z + dwdx.x);
123 double Syy = dvdx.y;
124 double Syz = 0.5 * (dvdx.z + dwdx.y);
125 double Szz = dwdx.z;
126 double Syx = Sxy, Szx = Sxz, Szy = Syz; // Symmetry
127
128 S_abs_arr[k][j][i] = sqrt( 2.0 * (Sxx*Sxx + Sxy*Sxy + Sxz*Sxz +
129 Syx*Syx + Syy*Syy + Syz*Syz +
130 Szx*Szx + Szy*Szy + Szz*Szz) );
131
132 // Store the full velocity gradient tensor for use in the next step
133 Sx_arr[k][j][i] = dudx; // Contains {du/dx, du/dy, du/dz}
134 Sy_arr[k][j][i] = dvdx; // Contains {dv/dx, dv/dy, dv/dz}
135 Sz_arr[k][j][i] = dwdx; // Contains {dw/dx, dw/dy, dw/dz}
136 }
137
138 // --- 2. Main loop to compute the dynamic coefficient ---
139 for (k=lzs; k<lze; k++)
140 for (j=lys; j<lye; j++)
141 for (i=lxs; i<lxe; i++) {
142 if(nvert[k][j][i] > 1.1) {
143 LM_arr[k][j][i] = 0.0;
144 MM_arr[k][j][i] = 0.0;
145 continue;
146 }
147
148 // --- 2a. Gather data from the 3x3x3 stencil around the current point (i,j,k) ---
149 double u[3][3][3], v[3][3][3], w[3][3][3];
150 double S_abs_stencil[3][3][3];
151 double S11[3][3][3], S12[3][3][3], S13[3][3][3];
152 double S22[3][3][3], S23[3][3][3], S33[3][3][3];
153 double weights[3][3][3];
154
155 for(r=-1; r<=1; r++)
156 for(q=-1; q<=1; q++)
157 for(p=-1; p<=1; p++) {
158 int R=r+1, Q=q+1, P=p+1; // Stencil array indices (0-2)
159 int KK=k+r, JJ=j+q, II=i+p; // Global grid indices
160
161 u[R][Q][P] = ucat[KK][JJ][II].x;
162 v[R][Q][P] = ucat[KK][JJ][II].y;
163 w[R][Q][P] = ucat[KK][JJ][II].z;
164
165 // Strain rate tensor components (Sij = 0.5 * (dui/dxj + duj/dxi))
166 S11[R][Q][P] = Sx_arr[KK][JJ][II].x; // du/dx
167 S12[R][Q][P] = 0.5 * (Sx_arr[KK][JJ][II].y + Sy_arr[KK][JJ][II].x); // 0.5 * (du/dy + dv/dx)
168 S13[R][Q][P] = 0.5 * (Sx_arr[KK][JJ][II].z + Sz_arr[KK][JJ][II].x); // 0.5 * (du/dz + dw/dx)
169 S22[R][Q][P] = Sy_arr[KK][JJ][II].y; // dv/dy
170 S23[R][Q][P] = 0.5 * (Sy_arr[KK][JJ][II].z + Sz_arr[KK][JJ][II].y); // 0.5 * (dv/dz + dw/dy)
171 S33[R][Q][P] = Sz_arr[KK][JJ][II].z; // dw/dz
172
173 S_abs_stencil[R][Q][P] = S_abs_arr[KK][JJ][II];
174 weights[R][Q][P] = aj[KK][JJ][II]; // Weight is Jacobian (1/volume)
175 }
176
177 // --- 2b. Apply the test filter to all required quantities ---
178 double u_filt = ApplyLESTestFilter(simCtx, u, weights);
179 double v_filt = ApplyLESTestFilter(simCtx, v, weights);
180 double w_filt = ApplyLESTestFilter(simCtx, w, weights);
181
182 double S11_filt = ApplyLESTestFilter(simCtx, S11, weights);
183 double S12_filt = ApplyLESTestFilter(simCtx, S12, weights);
184 double S13_filt = ApplyLESTestFilter(simCtx, S13, weights);
185 double S22_filt = ApplyLESTestFilter(simCtx, S22, weights);
186 double S23_filt = ApplyLESTestFilter(simCtx, S23, weights);
187 double S33_filt = ApplyLESTestFilter(simCtx, S33, weights);
188
189 double S_abs_filt = ApplyLESTestFilter(simCtx, S_abs_stencil, weights);
190
191 // Filter quadratic terms: <u*u>, <u*v>, etc.
192 double uu[3][3][3], uv[3][3][3], uw[3][3][3], vv[3][3][3], vw[3][3][3], ww[3][3][3];
193 for(r=0; r<3; r++) for(q=0; q<3; q++) for(p=0; p<3; p++) {
194 uu[r][q][p] = u[r][q][p] * u[r][q][p];
195 uv[r][q][p] = u[r][q][p] * v[r][q][p];
196 uw[r][q][p] = u[r][q][p] * w[r][q][p];
197 vv[r][q][p] = v[r][q][p] * v[r][q][p];
198 vw[r][q][p] = v[r][q][p] * w[r][q][p];
199 ww[r][q][p] = w[r][q][p] * w[r][q][p];
200 }
201 double uu_filt = ApplyLESTestFilter(simCtx, uu, weights);
202 double uv_filt = ApplyLESTestFilter(simCtx, uv, weights);
203 double uw_filt = ApplyLESTestFilter(simCtx, uw, weights);
204 double vv_filt = ApplyLESTestFilter(simCtx, vv, weights);
205 double vw_filt = ApplyLESTestFilter(simCtx, vw, weights);
206 double ww_filt = ApplyLESTestFilter(simCtx, ww, weights);
207
208 // --- 2c. Compute the Leonard stress tensor (L_ij = <u_i*u_j> - <u_i>*<u_j>) ---
209 double L11 = uu_filt - u_filt * u_filt;
210 double L12 = uv_filt - u_filt * v_filt;
211 double L13 = uw_filt - u_filt * w_filt;
212 double L22 = vv_filt - v_filt * v_filt;
213 double L23 = vw_filt - v_filt * w_filt;
214 double L33 = ww_filt - w_filt * w_filt;
215
216 // --- 2d. Compute the model tensor M_ij ---
217 double grid_filter_width, test_filter_width;
218 ierr = ComputeCellCharacteristicLengthScale(aj[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i], &grid_filter_width, &test_filter_width, &test_filter_width); CHKERRQ(ierr); // Simplified for now
219 grid_filter_width = pow(1.0 / aj[k][j][i], 1.0/3.0);
220 test_filter_width = 2.0 * grid_filter_width; // Standard box filter definition
221
222 double alpha = pow(test_filter_width / grid_filter_width, 2.0);
223
224 double M11 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S11_filt - S_abs_filt * S11_filt);
225 double M12 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S12_filt - S_abs_filt * S12_filt);
226 double M13 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S13_filt - S_abs_filt * S13_filt);
227 double M22 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S22_filt - S_abs_filt * S22_filt);
228 double M23 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S23_filt - S_abs_filt * S23_filt);
229 double M33 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S33_filt - S_abs_filt * S33_filt);
230
231 // --- 2e. Contract tensors to find Cs^2 (L_ij * M_ij / (M_kl * M_kl)) ---
232 LM_arr[k][j][i] = L11*M11 + 2*L12*M12 + 2*L13*M13 + L22*M22 + 2*L23*M23 + L33*M33;
233 MM_arr[k][j][i] = M11*M11 + 2*M12*M12 + 2*M13*M13 + M22*M22 + 2*M23*M23 + M33*M33;
234 }
235
236 // --- 3. Averaging and finalization ---
237
238 // At this point, LM_arr and MM_arr contain raw, unaveraged values.
239 // To stabilize the solution, these are often averaged over homogeneous directions or locally.
240 // The logic for homogeneous averaging would go here if simCtx->i_homo_filter is true.
241 // For this general version, we proceed with the local (unaveraged) values.
242
243 for (k=lzs; k<lze; k++)
244 for (j=lys; j<lye; j++)
245 for (i=lxs; i<lxe; i++) {
246 if(nvert[k][j][i] > 1.1) {
247 Cs_arr[k][j][i] = 0.0;
248 continue;
249 }
250
251 double C_sq = 0.0;
252 if (MM_arr[k][j][i] > LES_EPSILON) {
253 C_sq = LM_arr[k][j][i] / MM_arr[k][j][i];
254 }
255
256 // Final clipping and assignment: Cs must be positive and not excessively large.
257 double Cs_val = (C_sq > 0.0) ? sqrt(C_sq) : 0.0;
258 Cs_arr[k][j][i] = PetscMin(PetscMax(Cs_val, 0.0), simCtx->max_cs);
259 }
260
261 // --- 4. Cleanup ---
262
263 // Restore all PETSc arrays
264 ierr = DMDAVecRestoreArray(fda, user->lUcat, &ucat); CHKERRQ(ierr);
265 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, (Cmpnts***)&csi); CHKERRQ(ierr);
266 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, (Cmpnts***)&eta); CHKERRQ(ierr);
267 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, (Cmpnts***)&zet); CHKERRQ(ierr);
268 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
269 ierr = DMDAVecRestoreArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
270 ierr = DMDAVecRestoreArray(da, user->CS, &Cs_arr); CHKERRQ(ierr);
271 ierr = DMDAVecRestoreArray(fda, lSx, &Sx_arr); CHKERRQ(ierr);
272 ierr = DMDAVecRestoreArray(fda, lSy, &Sy_arr); CHKERRQ(ierr);
273 ierr = DMDAVecRestoreArray(fda, lSz, &Sz_arr); CHKERRQ(ierr);
274 ierr = DMDAVecRestoreArray(da, lS_abs, &S_abs_arr); CHKERRQ(ierr);
275 ierr = DMDAVecRestoreArray(da, lLM, &LM_arr); CHKERRQ(ierr);
276 ierr = DMDAVecRestoreArray(da, lMM, &MM_arr); CHKERRQ(ierr);
277
278 // Destroy temporary vectors to prevent memory leaks
279 ierr = VecDestroy(&lSx); CHKERRQ(ierr);
280 ierr = VecDestroy(&lSy); CHKERRQ(ierr);
281 ierr = VecDestroy(&lSz); CHKERRQ(ierr);
282 ierr = VecDestroy(&lS_abs); CHKERRQ(ierr);
283 ierr = VecDestroy(&lLM); CHKERRQ(ierr);
284 ierr = VecDestroy(&lMM); CHKERRQ(ierr);
285
286 // Communicate ghost point data for the newly computed Cs field
287 ierr = DMGlobalToLocalBegin(da, user->CS, INSERT_VALUES, user->lCs); CHKERRQ(ierr);
288 ierr = DMGlobalToLocalEnd(da, user->CS, INSERT_VALUES, user->lCs); CHKERRQ(ierr);
289
290 PetscReal max_norm;
291 ierr = VecMax(user->CS, NULL, &max_norm); CHKERRQ(ierr);
292 LOG_ALLOW(GLOBAL, LOG_INFO, " Max dynamic Smagorinsky constant (Cs) computed: %e (clip value: %e)\n", max_norm, simCtx->max_cs);
293
294 PROFILE_FUNCTION_END;
295 PetscFunctionReturn(0);
296}
ApplyLESTestFilter
double ApplyLESTestFilter(const SimCtx *simCtx, double values[3][3][3], double weights[3][3][3])
Applies a numerical "test filter" to a 3x3x3 stencil of data points.
Definition Filter.c:145
ComputeCellCharacteristicLengthScale
PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.
Definition Metric.c:313
LES_EPSILON
const double LES_EPSILON
Definition les.c:18
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
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:201
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
ComputeVectorFieldDerivatives
PetscErrorCode ComputeVectorFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
Computes the derivatives of a cell-centered vector field at a specific grid point.
Definition setup.c:2678
UserCtx::lNvert
Vec lNvert
Definition variables.h:688
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
UserCtx::lZet
Vec lZet
Definition variables.h:705
UserCtx::lCs
Vec lCs
Definition variables.h:712
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:548
Cmpnts::x
PetscScalar x
Definition variables.h:101
SimCtx::max_cs
PetscReal max_cs
Definition variables.h:611
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
Cmpnts::z
PetscScalar z
Definition variables.h:101
SimCtx::Const_CS
PetscReal Const_CS
Definition variables.h:611
UserCtx::lUcont
Vec lUcont
Definition variables.h:688
SimCtx::step
PetscInt step
Definition variables.h:546
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:705
SimCtx::les
PetscInt les
Definition variables.h:609
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:538
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeEddyViscosityLES()

PetscErrorCode ComputeEddyViscosityLES ( UserCtx *  user)

Computes the turbulent eddy viscosity (Nu_t) for the LES model.

Using the Smagorinsky constant (Cs) computed by ComputeSmagorinskyConstant, this function calculates the eddy viscosity at each grid point based on the local strain rate magnitude and filter width.

nu_t = (Cs * filter_width)^2 * |S|*

Parameters
userThe user context for a single computational block.
Returns
PetscErrorCode 0 on success.

Definition at line 301 of file les.c.

302{
303 PetscErrorCode ierr;
304 DM da = user->da, fda = user->fda;
305 DMDALocalInfo info;
306 PetscInt i, j, k;
307 PetscInt xs, xe, ys, ye, zs, ze;
308 PetscInt lxs, lxe, lys, lye, lzs, lze; // Local interior loop bounds
309 PetscInt mx, my, mz;
310
311 PetscFunctionBeginUser;
312 PROFILE_FUNCTION_BEGIN;
313
314 // Get grid dimensions and local bounds
315 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
316 mx = info.mx; my = info.my; mz = info.mz;
317 xs = info.xs; xe = xs + info.xm;
318 ys = info.ys; ye = ys + info.ym;
319 zs = info.zs; ze = zs + info.zm;
320
321 // Define loop bounds to exclude physical boundaries where stencils are invalid.
322 lxs = xs; lxe = xe; lys = ys; lye = ye; lzs = zs; lze = ze;
323 if (xs==0) lxs = xs+1;
324 if (ys==0) lys = ys+1;
325 if (zs==0) lzs = zs+1;
326 if (xe==mx) lxe = xe-1;
327 if (ye==my) lye = ye-1;
328 if (ze==mz) lze = ze-1;
329
330
331
332 // Get read/write access to PETSc data arrays
333 Cmpnts ***ucat;
334 PetscReal ***Cs_arr, ***nu_t_arr, ***nvert, ***aj;
335
336 ierr = DMDAVecGetArrayRead(fda, user->lUcat, (Cmpnts***)&ucat); CHKERRQ(ierr);
337 ierr = DMDAVecGetArrayRead(da, user->lCs, (PetscReal***)&Cs_arr); CHKERRQ(ierr);
338 ierr = DMDAVecGetArray(da, user->Nu_t, &nu_t_arr); CHKERRQ(ierr);
339 ierr = DMDAVecGetArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
340 ierr = DMDAVecGetArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
341
342 // Loop over the interior of the local domain
343 for (k=lzs; k<lze; k++)
344 for (j=lys; j<lye; j++)
345 for (i=lxs; i<lxe; i++) {
346 if(nvert[k][j][i] > 0.1) {
347 nu_t_arr[k][j][i] = 0.0;
348 continue;
349 }
350
351 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Computing eddy viscosity at point (%d,%d,%d)\n", i, j, k);
352 // 1. Compute the local strain rate magnitude |S|
353 Cmpnts dudx, dvdx, dwdx;
354 ierr = ComputeVectorFieldDerivatives(user, i, j, k, ucat, &dudx, &dvdx, &dwdx); CHKERRQ(ierr);
355
356 double Sxx = dudx.x;
357 double Sxy = 0.5 * (dudx.y + dvdx.x);
358 double Sxz = 0.5 * (dudx.z + dwdx.x);
359 double Syy = dvdx.y;
360 double Syz = 0.5 * (dvdx.z + dwdx.y);
361 double Szz = dwdx.z;
362 double strain_rate_mag = sqrt( 2.0 * (Sxx*Sxx + Sxy*Sxy + Sxz*Sxz + Sxy*Sxy + Syy*Syy + Syz*Syz + Sxz*Sxz + Syz*Syz + Szz*Szz) );
363
364 // 2. Determine the grid filter width, Delta = (cell volume)^(1/3)
365 double filter_width = pow( 1.0/aj[k][j][i], 1.0/3.0 );
366
367 // 3. Compute eddy viscosity: nu_t = (Cs * Delta)^2 * |S|
368 nu_t_arr[k][j][i] = pow(Cs_arr[k][j][i] * filter_width, 2.0) * strain_rate_mag;
369
370 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Cs=%.4e, Delta=%.4e, |S|=%.4e => nu_t=%.4e\n",
371 Cs_arr[k][j][i], filter_width, strain_rate_mag, nu_t_arr[k][j][i]);
372 }
373
374 // Restore PETSc data arrays
375 ierr = DMDAVecRestoreArrayRead(fda, user->lUcat, (Cmpnts***)&ucat); CHKERRQ(ierr);
376 ierr = DMDAVecRestoreArrayRead(da, user->lCs, (PetscReal***)&Cs_arr); CHKERRQ(ierr);
377 ierr = DMDAVecRestoreArray(da, user->Nu_t, &nu_t_arr); CHKERRQ(ierr);
378 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
379 ierr = DMDAVecRestoreArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
380
381 // Update ghost points for the newly computed eddy viscosity
382 ierr = DMGlobalToLocalBegin(da, user->Nu_t, INSERT_VALUES, user->lNu_t); CHKERRQ(ierr);
383 ierr = DMGlobalToLocalEnd(da, user->Nu_t, INSERT_VALUES, user->lNu_t); CHKERRQ(ierr);
384
385 PetscReal max_norm;
386 ierr = VecMax(user->Nu_t, NULL, &max_norm); CHKERRQ(ierr);
387 LOG_ALLOW(GLOBAL, LOG_INFO, " Max eddy viscosity (Nu_t) computed: %e\n", max_norm);
388
389 PROFILE_FUNCTION_END;
390 PetscFunctionReturn(0);
391}
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:712
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:712
Here is the call graph for this function:
Here is the caller graph for this function: