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