Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 7311) +++ /issm/trunk/src/c/Makefile.am (revision 7312) @@ -356,8 +356,5 @@ ./modules/GroundingLineMigrationx/GroundingLineMigrationx.h\ ./modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h\ - ./modules/GroundingLineMigrationx/CreateElementOnIceShelf.cpp\ - ./modules/GroundingLineMigrationx/CreateElementTouchingIceShelf.cpp\ - ./modules/GroundingLineMigrationx/CreateElementOnGroundingLine.cpp\ - ./modules/GroundingLineMigrationx/UpdateShelfStatus.cpp\ + ./modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp\ ./modules/ModelProcessorx/ModelProcessorx.h\ ./modules/ModelProcessorx/ModelProcessorx.cpp\ @@ -935,8 +932,5 @@ ./modules/GroundingLineMigrationx/GroundingLineMigrationx.h\ ./modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h\ - ./modules/GroundingLineMigrationx/CreateElementOnIceShelf.cpp\ - ./modules/GroundingLineMigrationx/CreateElementTouchingIceShelf.cpp\ - ./modules/GroundingLineMigrationx/CreateElementOnGroundingLine.cpp\ - ./modules/GroundingLineMigrationx/UpdateShelfStatus.cpp\ + ./modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp\ ./modules/ModelProcessorx/ModelProcessorx.h\ ./modules/ModelProcessorx/ModelProcessorx.cpp\ @@ -1180,4 +1174,5 @@ ./solutions/transient2d_core.cpp\ ./solutions/transient3d_core.cpp\ + ./solutions/groundinglinemigration2d_core.cpp\ ./solutions/steadystate_core.cpp\ ./solutions/steadystateconvergence.cpp\ Index: /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp (revision 7311) +++ /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp (revision 7312) @@ -5,45 +5,26 @@ #include "./GroundingLineMigrationx.h" #include "./GroundingLineMigrationxLocal.h" + #include "../../shared/shared.h" #include "../../include/include.h" #include "../../toolkits/toolkits.h" #include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../Container/Container.h" void GroundingLineMigrationx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ - int i,j; - Element* element=NULL; - double* element_touching_iceshelf=NULL; - int nflipped; + int migration_style; + _printf_(VerboseModule()," Migrating grounding line\n"); + + /*retrieve parameters: */ + parameters->FindParam(&migration_style,GroundingLineMigrationEnum); + + /*call different migration modules, according to user wishes: */ + if(migration_style==AgressiveMigrationEnum) AgressiveMigration(elements,nodes, vertices,loads,materials, parameters); + else if(migration_style==SoftMigrationEnum) AgressiveMigration(elements,nodes, vertices,loads,materials, parameters); + else if(migration_style==NoneEnum)_printf_("%s\n","NoneEnum supplied for migration style, doing nothing!"); + else _error_("%s not supported yet!",EnumToString(migration_style)); - /*Loop until no more nodes and elements change from grounded to floating: */ - nflipped=1; //arbitrary, get things started - - while(nflipped){ - - /*reset counter: */ - nflipped=0; - - /*Create vector holding all the elements that touch the ice shelf, by any node: */ - element_touching_iceshelf=CreateElementTouchingIceShelf(elements); - - /*Carry out grounding line migration for those elements: */ - for(i=0;iSize();i++){ - element=(Element*)elements->GetObjectByOffset(i); - if(element_touching_iceshelf[element->Sid()]) element->MigrateGroundingLine(); - } - - /*Now, update shelf flags in nodes and elements: */ - nflipped=UpdateShelfStatus(elements,nodes,parameters,element_touching_iceshelf); - _printf_(VerboseModule()," number of migrated nodes: %i\n",nflipped); - - /*avoid memory leaks: */ - xfree((void**)&element_touching_iceshelf); - } - - - /*free ressouces: */ - xfree((void**)&element_touching_iceshelf); } Index: /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h =================================================================== --- /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h (revision 7311) +++ /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h (revision 7312) @@ -11,8 +11,12 @@ /* local prototypes: */ + +void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); +void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); + double* CreateElementOnIceShelf(Elements* elements); bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf); double* CreateElementTouchingIceShelf(Elements* elements); -int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf); +int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf); #endif /* _GROUNDINGLINEMIGRATIONXLOCAL_H */ Index: /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp (revision 7312) +++ /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp (revision 7312) @@ -0,0 +1,223 @@ +/*!\file GroundingLineMigrationx + * \brief: migration grounding line position. + */ + +#include "./GroundingLineMigrationx.h" +#include "./GroundingLineMigrationxLocal.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + + /*Here, whatever grids inside the ice sheet want to unground, we allow -> instantaneous transmission of water through the bedrock: */ + + int i,j; + Element* element=NULL; + + _printf_(VerboseModule()," Migrating grounding line\n"); + + /*Carry out grounding line migration for those elements: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->AgressiveMigration(); + } + + /*Synchronize shelf status: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->AgressiveShelfSync(); + } + +} + +bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf){ + + int i; + int j; + Element *element = NULL; + bool *element_on_gl = NULL; + bool ongl=false; + int *neighboorsids = NULL; + int sid; + int shelf; + + /*Go through elements, and look for elements that can possibly have a grounding line migration. These + * are elements that touch the grounding line: */ + element_on_gl=(bool*)xmalloc(elements->Size()*sizeof(bool)); + + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + + sid=element->Sid(); + neighboorsids=element->GetHorizontalNeighboorSids(); + + shelf=(int)element_on_iceshelf[sid]; + ongl=false; + for(j=0;j<3;j++){ + if (neighboorsids[j]<0)continue; //this neighboor does not exist + if ((shelf==1) & (element_on_iceshelf[neighboorsids[j]]==1))continue; //both this element and this neighboor are on th ice shelf + if ((shelf==0) & (element_on_iceshelf[neighboorsids[j]]==0))continue; //both this element and this neighboor are on the ice sheet + ongl=true; //neighboor j is on a different location than us, ie we are touching the grounding line. + } + element_on_gl[i]=ongl; + } + + return element_on_gl; +} + +double* CreateElementOnIceShelf(Elements* elements){ + + int i; + Element* element=NULL; + Vec vec_element_on_iceshelf=NULL; + double* element_on_iceshelf=NULL; + + /*Create vector holding all the elements IsOnShelf flags: */ + vec_element_on_iceshelf=NewVec(elements->NumberOfElements(),true); + + /*Loop through elements, and fill vec_element_on_iceshelf: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + printf("%i\n",i); + VecSetValue(vec_element_on_iceshelf,element->Sid(),(int)element->IsOnShelf(),INSERT_VALUES); + } + + /*Assemble vector: */ + VecAssemblyBegin(vec_element_on_iceshelf); + VecAssemblyEnd(vec_element_on_iceshelf); + + /*Serialize vector: */ + VecToMPISerial(&element_on_iceshelf,vec_element_on_iceshelf); + + /*free ressouces: */ + VecFree(&vec_element_on_iceshelf); + + return element_on_iceshelf; +} + + +double* CreateElementTouchingIceShelf(Elements* elements){ + + int i; + Element* element=NULL; + Vec vec_element_touching_iceshelf=NULL; + double* element_touching_iceshelf=NULL; + + /*Create vector holding all the elements IsOnShelf flags: */ + vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true); + + /*Loop through elements, and fill vec_element_touching_iceshelf: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelf()?1.0:0.0,INSERT_VALUES); + } + + /*Assemble vector: */ + VecAssemblyBegin(vec_element_touching_iceshelf); + VecAssemblyEnd(vec_element_touching_iceshelf); + + /*Serialize vector: */ + VecToMPISerial(&element_touching_iceshelf,vec_element_touching_iceshelf); + + /*free ressouces: */ + VecFree(&vec_element_touching_iceshelf); + + return element_touching_iceshelf; +} + + +int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){ + + int i; + Element* element=NULL; + int analysis_type; + int numnods; + Vec vec_new_shelf_nodes=NULL; + double* new_shelf_nodes=NULL; + + /*output: */ + int local_nflipped=0; + int nflipped=0; + + /*recover parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */ + numnods=nodes->NumberOfNodes(analysis_type); + vec_new_shelf_nodes=NewVec(numnods); + + /*Ok, now go through the elements that have gone through grounding line migration, + * and update their flags: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + if(element_touching_iceshelf[element->Sid()]){ + local_nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes); + } + } + MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); + + /*Serialize vec_new_shelf_nodes: */ + VecToMPISerial(&new_shelf_nodes,vec_new_shelf_nodes); + + /*Now, go through ALL elements, and update the status of the nodes, to propagate what happened at the grounding line: */ + /*Carry out grounding line migration for those elements: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->UpdateShelfFlags(new_shelf_nodes); + } + + /*Free ressources: */ + VecFree(&vec_new_shelf_nodes); + xfree((void**)&new_shelf_nodes); + + return nflipped; +} + + +void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + + int i,j; + Element* element=NULL; + double* element_touching_iceshelf=NULL; + int nflipped; + + _printf_(VerboseModule()," Migrating grounding line\n"); + + /*Loop until no more nodes and elements change from grounded to floating: */ + nflipped=1; //arbitrary, get things started + + int count=0; + while(nflipped){ + + if (count==4)break; + + /*reset counter: */ + nflipped=0; + + /*Create vector holding all the elements that touch the ice shelf, by any node: */ + element_touching_iceshelf=CreateElementTouchingIceShelf(elements); + + + /*Carry out grounding line migration for those elements: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + if(element_touching_iceshelf[element->Sid()]) element->MigrateGroundingLine(); + } + + /*Now, update shelf flags in nodes and elements: */ + nflipped=UpdateShelfStatus(elements,nodes,parameters,element_touching_iceshelf); + //_printf_(VerboseModule()," number of migrated nodes: %i\n",nflipped); + extern int my_rank;if(my_rank==0)printf(" number of migrated nodes: %i\n",nflipped); + + /*avoid memory leaks: */ + xfree((void**)&element_touching_iceshelf); + count++; + } + + + /*free ressouces: */ + xfree((void**)&element_touching_iceshelf); + +}