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 42 of file les.c.

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}
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
FinalizeSmagorinskyConstantField
static PetscErrorCode FinalizeSmagorinskyConstantField(UserCtx *user)
Synchronizes the completed Smagorinsky coefficient field.
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: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
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
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
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::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
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 327 of file les.c.

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
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1755
UserCtx::lCs
Vec lCs
Definition variables.h:935
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:935
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8