/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Author
    Dr. Luca Mangani, HSLU email: luca.mangani@hslu.ch.  All rights reserved.

\*---------------------------------------------------------------------------*/


Info<< "Reading thermophysical properties\n" << endl;

autoPtr<psiThermo> thermo
(
    psiThermo::New(mesh)
);

Info<< "Calculating field rho\n" << endl;

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),        
    thermo->rho()
);

volScalarField& p = thermo->p();
volScalarField& h = thermo->he();

thermo->correct();

volScalarField gammaGas
(
    IOobject
    (
        "gammaGas",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo->Cp() /  thermo->Cv()
); 

volScalarField RGas
(
    IOobject
    (
        "RGas",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo->Cp() -  thermo->Cv()
); 


Info << "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


surfaceScalarField mdotf
(
    IOobject
    (
        "mdotf",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(rho*U) & mesh.Sf()
);

Info<< "Creating turbulence model\n" << endl;

autoPtr<compressible::RASModel> turbulence
(
    compressible::RASModel::New
    (
        rho,
        U,
        mdotf,
        thermo()
    )
);

const volScalarField::GeometricBoundaryField& pbf = p.boundaryField();

wordList pbt = pbf.types();

forAll(pbf, patchi)
{
    if (isType<totalPressureCompFvPatchScalarField>(p.boundaryField()[patchi]))
    {
        pbt[patchi] = totalPressureCorrectorCompFvPatchScalarField::typeName;
    }
}

volScalarField pp
(
    IOobject
    (

        "pp",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("zero", p.dimensions(), 0.0),
    pbt
);

// Set pcorr boundary values
forAll(pp.boundaryField(), patchi)
{
    if (isType<totalPressureCompFvPatchScalarField>(p.boundaryField()[patchi]))
    {           

         totalPressureCorrectorCompFvPatchScalarField& ppbound =
                         refCast<totalPressureCorrectorCompFvPatchScalarField>(pp.boundaryField()[patchi]);

               totalPressureCompFvPatchScalarField pbound =
                         refCast<totalPressureCompFvPatchScalarField>(p.boundaryField()[patchi]);

               ppbound.p0() = pbound.p0();
    }
    else if (isType<fixedValueFvPatchScalarField>(pp.boundaryField()[patchi]))
    {
       fixedValueFvPatchScalarField& ppbound =
                       refCast<fixedValueFvPatchScalarField>(pp.boundaryField()[patchi]);

       ppbound == scalarField(ppbound.size(),0.0);
    }
}

scalar alphaP = mesh.equationRelaxationFactor("pp");

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);

dimensionedScalar rhoMax(simple.dict().lookup("rhoMax"));
dimensionedScalar rhoMin(simple.dict().lookup("rhoMin"));

dimensionedScalar pMax(simple.dict().lookup("pMax"));
dimensionedScalar pMin(simple.dict().lookup("pMin"));
