Changeset 15015


Ignore:
Timestamp:
05/13/13 14:24:49 (12 years ago)
Author:
bdef
Message:

NEW: Adding the computation of the water transfer between layers

Location:
issm/trunk-jpl/src
Files:
3 added
16 edited

Legend:

Unmodified
Added
Removed

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

    r15006 r15015  
    7979
    8080                else if (hydrology_model==HydrologydcEnum){
    81                                 InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum,SedimentHeadOldEnum);
     81                        InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum,SedimentHeadOldEnum);
    8282                        femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     83
     84                        if(VerboseSolution()) _pprintLine_("   computing water transfer");
    8385                       
     86
    8487                        if(VerboseSolution()) _pprintLine_("   computing water head");
    8588                        solutionsequence_hydro_nonlinear(femmodel);

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r15012 r15015  
    136136                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
    137137                virtual void BasisIntegral(Vector<IssmDouble>* basisg)=0;
     138                virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
    138139                #endif
    139140};

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15012 r15015  
    93449344}
    93459345/*}}}*/
    9346 /*FUNCTION Tria::BasisIntegral {{{*/
     9346/*FUNCTION Penta::BasisIntegral {{{*/
    93479347void Penta::BasisIntegral(Vector<IssmDouble>* basisg){
    93489348        _error_("Hydrological stuff not suported in Penta");
    … …  
    93509350/*}}}*/
    93519351
     9352/*}}}*/
     9353/*FUNCTION Penta::GetHydrologyTransfer{{{*/
     9354void  Penta::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
     9355        _error_("Hydrological stuff not suported in Penta");
     9356}
     9357/*}}}*/
     9358       
     9359
    93529360#endif

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r15012 r15015  
    307307                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    308308                void    BasisIntegral(Vector<IssmDouble>* gbasis);
     309                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    309310                #endif
    310311                #ifdef _HAVE_THERMAL_

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r15012 r15015  
    64666466}
    64676467/*}}}*/
     6468/*FUNCTION Tria::GetHydrologyTransfer{{{*/
     6469void  Tria::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
     6470       
     6471        const int  numdof         = NDOF1 *NUMVERTICES;
     6472        int        *doflist       = NULL;
     6473        bool       isefficientlayer;
     6474        int        transfermethod;
     6475        IssmDouble sed_trans,sed_thick;
     6476        IssmDouble leakage;
     6477        IssmDouble wh_trans[numdof]={0.0};
     6478        IssmDouble storing[numdof];
     6479        IssmDouble epl_head[numdof],sed_head[numdof];
     6480
     6481
     6482        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     6483       
     6484        /*Get the flag to know if the efficient layer is present*/
     6485        this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6486       
     6487        if(!isefficientlayer){
     6488                transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL);
     6489                return;
     6490        }       
     6491        /*Also get the flag to the transfer method (TO DO)*/
     6492        this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
     6493       
     6494        /*Switch between the different transfer methods cases(TO DO)*/
     6495        switch(transfermethod){
     6496        case 0:
     6497                /*Just kipping the transfer to zero*/
     6498                break;
     6499        case 1:
     6500                this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum);
     6501               
     6502                sed_trans = matpar->GetSedimentTransmitivity();
     6503                sed_thick = matpar->GetSedimentThickness();
     6504               
     6505                GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
     6506                GetInputListOnVertices(&epl_head[0],EplHeadEnum);
     6507               
     6508                for(int i=0;i<numdof;i++){
     6509                        if(sed_head[i]>epl_head[i]){
     6510                                storing[i]=matpar->GetSedimentStoring();
     6511                                wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
     6512                        }
     6513                        else{
     6514                                storing[i]=matpar->GetEplStoring();
     6515                                wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
     6516                        }
     6517                }
     6518        default:
     6519                _error_("no case higher than 1 for the Transfer method");
     6520        }
     6521       
     6522        /*Assign output pointer*/
     6523        transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL);
     6524}
     6525
     6526/*}}}*/
    64686527
    64696528#endif

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r15012 r15015  
    254254                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    255255                void    BasisIntegral(Vector<IssmDouble>* gbasis);
     256                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    256257                #endif
    257258                #ifdef _HAVE_BALANCED_

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r14996 r15015  
    15941594}
    15951595/*}}}*/
     1596
     1597void FemModel::HydrologyTransferx(void){ /*{{{*/
     1598
     1599        Vector<IssmDouble>* transferg=NULL;
     1600       
     1601        /*Vector allocation*/
     1602        transferg=new Vector<IssmDouble>(nodes->NumberOfNodes());
     1603       
     1604        for (int i=0;i<elements->Size();i++){
     1605                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1606                element->GetHydrologyTransfer(transferg);
     1607        }
     1608        /*Assemble*/
     1609        transferg->Assemble();
     1610       
     1611        /*Update Inputs*/
     1612        InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,transferg,WaterTransferEnum,NodesEnum);
     1613       
     1614}
     1615/*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r15012 r15015  
    100100                void UpdateConstraintsx(void);
    101101                int  UpdateVertexPositionsx(void);
    102                 void BasisIntegralsx(void);
     102                void BasisIntegralsx(void);             
     103                void HydrologyTransferx(void);
    103104};
    104105

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r15012 r15015  
    370370}               
    371371/*}}}*/
     372/*FUNCTION Matpar::GetSedimentThickness {{{*/
     373IssmDouble Matpar::GetSedimentThickness(){
     374        return sediment_thickness;               
     375}               
     376/*}}}*/
    372377/*FUNCTION Matpar::GetEplTransitivity {{{*/
    373378IssmDouble Matpar::GetEplTransmitivity(){

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r14900 r15015  
    123123                IssmDouble GetEplStoring();
    124124                IssmDouble GetSedimentTransmitivity();
     125                IssmDouble GetSedimentThickness();
    125126                IssmDouble GetEplTransmitivity();
    126127                IssmDouble TMeltingPoint(IssmDouble pressure);

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r14990 r15015  
    104104        HydrologydcSedimentlimitFlagEnum,
    105105        HydrologydcSedimentlimitEnum,
     106        HydrologydcTransferFlagEnum,
     107        HydrologydcLeakageFactorEnum,
    106108        HydrologydcPenaltyFactorEnum,
    107109        HydrologyLayerEnum,
    … …  
    110112        HydrologySedimentKmaxEnum,
    111113        BasisIntegralEnum,
     114        WaterTransferEnum,
    112115        IndependentObjectEnum,
    113116        InversionControlParametersEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r14990 r15015  
    112112                case HydrologydcSedimentlimitFlagEnum : return "HydrologydcSedimentlimitFlag";
    113113                case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit";
     114                case HydrologydcTransferFlagEnum : return "HydrologydcTransferFlag";
     115                case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor";
    114116                case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor";
    115117                case HydrologyLayerEnum : return "HydrologyLayer";
    … …  
    118120                case HydrologySedimentKmaxEnum : return "HydrologySedimentKmax";
    119121                case BasisIntegralEnum : return "BasisIntegral";
     122                case WaterTransferEnum : return "WaterTransfer";
    120123                case IndependentObjectEnum : return "IndependentObject";
    121124                case InversionControlParametersEnum : return "InversionControlParameters";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r14990 r15015  
    112112              else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
    113113              else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
     114              else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
     115              else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
    114116              else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
    115117              else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
    … …  
    118120              else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
    119121              else if (strcmp(name,"BasisIntegral")==0) return BasisIntegralEnum;
     122              else if (strcmp(name,"WaterTransfer")==0) return WaterTransferEnum;
    120123              else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
    121124              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    … …  
    134137              else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
    135138              else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
    136               else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
    137               else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
    138               else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
     142              if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
     143              else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
     144              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
     145              else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
    143146              else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
    144147              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    … …  
    257260              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
    258261              else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
    259               else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
    260               else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
    261               else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
     265              if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
     266              else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
     267              else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
     268              else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
    266269              else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
    267270              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    … …  
    380383              else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
    381384              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    382               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    383               else if (strcmp(name,"Tria")==0) return TriaEnum;
    384               else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"Vertex")==0) return VertexEnum;
     388              if (strcmp(name,"StringParam")==0) return StringParamEnum;
     389              else if (strcmp(name,"Tria")==0) return TriaEnum;
     390              else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
     391              else if (strcmp(name,"Vertex")==0) return VertexEnum;
    389392              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
    390393              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    … …  
    503506              else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
    504507              else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
    505               else if (strcmp(name,"Step")==0) return StepEnum;
    506               else if (strcmp(name,"Time")==0) return TimeEnum;
    507               else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
     511              if (strcmp(name,"Step")==0) return StepEnum;
     512              else if (strcmp(name,"Time")==0) return TimeEnum;
     513              else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum;
     514              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    512515              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    513516              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;

  • issm/trunk-jpl/src/m/classes/hydrologydc.m

    r14838 r15015  
    1616                sedimentlimit_flag       = 0;
    1717                sedimentlimit            = 0;
     18                transfer_flag            = 0;
     19                leakage_factor           = 0;
    1820               
    1921                spcepl_head         = NaN;
    … …  
    4547                        obj.sedimentlimit_flag       = 0;
    4648                        obj.sedimentlimit            = 0;
     49                        obj.transfer_flag            = 0;
     50                        obj.leakage_factor           = 10.0;
    4751
    4852                        obj.epl_compressibility = 1.0e-08;
    … …  
    6872                        md = checkfield(md,'hydrology.penalty_factor','>',0,'numel',1);
    6973                        md = checkfield(md,'hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]);
    70 
     74                        md = checkfield(md,'hydrology.transfer_flag','numel',[1],'values',[0 1]);
     75                       
    7176                        if obj.sedimentlimit_flag==1,
    7277                                md = checkfield(md,'hydrology.sedimentlimit','>',0,'numel',1);
    7378            end
    74 
     79                       
     80                        if obj.transfer_flag==1,
     81                                md = checkfield(md,'hydrology.leakage_factor','>',0,'numel',1);
     82            end
     83                       
    7584                        if obj.isefficientlayer==1,
    7685                               
    … …  
    102111                                fielddisplay(obj,'sedimentlimit','user defined upper limit for the inefficient layer [m]');
    103112            end
     113                        fielddisplay(obj,'transfer_flag',['what kind of transfer method is applied between the layers']);
     114                        disp(sprintf('%55s  0: no transfer',' '));
     115                        disp(sprintf('%55s  1: constant leakage factor: %s',' ','leakage_factor'));
     116                        if obj.transfer_flag==1,
     117                                fielddisplay(obj,'leakage_factor','user defined leakage factor [m]');
     118            end
    104119
    105120                        if obj.isefficientlayer==1,
    … …  
    125140                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
    126141                        WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer');
     142                        WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer');
    127143
    128144                        if obj.sedimentlimit_flag==1,
    129145                                WriteData(fid,'object',obj,'fieldname','sedimentlimit','format','Double');
     146            end
     147
     148                        if obj.transfer_flag==1,
     149                                WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double');
    130150            end
    131151

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r14991 r15015  
    13391339        return StringToEnum('HydrologydcSedimentlimit')[0]
    13401340
     1341def HydrologydcTransferFlagEnum():
     1342        """
     1343        HYDROLOGYDCTRANSFERFLAGENUM - Enum of HydrologydcTransferFlag
     1344
     1345        WARNING: DO NOT MODIFY THIS FILE
     1346                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     1347                                Please read src/c/shared/Enum/README for more information
     1348
     1349           Usage:
     1350              macro=HydrologydcTransferFlagEnum()
     1351        """
     1352
     1353        return StringToEnum('HydrologydcTransferFlag')[0]
     1354
     1355def HydrologydcLeakageFactorEnum():
     1356        """
     1357        HYDROLOGYDCLEAKAGEFACTORENUM - Enum of HydrologydcLeakageFactor
     1358
     1359        WARNING: DO NOT MODIFY THIS FILE
     1360                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     1361                                Please read src/c/shared/Enum/README for more information
     1362
     1363           Usage:
     1364              macro=HydrologydcLeakageFactorEnum()
     1365        """
     1366
     1367        return StringToEnum('HydrologydcLeakageFactor')[0]
     1368
    13411369def HydrologydcPenaltyFactorEnum():
    13421370        """
    … …  
    14231451        return StringToEnum('BasisIntegral')[0]
    14241452
     1453def WaterTransferEnum():
     1454        """
     1455        WATERTRANSFERENUM - Enum of WaterTransfer
     1456
     1457        WARNING: DO NOT MODIFY THIS FILE
     1458                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     1459                                Please read src/c/shared/Enum/README for more information
     1460
     1461           Usage:
     1462              macro=WaterTransferEnum()
     1463        """
     1464
     1465        return StringToEnum('WaterTransfer')[0]
     1466
    14251467def IndependentObjectEnum():
    14261468        """
    … …  
    77217763        """
    77227764
    7723         return 550
    7724 
     7765        return 553
     7766

  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r14991 r15015  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=550;
     11macro=553;
Note: See TracChangeset for help on using the changeset viewer.