#include "fvCFD.H" // This file brings in the most fundamental tools for performing any finite volume calculation and includes a bunch of other files.

int main(int argc, char *argv[])
{
    #include "setRootCase.H" // Checks folder structure of the case
    #include "createTime.H" // Checks runtime according to the controlDict and initiates time variables
    #include "createMesh.H" // Defines mesh in the domain
    #include "createFields.H" // Creates fields (e.g. u, p, T) according to createFields.H
    #include "initContinuityErrs.H" // Declare and initialise the cumulative continuity error

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop()) // Start time loop
    {
        Info<< "Time = " << runTime.timeName() << nl << endl; // Print the current time

        #include "readPISOControls.H" // Read control parameters from fvSchemes
        #include "CourantNo.H" // Calculates and outputs the mean and maximum Courant Numbers

        // Set up the linear algebra for the momentum equation.
        // The flux of U, phi, is treated explicity
        // using the last known value of U.

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        // Solve using the last known value of p on the RHS.
        // This gives us a velocity field that is
        // not divergence free, but approximately satisfies momentum.
        // See Eqn. 7.31 of Ferziger & Peric (predicted velocity)

        solve(UEqn == -fvc::grad(p));

        // --- PISO loop

        for (int corr=0; corr<nCorr; corr++)
        {

            // From the last solution of velocity, extract the diag. term from the matrix and store the reciprocal
            // Note that the matrix coefficients are functions of U due to the non-linearity of convection.

            volScalarField rAU(1.0/UEqn.A());

            // take a Jacobi pass and update U.  See Hrv Jasak's thesis eqn. 3.137 and Henrik Rusche's thesis, eqn. 2.43
            // UEqn.H is the right-hand side of the UEqn minus the product of (the off-diagonal terms and U).
            // Note that since the pressure gradient is not included in the UEqn. above, this gives us U without
            // the pressure gradient.  Also note that UEqn.H() is a function of U.

            volVectorField HbyA("HbyA", U);

            // Calculate the fluxes by dotting the interpolated velocity (to cell faces) with face normals
            // The ddtPhiCorr term accounts for the divergence of the face velocity field by taking out the
            // difference between the interpolated velocity and the flux.

            HbyA = rAU*UEqn.H();
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

            // Adjusts the inlet and outlet fluxes to obey continuity, which is necessary for creating a well-posed
            // problem where a solution for pressure exists.

            adjustPhi(phiHbyA, U, p);

            // Iteratively correct for non-orthogonality.  The non-orthogonal part of the Laplacian is calculated from the most recent
            // solution for pressure, using a deferred-correction approach.

            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {

                // Set up the pressure equation
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                // In incompressible flow, only relative pressure matters.  Unless there is a pressure BC present,
                // one cell's pressure can be set arbitrarily to produce a unique pressure solution

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();

                // On the last non-orthogonality correction, correct the flux using the most up-to-date pressure
                // The .flux method includes contributions from all implicit terms of the pEqn (the Laplacian)

                if (nonOrth == nNonOrthCorr)
                {
                    phi = phiHbyA - pEqn.flux();
                }
            } // end of non-orthogonality looping

            #include "continuityErrs.H"

            // Add pressure gradient to interior velocity and BC's.
            // Note that this pressure is not just a small correction to a previous pressure,
            // but is the entire pressure field.
            // Contrast this to the use of p' in Ferziger & Peric, Eqn. 7.37.

            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    } // end of the time step loop

    Info<< "End\n" << endl;

    return 0;
    }

dimensionedScalar kappa
(
    transportProperties.lookup ("kappa")
);

dimensionedScalar hCoeff = kappa/( rho0 * Cv );

+ beta * g *( theta0 - theta )

fvScalarMatrix tempEqn
(
    fvm :: ddt (theta)
    + fvm :: div (phi, theta)
    - fvm :: laplacian (hCoeff, theta)
);
tempEqn.solve ();

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;          // 1= values given in meters, 0.001 = millimeters

vertices                    // it doesn't matter about order here
(
    (0 0 -0.05)
    (0.6 0 -0.05)
    (0 0.2 -0.05)
    (0.6 0.2 -0.05)
    (3 0.2 -0.05)
    (0 1 -0.05)
    (0.6 1 -0.05)
    (3 1 -0.05)
    (0 0 0.05)
    (0.6 0 0.05)
    (0 0.2 0.05)
    (0.6 0.2 0.05)
    (3 0.2 0.05)
    (0 1 0.05)
    (0.6 1 0.05)
    (3 1 0.05)
);

blocks                      // order must be right hand rule (x y z) grading
(
    hex (0 1 3 2 8 9 11 10) (25 10 1) simpleGrading (1 1 1)
    hex (2 3 6 5 10 11 14 13) (25 40 1) simpleGrading (1 1 1)
    hex (3 4 7 6 11 12 15 14) (100 40 1) simpleGrading (1 1 1)
);

edges
(
);

boundary                    // right hand rule again - normal vector must point outwards
(
    inlet
    {
        type patch;
        faces
        (
            (0 8 10 2)
            (2 10 13 5)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (4 7 15 12)
        );
    }
    bottom
    {
        type symmetryPlane;
        faces
        (
            (0 1 9 8)
        );
    }
    top
    {
        type symmetryPlane;
        faces
        (
            (5 13 14 6)
            (6 14 15 7)
        );
    }
    obstacle
    {
        type patch;
        faces
        (
            (1 3 11 9)
            (3 4 12 11)
        );
    }
);

mergePatchPairs
(
);

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

6
(
    inlet
    {
        type            patch;
        nFaces          50;
        startFace       10325;
    }
    outlet
    {
        type            patch;
        nFaces          40;
        startFace       10375;
    }
    bottom
    {
        type            symmetryPlane;
        inGroups        1(symmetryPlane);
        nFaces          25;
        startFace       10415;
    }
    top
    {
        type            symmetryPlane;
        inGroups        1(symmetryPlane);
        nFaces          125;
        startFace       10440;
    }
    obstacle
    {
        type            patch;
        nFaces          110;
        startFace       10565;
    }
    defaultFaces                        // unspecified front and back planes
    {
        type            empty;
        inGroups        1(empty);
        nFaces          10500;
        startFace       10675;
    }
)

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 1;                  // 1 Pa

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 1;          // 1 Pa
    }

    outlet
    {
        type            waveTransmissive;
        field           p;
        phi             phi;
        rho             rho;
        psi             thermo:psi;
        gamma           1.4;
        fieldInf        1;
        lInf            3;
        value           uniform 1;
    }

    bottom
    {
        type            symmetryPlane;
    }

    top
    {
        type            symmetryPlane;
    }

    obstacle
    {
        type            zeroGradient;       // Neumann boundary condition (derivative)
    }

    defaultFaces
    {
        type            empty;
    }
}

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 1;              // 1K(!)

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 1;      // 1K(!)
    }

    outlet
    {
        type            zeroGradient;   // Neumann
    }

    bottom
    {
        type            symmetryPlane;
    }

    top
    {
        type            symmetryPlane;
    }

    obstacle
    {
        type            zeroGradient;
    }

    defaultFaces
    {
        type            empty;
    }
}

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (3 0 0);            // 3 m/s

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (3 0 0);    // 3 m/s
    }

    outlet
    {
        type            zeroGradient;
    }

    bottom
    {
        type            symmetryPlane;
    }

    top
    {
        type            symmetryPlane;
    }

    obstacle
    {
        type            fixedValue;         // No slip
        value           uniform (0 0 0);
    }

    defaultFaces
    {
        type            empty;
    }
}

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType  laminar;


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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

// Note: these are the properties for a "normalised" inviscid gas
//       for which the speed of sound is 1 m/s at a temperature of 1K
//       and gamma = 7/5
mixture
{
    specie
    {
        nMoles          1;
        molWeight       11640.3;
    }
    thermodynamics
    {
        Cp              2.5;        // Specific heat
        Hf              0;
    }
    transport
    {
        mu              0;
        Pr              1;          // Prandtl number
    }
}


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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     sonicFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.002;

writeControl    runTime;

writeInterval   0.5;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

blocks                      // order must be right hand rule (x y z) grading
(
    hex (0 1 3 2 8 9 11 10) (25 10 1) simpleGrading (0.5 1 1)
    hex (2 3 6 5 10 11 14 13) (25 40 1) simpleGrading (0.5 1 1)
    hex (3 4 7 6 11 12 15 14) (100 40 1) simpleGrading (1 1 1)
);

blocks                      // order must be right hand rule (x y z) grading
(
    hex (0 1 3 2 8 9 11 10) (50 20 1) edgeGrading (1 1 1 1 2 2 2 2 1 1 1 1)
    hex (2 3 6 5 10 11 14 13) (25 40 1) simpleGrading (1 1 1)
    hex (3 4 7 6 11 12 15 14) (100 40 1) simpleGrading (1 1 1)
);

defaultFieldValues ( volVectorFieldValue U ( 0 0 0 ) volScalarFieldValue T 348.432 volScalarFieldValue p 100000 );

regions         (
// boxToCell { box ( 0 -1 -1 ) ( 5 1 1 ) ; fieldValues ( volScalarFieldValue T 278.746 volScalarFieldValue p 10000 ) ; }
);

defaultFieldValues ( volVectorFieldValue U ( 0 0 0 ) volScalarFieldValue T 348.432 volScalarFieldValue p 100000 );

regions         (
boxToCell { box ( 0 -1 -1 ) ( 5 1 1 ) ; fieldValues ( volScalarFieldValue T 278.746 volScalarFieldValue p 10000 ) ; }
);

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

interpolationScheme cellPoint;

setFormat       raw;

sets
(
    data
    {
        type    uniform;
        axis    x;
        start   (-4.995 0 0);
        end     (4.995 0 0);
        nPoints 10000;
    }
);

fields          (T magU p);

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