Add relaxation factors

[HendriceH]

Motivation

If we get performance improvements from kernel fusion or other stuff things. They will be more pronounced in a PIMPLE setup which solves the momentum and turbulence equation more often. The only missing piece for a neoPimpleFoam implementation are relaxation factors. This PR does a draft implementation of the functionality quickly added to dsl/solver.hpp
A neoPimpleFoam solver is implemented in https://github.com/exasim-project/NeoFOAM/pull/204
Status:

• Something is still off as the solver diverges, probably with the matrix relaxation.
• Is internalCoeffs_[patchi] = NeoN auxillary Boundary coeffs? I might need to use the internalCoeffs of the boundarycells?

Reference

The reference concept is the following

1. Calculate the sum-mag off-diagonal from the interior faces
2. Add the cmptMax(cmptMag() of the internalCoeffs of the boundary faces to the diagonal
3. Make the matrix diagonal dominant
4. Apply relaxation to the diagonal
5. remove the boundary contribution
6. add to the source term

template<class Type>
void Foam::fvMatrix<Type>::relax(const scalar alpha)
{
    if (alpha <= 0)
    {
        return;
    }

    DebugInFunction
        << "Relaxing " << psi_.name() << " by " << alpha << endl;

    Field<Type>& S = source();
    scalarField& D = diag();

    // Store the current unrelaxed diagonal for use in updating the source
    scalarField D0(D);

    // Calculate the sum-mag off-diagonal from the interior faces
    scalarField sumOff(D.size(), Zero);
    sumMagOffDiag(sumOff);

    // Handle the boundary contributions to the diagonal
    forAll(psi_.boundaryField(), patchi)
    {
        const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];

        if (ptf.size())
        {
            const labelUList& pa = lduAddr().patchAddr(patchi);
            Field<Type>& iCoeffs = internalCoeffs_[patchi];

            if (ptf.coupled())
            {
                const Field<Type>& pCoeffs = boundaryCoeffs_[patchi];

                // For coupled boundaries add the diagonal and
                // off-diagonal contributions
                forAll(pa, face)
                {
                    D[pa[face]] += component(iCoeffs[face], 0);
                    sumOff[pa[face]] += mag(component(pCoeffs[face], 0));
                }
            }
            else
            {
                // For non-coupled boundaries add the maximum magnitude diagonal
                // contribution to ensure stability
                forAll(pa, face)
                {
                    D[pa[face]] += cmptMax(cmptMag(iCoeffs[face]));
                }
            }
        }
    }

    // Ensure the matrix is diagonally dominant...
    // Assumes that the central coefficient is positive and ensures it is
    forAll(D, celli)
    {
        D[celli] = max(mag(D[celli]), sumOff[celli]);
    }

    // ... then relax
    D /= alpha;

    // Now remove the diagonal contribution from coupled boundaries
    forAll(psi_.boundaryField(), patchi)
    {
        const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];

        if (ptf.size())
        {
            const labelUList& pa = lduAddr().patchAddr(patchi);
            Field<Type>& iCoeffs = internalCoeffs_[patchi];

            if (ptf.coupled())
            {
                forAll(pa, face)
                {
                    D[pa[face]] -= component(iCoeffs[face], 0);
                }
            }
            else
            {
                forAll(pa, face)
                {
                    D[pa[face]] -= cmptMin(iCoeffs[face]);
                }
            }
        }
    }

    // Finally add the relaxation contribution to the source.
    S += (D - D0)*psi_.primitiveField();
}

with the helper:

template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::sumMagOffDiag
(
    Field<LUType>& sumOff
) const
{
    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();

    const labelUList& l = lduAddr().lowerAddr();
    const labelUList& u = lduAddr().upperAddr();

    for (label face = 0; face < l.size(); face++)
    {
        sumOff[u[face]] += cmptMag(Lower[face]);
        sumOff[l[face]] += cmptMag(Upper[face]);
    }
}

[github-actions]

Thank you for your PR, here are some useful tips:

• If this is your first PR, please add yourself to the AUTHORS.md file
• Use labels to change the workflow behavior. Set auto-fix to automatically apply fixes, full-ci for a comprehensive set of tests and sanitizer runs, ready-for-review to show that you're done and need a review.
• I deployed a test documentation to https://exasim-project.com/NeoN/Build_PR_436
• For full-ci builds benchmark results are available at https://github.com/exasim-project/NeoFOAM-BenchmarkData/436
• Also check if the GPU pipelines on https://gitlab-ce.lrz.de/greole/neon/-/jobs are successful.

[greole]

Needs a docstring what ithe purpose is.

[greole]

I think we already have this?

[greole]

Needs a docstring what this function does.

[greole]

Needs a docstring.

[greole]

Can we pull this out into a separate PR?