Changeset 21436


Ignore:
Timestamp:
12/09/16 06:03:42 (8 years ago)
Author:
bdef
Message:

NEW: extending the unconfined treatment to EPL

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

Legend:

Unmodified
Added
Removed

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

    r21434 r21436  
    207207        /* Intermediaries */
    208208        IssmDouble  D_scalar,Jdet,dt;
    209         IssmDouble  epl_thickness;
    210209        IssmDouble  transfer;
     210        IssmDouble  epl_transmitivity;
     211        IssmDouble  epl_storing;
    211212        IssmDouble *xyz_list = NULL;
    212213
    … …  
    225226
    226227        Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input);
    227         Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
    228         Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
    229         Input* thick_input     = basalelement->GetInput(ThicknessEnum);
    230         Input* base_input      = basalelement->GetInput(BaseEnum);
    231 
    232         IssmDouble epl_specificstoring   = EplSpecificStoring(basalelement);
    233         IssmDouble epl_conductivity      = basalelement->GetMaterialParameter(HydrologydcEplConductivityEnum);
     228        Input* epl_head_input   = basalelement->GetInput(EplHeadEnum);  _assert_(epl_head_input);
     229        Input* base_input                       = basalelement->GetInput(BaseEnum); _assert_(base_input);
    234230
    235231        /* Start  looping on the number of gaussian points: */
    236         Gauss* gauss=basalelement->NewGauss(2);
     232        Gauss* gauss                                    = basalelement->NewGauss(2);
    237233        for(int ig=gauss->begin();ig<gauss->end();ig++){
    238234                gauss           ->GaussPoint(ig);
    239235                basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
    240                 epl_thick_input ->GetInputValue(&epl_thickness,gauss);
     236
     237                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_head_input,epl_thick_input,base_input);
     238                epl_storing                             = EplStoring(basalelement,gauss,epl_head_input,epl_thick_input,base_input);
    241239
    242240                /*Diffusivity*/
    243                 D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet;
     241                D_scalar=epl_transmitivity*gauss->weight*Jdet;
    244242                if(dt!=0.) D_scalar=D_scalar*dt;
    245243                D[0][0]=D_scalar;
    … …  
    254252                if(dt!=0.){
    255253                        basalelement->NodalFunctions(&basis[0],gauss);
    256                         D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet;
     254                        D_scalar=epl_storing*gauss->weight*Jdet;
    257255                        TripleMultiply(basis,numnodes,1,0,
    258256                                                &D_scalar,1,1,0,
    … …  
    262260                       
    263261                        /*Transfer EPL part*/
    264                         transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input);
     262                        transfer=GetHydrologyKMatrixTransfer(basalelement);
    265263                        D_scalar=dt*transfer*gauss->weight*Jdet;
    266264                        TripleMultiply(basis,numnodes,1,0,
    … …  
    314312        IssmDouble dt,scalar,water_head;
    315313        IssmDouble water_load,transfer;
    316         IssmDouble epl_thickness;
     314        IssmDouble epl_storing;
    317315        IssmDouble Jdet;
    318316        IssmDouble residual,connectivity;
    … …  
    334332
    335333        Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input);
    336         Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
    337         Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
    338         Input* thick_input     = basalelement->GetInput(ThicknessEnum);
    339         Input* base_input      = basalelement->GetInput(BaseEnum);
     334        Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum); _assert_(sed_head_input);
     335        Input* epl_head_input    = basalelement->GetInput(EplHeadEnum); _assert_(epl_head_input);
    340336        Input* water_input               = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
    341         Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);    _assert_(residual_input);
    342         if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum);            _assert_(old_wh_input);}
    343 
    344         IssmDouble epl_specificstoring = EplSpecificStoring(basalelement);
    345 
     337        Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input);
     338        Input* base_input                        = basalelement->GetInput(BaseEnum); _assert_(base_input);
     339
     340        if(dt!= 0.){
     341                old_wh_input = basalelement->GetInput(EplHeadOldEnum);            _assert_(old_wh_input);
     342        }
    346343        /* Start  looping on the number of gaussian points: */
    347         Gauss* gauss=basalelement->NewGauss(2);
     344        Gauss* gauss           = basalelement->NewGauss(2);
    348345        for(int ig=gauss->begin();ig<gauss->end();ig++){
    349346                gauss->GaussPoint(ig);
    … …  
    351348                basalelement ->NodalFunctions(basis,gauss);
    352349
     350                epl_storing     = EplStoring(basalelement,gauss,epl_head_input,epl_thick_input,base_input);
    353351                /*Loading term*/
    354352                water_input->GetInputValue(&water_load,gauss);
    … …  
    362360                if(dt!=0.){
    363361                        old_wh_input    ->GetInputValue(&water_head,gauss);
    364                         epl_thick_input ->GetInputValue(&epl_thickness,gauss);
    365362                       
    366363                        /*Dealing with the epl part of the transfer term*/
    367                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input);
    368                         scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(dt*transfer));
     364                        transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input);
     365                        scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer));
    369366                        for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    370367                }
    … …  
    458455
    459456/*Intermediaries*/
    460 IssmDouble HydrologyDCEfficientAnalysis::EplSpecificStoring(Element* element){/*{{{*/
    461         IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    462         IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
    463         IssmDouble epl_porosity          = element->GetMaterialParameter(HydrologydcEplPorosityEnum);
    464         IssmDouble epl_compressibility   = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum);
    465         IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
    466         return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));                 
    467 }/*}}}*/
    468 
    469 IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/
    470         IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    471         IssmDouble g                        = element->GetMaterialParameter(ConstantsGEnum);
    472         IssmDouble sediment_porosity        = element->GetMaterialParameter(HydrologydcSedimentPorosityEnum);
    473         IssmDouble sediment_thickness       = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
    474         IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum);
    475         IssmDouble water_compressibility    = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
    476         return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));               
    477 }/*}}}*/
    478 
    479 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
     457IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/
    480458
    481459        int transfermethod;
    … …  
    483461
    484462        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    485 
    486463        /*Switch between the different transfer methods cases*/
    487464        switch(transfermethod){
    … …  
    500477}/*}}}*/
    501478
    502 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
     479IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input){/*{{{*/
    503480
    504481        int transfermethod;
    … …  
    507484
    508485        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    509 
    510486        /*Switch between the different transfer methods cases*/
    511487        switch(transfermethod){
    … …  
    537513       
    538514        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
     515        femmodel->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum);
     516        if(iseplthickcomp==0) return;
    539517
    540518        for(int j=0;j<femmodel->elements->Size();j++){
    541519               
    542520                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
    543                 element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum);
    544                 if(iseplthickcomp==0) return;
     521                /* element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); */
     522                /* if(iseplthickcomp==0) return; */
    545523               
    546524                switch(domaintype){
    … …  
    738716}
    739717/*}}}*/
    740 
    741 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
     718IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input){/*{{{*/
     719        IssmDouble epl_storing;
     720        IssmDouble storing,water_sheet;
     721        IssmDouble base_elev,prestep_head,epl_thick;
     722        IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     723        IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
     724        IssmDouble epl_porosity                                  = element->GetMaterialParameter(HydrologydcEplPorosityEnum);
     725        IssmDouble epl_compressibility   = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum);
     726        IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
     727
     728        base_input->GetInputValue(&base_elev,gauss);
     729        epl_head_input->GetInputValue(&prestep_head,gauss);
     730        epl_thick_input->GetInputValue(&epl_thick,gauss);
     731       
     732        water_sheet=max(0.0,(prestep_head-base_elev));
     733       
     734        storing=rho_freshwater*g*epl_porosity*epl_thick*(water_compressibility+(epl_compressibility/epl_porosity));
     735        //porosity for unconfined region
     736        if (water_sheet<=0.9*epl_thick){
     737                epl_storing=epl_porosity;
     738        }
     739        //continuity ramp
     740        else if((water_sheet<epl_thick) && (water_sheet>0.9*epl_thick)){
     741                epl_storing=(epl_thick-water_sheet)*(epl_porosity-storing)/(0.1*epl_thick)+storing;
     742        }
     743        //storing coefficient for confined
     744        else{
     745                epl_storing=storing;
     746        }
     747        return epl_storing;
     748}/*}}}*/
     749
     750IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input){/*{{{*/
     751        IssmDouble epl_transmitivity;
     752        IssmDouble base_elev,prestep_head,epl_thick;
     753        IssmDouble water_sheet;
     754        IssmDouble epl_conductivity      = element->GetMaterialParameter(HydrologydcEplConductivityEnum);
     755        base_input->GetInputValue(&base_elev,gauss);
     756        epl_head_input->GetInputValue(&prestep_head,gauss);
     757        epl_thick_input->GetInputValue(&epl_thick,gauss);
     758
     759        water_sheet=max(0.0,(prestep_head-base_elev));
     760        epl_transmitivity=epl_conductivity*min(water_sheet,epl_thick);
     761
     762        return epl_transmitivity;
     763}/*}}}*/
     764
     765void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){/*{{{*/
    742766        /*Constants*/
    743767        int      domaintype;
    … …  
    787811        xDelete<IssmDouble>(active);
    788812}
     813/*}}}*/
    789814
    790815void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/

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

    r19091 r21436  
    3636                /*Intermediaries*/
    3737                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    38                 IssmDouble EplSpecificStoring(Element* element);
    39                 IssmDouble SedimentStoring(Element* element);
    40                 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
    41                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
     38                IssmDouble GetHydrologyKMatrixTransfer(Element* element);
     39                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input);
     40                IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input);
     41                IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input);
    4242                void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element);
    4343                void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);

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

    r21433 r21436  
    191191        IssmDouble  D_scalar,Jdet,dt;
    192192        IssmDouble  sediment_transmitivity;
    193         IssmDouble  base_elev;
    194193        IssmDouble  transfer,sediment_storing;
    195194        IssmDouble *xyz_list  = NULL;
    … …  
    213212        Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
    214213        Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum);
    215         Input* epl_head_input = basalelement->GetInput(EplHeadEnum);
    216         Input* thick_input    = basalelement->GetInput(ThicknessEnum);
    217214        Input* base_input     = basalelement->GetInput(BaseEnum);
    218215
    … …  
    255252                                active_element_input->GetInputValue(&active_element);
    256253                                if(active_element){
    257                                         transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input);
     254                                        transfer=GetHydrologyKMatrixTransfer(basalelement);
    258255                                        basalelement->NodalFunctions(&basis[0],gauss);
    259256                                        D_scalar=dt*transfer*gauss->weight*Jdet;
    … …  
    319316        Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum);
    320317        Input* epl_head_input = basalelement->GetInput(EplHeadEnum);
    321         Input* thick_input        = basalelement->GetInput(ThicknessEnum);
    322318        Input* base_input                 = basalelement->GetInput(BaseEnum);
    323319        Input* water_input        = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
    324         if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                  _assert_(old_wh_input);}
     320        if(dt!= 0.){
     321                old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                  _assert_(old_wh_input);
     322        }
    325323        /*Transfer related Inputs*/
    326324        if(isefficientlayer){
    … …  
    365363                                active_element_input->GetInputValue(&active_element);
    366364                                if(active_element){
    367                                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input);
     365                                        transfer=GetHydrologyPVectorTransfer(basalelement,gauss,epl_head_input);
    368366                                }
    369367                                else{
    … …  
    514512}/*}}}*/
    515513
    516 /*Intermediaries*/
    517514IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/
    518515        IssmDouble sediment_storing;
    519         IssmDouble specific_porosity;
     516        IssmDouble storing;
    520517        IssmDouble base_elev,prestep_head,water_sheet;
    521518        IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    … …  
    528525        sed_head_input->GetInputValue(&prestep_head,gauss);
    529526        water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
    530         specific_porosity=sediment_porosity-rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
     527        storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
    531528        //porosity for unconfined region
    532529        if (water_sheet<=0.9*sediment_thickness){
    … …  
    535532        //continuity ramp
    536533        else if((water_sheet<sediment_thickness) && (water_sheet>0.9*sediment_thickness)){
    537                 sediment_storing=(sediment_thickness-water_sheet)*specific_porosity/(0.1*sediment_thickness)+
    538                         rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
     534                sediment_storing=(sediment_thickness-water_sheet)*(sediment_porosity-storing)/(0.1*sediment_thickness)+storing;
    539535        }
    540536        //storing coefficient for confined
    541537        else{
    542                 sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
     538                sediment_storing=storing;
    543539        }
    544540        return sediment_storing;
    … …  
    561557        }
    562558        return sediment_transmitivity;
    563 }/*}}}*/
    564 
    565 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss, Input* thick_input, Input* base_input){/*{{{*/
    566         int        hmax_flag;
    567         IssmDouble h_max;
    568         IssmDouble rho_ice,rho_water;
    569         IssmDouble thickness,bed;
    570         /*Get the flag to the limitation method*/
    571         element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    572        
    573         /*Switch between the different cases*/
    574         switch(hmax_flag){
    575         case 0:
    576                 h_max=1.0e+10;
    577                 break;
    578         case 1:
    579                 element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    580                 break;
    581         case 2:
    582                 rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    583                 rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
    584                 _assert_(thick_input);
    585                 _assert_(base_input);
    586                 /*Compute max*/
    587                 thick_input->GetInputValue(&thickness,gauss);
    588                 base_input->GetInputValue(&bed,gauss);
    589                 h_max=((rho_ice*thickness)/rho_water)+bed;
    590                 break;
    591         case 3:
    592                 _error_("Using normal stress  not supported yet");
    593                 break;
    594         default:
    595                 _error_("no case higher than 3 for SedimentlimitFlag");
    596         }
    597         return h_max;
    598559}/*}}}*/
    599560
    … …  
    634595/*}}}*/
    635596
    636 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
     597IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/
    637598
    638599        int transfermethod;
    … …  
    656617}/*}}}*/
    657618
    658 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
     619IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_head_input){/*{{{*/
    659620
    660621        int transfermethod;

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

    r21398 r21436  
    3636                IssmDouble SedimentStoring(Element* element, Gauss* gauss, Input* sed_head_input, Input* base_input);
    3737                IssmDouble SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input);
    38                 IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss, Input* thickness_input, Input* base_input);
    3938                void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
    40                 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
    41                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
     39                IssmDouble GetHydrologyKMatrixTransfer(Element* element);
     40                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_head_input);
    4241                void ElementizeEplMask(FemModel* femmodel);
    4342};
Note: See TracChangeset for help on using the changeset viewer.