setAlphaField crashes boundary condition related

Summary

setAlphaField crashes if used in combination with the inletOutlet boundary condition. Most likely also with others, but these are really common for alpha fields. This error does not exist in the v1812 and was hence introduced by some change. Tested in v2006

Steps to reproduce

Change one of the boundary conditions to inletOutlet

`Reading field alpha.water

--> FOAM FATAL ERROR:

request for surfaceScalarField phi from objectRegistry region0 failed
available objects of type surfaceScalarField are

1(mag(delta))

From const Type& Foam::objectRegistry::lookupObject(const Foam::word&, bool) const [with Type = Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>]
in file /home/bloerb/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 463.

FOAM aborting `

Many thanks @Goeke,

Could you please indicate which tutorial have you used for your test (so that we can first replicate the issue)?

assigned to @kuti

discInConstantFlow or one of the others with setAlphaFieldDict's inside

just replace some boundary condition with this:

    type            inletOutlet;
    inletValue      uniform 0;
    value           uniform 0;

@kuti @johan_roenby looking into it

In version of2006 and later, alpha1.correctBoundaryConditions() is called. inletOutlet requires phi to be set but only the alpha field is available in the object registry. setAlphaFields only updates the internalfield so this correctBoundaryConditions can be removed. This approach should have a lower risk than adding phi and U that might depend on additional variable.

The patch also removes a ununsed variable: 0001-BUG-setAlphaField-crashes-with-inletOutlet.patch

Please apply the patch to of2006 of2012 and ofdev

@johan_roenby What do you think?

I agree that correctBoundaryConditions() should be removed.

But there might be occasions where the value of alpha on boundary patches could depend on phi or U there.

Ideally, setAlphaField should read U first (if et exists) and set phi = U & Sf (or read it if it exists).

Then alpha1 can be set - first in internal cells (by the way, in the current implementation the forAll loops over the volScalarField alpha1, not the scalarField& alphaIn = alpha1.primitiveFieldRef()).

Next it should go through all the boundary pathces (like done here ) and set alpha to the isocut value for patches where fixedValue is false and to the fixed value specified by the BC (possibly depending on phi) if fixedValue is true.

Will you implement this, @Henning86 or should I?

on it

setAlphaField.C (updated version)

I refactored the code and the results are identical to the version from the patch

The boundary are set with the same approach as in this version

phi is only required because correctBoundaryConditions() was called and we can only set the value-field of a boundary condition with the the chosen approach

wall
{
        type            inletOutlet;
        inletValue      uniform 1;
        value           uniform 1; // set by setAlphaField but reset to inletValue after update 
}

So, adding U and phi would enable us to additionally modify the inletValue. But then we have to code it for every boundary condition

test-setAlpha.zip

I think you attached the old setAlphaField.C file, @Henning86 .

upload correct version in the old comment

Hi @Henning86 , @kuti and @Goeke

I now tested Henning's attached setAlphaField.C file and verified that 1) it compiles without errors, 2) runs without problems with the test-setAlpha.zip case above and 3) that the resulting alpha field looks right in paraview. So it is ready for a patch. I added minor style corrections and myself to the copyright notice: setAlphaField.C

Will you apply the fix, @kuti ?

Many thanks to three of you for your contributions.

@johan_roenby , could you please comment on the following, if possible?

facei and faceStatus are unused variables. Can we remove them?

            forAll(alphap, patchFacei)
            {
                const label facei = patchFacei + start;
                const label faceStatus = cutFace.calcSubFace(facei, 0.0);

                alphap[patchFacei] =
                    mag(cutFace.subFaceArea())/magSfp[patchFacei];
            }

            forAll(alphap, patchFacei)
            {
                const label facei = patchFacei + start;
                cutFace.calcSubFace(facei, 0.0); // removed faceStatus 

                alphap[patchFacei] =
                    mag(cutFace.subFaceArea())/magSfp[patchFacei];
            }

facei is used in calcSubFace not sure why you get an unused variable

@Henning86, The latest setAlphaField.C attachment contained faceStatus. Should I assume that I can remove it - or will you have another spin through with @johan_roenby ?

faceStatus can definitely be removed. facei not. Maybe the compiler thinks that since facei is used to calculate faceStatus which is not use, facei is not used either. I would guess, if you remove faceStatus, you also get rid of the warning that facei is not used.

The alpha field looks correct in my test. you can remove it

mentioned in commit de59afd7

mentioned in merge request !418 (merged)

mentioned in commit 54dfcf50

mentioned in commit 9f865914

closed