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).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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++){
Note: See TracChangeset for help on using the changeset viewer.