Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 14756) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 14757) @@ -90,4 +90,5 @@ HydrologydcEnum, SedimentHeadEnum, + EplHeadEnum, HydrologydcSpcsedimentHeadEnum, HydrologydcSedimentCompressibilityEnum, Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14756) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14757) @@ -5799,6 +5799,8 @@ ElementMatrix* Tria::CreateKMatrixHydrology(void){ - int hydrology_model; + int hydrology_model; + bool isefficientlayer; parameters->FindParam(&hydrology_model,HydrologyEnum); + switch(hydrology_model){ @@ -5806,5 +5808,9 @@ return CreateKMatrixHydrologyShreve(); case HydrologydcEnum: - return CreateKMatrixHydrologyDC(); + parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(isefficientlayer) + return CreateKMatrixHydrologyDCsediment(); + else + return CreateKMatrixHydrologyDCepl(); default: _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet"); @@ -5917,6 +5923,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrologyDC{{{*/ -ElementMatrix* Tria::CreateKMatrixHydrologyDC(void){ +/*FUNCTION Tria::CreateKMatrixHydrologyDCsediment{{{*/ +ElementMatrix* Tria::CreateKMatrixHydrologyDCsediment(void){ /*constants: */ @@ -5978,8 +5984,70 @@ } /*}}}*/ +/*FUNCTION Tria::CreateKMatrixHydrologyDCepl{{{*/ +ElementMatrix* Tria::CreateKMatrixHydrologyDCepl(void){ + + /*constants: */ + const int numdof=NDOF1*NUMVERTICES; + + /* Intermediaries */ + IssmDouble D_scalar,Jdet; + IssmDouble epl_transmitivity,dt; + IssmDouble epl_storing; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[2][numdof]; + IssmDouble L[NUMVERTICES]; + IssmDouble D[2][2]; + GaussTria *gauss = NULL; + + /*Initialize Element matrix*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + epl_storing = matpar->GetEplStoring(); + epl_transmitivity = matpar->GetEplTransmitivity(); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Diffusivity*/ + D_scalar=epl_transmitivity*gauss->weight*Jdet; + if(reCast(dt)) D_scalar=D_scalar*dt; + D[0][0]=D_scalar; D[0][1]=0.; + D[1][0]=0.; D[1][1]=D_scalar; + GetBHydro(&B[0][0],&xyz_list[0][0],gauss); + TripleMultiply(&B[0][0],2,numdof,1, + &D[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke->values[0],1); + + /*Transient*/ + if(reCast(dt)){ + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + D_scalar=epl_storing*gauss->weight*Jdet; + + TripleMultiply(&L[0],numdof,1,0, + &D_scalar,1,1,0, + &L[0],1,numdof,0, + &Ke->values[0],1); + } + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreatePVectorHydrology{{{*/ ElementVector* Tria::CreatePVectorHydrology(void){ - int hydrology_model; + int hydrology_model; + bool isefficientlayer; parameters->FindParam(&hydrology_model,HydrologyEnum); @@ -5988,5 +6056,9 @@ return CreatePVectorHydrologyShreve(); case HydrologydcEnum: - return CreatePVectorHydrologyDC(); + parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(isefficientlayer) + return CreatePVectorHydrologyDCsediment(); + else + return CreatePVectorHydrologyDCepl(); default: _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet"); @@ -6044,6 +6116,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrologyDC {{{*/ -ElementVector* Tria::CreatePVectorHydrologyDC(void){ +/*FUNCTION Tria::CreatePVectorHydrologyDCsediment {{{*/ +ElementVector* Tria::CreatePVectorHydrologyDCsediment(void){ /*Constants*/ @@ -6092,4 +6164,61 @@ old_wh_input->GetInputValue(&water_head,gauss); scalar = Jdet*gauss->weight*water_head*sediment_storing; + for(int i=0;ivalues[i]+=scalar*basis[i]; + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorHydrologyDCepl {{{*/ +ElementVector* Tria::CreatePVectorHydrologyDCepl(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dt,scalar,water_head; + IssmDouble water_load; + IssmDouble epl_storing; + IssmDouble basis[numdof]; + GaussTria* gauss=NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + epl_storing = matpar->GetEplStoring(); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input); + Input* old_wh_input=NULL; + + if(reCast(dt)){ + old_wh_input=inputs->GetInput(EplHeadEnum); _assert_(old_wh_input); + } + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(basis, gauss); + + /*Loading term*/ + water_input->GetInputValue(&water_load,gauss); + scalar = Jdet*gauss->weight*water_load; + if(reCast(dt)) scalar = scalar*dt; + for(int i=0;ivalues[i]+=scalar*basis[i]; + + /*Transient term*/ + if(reCast(dt)){ + old_wh_input->GetInputValue(&water_head,gauss); + scalar = Jdet*gauss->weight*water_head*epl_storing; for(int i=0;ivalues[i]+=scalar*basis[i]; } Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14756) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14757) @@ -240,8 +240,10 @@ ElementMatrix* CreateKMatrixHydrology(void); ElementMatrix* CreateKMatrixHydrologyShreve(void); - ElementMatrix* CreateKMatrixHydrologyDC(void); + ElementMatrix* CreateKMatrixHydrologyDCsediment(void); + ElementMatrix* CreateKMatrixHydrologyDCepl(void); ElementVector* CreatePVectorHydrology(void); ElementVector* CreatePVectorHydrologyShreve(void); - ElementVector* CreatePVectorHydrologyDC(void); + ElementVector* CreatePVectorHydrologyDCsediment(void); + ElementVector* CreatePVectorHydrologyDCepl(void); void CreateHydrologyWaterVelocityInput(void); void GetSolutionFromInputsHydrology(Vector* solution); Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp (revision 14756) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp (revision 14757) @@ -63,4 +63,9 @@ iomodel->Constant(&this->sediment_transmitivity,HydrologydcSedimentTransmitivityEnum); iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum); + + iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum); + iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum); + iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum); + iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum); } else{ @@ -355,7 +360,18 @@ } /*}}}*/ +/*FUNCTION Matpar::GetEplStoring {{{*/ +IssmDouble Matpar::GetEplStoring(){ + return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness* + ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity)); +} +/*}}}*/ /*FUNCTION Matpar::GetSedimentTransitivity {{{*/ IssmDouble Matpar::GetSedimentTransmitivity(){ return sediment_transmitivity; +} +/*}}}*/ +/*FUNCTION Matpar::GetEplTransitivity {{{*/ +IssmDouble Matpar::GetEplTransmitivity(){ + return epl_transmitivity; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h (revision 14756) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h (revision 14757) @@ -45,4 +45,9 @@ IssmDouble sediment_transmitivity; IssmDouble water_compressibility; + + IssmDouble epl_compressibility; + IssmDouble epl_porosity; + IssmDouble epl_thickness; + IssmDouble epl_transmitivity; /*gia: */ @@ -118,5 +123,7 @@ IssmDouble GetHydrologyN(); IssmDouble GetSedimentStoring(); + IssmDouble GetEplStoring(); IssmDouble GetSedimentTransmitivity(); + IssmDouble GetEplTransmitivity(); IssmDouble TMeltingPoint(IssmDouble pressure); IssmDouble PureIceEnthalpy(IssmDouble pressure); Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 14756) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 14757) @@ -95,4 +95,5 @@ case HydrologydcEnum : return "Hydrologydc"; case SedimentHeadEnum : return "SedimentHead"; + case EplHeadEnum : return "EplHead"; case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead"; case HydrologydcSedimentCompressibilityEnum : return "HydrologydcSedimentCompressibility"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp (revision 14756) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp (revision 14757) @@ -16,4 +16,5 @@ Parameters *parameters = NULL; int hydrology_model; + bool isefficientlayer; /*Get parameters: */ @@ -28,5 +29,6 @@ } else if(hydrology_model==HydrologydcEnum){ - //No parameters for now + iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum); + parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); } else{ Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 14756) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 14757) @@ -96,4 +96,5 @@ else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; + else if (strcmp(name,"EplHead")==0) return EplHeadEnum; else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum; else if (strcmp(name,"HydrologydcSedimentCompressibility")==0) return HydrologydcSedimentCompressibilityEnum; @@ -137,9 +138,9 @@ else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum; else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; - else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum; + if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum; + else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum; else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum; else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum; @@ -260,9 +261,9 @@ else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; - else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; + if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum; + else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum; else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; @@ -383,9 +384,9 @@ else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; else if (strcmp(name,"Adjointz")==0) return AdjointzEnum; - else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; + if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; + else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; else if (strcmp(name,"Boundary")==0) return BoundaryEnum; @@ -506,9 +507,9 @@ else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; else if (strcmp(name,"Incremental")==0) return IncrementalEnum; - else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"None")==0) return NoneEnum; + if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum; + else if (strcmp(name,"None")==0) return NoneEnum; else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; Index: /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp (revision 14756) +++ /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp (revision 14757) @@ -78,4 +78,5 @@ } } + else if (hydrology_model==HydrologydcEnum){ if(VerboseSolution()) _pprintLine_(" computing water head"); Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 14756) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 14757) @@ -6,13 +6,14 @@ classdef initialization properties (SetAccess=public) - vx = NaN; - vy = NaN; - vz = NaN; - vel = NaN; - pressure = NaN; - temperature = NaN; - waterfraction = NaN; - sediment_head = NaN; - watercolumn = NaN; + vx = NaN; + vy = NaN; + vz = NaN; + vel = NaN; + pressure = NaN; + temperature = NaN; + waterfraction = NaN; + sediment_head = NaN; + epl_head = NaN; + watercolumn = NaN; end methods @@ -59,4 +60,5 @@ if isa(md.hydrology,'hydrologydc'), md = checkfield(md,'initialization.sediment_head','NaN',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]); else md = checkfield(md,'initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices 1]); @@ -75,4 +77,5 @@ fielddisplay(obj,'waterfraction','fraction of water in the ice'); fielddisplay(obj,'sediment_head','sediment water head of subglacial system [m]'); + fielddisplay(obj,'epl_head','epl water head of subglacial system [m]'); fielddisplay(obj,'watercolumn','thickness of subglacial water [m]'); @@ -86,4 +89,5 @@ WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum); WriteData(fid,'data',obj.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum); + WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum); WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum); end % }}} Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 14756) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 14757) @@ -789,4 +789,14 @@ return StringToEnum('SedimentHead')[0] +def EplHeadEnum(): + """ + EPLHEADENUM - Enum of EplHead + + Usage: + macro=EplHeadEnum() + """ + + return StringToEnum('EplHead')[0] + def HydrologydcSpcsedimentHeadEnum(): """ @@ -5337,4 +5347,4 @@ """ - return 532 - + return 533 + Index: /issm/trunk-jpl/src/m/enum/EplHeadEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/EplHeadEnum.m (revision 14757) +++ /issm/trunk-jpl/src/m/enum/EplHeadEnum.m (revision 14757) @@ -0,0 +1,11 @@ +function macro=EplHeadEnum() +%EPLHEADENUM - Enum of EplHead +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=EplHeadEnum() + +macro=StringToEnum('EplHead'); Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 14756) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 14757) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=532; +macro=533;