Changeset 23937


Ignore:
Timestamp:
05/28/19 11:03:06 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on GlaDS

Location:
issm/trunk-jpl/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r23936 r23937  
    1414        if(hydrology_model!=HydrologyGlaDSEnum) return;
    1515
    16         _error_("Not implemented");
    17         //IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologyGlaDSAnalysisEnum,P1Enum);
     16        IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spcphi",HydrologyGlaDSAnalysisEnum,P1Enum);
    1817
    1918}/*}}}*/
    … …  
    2726        if(hydrology_model!=HydrologyGlaDSEnum) return;
    2827
    29         _error_("Not implemented");
    30 
    31         /*Create discrete loads for Moulins*/
    32         CreateSingleNodeToElementConnectivity(iomodel);
    33         for(int i=0;i<iomodel->numberofvertices;i++){
    34                 if (iomodel->domaintype!=Domain3DEnum){
    35                         /*keep only this partition's nodes:*/
    36                         if(iomodel->my_vertices[i]){
    37                                 loads->AddObject(new Moulin(i+1,i,iomodel,HydrologyGlaDSAnalysisEnum));
    38                         }
    39                 }
    40                 else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
    41                         if(iomodel->my_vertices[i]){
    42                                 loads->AddObject(new Moulin(i+1,i,iomodel,HydrologyGlaDSAnalysisEnum));
    43                         }       
    44                 }
    45         }
    46         iomodel->DeleteData(1,"md.mesh.vertexonbase");
    47 
    48         /*Deal with Neumann BC*/
    49         int M,N;
    50         int *segments = NULL;
    51         iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
    52 
    53         /*Check that the size seem right*/
    54         _assert_(N==3); _assert_(M>=3);
    55         for(int i=0;i<M;i++){
    56                 if(iomodel->my_elements[segments[i*3+2]-1]){
    57                         loads->AddObject(new Neumannflux(i+1,i,iomodel,segments,HydrologyGlaDSAnalysisEnum));
    58                 }
    59         }
    60 
    61         xDelete<int>(segments);
    62 
     28        /*No extra loads*/
    6329}/*}}}*/
    6430void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
    … …  
    7036        /*Now, do we really want GlaDS?*/
    7137        if(hydrology_model!=HydrologyGlaDSEnum) return;
    72 
    73         _error_("Not implemented");
    7438
    7539        if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    … …  
    10771        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    10872
    109         _error_("not impelemented yet");
     73        iomodel->FetchDataToInput(elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
    11074}/*}}}*/
    11175void HydrologyGlaDSAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    12185
    12286        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
    123         //parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
    124         _error_("not implemented");
     87        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.pressure_melt_coefficient",HydrologyPressureMeltCoefficientEnum));
     88        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.cavity_spacing",HydrologyCavitySpacingEnum));
     89        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.bump_height",HydrologyBumpHeightEnum));
     90        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.ischannels",HydrologyIschannelsEnum));
     91        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum));
     92        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_sheet_width",HydrologyChannelSheetWidthEnum));
     93        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
    12594
    12695  /*Requested outputs*/

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

    r23936 r23937  
    547547        HydrologydcMaskThawedNodeEnum,
    548548        HydrologydcSedimentTransmitivityEnum,
     549        HydrologySheetConductivityEnum,
     550        HydrologyPressureMeltCoefficientEnum,
     551        HydrologyCavitySpacingEnum,
     552        HydrologyIschannelsEnum,
     553        HydrologyChannelConductivityEnum,
     554        HydrologyChannelSheetWidthEnum,
     555        HydrologyEnglacialVoidRatioEnum,
    549556        HydrologyEnglacialInputEnum,
    550557        HydrologydcEplThicknessStackedEnum,

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

    r23936 r23937  
    553553                case HydrologydcMaskThawedNodeEnum : return "HydrologydcMaskThawedNode";
    554554                case HydrologydcSedimentTransmitivityEnum : return "HydrologydcSedimentTransmitivity";
     555                case HydrologySheetConductivityEnum : return "HydrologySheetConductivity";
     556                case HydrologyPressureMeltCoefficientEnum : return "HydrologyPressureMeltCoefficient";
     557                case HydrologyCavitySpacingEnum : return "HydrologyCavitySpacing";
     558                case HydrologyIschannelsEnum : return "HydrologyIschannels";
     559                case HydrologyChannelConductivityEnum : return "HydrologyChannelConductivity";
     560                case HydrologyChannelSheetWidthEnum : return "HydrologyChannelSheetWidth";
     561                case HydrologyEnglacialVoidRatioEnum : return "HydrologyEnglacialVoidRatio";
    555562                case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput";
    556563                case HydrologydcEplThicknessStackedEnum : return "HydrologydcEplThicknessStacked";

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

    r23936 r23937  
    565565              else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
    566566              else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
     567              else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
     568              else if (strcmp(name,"HydrologyPressureMeltCoefficient")==0) return HydrologyPressureMeltCoefficientEnum;
     569              else if (strcmp(name,"HydrologyCavitySpacing")==0) return HydrologyCavitySpacingEnum;
     570              else if (strcmp(name,"HydrologyIschannels")==0) return HydrologyIschannelsEnum;
     571              else if (strcmp(name,"HydrologyChannelConductivity")==0) return HydrologyChannelConductivityEnum;
     572              else if (strcmp(name,"HydrologyChannelSheetWidth")==0) return HydrologyChannelSheetWidthEnum;
     573              else if (strcmp(name,"HydrologyEnglacialVoidRatio")==0) return HydrologyEnglacialVoidRatioEnum;
    567574              else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
    568575              else if (strcmp(name,"HydrologydcEplThicknessStacked")==0) return HydrologydcEplThicknessStackedEnum;
    … …  
    622629              else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
    623630              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    624               else if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum;
    625635              else if (strcmp(name,"SealevelriseCumDeltathickness")==0) return SealevelriseCumDeltathicknessEnum;
    626636              else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
    … …  
    629639              else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
    630640              else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
     641              else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
    635642              else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
    636643              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
    … …  
    745752              else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
    746753              else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
    747               else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
    748758              else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
    749759              else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
    … …  
    752762              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    753763              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     764              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    758765              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    759766              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    … …  
    868875              else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
    869876              else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
    870               else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
    871881              else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
    872882              else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
    … …  
    875885              else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
    876886              else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
     887              else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
    881888              else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
    882889              else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
    … …  
    991998              else if (strcmp(name,"FileParam")==0) return FileParamEnum;
    992999              else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
    993               else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
    9941004              else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum;
    9951005              else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum;
    … …  
    9981008              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    9991009              else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
     1010              else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
    10041011              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    10051012              else if (strcmp(name,"Fset")==0) return FsetEnum;
    … …  
    11141121              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    11151122              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    1116               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"MinVy")==0) return MinVyEnum;
    11171127              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    11181128              else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
    … …  
    11211131              else if (strcmp(name,"Mpi")==0) return MpiEnum;
    11221132              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Mumps")==0) return MumpsEnum;
     1133              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    11271134              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    11281135              else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    … …  
    12371244              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    12381245              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    1239               else if (strcmp(name,"Tria")==0) return TriaEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Tria")==0) return TriaEnum;
    12401250              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    12411251              else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
    … …  
    12441254              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
    12451255              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"VertexLId")==0) return VertexLIdEnum;
     1256              else if (strcmp(name,"VertexLId")==0) return VertexLIdEnum;
    12501257              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    12511258              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;

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

    r23936 r23937  
    77        properties (SetAccess=public)
    88                %Sheet
    9                 pressure_melt_coefficient = 0.; %c_t
     9                pressure_melt_coefficient = 0.;
     10                sheet_conductivity        = NaN;
     11                cavity_spacing            = 0.;
     12                bump_height               = 0.;
    1013
    1114                %Channels
    12                 ischannels          = 0;
     15                ischannels           = 0;
     16                channel_conductivity = 0.;
     17                channel_sheet_width  = 0.;
    1318
    1419                %Other
    15                 requested_outputs = {};
     20                spcphi               = NaN;
     21                englacial_void_ratio = 0.;
     22                moulin_input         = NaN;
     23                requested_outputs    = {};
    1624        end
    1725        methods
    … …  
    3442                        %Sheet parameters
    3543                        self.pressure_melt_coefficient = 7.5e-8; %K/Pa (See table 1 in Erder et al. 2013)
     44                        self.cavity_spacing = 2.; %m
     45                        self.bump_height = .1; %m
    3646
    3747                        %Channel parameters
    3848                        self.ischannels=false;
     49                        self.channel_conductivity = 5.e-2; %Dow's default, Table uses 0.1
     50                        self.channel_sheet_width = 2.; %m
    3951
    4052                        %Other
     53                        self.englacial_void_ratio = 1.e-5;% Dow's default, Table from Werder et al. uses 1e-3;
    4154                        self.requested_outputs={'default'};
    4255                end % }}}
    … …  
    5063                        %Sheet
    5164                        md = checkfield(md,'fieldname','hydrology.pressure_melt_coefficient','numel',[1],'>=',0);       
     65                        md = checkfield(md,'fieldname','hydrology.sheet_conductivity','size',[md.mesh.numberofvertices 1],'>',0,'NaN',1,'Inf',1);       
     66                        md = checkfield(md,'fieldname','hydrology.cavity_spacing','numel',[1],'>',0);   
     67                        md = checkfield(md,'fieldname','hydrology.bump_height','numel',[1],'>',0);     
    5268
    5369                        %Channels
    5470                        md = checkfield(md,'fieldname','hydrology.ischannels','numel',[1],'values',[0 1]);
     71                        md = checkfield(md,'fieldname','hydrology.channel_conductivity','numel',[1],'>',0);     
     72                        md = checkfield(md,'fieldname','hydrology.channel_sheet_width','numel',[1],'>=',0);     
    5573
    5674                        %Other
     75                        md = checkfield(md,'fieldname','hydrology.spcphi','Inf',1,'timeseries',1);
     76                        md = checkfield(md,'fieldname','hydrology.englacial_void_ratio','numel',[1],'>=',0);
     77                        md = checkfield(md,'fieldname','hydrology.moulin_input','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1);
    5778                        md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1);
    5879                end % }}}
    … …  
    6182                        disp(sprintf('      SHEET'));
    6283                        fielddisplay(self,'pressure_melt_coefficient','Pressure melt coefficient (c_t) [K Pa^-1]');
     84                        fielddisplay(self,'sheet_conductivity','sheet conductivity (k) [m^(7/4) kg^(-1/2)]');
     85                        fielddisplay(self,'cavity_spacing','cavity spacing (l_r) [m]');
     86                        fielddisplay(self,'bump_height','typical bump height (h_r) [m]');
    6387                        disp(sprintf('      CHANNELS'));
    6488                        fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no');
     89                        fielddisplay(self,'channel_conductivity','channel conductivity (k_c) [m^(3/2) kg^(-1/2)]');
    6590                        disp(sprintf('      OTHER'));
     91                        fielddisplay(self,'spcphi','Hydraulic potential Dirichlet constraints [Pa]');
     92                        fielddisplay(self,'englacial_void_ratio','englacial void ratio (e_v)');
     93                        fielddisplay(self,'moulin_input','moulin input (Q_s) [m^3/s]');
    6694                        fielddisplay(self,'requested_outputs','additional outputs requested');
    6795                end % }}}
    … …  
    75103                        %Sheet
    76104                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','pressure_melt_coefficient','format','Double');
     105                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_conductivity','DoubleMat','mattype',1);
     106                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','cavity_spacing','format','Double');
     107                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_height','format','Double');
    77108
    78109                        %Channels
    79110                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','ischannels','format','Boolean');
     111                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_conductivity','format','Double');
     112                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_sheet_width','format','Double');
    80113
    81114                        %Others
     115                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spcphi','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     116                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_void_ratio','format','Double');
     117                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','moulin_input','format','DoubleMat','mattype',1);
    82118                        outputs = self.requested_outputs;
    83119                        pos  = find(ismember(outputs,'default'));
    … …  
    90126        end
    91127end
    92 
Note: See TracChangeset for help on using the changeset viewer.