Index: /issm/trunk-jpl/m4/analyses.m4 =================================================================== --- /issm/trunk-jpl/m4/analyses.m4 (revision 19719) +++ /issm/trunk-jpl/m4/analyses.m4 (revision 19720) @@ -478,4 +478,30 @@ dnl }}} +dnl with-HydrologySommers{{{ + +AC_ARG_WITH([HydrologySommers], + + AS_HELP_STRING([--with-HydrologySommers = YES], [compile with HydrologySommers capabilities (default is yes)]), + + [HYDROLOGYSOMMERS=$withval],[HYDROLOGYSOMMERS=yes]) + +AC_MSG_CHECKING(for HydrologySommers capability compilation) + + +HAVE_HYDROLOGYSOMMERS=no + +if test "x$HYDROLOGYSOMMERS" = "xyes"; then + + HAVE_HYDROLOGYSOMMERS=yes + + AC_DEFINE([_HAVE_HYDROLOGYSOMMERS_],[1],[with HydrologySommers capability]) + +fi + +AM_CONDITIONAL([HYDROLOGYSOMMERS], [test x$HAVE_HYDROLOGYSOMMERS = xyes]) + +AC_MSG_RESULT($HAVE_HYDROLOGYSOMMERS) + +dnl }}} dnl with-Melting{{{ Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 19719) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 19720) @@ -391,4 +391,7 @@ issm_sources += ./analyses/HydrologyShreveAnalysis.cpp endif +if HYDROLOGYSOMMERS +issm_sources += ./analyses/HydrologySommersAnalysis.cpp +endif if HYDROLOGYDCINEFFICIENT issm_sources += ./analyses/HydrologyDCInefficientAnalysis.cpp Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 19719) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 19720) @@ -68,4 +68,7 @@ case HydrologyDCEfficientAnalysisEnum : return new HydrologyDCEfficientAnalysis(); #endif + #ifdef _HAVE_HYDROLOGYSOMMERS_ + case HydrologySommersAnalysisEnum : return new HydrologySommersAnalysis(); + #endif #ifdef _HAVE_MELTING_ case MeltingAnalysisEnum : return new MeltingAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 19720) +++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 19720) @@ -0,0 +1,232 @@ +#include "./HydrologySommersAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void HydrologySommersAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + if(hydrology_model!=HydrologysommersEnum) return; + + IoModelToConstraintsx(constraints,iomodel,HydrologySpcheadEnum,HydrologySommersAnalysisEnum,P1Enum); + +}/*}}}*/ +void HydrologySommersAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ +void HydrologySommersAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want Sommers?*/ + if(hydrology_model!=HydrologysommersEnum) return; + + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(nodes,iomodel,HydrologySommersAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +int HydrologySommersAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void HydrologySommersAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want Sommers?*/ + if(hydrology_model!=HydrologysommersEnum) return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,BaseEnum); + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum); + iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + } + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum); + iomodel->FetchDataToInput(elements,HydrologyHeadEnum); + iomodel->FetchDataToInput(elements,HydrologyGapHeightEnum); + iomodel->FetchDataToInput(elements,HydrologyEnglacialInputEnum); + iomodel->FetchDataToInput(elements,HydrologyBumpSpacingEnum); + iomodel->FetchDataToInput(elements,HydrologyReynoldsEnum); +}/*}}}*/ +void HydrologySommersAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want Sommers?*/ + if(hydrology_model!=HydrologysommersEnum) return; + + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); +}/*}}}*/ + +/*Finite Element Analysis*/ +void HydrologySommersAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* HydrologySommersAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* HydrologySommersAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* HydrologySommersAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble conductivity; + IssmDouble Jdet; + IssmDouble gap,reynolds; + IssmDouble* xyz_list = NULL; + + /*Hard coded parameters*/ + IssmDouble omega = 0.001; // parameter controlling transition to nonlinear resistance in basal system (dimensionless) + IssmDouble nu = 1.787e-6; // kinematic water viscosity m^2/s + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* dbasis = xNew(2*numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + Input* reynolds_input = element->GetInput(HydrologyReynoldsEnum); _assert_(reynolds_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Compute conductivity*/ + reynolds_input->GetInputValue(&reynolds,gauss); + gap_input->GetInputValue(&gap,gauss); + conductivity = pow(gap,3)*g/(12.*nu*(1+omega*reynolds)); + + for(int i=0;ivalues[i*numnodes+j] += -conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); + } + } + + + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(dbasis); + delete gauss; + return Ke; +}/*}}}*/ +ElementVector* HydrologySommersAnalysis::CreatePVector(Element* element){/*{{{*/ + _error_("STOP"); + + /*Skip if water or ice shelf element*/ + if(element->IsFloating()) return NULL; + + /*Intermediaries */ + IssmDouble Jdet,dt; + IssmDouble mb,oldw; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input); + Input* oldw_input = element->GetInput(WaterColumnOldEnum); _assert_(oldw_input); + + /*Initialize mb_correction to 0, do not forget!:*/ + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + mb_input->GetInputValue(&mb,gauss); + oldw_input->GetInputValue(&oldw,gauss); + + if(dt!=0.){ + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i]; + } + else{ + for(int i=0;ivalues[i]+=Jdet*gauss->weight*mb*basis[i]; + } + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete gauss; + return pe; +}/*}}}*/ +void HydrologySommersAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum); +}/*}}}*/ +void HydrologySommersAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void HydrologySommersAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + /*Intermediary*/ + int* doflist = NULL; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Fetch dof list and allocate solution vector*/ + element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numnodes); + + /*Use the dof list to index into the solution vector: */ + for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + } + + /*Add input to the element: */ + element->AddInput(HydrologyHeadEnum,values,element->GetElementType()); + + /*Free ressources:*/ + xDelete(values); + xDelete(doflist); +}/*}}}*/ +void HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h (revision 19720) +++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h (revision 19720) @@ -0,0 +1,33 @@ +/*! \file HydrologySommersAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _HydrologySommersAnalysis_ +#define _HydrologySommersAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class HydrologySommersAnalysis: public Analysis{ + + public: + /*Model processing*/ + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 19719) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 19720) @@ -27,4 +27,5 @@ #include "./HydrologyDCInefficientAnalysis.h" #include "./HydrologyShreveAnalysis.h" +#include "./HydrologySommersAnalysis.h" #include "./LevelsetAnalysis.h" #include "./LsfReinitializationAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19719) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19720) @@ -586,4 +586,5 @@ if(ishydrology){ analyses_temp[numanalyses++]=HydrologyShreveAnalysisEnum; + analyses_temp[numanalyses++]=HydrologySommersAnalysisEnum; analyses_temp[numanalyses++]=HydrologyDCInefficientAnalysisEnum; analyses_temp[numanalyses++]=HydrologyDCEfficientAnalysisEnum; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 19719) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 19720) @@ -146,4 +146,7 @@ } else if(hydrology_model==HydrologyshreveEnum){ + /*Nothing to add*/ + } + else if(hydrology_model==HydrologysommersEnum){ /*Nothing to add*/ } Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 19719) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 19720) @@ -79,4 +79,24 @@ } } + + else if (hydrology_model==HydrologysommersEnum){ + if(VerboseSolution()) _printf0_(" computing water head\n"); + femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); + solutionsequence_nonlinear(femmodel,modify_loads); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results \n"); + int outputs[2] = {HydrologyHeadEnum,HydrologyGapHeightEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + + /*unload results*/ + if(VerboseSolution()) _printf0_(" saving temporary results\n"); + OutputResultsx(femmodel); + } + } + + else{ + _error_("Hydrology model "<< EnumToStringx(hydrology_model) <<" not supported yet"); + } } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19719) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19720) @@ -156,4 +156,12 @@ HydrologyEfficientEnum, HydrologySedimentKmaxEnum, + HydrologysommersEnum, + HydrologyHeadEnum, + HydrologyGapHeightEnum, + HydrologyBumpSpacingEnum, + HydrologyBumpHeightEnum, + HydrologyEnglacialInputEnum, + HydrologyReynoldsEnum, + HydrologySpcheadEnum, IndependentObjectEnum, InversionControlParametersEnum, @@ -481,4 +489,5 @@ HydrologyDCInefficientAnalysisEnum, HydrologyDCEfficientAnalysisEnum, + HydrologySommersAnalysisEnum, HydrologySolutionEnum, MeltingAnalysisEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19719) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19720) @@ -162,4 +162,12 @@ case HydrologyEfficientEnum : return "HydrologyEfficient"; case HydrologySedimentKmaxEnum : return "HydrologySedimentKmax"; + case HydrologysommersEnum : return "Hydrologysommers"; + case HydrologyHeadEnum : return "HydrologyHead"; + case HydrologyGapHeightEnum : return "HydrologyGapHeight"; + case HydrologyBumpSpacingEnum : return "HydrologyBumpSpacing"; + case HydrologyBumpHeightEnum : return "HydrologyBumpHeight"; + case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput"; + case HydrologyReynoldsEnum : return "HydrologyReynolds"; + case HydrologySpcheadEnum : return "HydrologySpchead"; case IndependentObjectEnum : return "IndependentObject"; case InversionControlParametersEnum : return "InversionControlParameters"; @@ -479,4 +487,5 @@ case HydrologyDCInefficientAnalysisEnum : return "HydrologyDCInefficientAnalysis"; case HydrologyDCEfficientAnalysisEnum : return "HydrologyDCEfficientAnalysis"; + case HydrologySommersAnalysisEnum : return "HydrologySommersAnalysis"; case HydrologySolutionEnum : return "HydrologySolution"; case MeltingAnalysisEnum : return "MeltingAnalysis"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19719) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19720) @@ -165,4 +165,12 @@ else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum; else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum; + else if (strcmp(name,"Hydrologysommers")==0) return HydrologysommersEnum; + else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum; + else if (strcmp(name,"HydrologyGapHeight")==0) return HydrologyGapHeightEnum; + else if (strcmp(name,"HydrologyBumpSpacing")==0) return HydrologyBumpSpacingEnum; + else if (strcmp(name,"HydrologyBumpHeight")==0) return HydrologyBumpHeightEnum; + else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum; + else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; + else if (strcmp(name,"HydrologySpchead")==0) return HydrologySpcheadEnum; else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum; else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; @@ -252,5 +260,8 @@ else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum; else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; - else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum; else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum; else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum; @@ -260,8 +271,5 @@ else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum; else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"MeshElements2d")==0) return MeshElements2dEnum; + else if (strcmp(name,"MeshElements2d")==0) return MeshElements2dEnum; else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; else if (strcmp(name,"MeshLowerelements")==0) return MeshLowerelementsEnum; @@ -375,5 +383,8 @@ else if (strcmp(name,"SmbP")==0) return SmbPEnum; else if (strcmp(name,"SmbSwf")==0) return SmbSwfEnum; - else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum; else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum; else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; @@ -383,8 +394,5 @@ else if (strcmp(name,"SmbDt")==0) return SmbDtEnum; else if (strcmp(name,"SmbDz")==0) return SmbDzEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; + else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum; @@ -488,4 +496,5 @@ else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; + else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum; else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; @@ -497,5 +506,8 @@ else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; - else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; @@ -506,8 +518,5 @@ else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum; + else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum; else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum; else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum; @@ -620,5 +629,8 @@ else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; - else if (strcmp(name,"Adjointz")==0) return AdjointzEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Adjointz")==0) return AdjointzEnum; else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; @@ -629,8 +641,5 @@ else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; else if (strcmp(name,"Friction")==0) return FrictionEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"Internal")==0) return InternalEnum; + else if (strcmp(name,"Internal")==0) return InternalEnum; else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; @@ -743,5 +752,8 @@ else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; - else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; @@ -752,8 +764,5 @@ else if (strcmp(name,"Outputdefinition1")==0) return Outputdefinition1Enum; else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"Outputdefinition3")==0) return Outputdefinition3Enum; + else if (strcmp(name,"Outputdefinition3")==0) return Outputdefinition3Enum; else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; @@ -866,5 +875,8 @@ else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum; else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum; - else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum; @@ -875,8 +887,5 @@ else if (strcmp(name,"MinVel")==0) return MinVelEnum; else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"MinVx")==0) return MinVxEnum; + else if (strcmp(name,"MinVx")==0) return MinVxEnum; else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; Index: /issm/trunk-jpl/src/m/classes/hydrologysommers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologysommers.m (revision 19720) +++ /issm/trunk-jpl/src/m/classes/hydrologysommers.m (revision 19720) @@ -0,0 +1,71 @@ +%HYDROLOGYSOMMERS class definition +% +% Usage: +% hydrologysommers=hydrologysommers(); + +classdef hydrologysommers + properties (SetAccess=public) + head = NaN; + gap_height = NaN; + bump_spacing = NaN; + bump_height = NaN; + englacial_input = NaN; + reynolds = NaN; + spchead = NaN; + end + methods + function self = extrude(self,md) % {{{ + end % }}} + function self = hydrologysommers(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + self=structtoobj(self,varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + %Early return + if ~ismember(HydrologySommersAnalysisEnum(),analyses) + return; + end + + md = checkfield(md,'fieldname','hydrology.head','size',[md.mesh.numberofvertices 1],'NaN',1); + md = checkfield(md,'fieldname','hydrology.gap_height','>=',0,'size',[md.mesh.numberofelements 1],'NaN',1); + md = checkfield(md,'fieldname','hydrology.bump_spacing','>',0,'size',[md.mesh.numberofelements 1],'NaN',1); + md = checkfield(md,'fieldname','hydrology.bump_height','>=',0,'size',[md.mesh.numberofelements 1],'NaN',1); + md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'size',[md.mesh.numberofvertices 1],'NaN',1); + md = checkfield(md,'fieldname','hydrology.reynolds','>',0,'size',[md.mesh.numberofelements 1],'NaN',1); + md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]); + end % }}} + function disp(self) % {{{ + disp(sprintf(' hydrologysommers solution parameters:')); + fielddisplay(self,'head','subglacial hydrology water head (m)'); + fielddisplay(self,'gap_height','height of gap separating ice to bed (m)'); + fielddisplay(self,'bump_spacing','characteristic bedrock bump spacing (m)'); + fielddisplay(self,'bump_height','characteristic bedrock bump height (m)'); + fielddisplay(self,'englacial_input','liquid water input from englacial to subglacial system (m/yr)'); + fielddisplay(self,'reynolds','Reynolds'' number'); + fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)'); + end % }}} + function marshall(self,md,fid) % {{{ + + yts=365.0*24.0*3600.0; + + WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologysommersEnum(),'format','Integer'); + WriteData(fid,'object',self,'class','hydrology','fieldname','head','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','hydrology','fieldname','gap_height','format','DoubleMat','mattype',2); + WriteData(fid,'object',self,'class','hydrology','fieldname','bump_spacing','format','DoubleMat','mattype',2); + WriteData(fid,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',2); + WriteData(fid,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts); + WriteData(fid,'object',self,'class','hydrology','fieldname','reynolds','format','DoubleMat','mattype',2); + WriteData(fid,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1); + end % }}} + end +end + Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 19719) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 19720) @@ -182,4 +182,9 @@ md.friction.C=project2d(md,md.friction.C,1); md.friction.m=project2d(md,md.friction.m,1); + elseif isa(md.friction,'frictionweertmantemp'), + md.friction.C=project2d(md,md.friction.C,1); + md.friction.m=project2d(md,md.friction.m,1); + else + disp('friction type not supported'); end Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19719) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19720) @@ -154,4 +154,12 @@ def HydrologyEfficientEnum(): return StringToEnum("HydrologyEfficient")[0] def HydrologySedimentKmaxEnum(): return StringToEnum("HydrologySedimentKmax")[0] +def HydrologysommersEnum(): return StringToEnum("Hydrologysommers")[0] +def HydrologyHeadEnum(): return StringToEnum("HydrologyHead")[0] +def HydrologyGapHeightEnum(): return StringToEnum("HydrologyGapHeight")[0] +def HydrologyBumpSpacingEnum(): return StringToEnum("HydrologyBumpSpacing")[0] +def HydrologyBumpHeightEnum(): return StringToEnum("HydrologyBumpHeight")[0] +def HydrologyEnglacialInputEnum(): return StringToEnum("HydrologyEnglacialInput")[0] +def HydrologyReynoldsEnum(): return StringToEnum("HydrologyReynolds")[0] +def HydrologySpcheadEnum(): return StringToEnum("HydrologySpchead")[0] def IndependentObjectEnum(): return StringToEnum("IndependentObject")[0] def InversionControlParametersEnum(): return StringToEnum("InversionControlParameters")[0] @@ -471,4 +479,5 @@ def HydrologyDCInefficientAnalysisEnum(): return StringToEnum("HydrologyDCInefficientAnalysis")[0] def HydrologyDCEfficientAnalysisEnum(): return StringToEnum("HydrologyDCEfficientAnalysis")[0] +def HydrologySommersAnalysisEnum(): return StringToEnum("HydrologySommersAnalysis")[0] def HydrologySolutionEnum(): return StringToEnum("HydrologySolution")[0] def MeltingAnalysisEnum(): return StringToEnum("MeltingAnalysis")[0] Index: /issm/trunk-jpl/src/m/enum/HydrologyBumpHeightEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologyBumpHeightEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologyBumpHeightEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologyBumpHeightEnum() +%HYDROLOGYBUMPHEIGHTENUM - Enum of HydrologyBumpHeight +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologyBumpHeightEnum() + +macro=StringToEnum('HydrologyBumpHeight'); Index: /issm/trunk-jpl/src/m/enum/HydrologyBumpSpacingEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologyBumpSpacingEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologyBumpSpacingEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologyBumpSpacingEnum() +%HYDROLOGYBUMPSPACINGENUM - Enum of HydrologyBumpSpacing +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologyBumpSpacingEnum() + +macro=StringToEnum('HydrologyBumpSpacing'); Index: /issm/trunk-jpl/src/m/enum/HydrologyEnglacialInputEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologyEnglacialInputEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologyEnglacialInputEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologyEnglacialInputEnum() +%HYDROLOGYENGLACIALINPUTENUM - Enum of HydrologyEnglacialInput +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologyEnglacialInputEnum() + +macro=StringToEnum('HydrologyEnglacialInput'); Index: /issm/trunk-jpl/src/m/enum/HydrologyGapHeightEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologyGapHeightEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologyGapHeightEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologyGapHeightEnum() +%HYDROLOGYGAPHEIGHTENUM - Enum of HydrologyGapHeight +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologyGapHeightEnum() + +macro=StringToEnum('HydrologyGapHeight'); Index: /issm/trunk-jpl/src/m/enum/HydrologyHeadEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologyHeadEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologyHeadEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologyHeadEnum() +%HYDROLOGYHEADENUM - Enum of HydrologyHead +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologyHeadEnum() + +macro=StringToEnum('HydrologyHead'); Index: /issm/trunk-jpl/src/m/enum/HydrologyReynoldsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologyReynoldsEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologyReynoldsEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologyReynoldsEnum() +%HYDROLOGYREYNOLDSENUM - Enum of HydrologyReynolds +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologyReynoldsEnum() + +macro=StringToEnum('HydrologyReynolds'); Index: /issm/trunk-jpl/src/m/enum/HydrologySommersAnalysisEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologySommersAnalysisEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologySommersAnalysisEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologySommersAnalysisEnum() +%HYDROLOGYSOMMERSANALYSISENUM - Enum of HydrologySommersAnalysis +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologySommersAnalysisEnum() + +macro=StringToEnum('HydrologySommersAnalysis'); Index: /issm/trunk-jpl/src/m/enum/HydrologySpcheadEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologySpcheadEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologySpcheadEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologySpcheadEnum() +%HYDROLOGYSPCHEADENUM - Enum of HydrologySpchead +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologySpcheadEnum() + +macro=StringToEnum('HydrologySpchead'); Index: /issm/trunk-jpl/src/m/enum/HydrologysommersEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologysommersEnum.m (revision 19720) +++ /issm/trunk-jpl/src/m/enum/HydrologysommersEnum.m (revision 19720) @@ -0,0 +1,11 @@ +function macro=HydrologysommersEnum() +%HYDROLOGYSOMMERSENUM - Enum of Hydrologysommers +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologysommersEnum() + +macro=StringToEnum('Hydrologysommers');