Changeset 22841


Ignore:
Timestamp:
06/13/18 02:10:04 (7 years ago)
Author:
bdef
Message:

NEW:Minor Hydro bug and fixing coupling with friction

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

Legend:

Unmodified
Added
Removed

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

    r22364 r22841  
    7474        bool   isefficientlayer;
    7575        int    hydrology_model;
    76        
     76
    7777        /*Fetch data needed: */
    7878        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
    … …  
    161161                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
    162162                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
    163                         }       
     163                        }
    164164                }
    165165        }
    … …  
    233233                active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    234234        }
    235        
     235
    236236        /* Start  looping on the number of gaussian points: */
    237237        Gauss* gauss=basalelement->NewGauss(2);
    … …  
    250250                D[0][0]=D_scalar;
    251251                D[1][1]=D_scalar;
    252                 GetB(B,basalelement,xyz_list,gauss); 
     252                GetB(B,basalelement,xyz_list,gauss);
    253253                TripleMultiply(B,2,numnodes,1,
    254254                                                                         &D[0][0],2,2,0,
    … …  
    265265                                                                                 basis,1,numnodes,0,
    266266                                                                                 &Ke->values[0],1);
    267                        
     267
    268268                        /*Transfer EPL part*/
    269269                        if(isefficientlayer){
    … …  
    288288        delete gauss;
    289289        if(domaintype!=Domain2DhorizontalEnum){
    290                 basalelement->DeleteMaterials(); 
     290                basalelement->DeleteMaterials();
    291291                delete basalelement;
    292292        }
    … …  
    412412        delete gauss;
    413413        if(domaintype!=Domain2DhorizontalEnum){
    414                 basalelement->DeleteMaterials(); 
     414                basalelement->DeleteMaterials();
    415415                delete basalelement;
    416416        }
    … …  
    419419
    420420void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    421         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
     421        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    422422         * For node i, Bi can be expressed in the actual coordinate system
    423          * by: 
     423         * by:
    424424         *       Bi=[ dN/dx ]
    425425         *          [ dN/dy ]
    … …  
    507507                IssmDouble rho_ice        = basalelement->GetMaterialParameter(MaterialsRhoIceEnum);
    508508                IssmDouble g              = basalelement->GetMaterialParameter(ConstantsGEnum);
    509                
     509
    510510                basalelement->GetInputListOnVertices(thickness,ThicknessEnum);
    511511                basalelement->GetInputListOnVertices(base,BaseEnum);
    … …  
    542542        xDelete<int>(doflist);
    543543        if(domaintype!=Domain2DhorizontalEnum){
    544                 basalelement->DeleteMaterials(); 
     544                basalelement->DeleteMaterials();
    545545                delete basalelement;
    546546        }
    … …  
    554554IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/
    555555        int unconf_scheme;
    556         IssmDouble expfac; 
     556        IssmDouble expfac;
    557557        IssmDouble sediment_storing;
    558558        IssmDouble storing,yield;
    … …  
    610610                sed_head_input->GetInputValue(&prestep_head,gauss);
    611611                water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
    612                
     612
    613613                //min definition of the if test
    614614                sediment_transmitivity=FullLayer_transmitivity*min(water_sheet/sediment_thickness,1.);
    … …  
    625625
    626626void  HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
    627        
     627
    628628        int        hmax_flag;
    629629        IssmDouble h_max;
    … …  
    632632        /*Get the flag to the limitation method*/
    633633        element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    634        
     634
    635635        /*Switch between the different cases*/
    636636        switch(hmax_flag){
    … …  
    674674        case 1:
    675675                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    676                 transfer=+leakage; 
     676                transfer=+leakage;
    677677                break;
    678678        default:
    … …  
    698698                _assert_(epl_head_input);
    699699                epl_head_input->GetInputValue(&epl_head,gauss);
    700                 if (element->Id()==42){
    701                         _printf_("epl head in sed Pvec transfer is "<<  epl_head <<"\n");
    702                 }
    703700                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    704701                transfer=+epl_head*leakage;
    … …  
    718715                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    719716                        Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
    720                
     717
    721718                if(node_mask_input->Max()>0.){
    722719                        element_active = true;

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

    r22047 r22841  
    1313#include "../classes.h"
    1414#include "shared/shared.h"
    15 /*}}}*/ 
     15/*}}}*/
    1616
    1717/*Constructors/destructors*/
    … …  
    8888
    8989        switch(CoupledFlag){
     90                //This should be removed at some point and the Enum uniformalized
    9091                case 0:
    9192                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    … …  
    9596                        break;
    9697                case 2:
    97                         element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     98                        element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    9899                        break;
    99100                default:
    … …  
    174175void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
    175176
    176         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p. 
    177          * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
     177        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
     178         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    178179         * alpha_complement= Neff ^r * vel ^s*/
    179180
    … …  
    207208                case 0:
    208209                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    209                         break; 
     210                        break;
    210211                case 1:
    211212                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    212213                        break;
    213214                case 2:
    214                         element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     215                        element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    215216                        break;
    216217                default:
    … …  
    285286void Friction::GetAlpha2Coulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    286287
    287         /*This routine calculates the basal friction coefficient 
     288        /*This routine calculates the basal friction coefficient
    288289          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
    289290
    … …  
    319320                case 0:
    320321                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    321                         break; 
     322                        break;
    322323                case 1:
    323324                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    324325                        break;
    325326                case 2:
    326                         element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     327                        element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    327328                        break;
    328329                default:
    … …  
    370371void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    371372
    372         /*This routine calculates the basal friction coefficient 
     373        /*This routine calculates the basal friction coefficient
    373374                Based on Gagliardini 2007, needs a good effective pressure computation
    374375                Not tested so far so use at your own risks
    … …  
    405406        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    406407        element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    407        
     408
    408409        switch(CoupledFlag){
    409410                case 0:
    410411                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    411                         break; 
     412                        break;
    412413                case 1:
    413414                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    414415                        break;
    415416                case 2:
    416                         element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     417                        element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    417418                        break;
    418419                default:
    … …  
    453454        Gamma=(Chi/(1. + alpha * pow(Chi,q_exp)));
    454455        /*Check to prevent dividing by zero if vmag==0*/
    455         if(vmag==0.) alpha2=0.; 
     456        if(vmag==0.) alpha2=0.;
    456457        else    if (Neff==0) alpha2=0.0;
    457458        else    alpha2=Neff * C_param * pow(Gamma,1./n) * 1/vmag;
    … …  
    570571void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    571572
    572         /*This routine calculates the basal friction coefficient 
     573        /*This routine calculates the basal friction coefficient
    573574          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
    574575
    … …  
    607608                        break;
    608609                case 2:
    609                         element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     610                        element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    610611                        break;
    611612                default:
    … …  
    644645void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    645646
    646         /*This routine calculates the basal friction coefficient 
     647        /*This routine calculates the basal friction coefficient
    647648          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
    648649

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

    r22372 r22841  
    11/*!\file: hydrology_core.cpp
    2  * \brief: core of the hydrology solution 
    3  */ 
     2 * \brief: core of the hydrology solution
     3 */
    44
    55#include "./cores.h"
    … …  
    2323        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    2424        femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
    25         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);       
     25        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    2626        femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum);
    2727        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum);
    … …  
    3636                femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
    3737                solutionsequence_nonlinear(femmodel,modify_loads);
    38                
     38
    3939                /*transfer water column thickness to old water column thickness: */
    4040                InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
    41                
     41
    4242        }
    4343
    … …  
    9090                                        /*save preceding timestep*/
    9191                                        InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
    92                                         InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
    93                                         InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
    9492                                        /*Proceed now to heads computations*/
    9593                                        solutionsequence_hydro_nonlinear(femmodel);
    … …  
    110108                }
    111109        }
    112         else if (hydrology_model==HydrologysommersEnum){       
    113                 femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum);       
    114       InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); 
    115                 solutionsequence_shakti_nonlinear(femmodel); 
     110        else if (hydrology_model==HydrologysommersEnum){
     111                femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum);
     112      InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum);
     113                solutionsequence_shakti_nonlinear(femmodel);
    116114                if(VerboseSolution()) _printf0_("   updating gap height\n");
    117115                HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis();
    … …  
    126124                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    127125        }
    128         /*Free ressources:*/   
     126        /*Free ressources:*/
    129127        if(numoutputs){
    130128                for(int i=0;i<numoutputs;i++){
    … …  
    134132        }
    135133}
    136 
Note: See TracChangeset for help on using the changeset viewer.