Changeset 7323


Ignore:
Timestamp:
02/03/11 18:09:00 (14 years ago)
Author:
Eric.Larour
Message:

Implemented soft version of grounding line migration. More testing needed.

Location:
issm/trunk/src/c
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp

    r7312 r7323  
    2424        /*call different migration modules, according to user wishes: */
    2525        if(migration_style==AgressiveMigrationEnum) AgressiveMigration(elements,nodes, vertices,loads,materials, parameters);
    26         else if(migration_style==SoftMigrationEnum) AgressiveMigration(elements,nodes, vertices,loads,materials, parameters);
     26        else if(migration_style==SoftMigrationEnum) SoftMigration(elements,nodes, vertices,loads,materials, parameters);
    2727        else if(migration_style==NoneEnum)_printf_("%s\n","NoneEnum supplied for migration style, doing nothing!");
    2828        else _error_("%s not supported yet!",EnumToString(migration_style));

  • issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h

    r7312 r7323  
    1111
    1212/* local prototypes: */
    13 
    14 void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
    15 void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
    16 
    17 double* CreateElementOnIceShelf(Elements* elements);
    18 bool*   CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf);
    19 double* CreateElementTouchingIceShelf(Elements* elements);
    20 int     UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf);
    21 
     13void       AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
     14void       SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
     15double*    PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters);
     16double*    PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding);
     17bool*      CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf);
     18double*    CreateElementOnIceShelf(Elements* elements);
     19double*    CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf);
     20int        UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf);
     21Vec        CreateNodesOnIceShelf(Nodes* nodes,int analysis_type);
    2222#endif  /* _GROUNDINGLINEMIGRATIONXLOCAL_H */

  • issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp

    r7312 r7323  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
     12void       AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ //{{{1
    1313
    1414        /*Here, whatever grids inside the ice sheet want to unground, we allow -> instantaneous transmission of water through the bedrock: */
    … …  
    2828        for(i=0;i<elements->Size();i++){
    2929                element=(Element*)elements->GetObjectByOffset(i);
    30                 element->AgressiveShelfSync();
    31         }
    32 
    33 }
    34 
    35 bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf){
     30                element->ShelfSync();
     31        }
     32
     33}
     34/*}}}*/
     35
     36void       SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ //{{{1
     37
     38
     39        /*intermediary: */
     40        int     i;
     41        double* potential_sheet_ungrounding=NULL;
     42        double* sheet_ungrounding=NULL;
     43        Element* element=NULL;
     44       
     45       
     46        _printf_(VerboseModule(),"   Migrating grounding line\n");
     47
     48        /*First, figure out which nodes that are on the ice sheet want to unground. Build a flags vector for this (size number of nodes: */
     49        potential_sheet_ungrounding=PotentialSheetUngrounding(elements,nodes,parameters);
     50
     51        /*Now, propagate ice shelf into connex areas of the ice sheet that potentially want to unground: */
     52        sheet_ungrounding=PropagateShelfIntoConnexIceSheet(elements,nodes,parameters,potential_sheet_ungrounding);
     53
     54        /*Now, use the sheet_ungrounding flags to unground the ice sheet (at the same time, take care of grounding elements of the ice shelf
     55         * that want to: */
     56        for(i=0;i<elements->Size();i++){
     57                element=(Element*)elements->GetObjectByOffset(i);
     58                element->SoftMigration(sheet_ungrounding);
     59        }
     60       
     61        /*Synchronize shelf status: */
     62        for(i=0;i<elements->Size();i++){
     63                element=(Element*)elements->GetObjectByOffset(i);
     64                element->ShelfSync();
     65        }
     66
     67        /*free ressouces: */
     68        xfree((void**)&potential_sheet_ungrounding);
     69        xfree((void**)&sheet_ungrounding);
     70
     71}
     72/*}}}*/
     73double*    PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters){ //{{{1
     74
     75        int      i;
     76        Element* element=NULL;
     77        Vec      vec_potential_sheet_ungrounding=NULL;
     78        double*  potential_sheet_ungrounding=NULL;
     79        int      numnods;
     80        int      analysis_type;
     81
     82        /*recover parameters: */
     83        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     84       
     85        /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */
     86        numnods=nodes->NumberOfNodes(analysis_type);
     87        vec_potential_sheet_ungrounding=NewVec(numnods);
     88
     89        /*Loop through elements, and fill vec_potential_sheet_ungrounding: */
     90        for(i=0;i<elements->Size();i++){
     91                element=(Element*)elements->GetObjectByOffset(i);
     92                element->PotentialSheetUngrounding(vec_potential_sheet_ungrounding);
     93        }
     94
     95        /*Assemble vector: */
     96        VecAssemblyBegin(vec_potential_sheet_ungrounding);
     97        VecAssemblyEnd(vec_potential_sheet_ungrounding);
     98
     99        /*Serialize vector: */
     100        VecToMPISerial(&potential_sheet_ungrounding,vec_potential_sheet_ungrounding);
     101
     102        /*free ressouces: */
     103        VecFree(&vec_potential_sheet_ungrounding);
     104
     105        return potential_sheet_ungrounding;
     106}
     107/*}}}*/
     108double*    PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding){ //{{{1
     109
     110        int      i;
     111        Element* element=NULL;
     112        int      numnods;
     113        int      analysis_type;
     114        Vec      vec_nodes_on_iceshelf=NULL;
     115        double*  nodes_on_iceshelf=NULL;
     116        int      nflipped,local_nflipped;
     117        double*  elements_touching_iceshelf=NULL;
     118
     119        /*recover parameters: */
     120        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     121
     122        /*recover vec_nodes_on_iceshelf: */
     123        vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,analysis_type);
     124
     125        nflipped=1; //bootstrap
     126        while(nflipped){
     127               
     128                /*get a list of potential elements that have nodes on ice shelf: */
     129                elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,vec_nodes_on_iceshelf);
     130
     131                /*now,  go through elements_touching_iceshelf, and if they have nodes inside potential_sheet_ungrounding,
     132                 * flag it: */
     133                local_nflipped=0;
     134
     135                /*serialize vec_nodes_on_iceshelf, needed by element's UpdatePotentialSheetUngrounding routine, to figure out if
     136                 * nodes have flipped from grounded to ungrounded: */
     137                VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
     138               
     139                for(i=0;i<elements->Size();i++){
     140                        element=(Element*)elements->GetObjectByOffset(i);
     141                        if(elements_touching_iceshelf[element->Sid()]){
     142                                nflipped+=element->UpdatePotentialSheetUngrounding(potential_sheet_ungrounding,vec_nodes_on_iceshelf,nodes_on_iceshelf);
     143                        }
     144                }
     145               
     146                MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
     147                //_printf_(VerboseConvergence(),"   number of grounded vertices  connected to grounding line: %i\n",nflipped);
     148                _printf_(1,"   number of grounded vertices  connected to grounding line: %i\n",nflipped);
     149
     150                /*Avoid leaks: */
     151                xfree((void**)&nodes_on_iceshelf);
     152        }
     153
     154        /*Serialize: */
     155        VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
     156
     157        /*Free ressources:*/
     158        VecFree(&vec_nodes_on_iceshelf);
     159        xfree((void**)&elements_touching_iceshelf);
     160
     161        return nodes_on_iceshelf;
     162}
     163/*}}}*/
     164
     165
     166bool*      CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf){ //{{{1
    36167
    37168        int      i;
    … …  
    67198        return element_on_gl;
    68199}
    69 
    70 double* CreateElementOnIceShelf(Elements* elements){
     200/*}}}*/
     201double*    CreateElementOnIceShelf(Elements* elements){ //{{{1
    71202
    72203        int i;
    … …  
    97228        return element_on_iceshelf;
    98229}
    99 
    100 
    101 double* CreateElementTouchingIceShelf(Elements* elements){
     230/*}}}*/
     231double*    CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf){ //{{{1
    102232
    103233        int i;
    … …  
    105235        Vec vec_element_touching_iceshelf=NULL;
    106236        double* element_touching_iceshelf=NULL;
     237        double* nodes_on_iceshelf=NULL;
     238
     239        /*serialize: */
     240        VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
    107241
    108242        /*Create  vector holding  all the elements IsOnShelf flags: */
    … …  
    112246        for(i=0;i<elements->Size();i++){
    113247                element=(Element*)elements->GetObjectByOffset(i);
    114                 VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelf()?1.0:0.0,INSERT_VALUES);
     248                VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES);
    115249        }
    116250
    … …  
    124258        /*free ressouces: */
    125259        VecFree(&vec_element_touching_iceshelf);
     260        xfree((void**)&nodes_on_iceshelf);
    126261
    127262        return element_touching_iceshelf;
    128263}
    129 
    130 
    131 int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){
    132        
     264/*}}}*/
     265int        UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){//{{{1
     266
    133267        int i;
    134268        Element* element=NULL;
    … …  
    175309        return nflipped;
    176310}
    177 
    178 
    179 void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
    180 
    181         int i,j;
    182         Element* element=NULL;
    183         double*  element_touching_iceshelf=NULL;
    184         int      nflipped;
    185 
    186         _printf_(VerboseModule(),"   Migrating grounding line\n");
    187 
    188         /*Loop until no more nodes and elements  change from grounded to floating: */
    189         nflipped=1; //arbitrary, get things started
    190 
    191         int count=0;
    192         while(nflipped){
    193                
    194                 if (count==4)break;
    195 
    196                 /*reset counter: */
    197                 nflipped=0;
    198 
    199                 /*Create  vector holding  all the elements that touch the ice shelf, by any node: */
    200                 element_touching_iceshelf=CreateElementTouchingIceShelf(elements);
    201 
    202 
    203                 /*Carry out grounding line migration for those elements: */
    204                 for(i=0;i<elements->Size();i++){
    205                         element=(Element*)elements->GetObjectByOffset(i);
    206                         if(element_touching_iceshelf[element->Sid()]) element->MigrateGroundingLine();
    207                 }
    208 
    209                 /*Now, update shelf flags in nodes and elements: */
    210                 nflipped=UpdateShelfStatus(elements,nodes,parameters,element_touching_iceshelf);
    211                 //_printf_(VerboseModule(),"      number of migrated nodes: %i\n",nflipped);
    212                 extern int  my_rank;if(my_rank==0)printf("      number of migrated nodes: %i\n",nflipped);
    213 
    214                 /*avoid memory leaks: */
    215                 xfree((void**)&element_touching_iceshelf);
    216                 count++;
    217         }
    218                
    219        
    220         /*free ressouces: */
    221         xfree((void**)&element_touching_iceshelf);
    222        
    223 }
     311/*}}}*/
     312Vec        CreateNodesOnIceShelf(Nodes* nodes,int analysis_type){ //{{{1
     313
     314        /*output: */
     315        Vec     nodes_on_iceshelf=NULL;
     316
     317        /*intermediary: */
     318        int     numnods;
     319        int     i;
     320        Node*   node=NULL;
     321
     322
     323        /*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
     324        numnods=nodes->NumberOfNodes(analysis_type);
     325        nodes_on_iceshelf=NewVec(numnods);
     326
     327        /*Loop through nodes, and fill nodes_on_iceshelf: */
     328        for(i=0;i<nodes->Size();i++){
     329                node=(Node*)nodes->GetObjectByOffset(i);
     330                if(node->IsOnShelf())VecSetValue(nodes_on_iceshelf,node->Sid()+1,1,INSERT_VALUES);
     331        }
     332
     333        /*Assemble vector: */
     334        VecAssemblyBegin(nodes_on_iceshelf);
     335        VecAssemblyEnd(nodes_on_iceshelf);
     336
     337        return nodes_on_iceshelf;
     338}
     339/*}}}*/

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

    r7319 r7323  
    3535                virtual bool   IsOnShelf()=0;
    3636                virtual bool   IsNodeOnShelf()=0;
     37                virtual bool   IsNodeOnShelfFromFlags(double* flags)=0;
    3738                virtual bool   IsOnBed()=0;
    3839                virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
    … …  
    8485                virtual double TimeAdapt()=0;
    8586                virtual void   AgressiveMigration()=0;
    86                 virtual void   AgressiveShelfSync()=0;
     87                virtual void   SoftMigration(double* sheet_ungrounding)=0;
     88                virtual void   ShelfSync()=0;
     89                virtual void   PotentialSheetUngrounding(Vec potential_sheet_ungrounding)=0;
    8790                virtual void   MigrateGroundingLine()=0;
    8891                virtual int    UpdateShelfStatus(Vec new_shelf_nodes)=0;
    8992                virtual void   UpdateShelfFlags(double* new_shelf_nodes)=0;
     93                virtual int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0;
    9094
    9195                /*Implementation: */

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

    r7319 r7323  
    309309}
    310310/*}}}*/
    311 /*FUNCTION Penta::AgressiveShelfSync{{{1*/
    312 void  Penta::AgressiveShelfSync(void){
     311/*FUNCTION Penta::SoftMigration{{{1*/
     312void  Penta::SoftMigration(double* sheet_ungrounding){
     313        _error_("not supported yet!");
     314}
     315/*}}}*/
     316/*FUNCTION Penta::ShelfSync{{{1*/
     317void  Penta::ShelfSync(void){
    313318        _error_("not supported yet!");
    314319}
    … …  
    52535258}
    52545259/*}}}*/
     5260/*FUNCTION Penta::IsNodeOnShelf {{{1*/
     5261bool   Penta::IsNodeOnShelfFromFlags(double* flags){
     5262
     5263        int  i;
     5264        bool shelf=false;
     5265
     5266        for(i=0;i<6;i++){
     5267                if (flags[nodes[i]->Sid()]){
     5268                        shelf=true;
     5269                        break;
     5270                }
     5271        }
     5272        return shelf;
     5273}
     5274/*}}}*/
    52555275/*FUNCTION Penta::IsOnSurface{{{1*/
    52565276bool Penta::IsOnSurface(void){
    … …  
    55005520        *pnumvertices=NUMVERTICES;
    55015521        *pnumnodes=numnodes;
     5522}
     5523/*}}}*/
     5524/*FUNCTION Penta::PotentialSheetUngrounding{{{1*/
     5525void  Penta::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){
     5526        _error_("not supported yet!");
    55025527}
    55035528/*}}}*/
    … …  
    61506175}
    61516176/*}}}*/
     6177/*FUNCTION Penta::UpdatePotentialSheetUngrounding{{{1*/
     6178int Penta::UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){
     6179        _error_("Not implemented yet");
     6180}
     6181/*}}}*/

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

    r7319 r7323  
    110110                void   MaxVz(double* pmaxvz, bool process_units);
    111111                void   AgressiveMigration();
    112                 void   AgressiveShelfSync();
     112                void   SoftMigration(double* sheet_ungrounding);
     113                void   PotentialSheetUngrounding(Vec potential_sheet_ungrounding);
     114                void   ShelfSync();
    113115                void   MigrateGroundingLine();
    114116                void   MinVel(double* pminvel, bool process_units);
    … …  
    129131                int    UpdateShelfStatus(Vec new_shelf_nodes);
    130132                void   UpdateShelfFlags(double* new_shelf_nodes);
     133                int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf);
    131134                double TimeAdapt();
    132135                int*   GetHorizontalNeighboorSids(void);
    … …  
    239242                bool    IsOnShelf(void);
    240243                bool    IsNodeOnShelf(void);
     244                bool    IsNodeOnShelfFromFlags(double* flags);
    241245                bool    IsOnWater(void);
    242246                double  MinEdgeLength(double xyz_list[6][3]);

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

    r7319 r7323  
    335335
    336336        }
    337 /*}}}*/
    338 /*FUNCTION Tria::AgressiveShelfSync{{{1*/
    339 void  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 }
    388337/*}}}*/
    389338/*FUNCTION Tria::AverageOntoPartition {{{1*/
    … …  
    43484297}
    43494298/*}}}*/
     4299/*FUNCTION Tria::IsNodeOnShelf {{{1*/
     4300bool   Tria::IsNodeOnShelfFromFlags(double* flags){
     4301
     4302        int  i;
     4303        bool shelf=false;
     4304
     4305        for(i=0;i<3;i++){
     4306                if (flags[nodes[i]->Sid()]){
     4307                        shelf=true;
     4308                        break;
     4309                }
     4310        }
     4311        return shelf;
     4312}
     4313/*}}}*/
    43504314/*FUNCTION Tria::IsOnWater {{{1*/
    43514315bool   Tria::IsOnWater(){
    … …  
    46954659}
    46964660/*}}}*/
     4661/*FUNCTION Tria::PotentialSheetUngrounding{{{1*/
     4662void  Tria::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){
     4663
     4664
     4665        double *values         = NULL;
     4666        double  h[3],s[3],b[3],ba[3];
     4667        double  bed_hydro;
     4668        double  rho_water,rho_ice,density;
     4669        int     i;
     4670        bool    elementonshelf = false;
     4671
     4672        /*Recover info at the vertices: */
     4673        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     4674        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     4675        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     4676
     4677        GetParameterListOnVertices(&h[0],ThicknessEnum);
     4678        GetParameterListOnVertices(&s[0],SurfaceEnum);
     4679        GetParameterListOnVertices(&b[0],BedEnum);
     4680        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     4681
     4682        /*material parameters: */
     4683        rho_water=matpar->GetRhoWater();
     4684        rho_ice=matpar->GetRhoIce();
     4685        density=rho_ice/rho_water;
     4686
     4687        /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
     4688        for(i=0;i<3;i++){
     4689                if (!nodes[i]->IsOnShelf()){
     4690                       
     4691                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     4692                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     4693                        bed_hydro=-density*h[i];
     4694                        if (bed_hydro>ba[i]){
     4695                                /*ok, this node wants to unground. flag it: */
     4696                                VecSetValue(potential_sheet_ungrounding,nodes[i]->Sid(),1,INSERT_VALUES);
     4697                        }
     4698                }
     4699        }
     4700}
     4701/*}}}*/
    46974702/*FUNCTION Tria::ProcessResultsUnits{{{1*/
    46984703void  Tria::ProcessResultsUnits(void){
    … …  
    47994804
    48004805}
     4806/*}}}*/
     4807/*FUNCTION Tria::ShelfSync{{{1*/
     4808void  Tria::ShelfSync(void){
     4809
     4810
     4811        double *values         = NULL;
     4812        double  h[3],s[3],b[3],ba[3];
     4813        double  bed_hydro;
     4814        double  rho_water,rho_ice,density;
     4815        int     i;
     4816        bool    elementonshelf = false;
     4817
     4818        /*Recover info at the vertices: */
     4819        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     4820        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     4821        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     4822
     4823        GetParameterListOnVertices(&h[0],ThicknessEnum);
     4824        GetParameterListOnVertices(&s[0],SurfaceEnum);
     4825        GetParameterListOnVertices(&b[0],BedEnum);
     4826        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     4827
     4828        /*material parameters: */
     4829        rho_water=matpar->GetRhoWater();
     4830        rho_ice=matpar->GetRhoIce();
     4831        density=rho_ice/rho_water;
     4832
     4833        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     4834        for(i=0;i<3;i++){
     4835                if(b[i]==ba[i]){
     4836                               
     4837                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false));
     4838                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true));
     4839                }
     4840                else{
     4841                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true));
     4842                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false));
     4843
     4844                }
     4845        }
     4846
     4847        /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
     4848        elementonshelf=false;
     4849        for(i=0;i<3;i++){
     4850                if(nodes[i]->IsOnShelf()){
     4851                        elementonshelf=true;
     4852                        break;
     4853                }
     4854        }
     4855    this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf));
     4856}
     4857/*}}}*/
     4858/*FUNCTION Tria::SoftMigration{{{1*/
     4859void  Tria::SoftMigration(double* sheet_ungrounding){
     4860
     4861
     4862        double *values         = NULL;
     4863        double  h[3],s[3],b[3],ba[3];
     4864        double  bed_hydro;
     4865        double  rho_water,rho_ice,density;
     4866        int     i;
     4867        bool    elementonshelf = false;
     4868
     4869        /*Recover info at the vertices: */
     4870        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     4871        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     4872        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     4873
     4874        GetParameterListOnVertices(&h[0],ThicknessEnum);
     4875        GetParameterListOnVertices(&s[0],SurfaceEnum);
     4876        GetParameterListOnVertices(&b[0],BedEnum);
     4877        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     4878
     4879        /*material parameters: */
     4880        rho_water=matpar->GetRhoWater();
     4881        rho_ice=matpar->GetRhoIce();
     4882        density=rho_ice/rho_water;
     4883
     4884        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     4885        for(i=0;i<3;i++){
     4886                if (nodes[i]->IsOnShelf()){
     4887                        /*This node is on the shelf. See if its bed is going under the bathymetry: */
     4888                        if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
     4889                                /*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: */
     4890                                b[i]=ba[i];
     4891                                s[i]=b[i]+h[i];
     4892                        }
     4893                        else{
     4894                                /*do nothing, we are still floating.*/
     4895                        }
     4896                }
     4897                else{
     4898                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     4899                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     4900                        bed_hydro=-density*h[i];
     4901                        if (bed_hydro>ba[i]){
     4902
     4903                                /*Now, are we connected to the grounding line, if so, go forward, otherwise, bail out: */
     4904                                if(sheet_ungrounding[nodes[i]->Sid()]){
     4905                                        /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
     4906                                        s[i]=(1-density)*h[i];
     4907                                        b[i]=-density*h[i];
     4908                                }
     4909                        }
     4910                        else{
     4911                                /*do nothing, we are still grounded.*/
     4912                        }
     4913                }
     4914        }
     4915
     4916        /*Surface and bed are updated. Update inputs:*/
     4917        surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
     4918        bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
     4919
     4920        }
    48014921/*}}}*/
    48024922/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/
    … …  
    55485668}
    55495669/*}}}*/
     5670/*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/
     5671int Tria::UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){
     5672
     5673        /*intermediary: */
     5674        int i;
     5675        int nflipped=0;
     5676
     5677        /*Ok, go through our 3 nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */
     5678        for(i=0;i<3;i++){
     5679                if (potential_sheet_ungrounding[nodes[i]->Sid()]){
     5680                        VecSetValue(vec_nodes_on_iceshelf,nodes[i]->Sid(),1,INSERT_VALUES);
     5681                }
     5682                /*Figure out if we flipped: */
     5683                if (potential_sheet_ungrounding[nodes[i]->Sid()] != nodes_on_iceshelf[nodes[i]->Sid()])nflipped++;
     5684        }
     5685        return nflipped;
     5686}
     5687/*}}}*/

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

    r7319 r7323  
    8484                bool   IsOnShelf();
    8585                bool   IsNodeOnShelf();
     86                bool   IsNodeOnShelfFromFlags(double* flags);
    8687                bool   IsOnWater();
    8788                void   GetSolutionFromInputs(Vec solution);
    … …  
    114115                void   MaxVz(double* pmaxvz, bool process_units);
    115116                void   AgressiveMigration();
    116                 void   AgressiveShelfSync();
     117                void   SoftMigration(double* sheet_ungrounding);
     118                void   ShelfSync();
     119                void   PotentialSheetUngrounding(Vec potential_sheet_ungrounding);
    117120                void   MigrateGroundingLine();
    118121                void   MinVel(double* pminvel, bool process_units);
    … …  
    132135                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    133136                int    UpdateShelfStatus(Vec new_shelf_nodes);
     137                int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf);
    134138                void   UpdateShelfFlags(double* new_shelf_nodes);
    135139                double TimeAdapt();
Note: See TracChangeset for help on using the changeset viewer.