Changeset 19740


Ignore:
Timestamp:
11/17/15 15:03:35 (9 years ago)
Author:
Mathieu Morlighem
Message:

NEW: done with gap height update

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

Legend:

Unmodified
Added
Removed

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

    r19731 r19740  
    4343
    4444        /*Fetch data needed: */
    45         int    hydrology_model;
     45        int    hydrology_model,frictionlaw;
    4646        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
    4747
    … …  
    7373        iomodel->FetchDataToInput(elements,HydrologyEnglacialInputEnum);
    7474        iomodel->FetchDataToInput(elements,HydrologyBumpSpacingEnum);
     75        iomodel->FetchDataToInput(elements,HydrologyBumpHeightEnum);
    7576        iomodel->FetchDataToInput(elements,HydrologyReynoldsEnum);
     77        iomodel->FetchDataToInput(elements,VxEnum);
     78        iomodel->FetchDataToInput(elements,VyEnum);
     79
     80        iomodel->Constant(&frictionlaw,FrictionLawEnum);
     81        /*Friction law variables*/
     82        switch(frictionlaw){
     83                case 8:
     84                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     85                        break;
     86                default:
     87                        _error_("Friction law "<< frictionlaw <<" not supported");
     88        }
    7689}/*}}}*/
    7790void HydrologySommersAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    8598
    8699        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
     100        parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
    87101}/*}}}*/
    88102
    … …  
    145159        /*Intermediaries */
    146160        IssmDouble  Jdet,meltrate,G,dh[2],B,A,n;
    147         IssmDouble  gap,bed,thickness,head;
     161        IssmDouble  gap,bed,thickness,head,ieb;
     162        IssmDouble  lr,br,vx,vy,beta;
     163        IssmDouble  alpha2,frictionheat;
    148164        IssmDouble* xyz_list = NULL;
    149165
    … …  
    168184        Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    169185        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     186        Input* englacial_input      = element->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
     187        Input* vx_input             = element->GetInput(VxEnum);                         _assert_(vx_input);
     188        Input* vy_input             = element->GetInput(VyEnum);                         _assert_(vy_input);
     189        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
     190        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    170191
    171192        /*Get conductivity from inputs*/
    172193        IssmDouble conductivity = GetConductivity(element);
     194
     195        /*Build friction element, needed later: */
     196        Friction* friction=new Friction(element,3);
    173197
    174198        /* Start  looping on the number of gaussian points: */
    … …  
    185209                head_input->GetInputValue(&head,gauss);
    186210                head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
     211                englacial_input->GetInputValue(&ieb,gauss);
     212                lr_input->GetInputValue(&lr,gauss);
     213                br_input->GetInputValue(&br,gauss);
     214                vx_input->GetInputValue(&vx,gauss);
     215                vy_input->GetInputValue(&vy,gauss);
    187216
    188217                /*Get ice A parameter*/
    … …  
    190219                n_input->GetInputValue(&n,gauss);
    191220                A=pow(B,-n);
     221               
     222                /*Compute beta term*/
     223                if(gap<br)
     224                 beta = (br-gap)/lr;
     225                else
     226                 beta = 0.;
     227
     228                /*Compute frictional heat flux*/
     229                friction->GetAlpha2(&alpha2,gauss);
     230                vx_input->GetInputValue(&vx,gauss);
     231                vy_input->GetInputValue(&vy,gauss);
     232                frictionheat=alpha2*(vx*vx+vy*vy);
    192233
    193234                /*Get water and ice pressures*/
    … …  
    196237                _assert_(pressure_water<=pressure_ice);
    197238
    198                 meltrate = 1/latentheat*(G+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
     239                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
    199240                _assert_(meltrate>0.);
    200241
    … …  
    203244                  meltrate*(1/rho_water-1/rho_ice)
    204245                  +A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
     246                  -beta*sqrt(vx*vx+vy*vy)
     247                  +ieb
    205248                  )*basis[i];
    206249        }
    … …  
    209252        xDelete<IssmDouble>(xyz_list);
    210253        xDelete<IssmDouble>(basis);
     254        delete friction;
    211255        delete gauss;
    212256        return pe;
    … …  
    297341        return conductivity;
    298342}/*}}}*/
     343void HydrologySommersAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/
     344
     345
     346        for(int j=0;j<femmodel->elements->Size();j++){
     347                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
     348                UpdateGapHeight(element);
     349        }
     350
     351}/*}}}*/
     352void HydrologySommersAnalysis::UpdateGapHeight(Element* element){/*{{{*/
     353
     354        /*Skip if water or ice shelf element*/
     355        if(element->IsFloating()) return;
     356
     357        /*Intermediaries */
     358        IssmDouble newgap = 0.;
     359        IssmDouble  Jdet,meltrate,G,dh[2],B,A,n,dt;
     360        IssmDouble  gap,bed,thickness,head,ieb;
     361        IssmDouble  lr,br,vx,vy,beta;
     362        IssmDouble  alpha2,frictionheat;
     363        IssmDouble* xyz_list = NULL;
     364
     365
     366        /*Retrieve all inputs and parameters*/
     367        element->GetVerticesCoordinates(&xyz_list);
     368        element->FindParam(&dt,TimesteppingTimeStepEnum);
     369        IssmDouble  latentheat      = element->GetMaterialParameter(MaterialsLatentheatEnum);
     370        IssmDouble  g               = element->GetMaterialParameter(ConstantsGEnum);
     371        IssmDouble  rho_ice         = element->GetMaterialParameter(MaterialsRhoIceEnum);
     372        IssmDouble  rho_water       = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     373        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
     374        Input* head_input           = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
     375        Input* gap_input            = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
     376        Input* thickness_input      = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
     377        Input* base_input           = element->GetInput(BaseEnum);                       _assert_(base_input);
     378        Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
     379        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     380        Input* englacial_input      = element->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
     381        Input* vx_input             = element->GetInput(VxEnum);                         _assert_(vx_input);
     382        Input* vy_input             = element->GetInput(VyEnum);                         _assert_(vy_input);
     383        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
     384        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
     385
     386        /*Get conductivity from inputs*/
     387        IssmDouble conductivity = GetConductivity(element);
     388
     389        /*Build friction element, needed later: */
     390        Friction* friction=new Friction(element,3);
     391
     392        /*Keep track of weights*/
     393        IssmDouble totalweights=0.;
     394
     395        /* Start  looping on the number of gaussian points: */
     396        Gauss* gauss=element->NewGauss(2);
     397        for(int ig=gauss->begin();ig<gauss->end();ig++){
     398                gauss->GaussPoint(ig);
     399
     400                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     401
     402                geothermalflux_input->GetInputValue(&G,gauss);
     403                base_input->GetInputValue(&bed,gauss);
     404                thickness_input->GetInputValue(&thickness,gauss);
     405                gap_input->GetInputValue(&gap,gauss);
     406                head_input->GetInputValue(&head,gauss);
     407                head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
     408                englacial_input->GetInputValue(&ieb,gauss);
     409                lr_input->GetInputValue(&lr,gauss);
     410                br_input->GetInputValue(&br,gauss);
     411                vx_input->GetInputValue(&vx,gauss);
     412                vy_input->GetInputValue(&vy,gauss);
     413
     414                /*Get ice A parameter*/
     415                B_input->GetInputValue(&B,gauss);
     416                n_input->GetInputValue(&n,gauss);
     417                A=pow(B,-n);
     418
     419                /*Compute beta term*/
     420                if(gap<br)
     421                 beta = (br-gap)/lr;
     422                else
     423                 beta = 0.;
     424
     425                /*Compute frictional heat flux*/
     426                friction->GetAlpha2(&alpha2,gauss);
     427                vx_input->GetInputValue(&vx,gauss);
     428                vy_input->GetInputValue(&vy,gauss);
     429                frictionheat=alpha2*(vx*vx+vy*vy);
     430
     431                /*Get water and ice pressures*/
     432                IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
     433                IssmDouble pressure_water = rho_water*g*(head-bed);
     434                _assert_(pressure_water<=pressure_ice);
     435
     436                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
     437                _assert_(meltrate>0.);
     438
     439                newgap += gauss->weight*Jdet*(gap+dt*(
     440                                        meltrate/rho_ice
     441                                        -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
     442                                        +beta*sqrt(vx*vx+vy*vy)
     443                                        ));
     444                totalweights +=gauss->weight*Jdet;
     445        }
     446
     447        /*Divide by connectivity*/
     448        newgap = newgap/totalweights;
     449
     450        /*Add new gap as an input*/
     451        element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
     452
     453        /*Clean up and return*/
     454        xDelete<IssmDouble>(xyz_list);
     455        delete friction;
     456        delete gauss;
     457}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h

    r19725 r19740  
    3333                /*Intermediaries*/
    3434                IssmDouble GetConductivity(Element* element);
     35                void UpdateGapHeight(FemModel* femmodel);
     36                void UpdateGapHeight(Element* element);
    3537};
    3638#endif

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

    r19488 r19740  
    214214                        GetAlpha2Coulomb(palpha2,gauss);
    215215                        break;
     216                case 8:
     217                        GetAlpha2Sommers(palpha2,gauss);
     218                        break;
    216219          default:
    217                         _error_("not supported");
     220                        _error_("Friction law "<< this->law <<" not supported");
    218221        }
    219222
    … …  
    588591        *palpha2=alpha2;
    589592}/*}}}*/
     593void Friction::GetAlpha2Sommers(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     594
     595        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-bed)*/
     596
     597        /*diverse: */
     598        IssmDouble  pressure_ice,pressure_water;
     599        IssmDouble  Neff;
     600        IssmDouble  drag_coefficient;
     601        IssmDouble  bed,thickness,head;
     602        IssmDouble  alpha2;
     603
     604        /*Recover parameters: */
     605        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     606        element->GetInputValue(&bed, gauss,BaseEnum);
     607        element->GetInputValue(&head, gauss,HydrologyHeadEnum);
     608        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
     609        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     610        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     611        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
     612
     613        //From bed and thickness, compute effective pressure when drag is viscous:
     614        pressure_ice   = rho_ice*gravity*thickness;
     615        pressure_water = rho_water*gravity*(head-bed);
     616        Neff=pressure_ice-pressure_water;
     617        if(Neff<0.) Neff=0.;
     618
     619        alpha2=drag_coefficient*drag_coefficient*Neff;
     620        _assert_(!xIsNan<IssmDouble>(alpha2));
     621
     622        /*Assign output pointers:*/
     623        *palpha2=alpha2;
     624}
     625/*}}}*/

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

    r19479 r19740  
    4040                void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
    4141                void  GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss);
     42                void  GetAlpha2Sommers(IssmDouble* palpha2,Gauss* gauss);
    4243};
    4344

  • issm/trunk-jpl/src/c/cores/hydrology_core.cpp

    r19720 r19740  
    1313
    1414        /*intermediary*/
    15         int        hydrology_model;
    16         bool       save_results;
    17         bool       modify_loads=true;
    18         bool       isefficientlayer;
    19        
     15        int  hydrology_model;
     16        bool save_results;
     17        bool modify_loads=true;
     18        bool isefficientlayer;
    2019
    2120        /*first recover parameters common to all solutions*/
    … …  
    8483                femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum);
    8584                solutionsequence_nonlinear(femmodel,modify_loads);
     85                if(VerboseSolution()) _printf0_("   updating gap height\n");
     86                HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis();
     87                analysis->UpdateGapHeight(femmodel);
     88                delete analysis;       
    8689               
    8790                if(save_results){

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp

    r19150 r19740  
    5757        //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
    5858        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    59         InputUpdateFromSolutionx(femmodel,ug);
     59        //InputUpdateFromSolutionx(femmodel,ug);
    6060
    6161        for(;;){
Note: See TracChangeset for help on using the changeset viewer.