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

testInterpolator.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 */
00036 #include <mpi.h>
00037 #include <StGermain/StGermain.h>
00038 #ifdef SDL_INSTALLED
00039     #include <SDL/SDL.h>
00040 #endif
00041 #include "SPModel/SPModel.h"
00042 
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <string.h>
00046 #include <assert.h>
00047 
00048 float *erosion;
00049 float *sedimentHistory;
00050 
00051 float interpolateFunction ( int meshLocalIndex, void *args )
00052 {
00053     SPModel_Context *context = (SPModel_Context*)args;
00054     float result = 0.0f;
00055 
00056     assert( context );
00057 
00058     result = 0.0f;
00059     
00060     if( erosion ){
00061         result = -erosion[meshLocalIndex] * context->bedRockDensity*context->gravity;
00062     }
00063     if( sedimentHistory ){
00064         result += sedimentHistory[meshLocalIndex] * context->sedimentDensity*context->gravity;
00065     }
00066     
00067     return result;
00068 }
00069 
00070 int main (int argc, char **argv)
00071 {
00072     MPI_Comm            CommWorld;
00073     int             rank;
00074     int             numProcessors;
00075     int             procToWatch;
00076     Dictionary*         dictionary;
00077     XML_IO_Handler*         ioHandler;
00078     char*               filename;
00079     SPModel_Context*        self;
00080     LinearInterpolator      *interpolator;
00081     int i = 0, j = 0;
00082     Stream *stream;
00083     double expectedValue = 0.0f, sqSum = 0.0f, deviation = 0.0f;
00084     const double epsilon = 1e-2;
00085     
00086     /* Initialise MPI, get world info */
00087     MPI_Init( &argc, &argv );
00088     MPI_Comm_dup( MPI_COMM_WORLD, &CommWorld );
00089     MPI_Comm_size( CommWorld, &numProcessors );
00090     MPI_Comm_rank( CommWorld, &rank );
00091     if( argc >= 3 ) {
00092         procToWatch = atoi( argv[1] );
00093     }
00094     else {
00095         procToWatch = 0;
00096     }
00097     if( rank == procToWatch ) printf( "Watching rank: %i\n", rank );
00098     
00099     
00100     if (!StGermain_Init( &argc, &argv )) {
00101         fprintf(stderr, "Error initialising StGermain, exiting.\n" );
00102         exit(EXIT_FAILURE);
00103     }
00104     if (!SPModel_Init( &argc, &argv )) {
00105         fprintf(stderr, "Error initialising SPModel, exiting.\n" );
00106         exit(EXIT_FAILURE);
00107     }
00108     
00109     /* Create the dictionary, and some fixed values */
00110     dictionary = Dictionary_New();
00111     Dictionary_Add( dictionary, "rank", Dictionary_Entry_Value_FromUnsignedInt( rank ) );
00112     Dictionary_Add( dictionary, "numProcessors", Dictionary_Entry_Value_FromUnsignedInt( numProcessors ) );
00113     
00114     /* Read input */
00115     ioHandler = XML_IO_Handler_New();
00116     if( argc >= 3 ) {
00117         filename = strdup( argv[2] );
00118     }
00119     else {
00120         filename = strdup( "input.xml" );
00121     }
00122     if ( False == IO_Handler_ReadAllFromFile( ioHandler, filename, dictionary ) )
00123     {
00124         fprintf( stderr, "Error: SPModel couldn't find specified input file %s. Exiting.\n", filename );
00125         exit( EXIT_FAILURE );
00126     }
00127     Journal_ReadFromDictionary( dictionary );
00128     
00129     self = SPModel_Context_New( "SPModel_Context", 0.0f, 0.0f, CommWorld, dictionary );
00130     
00131     if( self->rank == MASTER_PROC ){
00132         
00133         Stg_Component_Build( self->globalMesh, NULL, False );
00134         Stg_Component_Initialise( self->globalMesh, NULL, False );
00135         Stg_Component_Build( self->catchmentList, NULL, False );
00136         Stg_Component_Build( self->meshDecomp, NULL, False );       
00137     }
00138 
00139     Stg_Component_Execute( (void*)self->meshDecomp, (void*)self->localMesh, False );
00140     Stg_Component_Initialise( self->localMesh, NULL, False );
00141     
00142     Stg_Component_Build( self->regularMesh, NULL, False );
00143     Stg_Component_Build( self->interpolator, self->globalMesh, False );
00144 
00145     stream = Journal_Register( Info_Type, "stream" );
00146 
00147     if( self->rank == MASTER_PROC ){
00148         Journal_Printf( stream, "Printing the Interpolator on the master process\n" );
00149         Stg_Class_Print( self->interpolator, stream );
00150     }
00151 
00152     interpolator = (LinearInterpolator*)self->interpolator;
00153     for( i=0; i<self->localMesh->myLoad; i++ ){
00154         for( j=0; j<4; j++ ){
00155             if( interpolator->meshWeights[self->localMesh->id[i]-1][j] > 1.0f ){
00156                 Journal_Printf( stream, "Invalid meshWeight in Interpolator of processor %d\n", self->rank );
00157                 assert( 0 );
00158             }
00159         }
00160     }
00161 
00162     erosion = malloc( sizeof( float ) * self->localMesh->myLoad );
00163     sedimentHistory = malloc( sizeof( float ) * self->localMesh->myLoad );
00164 
00165     for( i=0; i<self->localMesh->myLoad; i++ ){
00166         erosion[i] = 0.1f;
00167         sedimentHistory[i] = 0.01f;
00168     }
00169 
00170     _Interpolator_InterpolateFromMeshToGrid( (_Interpolator*)self->interpolator, interpolateFunction, self );
00171     
00172     SurfaceRegularMesh_GatherData( self->regularMesh, self->interpolator );
00173 
00174     /*if( self->rank == MASTER_PROC ){
00175         for( i=0; i<self->regularMesh->nx; i++ ){
00176             for( j=0; j<self->regularMesh->ny; j++ ){
00177                 printf( "flexureArray[%d][%d] = %lf-->\n", i, j, self->regularMesh->flexureArray[i][j] );
00178             }
00179         }
00180     }*/
00181 
00182     expectedValue = -(0.1f * self->bedRockDensity * self->gravity) + ( 0.01 * self->sedimentDensity * self->gravity );
00183     
00184     sqSum = 0.0f;
00185     if( self->rank == MASTER_PROC ){
00186         Journal_Printf( stream, "\n\nChecking Interpolation weights.. Done\n" );
00187         Journal_Printf( stream, "Checking parallel interpolation.. Done\n" );
00188 
00189         for( i=0; i<self->regularMesh->nx; i++ ){
00190             for( j=0; j<self->regularMesh->nx; j++ ){
00191                 sqSum += pow( expectedValue - self->regularMesh->dataArray[i][j], 2.0f );
00192             }
00193         }
00194         
00195         deviation = sqrt(sqSum / (double)(self->regularMesh->nx * self->regularMesh->ny));
00196 
00197         if( deviation > epsilon ){
00198             Journal_Printf( stream, "Standard Deviation too high.. aborting..!\n" );
00199             assert( 0 );
00200         }
00201         Journal_Printf( stream, "Verifying parallel assembly.. Done\n" );
00202     }
00203 
00204     SurfaceRegularMesh_BroadcastData( self->regularMesh );
00205     if( self->rank == MASTER_PROC ){
00206     }
00207     else{
00208         for( i=0; i<self->regularMesh->nx; i++ ){
00209             for( j=0; j<self->regularMesh->nx; j++ ){
00210                 self->regularMesh->dataArray[i][j] = 0.0f;
00211             }
00212         }
00213     
00214         /*for( i=0; i<self->regularMesh->nx; i++ ){
00215             for( j=0; j<self->regularMesh->nx; j++ ){
00216                 printf( "-->flexureArray[%d][%d] = %lf\n", i, j, self->regularMesh->flexureArray[i][j] );
00217             }
00218         }*/
00219     }
00220 
00221     StGermain_Finalise();
00222     MPI_Finalize();
00223     
00224     return 0; /* success */
00225 }
00226