Changeset 18633


Ignore:
Timestamp:
10/14/14 12:15:04 (10 years ago)
Author:
bdef
Message:

CHG: fixing hydrological transfer in DC

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

Legend:

Unmodified
Added
Removed

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

    r18581 r18633  
    158158        basalelement->GetVerticesCoordinates(&xyz_list);
    159159        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    160         Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum);          _assert_(thickness_input);
     160
     161        Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input);
    161162        Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
    162163        Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
    163         Input* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum);
    164164        Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);
    165165
    … …  
    172172                gauss           ->GaussPoint(ig);
    173173                basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
    174                 thickness_input ->GetInputValue(&epl_thickness,gauss);
     174                epl_thick_input ->GetInputValue(&epl_thickness,gauss);
    175175
    176176                /*Diffusivity*/
    … …  
    195195                       
    196196                        /*Transfer EPL part*/
    197                         transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,thickness_input,sed_head_input,epl_head_input,sed_trans_input,residual_input);
     197                        transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input);
    198198                        D_scalar=transfer*gauss->weight*Jdet*dt;
    199199                        TripleMultiply(basis,numnodes,1,0,
    … …  
    266266        basalelement->FindParam(&dt,TimesteppingTimeStepEnum); 
    267267
    268         Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
    269         Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
    270         Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
    271         Input* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum);
    272         Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);    _assert_(residual_input);
     268        Input* epl_thick_input   = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input);
     269        Input* sed_head_input    = basalelement->GetInput(SedimentHeadEnum);
     270        Input* epl_head_input    = basalelement->GetInput(EplHeadEnum);
     271        Input* residual_input    = basalelement->GetInput(SedimentHeadResidualEnum);    _assert_(residual_input);
    273272        if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum);            _assert_(old_wh_input);}
    274273
    … …  
    286285                if(dt!=0.){
    287286                        old_wh_input    ->GetInputValue(&water_head,gauss);
    288                         thickness_input ->GetInputValue(&epl_thickness,gauss);
     287                        epl_thick_input ->GetInputValue(&epl_thickness,gauss);
    289288                       
    290289                        /*Dealing with the epl part of the transfer term*/
    291                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,thickness_input,sed_head_input,epl_head_input,sed_trans_input,residual_input);
     290                        transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input);
    292291                        scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(transfer*dt));
    293292                        for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    … …  
    322321void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    323322
     323        bool active_element;
    324324        int domaintype,i;
    325325        Element*   basalelement=NULL;
    … …  
    341341        int numnodes = basalelement->GetNumberOfNodes();
    342342
     343
    343344        /*Fetch dof list and allocate solution vector*/
     345        IssmDouble* sedhead     = xNew<IssmDouble>(numnodes);
     346        IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
     347
    344348        basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    345         IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
    346 
    347         /*Use the dof list to index into the solution vector: */
     349        basalelement->GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
     350
     351        Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
     352        active_element_input->GetInputValue(&active_element);
     353
     354        /* if(!active_element){ */
     355        /*      /\*Keeping thickness to initial value if EPL is not active*\/ */
     356        /*      for(i=0;i<numnodes;i++){ */
     357        /*              eplHeads[i]=sedhead[i]; */
     358        /*      } */
     359        /* } */
     360        /* else{ */
     361                /*Use the dof list to index into the solution vector: */
    348362        for(i=0;i<numnodes;i++){
    349363                eplHeads[i]=solution[doflist[i]];
    350364                if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
    351365        }
     366                //      }
    352367        /*Add input to the element: */
    353368        element->AddBasalInput(EplHeadEnum,eplHeads,P1Enum);
    … …  
    355370        /*Free ressources:*/
    356371        xDelete<IssmDouble>(eplHeads);
     372        xDelete<IssmDouble>(sedhead);
    357373        xDelete<int>(doflist);
    358374        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    … …  
    381397        return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));               
    382398}/*}}}*/
    383 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input){/*{{{*/
     399
     400IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/
    384401       
     402
    385403        int transfermethod;
    386         IssmDouble epl_thickness;
    387         IssmDouble epl_head,sed_head;
    388         IssmDouble sediment_transmitivity;
    389         IssmDouble leakage,residual,transfer;
    390 
    391         IssmDouble sediment_thickness  = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
    392         IssmDouble sediment_storing    = SedimentStoring(element);
    393         IssmDouble epl_specificstoring = EplSpecificStoring(element);           
     404        IssmDouble epl_thickness,sed_residual;
     405        IssmDouble epl_head,sediment_head;
     406        IssmDouble leakage,transfer;
    394407
    395408        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
     409
    396410        /*Switch between the different transfer methods cases*/
    397411        switch(transfermethod){
    … …  
    401415                break;
    402416        case 1:
    403                 _assert_(epl_thick_input);
    404                 _assert_(sed_head_input);
    405                 _assert_(epl_head_input);
    406                 _assert_(sed_trans_input);
    407                 _assert_(residual_input);
    408                 /* get input */
    409                 epl_thick_input->GetInputValue(&epl_thickness,gauss);
    410                 sed_head_input->GetInputValue(&sed_head,gauss);
     417                _assert_(sed_head_input);
     418                _assert_(epl_head_input);
     419                _assert_(residual_input);
     420               
     421                sed_head_input->GetInputValue(&sediment_head,gauss);
    411422                epl_head_input->GetInputValue(&epl_head,gauss);
    412                 sed_trans_input->GetInputValue(&sediment_transmitivity,gauss);
    413                 residual_input->GetInputValue(&residual,gauss);
     423                residual_input->GetInputValue(&sed_residual,gauss);
    414424                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    415                 transfer=(sediment_transmitivity)/(sediment_thickness*leakage);
     425               
     426                if(epl_head>sediment_head){ 
     427                        if(sed_residual>0.0){
     428                                transfer=0.0;
     429                        }
     430                        else{
     431                                transfer=(leakage);
     432                        }
     433                }
     434                else{
     435                        transfer=(leakage);
     436                }
    416437                break;
    417438        default:
    418439                _error_("no case higher than 1 for the Transfer method");
    419440        }
    420        
     441        if(element->Id()==42){
     442                printf("Transferefficient Kmat %e, %e, %e, %i \n",transfer,sediment_head,epl_head,element->Id());
     443        }
    421444        return transfer;
    422445}/*}}}*/
    423 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input){/*{{{*/
     446
     447IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/
    424448
    425449        int transfermethod;
    426         IssmDouble epl_thickness;
     450        IssmDouble sed_residual;
    427451        IssmDouble epl_head,sediment_head;
    428         IssmDouble sediment_transmitivity;
    429         IssmDouble leakage,residual,transfer;
    430 
    431         IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
    432         IssmDouble sediment_storing   = SedimentStoring(element);
    433         IssmDouble epl_specificstoring = EplSpecificStoring(element);           
     452        IssmDouble leakage,transfer;
    434453
    435454        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
     455
    436456        /*Switch between the different transfer methods cases*/
    437457        switch(transfermethod){
    … …  
    441461                break;
    442462        case 1:
    443                 _assert_(epl_thick_input);
    444                 _assert_(sed_head_input);
    445                 _assert_(epl_head_input);
    446                 _assert_(sed_trans_input);
    447                 _assert_(residual_input);
    448                 /* get input */
    449                 epl_thick_input->GetInputValue(&epl_thickness,gauss);
     463                _assert_(sed_head_input);
     464                _assert_(epl_head_input);
     465                _assert_(residual_input);
     466               
    450467                sed_head_input->GetInputValue(&sediment_head,gauss);
    451468                epl_head_input->GetInputValue(&epl_head,gauss);
    452                 sed_trans_input->GetInputValue(&sediment_transmitivity,gauss);
    453                 residual_input->GetInputValue(&residual,gauss);
     469                residual_input->GetInputValue(&sed_residual,gauss);
    454470                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    455471
    456                 transfer=(sediment_transmitivity*sediment_head)/(sediment_thickness*leakage);
     472                if(epl_head>sediment_head){ 
     473                        if(sed_residual>0.0){
     474                                transfer=0.0;
     475                        }
     476                        else{
     477                                transfer=(sediment_head*leakage);
     478                        }
     479                }
     480                else{
     481                        transfer=(sediment_head*leakage);
     482                }
    457483                break;
    458484        default:
    459485                _error_("no case higher than 1 for the Transfer method");
     486        }
     487        if(element->Id()==42){
     488                printf("Transferefficient Pvec %e, %e, %e\n",transfer,sediment_head,epl_head);
    460489        }
    461490        return transfer;
    … …  
    519548                element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
    520549                element->GetInputListOnVertices(&bed[0],BaseEnum);
    521                        
     550               
    522551                if(!active_element){
    523552                       

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

    r18057 r18633  
    3636                IssmDouble EplSpecificStoring(Element* element);
    3737                IssmDouble SedimentStoring(Element* element);
    38                 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input);
    39                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input);
     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);
    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

    r18576 r18633  
    197197        basalelement ->FindParam(&dt,TimesteppingTimeStepEnum);
    198198        basalelement ->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    199 
    200         Input* epl_thick_input   = basalelement->GetInput(HydrologydcEplThicknessEnum);
     199        Input* SedTrans_input    = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
    201200        Input* sed_head_input    = basalelement->GetInput(SedimentHeadEnum);
    202201        Input* epl_head_input    = basalelement->GetInput(EplHeadEnum);
    203         Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
    204         Input* base_input         = basalelement->GetInput(BaseEnum);
    205         Input* SedTrans_input    = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
     202        Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);
     203
    206204        IssmDouble sediment_storing = SedimentStoring(basalelement);
    207205        /*Transfer related Inputs*/
    … …  
    240238                                active_element_input->GetInputValue(&active_element);
    241239                                if(active_element){
    242                                         transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,epl_thick_input,sed_head_input,epl_head_input,SedTrans_input,thickness_input,base_input);
     240                                        transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input);
    243241                                        basalelement->NodalFunctions(&basis[0],gauss);
    244242                                        D_scalar=transfer*gauss->weight*Jdet*dt;
    … …  
    296294        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    297295
     296
    298297        /*Retrieve all inputs and parameters*/
    299298        basalelement->GetVerticesCoordinates(&xyz_list);
    … …  
    301300        basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    302301
    303         Input* epl_thick_input   = basalelement->GetInput(HydrologydcEplThicknessEnum);
    304302        Input* sed_head_input    = basalelement->GetInput(SedimentHeadEnum);
    305303        Input* epl_head_input    = basalelement->GetInput(EplHeadEnum);
    306         Input* sed_trans_input   = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum);
    307         Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
    308         Input* base_input        = basalelement->GetInput(BaseEnum);
    309304        Input* water_input       = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
     305        Input* residual_input    = basalelement->GetInput(SedimentHeadResidualEnum);
    310306        if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                     _assert_(old_wh_input);}
    311307
    … …  
    342338                                active_element_input->GetInputValue(&active_element);
    343339                                if(active_element){
    344                                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,epl_thick_input,sed_head_input,epl_head_input,sed_trans_input,thickness_input,base_input);
     340                                        transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input);
    345341                                }
    346342                                else{
    … …  
    571567}
    572568/*}}}*/
    573 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input){/*{{{*/
     569
     570IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/
    574571
    575572        int transfermethod;
    576         IssmDouble epl_thickness;
    577         IssmDouble epl_head,sed_head;
    578         IssmDouble sediment_transmitivity;
    579         IssmDouble leakage,h_max,transfer;
    580 
    581         IssmDouble sediment_thickness  = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
    582         IssmDouble sediment_storing    = SedimentStoring(element);
    583         IssmDouble epl_specificstoring = EplSpecificStoring(element);           
     573        IssmDouble sed_residual;
     574        IssmDouble epl_head,sediment_head;
     575        IssmDouble leakage,transfer;
    584576
    585577        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
     578
    586579        /*Switch between the different transfer methods cases*/
    587580        switch(transfermethod){
    … …  
    591584                break;
    592585        case 1:
    593                
    594                 _assert_(epl_thick_input);
    595586                _assert_(sed_head_input);
    596587                _assert_(epl_head_input);
    597                 _assert_(sed_trans_input);
    598                 _assert_(thickness_input);
    599                 _assert_(base_input);
    600 
    601                 epl_thick_input->GetInputValue(&epl_thickness,gauss);
    602                 sed_head_input->GetInputValue(&sed_head,gauss);
     588                _assert_(residual_input);
     589               
     590                sed_head_input->GetInputValue(&sediment_head,gauss);
    603591                epl_head_input->GetInputValue(&epl_head,gauss);
    604                 sed_trans_input->GetInputValue(&sediment_transmitivity,gauss);
     592                residual_input->GetInputValue(&sed_residual,gauss);
    605593                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    606594               
    607                 transfer=(sediment_transmitivity)/(sediment_thickness*leakage);
     595                if(epl_head>sediment_head){ 
     596                        if(sediment_head>=400.0){
     597                                transfer=0.0;
     598                        }
     599                        else{
     600                                transfer=(leakage);
     601                        }
     602                }
     603                else{
     604                        transfer=(leakage);
     605                }               
    608606                break;
    609607        default:
    610608                _error_("no case higher than 1 for the Transfer method");
    611609        }
    612        
     610        if(element->Id()==42){
     611                printf("TransferInfficient Kmat %e, %e, %e, %i\n",transfer,sediment_head,epl_head,element->Id());
     612        }
    613613        return transfer;
    614614}/*}}}*/
    615615
    616 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input){/*{{{*/
     616IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/
    617617
    618618        int transfermethod;
    619         IssmDouble epl_thickness;
     619        IssmDouble sed_residual;
    620620        IssmDouble epl_head,sediment_head;
    621         IssmDouble sediment_transmitivity;
    622         IssmDouble leakage,h_max,transfer;
    623 
    624         IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
    625         IssmDouble sediment_storing   = SedimentStoring(element);
    626         IssmDouble epl_specificstoring = EplSpecificStoring(element);           
     621        IssmDouble leakage,transfer;
    627622
    628623        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
     624
    629625        /*Switch between the different transfer methods cases*/
    630626        switch(transfermethod){
    … …  
    634630                break;
    635631        case 1:
    636                
    637                 _assert_(epl_thick_input);
    638632                _assert_(sed_head_input);
    639633                _assert_(epl_head_input);
    640                 _assert_(sed_trans_input);
    641                 _assert_(thickness_input);
    642                 _assert_(base_input);
    643 
    644                 epl_thick_input->GetInputValue(&epl_thickness,gauss);
     634                _assert_(residual_input);
     635               
    645636                sed_head_input->GetInputValue(&sediment_head,gauss);
    646637                epl_head_input->GetInputValue(&epl_head,gauss);
    647                 sed_trans_input->GetInputValue(&sediment_transmitivity,gauss);
     638                residual_input->GetInputValue(&sed_residual,gauss);
    648639                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    649640
    650                 transfer=(sediment_transmitivity*epl_head)/(sediment_thickness*leakage);
     641                if(epl_head>sediment_head){ 
     642                        if(sediment_head>=400.0){
     643                                transfer=0.0;
     644                        }
     645                        else{
     646                                transfer=(epl_head*leakage);
     647                        }
     648                }
     649                else{
     650                        transfer=(epl_head*leakage);
     651                }
    651652                break;
    652653        default:
    653654                _error_("no case higher than 1 for the Transfer method");
     655        }
     656        if(element->Id()==42){
     657                printf("TransferInfficient Pvec %e, %e, %e, %i\n",transfer,sediment_head,epl_head,element->Id());
    654658        }
    655659        return transfer;

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

    r18057 r18633  
    3838                IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss, Input* thickness_input, Input* base_input);
    3939                void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
    40                 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input);
    41                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input);
     40                IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input);
     41                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input);
    4242                void ElementizeEplMask(FemModel* femmodel);
    4343};
Note: See TracChangeset for help on using the changeset viewer.