/*
 * Cooldt cool-down time calculator for multi-layers pipe
 *
 * Copyright 2010, 2013, 2014, Benjamin DEGLO DE BESSES. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification, are
 * permitted provided that the following conditions are met:
 *
 *    1. Redistributions of source code must retain the above copyright notice, this list of
 *       conditions and the following disclaimer.
 *
 *    2. Redistributions in binary form must reproduce the above copyright notice, this list
 *       of conditions and the following disclaimer in the documentation and/or other materials
 *       provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY Benjamin DEGLO DE BESSES "AS IS" AND ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
 * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Benjamin DEGLO DE BESSES OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include "cooldt.h"

/*
 * Discretize the fluid and pipe layers
 * Calculate interface conductivity km and kp
 * Calculate cells surface
 */
void MeshPipe(PipeProblem *PP, PipeDicretize *PZ)
{
    int l ;
    int i ;

    /*
     * For each layer
     */
    for( l = 0 ; l <= PP->N ; l++ )
    {
        /*
         * Special case : the Pipe/fluid (layer 0)
         */
        if (l == 0)
        {
            double HalfInterval = PP->r[l] / ( 2.0 * (double)PP->NbCell[l] - 1.0 ) ;

            Size0PipeDiscretize(PZ) ;

            PZ->L[PZ->N] = l ;

            PZ->rm[PZ->N] = 0.0 ;
            PZ->rc[PZ->N] = 0.0 ;
            PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;

            PZ->Cp[PZ->N] = PP->Cp[l] ;
            PZ->d[PZ->N] = PP->d[l] ;

            PZ->kc[PZ->N] = PP->k[l] ;

            /*
             * Dicretize the others intervals
             */
            for( i = 2 ; i <= PP->NbCell[l] ; i++ )
            {
                IncrementSizePipeDiscretize(PZ) ;

                PZ->L[PZ->N] = l ;

                PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
                PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
                PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;

                PZ->Cp[PZ->N] = PP->Cp[l] ;
                PZ->d[PZ->N] = PP->d[l] ;

                PZ->kc[PZ->N] = PP->k[l] ;
            }
        }

        /*
         * Special case : the last layer
         */
        else if ( l == PP->N )
        {
            double HalfInterval = ( PP->r[l] - PP->r[l-1] ) / ( 2.0 * (double)PP->NbCell[l] - 1.0 ) ;

            for( i = 2 ; i <= PP->NbCell[l] ; i++ )
            {
                IncrementSizePipeDiscretize(PZ) ;

                PZ->L[PZ->N] = l ;

                PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
                PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
                PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;

                PZ->Cp[PZ->N] = PP->Cp[l] ;
                PZ->d[PZ->N] = PP->d[l] ;

                PZ->kc[PZ->N] = PP->k[l] ;
            }

            IncrementSizePipeDiscretize(PZ) ;

            PZ->L[PZ->N] = l ;

            PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
            PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
            PZ->rp[PZ->N] = PZ->rc[PZ->N] ;

            PZ->Cp[PZ->N] = PP->Cp[l] ;
            PZ->d[PZ->N] = PP->d[l] ;

            PZ->kc[PZ->N] = PP->k[l] ;
        }

        /*
         * Normal case: the other layers
         */
        else
        {
            double HalfInterval = ( PP->r[l] - PP->r[l-1] ) / ( 2.0 * (double)PP->NbCell[l] ) ;

            for( i = 1 ; i <= PP->NbCell[l] ; i++ )
            {
                IncrementSizePipeDiscretize(PZ) ;

                PZ->L[PZ->N] = l ;

                PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
                PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
                PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;

                PZ->Cp[PZ->N] = PP->Cp[l] ;
                PZ->d[PZ->N] = PP->d[l] ;

                PZ->kc[PZ->N] = PP->k[l] ;

            }
        }
    }

    /*
     * Calculate k[i-1/2] and k[i+1/2] conductivity (resp. km and kp)
     */
    i = 0 ;
    PZ->km[i] = PZ->kc[i] ;
    PZ->kp[i] = PZ->kc[i] ;

    for( i = 1 ; i <= (PZ->N-1) ; i++ )
    {
        PZ->km[i] = PZ->kp[i-1] ;
        PZ->kp[i] = PZ->kc[i] * PZ->kc[i+1] * log(PZ->rc[i+1]/PZ->rc[i])
                    / ( PZ->kc[i+1] * log(PZ->rp[i]/PZ->rc[i]) + PZ->kc[i] * log(PZ->rc[i+1]/PZ->rp[i]) ) ;
    }

    i = PZ->N ;
    PZ->km[i] = PZ->kp[i-1] ;
    PZ->kp[i] = PZ->kc[i] ;

    /*
     * Calculate V[i] = Cells Volume / unit length = Surface
     */
    for( i = 0 ; i <= PZ->N ; i++ )
        PZ->S[i] = C_Pi * ( PZ->rp[i]*PZ->rp[i] - PZ->rm[i]*PZ->rm[i] ) ;
}

/*
 * Alpha, Beta & Gamma coefficients calculation
 */
void CalcModelCoefs(PipeProblem *PP, PipeDicretize *PZ)
{
    /*
     * At the pipe center, i=0
     */
    int i = 0 ;
    PZ->alpha[i] = 0.0 ;
    PZ->beta[i]  = PZ->S[i] * PZ->d[i] * PZ->Cp[i] ;
    PZ->gamma[i] = 2.0 * C_Pi * PZ->rp[i] * PZ->kp[i] / ( PZ->rc[i+1] - PZ->rc[i] ) ;

    /*
     * In the pipe layers, 2<=i<=N-1
     */
    for( i = 1 ; i < PZ->N ; i++ )
    {
        PZ->alpha[i] = PZ->gamma[i-1] ;
        PZ->beta[i]  = PZ->S[i] * PZ->d[i] * PZ->Cp[i] ;
        PZ->gamma[i] = 2.0 * C_Pi * PZ->rp[i] * PZ->kp[i] / ( PZ->rc[i+1] - PZ->rc[i] ) ;
    }

    /*
     * At the ambiant fluid / outer wall interface, i=N
     */
    i = PZ->N ;
    PZ->alpha[i] = PZ->gamma[i-1] ;
    PZ->beta[i]  = PZ->S[i] * PZ->d[i] * PZ->Cp[i] ;
    PZ->gamma[i] = 2.0 * C_Pi * PZ->rc[i] * PP->Hout ;

    /*
     * AlphaSharp Coef for fluid/solid temperature calculation
     */
    i = PP->NbCell[0] ;
    PZ->alphaSharp = 2.0 * C_Pi * PZ->rc[i] * PZ->kc[i] / ( PZ->rc[i] - PZ->rm[i] ) ;
}

/*
 * Build the Pipe discrete model
 */
void BuildModel(PipeProblem *PP, PipeDicretize *PZ)
{
    int l ;

    /*
     * Apply discretization limitations
     */

    /* For each layer */
    for( l = 0 ; l <= PP->N ; l++ )
    {
        if ( PP->NbCell[l] < 1 )  PP->NbCell[l] = 1 ;
        if ( PP->NbCell[l] > MaxCellNumberPerLayer )  PP->NbCell[l] = MaxCellNumberPerLayer ;
    }

    /* For time step */
    if (PP->dt < MinTimeStep)  PP->dt = MinTimeStep ;

    /*
     * Mesh the pipe and calculate alpha, beta & gamma coefs
     */
    MeshPipe(PP, PZ) ;
    CalcModelCoefs(PP, PZ) ;
}

/* End of mesh.c */