Cooldt Code

/*
 * 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"

/*
 * Calculate transient problem solution with automatic accurancy control
 */
void aaControl(FILE *fp, PipeProblem *PP)
{
    int* aacNbCell = (int*)malloc((PP->N + 1)*sizeof(int)) ;
    double aacdt ;
    double* aacCDT = (double*)malloc((PP->N + 2)*sizeof(double)) ;
    double* aacError = (double*)malloc((PP->N + 2)*sizeof(double)) ;
    double CDT ;
    double MaxCdtError = 1.0 ;
    int i ;
    int k = 0 ;

    /*
     * Write the header of the aac report
     */
    fprintf(fp, "#\tIteration\tError\tCDT\tTime step") ;

    for( i=0 ; i <=PP->N ; i++ )
        fprintf(fp, "\tLayer %i", i);

    fprintf(fp, "\n#\t(-)\t(%%)\t(h)\t(s)") ;

    for( i=0 ; i <=PP->N ; i++ )
        fprintf(fp, "\tCell #");

    fprintf(fp, "\n") ;

    /*
     * Start with two cells for fluid, one for others layers
     */
    aacNbCell[0] = PP->NbCell[0] = aacInitialFluidCellNumber ;
    for( i=1 ; i <=PP->N ; i++ )
        aacNbCell[i] = PP->NbCell[i] = aacInitialSolidCellNumber ;

    /*
     * Start with a large time step
     */
    aacdt = PP->dt = aacInitialTimeStep ;

    /*
     * Build the model, Calculate Steady State then Transient
     */
    CDT = SolutionSimpleOutput(PP) ;

    /*
     * Improves Parameter while the CDT has not converged
     */
    while( MaxCdtError > aacCDTconverge )
    {
        /*
         * Print the aac report for the reference case
         */
        fprintf(fp, "#\t%i\t%f\t%f\t%.0e", k+1, 100.0 * MaxCdtError, CDT/3600.0, aacdt) ;

        for( i=0 ; i <=PP->N ; i++ )
            fprintf(fp, "\t%i", aacNbCell[i]) ;

        fprintf(fp, "\n") ;

        /*
         * Maximum searching initialization
         */
        int IndexMaxCdtError = 0 ;
        MaxCdtError = -1.0 ;

        /*
         * Calculate the CDTs for each parameter improvement
         */
        for( i=0 ; i <=(PP->N+1) ; i++ )
        {
            int j ;

            /*
             * Copy parameters to PP structure
             */
            for( j=0 ; j <=PP->N ; j++ )
            {
                if( i == j) /* Increase the number of cell */
                    PP->NbCell[j] =  aacNbCellGroth * aacNbCell[j] ;
                else
                    PP->NbCell[j] = aacNbCell[j] ;
            }

            if( i == PP->N+1) /* Reduce the time step */
                PP->dt = aacTimeStepDecay * aacdt ;
            else
                PP->dt = aacdt ;

            /*
             * Build the model, Calculate Steady State then Transient
             */
            aacCDT[i] = SolutionSimpleOutput(PP) ;

            /*
             * Calculate the relative error / reference case
             */
            aacError[i] = fabs( (aacCDT[i] - CDT) / aacCDT[i] ) ;

            /*
             * Find for the Best Parameter / Solution improvement
             */
            if ( (aacError[i] > MaxCdtError) )
            {
                MaxCdtError = aacError[i] ;
                IndexMaxCdtError = i ;
            }
        }

        /*
         * Keep only the Best Parameter / Solution Improvement
         */
        CDT = aacCDT[IndexMaxCdtError] ;

        /*
         * Save the parameter set corresponding to the best solution improvement
         */
        if(IndexMaxCdtError == (PP->N+1) )
            aacdt = aacTimeStepDecay * aacdt ;
        else
            aacNbCell[IndexMaxCdtError] = aacNbCellGroth * aacNbCell[IndexMaxCdtError] ;

        /*
         * Increment the iteration counter
         */
        k++ ;
    }

    /*
     * Save the final parameter set to the PP structure
     */
    for( i=0 ; i <=PP->N ; i++ )
        PP->NbCell[i] = aacNbCell[i] ;

    PP->dt = aacdt ;

    /*
     * free the local pointers
     */
    free(aacError) ;
    free(aacNbCell) ;
    free(aacCDT) ;

    /*
     * Print the aac report after convergence
     */
    fprintf(fp, "#\t%i\t%f\t%f\t%.0e", k+1, 100.0 * MaxCdtError, CDT/3600.0, PP->dt) ;

    for( i=0 ; i <=PP->N ; i++ )
        fprintf(fp, "\t%i", PP->NbCell[i])  ;

    fprintf(fp, "\n#\n") ;

}

/* End of automatic.c */