PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
les.c
Go to the documentation of this file.
1/**
2 * @file les.c
3 * @brief Implements the Large Eddy Simulation (LES) turbulence models.
4 *
5 * This file contains the core logic for the dynamic Smagorinsky LES model.
6 * It is responsible for two primary tasks, orchestrated by the main FlowSolver:
7 * 1. `ComputeSmagorinskyConstant`: Dynamically calculates the Smagorinsky coefficient (Cs)
8 * at each grid point by applying a test filter to the resolved flow field. This
9 * is computationally intensive and is typically not performed every time step.
10 * 2. `ComputeEddyViscosityLES`: Uses the calculated Cs to compute the turbulent
11 * eddy viscosity (nu_t), which is then added to the molecular viscosity in the
12 * momentum equations to model the dissipative effects of sub-grid scale turbulence.
13 */
14
15#include "les.h"
16
17// A small constant to prevent division by zero in sensitive calculations.
18const double LES_EPSILON = 1.0e-12;
19
20
21#undef __FUNCT__
22#define __FUNCT__ "ComputeSmagorinskyConstant"
23PetscErrorCode ComputeSmagorinskyConstant(UserCtx *user)
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}
297
298
299#undef __FUNCT__
300#define __FUNCT__ "ComputeEddyViscosityLES"
301PetscErrorCode ComputeEddyViscosityLES(UserCtx *user)
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}
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
ComputeEddyViscosityLES
PetscErrorCode ComputeEddyViscosityLES(UserCtx *user)
Computes the turbulent eddy viscosity (Nu_t) for the LES model.
Definition les.c:301
ComputeSmagorinskyConstant
PetscErrorCode ComputeSmagorinskyConstant(UserCtx *user)
Computes the dynamic Smagorinsky constant (Cs) for the LES model.
Definition les.c:23
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
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
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::lNu_t
Vec lNu_t
Definition variables.h:712
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:712
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
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:661