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])
	(Internal) Applies a 2D homogeneous test filter in the i-k plane using Simpson's rule.
static double	ApplyVolumeWeightedBoxFilter (double values[3][3][3], double weights[3][3][3])
	(Internal) Applies a volume-weighted 3D box filter to a 3x3x3 stencil.
double	ApplyLESTestFilter (const SimCtx *simCtx, double values[3][3][3], double weights[3][3][3])
	Public interface for applying the LES test filter.

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 141 of file Filter.c.

Function Documentation

ApplySimpsonRuleHomogeneousFilter()

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

static

(Internal) Applies a 2D homogeneous test filter in the i-k plane using Simpson's rule.

This is a specialized, high-accuracy filter designed for flow cases that are statistically homogeneous in two directions (e.g., the streamwise and spanwise directions of a channel flow). It uses a fixed-weight stencil based on Simpson's rule for numerical integration, which offers better spectral properties than a simple box filter. The weights are constant and do not depend on cell volume, assuming a uniform grid in the homogeneous plane.

Morinishi, Kera, & Vasilyev, "Skew-symmetric formulation for large eddy simulation of compressible turbulent flows," Journal of Computational Physics, 2004.

Parameters

values

The 3x3x3 array of scalar values from the local stencil. For this 2D filter, only the central j-plane (j-index = 1) is utilized.

Returns

The filtered scalar value.

Definition at line 34 of file Filter.c.

{
    // The stencil only uses the central j-plane (j-index = 1, corresponding to the y-direction).
    // The formula is a weighted sum of the 9 points in that plane.
    const double corners = values[0][1][0] + values[2][1][0] + values[0][1][2] + values[2][1][2]; // 4 corner points
    const double edges   = values[0][1][1] + values[1][1][0] + values[2][1][1] + values[1][1][2]; // 4 edge-center points
    const double center  = values[1][1][1];                                                        // 1 center point
 
    // The weights (1, 4, 16) and normalization factor (36) come from the 2D Simpson's rule.
    return (corners + 4.0 * edges + 16.0 * center) / 36.0;
}

Here is the caller graph for this function:

ApplyVolumeWeightedBoxFilter()

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

static

(Internal) Applies a volume-weighted 3D box filter to a 3x3x3 stencil.

This is the general-purpose filter for non-uniform, curvilinear grids. To correctly average a quantity over a space with varying cell sizes, the contribution of each point must be weighted by its associated volume (or, equivalently, 1/Jacobian).

The logic calculates a weighted average over eight 2x2x2 sub-cubes (or "octants") that surround the central node of the 3x3x3 stencil. This produces a filtered value representative of the larger volume covered by the test filter.

Parameters

values	The 3x3x3 array of scalar values at the stencil points.
weights	The 3x3x3 array of weights, where weight = 1.0 / cell_volume (i.e., the Jacobian).

Returns

The volume-weighted filtered scalar value.

Definition at line 62 of file Filter.c.

{
    // v1...v8 store the sum of (value * weight) for each of the 8 sub-cubes.
    // w1...w8 store the sum of (weight) for each of the 8 sub-cubes.
    double v1, v2, v3, v4, v5, v6, v7, v8;
    double w1, w2, w3, w4, w5, w6, w7, w8;
 
    // --- Calculations for the 4 sub-cubes on the bottom layer (k-indices 0 and 1) ---
 
    // Bottom-Back-Left sub-cube (i-indices: 0,1; j-indices: 0,1; k-indices: 0,1)
    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] +
           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] );
    w1 = ( weights[0][0][0] + weights[1][0][0] + weights[0][1][0] + weights[1][1][0] +
           weights[0][0][1] + weights[1][0][1] + weights[0][1][1] + weights[1][1][1] );
 
    // Bottom-Back-Right sub-cube (i-indices: 1,2; j-indices: 0,1; k-indices: 0,1)
    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] +
           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] );
    w2 = ( weights[1][0][0] + weights[2][0][0] + weights[1][1][0] + weights[2][1][0] +
           weights[1][0][1] + weights[2][0][1] + weights[1][1][1] + weights[2][1][1] );
 
    // Bottom-Front-Left sub-cube (i-indices: 0,1; j-indices: 1,2; k-indices: 0,1)
    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] +
           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] );
    w3 = ( weights[0][1][0] + weights[1][1][0] + weights[0][2][0] + weights[1][2][0] +
           weights[0][1][1] + weights[1][1][1] + weights[0][2][1] + weights[1][2][1] );
 
    // Bottom-Front-Right sub-cube (i-indices: 1,2; j-indices: 1,2; k-indices: 0,1)
    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] +
           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] );
    w4 = ( weights[1][1][0] + weights[2][1][0] + weights[1][2][0] + weights[2][2][0] +
           weights[1][1][1] + weights[2][1][1] + weights[1][2][1] + weights[2][2][1] );
 
 
    // --- Calculations for the 4 sub-cubes on the top layer (k-indices 1 and 2) ---
 
    // Top-Back-Left sub-cube (i-indices: 0,1; j-indices: 0,1; k-indices: 1,2)
    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] +
           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] );
    w5 = ( weights[0][0][1] + weights[1][0][1] + weights[0][1][1] + weights[1][1][1] +
           weights[0][0][2] + weights[1][0][2] + weights[0][1][2] + weights[1][1][2] );
 
    // Top-Back-Right sub-cube (i-indices: 1,2; j-indices: 0,1; k-indices: 1,2)
    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] +
           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] );
    w6 = ( weights[1][0][1] + weights[2][0][1] + weights[1][1][1] + weights[2][1][1] +
           weights[1][0][2] + weights[2][0][2] + weights[1][1][2] + weights[2][1][2] );
 
    // Top-Front-Left sub-cube (i-indices: 0,1; j-indices: 1,2; k-indices: 1,2)
    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] +
           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] );
    w7 = ( weights[0][1][1] + weights[1][1][1] + weights[0][2][1] + weights[1][2][1] +
           weights[0][1][2] + weights[1][1][2] + weights[0][2][2] + weights[1][2][2] );
 
    // Top-Front-Right sub-cube (i-indices: 1,2; j-indices: 1,2; k-indices: 1,2)
    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] +
           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] );
    w8 = ( weights[1][1][1] + weights[2][1][1] + weights[1][2][1] + weights[2][2][1] +
           weights[1][1][2] + weights[2][1][2] + weights[1][2][2] + weights[2][2][2] );
 
    // Sum the contributions from all 8 octants.
    double total_weighted_value = v1+v2+v3+v4+v5+v6+v7+v8;
    double total_weight         = w1+w2+w3+w4+w5+w6+w7+w8;
 
    // Production safety check: avoid division by zero if all weights are zero
    // (e.g., if the stencil is entirely inside a solid body).
    if (fabs(total_weight) < 1.0e-12) {
        return 0.0;
    }
 
    return total_weighted_value / total_weight;
}

Here is the caller graph for this function:

ApplyLESTestFilter()

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

Public interface for applying the LES test filter.

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

Definition at line 145 of file Filter.c.

{
    // This function acts as a dispatcher. It reads the simulation configuration
    // and calls the appropriate internal filtering routine.
    if (simCtx->testfilter_ik) {
        // If the "-testfilter_ik" flag is set, it indicates that the simulation
        // is homogeneous in the i and k directions, so we use the specialized,
        // more accurate Simpson's rule filter. The `weights` are ignored in this case.
        return ApplySimpsonRuleHomogeneousFilter(values);
    } else {
        // This is the default case for general, non-uniform, curvilinear grids.
        // The volume-weighted box filter is used to ensure correct averaging.
        return ApplyVolumeWeightedBoxFilter(values, weights);
    }
}

Here is the call graph for this function:

Here is the caller graph for this function: