/*--------------------------------*- C++ -*----------------------------------*\
  ==  == ====== ====   ====    |
                    \\     ||  | Multiphase Code Repository by HZDR
  ======   //   ||  || ===//   | Website: https://doi.org/10.14278/rodare.767
  ||  ||  //    ||  // || \\   | License: GPL-3.0-or-later
  ==  == ====== ====   ==  ==  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      generateAlphas;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

type            coded;

libs            ("libutilityFunctionObjects.so");

name            generateAlphas;

codeWrite
#{
    const dimensionedVector dy(dimless/dimLength, vector(0, 1, 0));

    const volScalarField y(mesh().C() & dy);

    const scalar centre=0.01;

    volScalarField alphaAir
    (
        IOobject
        (
            IOobject::groupName("alpha", "air"),
            mesh().time().name(),
            mesh(),
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh(),
        dimensionedScalar(dimless, 1)
    );

    volScalarField alphaWater
    (
        IOobject
        (
            IOobject::groupName("alpha", "water"),
            mesh().time().name(),
            mesh(),
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh(),
        dimensionedScalar(dimless, 0)
    );

    // linear
    // const scalar delta=0.006, y0=centre-0.5*delta, y1=centre+0.5*delta;
    // const volScalarField f
    // (
    //                 neg0(y - y0)*(0                 )
    //   + pos(y - y0)*neg0(y - y1)*((y - y0)/(y1 - y0))
    //   + pos(y - y1)             *(1                 )
    // );

    // tanh
    // const scalar gradAlpha=800;
    // const volScalarField f
    // (
    //     0.5*(tanh(gradAlpha*(y - centre)) + 1.0)
    // );

    // polynomial
    const scalar gradAlpha=600, y0=-1.5/gradAlpha, y1=1.5/gradAlpha;
    const volScalarField yTrafo(y - centre);
    const volScalarField g
    (
                         neg0(yTrafo - y0)*(-1)
      - pos(yTrafo - y0)*neg0(yTrafo - y1)*(4.0/27.0*pow3(gradAlpha*yTrafo)
                                            - gradAlpha*yTrafo)
      + pos(yTrafo - y1)                  *(1)
    );
    const volScalarField f(0.5*(g + 1));

    alphaAir = f;
    alphaWater = 1 - f;

    alphaAir.write();
    alphaWater.write();
#};

// ************************************************************************* //
