Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23936) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23937) @@ -14,6 +14,5 @@ if(hydrology_model!=HydrologyGlaDSEnum) return; - _error_("Not implemented"); - //IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologyGlaDSAnalysisEnum,P1Enum); + IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spcphi",HydrologyGlaDSAnalysisEnum,P1Enum); }/*}}}*/ @@ -27,38 +26,5 @@ if(hydrology_model!=HydrologyGlaDSEnum) return; - _error_("Not implemented"); - - /*Create discrete loads for Moulins*/ - CreateSingleNodeToElementConnectivity(iomodel); - for(int i=0;inumberofvertices;i++){ - if (iomodel->domaintype!=Domain3DEnum){ - /*keep only this partition's nodes:*/ - if(iomodel->my_vertices[i]){ - loads->AddObject(new Moulin(i+1,i,iomodel,HydrologyGlaDSAnalysisEnum)); - } - } - else if(reCast(iomodel->Data("md.mesh.vertexonbase")[i])){ - if(iomodel->my_vertices[i]){ - loads->AddObject(new Moulin(i+1,i,iomodel,HydrologyGlaDSAnalysisEnum)); - } - } - } - iomodel->DeleteData(1,"md.mesh.vertexonbase"); - - /*Deal with Neumann BC*/ - int M,N; - int *segments = NULL; - iomodel->FetchData(&segments,&M,&N,"md.mesh.segments"); - - /*Check that the size seem right*/ - _assert_(N==3); _assert_(M>=3); - for(int i=0;imy_elements[segments[i*3+2]-1]){ - loads->AddObject(new Neumannflux(i+1,i,iomodel,segments,HydrologyGlaDSAnalysisEnum)); - } - } - - xDelete(segments); - + /*No extra loads*/ }/*}}}*/ void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ @@ -70,6 +36,4 @@ /*Now, do we really want GlaDS?*/ if(hydrology_model!=HydrologyGlaDSEnum) return; - - _error_("Not implemented"); if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); @@ -107,5 +71,5 @@ iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); - _error_("not impelemented yet"); + iomodel->FetchDataToInput(elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum); }/*}}}*/ void HydrologyGlaDSAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -121,6 +85,11 @@ parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); - //parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum)); - _error_("not implemented"); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.pressure_melt_coefficient",HydrologyPressureMeltCoefficientEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.cavity_spacing",HydrologyCavitySpacingEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.bump_height",HydrologyBumpHeightEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.ischannels",HydrologyIschannelsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_sheet_width",HydrologyChannelSheetWidthEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum)); /*Requested outputs*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23936) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23937) @@ -547,4 +547,11 @@ HydrologydcMaskThawedNodeEnum, HydrologydcSedimentTransmitivityEnum, + HydrologySheetConductivityEnum, + HydrologyPressureMeltCoefficientEnum, + HydrologyCavitySpacingEnum, + HydrologyIschannelsEnum, + HydrologyChannelConductivityEnum, + HydrologyChannelSheetWidthEnum, + HydrologyEnglacialVoidRatioEnum, HydrologyEnglacialInputEnum, HydrologydcEplThicknessStackedEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23936) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23937) @@ -553,4 +553,11 @@ case HydrologydcMaskThawedNodeEnum : return "HydrologydcMaskThawedNode"; case HydrologydcSedimentTransmitivityEnum : return "HydrologydcSedimentTransmitivity"; + case HydrologySheetConductivityEnum : return "HydrologySheetConductivity"; + case HydrologyPressureMeltCoefficientEnum : return "HydrologyPressureMeltCoefficient"; + case HydrologyCavitySpacingEnum : return "HydrologyCavitySpacing"; + case HydrologyIschannelsEnum : return "HydrologyIschannels"; + case HydrologyChannelConductivityEnum : return "HydrologyChannelConductivity"; + case HydrologyChannelSheetWidthEnum : return "HydrologyChannelSheetWidth"; + case HydrologyEnglacialVoidRatioEnum : return "HydrologyEnglacialVoidRatio"; case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput"; case HydrologydcEplThicknessStackedEnum : return "HydrologydcEplThicknessStacked"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23936) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23937) @@ -565,4 +565,11 @@ else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum; else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum; + else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; + else if (strcmp(name,"HydrologyPressureMeltCoefficient")==0) return HydrologyPressureMeltCoefficientEnum; + else if (strcmp(name,"HydrologyCavitySpacing")==0) return HydrologyCavitySpacingEnum; + else if (strcmp(name,"HydrologyIschannels")==0) return HydrologyIschannelsEnum; + else if (strcmp(name,"HydrologyChannelConductivity")==0) return HydrologyChannelConductivityEnum; + else if (strcmp(name,"HydrologyChannelSheetWidth")==0) return HydrologyChannelSheetWidthEnum; + else if (strcmp(name,"HydrologyEnglacialVoidRatio")==0) return HydrologyEnglacialVoidRatioEnum; else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum; else if (strcmp(name,"HydrologydcEplThicknessStacked")==0) return HydrologydcEplThicknessStackedEnum; @@ -622,5 +629,8 @@ else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; - else if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum; else if (strcmp(name,"SealevelriseCumDeltathickness")==0) return SealevelriseCumDeltathicknessEnum; else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; @@ -629,8 +639,5 @@ else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; + else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum; @@ -745,5 +752,8 @@ else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; - else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; @@ -752,8 +762,5 @@ else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; + else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; - else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; @@ -875,8 +885,5 @@ else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; + else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; @@ -991,5 +998,8 @@ else if (strcmp(name,"FileParam")==0) return FileParamEnum; else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; - else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum; else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; + else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; @@ -1114,5 +1121,8 @@ else if (strcmp(name,"MinVel")==0) return MinVelEnum; else if (strcmp(name,"MinVx")==0) return MinVxEnum; - else if (strcmp(name,"MinVy")==0) return MinVyEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; @@ -1121,8 +1131,5 @@ else if (strcmp(name,"Mpi")==0) return MpiEnum; else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Mumps")==0) return MumpsEnum; + else if (strcmp(name,"Mumps")==0) return MumpsEnum; else if (strcmp(name,"Nodal")==0) return NodalEnum; else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; @@ -1237,5 +1244,8 @@ else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; - else if (strcmp(name,"Tria")==0) return TriaEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Tria")==0) return TriaEnum; else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; @@ -1244,8 +1254,5 @@ else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"VertexLId")==0) return VertexLIdEnum; + else if (strcmp(name,"VertexLId")==0) return VertexLIdEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; Index: /issm/trunk-jpl/src/m/classes/hydrologyglads.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyglads.m (revision 23936) +++ /issm/trunk-jpl/src/m/classes/hydrologyglads.m (revision 23937) @@ -7,11 +7,19 @@ properties (SetAccess=public) %Sheet - pressure_melt_coefficient = 0.; %c_t + pressure_melt_coefficient = 0.; + sheet_conductivity = NaN; + cavity_spacing = 0.; + bump_height = 0.; %Channels - ischannels = 0; + ischannels = 0; + channel_conductivity = 0.; + channel_sheet_width = 0.; %Other - requested_outputs = {}; + spcphi = NaN; + englacial_void_ratio = 0.; + moulin_input = NaN; + requested_outputs = {}; end methods @@ -34,9 +42,14 @@ %Sheet parameters self.pressure_melt_coefficient = 7.5e-8; %K/Pa (See table 1 in Erder et al. 2013) + self.cavity_spacing = 2.; %m + self.bump_height = .1; %m %Channel parameters self.ischannels=false; + self.channel_conductivity = 5.e-2; %Dow's default, Table uses 0.1 + self.channel_sheet_width = 2.; %m %Other + self.englacial_void_ratio = 1.e-5;% Dow's default, Table from Werder et al. uses 1e-3; self.requested_outputs={'default'}; end % }}} @@ -50,9 +63,17 @@ %Sheet md = checkfield(md,'fieldname','hydrology.pressure_melt_coefficient','numel',[1],'>=',0); + md = checkfield(md,'fieldname','hydrology.sheet_conductivity','size',[md.mesh.numberofvertices 1],'>',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','hydrology.cavity_spacing','numel',[1],'>',0); + md = checkfield(md,'fieldname','hydrology.bump_height','numel',[1],'>',0); %Channels md = checkfield(md,'fieldname','hydrology.ischannels','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','hydrology.channel_conductivity','numel',[1],'>',0); + md = checkfield(md,'fieldname','hydrology.channel_sheet_width','numel',[1],'>=',0); %Other + md = checkfield(md,'fieldname','hydrology.spcphi','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','hydrology.englacial_void_ratio','numel',[1],'>=',0); + md = checkfield(md,'fieldname','hydrology.moulin_input','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); end % }}} @@ -61,7 +82,14 @@ disp(sprintf(' SHEET')); fielddisplay(self,'pressure_melt_coefficient','Pressure melt coefficient (c_t) [K Pa^-1]'); + fielddisplay(self,'sheet_conductivity','sheet conductivity (k) [m^(7/4) kg^(-1/2)]'); + fielddisplay(self,'cavity_spacing','cavity spacing (l_r) [m]'); + fielddisplay(self,'bump_height','typical bump height (h_r) [m]'); disp(sprintf(' CHANNELS')); fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no'); + fielddisplay(self,'channel_conductivity','channel conductivity (k_c) [m^(3/2) kg^(-1/2)]'); disp(sprintf(' OTHER')); + fielddisplay(self,'spcphi','Hydraulic potential Dirichlet constraints [Pa]'); + fielddisplay(self,'englacial_void_ratio','englacial void ratio (e_v)'); + fielddisplay(self,'moulin_input','moulin input (Q_s) [m^3/s]'); fielddisplay(self,'requested_outputs','additional outputs requested'); end % }}} @@ -75,9 +103,17 @@ %Sheet WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','pressure_melt_coefficient','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_conductivity','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','cavity_spacing','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_height','format','Double'); %Channels WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','ischannels','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_conductivity','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_sheet_width','format','Double'); %Others + 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); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_void_ratio','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','moulin_input','format','DoubleMat','mattype',1); outputs = self.requested_outputs; pos = find(ismember(outputs,'default')); @@ -90,3 +126,2 @@ end end -