PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
Filter.c File Reference

Implements numerical filtering schemes for Large Eddy Simulation (LES). More...

#include "Filter.h"
#include <math.h>
Include dependency graph for Filter.c:

Go to the source code of this file.

Macros

#define __FUNCT__   "ApplyLESTestFilter"
 

Functions

static double ApplySimpsonRuleHomogeneousFilter (double values[3][3][3])
 Apply the two-dimensional Simpson stencil on the central homogeneous-plane samples.
 
static double ApplyVolumeWeightedBoxFilter (double values[3][3][3], double weights[3][3][3])
 Average a 3x3x3 stencil with cell-volume weights while excluding solid cells.
 
double ApplyLESTestFilter (const SimCtx *simCtx, double values[3][3][3], double weights[3][3][3])
 Internal helper implementation: ApplyLESTestFilter().
 

Detailed Description

Implements numerical filtering schemes for Large Eddy Simulation (LES).

This file contains the functions necessary to apply a "test filter" to the resolved velocity field, which is a core component of the dynamic Smagorinsky turbulence model. The choice of filter (e.g., a general 3D box filter or a specialized 2D homogeneous filter) is determined by the simulation's configuration.

Definition in file Filter.c.

Macro Definition Documentation

◆ __FUNCT__

#define __FUNCT__   "ApplyLESTestFilter"

Definition at line 116 of file Filter.c.

Function Documentation

◆ ApplySimpsonRuleHomogeneousFilter()

static double ApplySimpsonRuleHomogeneousFilter ( double  values[3][3][3])
static

Apply the two-dimensional Simpson stencil on the central homogeneous-plane samples.

Definition at line 21 of file Filter.c.

22{
23 // The stencil only uses the central j-plane (j-index = 1, corresponding to the y-direction).
24 // The formula is a weighted sum of the 9 points in that plane.
25 const double corners = values[0][1][0] + values[2][1][0] + values[0][1][2] + values[2][1][2]; // 4 corner points
26 const double edges = values[0][1][1] + values[1][1][0] + values[2][1][1] + values[1][1][2]; // 4 edge-center points
27 const double center = values[1][1][1]; // 1 center point
28
29 // The weights (1, 4, 16) and normalization factor (36) come from the 2D Simpson's rule.
30 return (corners + 4.0 * edges + 16.0 * center) / 36.0;
31}
Here is the caller graph for this function:

◆ ApplyVolumeWeightedBoxFilter()

static double ApplyVolumeWeightedBoxFilter ( double  values[3][3][3],
double  weights[3][3][3] 
)
static

Average a 3x3x3 stencil with cell-volume weights while excluding solid cells.

Definition at line 37 of file Filter.c.

38{
39 // v1...v8 store the sum of (value * weight) for each of the 8 sub-cubes.
40 // w1...w8 store the sum of (weight) for each of the 8 sub-cubes.
41 double v1, v2, v3, v4, v5, v6, v7, v8;
42 double w1, w2, w3, w4, w5, w6, w7, w8;
43
44 // --- Calculations for the 4 sub-cubes on the bottom layer (k-indices 0 and 1) ---
45
46 // Bottom-Back-Left sub-cube (i-indices: 0,1; j-indices: 0,1; k-indices: 0,1)
47 v1 = ( values[0][0][0]*weights[0][0][0] + values[1][0][0]*weights[1][0][0] + values[0][1][0]*weights[0][1][0] + values[1][1][0]*weights[1][1][0] +
48 values[0][0][1]*weights[0][0][1] + values[1][0][1]*weights[1][0][1] + values[0][1][1]*weights[0][1][1] + values[1][1][1]*weights[1][1][1] );
49 w1 = ( weights[0][0][0] + weights[1][0][0] + weights[0][1][0] + weights[1][1][0] +
50 weights[0][0][1] + weights[1][0][1] + weights[0][1][1] + weights[1][1][1] );
51
52 // Bottom-Back-Right sub-cube (i-indices: 1,2; j-indices: 0,1; k-indices: 0,1)
53 v2 = ( values[1][0][0]*weights[1][0][0] + values[2][0][0]*weights[2][0][0] + values[1][1][0]*weights[1][1][0] + values[2][1][0]*weights[2][1][0] +
54 values[1][0][1]*weights[1][0][1] + values[2][0][1]*weights[2][0][1] + values[1][1][1]*weights[1][1][1] + values[2][1][1]*weights[2][1][1] );
55 w2 = ( weights[1][0][0] + weights[2][0][0] + weights[1][1][0] + weights[2][1][0] +
56 weights[1][0][1] + weights[2][0][1] + weights[1][1][1] + weights[2][1][1] );
57
58 // Bottom-Front-Left sub-cube (i-indices: 0,1; j-indices: 1,2; k-indices: 0,1)
59 v3 = ( values[0][1][0]*weights[0][1][0] + values[1][1][0]*weights[1][1][0] + values[0][2][0]*weights[0][2][0] + values[1][2][0]*weights[1][2][0] +
60 values[0][1][1]*weights[0][1][1] + values[1][1][1]*weights[1][1][1] + values[0][2][1]*weights[0][2][1] + values[1][2][1]*weights[1][2][1] );
61 w3 = ( weights[0][1][0] + weights[1][1][0] + weights[0][2][0] + weights[1][2][0] +
62 weights[0][1][1] + weights[1][1][1] + weights[0][2][1] + weights[1][2][1] );
63
64 // Bottom-Front-Right sub-cube (i-indices: 1,2; j-indices: 1,2; k-indices: 0,1)
65 v4 = ( values[1][1][0]*weights[1][1][0] + values[2][1][0]*weights[2][1][0] + values[1][2][0]*weights[1][2][0] + values[2][2][0]*weights[2][2][0] +
66 values[1][1][1]*weights[1][1][1] + values[2][1][1]*weights[2][1][1] + values[1][2][1]*weights[1][2][1] + values[2][2][1]*weights[2][2][1] );
67 w4 = ( weights[1][1][0] + weights[2][1][0] + weights[1][2][0] + weights[2][2][0] +
68 weights[1][1][1] + weights[2][1][1] + weights[1][2][1] + weights[2][2][1] );
69
70
71 // --- Calculations for the 4 sub-cubes on the top layer (k-indices 1 and 2) ---
72
73 // Top-Back-Left sub-cube (i-indices: 0,1; j-indices: 0,1; k-indices: 1,2)
74 v5 = ( values[0][0][1]*weights[0][0][1] + values[1][0][1]*weights[1][0][1] + values[0][1][1]*weights[0][1][1] + values[1][1][1]*weights[1][1][1] +
75 values[0][0][2]*weights[0][0][2] + values[1][0][2]*weights[1][0][2] + values[0][1][2]*weights[0][1][2] + values[1][1][2]*weights[1][1][2] );
76 w5 = ( weights[0][0][1] + weights[1][0][1] + weights[0][1][1] + weights[1][1][1] +
77 weights[0][0][2] + weights[1][0][2] + weights[0][1][2] + weights[1][1][2] );
78
79 // Top-Back-Right sub-cube (i-indices: 1,2; j-indices: 0,1; k-indices: 1,2)
80 v6 = ( values[1][0][1]*weights[1][0][1] + values[2][0][1]*weights[2][0][1] + values[1][1][1]*weights[1][1][1] + values[2][1][1]*weights[2][1][1] +
81 values[1][0][2]*weights[1][0][2] + values[2][0][2]*weights[2][0][2] + values[1][1][2]*weights[1][1][2] + values[2][1][2]*weights[2][1][2] );
82 w6 = ( weights[1][0][1] + weights[2][0][1] + weights[1][1][1] + weights[2][1][1] +
83 weights[1][0][2] + weights[2][0][2] + weights[1][1][2] + weights[2][1][2] );
84
85 // Top-Front-Left sub-cube (i-indices: 0,1; j-indices: 1,2; k-indices: 1,2)
86 v7 = ( values[0][1][1]*weights[0][1][1] + values[1][1][1]*weights[1][1][1] + values[0][2][1]*weights[0][2][1] + values[1][2][1]*weights[1][2][1] +
87 values[0][1][2]*weights[0][1][2] + values[1][1][2]*weights[1][1][2] + values[0][2][2]*weights[0][2][2] + values[1][2][2]*weights[1][2][2] );
88 w7 = ( weights[0][1][1] + weights[1][1][1] + weights[0][2][1] + weights[1][2][1] +
89 weights[0][1][2] + weights[1][1][2] + weights[0][2][2] + weights[1][2][2] );
90
91 // Top-Front-Right sub-cube (i-indices: 1,2; j-indices: 1,2; k-indices: 1,2)
92 v8 = ( values[1][1][1]*weights[1][1][1] + values[2][1][1]*weights[2][1][1] + values[1][2][1]*weights[1][2][1] + values[2][2][1]*weights[2][2][1] +
93 values[1][1][2]*weights[1][1][2] + values[2][1][2]*weights[2][1][2] + values[1][2][2]*weights[1][2][2] + values[2][2][2]*weights[2][2][2] );
94 w8 = ( weights[1][1][1] + weights[2][1][1] + weights[1][2][1] + weights[2][2][1] +
95 weights[1][1][2] + weights[2][1][2] + weights[1][2][2] + weights[2][2][2] );
96
97 // Sum the contributions from all 8 octants.
98 double total_weighted_value = v1+v2+v3+v4+v5+v6+v7+v8;
99 double total_weight = w1+w2+w3+w4+w5+w6+w7+w8;
100
101 // Production safety check: avoid division by zero if all weights are zero
102 // (e.g., if the stencil is entirely inside a solid body).
103 if (fabs(total_weight) < 1.0e-12) {
104 return 0.0;
105 }
106
107 return total_weighted_value / total_weight;
108}
Here is the caller graph for this function:

◆ ApplyLESTestFilter()

double ApplyLESTestFilter ( const SimCtx *  simCtx,
double  values[3][3][3],
double  weights[3][3][3] 
)

Internal helper implementation: ApplyLESTestFilter().

Applies a numerical "test filter" to a 3x3x3 stencil of data points.

Local to this translation unit.

Definition at line 121 of file Filter.c.

122{
123 // This function acts as a dispatcher. It reads the simulation configuration
124 // and calls the appropriate internal filtering routine.
125 if (simCtx->testfilter_ik) {
126 // If the "-testfilter_ik" flag is set, it indicates that the simulation
127 // is homogeneous in the i and k directions, so we use the specialized,
128 // more accurate Simpson's rule filter. The `weights` are ignored in this case.
129 return ApplySimpsonRuleHomogeneousFilter(values);
130 } else {
131 // This is the default case for general, non-uniform, curvilinear grids.
132 // The volume-weighted box filter is used to ensure correct averaging.
133 return ApplyVolumeWeightedBoxFilter(values, weights);
134 }
135}
ApplyVolumeWeightedBoxFilter
static double ApplyVolumeWeightedBoxFilter(double values[3][3][3], double weights[3][3][3])
Average a 3x3x3 stencil with cell-volume weights while excluding solid cells.
Definition Filter.c:37
ApplySimpsonRuleHomogeneousFilter
static double ApplySimpsonRuleHomogeneousFilter(double values[3][3][3])
Apply the two-dimensional Simpson stencil on the central homogeneous-plane samples.
Definition Filter.c:21
SimCtx::testfilter_ik
PetscInt testfilter_ik
Definition variables.h:824
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8