v2412: New and improved numerics

New fused discretisation

TOP

New fused discretisation methods avoid intermediate interpolation steps and fuse with the face-based integration. The functionality is equivalent to the usual discretisation:

  • gradSchemes : fusedGauss
  • divSchemes : fusedGauss
  • laplacianSchemes : fusedGauss, fusedLeastSquares

These schemes avoid memory operations, benefiting cases with large numbers of patches or on nodes with slow memory access, e.g. GPU accessing CPU data. Comparing execution time vs. simulation time for a large rotating wheel assembly case where 'logs-gauss' show the default behaviour and 'logs-fused' for the drop-in replacement fused schemes:

Note that this is a preliminary work and only considers explicit operations, and the order of operations is different leading to different truncation errors.

Tutorials

Source code

Merge request

New GAMG agglomeration caching

TOP

Moving mesh cases automatically rebuild the agglomeration for every time step (or even every pressure corrector). For most mesh motion cases, however, the GAMG agglomeration table can be kept for some iterations before needing to be recreated. For solid-body motion there is no need to update the table since the agglomeration logic only depends on the matrix coefficients.

Caching is enabled via the new updateInterval entry. This can be a very large number for solid-body simulations; for general mesh motion this should be based on the expected mesh Courant number.

p
{
    solver              GAMG;
    ..
    cacheAgglomeration  yes;
    updateInterval      100000;
}

For verbose feedback on the agglomeration, activate the debug switches in the system/controlDict

DebugSwitches
{
    // Print agglomeration
    GAMGAgglomeration   1;
}

faceAreaWeight optimisation

The default agglomeration method employs a pair-wise combination of 'cells' based on the magnitude of the face weight connecting the cells. Currently this uses a face weight derived from face area vectors, perturbed to avoid jitter in the agglomeration boundary on an axis-aligned mesh.

When the updateInterval is not 1, the weights are calculated using unscaled, unperturbed face area magnitudes to be more consistent with Gaussian-type discretisation. This can yield a slight improvement in the number of GAMG cycles, e.g. the incompressible/pisoFoam/LES/motorBike tutorial using 10 iterations of simpleFoam showed an average reduction in cycles of around 10%:

Note that algebraic agglomeration can run with updateInterval; but it

  • does not update cyclicA(C)MI and will do a bad job of the interpolation since it uses an old stencil at the coarser levels
  • is based on the matrix coefficients which are changing during the simulation so would not be constant even during solid-body motion

Tutorials

Source code

Merge request

Improved lduMatrix, lduAddressing

TOP

The scalar matrix class lduMatrix has been extended to make it easier to work with cell-based algorithms.

The addressing lduAddressing includes a lowerCSRAddr accessor that returns the lower addressing, avoiding the additional indirection using loSortAddressing:

const labelUList& lduAddressing::lowerCSRAddr() const

This is constructed on-the-fly from the lowerAddr() addressing.

Additionally, the matrix has been extended with a lowerCSR() accessor that returns the corresponding reordering of the lower() coefficients, or can be constructed directly:

const scalarField& lduMatrix::lowerCSR() const

Availability of the addressing can be checked using the hasLowerCSR() function.

An example usage is given inside the lduMatrix::Amul routine:

if (hasLowerCSR())
{
    // Note: lowerCSR constructed from lower if available, upper otherwise
    //       so is handling symmetric()
    const scalar* const __restrict__ lowercsrPtr = lowerCSR().begin();

    for (label cell=0; cell<nCells; cell++)
    {
        scalar& val = ApsiPtr[cell];

        val = diagPtr[cell]*psiPtr[cell];

        // Add lower contributions
        {
            const label start = loStartPtr[cell];
            const label end = loStartPtr[cell+1];

            for (label i = start; i < end; i++)
            {
                const label nbrCell = lcsrPtr[i];
                val += lowercsrPtr[i]*psiPtr[nbrCell];
            }
        }
        // Add upper contributions
        {
            const label start = oStartPtr[cell];
            const label end = oStartPtr[cell+1];

            for (label i = start; i < end; i++)
            {
                const label nbrCell = uPtr[i];
                val += upperPtr[i]*psiPtr[nbrCell];
            }
        }
    }
}

Source code

Improved wall distance

TOP

The near-wall distance field used in e.g. turbulence wall functions, has been updated for wall faces, changing the behaviour of cells with faces on multiple wall patches.

The nearest distance is calculated as the minimum distance from the cell centre to the local wall face or any local edge- or point-connected face. Note that this

  • is not parallel consistent
  • only considers one face

For more consistent behaviour the exact wall distance method should be employed.

Previous behaviour (OpenFOAM v2406 and earlier) can be obtained by deactivating the useCombinedWallPatch info switch in etc/controlDict or the local system/controlDict:


OptimisationSwitches
{
    useCombinedWallPatch    0;
}

With a skewed/tetified version of the cavity mesh, v2406 wall distance:

The new behaviour looks across the side patches as well as the top patch:

Source code

Merge request

Issue

Community contribution: Improved adjoint optimisation

TOP

OpenFOAM v2412 adds some quality-of-life changes to the adjoint library.

The initial values of the design variables in topology optimisation can now be read through the 0/alpha field, which can be easily set using the topoSet/setFields utilities. Previously, this was only possible through 0/uniform/topOVars, which is harder to manipulate.

The multiplier of the term added to the mathematical optimisation problem of ISQP and MMA to guarantee feasibility (named 'c') can now be given by a Function1, to allow its manipulation throughout the optimisation loop. A typical example would be setting it to a small value in the first few optimisation cycles, to focus on optimality rather than feasibility, and gradually increase it.

Added the option to disable damping of the approximate Hessian in ISQP. In the general case, this is to be avoided but can accelerate the algorithm in some fringe cases.

Tutorials

Source code

Merge request

New zoneBlended scheme

TOP

The new zoneBlended scheme enables users to apply differencing schemes per-face-zone. Schemes are set in dictionary format according to:

divSchemes
{
    .
    .
    div(phi,U)      Gauss zoneBlended
    {
        default         defaultScheme;
        faceZone1       scheme1;
        faceZone2       scheme2;
        ...
        faceZoneN       schemeN;
    }
    .
    .
}

The default entry specifies the background scheme; additional schemes can be set per faceZone, e.g. scheme1 is applied to faceZone1, scheme2 is applied to faceZone2 etc.

Source code