Changeset 18645


Ignore:
Timestamp:
10/14/14 16:36:18 (10 years ago)
Author:
bdef
Message:

CHG:fixing transfer in DC

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

Legend:

Unmodified
Added
Removed

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

    r18641 r18645  
    162162        Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
    163163        Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
    164         Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);
     164        Input* thick_input     = basalelement->GetInput(ThicknessEnum);
     165        Input* base_input      = basalelement->GetInput(BaseEnum);
    165166
    166167        IssmDouble epl_specificstoring   = EplSpecificStoring(basalelement);
    … …  
    195196                       
    196197                        /*Transfer EPL part*/
    197                         transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input);
     198                        transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input);
    198199                        D_scalar=transfer*gauss->weight*Jdet*dt;
    199200                        TripleMultiply(basis,numnodes,1,0,
    … …  
    269270        Input* sed_head_input    = basalelement->GetInput(SedimentHeadEnum);
    270271        Input* epl_head_input    = basalelement->GetInput(EplHeadEnum);
     272        Input* thick_input       = basalelement->GetInput(ThicknessEnum);
     273        Input* base_input        = basalelement->GetInput(BaseEnum);
    271274        Input* residual_input    = basalelement->GetInput(SedimentHeadResidualEnum);    _assert_(residual_input);
    272275        if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum);            _assert_(old_wh_input);}
    … …  
    288291                       
    289292                        /*Dealing with the epl part of the transfer term*/
    290                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input);
     293                        transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input);
    291294                        scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(transfer*dt));
    292295                        for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    … …  
    390393}/*}}}*/
    391394
    392 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/
    393        
     395IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
    394396
    395397        int transfermethod;
    396         IssmDouble epl_thickness,sed_residual;
     398        IssmDouble hmax;
    397399        IssmDouble epl_head,sediment_head;
    398400        IssmDouble leakage,transfer;
     401        HydrologyDCInefficientAnalysis* inefanalysis = NULL;
    399402
    400403        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    … …  
    409412                _assert_(sed_head_input);
    410413                _assert_(epl_head_input);
    411                 _assert_(residual_input);
    412                
     414               
     415                inefanalysis = new HydrologyDCInefficientAnalysis();
     416                hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);
     417                delete inefanalysis;
     418
    413419                sed_head_input->GetInputValue(&sediment_head,gauss);
    414420                epl_head_input->GetInputValue(&epl_head,gauss);
    415                 residual_input->GetInputValue(&sed_residual,gauss);
    416421                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    417422               
    418423                if(epl_head>sediment_head){ 
    419                         if(sed_residual>0.0){
     424                        if(sediment_head>=hmax){
    420425                                transfer=0.0;
    421426                        }
    … …  
    431436                _error_("no case higher than 1 for the Transfer method");
    432437        }
    433         /* if(element->Id()==42){ */
    434         /*      printf("Transferefficient Kmat %e, %e, %e, %e \n",transfer,sediment_head,epl_head,sed_residual); */
    435         /* } */
    436438        return transfer;
    437439}/*}}}*/
    438440
    439 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/
     441IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
    440442
    441443        int transfermethod;
    442         IssmDouble sed_residual;
     444        IssmDouble hmax;
    443445        IssmDouble epl_head,sediment_head;
    444446        IssmDouble leakage,transfer;
     447
     448        HydrologyDCInefficientAnalysis* inefanalysis = NULL;
    445449
    446450        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    … …  
    455459                _assert_(sed_head_input);
    456460                _assert_(epl_head_input);
    457                 _assert_(residual_input);
     461
     462                inefanalysis = new HydrologyDCInefficientAnalysis();
     463                hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);
     464                delete inefanalysis;
    458465               
    459466                sed_head_input->GetInputValue(&sediment_head,gauss);
    460467                epl_head_input->GetInputValue(&epl_head,gauss);
    461                 residual_input->GetInputValue(&sed_residual,gauss);
    462468                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    463469
    464470                if(epl_head>sediment_head){ 
    465                         if(sed_residual>0.0){
     471                        if(sediment_head>=hmax){
    466472                                transfer=0.0;
    467473                        }
    … …  
    477483                _error_("no case higher than 1 for the Transfer method");
    478484        }
    479         /* if(element->Id()==42){ */
    480         /*      printf("Transferefficient Pvec %e, %e, %e, %e\n",transfer,sediment_head,epl_head,sed_residual); */
    481         /* } */
    482485        return transfer;
    483486}/*}}}*/

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

    r18633 r18645  
    3636                IssmDouble EplSpecificStoring(Element* element);
    3737                IssmDouble SedimentStoring(Element* element);
    38                 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input);
    39                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input);
     38                IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
     39                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
    4040                void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element);
    4141                void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);

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

    r18643 r18645  
    200200        Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum);
    201201        Input* epl_head_input = basalelement->GetInput(EplHeadEnum);
     202        Input* thick_input    = basalelement->GetInput(ThicknessEnum);
    202203        Input* base_input     = basalelement->GetInput(BaseEnum);
    203         Input* thick_input    = basalelement->GetInput(ThicknessEnum);
    204204
    205205        IssmDouble sediment_storing = SedimentStoring(basalelement);
    … …  
    304304        Input* epl_head_input = basalelement->GetInput(EplHeadEnum);
    305305        Input* water_input    = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
     306        Input* thick_input    = basalelement->GetInput(ThicknessEnum);
    306307        Input* base_input     = basalelement->GetInput(BaseEnum);
    307         Input* thick_input    = basalelement->GetInput(ThicknessEnum);
    308308        if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                     _assert_(old_wh_input);}
    309309
    … …  
    599599                                transfer=0.0;
    600600                        }
    601                         else if(sediment_head>=hmax-50.0){
    602                                 transfer=leakage*1.0e-2;
    603                         }
    604601                        else{
    605602                                transfer=(leakage);
    … …  
    613610                _error_("no case higher than 1 for the Transfer method");
    614611        }
    615         /* if(element->Id()==42){ */
    616         /*      printf("TransferInfficient Kmat %e, %e, %e\n",transfer,sediment_head,epl_head); */
    617         /* } */
    618612        return transfer;
    619613}/*}}}*/
    … …  
    648642                                transfer=0.0;
    649643                        }
    650                         else if(sediment_head>=hmax-50.0){
    651                                 transfer=leakage*1.0e-2;
    652                         }
    653644                        else{
    654645                                transfer=(epl_head*leakage);
    … …  
    662653                _error_("no case higher than 1 for the Transfer method");
    663654        }
    664         /* if(element->Id()==42){ */
    665         /*      printf("TransferInfficient Pvec %e, %e, %e\n",transfer,sediment_head,epl_head); */
    666         /* } */
    667655        return transfer;
    668656}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.