SPModel: plugins/diffusionErosion/DiffusionErosion.c Source File
VPAC - Computational Software Development
Main | SPModel | StGermain FrameWork |
Main Page | Alphabetical List | Class List | Directories | File List | Class Members | File Members

DiffusionErosion.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: Erosion.c 225 2005-12-22 00:01:19Z AlanLo $
00021 **
00022 **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
00023 
00024 
00025 #include <mpi.h>
00026 #include <StGermain/StGermain.h>
00027 #include "SPModel/SPModel.h"
00028 #include <stdlib.h>
00029 #include <string.h>
00030 #include <assert.h>
00031 #include "DiffusionErosion.h"
00032 
00033 const Name SPModelFlowHierarchyViewerPluginName = "SPModelFlowHierarchyViewer";
00034 
00035 ExtensionInfo_Index SPModel_Simulation_ContextExtHandle;
00036 
00037 /* Textual name of this module */
00038 const Type SPModelDiffusionErosion_Type = "SPModelDiffusionErosion";
00039 
00040 void* _SPModelDiffusionErosion_DefaultNew( Name name ) {
00041     return Codelet_New(
00042             SPModelDiffusionErosion_Type,
00043             _SPModelDiffusionErosion_DefaultNew,
00044             _SPModelDiffusionErosion_Construct,
00045             _Codelet_Build,
00046             _Codelet_Initialise,
00047             _Codelet_Execute,
00048             _Codelet_Destroy,
00049             name );
00050 }
00051 
00052 void _SPModelDiffusionErosion_Construct( void* component, Stg_ComponentFactory* data ) {    
00053     SPModel_Context*                context;
00054     
00055     context = (SPModel_Context*)Stg_ComponentFactory_ConstructByName( data, "context", SPModel_Context, True );
00056     
00057     _SPModelDiffusionErosion_Init( context );
00058     
00059     Journal_DPrintf( context->debug, "In %s():\n", __func__ );
00060     ContextEP_Append( context, AbstractContext_EP_Solve, SPModelDiffusionErosion_Solve );
00061 }
00062 
00063 void _SPModelDiffusionErosion_Init( void* _context ) {
00064     SPModel_Context *context = (SPModel_Context*)_context;
00065     
00066     assert( context );
00067 
00068     SPModel_Simulation_ContextExtHandle = ExtensionManager_GetHandle( context->extensionMgr, "Simulation" );
00069 }
00070 
00071 void SPModelDiffusionErosion_Solve( void *_context )
00072 {
00073     int i, j;
00074     int localIndex = 0;
00075     int neighbour, neighbourLocalIndex;
00076     float dx, dy, dr, dh, dhdr;
00077     float minSediment;
00078     SPModel_Context *context;
00079     SurfaceMesh* mesh = NULL;
00080     float diffusivity;
00081     double material = 0.0f;
00082     double penalty = 0.0f;
00083     double factor = 0.0f;
00084     SPModelSimulationContextExtension *simulationExt= NULL;
00085     
00086     context = (SPModel_Context*)_context;
00087     assert( context );
00088 
00089     simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
00090                 ( context->extensionMgr, context, SPModel_Simulation_ContextExtHandle );
00091     assert( simulationExt );
00092 
00093     mesh = context->localMesh;
00094 
00095     minSediment = 0.1f;
00096     for( i=0; i<mesh->myLoad; i++ ){
00097         localIndex = mesh->mapGlobalToLocal[mesh->sortedId[i]-1];
00098         assert( localIndex != mesh->myLoad );
00099 
00100         
00101         for (j=0; j<mesh->numNeigh[mesh->sortedId[i]-1]; j++) {
00102             neighbour = mesh->nodeNeighbours[localIndex][j];
00103             
00104             factor = 1.0f;
00105             if( ((mesh->boundaryConditions[mesh->sortedId[i]-1])==2) && (!(mesh->boundaryConditions[neighbour-1])) ){
00106                 factor = 0.5f;
00107             }
00108             else if( (!(mesh->boundaryConditions[mesh->sortedId[i]-1])) && ((mesh->boundaryConditions[neighbour-1])==2) ){
00109                 factor = 2.0f;
00110             }
00111 
00112             neighbourLocalIndex = mesh->mapGlobalToLocal[neighbour-1];
00113     
00114             if( simulationExt->sedimentHistory[localIndex] > minSediment ){
00115                 if( (mesh->h[localIndex] < context->seaLevel->value) ||  (mesh->h[neighbourLocalIndex] < context->seaLevel->value) ){
00116                     diffusivity = context->submarineSedimentDiffusivity;
00117                 }
00118                 else{
00119                     diffusivity = context->sedimentDiffusivity;
00120                 }
00121             }
00122             else{
00123                 if( (mesh->h[localIndex] < context->seaLevel->value) ||  (mesh->h[neighbourLocalIndex] < context->seaLevel->value) ){
00124                     diffusivity = context->submarineBedRockDiffusivity;
00125                 }
00126                 else{
00127                     diffusivity = context->bedRockDiffusivity;
00128                 }
00129             }
00130 
00131         
00132             if( neighbourLocalIndex != mesh->myLoad ){
00133                 if( mesh->h[neighbourLocalIndex] < mesh->h[localIndex] ){
00134                     dx = mesh->x[localIndex] - mesh->x[neighbourLocalIndex];
00135                     dy = mesh->y[localIndex] - mesh->y[neighbourLocalIndex];
00136                     dh = mesh->h[localIndex] - mesh->h[neighbourLocalIndex];
00137                 
00138                     dr = sqrt( dx*dx + dy*dy );
00139                     dhdr = dh/dr;
00140                     
00141                     material = (2*sqrt(PI*mesh->surface[localIndex]) * diffusivity * dhdr * context->dt) / 
00142                                 ((float)(mesh->numNeigh[localIndex]));
00143                     
00144                     penalty = (mesh->surface[localIndex] * mesh->surface[neighbourLocalIndex]) / 
00145                                 (mesh->surface[localIndex] + mesh->surface[neighbourLocalIndex]) *
00146                                 (mesh->h[localIndex] - mesh->h[neighbourLocalIndex]);
00147                     
00148                     if( material > penalty ){
00149                         material = penalty;
00150                         printf("rank %d, localIndex %d, material %8.20f, penalty %8.20f\n", mesh->rank, localIndex, material, penalty);
00151                         fprintf( stderr, "incurred penalty in diffusionErosion routine\n" );
00152                     }
00153 
00154                     mesh->h[localIndex] -= material / mesh->surface[localIndex]*factor;
00155                     mesh->h[neighbourLocalIndex] += material / mesh->surface[neighbourLocalIndex]*factor;
00156                     
00157                     simulationExt->sedimentHistory[localIndex] -= material / mesh->surface[localIndex]*factor;
00158                     if( simulationExt->sedimentHistory[localIndex] < 0.0f ){
00159                         simulationExt->erosion[localIndex] -= simulationExt->sedimentHistory[localIndex]*factor;
00160                         simulationExt->sedimentHistory[localIndex] = 0.0f;
00161                     }
00162                     simulationExt->sedimentHistory[neighbourLocalIndex] += material / mesh->surface[neighbourLocalIndex]*factor;
00163                 }
00164             }
00165             else{
00166 
00167                 /*printf( "foreign nodes not handled at the moment\n" );*/
00168             }
00169         }
00170     }
00171 }
00172 
00173 Index SPModelDiffusionErosion_Register( PluginsManager* pluginsManager )
00174 {
00175     return PluginsManager_Submit( pluginsManager, SPModelDiffusionErosion_Type, "0", _SPModelDiffusionErosion_DefaultNew );
00176 }
00177