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 "Boundaries.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.

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

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

See also
ComputeSmagorinskyConstant()

Definition at line 30 of file les.c.

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}
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
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
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::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
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.

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

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

See also
ComputeEddyViscosityLES()

Definition at line 315 of file les.c.

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
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:865
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:865
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8