VPAC - Computational Software Development

Main | SPModel | StGermain FrameWork |

Context.c

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
**
** Copyright (C), 2004, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
**
** Authors:
**  Ogar R. Widjaja, Computational Scientist, VPAC.
**  Raquibul Hassan, Software Engineer, VPAC. (raq@vpac.org)
**  Keith Hsuan, Computational Scientist, VPAC (keith@vpac.org)
**  William F. Appelbe, Director, VPAC. (bill@vpac.org)
**  Stevan M. Quenette, Senior Software Engineer, VPAC. (steve@vpac.org)
**  Patrick D. Sunter, Software Engineer, VPAC. (patrick@vpac.org)
**
** This file may be distributed under the terms of the VPAC Public License
** as defined by VPAC of Australia and appearing in the file
** LICENSE.VPL included in the packaging of this file.
**
** This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
** WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
**
** $Id: Context.c 285 2006-03-20 07:26:25Z RaquibulHassan $
**
**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

#include <mpi.h>
#include <StGermain/StGermain.h>

#include "types.h"
#include "Context.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <assert.h>
#include "SurfaceMesh.h"
#include "SurfaceRegularMesh.h"
#include "SurfaceMeshDecomp.h"
#include "SurfaceMeshIrregularDecomp.h"
#include "SurfaceMeshRegularDecomp.h"
#include "Catchment.h"
#include "_Interpolator.h"
#include "LinearInterpolator.h"
#include "SplineInterpolator.h"
#include "SurfaceMeshSmoother.h"
#include "ParameterTimeSeries.h"
#include "SurfaceMeshCyclicBC.h"
#include "SurfaceMeshRectangularCyclicBC.h"
#include "CommHandler.h"

/* Textual name of this class */
const Type SPModel_Context_Type = "SPModel_Context";

ExtensionInfo_Index SPModel_Simulation_ContextExtHandle;

SPModel_Context* SPModel_Context_DefaultNew( Name name )
{
    return _SPModel_Context_New( sizeof(SPModel_Context), SPModel_Context_Type,
            _SPModel_Context_Delete, _SPModel_Context_Print, NULL, (void*)SPModel_Context_DefaultNew,
            _SPModel_Context_ComponentConstruct, _SPModel_Context_ComponentBuild, _SPModel_Context_ComponentInitialise,
            _SPModel_Context_ComponentExecute, _SPModel_Context_ComponentDestroy, name, False,
        _SPModel_Context_SetDt, 0, 0, 0, NULL );
}

SPModel_Context* SPModel_Context_New(
        Name                    name,
        double                  startTime,
        double                  stopTime,
        MPI_Comm                communicator,
        Dictionary*             dictionary )
{
    return _SPModel_Context_New( sizeof(SPModel_Context), SPModel_Context_Type,
            _SPModel_Context_Delete, _SPModel_Context_Print, NULL, (void*)SPModel_Context_DefaultNew,
            _SPModel_Context_ComponentConstruct, _SPModel_Context_ComponentBuild, _SPModel_Context_ComponentInitialise,
            _SPModel_Context_ComponentExecute, _SPModel_Context_ComponentDestroy, name, True,
        _SPModel_Context_SetDt, startTime, stopTime, communicator, dictionary );
}


SPModel_Context* _SPModel_Context_New(
        SizeT                   _sizeOfSelf,
        Type                    type,
        Stg_Class_DeleteFunction*               _delete,
        Stg_Class_PrintFunction*                _print,
        Stg_Class_CopyFunction*             _copy, 
        Stg_Component_DefaultConstructorFunction*   _defaultConstructor,
        Stg_Component_ConstructFunction*            _construct,
        Stg_Component_BuildFunction*        _build,
        Stg_Component_InitialiseFunction*       _initialise,
        Stg_Component_ExecuteFunction*      _execute,
        Stg_Component_DestroyFunction*      _destroy,
        Name                            name,
        Bool                            initFlag,
        AbstractContext_SetDt*          _setDt,
        double                  startTime,
        double                  stopTime,
        MPI_Comm                communicator,
        Dictionary*             dictionary )
{
    SPModel_Context*    self;

    /* Allocate memory */
    self = (SPModel_Context*)_AbstractContext_New( _sizeOfSelf, type, _delete, _print, NULL, 
            _defaultConstructor, _construct, _build, _initialise, _execute, _destroy, name, 
            initFlag, _setDt, startTime, stopTime, communicator, dictionary );

    /* General info */

    /* Virtual info */

    /* SPModel_Context info */
    if( initFlag ){
        _SPModel_Context_Init( self );
    }

    return self;
}


void _SPModel_Context_Init( SPModel_Context* self ) {

    Dictionary *fileOutputDict = NULL;
    Dictionary *cyclicBCDict = NULL;
    char *interpolationOption = NULL;
    
    if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );

    /* construct the mesh */
    self->globalMesh = SurfaceMesh_New( "globalMesh", self->nproc, self->rank, self->extensionMgr_Register, self->dictionary );
    self->localMesh = SurfaceMesh_New( "localMesh", self->nproc, self->rank, self->extensionMgr_Register, self->dictionary );
    self->catchmentList = CatchmentList_New( "catchments", self->globalMesh, self->nproc, self->dictionary );
    self->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshIrregularDecomp_New( "meshDecomp", self->catchmentList, self->globalMesh, self->dictionary );
    /*self->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshRegularDecomp_New( "meshDecomp", self->globalMesh, self->dictionary );*/
    self->meshSmoother = SurfaceMeshSmoother_New( "surfaceMeshSmoother", self->dictionary );
    self->regularMesh = SurfaceRegularMesh_New( "regularMesh", self->dictionary, self->localMesh );

    interpolationOption = Dictionary_GetString_WithDefault( self->dictionary, "interpolation", "Linear" );

    if( !strcmp(interpolationOption, "Linear") ){
        self->interpolator = (_Interpolator*)LinearInterpolator_New( "linearInterpolator", self->dictionary, self->localMesh, self->regularMesh );
    }
    else if ( !strcmp(interpolationOption, "Cubic") ){
        self->interpolator = (_Interpolator*)SplineInterpolator_New( "splineInterpolator", self->dictionary, self->localMesh, self->regularMesh );      
    }

    if ( (cyclicBCDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "CyclicBoundaryConditions" ) ))){
        self->cyclicBC = (_SurfaceMeshCyclicBC*)SurfaceMeshRectangularCyclicBC_New
                            ( "cyclicBC", self->globalMesh, self->dictionary );
    }
    
    /* defaults of physical parameters */
    self->gravity = Dictionary_GetDouble_WithDefault( self->dictionary, "gravity", 9.81f );
    self->flexureInterval = Dictionary_GetDouble_WithDefault( self->dictionary, "flexureInterval", 100 );
    
    self->seaLevel = ParameterTimeSeries_New( "seaLevel", 0.0f, self->dictionary );
    
    self->upliftRate = Dictionary_GetFloat_WithDefault( self->dictionary, "upliftRate", 1.0e-5 );
    
    self->oroLength = Dictionary_GetFloat_WithDefault( self->dictionary, "oroLength", 0 );
    self->oroRate = Dictionary_GetFloat_WithDefault( self->dictionary, "oroRate", 5.0e-1 );
    self->fluvialConstant = Dictionary_GetFloat_WithDefault( self->dictionary, "fluvialConstant",3.0e-2 );
    
    self->bedRockDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "bedRockDiffusivity", 1.0e-3 );
    self->sedimentDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "sedimentDiffusivity", 1.0e3 );
    self->submarineBedRockDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "submarineBedRockDiffusivity", 1.0e-3 );
    self->submarineSedimentDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "submarineSedimentDiffusivity", 1.0e3 );
    
    self->bedRockDensity = Dictionary_GetFloat_WithDefault( self->dictionary, "bedRockDensity", 2.8e3 );
    self->sedimentDensity = Dictionary_GetFloat_WithDefault( self->dictionary, "sedimentDensity", 2.8e3 );
    self->asthenosphereDensity = Dictionary_GetFloat_WithDefault( self->dictionary, "asthenosphereDensity", 3.3e3 );
    
    self->flexuralRigidity = Dictionary_GetFloat_WithDefault( self->dictionary, "flexuralRigidity", 3.0e22 );
    self->elasticPlateLength = Dictionary_GetFloat_WithDefault( self->dictionary, "elasticPlateLength", 1000000.0f );
    
    self->bedRockErosionLengthScale = Dictionary_GetFloat_WithDefault( self->dictionary,
        "bedRockErosionLengthScale", 5.0e3 );
    self->sedimentErosionLengthScale = Dictionary_GetFloat_WithDefault( self->dictionary,
        "sedimentErosionLengthScale", 1.0e3 );

    self->staticDt = Dictionary_GetDouble_WithDefault( self->dictionary, "staticDt", 2.0e0 );
    self->dumpEvery = Dictionary_GetDouble_WithDefault( self->dictionary, "dumpEvery", 500 );

    /* File output options */

    fileOutputDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "fileOutput" ) );
    if( fileOutputDict ){
        self->outputFormat = (!strcmp(Dictionary_GetString_WithDefault( fileOutputDict, "outputFormat", "text" ), "binary")?BINARY:TEXT);
        self->outputPath = Dictionary_GetString_WithDefault( fileOutputDict, "outputDirectory", "/." );
        self->itopography = Dictionary_GetBool_WithDefault( fileOutputDict, "itopography", True );
        self->ierosion = Dictionary_GetBool_WithDefault( fileOutputDict, "ierosion", False );
        self->iwater = Dictionary_GetBool_WithDefault( fileOutputDict, "iwater", False );
        self->iuplift = Dictionary_GetBool_WithDefault( fileOutputDict, "iuplift", False );
        self->isedimentation = Dictionary_GetBool_WithDefault( fileOutputDict, "isedimentation", False );
        self->iisostacy = Dictionary_GetBool_WithDefault( fileOutputDict, "iisostacy", False );
        self->idaughterNodes = Dictionary_GetBool_WithDefault( fileOutputDict, "idaughterNodes", False );
    }

    /* Add hooks to exisiting entry points */
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_Build ),
        "SPModelBuild",
        _SPModel_Context_Build,
        SPModel_Context_Type );
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_IC ),
        "SPModelIC",
        _SPModel_Context_InitialConditions,
        SPModel_Context_Type );
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_Dt ),
        "SPModelDt",
        _SPModel_Context_Dt,
        SPModel_Context_Type );
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_BC ),
        "SPModelBC",
        _SPModel_Context_BoundaryConditions,
        SPModel_Context_Type );
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_Solve ),
        "SPModelSolve",
        _SPModel_Context_Solve,
        SPModel_Context_Type );
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_Sync ),
        "SPModelSync",
        _SPModel_Context_Sync,
        SPModel_Context_Type );
    EntryPoint_Append(
        Context_GetEntryPoint( self, AbstractContext_EP_Dump ),
        "SPModelDump",
        _SPModel_Context_Dump,
        SPModel_Context_Type );
    
    SPModel_Simulation_ContextExtHandle = ExtensionManager_Add( self->extensionMgr, "Simulation",
                                        sizeof(SPModelSimulationContextExtension) );
}

void _SPModel_Context_Delete( void* context ) {
    SPModel_Context*    self = (SPModel_Context*)context;

    if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );
    
    Stg_Class_Delete( self->globalMesh );
    Stg_Class_Delete( self->localMesh );

    Stg_Class_Delete( self->catchmentList );
    Stg_Class_Delete( self->meshDecomp );
    Stg_Class_Delete( self->meshSmoother );
    Stg_Class_Delete( self->regularMesh );
    Stg_Class_Delete( self->interpolator );

    if ( self->cyclicBC ){
        Stg_Class_Delete( self->cyclicBC );
    }
    
    
    Stg_Class_Delete( self->seaLevel );
    
    /* Stg_Class_Delete parent class */
    _AbstractContext_Delete( context );
}



void _SPModel_Context_Print( void* context, Stream* stream ) {
    SPModel_Context* self = (SPModel_Context*)context;

    /* General info */
    Journal_Printf( stream, "SPModel_Context (%p):\n", self );

    /* Virtual info */

    /* SPModel_Context info */
    _Stg_Component_Print( self->globalMesh, stream );

    /* Parent class info */
    _AbstractContext_Print( context, stream );
}

void _SPModel_SimulationContextExtension_AllocateMemory( void *context, MemoryFlag mf )
{
    SPModelSimulationContextExtension *simulationExt = NULL;
    SPModel_Context* self = (SPModel_Context*)context;
    
    simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
                    ( self->extensionMgr, self, SPModel_Simulation_ContextExtHandle );

    assert( simulationExt );
    if( (self->rank == MASTER_PROC) && (mf == GLOBAL_EXTENSION) ){
        simulationExt->globalWater = malloc( sizeof( float ) * self->localMesh->numNodes );
        memset( simulationExt->globalWater, 0, sizeof( float ) * self->localMesh->numNodes );
        
        simulationExt->globalErosion = malloc( sizeof( float ) * self->localMesh->numNodes );
        memset( simulationExt->globalErosion, 0, sizeof( float ) * self->localMesh->numNodes );
        
        simulationExt->globalSedimentHistory = malloc( sizeof( float ) * self->localMesh->numNodes );
        memset( simulationExt->globalSedimentHistory, 0, sizeof( float ) * self->localMesh->numNodes );
        
        simulationExt->globalUplift = malloc( sizeof( float ) * self->localMesh->numNodes );
        memset( simulationExt->globalUplift, 0, sizeof( float ) * self->localMesh->numNodes );
    }

    simulationExt->water = malloc( sizeof( float ) * self->localMesh->myLoad );
    memset( simulationExt->water, 0, sizeof( float ) * self->localMesh->myLoad );
    
    simulationExt->erosion = malloc( sizeof( float ) * self->localMesh->myLoad );
    memset( simulationExt->erosion, 0, sizeof( float ) * self->localMesh->myLoad );
    
    simulationExt->sediment = malloc( sizeof( float ) * self->localMesh->myLoad );
    memset( simulationExt->sediment, 0, sizeof( float ) * self->localMesh->myLoad );
    
    simulationExt->sedimentHistory = malloc( sizeof( float ) * self->localMesh->myLoad );
    memset( simulationExt->sedimentHistory, 0, sizeof( float ) * self->localMesh->myLoad );
    
    simulationExt->uplift = malloc( sizeof( float ) * self->localMesh->myLoad );
    memset( simulationExt->uplift, 0, sizeof( float ) * self->localMesh->myLoad );
    
    simulationExt->hiso = malloc( sizeof( float ) * self->localMesh->numNodes );
    memset( simulationExt->hiso, 0, sizeof( float ) * self->localMesh->numNodes );
    
    simulationExt->hisoPrev = malloc( sizeof( float ) * self->localMesh->numNodes );
    memset( simulationExt->hisoPrev, 0, sizeof( float ) * self->localMesh->numNodes );
    
    simulationExt->hisotot = malloc( sizeof( float ) * self->localMesh->numNodes );
    memset( simulationExt->hisotot, 0, sizeof( float ) * self->localMesh->numNodes );
}

void _SPModel_SimulationContextExtension_ReleaseMemory( void *context, MemoryFlag mf )
{
    SPModelSimulationContextExtension *simulationExt = NULL;
    SPModel_Context* self = (SPModel_Context*)context;
    
    simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
                    ( self->extensionMgr, self, SPModel_Simulation_ContextExtHandle );

    assert( simulationExt );
    
    if( (self->rank == MASTER_PROC) && (mf == GLOBAL_EXTENSION) ){
        free( simulationExt->globalWater );
        free( simulationExt->globalErosion );
        free( simulationExt->globalSedimentHistory );
        free( simulationExt->globalUplift );
    }

    free( simulationExt->water );
    free( simulationExt->erosion );
    free( simulationExt->sediment );
    free( simulationExt->sedimentHistory );
    free( simulationExt->uplift );
    free( simulationExt->hiso );
    free( simulationExt->hisoPrev );
    free( simulationExt->hisotot );
}


void _SPModel_Context_Build( void* context ) {
    SPModel_Context*    self = (SPModel_Context*)context;
    int smoothMesh = 0;
    Dictionary *meshDict = NULL;
    
    if( self->rank == MASTER_PROC ) Journal_DPrintf( self->debug, "In: %s\n", __func__ );

    /*Building the parameters*/
    Stg_Component_Build( self->seaLevel, NULL, False );
    
    if( self->rank == MASTER_PROC ){
        
        Stg_Component_Build( self->globalMesh, NULL, False );
        if( self->cyclicBC ){
            Stg_Component_Build( self->cyclicBC, NULL, False );
            Stg_Component_Execute( self->cyclicBC, NULL, False );
        }
        Stg_Component_Build( self->catchmentList, NULL, False );
        Stg_Component_Build( self->meshSmoother, NULL, False );
        
#ifdef USE_DEBUG
        Print( self->interpolator, Journal_Register( "stream", Info_Type ) );
        Print( self->regularMesh, Journal_Register( "stream", Info_Type ) );
        Print( self->globalMesh, Journal_Register( Info_Type, "stream" ) );
        Print( self->catchmentList, Journal_Register( Info_Type, "catchmentStream" ) );
#endif
        
        Stg_Component_Build( self->meshDecomp, NULL, False );
        
        meshDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "mesh" ) );
        smoothMesh = Dictionary_GetBool_WithDefault( meshDict, "smoothMesh", False );

        if( smoothMesh ){
            Journal_Printf( Journal_Register( "stream", Info_Type ), "Smoothing Mesh\n" );
            Stg_Component_Execute( self->meshSmoother, self->globalMesh, False );       
        }
    }
    
    Stg_Component_Execute( (void*)self->meshDecomp, (void*)self->localMesh, False );
    
    Stg_Component_Build( self->regularMesh, NULL, False );
    Stg_Component_Build( self->interpolator, self->globalMesh, False );
    
    _SPModel_SimulationContextExtension_AllocateMemory( self, GLOBAL_EXTENSION );
}

void _SPModel_Context_InitialConditions( void* context )
{
    SPModel_Context*    self = (SPModel_Context*)context;

    self->dt = self->staticDt;
}

double _SPModel_Context_Dt( void* context ) {
    SPModel_Context*    self = (SPModel_Context*)context;

    if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );

    return self->dt;
}

void _SPModel_Context_SetDt( void* context, double dt ) {
    SPModel_Context*    self = (SPModel_Context*)context;

    if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );

    self->dt = dt;
}


void _SPModel_Context_BoundaryConditions( void* context ) 
{
    SPModel_Context *spmContext = NULL;
    
    assert( context );

    spmContext = (SPModel_Context*)context;
    
    SurfaceMesh_BoundaryConditions( spmContext->localMesh );
}

void _SPModel_Context_Solve( void* context )
{
    int globalMaxComm = 0;
    int localMaxComm = 0;
    SPModel_Context*        self = (SPModel_Context*)context;
    
    ParameterTimeSeries_Interpolate( self->seaLevel, self );
    
        /* Perform the SPModel solve loop */
    SurfaceMesh_OrderNodes( self->localMesh );

    localMaxComm = (self->itemsSent > self->itemsReceived)?self->itemsSent:self->itemsReceived;
    MPI_Allreduce( &(localMaxComm), &globalMaxComm, 1, MPI_INT, MPI_MAX, self->communicator );

#ifdef DEBUG
    printf( "globalMaxComm -> %d\n", globalMaxComm );
#endif

    self->redistributeNodesFlag = False;
    
    /*TODO dynamically decide to redistribute */
    /*if( globalMaxComm > 50 ){
        _SPModel_Context_RedistributeNodes( self );
        self->redistributeNodesFlag = True;
    }*/
}

void _SPModel_Context_Sync( void* context )
{
    SPModel_Context *self = NULL;
    _SurfaceMeshDecomp *meshDecomp = NULL;
    
    assert( context );
    
    self = (SPModel_Context*)context;

    meshDecomp = self->meshDecomp;

    assert( meshDecomp );
    
    meshDecomp->syncMesh( meshDecomp, self->localMesh );
    
    SurfaceMesh_UpdateFlows( self->localMesh );
}

void _SPModel_Context_Gather( SPModel_Context* self )
{
    typedef struct Package_t{
        float h;
        int id, receiver;
    }Package; 

    int i, j;
    Package package;
    MPI_Status status;

    SPModel_Context *context;
    SurfaceMesh *localMesh, *globalMesh;
    _SurfaceMeshDecomp *meshDecomp;

    assert( self );
    
    context = (SPModel_Context*)self;
    
    localMesh = self->localMesh;
    globalMesh = self->globalMesh;
    meshDecomp = context->meshDecomp;
    
    if( self->rank == MASTER_PROC ){
        for( i=0; i<context->nproc; i++ ){
            for( j=0; j<meshDecomp->processorLoad[i]; j++ ){
                if( i == MASTER_PROC ){
                    assert( globalMesh->id[localMesh->id[j]-1] == localMesh->id[j] );
                    
                    globalMesh->h[localMesh->id[j]-1] = localMesh->h[j];
                    globalMesh->receiver[localMesh->id[j]-1] = localMesh->receiver[j];
                }
                else{
                    MPI_Recv( &package, sizeof(Package), MPI_BYTE, i,
                        PACKAGE_TAG, MPI_COMM_WORLD, &status );
                    
                    assert( globalMesh->id[package.id-1] == package.id );

                    globalMesh->h[package.id-1] = package.h;
                    globalMesh->receiver[package.id-1] = package.receiver;
                }
            }
        }
    }
    else{
        for( i=0; i<localMesh->myLoad; i++ ){
            package.id = localMesh->id[i];
            package.h = localMesh->h[i];
            package.receiver = localMesh->receiver[i];

            MPI_Send( &package, sizeof(Package), MPI_BYTE, MASTER_PROC,
                    PACKAGE_TAG, MPI_COMM_WORLD );
        }
    }
}

void _SPModel_Context_RedistributeNodes( SPModel_Context *context )
{
    char *interpolationOption;
    
    assert( context );
    
    _SPModel_Context_Gather( context );
    SPModelSimulationContextExtension_Gather( context );

    Stg_Class_Delete( context->meshDecomp );

    if( context->rank == MASTER_PROC ){
        SurfaceMesh_BuildRiverNetwork( context->globalMesh );
        
        Stg_Class_Delete( context->catchmentList );
        
        context->catchmentList = CatchmentList_New( "catchments", context->globalMesh, context->nproc, context->dictionary );
        
        Stg_Component_Build( context->catchmentList, NULL, False );
        
        context->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshIrregularDecomp_New( "meshDecomp", 
                            context->catchmentList, context->globalMesh, context->dictionary );
        
        Stg_Component_Build( context->meshDecomp, NULL, False );
    }
    else{
        context->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshIrregularDecomp_New( "meshDecomp", 
                                context->catchmentList, context->globalMesh, context->dictionary );
        
    }
    
    SurfaceMesh_ReleaseMemory( context->localMesh );
    Stg_Component_Execute( (void*)context->meshDecomp, (void*)context->localMesh, False );
    
    Stg_Class_Delete( context->interpolator );

    interpolationOption = Dictionary_GetString_WithDefault( context->dictionary, "interpolation", "Linear" );

    if( !strcmp(interpolationOption, "Linear") ){
        context->interpolator = (_Interpolator*)LinearInterpolator_New( "linearInterpolator", context->dictionary, 
                                            context->localMesh, context->regularMesh );
    }
    else if ( !strcmp(interpolationOption, "Cubic") ){
        context->interpolator = (_Interpolator*)LinearInterpolator_New( "cubicInterpolator", context->dictionary, 
                                            context->localMesh, context->regularMesh );
    }

    Stg_Component_Build( context->interpolator, context->globalMesh, False );
    
    _SPModel_SimulationContextExtension_ReleaseMemory( context, LOCAL_EXTENSION );
    _SPModel_SimulationContextExtension_AllocateMemory( context, LOCAL_EXTENSION );
    
    SPModelSimulationContextExtension_Scatter( context );
}

void _SPModel_Context_Dump( SPModel_Context* self )
{
    SurfaceMesh *mesh = NULL;
    int i = 0;
    
    assert( self );

    _SPModel_Context_Gather( self );
    
    if( self->rank != MASTER_PROC ){
        return;
    }

    mesh = self->globalMesh;

    if( self->itopography ){
        char fileName[500] = {0};
        FILE *topo = NULL;
        
        if( self->outputFormat == TEXT ){
            memset( fileName, 0, sizeof(fileName) );
            sprintf( fileName, "%s/topography%d.txt", self->outputPath, (int)self->currentTime );
            topo = fopen( fileName, "w+" );
            assert( topo );

            fprintf( topo, "Topography\n" );
            fprintf( topo, "Column Description\n" );
            fprintf( topo, "x\ty\th\tboundaryCondition\n" );
            for( i=0; i<mesh->numNodes; i++ ){
                fprintf( topo, "%lf %lf %lf %lf\n", mesh->x[i],
                        mesh->y[i], mesh->h[i], mesh->boundaryConditions[i] );
            }
            fclose( topo );
        }
        else if( self->outputFormat == BINARY ){
            typedef struct outputPackage_t{
                float x, y, h, bc;
            }outputPackage;

            outputPackage output;
            
            memset( fileName, 0, sizeof(fileName) );
            sprintf( fileName, "%s/topography%d", self->outputPath, (int)self->currentTime );
            topo = fopen( fileName, "w+" );
            assert( topo );
            
            for( i=0; i<mesh->numNodes; i++ ){
                output.x = mesh->x[i];
                output.y = mesh->y[i];
                output.h = mesh->h[i];
                output.bc = mesh->boundaryConditions[i];
                fwrite( &output, sizeof( outputPackage ), 1, topo );
            }
            fclose( topo );
        }
    }
    
    if( self->idaughterNodes ){
        char fileName[500] = {0};
        FILE *topo = NULL;
        int j = 0, k = 0, value;

        if( self->outputFormat == TEXT ){
            memset( fileName, 0, sizeof(fileName) );
            sprintf( fileName, "%s/daughterNodes%d.txt", self->outputPath, (int)self->currentTime );
            topo = fopen( fileName, "w+" );
            assert( topo );

            fprintf( topo, "Daughter Nodes\n" );
            fprintf( topo, "Column Description\n" );
            fprintf( topo, "Immediate Daughter Node of ith node(-1 = boundary node, 0 = neighbour of itself)\n" );
            for( i=0; i<mesh->numNodes; i++ ){
                if( !(mesh->boundaryConditions[mesh->id[i]-1]) ){
                    value = -1;
                    if( mesh->receiver[mesh->id[i]-1] != mesh->id[i] ){
                        value = mesh->receiver[mesh->id[i]-1];
                    }
                    fprintf( topo, "%d\n", value );
                }
                else{
                    if( mesh->receiver[mesh->id[i]-1] == mesh->id[i] ){
                        value = 0;
                        fprintf( topo, "%d\n", value );
                    }
                    else{
                        k = mesh->id[i];
                        j = mesh->receiver[mesh->id[i]-1];
                        value = j;
                        fprintf( topo, "%d\n", value );
                    }
                }
            }
            fclose( topo );
        }
        else if( self->outputFormat == BINARY ){
            memset( fileName, 0, sizeof(fileName) );
            sprintf( fileName, "%s/daughterNodes%d", self->outputPath, (int)self->currentTime );
            topo = fopen( fileName, "w+" );
            assert( topo );

            for( i=0; i<mesh->numNodes; i++ ){
                if( !(mesh->boundaryConditions[mesh->id[i]-1]) ){
                    value = -1;
                    if( mesh->receiver[mesh->id[i]-1] != mesh->id[i] ){
                        value = mesh->receiver[mesh->id[i]-1];
                    }
                    fwrite( &value, sizeof( int ), 1, topo );
                }
                else{
                    if( mesh->receiver[mesh->id[i]-1] == mesh->id[i] ){
                        value = 0;
                        fwrite( &value, sizeof( int ), 1, topo );
                    }
                    else{
                        k = mesh->id[i];
                        j = mesh->receiver[mesh->id[i]-1];
                        value = j;
                        fwrite( &value, sizeof( int ), 1, topo );
                    }
                }
            }
            fclose( topo );
        }
    }
}

void _SPModel_Context_ComponentConstruct( void* context, Stg_ComponentFactory* cf )
{
}

void _SPModel_Context_ComponentBuild( void* context, void* data )
{
}

void _SPModel_Context_ComponentInitialise( void* context, void* data )
{
}

void _SPModel_Context_ComponentExecute( void* context, void* data )
{
}

void _SPModel_Context_ComponentDestroy( void* context, void* data )
{
}

SPModel documentation generated by doxygen at Fri Mar 24 01:41:45 2006