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"
23/**
24 * @brief Implementation of \ref ComputeSmagorinskyConstant().
25 * @details Full API contract (arguments, ownership, side effects) is documented with
26 * the header declaration in `include/les.h`.
27 * @see ComputeSmagorinskyConstant()
28 */
29
30PetscErrorCode ComputeSmagorinskyConstant(UserCtx *user)
31{
32 PetscErrorCode ierr;
33 SimCtx *simCtx = user->simCtx; // Get the global context for simulation parameters
34
35 PetscFunctionBeginUser;
36 PROFILE_FUNCTION_BEGIN;
37
38 // --- Initial Condition / Simple Model Handling ---
39
40 // For the first couple of steps of a simulation from t=0, the flow is not developed
41 // enough for the dynamic procedure to be stable. Set Cs to zero.
42 if (simCtx->step < 2 && simCtx->StartStep == 0) {
43 ierr = VecSet(user->CS, 0.0); CHKERRQ(ierr);
44 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Setting Smagorinsky coefficient Cs=0.0 for initial steps (step=%d)\n", simCtx->step);
45 PROFILE_FUNCTION_END;
46 PetscFunctionReturn(0);
47 }
48
49 // If the user requests the non-dynamic, constant-coefficient Smagorinsky model (les=1),
50 // set a constant value and exit.
51 if (simCtx->les == CONSTANT_SMAGORINSKY) {
52 LOG_ALLOW(GLOBAL,LOG_INFO,"Using constant-coefficient Smagorinsky model with Cs=%.4f \n", simCtx->Const_CS);
53 ierr = VecSet(user->CS, simCtx->Const_CS); CHKERRQ(ierr); // A typical constant value
54 PROFILE_FUNCTION_END;
55 PetscFunctionReturn(0);
56 }
57
58 // --- Setup for Dynamic Procedure ---
59
60 DM da = user->da, fda = user->fda;
61 DMDALocalInfo info;
62 PetscInt i, j, k, p, q, r;
63 PetscInt xs, xe, ys, ye, zs, ze; // Local grid corner indices
64 PetscInt mx, my, mz; // Global grid dimensions
65 PetscInt lxs, lxe, lys, lye, lzs, lze; // Local interior loop bounds
66
67 // Temporary PETSc vectors required for the calculation. These are allocated
68 // here and destroyed at the end to keep the UserCtx clean.
69 Vec lSx, lSy, lSz, lS_abs;
70 Vec lLM, lMM; // Leonard & cross-term stress tensors
71
72 // Get grid dimensions and local bounds
73 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
74 mx = info.mx; my = info.my; mz = info.mz;
75 xs = info.xs; xe = xs + info.xm;
76 ys = info.ys; ye = ys + info.ym;
77 zs = info.zs; ze = zs + info.zm;
78
79 // Define loop bounds to exclude physical boundaries where stencils are invalid.
80 lxs = xs; lxe = xe; lys = ys; lye = ye; lzs = zs; lze = ze;
81 if (xs==0) lxs = xs+1;
82 if (ys==0) lys = ys+1;
83 if (zs==0) lzs = zs+1;
84 if (xe==mx) lxe = xe-1;
85 if (ye==my) lye = ye-1;
86 if (ze==mz) lze = ze-1;
87
88 // Allocate temporary vectors
89 ierr = VecDuplicate(user->lUcont, &lSx); CHKERRQ(ierr);
90 ierr = VecDuplicate(user->lUcont, &lSy); CHKERRQ(ierr);
91 ierr = VecDuplicate(user->lUcont, &lSz); CHKERRQ(ierr);
92 ierr = VecDuplicate(user->lP, &lS_abs); CHKERRQ(ierr);
93 ierr = VecDuplicate(user->lP, &lLM); CHKERRQ(ierr);
94 ierr = VecDuplicate(user->lP, &lMM); CHKERRQ(ierr);
95 ierr = VecSet(lLM, 0.0); CHKERRQ(ierr);
96 ierr = VecSet(lMM, 0.0); CHKERRQ(ierr);
97
98 // Get read/write access to PETSc vector data as multidimensional arrays
99 Cmpnts ***ucat, ***csi, ***eta, ***zet;
100 PetscReal ***nvert, ***Cs_arr, ***aj, ***S_abs_arr, ***LM_arr, ***MM_arr;
101 Cmpnts ***Sx_arr, ***Sy_arr, ***Sz_arr;
102
103 ierr = DMDAVecGetArray(fda, user->lUcat, &ucat); CHKERRQ(ierr);
104 ierr = DMDAVecGetArrayRead(fda, user->lCsi, (Cmpnts***)&csi); CHKERRQ(ierr);
105 ierr = DMDAVecGetArrayRead(fda, user->lEta, (Cmpnts***)&eta); CHKERRQ(ierr);
106 ierr = DMDAVecGetArrayRead(fda, user->lZet, (Cmpnts***)&zet); CHKERRQ(ierr);
107 ierr = DMDAVecGetArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
108 ierr = DMDAVecGetArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
109 ierr = DMDAVecGetArray(da, user->CS, &Cs_arr); CHKERRQ(ierr);
110
111 ierr = DMDAVecGetArray(fda, lSx, &Sx_arr); CHKERRQ(ierr);
112 ierr = DMDAVecGetArray(fda, lSy, &Sy_arr); CHKERRQ(ierr);
113 ierr = DMDAVecGetArray(fda, lSz, &Sz_arr); CHKERRQ(ierr);
114 ierr = DMDAVecGetArray(da, lS_abs, &S_abs_arr); CHKERRQ(ierr);
115 ierr = DMDAVecGetArray(da, lLM, &LM_arr); CHKERRQ(ierr);
116 ierr = DMDAVecGetArray(da, lMM, &MM_arr); CHKERRQ(ierr);
117
118 // --- 1. Compute and store the strain rate tensor |S| and velocity gradients at all points ---
119 for (k=lzs; k<lze; k++)
120 for (j=lys; j<lye; j++)
121 for (i=lxs; i<lxe; i++) {
122 if( nvert[k][j][i] > 1.1) continue; // Skip points inside solid bodies
123
124 Cmpnts dudx, dvdx, dwdx;
125 ierr = ComputeVectorFieldDerivatives(user, i, j, k, ucat, &dudx, &dvdx, &dwdx); CHKERRQ(ierr);
126
127 double Sxx = dudx.x;
128 double Sxy = 0.5 * (dudx.y + dvdx.x);
129 double Sxz = 0.5 * (dudx.z + dwdx.x);
130 double Syy = dvdx.y;
131 double Syz = 0.5 * (dvdx.z + dwdx.y);
132 double Szz = dwdx.z;
133 double Syx = Sxy, Szx = Sxz, Szy = Syz; // Symmetry
134
135 S_abs_arr[k][j][i] = sqrt( 2.0 * (Sxx*Sxx + Sxy*Sxy + Sxz*Sxz +
136 Syx*Syx + Syy*Syy + Syz*Syz +
137 Szx*Szx + Szy*Szy + Szz*Szz) );
138
139 // Store the full velocity gradient tensor for use in the next step
140 Sx_arr[k][j][i] = dudx; // Contains {du/dx, du/dy, du/dz}
141 Sy_arr[k][j][i] = dvdx; // Contains {dv/dx, dv/dy, dv/dz}
142 Sz_arr[k][j][i] = dwdx; // Contains {dw/dx, dw/dy, dw/dz}
143 }
144
145 // --- 2. Main loop to compute the dynamic coefficient ---
146 for (k=lzs; k<lze; k++)
147 for (j=lys; j<lye; j++)
148 for (i=lxs; i<lxe; i++) {
149 if(nvert[k][j][i] > 1.1) {
150 LM_arr[k][j][i] = 0.0;
151 MM_arr[k][j][i] = 0.0;
152 continue;
153 }
154
155 // --- 2a. Gather data from the 3x3x3 stencil around the current point (i,j,k) ---
156 double u[3][3][3], v[3][3][3], w[3][3][3];
157 double S_abs_stencil[3][3][3];
158 double S11[3][3][3], S12[3][3][3], S13[3][3][3];
159 double S22[3][3][3], S23[3][3][3], S33[3][3][3];
160 double weights[3][3][3];
161
162 for(r=-1; r<=1; r++)
163 for(q=-1; q<=1; q++)
164 for(p=-1; p<=1; p++) {
165 int R=r+1, Q=q+1, P=p+1; // Stencil array indices (0-2)
166 int KK=k+r, JJ=j+q, II=i+p; // Global grid indices
167
168 u[R][Q][P] = ucat[KK][JJ][II].x;
169 v[R][Q][P] = ucat[KK][JJ][II].y;
170 w[R][Q][P] = ucat[KK][JJ][II].z;
171
172 // Strain rate tensor components (Sij = 0.5 * (dui/dxj + duj/dxi))
173 S11[R][Q][P] = Sx_arr[KK][JJ][II].x; // du/dx
174 S12[R][Q][P] = 0.5 * (Sx_arr[KK][JJ][II].y + Sy_arr[KK][JJ][II].x); // 0.5 * (du/dy + dv/dx)
175 S13[R][Q][P] = 0.5 * (Sx_arr[KK][JJ][II].z + Sz_arr[KK][JJ][II].x); // 0.5 * (du/dz + dw/dx)
176 S22[R][Q][P] = Sy_arr[KK][JJ][II].y; // dv/dy
177 S23[R][Q][P] = 0.5 * (Sy_arr[KK][JJ][II].z + Sz_arr[KK][JJ][II].y); // 0.5 * (dv/dz + dw/dy)
178 S33[R][Q][P] = Sz_arr[KK][JJ][II].z; // dw/dz
179
180 S_abs_stencil[R][Q][P] = S_abs_arr[KK][JJ][II];
181 weights[R][Q][P] = aj[KK][JJ][II]; // Weight is Jacobian (1/volume)
182 }
183
184 // --- 2b. Apply the test filter to all required quantities ---
185 double u_filt = ApplyLESTestFilter(simCtx, u, weights);
186 double v_filt = ApplyLESTestFilter(simCtx, v, weights);
187 double w_filt = ApplyLESTestFilter(simCtx, w, weights);
188
189 double S11_filt = ApplyLESTestFilter(simCtx, S11, weights);
190 double S12_filt = ApplyLESTestFilter(simCtx, S12, weights);
191 double S13_filt = ApplyLESTestFilter(simCtx, S13, weights);
192 double S22_filt = ApplyLESTestFilter(simCtx, S22, weights);
193 double S23_filt = ApplyLESTestFilter(simCtx, S23, weights);
194 double S33_filt = ApplyLESTestFilter(simCtx, S33, weights);
195
196 double S_abs_filt = ApplyLESTestFilter(simCtx, S_abs_stencil, weights);
197
198 // Filter quadratic terms: <u*u>, <u*v>, etc.
199 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];
200 for(r=0; r<3; r++) for(q=0; q<3; q++) for(p=0; p<3; p++) {
201 uu[r][q][p] = u[r][q][p] * u[r][q][p];
202 uv[r][q][p] = u[r][q][p] * v[r][q][p];
203 uw[r][q][p] = u[r][q][p] * w[r][q][p];
204 vv[r][q][p] = v[r][q][p] * v[r][q][p];
205 vw[r][q][p] = v[r][q][p] * w[r][q][p];
206 ww[r][q][p] = w[r][q][p] * w[r][q][p];
207 }
208 double uu_filt = ApplyLESTestFilter(simCtx, uu, weights);
209 double uv_filt = ApplyLESTestFilter(simCtx, uv, weights);
210 double uw_filt = ApplyLESTestFilter(simCtx, uw, weights);
211 double vv_filt = ApplyLESTestFilter(simCtx, vv, weights);
212 double vw_filt = ApplyLESTestFilter(simCtx, vw, weights);
213 double ww_filt = ApplyLESTestFilter(simCtx, ww, weights);
214
215 // --- 2c. Compute the Leonard stress tensor (L_ij = <u_i*u_j> - <u_i>*<u_j>) ---
216 double L11 = uu_filt - u_filt * u_filt;
217 double L12 = uv_filt - u_filt * v_filt;
218 double L13 = uw_filt - u_filt * w_filt;
219 double L22 = vv_filt - v_filt * v_filt;
220 double L23 = vw_filt - v_filt * w_filt;
221 double L33 = ww_filt - w_filt * w_filt;
222
223 // --- 2d. Compute the model tensor M_ij ---
224 double grid_filter_width, test_filter_width;
225 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
226 grid_filter_width = pow(1.0 / aj[k][j][i], 1.0/3.0);
227 test_filter_width = 2.0 * grid_filter_width; // Standard box filter definition
228
229 double alpha = pow(test_filter_width / grid_filter_width, 2.0);
230
231 double M11 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S11_filt - S_abs_filt * S11_filt);
232 double M12 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S12_filt - S_abs_filt * S12_filt);
233 double M13 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S13_filt - S_abs_filt * S13_filt);
234 double M22 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S22_filt - S_abs_filt * S22_filt);
235 double M23 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S23_filt - S_abs_filt * S23_filt);
236 double M33 = -2.0 * grid_filter_width * grid_filter_width * (alpha * S_abs_filt * S33_filt - S_abs_filt * S33_filt);
237
238 // --- 2e. Contract tensors to find Cs^2 (L_ij * M_ij / (M_kl * M_kl)) ---
239 LM_arr[k][j][i] = L11*M11 + 2*L12*M12 + 2*L13*M13 + L22*M22 + 2*L23*M23 + L33*M33;
240 MM_arr[k][j][i] = M11*M11 + 2*M12*M12 + 2*M13*M13 + M22*M22 + 2*M23*M23 + M33*M33;
241 }
242
243 // --- 3. Averaging and finalization ---
244
245 // At this point, LM_arr and MM_arr contain raw, unaveraged values.
246 // To stabilize the solution, these are often averaged over homogeneous directions or locally.
247 // The logic for homogeneous averaging would go here if simCtx->i_homo_filter is true.
248 // For this general version, we proceed with the local (unaveraged) values.
249
250 for (k=lzs; k<lze; k++)
251 for (j=lys; j<lye; j++)
252 for (i=lxs; i<lxe; i++) {
253 if(nvert[k][j][i] > 1.1) {
254 Cs_arr[k][j][i] = 0.0;
255 continue;
256 }
257
258 double C_sq = 0.0;
259 if (MM_arr[k][j][i] > LES_EPSILON) {
260 C_sq = LM_arr[k][j][i] / MM_arr[k][j][i];
261 }
262
263 // Final clipping and assignment: Cs must be positive and not excessively large.
264 double Cs_val = (C_sq > 0.0) ? sqrt(C_sq) : 0.0;
265 Cs_arr[k][j][i] = PetscMin(PetscMax(Cs_val, 0.0), simCtx->max_cs);
266 }
267
268 // --- 4. Cleanup ---
269
270 // Restore all PETSc arrays
271 ierr = DMDAVecRestoreArray(fda, user->lUcat, &ucat); CHKERRQ(ierr);
272 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, (Cmpnts***)&csi); CHKERRQ(ierr);
273 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, (Cmpnts***)&eta); CHKERRQ(ierr);
274 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, (Cmpnts***)&zet); CHKERRQ(ierr);
275 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
276 ierr = DMDAVecRestoreArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
277 ierr = DMDAVecRestoreArray(da, user->CS, &Cs_arr); CHKERRQ(ierr);
278 ierr = DMDAVecRestoreArray(fda, lSx, &Sx_arr); CHKERRQ(ierr);
279 ierr = DMDAVecRestoreArray(fda, lSy, &Sy_arr); CHKERRQ(ierr);
280 ierr = DMDAVecRestoreArray(fda, lSz, &Sz_arr); CHKERRQ(ierr);
281 ierr = DMDAVecRestoreArray(da, lS_abs, &S_abs_arr); CHKERRQ(ierr);
282 ierr = DMDAVecRestoreArray(da, lLM, &LM_arr); CHKERRQ(ierr);
283 ierr = DMDAVecRestoreArray(da, lMM, &MM_arr); CHKERRQ(ierr);
284
285 // Destroy temporary vectors to prevent memory leaks
286 ierr = VecDestroy(&lSx); CHKERRQ(ierr);
287 ierr = VecDestroy(&lSy); CHKERRQ(ierr);
288 ierr = VecDestroy(&lSz); CHKERRQ(ierr);
289 ierr = VecDestroy(&lS_abs); CHKERRQ(ierr);
290 ierr = VecDestroy(&lLM); CHKERRQ(ierr);
291 ierr = VecDestroy(&lMM); CHKERRQ(ierr);
292
293 // Communicate ghost point data for the newly computed Cs field
294 ierr = DMGlobalToLocalBegin(da, user->CS, INSERT_VALUES, user->lCs); CHKERRQ(ierr);
295 ierr = DMGlobalToLocalEnd(da, user->CS, INSERT_VALUES, user->lCs); CHKERRQ(ierr);
296
297 PetscReal max_norm;
298 ierr = VecMax(user->CS, NULL, &max_norm); CHKERRQ(ierr);
299 LOG_ALLOW(GLOBAL, LOG_INFO, " Max dynamic Smagorinsky constant (Cs) computed: %e (clip value: %e)\n", max_norm, simCtx->max_cs);
300
301 PROFILE_FUNCTION_END;
302 PetscFunctionReturn(0);
303}
304
305
306#undef __FUNCT__
307#define __FUNCT__ "ComputeEddyViscosityLES"
308/**
309 * @brief Implementation of \ref ComputeEddyViscosityLES().
310 * @details Full API contract (arguments, ownership, side effects) is documented with
311 * the header declaration in `include/les.h`.
312 * @see ComputeEddyViscosityLES()
313 */
314
315PetscErrorCode ComputeEddyViscosityLES(UserCtx *user)
316{
317 PetscErrorCode ierr;
318 DM da = user->da, fda = user->fda;
319 DMDALocalInfo info;
320 PetscInt i, j, k;
321 PetscInt xs, xe, ys, ye, zs, ze;
322 PetscInt lxs, lxe, lys, lye, lzs, lze; // Local interior loop bounds
323 PetscInt mx, my, mz;
324
325 PetscFunctionBeginUser;
326 PROFILE_FUNCTION_BEGIN;
327
328 // Get grid dimensions and local bounds
329 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
330 mx = info.mx; my = info.my; mz = info.mz;
331 xs = info.xs; xe = xs + info.xm;
332 ys = info.ys; ye = ys + info.ym;
333 zs = info.zs; ze = zs + info.zm;
334
335 // Define loop bounds to exclude physical boundaries where stencils are invalid.
336 lxs = xs; lxe = xe; lys = ys; lye = ye; lzs = zs; lze = ze;
337 if (xs==0) lxs = xs+1;
338 if (ys==0) lys = ys+1;
339 if (zs==0) lzs = zs+1;
340 if (xe==mx) lxe = xe-1;
341 if (ye==my) lye = ye-1;
342 if (ze==mz) lze = ze-1;
343
344
345
346 // Get read/write access to PETSc data arrays
347 Cmpnts ***ucat;
348 PetscReal ***Cs_arr, ***nu_t_arr, ***nvert, ***aj;
349
350 ierr = DMDAVecGetArrayRead(fda, user->lUcat, (Cmpnts***)&ucat); CHKERRQ(ierr);
351 ierr = DMDAVecGetArrayRead(da, user->lCs, (PetscReal***)&Cs_arr); CHKERRQ(ierr);
352 ierr = DMDAVecGetArray(da, user->Nu_t, &nu_t_arr); CHKERRQ(ierr);
353 ierr = DMDAVecGetArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
354 ierr = DMDAVecGetArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
355
356 // Loop over the interior of the local domain
357 for (k=lzs; k<lze; k++)
358 for (j=lys; j<lye; j++)
359 for (i=lxs; i<lxe; i++) {
360 if(nvert[k][j][i] > 0.1) {
361 nu_t_arr[k][j][i] = 0.0;
362 continue;
363 }
364
365 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Computing eddy viscosity at point (%d,%d,%d)\n", i, j, k);
366 // 1. Compute the local strain rate magnitude |S|
367 Cmpnts dudx, dvdx, dwdx;
368 ierr = ComputeVectorFieldDerivatives(user, i, j, k, ucat, &dudx, &dvdx, &dwdx); CHKERRQ(ierr);
369
370 double Sxx = dudx.x;
371 double Sxy = 0.5 * (dudx.y + dvdx.x);
372 double Sxz = 0.5 * (dudx.z + dwdx.x);
373 double Syy = dvdx.y;
374 double Syz = 0.5 * (dvdx.z + dwdx.y);
375 double Szz = dwdx.z;
376 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) );
377
378 // 2. Determine the grid filter width, Delta = (cell volume)^(1/3)
379 double filter_width = pow( 1.0/aj[k][j][i], 1.0/3.0 );
380
381 // 3. Compute eddy viscosity: nu_t = (Cs * Delta)^2 * |S|
382 nu_t_arr[k][j][i] = pow(Cs_arr[k][j][i] * filter_width, 2.0) * strain_rate_mag;
383
384 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Cs=%.4e, Delta=%.4e, |S|=%.4e => nu_t=%.4e\n",
385 Cs_arr[k][j][i], filter_width, strain_rate_mag, nu_t_arr[k][j][i]);
386 }
387
388 // Restore PETSc data arrays
389 ierr = DMDAVecRestoreArrayRead(fda, user->lUcat, (Cmpnts***)&ucat); CHKERRQ(ierr);
390 ierr = DMDAVecRestoreArrayRead(da, user->lCs, (PetscReal***)&Cs_arr); CHKERRQ(ierr);
391 ierr = DMDAVecRestoreArray(da, user->Nu_t, &nu_t_arr); CHKERRQ(ierr);
392 ierr = DMDAVecRestoreArrayRead(da, user->lNvert, (PetscReal***)&nvert); CHKERRQ(ierr);
393 ierr = DMDAVecRestoreArrayRead(da, user->lAj, (PetscReal***)&aj); CHKERRQ(ierr);
394
395 // Update ghost points for the newly computed eddy viscosity
396 ierr = DMGlobalToLocalBegin(da, user->Nu_t, INSERT_VALUES, user->lNu_t); CHKERRQ(ierr);
397 ierr = DMGlobalToLocalEnd(da, user->Nu_t, INSERT_VALUES, user->lNu_t); CHKERRQ(ierr);
398
399 ierr = TransferPeriodicField(user,"Eddy Viscosity");CHKERRQ(ierr);
400
401 PetscReal max_norm;
402 ierr = VecMax(user->Nu_t, NULL, &max_norm); CHKERRQ(ierr);
403 LOG_ALLOW(GLOBAL, LOG_INFO, " Max eddy viscosity (Nu_t) computed: %e\n", max_norm);
404
405 PROFILE_FUNCTION_END;
406 PetscFunctionReturn(0);
407}
TransferPeriodicField
PetscErrorCode TransferPeriodicField(UserCtx *user, const char *field_name)
(Orchestrator) Applies periodic transfer for one field across all i/j/k directions.
Definition Boundaries.c:1656
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:123
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:283
ComputeEddyViscosityLES
PetscErrorCode ComputeEddyViscosityLES(UserCtx *user)
Implementation of ComputeEddyViscosityLES().
Definition les.c:315
ComputeSmagorinskyConstant
PetscErrorCode ComputeSmagorinskyConstant(UserCtx *user)
Implementation of ComputeSmagorinskyConstant().
Definition les.c:30
LES_EPSILON
const double LES_EPSILON
Definition les.c:18
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
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:199
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
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:770
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:2900
CONSTANT_SMAGORINSKY
@ CONSTANT_SMAGORINSKY
Definition variables.h:490
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
UserCtx::lZet
Vec lZet
Definition variables.h:858
UserCtx::lCs
Vec lCs
Definition variables.h:865
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:653
Cmpnts::x
PetscScalar x
Definition variables.h:101
SimCtx::max_cs
PetscReal max_cs
Definition variables.h:734
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:865
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:865
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
Cmpnts::z
PetscScalar z
Definition variables.h:101
SimCtx::Const_CS
PetscReal Const_CS
Definition variables.h:734
UserCtx::lUcont
Vec lUcont
Definition variables.h:837
SimCtx::step
PetscInt step
Definition variables.h:651
UserCtx::lAj
Vec lAj
Definition variables.h:858
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:858
SimCtx::les
PetscInt les
Definition variables.h:732
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:643
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:811
Generated by doxygen 1.9.8