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

SurfaceProcesses.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: SurfaceProcesses.c 262 2006-02-16 02:51:37Z RaquibulHassan $
00021 **
00022 **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
00023 
00024 #include <mpi.h>
00025 #include <StGermain/StGermain.h>
00026 #include <Cascade/cascade.h>
00027 
00028 
00029 #include "types.h"
00030 #include "SurfaceMesh.h"
00031 
00032 #include <stdio.h>
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include <assert.h>
00036 #include <math.h>
00037 #include "SPModel.h"
00038 
00039 void fluvial( SPModel_Context *context, SPModelSimulationContextExtension *simulationExt, TransportEquation *transportEq )
00040 {
00041     typedef struct package_t{
00042         int id;
00043         float water;
00044         float sediment;
00045     }package;
00046     
00047     double Qe, lFac, temp, minSediment, eFrac, sFrac;
00048     SurfaceMesh *mesh = NULL;
00049     int i;
00050     int receiver;
00051     float dr, highestNeighbourHeight;
00052     float receiverHeight, slope;
00053     int haloNodeProc, haloNodeIndex;
00054     int tempHaloNodeProc, tempHaloNodeIndex;
00055     int providerHaloNodeProc, providerHaloNodeIndex;
00056     int foreignReceiver = 0;
00057     int localIndex = 0;
00058     int receiverLocalIndex = 0, highestNeighbourLocalIndex = 0;
00059     package waterSed;
00060     LinkedListIterator *iterator = NULL;
00061     int *data = NULL;
00062     double factor = 1.0f;
00063  
00064     minSediment = 1.0; /* maximum detectable sediment set to 1m */
00065 
00066     eFrac = 0.5; /* maximum erosional fraction */
00067 
00068     sFrac = 1.0; /* maximum sedimentational fraction */
00069     
00070     assert( context );
00071 
00072     assert( transportEq );
00073     
00074     assert( simulationExt );
00075 
00076     mesh = context->localMesh;
00077     
00078     iterator = LinkedListIterator_New( mesh->nodeProviders[0] );
00079     
00080     context->itemsSent = 0;
00081     context->itemsReceived = 0;
00082     
00083     for( i=0; i<mesh->myLoad; i++ )
00084     {
00085         localIndex = mesh->mapGlobalToLocal[mesh->sortedId[i] - 1];
00086         assert( localIndex != mesh->myLoad );
00087         
00088         /* Only process nodes that are above the sea level */
00089         if( mesh->h[localIndex] > context->seaLevel->value ){
00090             foreignReceiver = FALSE;
00091 
00092             receiver = mesh->receiver[localIndex];
00093             receiverLocalIndex = mesh->mapGlobalToLocal[receiver-1];
00094         
00095             factor = 1.0f;
00096             if( ((mesh->boundaryConditions[mesh->sortedId[i]-1])==2) && (!(mesh->boundaryConditions[receiver-1])) ){
00097                 factor = 0.5f;
00098             }
00099             else if( (!(mesh->boundaryConditions[mesh->sortedId[i]-1])) && ((mesh->boundaryConditions[receiver-1])==2) ){
00100                 factor = 2.0f;
00101             }
00102             
00103             /* Checking if the receiver is on a different processor */
00104             if( receiverLocalIndex == mesh->myLoad ){
00105                 foreignReceiver = TRUE;
00106 
00107                 SurfaceMesh_FindHaloNode( mesh, receiver, &haloNodeProc, &haloNodeIndex );
00108             }
00109             
00110             /* Finding the receivers height */
00111             if( receiverLocalIndex == mesh->myLoad ){
00112                 receiverHeight = mesh->haloNodes[haloNodeProc].h[haloNodeIndex];
00113             }
00114             else{
00115                 receiverHeight = mesh->h[receiverLocalIndex];
00116             }
00117             
00118             if( receiverHeight < mesh->h[localIndex] ){
00119                 iterator->list = mesh->nodeProviders[localIndex];
00120                 for( data=(int*)LinkedListIterator_First( iterator ); data; 
00121                         data = (int*)LinkedListIterator_Next( iterator ) ){
00122                 
00123                     if( mesh->mapGlobalToLocal[*data-1] == mesh->myLoad ){
00124                         SurfaceMesh_FindHaloNode( mesh, *data, &providerHaloNodeProc, &providerHaloNodeIndex );
00125 
00126                         MPI_Recv( &waterSed, sizeof( package ), MPI_BYTE, providerHaloNodeProc,
00127                                *data, MPI_COMM_WORLD, &(mesh->statusTable[0]) );
00128                         context->itemsReceived++;
00129                         assert( waterSed.id == *data );
00130                         simulationExt->water[localIndex] += waterSed.water;
00131                         simulationExt->sediment[localIndex] += waterSed.sediment;
00132                     
00133                         #ifdef USE_DEBUG
00134                         printf( "receiving water %lf and sediment %lf, on proc %d from proc %d from node %d\n",
00135                                 waterSed.water, waterSed.sediment, mesh->rank, haloNodeProc, *data );
00136                         #endif
00137                     }
00138                 }
00139 
00140                 slope = -mesh->slope[localIndex];
00141 
00142                 /* Transport equation */
00143                 Qe = transportEq( localIndex, mesh, simulationExt );
00144                 
00145                 
00146                 /* Finding the highest neighbours height */
00147                 highestNeighbourLocalIndex = mesh->mapGlobalToLocal[mesh->highestNeighbour[localIndex]-1];
00148                 if( highestNeighbourLocalIndex == mesh->myLoad ){
00149                     SurfaceMesh_FindHaloNode( mesh, mesh->highestNeighbour[localIndex], &tempHaloNodeProc, &tempHaloNodeIndex );
00150                     highestNeighbourHeight = mesh->haloNodes[tempHaloNodeProc].h[tempHaloNodeIndex];
00151                 }
00152                 else{
00153                     highestNeighbourHeight = mesh->h[highestNeighbourLocalIndex];
00154                 }
00155                 
00156                 
00157                 /* sedimentation occuring */
00158                 if( simulationExt->sediment[localIndex] > Qe ){
00159                     temp = ( simulationExt->sediment[localIndex] - Qe ) / mesh->surface[localIndex];
00160                         
00161                     /* Penalty check */
00162                     if( temp > sFrac * ( highestNeighbourHeight - receiverHeight ) ){
00163                         temp = sFrac * ( highestNeighbourHeight - receiverHeight );
00164                         fprintf( stderr, "penalty incurred in  transport routine while depositing\n" );
00165                     }
00166                     
00167                     mesh->h[localIndex] += temp*factor;
00168                     simulationExt->sedimentHistory[localIndex] += temp*factor;
00169                     simulationExt->sediment[localIndex] -= temp * mesh->surface[localIndex]*factor;
00170                 }
00171                 /* erosion occuring */
00172                 else{
00173                     dr = fabs( slope / (mesh->h[localIndex] - receiverHeight) );
00174                     if( simulationExt->sedimentHistory[localIndex] > minSediment ){ /* sediments are present */
00175                         lFac = dr / context->sedimentErosionLengthScale;
00176                     }
00177                     else{ /* no deposited sediments are present */
00178                         lFac = dr / context->bedRockErosionLengthScale;
00179                     }
00180 
00181                     /* Penalty check on length factor */
00182                     if( lFac > 1.0f ){
00183                         lFac = 1.0f;
00184                     }
00185 
00186                     temp = ( (simulationExt->sediment[localIndex] - Qe)/mesh->surface[localIndex] * lFac );
00187                     
00188                     /* penalty check */
00189                     if( temp > eFrac*( mesh->h[localIndex] - receiverHeight ) ){
00190                         temp = eFrac*( mesh->h[localIndex] - receiverHeight );
00191                         fprintf( stderr, "penalty incurred in  transport routine while eroding\n" );
00192                     }
00193                     
00194                     mesh->h[localIndex] += temp*factor;
00195                     simulationExt->sedimentHistory[localIndex] += temp*factor;
00196                     if( simulationExt->sedimentHistory[localIndex] < 0.0f ){
00197                         simulationExt->erosion[localIndex] -= simulationExt->sedimentHistory[localIndex]*factor;
00198                         simulationExt->sedimentHistory[localIndex] = 0.0f;
00199                     }
00200                     simulationExt->sediment[localIndex] -= temp * mesh->surface[localIndex]*factor;
00201                 }
00202         
00203                 if( foreignReceiver == FALSE ){
00204                     simulationExt->water[receiverLocalIndex] += simulationExt->water[localIndex] * mesh->surface[localIndex]
00205                                                 / mesh->surface[receiverLocalIndex]
00206                                                 * factor;
00207                     simulationExt->sediment[receiverLocalIndex] += simulationExt->sediment[localIndex] * factor;
00208                 }
00209                 else{
00210                     waterSed.id = mesh->id[localIndex];
00211                     waterSed.water = simulationExt->water[localIndex] * mesh->surface[localIndex] 
00212                                         / mesh->haloNodes[haloNodeProc].surface[haloNodeIndex]
00213                                         * factor;
00214                     waterSed.sediment = simulationExt->sediment[localIndex] * factor;
00215 
00216                     MPI_Isend( &waterSed, sizeof( package ), MPI_BYTE, haloNodeProc,
00217                             mesh->id[localIndex], MPI_COMM_WORLD, &(mesh->requestTable[localIndex]) );
00218                     context->itemsSent++;
00219                     #ifdef USE_DEBUG
00220                     printf( "sending water %lf and sediment %lf, from proc %d to proc %d from node %d\n",
00221                         waterSed.water, waterSed.sediment, mesh->rank, haloNodeProc, mesh->id[localIndex] );
00222                     #endif
00223                 }
00224                 simulationExt->sediment[localIndex] = 0.0f;
00225                 simulationExt->water[localIndex] = 0.0f;
00226             }
00227             else{
00228                 if( foreignReceiver ){
00229                     memset( &waterSed, 0, sizeof( package ) );
00230                     
00231                     waterSed.id = mesh->id[localIndex];
00232                     
00233                     MPI_Isend( &waterSed, sizeof( package ), MPI_BYTE, haloNodeProc,
00234                         mesh->id[localIndex], MPI_COMM_WORLD, &(mesh->requestTable[localIndex]) );
00235                     context->itemsSent++;
00236                     #ifdef USE_DEBUG
00237                     printf( "sending dummy values to node %d on proc %d\n", mesh->id[localIndex], haloNodeProc );
00238                     #endif
00239                 }
00240                 
00241                 temp = simulationExt->sediment[localIndex]/mesh->surface[localIndex];
00242                 mesh->h[localIndex] += temp*factor;
00243                 simulationExt->sedimentHistory[localIndex] += temp*factor;
00244                 simulationExt->sediment[localIndex] = 0.0f;
00245                 simulationExt->water[localIndex] = 0.0f; /* water is lost to the subsurface */
00246             }
00247         }
00248         else{
00249             temp = simulationExt->sediment[localIndex] / mesh->surface[localIndex];
00250             mesh->h[localIndex] += temp*factor;
00251             simulationExt->sedimentHistory[localIndex] += temp*factor;
00252             simulationExt->sediment[localIndex] = 0.0f;
00253         }
00254     }
00255     #ifdef USE_DEBUG
00256     printf( "sent %d, received %d, leaving  erosion\n", context->itemsSent, context->itemsReceived );
00257     #endif
00258     Stg_Class_Delete( iterator );
00259 }