Changeset 20156


Ignore:
Timestamp:
02/13/16 22:17:14 (9 years ago)
Author:
Eric.Larour
Message:

CHG: introduced sealevel in the friction (Neff), the neumann boundary conditions (StressbalanceAnalysis), the grouding
line dynamics (MigrateGroundingLIne of Element.cpp) and in the update of free surfaces (MasstransportAnalysis).

Location:
issm/trunk-jpl/src/c
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r20155 r20156  
    673673        IssmDouble* newthickness   = xNew<IssmDouble>(numnodes);
    674674        IssmDouble* deltathickness = xNew<IssmDouble>(numnodes);
    675         IssmDouble* newbed         = xNew<IssmDouble>(numnodes);
     675        IssmDouble* newbase        = xNew<IssmDouble>(numnodes);
    676676        IssmDouble* newsurface     = xNew<IssmDouble>(numnodes);
    677677        IssmDouble* oldthickness   = xNew<IssmDouble>(numnodes);
    678         IssmDouble* oldbed         = xNew<IssmDouble>(numnodes);
     678        IssmDouble* oldbase        = xNew<IssmDouble>(numnodes);
    679679        IssmDouble* oldsurface     = xNew<IssmDouble>(numnodes);
    680680        IssmDouble* phi            = xNew<IssmDouble>(numnodes);
     681        IssmDouble* sealevel       = xNew<IssmDouble>(numnodes);
    681682
    682683        /*Use the dof list to index into the solution vector: */
    … …  
    689690        }
    690691
    691         /*Get previous bed, thickness and surface*/
    692         basalelement->GetInputListOnNodes(&oldbed[0],BaseEnum);
     692        /*Get previous base, thickness, surfac and current sealevel:*/
     693        basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum);
    693694        basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
    694695        basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
    695696        basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
    696 
    697        
    698         /*What is the delta thickness, useful for Sea-level rise:*/
     697        basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum);
     698
     699        /*What is the delta thickness forcing the sea-level rise core: */
    699700        for(i=0;i<numnodes;i++) deltathickness[i]=newthickness[i]-oldthickness[i];
    700701
    … …  
    705706
    706707        for(i=0;i<numnodes;i++) {
    707                 /*If shelf: hydrostatic equilibrium*/
    708                 if (phi[i]>0.){
    709                         newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
    710                         newbed[i]     = oldbed[i];                 //same bed: do nothing
    711                 }
    712                 else{ //this is an ice shelf
     708                if (phi[i]>0.){ //this is an ice sheet: just add thickness to base.
     709                        newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness
     710                        newbase[i]     = oldbase[i];                 //same base: do nothing
     711                }
     712                else{ //this is an ice shelf: hydrostatic equilibrium*/
    713713                        if(hydroadjustment==AbsoluteEnum){
    714                                 newsurface[i] = newthickness[i]*(1-rho_ice/rho_water);
    715                                 newbed[i]     = newthickness[i]*(-rho_ice/rho_water);
     714                                newsurface[i] = newthickness[i]*(1-rho_ice/rho_water)+sealevel[i];
     715                                newbase[i]     = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
    716716                        }
    717717                        else if(hydroadjustment==IncrementalEnum){
    718                                 newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
    719                                 newbed[i]     = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed               = oldbed + di * dH
     718                                newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH
     719                                newbase[i]     = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base               = oldbed + di * dH
    720720                        }
    721721                        else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
    … …  
    727727        element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
    728728        element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
    729         element->AddBasalInput(BaseEnum,newbed,P1Enum);
     729        element->AddBasalInput(BaseEnum,newbase,P1Enum);
    730730
    731731        /*Free ressources:*/
    732732        xDelete<IssmDouble>(newthickness);
    733         xDelete<IssmDouble>(newbed);
     733        xDelete<IssmDouble>(newbase);
    734734        xDelete<IssmDouble>(newsurface);
    735735        xDelete<IssmDouble>(oldthickness);
    736736        xDelete<IssmDouble>(deltathickness);
    737         xDelete<IssmDouble>(oldbed);
     737        xDelete<IssmDouble>(oldbase);
    738738        xDelete<IssmDouble>(oldsurface);
    739739        xDelete<IssmDouble>(phi);

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r20155 r20156  
    16101610        /*Intermediaries*/
    16111611        int         dim,domaintype;
    1612         IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
     1612        IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
    16131613        IssmDouble  surface_under_water,base_under_water,pressure;
    16141614        IssmDouble *xyz_list = NULL;
    … …  
    16351635        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
    16361636        Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
     1637        Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
    16371638        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    16381639        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    16481649                gauss->GaussPoint(ig);
    16491650                thickness_input->GetInputValue(&thickness,gauss);
    1650                 base_input->GetInputValue(&bed,gauss);
     1651                sealevel_input->GetInputValue(&sealevel,gauss);
     1652                base_input->GetInputValue(&base,gauss);
    16511653                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    16521654                element->NodalFunctions(basis,gauss);
    16531655
    1654                 surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
    1655                 base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
     1656                surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level
     1657                base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
    16561658                water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
    16571659                ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
    … …  
    20702072
    20712073        /*Intermediaries*/
    2072         IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
     2074        IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
    20732075        IssmDouble  surface_under_water,base_under_water,pressure;
    20742076        IssmDouble *xyz_list = NULL;
    … …  
    20862088        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
    20872089        Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
     2090        Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
    20882091        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    20892092        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    21002103                thickness_input->GetInputValue(&thickness,gauss);
    21012104                base_input->GetInputValue(&bed,gauss);
     2105                sealevel_input->GetInputValue(&sealevel,gauss);
    21022106                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    21032107                element->NodalFunctions(basis,gauss);
    21042108
    2105                 surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
    2106                 base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
     2109                surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
     2110                base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
    21072111                water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
    21082112                ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
    … …  
    25842588        /*Intermediaries*/
    25852589        int         dim;
    2586         IssmDouble  Jdet,surface,z,water_pressure,ice_pressure;
     2590        IssmDouble  Jdet,surface,sealevel,z,water_pressure,ice_pressure;
    25872591        IssmDouble  surface_under_water,base_under_water,pressure;
    25882592        IssmDouble* xyz_list       = NULL;
    … …  
    26042608        /*Retrieve all inputs and parameters*/
    26052609        Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
     2610        Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
    26062611        IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    26072612        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    26222627                gauss->GaussPoint(ig);
    26232628                surface_input->GetInputValue(&surface,gauss);
     2629                sealevel_input->GetInputValue(&sealevel,gauss);
    26242630                if(dim==3) z=element->GetZcoord(xyz_list,gauss);
    26252631                else       z=element->GetYcoord(xyz_list,gauss);
    … …  
    26272633                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    26282634
    2629                 water_pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
     2635                water_pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level
    26302636                ice_pressure   = rho_ice*gravity*(surface-z);
    26312637                pressure       = ice_pressure + water_pressure;
    … …  
    38113817        /*Intermediaries*/
    38123818        int         i,dim;
    3813         IssmDouble  Jdet,pressure,surface,z;
     3819        IssmDouble  Jdet,pressure,surface,sealevel,z;
    38143820        IssmDouble      normal[3];
    38153821        IssmDouble *xyz_list       = NULL;
    … …  
    38433849        element->NormalSection(&normal[0],xyz_list_front);
    38443850        Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
     3851        Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
    38453852        IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    38463853        IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
    … …  
    38593866                element->NodalFunctionsVelocity(vbasis,gauss);
    38603867                surface_input->GetInputValue(&surface,gauss);
     3868                sealevel_input->GetInputValue(&sealevel,gauss);
    38613869                if(dim==3) z=element->GetZcoord(xyz_list,gauss);
    38623870                else       z=element->GetYcoord(xyz_list,gauss);
    3863                 pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
     3871                pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level
    38643872
    38653873                for (int i=0;i<vnumnodes;i++){

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

    r20154 r20156  
    16531653        IssmDouble* b       = xNew<IssmDouble>(numvertices);
    16541654        IssmDouble* r       = xNew<IssmDouble>(numvertices);
     1655        IssmDouble* sl      = xNew<IssmDouble>(numvertices);
    16551656
    16561657        /*Recover info at the vertices: */
    … …  
    16611662        GetInputListOnVertices(&b[0],BaseEnum);
    16621663        GetInputListOnVertices(&r[0],BedEnum);
     1664        GetInputListOnVertices(&sl[0],SealevelEnum);
    16631665        GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
    16641666        rho_water   = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    … …  
    16841686                /*Change only if AggressiveMigration or if the ice sheet is in contact with the ocean*/
    16851687                else{ // phi>0
    1686                         bed_hydro=-density*h[i];
     1688                        bed_hydro=-density*h[i]+sl[i];
    16871689                        if (bed_hydro>r[i]){
    16881690                                /*Unground only if the element is connected to the ice shelf*/
    16891691                                if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
    1690                                         s[i]        = (1-density)*h[i];
    1691                                         b[i]        = -density*h[i];
     1692                                        s[i]        = (1-density)*h[i]+sl[i];
     1693                                        b[i]        = -density*h[i]+sl[i];
    16921694                                }
    16931695                                else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){
    1694                                         s[i]        = (1-density)*h[i];
    1695                                         b[i]        = -density*h[i];
     1696                                        s[i]        = (1-density)*h[i]+sl[i];
     1697                                        b[i]        = -density*h[i]+sl[i];
    16961698                                }
    16971699                                else{
    … …  
    17051707        for(i=0;i<numvertices;i++){
    17061708                if(migration_style==SoftMigrationEnum){
    1707                         bed_hydro=-density*h[i];
     1709                        bed_hydro=-density*h[i]+sl[i];
    17081710                        if(phi[i]<0. || bed_hydro<=r[i] || phi_ungrounding[vertices[i]->Pid()]<0.){
    1709                                 phi[i]=h[i]+r[i]/density;
     1711                                phi[i]=h[i]+(r[i]-sl[i])/density;
    17101712                        }
    17111713                }
    1712                 else if(migration_style!=ContactEnum) phi[i]=h[i]+r[i]/density;
     1714                else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density;
    17131715                else{
    17141716                        /*do nothing*/
    … …  
    17271729        xDelete<IssmDouble>(b);
    17281730        xDelete<IssmDouble>(s);
     1731        xDelete<IssmDouble>(sl);
    17291732        xDelete<IssmDouble>(h);
    17301733

  • issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r19740 r20156  
    6767        IssmDouble  Neff;
    6868        IssmDouble  drag_coefficient;
    69         IssmDouble  bed,thickness;
     69        IssmDouble  bed,thickness,sealevel;
    7070        IssmDouble  alpha_complement;
    7171
    … …  
    7575        element->GetInputValue(&thickness, gauss,ThicknessEnum);
    7676        element->GetInputValue(&bed, gauss,BaseEnum);
     77        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    7778        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    7879        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    … …  
    8586
    8687        //From bed and thickness, compute effective pressure when drag is viscous:
    87         Neff=gravity*(rho_ice*thickness+rho_water*bed);
     88        Neff=gravity*(rho_ice*thickness+rho_water*(bed-sealevel));
    8889        if(Neff<0)Neff=0;
    8990
    … …  
    231232        IssmDouble  drag_p, drag_q;
    232233        IssmDouble  Neff;
    233         IssmDouble  thickness,base,bed,floatation_thickness;
     234        IssmDouble  thickness,base,bed,floatation_thickness,sealevel;
    234235        IssmDouble  vx,vy,vz,vmag;
    235236        IssmDouble  drag_coefficient,drag_coefficient_coulomb;
    … …  
    241242        element->GetInputValue(&thickness, gauss,ThicknessEnum);
    242243        element->GetInputValue(&base, gauss,BaseEnum);
     244        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    243245        element->GetInputValue(&bed, gauss,BedEnum);
    244246        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    … …  
    253255
    254256        //From base and thickness, compute effective pressure when drag is viscous:
    255         Neff=gravity*(rho_ice*thickness+rho_water*base);
     257        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    256258        if(Neff<0)Neff=0;
    257259
    … …  
    407409        IssmDouble  drag_p, drag_q;
    408410        IssmDouble  Neff;
    409         IssmDouble  thickness,bed;
     411        IssmDouble  thickness,base,sealevel;
    410412        IssmDouble  vx,vy,vz,vmag;
    411413        IssmDouble  drag_coefficient;
    … …  
    416418        element->GetInputValue(&drag_q,FrictionQEnum);
    417419        element->GetInputValue(&thickness, gauss,ThicknessEnum);
    418         element->GetInputValue(&bed, gauss,BaseEnum);
     420        element->GetInputValue(&base, gauss,BaseEnum);
     421        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    419422        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    420423        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    … …  
    426429        s=1./drag_p;
    427430
    428         //From bed and thickness, compute effective pressure when drag is viscous:
    429         Neff=gravity*(rho_ice*thickness+rho_water*bed);
     431        //From base and thickness, compute effective pressure when drag is viscous:
     432        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    430433        if(Neff<0)Neff=0;
    431434
    … …  
    467470        IssmDouble  drag_p, drag_q;
    468471        IssmDouble  Neff,F;
    469         IssmDouble  thickness,bed;
     472        IssmDouble  thickness,bed,sealevel;
    470473        IssmDouble  vx,vy,vz,vmag;
    471474        IssmDouble  drag_coefficient,water_layer;
    … …  
    478481        element->GetInputValue(&thickness, gauss,ThicknessEnum);
    479482        element->GetInputValue(&bed, gauss,BaseEnum);
     483        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    480484        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    481485        element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
    … …  
    490494        //From bed and thickness, compute effective pressure when drag is viscous:
    491495        if(bed>0) bed=0;
    492         if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*bed;
     496        if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(bed-sealevel);
    493497        else if(water_layer>0) Neff=gravity*rho_ice*thickness*F;
    494498        else _error_("negative water layer thickness");
    … …  
    599603        IssmDouble  Neff;
    600604        IssmDouble  drag_coefficient;
    601         IssmDouble  bed,thickness,head;
     605        IssmDouble  bed,thickness,head,sealevel;
    602606        IssmDouble  alpha2;
    603607
    … …  
    606610        element->GetInputValue(&bed, gauss,BaseEnum);
    607611        element->GetInputValue(&head, gauss,HydrologyHeadEnum);
     612        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    608613        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    609614        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    … …  
    613618        //From bed and thickness, compute effective pressure when drag is viscous:
    614619        pressure_ice   = rho_ice*gravity*thickness;
    615         pressure_water = rho_water*gravity*(head-bed);
     620        pressure_water = rho_water*gravity*(head-bed+sealevel);
    616621        Neff=pressure_ice-pressure_water;
    617622        if(Neff<0.) Neff=0.;
Note: See TracChangeset for help on using the changeset viewer.