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

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

internalField   uniform (0 0 0);

boundaryField
{
    roof
    {
        type          noSlip;
    }

    outlet
    {          
        type 		pressureInletOutletVelocity;
        value		$internalField;
    }

  "(bottom.*)"
    {
        name     surfaceEvaporationVelocity;
        type     codedFixedValue;
        value    uniform (0 0 0);
        code
        #{
            // Parameters

            // Read thermophysical properties dictonary
            const auto thermoPhys = this->db().objectRegistry::lookupObject<IOdictionary>("thermophysicalProperties");

            // Extract tank geometry and heat transfer properties
            const scalar tankHeight = readScalar(thermoPhys.subDict("tank").lookup("tankHeight"));
            const scalar valve_radius = readScalar(thermoPhys.subDict("tank").lookup("valve_radius"));
            const scalar U_L = readScalar(thermoPhys.subDict("tank").lookup("U")); // Liquid overall heat transfer coefficient W/m^2K
            const scalar R_T = readScalar(thermoPhys.subDict("tank").lookup("R_T")); // Tank radius [m]
            const scalar radius_ext =  readScalar(thermoPhys.subDict("tank").lookup("radius_ext")); // Tank radius [m]
            const scalar T_air = readScalar(thermoPhys.subDict("tank").lookup("T_air")); //288.1500; // Air temperature [K]
            const scalar T_sat = readScalar(thermoPhys.subDict("mixture").subDict("thermodynamics").lookup("T_sat"));
          // Thermophysical

            const scalar Cp = readScalar(thermoPhys.subDict("mixture").subDict("thermodynamics").lookup("Cp"));
            const scalar Pr = readScalar(thermoPhys.subDict("mixture").subDict("transport").lookup("Pr"));
            const scalar mu = readScalar(thermoPhys.subDict("mixture").subDict("transport").lookup("mu"));

            const scalar k_V = mu*Cp/Pr; // Thermal conductivity
            const scalar rho_V = readScalar(thermoPhys.subDict("mixture").subDict("thermodynamics").lookup("rho_Tsat"));
            const scalar dH_LV = readScalar(thermoPhys.subDict("mixture").subDict("thermodynamics").lookup("dH_LV"));

            // Access to the patch and the mesh
            const fvMesh& mesh = patch().boundaryMesh().mesh();
            const scalarField& V = mesh.V(); // Get the volume of the mesh
            const scalar V_vap = gSum(V);

            // Get the area of the bottom patch
            label interface_1 = mesh.boundaryMesh().findPatchID("bottom_1");
            label interface_2 = mesh.boundaryMesh().findPatchID("bottom_2");
            const scalarField& interfaceAreas1 = mesh.boundary()[interface_1].magSf();
            const scalarField& interfaceAreas2 = mesh.boundary()[interface_2].magSf();

            const scalar A1 = gSum(interfaceAreas1);
            const scalar A2 = gSum(interfaceAreas2);

            // Convert wedge area to cylindrical equivalent
            const scalar cyl_over_wedge = sqrt((2.5/180) / (pow(cos(2.5/180),3)*sin(2.5/180) ));
            
            // Calculate the liquid height
            const scalar liquidHeight = tankHeight - V_vap/(A1+A2) * cyl_over_wedge;

            // Calculate heat transfer areas
            const scalar A_wall = 2*constant::mathematical::pi*radius_ext*liquidHeight;
            const scalar A_int = constant::mathematical::pi*pow(R_T,2); // Interfacial Area : based in internal Area
            const scalar A_bottom = constant::mathematical::pi*pow(radius_ext,2);

            // Access to the liquid temperature
            const volScalarField& T = this->db().objectRegistry::lookupObject<volScalarField>  ("T");
            const fvPatchField<double>& T_int_1 = T.boundaryField()[interface_1];
            const fvPatchField<double>& T_int_2 = T.boundaryField()[interface_2];

            const scalarField grad_T_int_1 = T_int_1.snGrad();
            const scalarField grad_T_int_2 = T_int_2.snGrad();

            // Recall that the wedge mesh is divided in 360/2.5 parts
            // We will calculate the average gradient at the interface
            // And then multiply it by the area of the cylinder

            const scalar A1_exact = constant::mathematical::pi * pow((valve_radius*R_T),2);
            const scalar A2_exact = A_int - A1_exact;

            const double Q_VL_1 = - k_V * A1_exact * sumProd(grad_T_int_1, interfaceAreas1)/A1;
            const double Q_VL_2 = - k_V * A2_exact * sumProd(grad_T_int_2, interfaceAreas2)/A2;

            const scalar T_interface = T_sat;

            // Calculate wall heat ingress
            // ASSUMPTION: THE BOTTOM HEAT INGRESS WILL BE NEGLECTED
            const scalar Q_L = U_L*(A_wall)*(T_air-T_interface);
            // Calculate vapour to liquid heat ingress
            const scalar Q_VL = Q_VL_1 + Q_VL_2;
            // Calculate total heat ingress (W)
            const scalar Q_L_tot = Q_VL + Q_L;

            // Calculate inlet velocity
            const scalar v_evap = Q_L_tot/(dH_LV*rho_V*A_int);
            // Info << "v_inlet = " << v_evap << endl;
            const vector axis (0,1,0);
            vectorField v(this->patch().Cf());

            v.replace(vector::X, 0);
            v.replace(vector::Y, v_evap);
            v.replace(vector::Z, 0);

             operator==(v);

             bool debug = false;
             if(debug)
             {
                 Info << "Printing velocity" << endl;
                 Info << v << endl;
                 Info << "Tank height = " << tankHeight << " m " << endl;
                 Info << "Liquid height = " << liquidHeight << " m " << endl;
                 Info << "A_wall = " << A_wall << " m^2 " << endl;
                 Info << "Interface area = " << A1 + A2 << " m^2 " <<  endl;
                 Info << "Wall wet area = " << A_wall << " m^2 " << endl;;
                 Info << "Bottom area = " << A_bottom<<" m^2 " << endl;


                 Info << "T_interface = " << T_interface << " K " << endl;
                 Info << "Vapour volume = " << V_vap << " m^3 " << endl;
                 Info << "valve_radius = " << valve_radius << " m " << endl;
                 Info << "U_L = " << U_L << " W/m^2K" << endl;
                 Info << "R_T = " << R_T << " m" << endl;
                 Info << "External radius = " << radius_ext << " m " << endl;
                 Info << "T_air = " << T_air << " K" << endl;
                 Info << "Cp = " << Cp << " J/kgK" <<endl;
                 Info << "Pr = " << Pr << endl;
                 Info << "mu = " << mu << " kg/(m^2*s)" <<endl;
                 Info << "k_V = " << k_V << " W/mK " << endl;
                 Info << "rho_V_Tsat = " << rho_V << " kg/m^3" <<endl;
                 Info << "dH_LV = " << dH_LV << " J/kg" <<endl;

                 Info << "Q_VL = " << Q_VL << " W " << endl;
                 Info << "Q_L = " << Q_L << " W " << endl;   
                 Info << "Liquid heat ingress = " << Q_L_tot << " W " << endl;
             }
            
        #};
        value    $internalField;
    }

    tank_wall
    {
        type            noSlip;
    }

    valve_wall
    {
        type            noSlip;
    }

    "(wedge.*)"
    {
        type            wedge;
    }
}

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