Changeset 7319


Ignore:
Timestamp:
02/03/11 16:23:38 (14 years ago)
Author:
Eric.Larour
Message:

Grounding line migration routines

Location:
issm/trunk/src/c/objects/Elements
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Element.h

    r7288 r7319  
    8383                virtual int*   GetHorizontalNeighboorSids(void)=0;
    8484                virtual double TimeAdapt()=0;
     85                virtual void   AgressiveMigration()=0;
     86                virtual void   AgressiveShelfSync()=0;
    8587                virtual void   MigrateGroundingLine()=0;
    8688                virtual int    UpdateShelfStatus(Vec new_shelf_nodes)=0;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r7288 r7319  
    302302void  Penta::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){
    303303        _error_("Not supported yet!");
     304}
     305/*}}}*/
     306/*FUNCTION Penta::AgressiveMigration{{{1*/
     307void  Penta::AgressiveMigration(void){
     308        _error_("not supported yet!");
     309}
     310/*}}}*/
     311/*FUNCTION Penta::AgressiveShelfSync{{{1*/
     312void  Penta::AgressiveShelfSync(void){
     313        _error_("not supported yet!");
    304314}
    305315/*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r7288 r7319  
    109109                void   MaxVy(double* pmaxvy, bool process_units);
    110110                void   MaxVz(double* pmaxvz, bool process_units);
     111                void   AgressiveMigration();
     112                void   AgressiveShelfSync();
    111113                void   MigrateGroundingLine();
    112114                void   MinVel(double* pminvel, bool process_units);

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r7288 r7319  
    276276
    277277/*Other*/
     278/*FUNCTION Tria::AgressiveMigration{{{1*/
     279void  Tria::AgressiveMigration(void){
     280
     281
     282        double *values         = NULL;
     283        double  h[3],s[3],b[3],ba[3];
     284        double  bed_hydro;
     285        double  rho_water,rho_ice,density;
     286        int     i;
     287        bool    elementonshelf = false;
     288
     289        /*Recover info at the vertices: */
     290        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     291        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     292        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     293
     294        GetParameterListOnVertices(&h[0],ThicknessEnum);
     295        GetParameterListOnVertices(&s[0],SurfaceEnum);
     296        GetParameterListOnVertices(&b[0],BedEnum);
     297        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     298
     299        /*material parameters: */
     300        rho_water=matpar->GetRhoWater();
     301        rho_ice=matpar->GetRhoIce();
     302        density=rho_ice/rho_water;
     303
     304        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     305        for(i=0;i<3;i++){
     306                if (nodes[i]->IsOnShelf()){
     307                        /*This node is on the shelf. See if its bed is going under the bathymetry: */
     308                        if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
     309                                /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
     310                                b[i]=ba[i];
     311                                s[i]=b[i]+h[i];
     312                        }
     313                        else{
     314                                /*do nothing, we are still floating.*/
     315                        }
     316                }
     317                else{
     318                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     319                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     320                        bed_hydro=-density*h[i];
     321                        if (bed_hydro>ba[i]){
     322                                /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
     323                                s[i]=(1-density)*h[i];
     324                                b[i]=-density*h[i];
     325                        }
     326                        else{
     327                                /*do nothing, we are still grounded.*/
     328                        }
     329                }
     330        }
     331
     332        /*Surface and bed are updated. Update inputs:*/
     333        surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
     334        bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
     335
     336        }
     337/*}}}*/
     338/*FUNCTION Tria::AgressiveShelfSync{{{1*/
     339void  Tria::AgressiveShelfSync(void){
     340
     341
     342        double *values         = NULL;
     343        double  h[3],s[3],b[3],ba[3];
     344        double  bed_hydro;
     345        double  rho_water,rho_ice,density;
     346        int     i;
     347        bool    elementonshelf = false;
     348
     349        /*Recover info at the vertices: */
     350        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     351        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     352        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     353
     354        GetParameterListOnVertices(&h[0],ThicknessEnum);
     355        GetParameterListOnVertices(&s[0],SurfaceEnum);
     356        GetParameterListOnVertices(&b[0],BedEnum);
     357        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     358
     359        /*material parameters: */
     360        rho_water=matpar->GetRhoWater();
     361        rho_ice=matpar->GetRhoIce();
     362        density=rho_ice/rho_water;
     363
     364        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     365        for(i=0;i<3;i++){
     366                if(b[i]==ba[i]){
     367                               
     368                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false));
     369                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true));
     370                }
     371                else{
     372                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true));
     373                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false));
     374
     375                }
     376        }
     377
     378        /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
     379        elementonshelf=false;
     380        for(i=0;i<3;i++){
     381                if(nodes[i]->IsOnShelf()){
     382                        elementonshelf=true;
     383                        break;
     384                }
     385        }
     386    this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf));
     387}
     388/*}}}*/
    278389/*FUNCTION Tria::AverageOntoPartition {{{1*/
    279390void  Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){
    … …  
    35993710                this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs));
    36003711        }
    3601         if(iomodel->gl_migration){
     3712        if(iomodel->gl_migration!=NoneEnum){
    36023713                for(i=0;i<3;i++)nodeinputs[i]=iomodel->bathymetry[tria_vertex_ids[i]-1];
    36033714                this->inputs->AddInput(new TriaVertexInput(BathymetryEnum,nodeinputs));
    … …  
    44224533        GetParameterListOnVertices(&b[0],BedEnum);
    44234534        GetParameterListOnVertices(&ba[0],BathymetryEnum);
    4424         for(i=0;i<3;i++)isonshelf[i]=nodes[i]->IsOnShelf();
     4535        for(i=0;i<3;i++){
     4536                isonshelf[i]=nodes[i]->IsOnShelf();
     4537                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
     4538        }
    44254539
    44264540        /*material parameters: */
    … …  
    44374551                                b[i]=ba[i];
    44384552                                s[i]=b[i]+h[i];
    4439                                 //printf("Node id %i is starting to ground\n",nodes[i]->Id());
     4553                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    44404554                        }
    44414555                        else{
    44424556                                /*do nothing, we are still floating.*/
     4557                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    44434558                        }
    44444559                }
    … …  
    44514566                                s[i]=(1-density)*h[i];
    44524567                                b[i]=-density*h[i];
    4453                                 //printf("Node id %i is starting to float\n",nodes[i]->Id());
     4568                                printf("%i\n",nodes[i]->Sid()+1);
     4569                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    44544570                        }
    44554571                        else{
    44564572                                /*do nothing, we are still grounded.*/
     4573                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still grounded %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    44574574                        }
    44584575                }
    … …  
    44624579        surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
    44634580        bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
    4464 
     4581       
     4582        for(i=0;i<3;i++){
     4583                isonshelf[i]=nodes[i]->IsOnShelf();
     4584                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i second shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
     4585        }
    44654586}
    44664587/*}}}*/
    … …  
    53445465        bool    elementonshelf = false;
    53455466        int     flipped=0;
    5346         bool    shelfstatus[3];
     5467        int     shelfstatus[3];
    53475468
    53485469
    … …  
    53645485        /*Initialize current status of nodes: */
    53655486        for(i=0;i<3;i++){
    5366                 shelfstatus[i]=nodes[i]->inputs->GetInput(NodeOnIceShelfEnum);
     5487                shelfstatus[i]=nodes[i]->IsOnShelf();
     5488                if((nodes[i]->Sid()+1)==36) printf("UpdateShelfStatus: El %i Node %i shelf status %i\n",this->Id(),nodes[i]->Sid()+1,shelfstatus[i]);
    53675489        }
    53685490       
    … …  
    53705492        flipped=0;
    53715493        for(i=0;i<3;i++){
    5372                 if(b[i]<=ba[i]){
     5494                if(b[i]<=ba[i]){ //the = will lead to oscillations.
    53735495                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false));
    53745496                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true));
    5375                         if(shelfstatus[i])flipped++;
    5376                         //printf("node %i flipping to sheet\n",nodes[i]->Sid());
     5497                        if(shelfstatus[i]){
     5498                                flipped++;
     5499                                if((nodes[i]->Sid()+1)==36)printf("UpdateShelfStatus: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     5500                        }
    53775501                        VecSetValue(new_shelf_nodes,nodes[i]->Sid(),0,INSERT_VALUES);
    53785502                }
    … …  
    53805504                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true));
    53815505                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false));
    5382                         if(!shelfstatus[i])flipped++;
    5383                         //printf("node %i flipping to shelf\n",nodes[i]->Sid());
     5506                        if(!shelfstatus[i]){
     5507                                flipped++;
     5508                                if((nodes[i]->Sid()+1)==36)printf("UpdateShelfStatus: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     5509                        }
    53845510                        VecSetValue(new_shelf_nodes,nodes[i]->Sid(),1,INSERT_VALUES);
    53855511                }
    … …  
    54055531        bool flag;
    54065532        int  i;
     5533        double  h[3];
     5534        double  s[3];
     5535        double  b[3];
     5536        double  ba[3];
     5537
     5538        GetParameterListOnVertices(&h[0],ThicknessEnum);
     5539        GetParameterListOnVertices(&s[0],SurfaceEnum);
     5540        GetParameterListOnVertices(&b[0],BedEnum);
     5541        GetParameterListOnVertices(&ba[0],BathymetryEnum);
    54075542
    54085543        for(i=0;i<3;i++){

  • issm/trunk/src/c/objects/Elements/Tria.h

    r7288 r7319  
    113113                void   MaxVy(double* pmaxvy, bool process_units);
    114114                void   MaxVz(double* pmaxvz, bool process_units);
     115                void   AgressiveMigration();
     116                void   AgressiveShelfSync();
    115117                void   MigrateGroundingLine();
    116118                void   MinVel(double* pminvel, bool process_units);
Note: See TracChangeset for help on using the changeset viewer.