Changeset 14757


Ignore:
Timestamp:
04/25/13 13:07:40 (12 years ago)
Author:
bdef
Message:

CHG: Modification to inclue a second layer in hydrology DC

Location:
issm/trunk-jpl/src
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

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

    r14744 r14757  
    9090        HydrologydcEnum,
    9191        SedimentHeadEnum,
     92        EplHeadEnum,
    9293        HydrologydcSpcsedimentHeadEnum,
    9394        HydrologydcSedimentCompressibilityEnum,

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

    r14746 r14757  
    57995799ElementMatrix* Tria::CreateKMatrixHydrology(void){
    58005800
    5801         int hydrology_model;
     5801        int  hydrology_model;
     5802        bool isefficientlayer;
    58025803        parameters->FindParam(&hydrology_model,HydrologyEnum);
     5804       
    58035805
    58045806        switch(hydrology_model){
    … …  
    58065808                        return CreateKMatrixHydrologyShreve();
    58075809                case HydrologydcEnum:
    5808                         return CreateKMatrixHydrologyDC();
     5810                        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     5811                        if(isefficientlayer)
     5812                                return CreateKMatrixHydrologyDCsediment();
     5813                        else
     5814                                return CreateKMatrixHydrologyDCepl();
    58095815                default:
    58105816                        _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
    … …  
    59175923}
    59185924/*}}}*/
    5919 /*FUNCTION Tria::CreateKMatrixHydrologyDC{{{*/
    5920 ElementMatrix* Tria::CreateKMatrixHydrologyDC(void){
     5925/*FUNCTION Tria::CreateKMatrixHydrologyDCsediment{{{*/
     5926ElementMatrix* Tria::CreateKMatrixHydrologyDCsediment(void){
    59215927
    59225928        /*constants: */
    … …  
    59785984}
    59795985/*}}}*/
     5986/*FUNCTION Tria::CreateKMatrixHydrologyDCepl{{{*/
     5987ElementMatrix* Tria::CreateKMatrixHydrologyDCepl(void){
     5988
     5989        /*constants: */
     5990        const int    numdof=NDOF1*NUMVERTICES;
     5991
     5992        /* Intermediaries */
     5993        IssmDouble  D_scalar,Jdet;
     5994        IssmDouble      epl_transmitivity,dt;
     5995        IssmDouble  epl_storing;
     5996        IssmDouble  xyz_list[NUMVERTICES][3];
     5997        IssmDouble  B[2][numdof];
     5998        IssmDouble  L[NUMVERTICES];
     5999        IssmDouble  D[2][2];
     6000        GaussTria   *gauss = NULL;
     6001
     6002        /*Initialize Element matrix*/
     6003        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     6004
     6005        /*Retrieve all inputs and parameters*/
     6006        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
     6007        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6008        epl_storing       = matpar->GetEplStoring();
     6009        epl_transmitivity = matpar->GetEplTransmitivity();
     6010
     6011        /* Start looping on the number of gaussian points: */
     6012        gauss=new GaussTria(2);
     6013        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6014
     6015                gauss->GaussPoint(ig);
     6016
     6017                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6018
     6019                /*Diffusivity*/
     6020                D_scalar=epl_transmitivity*gauss->weight*Jdet;
     6021                if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;
     6022                D[0][0]=D_scalar; D[0][1]=0.;
     6023                D[1][0]=0.;       D[1][1]=D_scalar;
     6024                GetBHydro(&B[0][0],&xyz_list[0][0],gauss);
     6025                TripleMultiply(&B[0][0],2,numdof,1,
     6026                                                                         &D[0][0],2,2,0,
     6027                                                                         &B[0][0],2,numdof,0,
     6028                                                                         &Ke->values[0],1);
     6029
     6030                /*Transient*/
     6031                if(reCast<bool,IssmDouble>(dt)){
     6032                        GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     6033                        D_scalar=epl_storing*gauss->weight*Jdet;
     6034                       
     6035                        TripleMultiply(&L[0],numdof,1,0,
     6036                                                                                 &D_scalar,1,1,0,
     6037                                                                                 &L[0],1,numdof,0,
     6038                                                                                 &Ke->values[0],1);
     6039                }
     6040        }
     6041
     6042        /*Clean up and return*/
     6043        delete gauss;
     6044        return Ke;
     6045}
     6046/*}}}*/
    59806047/*FUNCTION Tria::CreatePVectorHydrology{{{*/
    59816048ElementVector* Tria::CreatePVectorHydrology(void){
    59826049
    5983         int hydrology_model;
     6050        int  hydrology_model;
     6051        bool isefficientlayer;
    59846052        parameters->FindParam(&hydrology_model,HydrologyEnum);
    59856053
    … …  
    59886056                        return CreatePVectorHydrologyShreve();
    59896057                case HydrologydcEnum:
    5990                         return CreatePVectorHydrologyDC();
     6058                        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6059                        if(isefficientlayer)
     6060                                return CreatePVectorHydrologyDCsediment();
     6061                        else
     6062                                return CreatePVectorHydrologyDCepl();
    59916063                default:
    59926064                        _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
    … …  
    60446116}
    60456117/*}}}*/
    6046 /*FUNCTION Tria::CreatePVectorHydrologyDC {{{*/
    6047 ElementVector* Tria::CreatePVectorHydrologyDC(void){
     6118/*FUNCTION Tria::CreatePVectorHydrologyDCsediment {{{*/
     6119ElementVector* Tria::CreatePVectorHydrologyDCsediment(void){
    60486120
    60496121        /*Constants*/
    … …  
    60926164                        old_wh_input->GetInputValue(&water_head,gauss);
    60936165                        scalar = Jdet*gauss->weight*water_head*sediment_storing;
     6166                        for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     6167                }
     6168        }
     6169
     6170        /*Clean up and return*/
     6171        delete gauss;
     6172        return pe;
     6173}
     6174/*}}}*/
     6175/*FUNCTION Tria::CreatePVectorHydrologyDCepl {{{*/
     6176ElementVector* Tria::CreatePVectorHydrologyDCepl(void){
     6177
     6178        /*Constants*/
     6179        const int    numdof=NDOF1*NUMVERTICES;
     6180
     6181        /*Intermediaries */
     6182        IssmDouble Jdet;
     6183        IssmDouble xyz_list[NUMVERTICES][3];
     6184        IssmDouble dt,scalar,water_head;
     6185        IssmDouble water_load;
     6186        IssmDouble epl_storing;
     6187        IssmDouble basis[numdof];
     6188        GaussTria* gauss=NULL;
     6189
     6190        /*Initialize Element vector*/
     6191        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     6192
     6193        /*Retrieve all inputs and parameters*/
     6194        epl_storing = matpar->GetEplStoring();
     6195        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     6196        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6197        Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
     6198        Input* old_wh_input=NULL;
     6199       
     6200        if(reCast<bool,IssmDouble>(dt)){
     6201                old_wh_input=inputs->GetInput(EplHeadEnum); _assert_(old_wh_input);
     6202        }
     6203       
     6204        /* Start  looping on the number of gaussian points: */
     6205        gauss=new GaussTria(2);
     6206        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6207
     6208                gauss->GaussPoint(ig);
     6209
     6210                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6211                GetNodalFunctions(basis, gauss);
     6212
     6213                /*Loading term*/
     6214                water_input->GetInputValue(&water_load,gauss);
     6215                scalar = Jdet*gauss->weight*water_load;
     6216                if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
     6217                for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     6218
     6219                /*Transient term*/
     6220                if(reCast<bool,IssmDouble>(dt)){
     6221                        old_wh_input->GetInputValue(&water_head,gauss);
     6222                        scalar = Jdet*gauss->weight*water_head*epl_storing;
    60946223                        for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    60956224                }

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

    r14695 r14757  
    240240                ElementMatrix* CreateKMatrixHydrology(void);
    241241                ElementMatrix* CreateKMatrixHydrologyShreve(void);
    242                 ElementMatrix* CreateKMatrixHydrologyDC(void);
     242                ElementMatrix* CreateKMatrixHydrologyDCsediment(void);
     243                ElementMatrix* CreateKMatrixHydrologyDCepl(void);
    243244                ElementVector* CreatePVectorHydrology(void);
    244245                ElementVector* CreatePVectorHydrologyShreve(void);
    245                 ElementVector* CreatePVectorHydrologyDC(void);
     246                ElementVector* CreatePVectorHydrologyDCsediment(void);
     247                ElementVector* CreatePVectorHydrologyDCepl(void);
    246248                void    CreateHydrologyWaterVelocityInput(void);
    247249                void      GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution);

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

    r14743 r14757  
    6363                iomodel->Constant(&this->sediment_transmitivity,HydrologydcSedimentTransmitivityEnum);
    6464                iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
     65
     66                iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
     67                iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
     68                iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum);
     69                iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum);
    6570        }
    6671        else{
    … …  
    355360}               
    356361/*}}}*/
     362/*FUNCTION Matpar::GetEplStoring {{{*/
     363IssmDouble Matpar::GetEplStoring(){
     364        return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness*
     365    ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity));             
     366}               
     367/*}}}*/
    357368/*FUNCTION Matpar::GetSedimentTransitivity {{{*/
    358369IssmDouble Matpar::GetSedimentTransmitivity(){
    359370        return sediment_transmitivity;           
     371}               
     372/*}}}*/
     373/*FUNCTION Matpar::GetEplTransitivity {{{*/
     374IssmDouble Matpar::GetEplTransmitivity(){
     375        return epl_transmitivity;               
    360376}               
    361377/*}}}*/                 

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

    r14735 r14757  
    4545                IssmDouble  sediment_transmitivity;     
    4646                IssmDouble  water_compressibility;
     47
     48                IssmDouble  epl_compressibility;
     49                IssmDouble  epl_porosity;       
     50                IssmDouble  epl_thickness;
     51                IssmDouble  epl_transmitivity;   
    4752
    4853                /*gia: */
    … …  
    118123                IssmDouble GetHydrologyN();
    119124                IssmDouble GetSedimentStoring();
     125                IssmDouble GetEplStoring();
    120126                IssmDouble GetSedimentTransmitivity();
     127                IssmDouble GetEplTransmitivity();
    121128                IssmDouble TMeltingPoint(IssmDouble pressure);
    122129                IssmDouble PureIceEnthalpy(IssmDouble pressure);

  • issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r14744 r14757  
    9595                case HydrologydcEnum : return "Hydrologydc";
    9696                case SedimentHeadEnum : return "SedimentHead";
     97                case EplHeadEnum : return "EplHead";
    9798                case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead";
    9899                case HydrologydcSedimentCompressibilityEnum : return "HydrologydcSedimentCompressibility";

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp

    r14562 r14757  
    1616        Parameters *parameters = NULL;
    1717        int         hydrology_model;
     18        bool        isefficientlayer;
    1819
    1920        /*Get parameters: */
    … …  
    2829        }
    2930        else if(hydrology_model==HydrologydcEnum){
    30                 //No parameters for now
     31                iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     32                        parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
    3133        }
    3234        else{

  • issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r14744 r14757  
    9696              else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
    9797              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     98              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    9899              else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum;
    99100              else if (strcmp(name,"HydrologydcSedimentCompressibility")==0) return HydrologydcSedimentCompressibilityEnum;
    … …  
    137138              else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
    138139              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    139               else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
    140140         else stage=2;
    141141   }
    142142   if(stage==2){
    143               if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
     143              if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
     144              else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
    144145              else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
    145146              else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
    … …  
    260261              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    261262              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    262               else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
    263263         else stage=3;
    264264   }
    265265   if(stage==3){
    266               if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
     266              if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
     267              else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
    267268              else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum;
    268269              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    … …  
    383384              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
    384385              else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
    385               else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
    386386         else stage=4;
    387387   }
    388388   if(stage==4){
    389               if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
     389              if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
     390              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    390391              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    391392              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    … …  
    506507              else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
    507508              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    508               else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
    509509         else stage=5;
    510510   }
    511511   if(stage==5){
    512               if (strcmp(name,"None")==0) return NoneEnum;
     512              if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
     513              else if (strcmp(name,"None")==0) return NoneEnum;
    513514              else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
    514515              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;

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

    r14591 r14757  
    7878                        }
    7979                }
     80
    8081                else if (hydrology_model==HydrologydcEnum){
    8182                        if(VerboseSolution()) _pprintLine_("   computing water head");

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

    r14643 r14757  
    66classdef initialization
    77        properties (SetAccess=public)
    8                 vx             = NaN;
    9                 vy             = NaN;
    10                 vz             = NaN;
    11                 vel            = NaN;
    12                 pressure       = NaN;
    13                 temperature    = NaN;
    14                 waterfraction  = NaN;
    15                 sediment_head  = NaN;
    16                 watercolumn    = NaN;
     8                vx            = NaN;
     9                vy            = NaN;
     10                vz            = NaN;
     11                vel           = NaN;
     12                pressure      = NaN;
     13                temperature   = NaN;
     14                waterfraction = NaN;
     15                sediment_head = NaN;
     16                epl_head      = NaN;
     17                watercolumn   = NaN;
    1718        end
    1819        methods
    … …  
    5960                                if isa(md.hydrology,'hydrologydc'),
    6061                                        md = checkfield(md,'initialization.sediment_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
     62                                        md = checkfield(md,'initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
    6163                                else
    6264                                        md = checkfield(md,'initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices 1]);
    … …  
    7577                        fielddisplay(obj,'waterfraction','fraction of water in the ice');
    7678                        fielddisplay(obj,'sediment_head','sediment water head of subglacial system [m]');
     79                        fielddisplay(obj,'epl_head','epl water head of subglacial system [m]');
    7780                        fielddisplay(obj,'watercolumn','thickness of subglacial water [m]');
    7881
    … …  
    8689                        WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum);
    8790                        WriteData(fid,'data',obj.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum);
     91                        WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum);
    8892                        WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);
    8993                end % }}}

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

    r14744 r14757  
    789789        return StringToEnum('SedimentHead')[0]
    790790
     791def EplHeadEnum():
     792        """
     793        EPLHEADENUM - Enum of EplHead
     794
     795           Usage:
     796              macro=EplHeadEnum()
     797        """
     798
     799        return StringToEnum('EplHead')[0]
     800
    791801def HydrologydcSpcsedimentHeadEnum():
    792802        """
    … …  
    53375347        """
    53385348
    5339         return 532
    5340 
     5349        return 533
     5350

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

    r14744 r14757  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=532;
     11macro=533;
Note: See TracChangeset for help on using the changeset viewer.