Ignore:
Timestamp:
07/01/13 23:32:16 (12 years ago)
Author:
Eric.Larour
Message:

CHG: moved code around, to decrease the scope, and hopefully increase AD mode computations, which have too many actives.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15377 r15386  
    142142
    143143/*Other*/
    144 /*FUNCTION Penta::AverageOntoPartition {{{*/
    145 void  Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
    146         _error_("Not supported yet!");
    147 }
    148 /*}}}*/
    149144/*FUNCTION Penta::BedNormal {{{*/
    150145void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){
    … …  
    23182313
    23192314}/*}}}*/
    2320 /*FUNCTION Penta::MigrateGroundingLine{{{*/
    2321 void  Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
    2322 
    2323         int     i,migration_style;
    2324         bool    floatingelement = false;
    2325         bool    groundedelement = false;
    2326         IssmDouble  bed_hydro,yts,gl_melting_rate;
    2327         IssmDouble  rho_water,rho_ice,density;
    2328         IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];
    2329         IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
    2330 
    2331         if(!IsOnBed()) return;
    2332 
    2333         /*Recover info at the vertices: */
    2334         parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    2335         parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
    2336         parameters->FindParam(&yts,ConstantsYtsEnum);
    2337         GetInputListOnVertices(&h[0],ThicknessEnum);
    2338         GetInputListOnVertices(&s[0],SurfaceEnum);
    2339         GetInputListOnVertices(&b[0],BedEnum);
    2340         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2341         if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
    2342         rho_water=matpar->GetRhoWater();
    2343         rho_ice=matpar->GetRhoIce();
    2344         density=rho_ice/rho_water;
    2345 
    2346         /*go through vertices, and update inputs, considering them to be PentaVertex type: */
    2347         for(i=0;i<NUMVERTICES;i++){
    2348                 /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
    2349                 if(reCast<bool,IssmDouble>(old_floating_ice[nodes[i]->Sid()])){
    2350                         if(b[i]<=ba[i]){
    2351                                 b[i]=ba[i];
    2352                                 s[i]=b[i]+h[i];
    2353                                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
    2354                                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
    2355                         }
    2356                 }
    2357                 /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
    2358                 /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
    2359                 else{
    2360                         bed_hydro=-density*h[i];
    2361                         if (bed_hydro>ba[i]){
    2362                                 /*Unground only if the element is connected to the ice shelf*/
    2363                                 if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
    2364                                         s[i]=(1-density)*h[i];
    2365                                         b[i]=-density*h[i];
    2366                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2367                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2368                                 }
    2369                                 else if(migration_style==SoftMigrationEnum && reCast<int,IssmDouble>(sheet_ungrounding[nodes[i]->Sid()])){
    2370                                         s[i]=(1-density)*h[i];
    2371                                         b[i]=-density*h[i];
    2372                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2373                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2374                                 }
    2375                                 else{
    2376                                         if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
    2377                                 }
    2378                         }
    2379                 }
    2380         }
    2381 
    2382         /*SubelementMigrationEnum: if one grounded, all grounded*/
    2383         if(migration_style==SubelementMigrationEnum){
    2384                 for(i=0;i<NUMVERTICES;i++){
    2385                         if(nodes[i]->IsGrounded()){
    2386                                 groundedelement=true;
    2387                                 break;
    2388                         }
    2389                 }
    2390                 floatingelement=!groundedelement;
    2391         }
    2392         else{
    2393                 for(i=0;i<NUMVERTICES;i++){
    2394                         if(nodes[i]->IsFloating()){
    2395                                 floatingelement=true;
    2396                                 break;
    2397                         }
    2398                 }
    2399         }
    2400 
    2401    /*Add basal melting rate if element just ungrounded*/
    2402         if(!this->IsFloating() && floatingelement==true){
    2403                 for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
    2404                 this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
    2405         }
    2406 
    2407         /*Update inputs*/
    2408         this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
    2409         this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
    2410    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
    2411 
    2412         /*Recalculate phi*/
    2413         if(migration_style==SubelementMigrationEnum){
    2414                 for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
    2415                 this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
    2416                 this->InputExtrude(GLlevelsetEnum,ElementEnum);
    2417         }
    2418 
    2419         /*Extrude inputs*/
    2420         this->InputExtrude(SurfaceEnum,ElementEnum);
    2421         this->InputExtrude(BedEnum,ElementEnum);
    2422         this->InputExtrude(MaskElementonfloatingiceEnum,ElementEnum);
    2423         this->InputExtrude(MaskVertexonfloatingiceEnum,NodeEnum);
    2424         this->InputExtrude(MaskVertexongroundediceEnum,NodeEnum);
    2425 }
    2426 /*}}}*/
    24272315/*FUNCTION Penta::MinEdgeLength{{{*/
    24282316IssmDouble Penta::MinEdgeLength(IssmDouble xyz_list[6][3]){
    … …  
    25812469        /*clean-up*/
    25822470        delete gauss;
    2583 }
    2584 /*}}}*/
    2585 /*FUNCTION Penta::PotentialUngrounding{{{*/
    2586 void  Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
    2587 
    2588         int     i;
    2589         IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
    2590         IssmDouble  bed_hydro;
    2591         IssmDouble  rho_water,rho_ice,density;
    2592 
    2593         /*material parameters: */
    2594         rho_water=matpar->GetRhoWater();
    2595         rho_ice=matpar->GetRhoIce();
    2596         density=rho_ice/rho_water;
    2597         GetInputListOnVertices(&h[0],ThicknessEnum);
    2598         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2599 
    2600         /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
    2601         for(i=0;i<NUMVERTICES;i++){
    2602                 /*Find if grounded vertices want to start floating*/
    2603                 if (!nodes[i]->IsFloating()){
    2604                         bed_hydro=-density*h[i];
    2605                         if (bed_hydro>ba[i]){
    2606                                 /*Vertex that could potentially unground, flag it*/
    2607                                 potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
    2608                         }
    2609                 }
    2610         }
    26112471}
    26122472/*}}}*/
    … …  
    32063066}
    32073067/*}}}*/
    3208 /*FUNCTION Penta::UpdatePotentialUngrounding{{{*/
    3209 int Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
    3210 
    3211         int i;
    3212         int nflipped=0;
    3213 
    3214         /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
    3215         for(i=0;i<NUMVERTICES;i++){
    3216                 if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
    3217                         vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
    3218 
    3219                         /*If node was not on ice shelf, we flipped*/
    3220                         if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
    3221                                 nflipped++;
    3222                         }
    3223                 }
    3224         }
    3225         return nflipped;
    3226 }
    3227 /*}}}*/
    32283068/*FUNCTION Penta::UpdateConstraints{{{*/
    32293069void Penta::UpdateConstraints(void){
    … …  
    94259265        xDelete<int>(doflist);
    94269266}
    9427 /*{{{*/
     9267/*}}}*/
    94289268/*FUNCTION Penta::HydrologyEPLGetActive {{{*/
    94299269void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
    … …  
    95169356/*}}}*/
    95179357#endif
     9358
     9359#ifdef _HAVE_GROUNDINGLINE_
     9360/*FUNCTION Penta::MigrateGroundingLine{{{*/
     9361void  Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
     9362
     9363        int     i,migration_style;
     9364        bool    floatingelement = false;
     9365        bool    groundedelement = false;
     9366        IssmDouble  bed_hydro,yts,gl_melting_rate;
     9367        IssmDouble  rho_water,rho_ice,density;
     9368        IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];
     9369        IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
     9370
     9371        if(!IsOnBed()) return;
     9372
     9373        /*Recover info at the vertices: */
     9374        parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
     9375        parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
     9376        parameters->FindParam(&yts,ConstantsYtsEnum);
     9377        GetInputListOnVertices(&h[0],ThicknessEnum);
     9378        GetInputListOnVertices(&s[0],SurfaceEnum);
     9379        GetInputListOnVertices(&b[0],BedEnum);
     9380        GetInputListOnVertices(&ba[0],BathymetryEnum);
     9381        if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
     9382        rho_water=matpar->GetRhoWater();
     9383        rho_ice=matpar->GetRhoIce();
     9384        density=rho_ice/rho_water;
     9385
     9386        /*go through vertices, and update inputs, considering them to be PentaVertex type: */
     9387        for(i=0;i<NUMVERTICES;i++){
     9388                /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
     9389                if(reCast<bool,IssmDouble>(old_floating_ice[nodes[i]->Sid()])){
     9390                        if(b[i]<=ba[i]){
     9391                                b[i]=ba[i];
     9392                                s[i]=b[i]+h[i];
     9393                                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
     9394                                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
     9395                        }
     9396                }
     9397                /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
     9398                /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
     9399                else{
     9400                        bed_hydro=-density*h[i];
     9401                        if (bed_hydro>ba[i]){
     9402                                /*Unground only if the element is connected to the ice shelf*/
     9403                                if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
     9404                                        s[i]=(1-density)*h[i];
     9405                                        b[i]=-density*h[i];
     9406                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     9407                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     9408                                }
     9409                                else if(migration_style==SoftMigrationEnum && reCast<int,IssmDouble>(sheet_ungrounding[nodes[i]->Sid()])){
     9410                                        s[i]=(1-density)*h[i];
     9411                                        b[i]=-density*h[i];
     9412                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     9413                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     9414                                }
     9415                                else{
     9416                                        if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
     9417                                }
     9418                        }
     9419                }
     9420        }
     9421
     9422        /*SubelementMigrationEnum: if one grounded, all grounded*/
     9423        if(migration_style==SubelementMigrationEnum){
     9424                for(i=0;i<NUMVERTICES;i++){
     9425                        if(nodes[i]->IsGrounded()){
     9426                                groundedelement=true;
     9427                                break;
     9428                        }
     9429                }
     9430                floatingelement=!groundedelement;
     9431        }
     9432        else{
     9433                for(i=0;i<NUMVERTICES;i++){
     9434                        if(nodes[i]->IsFloating()){
     9435                                floatingelement=true;
     9436                                break;
     9437                        }
     9438                }
     9439        }
     9440
     9441   /*Add basal melting rate if element just ungrounded*/
     9442        if(!this->IsFloating() && floatingelement==true){
     9443                for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
     9444                this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
     9445        }
     9446
     9447        /*Update inputs*/
     9448        this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
     9449        this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
     9450   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
     9451
     9452        /*Recalculate phi*/
     9453        if(migration_style==SubelementMigrationEnum){
     9454                for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
     9455                this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
     9456                this->InputExtrude(GLlevelsetEnum,ElementEnum);
     9457        }
     9458
     9459        /*Extrude inputs*/
     9460        this->InputExtrude(SurfaceEnum,ElementEnum);
     9461        this->InputExtrude(BedEnum,ElementEnum);
     9462        this->InputExtrude(MaskElementonfloatingiceEnum,ElementEnum);
     9463        this->InputExtrude(MaskVertexonfloatingiceEnum,NodeEnum);
     9464        this->InputExtrude(MaskVertexongroundediceEnum,NodeEnum);
     9465}
     9466/*}}}*/
     9467#endif
     9468
     9469#ifdef _HAVE_WRAPPERS_
     9470/*FUNCTION Penta::PotentialUngrounding{{{*/
     9471void  Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
     9472
     9473        int     i;
     9474        IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
     9475        IssmDouble  bed_hydro;
     9476        IssmDouble  rho_water,rho_ice,density;
     9477
     9478        /*material parameters: */
     9479        rho_water=matpar->GetRhoWater();
     9480        rho_ice=matpar->GetRhoIce();
     9481        density=rho_ice/rho_water;
     9482        GetInputListOnVertices(&h[0],ThicknessEnum);
     9483        GetInputListOnVertices(&ba[0],BathymetryEnum);
     9484
     9485        /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
     9486        for(i=0;i<NUMVERTICES;i++){
     9487                /*Find if grounded vertices want to start floating*/
     9488                if (!nodes[i]->IsFloating()){
     9489                        bed_hydro=-density*h[i];
     9490                        if (bed_hydro>ba[i]){
     9491                                /*Vertex that could potentially unground, flag it*/
     9492                                potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
     9493                        }
     9494                }
     9495        }
     9496}
     9497/*}}}*/
     9498/*FUNCTION Penta::AverageOntoPartition {{{*/
     9499void  Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
     9500        _error_("Not supported yet!");
     9501}
     9502/*}}}*/
     9503/*FUNCTION Penta::UpdatePotentialUngrounding{{{*/
     9504int Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
     9505
     9506        int i;
     9507        int nflipped=0;
     9508
     9509        /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
     9510        for(i=0;i<NUMVERTICES;i++){
     9511                if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
     9512                        vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
     9513
     9514                        /*If node was not on ice shelf, we flipped*/
     9515                        if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
     9516                                nflipped++;
     9517                        }
     9518                }
     9519        }
     9520        return nflipped;
     9521}
     9522/*}}}*/
     9523#endif
Note: See TracChangeset for help on using the changeset viewer.