SPModel: libSPModel/src/Context.c Source File
VPAC - Computational Software Development
Main | SPModel | StGermain FrameWork |
Main Page | Alphabetical List | Class List | Directories | File List | Class Members | File Members

Context.c

Go to the documentation of this file.
00001 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00002 **
00003 ** Copyright (C), 2004, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
00004 **
00005 ** Authors:
00006 **  Ogar R. Widjaja, Computational Scientist, VPAC.
00007 **  Raquibul Hassan, Software Engineer, VPAC. (raq@vpac.org)
00008 **  Keith Hsuan, Computational Scientist, VPAC (keith@vpac.org)
00009 **  William F. Appelbe, Director, VPAC. (bill@vpac.org)
00010 **  Stevan M. Quenette, Senior Software Engineer, VPAC. (steve@vpac.org)
00011 **  Patrick D. Sunter, Software Engineer, VPAC. (patrick@vpac.org)
00012 **
00013 ** This file may be distributed under the terms of the VPAC Public License
00014 ** as defined by VPAC of Australia and appearing in the file
00015 ** LICENSE.VPL included in the packaging of this file.
00016 **
00017 ** This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
00018 ** WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00019 **
00020 ** $Id: Context.c 285 2006-03-20 07:26:25Z RaquibulHassan $
00021 **
00022 **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
00023 
00024 #include <mpi.h>
00025 #include <StGermain/StGermain.h>
00026 
00027 #include "types.h"
00028 #include "Context.h"
00029 
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <string.h>
00033 #include <limits.h>
00034 #include <assert.h>
00035 #include "SurfaceMesh.h"
00036 #include "SurfaceRegularMesh.h"
00037 #include "SurfaceMeshDecomp.h"
00038 #include "SurfaceMeshIrregularDecomp.h"
00039 #include "SurfaceMeshRegularDecomp.h"
00040 #include "Catchment.h"
00041 #include "_Interpolator.h"
00042 #include "LinearInterpolator.h"
00043 #include "SplineInterpolator.h"
00044 #include "SurfaceMeshSmoother.h"
00045 #include "ParameterTimeSeries.h"
00046 #include "SurfaceMeshCyclicBC.h"
00047 #include "SurfaceMeshRectangularCyclicBC.h"
00048 #include "CommHandler.h"
00049 
00050 /* Textual name of this class */
00051 const Type SPModel_Context_Type = "SPModel_Context";
00052 
00053 ExtensionInfo_Index SPModel_Simulation_ContextExtHandle;
00054 
00055 SPModel_Context* SPModel_Context_DefaultNew( Name name )
00056 {
00057     return _SPModel_Context_New( sizeof(SPModel_Context), SPModel_Context_Type,
00058             _SPModel_Context_Delete, _SPModel_Context_Print, NULL, (void*)SPModel_Context_DefaultNew,
00059             _SPModel_Context_ComponentConstruct, _SPModel_Context_ComponentBuild, _SPModel_Context_ComponentInitialise,
00060             _SPModel_Context_ComponentExecute, _SPModel_Context_ComponentDestroy, name, False,
00061         _SPModel_Context_SetDt, 0, 0, 0, NULL );
00062 }
00063 
00064 SPModel_Context* SPModel_Context_New(
00065         Name                    name,
00066         double                  startTime,
00067         double                  stopTime,
00068         MPI_Comm                communicator,
00069         Dictionary*             dictionary )
00070 {
00071     return _SPModel_Context_New( sizeof(SPModel_Context), SPModel_Context_Type,
00072             _SPModel_Context_Delete, _SPModel_Context_Print, NULL, (void*)SPModel_Context_DefaultNew,
00073             _SPModel_Context_ComponentConstruct, _SPModel_Context_ComponentBuild, _SPModel_Context_ComponentInitialise,
00074             _SPModel_Context_ComponentExecute, _SPModel_Context_ComponentDestroy, name, True,
00075         _SPModel_Context_SetDt, startTime, stopTime, communicator, dictionary );
00076 }
00077 
00078 
00079 SPModel_Context* _SPModel_Context_New(
00080         SizeT                   _sizeOfSelf,
00081         Type                    type,
00082         Stg_Class_DeleteFunction*               _delete,
00083         Stg_Class_PrintFunction*                _print,
00084         Stg_Class_CopyFunction*             _copy, 
00085         Stg_Component_DefaultConstructorFunction*   _defaultConstructor,
00086         Stg_Component_ConstructFunction*            _construct,
00087         Stg_Component_BuildFunction*        _build,
00088         Stg_Component_InitialiseFunction*       _initialise,
00089         Stg_Component_ExecuteFunction*      _execute,
00090         Stg_Component_DestroyFunction*      _destroy,
00091         Name                            name,
00092         Bool                            initFlag,
00093         AbstractContext_SetDt*          _setDt,
00094         double                  startTime,
00095         double                  stopTime,
00096         MPI_Comm                communicator,
00097         Dictionary*             dictionary )
00098 {
00099     SPModel_Context*    self;
00100 
00101     /* Allocate memory */
00102     self = (SPModel_Context*)_AbstractContext_New( _sizeOfSelf, type, _delete, _print, NULL, 
00103             _defaultConstructor, _construct, _build, _initialise, _execute, _destroy, name, 
00104             initFlag, _setDt, startTime, stopTime, communicator, dictionary );
00105 
00106     /* General info */
00107 
00108     /* Virtual info */
00109 
00110     /* SPModel_Context info */
00111     if( initFlag ){
00112         _SPModel_Context_Init( self );
00113     }
00114 
00115     return self;
00116 }
00117 
00118 
00119 void _SPModel_Context_Init( SPModel_Context* self ) {
00120 
00121     Dictionary *fileOutputDict = NULL;
00122     Dictionary *cyclicBCDict = NULL;
00123     char *interpolationOption = NULL;
00124     
00125     if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );
00126 
00127     /* construct the mesh */
00128     self->globalMesh = SurfaceMesh_New( "globalMesh", self->nproc, self->rank, self->extensionMgr_Register, self->dictionary );
00129     self->localMesh = SurfaceMesh_New( "localMesh", self->nproc, self->rank, self->extensionMgr_Register, self->dictionary );
00130     self->catchmentList = CatchmentList_New( "catchments", self->globalMesh, self->nproc, self->dictionary );
00131     self->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshIrregularDecomp_New( "meshDecomp", self->catchmentList, self->globalMesh, self->dictionary );
00132     /*self->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshRegularDecomp_New( "meshDecomp", self->globalMesh, self->dictionary );*/
00133     self->meshSmoother = SurfaceMeshSmoother_New( "surfaceMeshSmoother", self->dictionary );
00134     self->regularMesh = SurfaceRegularMesh_New( "regularMesh", self->dictionary, self->localMesh );
00135 
00136     interpolationOption = Dictionary_GetString_WithDefault( self->dictionary, "interpolation", "Linear" );
00137 
00138     if( !strcmp(interpolationOption, "Linear") ){
00139         self->interpolator = (_Interpolator*)LinearInterpolator_New( "linearInterpolator", self->dictionary, self->localMesh, self->regularMesh );
00140     }
00141     else if ( !strcmp(interpolationOption, "Cubic") ){
00142         self->interpolator = (_Interpolator*)SplineInterpolator_New( "splineInterpolator", self->dictionary, self->localMesh, self->regularMesh );      
00143     }
00144 
00145     if ( (cyclicBCDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "CyclicBoundaryConditions" ) ))){
00146         self->cyclicBC = (_SurfaceMeshCyclicBC*)SurfaceMeshRectangularCyclicBC_New
00147                             ( "cyclicBC", self->globalMesh, self->dictionary );
00148     }
00149     
00150     /* defaults of physical parameters */
00151     self->gravity = Dictionary_GetDouble_WithDefault( self->dictionary, "gravity", 9.81f );
00152     self->flexureInterval = Dictionary_GetDouble_WithDefault( self->dictionary, "flexureInterval", 100 );
00153     
00154     self->seaLevel = ParameterTimeSeries_New( "seaLevel", 0.0f, self->dictionary );
00155     
00156     self->upliftRate = Dictionary_GetFloat_WithDefault( self->dictionary, "upliftRate", 1.0e-5 );
00157     
00158     self->oroLength = Dictionary_GetFloat_WithDefault( self->dictionary, "oroLength", 0 );
00159     self->oroRate = Dictionary_GetFloat_WithDefault( self->dictionary, "oroRate", 5.0e-1 );
00160     self->fluvialConstant = Dictionary_GetFloat_WithDefault( self->dictionary, "fluvialConstant",3.0e-2 );
00161     
00162     self->bedRockDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "bedRockDiffusivity", 1.0e-3 );
00163     self->sedimentDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "sedimentDiffusivity", 1.0e3 );
00164     self->submarineBedRockDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "submarineBedRockDiffusivity", 1.0e-3 );
00165     self->submarineSedimentDiffusivity = Dictionary_GetFloat_WithDefault( self->dictionary, "submarineSedimentDiffusivity", 1.0e3 );
00166     
00167     self->bedRockDensity = Dictionary_GetFloat_WithDefault( self->dictionary, "bedRockDensity", 2.8e3 );
00168     self->sedimentDensity = Dictionary_GetFloat_WithDefault( self->dictionary, "sedimentDensity", 2.8e3 );
00169     self->asthenosphereDensity = Dictionary_GetFloat_WithDefault( self->dictionary, "asthenosphereDensity", 3.3e3 );
00170     
00171     self->flexuralRigidity = Dictionary_GetFloat_WithDefault( self->dictionary, "flexuralRigidity", 3.0e22 );
00172     self->elasticPlateLength = Dictionary_GetFloat_WithDefault( self->dictionary, "elasticPlateLength", 1000000.0f );
00173     
00174     self->bedRockErosionLengthScale = Dictionary_GetFloat_WithDefault( self->dictionary,
00175         "bedRockErosionLengthScale", 5.0e3 );
00176     self->sedimentErosionLengthScale = Dictionary_GetFloat_WithDefault( self->dictionary,
00177         "sedimentErosionLengthScale", 1.0e3 );
00178 
00179     self->staticDt = Dictionary_GetDouble_WithDefault( self->dictionary, "staticDt", 2.0e0 );
00180     self->dumpEvery = Dictionary_GetDouble_WithDefault( self->dictionary, "dumpEvery", 500 );
00181 
00182     /* File output options */
00183 
00184     fileOutputDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "fileOutput" ) );
00185     if( fileOutputDict ){
00186         self->outputFormat = (!strcmp(Dictionary_GetString_WithDefault( fileOutputDict, "outputFormat", "text" ), "binary")?BINARY:TEXT);
00187         self->outputPath = Dictionary_GetString_WithDefault( fileOutputDict, "outputDirectory", "/." );
00188         self->itopography = Dictionary_GetBool_WithDefault( fileOutputDict, "itopography", True );
00189         self->ierosion = Dictionary_GetBool_WithDefault( fileOutputDict, "ierosion", False );
00190         self->iwater = Dictionary_GetBool_WithDefault( fileOutputDict, "iwater", False );
00191         self->iuplift = Dictionary_GetBool_WithDefault( fileOutputDict, "iuplift", False );
00192         self->isedimentation = Dictionary_GetBool_WithDefault( fileOutputDict, "isedimentation", False );
00193         self->iisostacy = Dictionary_GetBool_WithDefault( fileOutputDict, "iisostacy", False );
00194         self->idaughterNodes = Dictionary_GetBool_WithDefault( fileOutputDict, "idaughterNodes", False );
00195     }
00196 
00197     /* Add hooks to exisiting entry points */
00198     EntryPoint_Append(
00199         Context_GetEntryPoint( self, AbstractContext_EP_Build ),
00200         "SPModelBuild",
00201         _SPModel_Context_Build,
00202         SPModel_Context_Type );
00203     EntryPoint_Append(
00204         Context_GetEntryPoint( self, AbstractContext_EP_IC ),
00205         "SPModelIC",
00206         _SPModel_Context_InitialConditions,
00207         SPModel_Context_Type );
00208     EntryPoint_Append(
00209         Context_GetEntryPoint( self, AbstractContext_EP_Dt ),
00210         "SPModelDt",
00211         _SPModel_Context_Dt,
00212         SPModel_Context_Type );
00213     EntryPoint_Append(
00214         Context_GetEntryPoint( self, AbstractContext_EP_BC ),
00215         "SPModelBC",
00216         _SPModel_Context_BoundaryConditions,
00217         SPModel_Context_Type );
00218     EntryPoint_Append(
00219         Context_GetEntryPoint( self, AbstractContext_EP_Solve ),
00220         "SPModelSolve",
00221         _SPModel_Context_Solve,
00222         SPModel_Context_Type );
00223     EntryPoint_Append(
00224         Context_GetEntryPoint( self, AbstractContext_EP_Sync ),
00225         "SPModelSync",
00226         _SPModel_Context_Sync,
00227         SPModel_Context_Type );
00228     EntryPoint_Append(
00229         Context_GetEntryPoint( self, AbstractContext_EP_Dump ),
00230         "SPModelDump",
00231         _SPModel_Context_Dump,
00232         SPModel_Context_Type );
00233     
00234     SPModel_Simulation_ContextExtHandle = ExtensionManager_Add( self->extensionMgr, "Simulation",
00235                                         sizeof(SPModelSimulationContextExtension) );
00236 }
00237 
00238 void _SPModel_Context_Delete( void* context ) {
00239     SPModel_Context*    self = (SPModel_Context*)context;
00240 
00241     if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );
00242     
00243     Stg_Class_Delete( self->globalMesh );
00244     Stg_Class_Delete( self->localMesh );
00245 
00246     Stg_Class_Delete( self->catchmentList );
00247     Stg_Class_Delete( self->meshDecomp );
00248     Stg_Class_Delete( self->meshSmoother );
00249     Stg_Class_Delete( self->regularMesh );
00250     Stg_Class_Delete( self->interpolator );
00251 
00252     if ( self->cyclicBC ){
00253         Stg_Class_Delete( self->cyclicBC );
00254     }
00255     
00256     
00257     Stg_Class_Delete( self->seaLevel );
00258     
00259     /* Stg_Class_Delete parent class */
00260     _AbstractContext_Delete( context );
00261 }
00262 
00263 
00264 
00265 void _SPModel_Context_Print( void* context, Stream* stream ) {
00266     SPModel_Context* self = (SPModel_Context*)context;
00267 
00268     /* General info */
00269     Journal_Printf( stream, "SPModel_Context (%p):\n", self );
00270 
00271     /* Virtual info */
00272 
00273     /* SPModel_Context info */
00274     _Stg_Component_Print( self->globalMesh, stream );
00275 
00276     /* Parent class info */
00277     _AbstractContext_Print( context, stream );
00278 }
00279 
00280 void _SPModel_SimulationContextExtension_AllocateMemory( void *context, MemoryFlag mf )
00281 {
00282     SPModelSimulationContextExtension *simulationExt = NULL;
00283     SPModel_Context* self = (SPModel_Context*)context;
00284     
00285     simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
00286                     ( self->extensionMgr, self, SPModel_Simulation_ContextExtHandle );
00287 
00288     assert( simulationExt );
00289     if( (self->rank == MASTER_PROC) && (mf == GLOBAL_EXTENSION) ){
00290         simulationExt->globalWater = malloc( sizeof( float ) * self->localMesh->numNodes );
00291         memset( simulationExt->globalWater, 0, sizeof( float ) * self->localMesh->numNodes );
00292         
00293         simulationExt->globalErosion = malloc( sizeof( float ) * self->localMesh->numNodes );
00294         memset( simulationExt->globalErosion, 0, sizeof( float ) * self->localMesh->numNodes );
00295         
00296         simulationExt->globalSedimentHistory = malloc( sizeof( float ) * self->localMesh->numNodes );
00297         memset( simulationExt->globalSedimentHistory, 0, sizeof( float ) * self->localMesh->numNodes );
00298         
00299         simulationExt->globalUplift = malloc( sizeof( float ) * self->localMesh->numNodes );
00300         memset( simulationExt->globalUplift, 0, sizeof( float ) * self->localMesh->numNodes );
00301     }
00302 
00303     simulationExt->water = malloc( sizeof( float ) * self->localMesh->myLoad );
00304     memset( simulationExt->water, 0, sizeof( float ) * self->localMesh->myLoad );
00305     
00306     simulationExt->erosion = malloc( sizeof( float ) * self->localMesh->myLoad );
00307     memset( simulationExt->erosion, 0, sizeof( float ) * self->localMesh->myLoad );
00308     
00309     simulationExt->sediment = malloc( sizeof( float ) * self->localMesh->myLoad );
00310     memset( simulationExt->sediment, 0, sizeof( float ) * self->localMesh->myLoad );
00311     
00312     simulationExt->sedimentHistory = malloc( sizeof( float ) * self->localMesh->myLoad );
00313     memset( simulationExt->sedimentHistory, 0, sizeof( float ) * self->localMesh->myLoad );
00314     
00315     simulationExt->uplift = malloc( sizeof( float ) * self->localMesh->myLoad );
00316     memset( simulationExt->uplift, 0, sizeof( float ) * self->localMesh->myLoad );
00317     
00318     simulationExt->hiso = malloc( sizeof( float ) * self->localMesh->numNodes );
00319     memset( simulationExt->hiso, 0, sizeof( float ) * self->localMesh->numNodes );
00320     
00321     simulationExt->hisoPrev = malloc( sizeof( float ) * self->localMesh->numNodes );
00322     memset( simulationExt->hisoPrev, 0, sizeof( float ) * self->localMesh->numNodes );
00323     
00324     simulationExt->hisotot = malloc( sizeof( float ) * self->localMesh->numNodes );
00325     memset( simulationExt->hisotot, 0, sizeof( float ) * self->localMesh->numNodes );
00326 }
00327 
00328 void _SPModel_SimulationContextExtension_ReleaseMemory( void *context, MemoryFlag mf )
00329 {
00330     SPModelSimulationContextExtension *simulationExt = NULL;
00331     SPModel_Context* self = (SPModel_Context*)context;
00332     
00333     simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
00334                     ( self->extensionMgr, self, SPModel_Simulation_ContextExtHandle );
00335 
00336     assert( simulationExt );
00337     
00338     if( (self->rank == MASTER_PROC) && (mf == GLOBAL_EXTENSION) ){
00339         free( simulationExt->globalWater );
00340         free( simulationExt->globalErosion );
00341         free( simulationExt->globalSedimentHistory );
00342         free( simulationExt->globalUplift );
00343     }
00344 
00345     free( simulationExt->water );
00346     free( simulationExt->erosion );
00347     free( simulationExt->sediment );
00348     free( simulationExt->sedimentHistory );
00349     free( simulationExt->uplift );
00350     free( simulationExt->hiso );
00351     free( simulationExt->hisoPrev );
00352     free( simulationExt->hisotot );
00353 }
00354 
00355 
00356 void _SPModel_Context_Build( void* context ) {
00357     SPModel_Context*    self = (SPModel_Context*)context;
00358     int smoothMesh = 0;
00359     Dictionary *meshDict = NULL;
00360     
00361     if( self->rank == MASTER_PROC ) Journal_DPrintf( self->debug, "In: %s\n", __func__ );
00362 
00363     /*Building the parameters*/
00364     Stg_Component_Build( self->seaLevel, NULL, False );
00365     
00366     if( self->rank == MASTER_PROC ){
00367         
00368         Stg_Component_Build( self->globalMesh, NULL, False );
00369         if( self->cyclicBC ){
00370             Stg_Component_Build( self->cyclicBC, NULL, False );
00371             Stg_Component_Execute( self->cyclicBC, NULL, False );
00372         }
00373         Stg_Component_Build( self->catchmentList, NULL, False );
00374         Stg_Component_Build( self->meshSmoother, NULL, False );
00375         
00376 #ifdef USE_DEBUG
00377         Print( self->interpolator, Journal_Register( "stream", Info_Type ) );
00378         Print( self->regularMesh, Journal_Register( "stream", Info_Type ) );
00379         Print( self->globalMesh, Journal_Register( Info_Type, "stream" ) );
00380         Print( self->catchmentList, Journal_Register( Info_Type, "catchmentStream" ) );
00381 #endif
00382         
00383         Stg_Component_Build( self->meshDecomp, NULL, False );
00384         
00385         meshDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "mesh" ) );
00386         smoothMesh = Dictionary_GetBool_WithDefault( meshDict, "smoothMesh", False );
00387 
00388         if( smoothMesh ){
00389             Journal_Printf( Journal_Register( "stream", Info_Type ), "Smoothing Mesh\n" );
00390             Stg_Component_Execute( self->meshSmoother, self->globalMesh, False );       
00391         }
00392     }
00393     
00394     Stg_Component_Execute( (void*)self->meshDecomp, (void*)self->localMesh, False );
00395     
00396     Stg_Component_Build( self->regularMesh, NULL, False );
00397     Stg_Component_Build( self->interpolator, self->globalMesh, False );
00398     
00399     _SPModel_SimulationContextExtension_AllocateMemory( self, GLOBAL_EXTENSION );
00400 }
00401 
00402 void _SPModel_Context_InitialConditions( void* context )
00403 {
00404     SPModel_Context*    self = (SPModel_Context*)context;
00405 
00406     self->dt = self->staticDt;
00407 }
00408 
00409 double _SPModel_Context_Dt( void* context ) {
00410     SPModel_Context*    self = (SPModel_Context*)context;
00411 
00412     if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );
00413 
00414     return self->dt;
00415 }
00416 
00417 void _SPModel_Context_SetDt( void* context, double dt ) {
00418     SPModel_Context*    self = (SPModel_Context*)context;
00419 
00420     if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );
00421 
00422     self->dt = dt;
00423 }
00424 
00425 
00426 void _SPModel_Context_BoundaryConditions( void* context ) 
00427 {
00428     SPModel_Context *spmContext = NULL;
00429     
00430     assert( context );
00431 
00432     spmContext = (SPModel_Context*)context;
00433     
00434     SurfaceMesh_BoundaryConditions( spmContext->localMesh );
00435 }
00436 
00437 void _SPModel_Context_Solve( void* context )
00438 {
00439     int globalMaxComm = 0;
00440     int localMaxComm = 0;
00441     SPModel_Context*        self = (SPModel_Context*)context;
00442     
00443     ParameterTimeSeries_Interpolate( self->seaLevel, self );
00444     
00445         /* Perform the SPModel solve loop */
00446     SurfaceMesh_OrderNodes( self->localMesh );
00447 
00448     localMaxComm = (self->itemsSent > self->itemsReceived)?self->itemsSent:self->itemsReceived;
00449     MPI_Allreduce( &(localMaxComm), &globalMaxComm, 1, MPI_INT, MPI_MAX, self->communicator );
00450 
00451 #ifdef DEBUG
00452     printf( "globalMaxComm -> %d\n", globalMaxComm );
00453 #endif
00454 
00455     self->redistributeNodesFlag = False;
00456     
00457     /*TODO dynamically decide to redistribute */
00458     /*if( globalMaxComm > 50 ){
00459         _SPModel_Context_RedistributeNodes( self );
00460         self->redistributeNodesFlag = True;
00461     }*/
00462 }
00463 
00464 void _SPModel_Context_Sync( void* context )
00465 {
00466     SPModel_Context *self = NULL;
00467     _SurfaceMeshDecomp *meshDecomp = NULL;
00468     
00469     assert( context );
00470     
00471     self = (SPModel_Context*)context;
00472 
00473     meshDecomp = self->meshDecomp;
00474 
00475     assert( meshDecomp );
00476     
00477     meshDecomp->syncMesh( meshDecomp, self->localMesh );
00478     
00479     SurfaceMesh_UpdateFlows( self->localMesh );
00480 }
00481 
00482 void _SPModel_Context_Gather( SPModel_Context* self )
00483 {
00484     typedef struct Package_t{
00485         float h;
00486         int id, receiver;
00487     }Package; 
00488 
00489     int i, j;
00490     Package package;
00491     MPI_Status status;
00492 
00493     SPModel_Context *context;
00494     SurfaceMesh *localMesh, *globalMesh;
00495     _SurfaceMeshDecomp *meshDecomp;
00496 
00497     assert( self );
00498     
00499     context = (SPModel_Context*)self;
00500     
00501     localMesh = self->localMesh;
00502     globalMesh = self->globalMesh;
00503     meshDecomp = context->meshDecomp;
00504     
00505     if( self->rank == MASTER_PROC ){
00506         for( i=0; i<context->nproc; i++ ){
00507             for( j=0; j<meshDecomp->processorLoad[i]; j++ ){
00508                 if( i == MASTER_PROC ){
00509                     assert( globalMesh->id[localMesh->id[j]-1] == localMesh->id[j] );
00510                     
00511                     globalMesh->h[localMesh->id[j]-1] = localMesh->h[j];
00512                     globalMesh->receiver[localMesh->id[j]-1] = localMesh->receiver[j];
00513                 }
00514                 else{
00515                     MPI_Recv( &package, sizeof(Package), MPI_BYTE, i,
00516                         PACKAGE_TAG, MPI_COMM_WORLD, &status );
00517                     
00518                     assert( globalMesh->id[package.id-1] == package.id );
00519 
00520                     globalMesh->h[package.id-1] = package.h;
00521                     globalMesh->receiver[package.id-1] = package.receiver;
00522                 }
00523             }
00524         }
00525     }
00526     else{
00527         for( i=0; i<localMesh->myLoad; i++ ){
00528             package.id = localMesh->id[i];
00529             package.h = localMesh->h[i];
00530             package.receiver = localMesh->receiver[i];
00531 
00532             MPI_Send( &package, sizeof(Package), MPI_BYTE, MASTER_PROC,
00533                     PACKAGE_TAG, MPI_COMM_WORLD );
00534         }
00535     }
00536 }
00537 
00538 void _SPModel_Context_RedistributeNodes( SPModel_Context *context )
00539 {
00540     char *interpolationOption;
00541     
00542     assert( context );
00543     
00544     _SPModel_Context_Gather( context );
00545     SPModelSimulationContextExtension_Gather( context );
00546 
00547     Stg_Class_Delete( context->meshDecomp );
00548 
00549     if( context->rank == MASTER_PROC ){
00550         SurfaceMesh_BuildRiverNetwork( context->globalMesh );
00551         
00552         Stg_Class_Delete( context->catchmentList );
00553         
00554         context->catchmentList = CatchmentList_New( "catchments", context->globalMesh, context->nproc, context->dictionary );
00555         
00556         Stg_Component_Build( context->catchmentList, NULL, False );
00557         
00558         context->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshIrregularDecomp_New( "meshDecomp", 
00559                             context->catchmentList, context->globalMesh, context->dictionary );
00560         
00561         Stg_Component_Build( context->meshDecomp, NULL, False );
00562     }
00563     else{
00564         context->meshDecomp = (_SurfaceMeshDecomp*)SurfaceMeshIrregularDecomp_New( "meshDecomp", 
00565                                 context->catchmentList, context->globalMesh, context->dictionary );
00566         
00567     }
00568     
00569     SurfaceMesh_ReleaseMemory( context->localMesh );
00570     Stg_Component_Execute( (void*)context->meshDecomp, (void*)context->localMesh, False );
00571     
00572     Stg_Class_Delete( context->interpolator );
00573 
00574     interpolationOption = Dictionary_GetString_WithDefault( context->dictionary, "interpolation", "Linear" );
00575 
00576     if( !strcmp(interpolationOption, "Linear") ){
00577         context->interpolator = (_Interpolator*)LinearInterpolator_New( "linearInterpolator", context->dictionary, 
00578                                             context->localMesh, context->regularMesh );
00579     }
00580     else if ( !strcmp(interpolationOption, "Cubic") ){
00581         context->interpolator = (_Interpolator*)LinearInterpolator_New( "cubicInterpolator", context->dictionary, 
00582                                             context->localMesh, context->regularMesh );
00583     }
00584 
00585     Stg_Component_Build( context->interpolator, context->globalMesh, False );
00586     
00587     _SPModel_SimulationContextExtension_ReleaseMemory( context, LOCAL_EXTENSION );
00588     _SPModel_SimulationContextExtension_AllocateMemory( context, LOCAL_EXTENSION );
00589     
00590     SPModelSimulationContextExtension_Scatter( context );
00591 }
00592 
00593 void _SPModel_Context_Dump( SPModel_Context* self )
00594 {
00595     SurfaceMesh *mesh = NULL;
00596     int i = 0;
00597     
00598     assert( self );
00599 
00600     _SPModel_Context_Gather( self );
00601     
00602     if( self->rank != MASTER_PROC ){
00603         return;
00604     }
00605 
00606     mesh = self->globalMesh;
00607 
00608     if( self->itopography ){
00609         char fileName[500] = {0};
00610         FILE *topo = NULL;
00611         
00612         if( self->outputFormat == TEXT ){
00613             memset( fileName, 0, sizeof(fileName) );
00614             sprintf( fileName, "%s/topography%d.txt", self->outputPath, (int)self->currentTime );
00615             topo = fopen( fileName, "w+" );
00616             assert( topo );
00617 
00618             fprintf( topo, "Topography\n" );
00619             fprintf( topo, "Column Description\n" );
00620             fprintf( topo, "x\ty\th\tboundaryCondition\n" );
00621             for( i=0; i<mesh->numNodes; i++ ){
00622                 fprintf( topo, "%lf %lf %lf %lf\n", mesh->x[i],
00623                         mesh->y[i], mesh->h[i], mesh->boundaryConditions[i] );
00624             }
00625             fclose( topo );
00626         }
00627         else if( self->outputFormat == BINARY ){
00628             typedef struct outputPackage_t{
00629                 float x, y, h, bc;
00630             }outputPackage;
00631 
00632             outputPackage output;
00633             
00634             memset( fileName, 0, sizeof(fileName) );
00635             sprintf( fileName, "%s/topography%d", self->outputPath, (int)self->currentTime );
00636             topo = fopen( fileName, "w+" );
00637             assert( topo );
00638             
00639             for( i=0; i<mesh->numNodes; i++ ){
00640                 output.x = mesh->x[i];
00641                 output.y = mesh->y[i];
00642                 output.h = mesh->h[i];
00643                 output.bc = mesh->boundaryConditions[i];
00644                 fwrite( &output, sizeof( outputPackage ), 1, topo );
00645             }
00646             fclose( topo );
00647         }
00648     }
00649     
00650     if( self->idaughterNodes ){
00651         char fileName[500] = {0};
00652         FILE *topo = NULL;
00653         int j = 0, k = 0, value;
00654 
00655         if( self->outputFormat == TEXT ){
00656             memset( fileName, 0, sizeof(fileName) );
00657             sprintf( fileName, "%s/daughterNodes%d.txt", self->outputPath, (int)self->currentTime );
00658             topo = fopen( fileName, "w+" );
00659             assert( topo );
00660 
00661             fprintf( topo, "Daughter Nodes\n" );
00662             fprintf( topo, "Column Description\n" );
00663             fprintf( topo, "Immediate Daughter Node of ith node(-1 = boundary node, 0 = neighbour of itself)\n" );
00664             for( i=0; i<mesh->numNodes; i++ ){
00665                 if( !(mesh->boundaryConditions[mesh->id[i]-1]) ){
00666                     value = -1;
00667                     if( mesh->receiver[mesh->id[i]-1] != mesh->id[i] ){
00668                         value = mesh->receiver[mesh->id[i]-1];
00669                     }
00670                     fprintf( topo, "%d\n", value );
00671                 }
00672                 else{
00673                     if( mesh->receiver[mesh->id[i]-1] == mesh->id[i] ){
00674                         value = 0;
00675                         fprintf( topo, "%d\n", value );
00676                     }
00677                     else{
00678                         k = mesh->id[i];
00679                         j = mesh->receiver[mesh->id[i]-1];
00680                         value = j;
00681                         fprintf( topo, "%d\n", value );
00682                     }
00683                 }
00684             }
00685             fclose( topo );
00686         }
00687         else if( self->outputFormat == BINARY ){
00688             memset( fileName, 0, sizeof(fileName) );
00689             sprintf( fileName, "%s/daughterNodes%d", self->outputPath, (int)self->currentTime );
00690             topo = fopen( fileName, "w+" );
00691             assert( topo );
00692 
00693             for( i=0; i<mesh->numNodes; i++ ){
00694                 if( !(mesh->boundaryConditions[mesh->id[i]-1]) ){
00695                     value = -1;
00696                     if( mesh->receiver[mesh->id[i]-1] != mesh->id[i] ){
00697                         value = mesh->receiver[mesh->id[i]-1];
00698                     }
00699                     fwrite( &value, sizeof( int ), 1, topo );
00700                 }
00701                 else{
00702                     if( mesh->receiver[mesh->id[i]-1] == mesh->id[i] ){
00703                         value = 0;
00704                         fwrite( &value, sizeof( int ), 1, topo );
00705                     }
00706                     else{
00707                         k = mesh->id[i];
00708                         j = mesh->receiver[mesh->id[i]-1];
00709                         value = j;
00710                         fwrite( &value, sizeof( int ), 1, topo );
00711                     }
00712                 }
00713             }
00714             fclose( topo );
00715         }
00716     }
00717 }
00718 
00719 void _SPModel_Context_ComponentConstruct( void* context, Stg_ComponentFactory* cf )
00720 {
00721 }
00722 
00723 void _SPModel_Context_ComponentBuild( void* context, void* data )
00724 {
00725 }
00726 
00727 void _SPModel_Context_ComponentInitialise( void* context, void* data )
00728 {
00729 }
00730 
00731 void _SPModel_Context_ComponentExecute( void* context, void* data )
00732 {
00733 }
00734 
00735 void _SPModel_Context_ComponentDestroy( void* context, void* data )
00736 {
00737 }
00738