/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1906                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      binary;
    class       volScalarField;
    arch        "LSB;label=32;scalar=64";
    location    "0";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField  uniform 111.508;

boundaryField
{
    bottom_1
    {
        type            fixedValue;
        value           uniform 111.508;
    }
    outlet
    {
        type            zeroGradient;
    }
    wedge_0
    {
        type            wedge;
    }
    wedge_1
    {
        type            wedge;
    }
    wedge_2
    {
        type            wedge;
    }
    bottom_2
    {
        type            fixedValue;
        value           uniform 111.508;
    }
    tank_wall
    {
        type            codedMixed;
        refValue        uniform 288.15;
        refGradient     uniform 0;
        valueFraction   uniform 0;
        value           uniform 111.508;
        name            convection_Conduction;
        codeInclude     #{
            #include "fvCFD.H"
            #include "dynamicFvMesh.H"
            #include "rhoThermo.H"
            #include "turbulentFluidThermoModel.H"
            #include "radiationModel.H"
            #include "fvOptions.H"
            #include "pimpleControl.H"
            #include "pressureControl.H"
        #};
    
        codeOptions     #{
        -I$(FOAM_SOLVERS)/heatTransfer/buoyantPimpleFoam \
        -I$(FOAM_SOLVERS)/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam \
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/sampling/lnInclude \
        -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
        -I$(LIB_SRC)/dynamicMesh/lnInclude \
        -I$(LIB_SRC)/overset/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude \
        -I$(LIB_SRC)/transportModels/compressible/lnInclude \
        -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
        -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
        -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
        -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
               #};

        code            #{
            const fvPatch& boundaryPatch = patch();     // Access to the patch
            const fvMesh& mesh = patch().boundaryMesh().mesh(); // Access to the mesh
            const scalarField& deltaCoeffs = boundaryPatch.deltaCoeffs();

            const compressibleTurbulenceModel& tur = this->db().objectRegistry::lookupObject<compressibleTurbulenceModel>("turbulenceProperties");
            // Read turbulent thermal difussivity
            const volScalarField& alphat = this->db().objectRegistry::lookupObject<volScalarField>("alphat");
            const volScalarField& rho = this->db().objectRegistry::lookupObject<volScalarField>("rho");
            const auto thermoPhys = this->db().objectRegistry::lookupObject<IOdictionary>("thermophysicalProperties");

            // Read laminar Prandtl number from thermophysicalProperties dictionary
            const scalar Pr = readScalar(thermoPhys.subDict("mixture").subDict("transport").lookup("Pr"));
            const scalar Cp = readScalar(thermoPhys.subDict("mixture").subDict("thermodynamics").lookup("Cp"));
            // Get laminar viscosity
            const volScalarField nu  = tur.nu();///rho; // This is the dynamic viscosity
            const volScalarField alphaEff = alphat/rho + nu/Pr; // Alphat is in kg/(m*s) where it should be m^2/s.

            // Find Patch ID of the tank wall
            label tankWall = mesh.boundaryMesh().findPatchID("tank_wall");
            // Wall thermal diffusivity
            const fvPatchField<double>& alphatWall = alphaEff.boundaryField()[tankWall];
            const fvPatchField<double>& rhoWall = rho.boundaryField()[tankWall];
            const scalarField k_eff = alphatWall*rhoWall*Cp;

            const scalar U = readScalar(thermoPhys.subDict("tank").lookup("U"));
            const scalar T_air = readScalar(thermoPhys.subDict("tank").lookup("T_air"));

            bool debug = false;

            if(debug)
            {
            Info << "alpha =" << nu[0]/Pr << endl;
            Info << "alphaEff =" << alphatWall[0] << endl;
            Info << "rhoWall =" << rhoWall << endl;
            Info << "alphat =" << alphat.boundaryField()[tankWall] << endl;
            //Info << "alphat 2 =" << tur.nut().boundaryField()[tankWall]/Prt << endl;
            Info << "kEff = " << k_eff << endl;
            Info << "k_lam =" << nu[0]/Pr*rho[0]*Cp << endl;
            Info << Pr << endl;
            }
            Info << "U = " << U << endl;
            Info << "T_air = " << T_air << endl;

            // Update boundary conditions
            this->refValue() = T_air;
            this->refGrad() = 0;       // 0 for convection-conduction bc
            this->valueFraction() = 1.0/(1.0 + k_eff/(U*(1/deltaCoeffs)) ); // Robin BC
        #};

    }
    roof
    {
        type            zeroGradient;
    }
    wedge_3
    {
        type            wedge;
    }
}


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