[ChrisA]

I'm running a modified rhoCentralFoam which is currently based on openfoam 2.2.1 and has been modified to run with binary gas mixtures. The issue I'm having is with a modified version of the maxwellSlip boundary condition which I've changed around to use a different mean free path definition. The issue is when I try to run in parallel the code crashes with a floating point exception. It seems to run fine on a single core.

The behavior is very similar to that explained here but since it should have been fixed in this version and I have thermalCreep off (commented out in this iteration of code JUST to be sure) that can't be the issue. The standard maxwellSlip boundary condition works fine as well.

The modifications made to the maxwellSlip boundary are as follows:
In the part where the valueFraction is defined, the code has been modified to accept a lambda value which is passed from the main solver to the boundary condition:

    const fvPatchField<scalar>& plambda =
        patch().lookupPatchField<volScalarField, scalar>(lambdaName_);

    Field<scalar> C1
    (
        (2.0 - accommodationCoeff_)/accommodationCoeff_
    );

    Field<scalar> pnu(pmu/prho);
    valueFraction() = (1.0/(1.0 + patch().deltaCoeffs()*C1*plambda));

    refValue() = Uwall_;

lambda is the mean free path and is defined in the rhoCentralFoam code because it depends on quite a few variables not immediately/easily accessible by the boundary condition. Lambda is declared as follows:

volScalarField lambda
(
        IOobject
        (
             "lambda",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
    mesh,
    dimensionedScalar("lambda", dimensionSet(0,1,0,0,0), 0.0)
);

and given a value with the following:

lambda = 1/(3.14*(n1*pow(dsi,2)*Foam::sqrt(1+(m1/m2))+n2*pow(dsi,2)*Foam::sqrt(1+(m2/m1)))); //mean free path

        Info << "lambda = " << max(lambda) << endl;


forAll(lambda.boundaryField(),patchi)
{
     lambda.boundaryField()[patchi] == lambda.boundaryField()[patchi].patchInternalField();
}
    lambda.correctBoundaryConditions();

I'm not sure if the forAll loop and the line after it are required, but I've had issues with it before so it's there as one of the things I've been trying. I'm stumped on this one. The definitions don't rely on any adjacent cells and the only thing that's changed is the mean free path. I should also note that when I had a mistake in my code and the valueFraction was essentially hard coded to a very small number (1e-6).

I'm currently running the simulation on a flat plate.

Edit: turns out it crashes in serial too, it just progresses further in the simulation before this happens. It could be possible I made a mistake when I interpreted the valueFraction expected value and the place where mean free path pops up? The code diverges when it tries to solve for Ux (tangential plate velocity).