Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 23019) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 23020) @@ -419,6 +419,9 @@ issm_sources += ./analyses/HydrologyShreveAnalysis.cpp endif -if HYDROLOGYSOMMERS -issm_sources += ./analyses/HydrologySommersAnalysis.cpp +if HYDROLOGYSHAKTI +issm_sources += ./analyses/HydrologyShaktiAnalysis.cpp +endif +if HYDROLOGYPISM +issm_sources += ./analyses/HydrologyPismAnalysis.cpp endif if HYDROLOGYDCINEFFICIENT Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 23020) @@ -77,6 +77,9 @@ case HydrologyShreveAnalysisEnum : return new HydrologyShreveAnalysis(); #endif - #ifdef _HAVE_HYDROLOGYSOMMERS_ - case HydrologySommersAnalysisEnum : return new HydrologySommersAnalysis(); + #ifdef _HAVE_HYDROLOGYSHAKTI_ + case HydrologyShaktiAnalysisEnum : return new HydrologyShaktiAnalysis(); + #endif + #ifdef _HAVE_HYDROLOGYPISM_ + case HydrologyPismAnalysisEnum : return new HydrologyPismAnalysis(); #endif #ifdef _HAVE_L2PROJECTIONBASE_ Index: /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23020) +++ /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23020) @@ -0,0 +1,119 @@ +#include "./HydrologyPismAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void HydrologyPismAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + return; + +}/*}}}*/ +void HydrologyPismAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + return; + +}/*}}}*/ +void HydrologyPismAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + return; +}/*}}}*/ +int HydrologyPismAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 0; +}/*}}}*/ +void HydrologyPismAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + int hydrology_model,frictionlaw; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Pism?*/ + if(hydrology_model!=HydrologypismEnum) return; + + /*Add input to elements*/ + iomodel->FetchDataToInput(elements,"md.hydrology.drainage_rate",HydrologyDrainageRateEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.watercolumn",WatercolumnEnum,0.); +}/*}}}*/ +void HydrologyPismAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + int numoutputs; + char** requestedoutputs = NULL; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Pism?*/ + if(hydrology_model!=HydrologypismEnum) return; + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + + /*Nothing else to add for now*/ +}/*}}}*/ + +/*Finite Element Analysis*/ +void HydrologyPismAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* HydrologyPismAnalysis::CreateDVector(Element* element){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementMatrix* HydrologyPismAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* HydrologyPismAnalysis::CreateKMatrix(Element* element){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* HydrologyPismAnalysis::CreatePVector(Element* element){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +void HydrologyPismAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +void HydrologyPismAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void HydrologyPismAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +void HydrologyPismAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ + +/*Additional methods*/ +void HydrologyPismAnalysis::UpdateWaterColumn(FemModel* femmodel){/*{{{*/ + + for(int j=0;jelements->Size();j++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + this->UpdateWaterColumn(element); + } + +}/*}}}*/ +void HydrologyPismAnalysis::UpdateWaterColumn(Element* element){/*{{{*/ + + /*Skip if water or ice shelf element*/ + if(element->IsFloating()) return; + + /*Intermediaries */ + IssmDouble dt,drainage_rate,water_column; + + /*Retrieve all inputs and parameters*/ + element->FindParam(&dt,TimesteppingTimeStepEnum); + + /*Get water column and drainage rate*/ + const int numvertices= element->GetNumberOfVertices(); + IssmDouble* watercolumn = xNew(numvertices); + IssmDouble* drainagerate = xNew(numvertices); + IssmDouble* meltingrate = xNew(numvertices); + element->GetInputListOnVertices(&watercolumn[0],WatercolumnEnum); + element->GetInputListOnVertices(&drainagerate[0],HydrologyDrainageRateEnum); + element->GetInputListOnVertices(&meltingrate[0],BasalforcingsGroundediceMeltingRateEnum); + + /*Add water*/ + for(int i=0;iAddInput(WatercolumnEnum,&watercolumn[0],P1Enum); + + /*Clean up and return*/ + xDelete(watercolumn); + xDelete(drainagerate); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.h (revision 23020) +++ /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.h (revision 23020) @@ -0,0 +1,37 @@ +/*! \file HydrologyPismAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _HydrologyPismAnalysis_ +#define _HydrologyPismAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class HydrologyPismAnalysis: 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); + + /*Intermediaries*/ + void UpdateWaterColumn(FemModel* femmodel); + void UpdateWaterColumn(Element* element); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23020) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23020) @@ -0,0 +1,596 @@ +#include "./HydrologyShaktiAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Define 2 hardcoded parameters*/ +#define OMEGA 0.001 // parameter controlling transition to nonlinear resistance in basal system (dimensionless) +#define NU 1.787e-6 //kinematic water viscosity m^2/s +#define CT 7.5e-8 // Clapeyron slope (K/Pa) +#define CW 4.22e3 // specific heat capacity of water (J/kg/K) + +/*Model processing*/ +void HydrologyShaktiAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + if(hydrology_model!=HydrologyshaktiEnum) return; + + IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologyShaktiAnalysisEnum,P1Enum); + +}/*}}}*/ +void HydrologyShaktiAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Shakti?*/ + if(hydrology_model!=HydrologyshaktiEnum) return; + + /*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(iomodel->loadcounter+i+1,i,iomodel,HydrologyShaktiAnalysisEnum)); + } + } + else if(reCast(iomodel->Data("md.mesh.vertexonbase")[i])){ + if(iomodel->my_vertices[i]){ + loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyShaktiAnalysisEnum)); + } + } + } + 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(iomodel->loadcounter+i+1,i,iomodel,segments,HydrologyShaktiAnalysisEnum)); + } + } + + xDelete(segments); + +}/*}}}*/ +void HydrologyShaktiAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Shakti?*/ + if(hydrology_model!=HydrologyshaktiEnum) return; + + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); + ::CreateNodes(nodes,iomodel,HydrologyShaktiAnalysisEnum,P1Enum); + iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); +}/*}}}*/ +int HydrologyShaktiAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void HydrologyShaktiAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + int hydrology_model,frictionlaw; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Shakti?*/ + if(hydrology_model!=HydrologyshaktiEnum) 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,"md.geometry.thickness",ThicknessEnum); + iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); + iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); + } + iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.head",HydrologyHeadEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.gap_height",HydrologyGapHeightEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.englacial_input",HydrologyEnglacialInputEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.bump_spacing",HydrologyBumpSpacingEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.reynolds",HydrologyReynoldsEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); + iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); + iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); + iomodel->FindConstant(&frictionlaw,"md.friction.law"); + + /*Friction law variables*/ + switch(frictionlaw){ + case 1: + iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); + iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); + break; + case 8: + iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); + break; + default: + _error_("Friction law "<< frictionlaw <<" not supported"); + } +}/*}}}*/ +void HydrologyShaktiAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + int numoutputs; + char** requestedoutputs = NULL; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Shakti?*/ + if(hydrology_model!=HydrologyshaktiEnum) return; + + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum)); + + /*Requested outputs*/ + iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs"); + parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs"); +}/*}}}*/ + +/*Finite Element Analysis*/ +void HydrologyShaktiAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* HydrologyShaktiAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* HydrologyShaktiAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* HydrologyShaktiAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble* xyz_list = NULL; + + /*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); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Get conductivity from inputs*/ + IssmDouble conductivity = GetConductivity(element); + + /*Get englacial storage coefficient*/ + IssmDouble storage,dt; + element->FindParam(&storage,HydrologyStorageEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(1); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + 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]) + + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]; + } + } + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(dbasis); + delete gauss; + return Ke; +}/*}}}*/ +ElementVector* HydrologyShaktiAnalysis::CreatePVector(Element* element){/*{{{*/ + + /*Skip if water or ice shelf element*/ + if(element->IsFloating()) return NULL; + + /*Intermediaries */ + IssmDouble Jdet,meltrate,G,dh[2],B,A,n; + IssmDouble gap,bed,thickness,head,ieb,head_old; + IssmDouble lr,br,vx,vy,beta; + IssmDouble alpha2,frictionheat; + IssmDouble PMPheat,dpressure_water[2],dbed[2]; + 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); + IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); + Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); + Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); + Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); + + /*Get conductivity from inputs*/ + IssmDouble conductivity = GetConductivity(element); + + /*Get englacial storage coefficient*/ + IssmDouble storage,dt; + element->FindParam(&storage,HydrologyStorageEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,2); + + /* 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); + geothermalflux_input->GetInputValue(&G,gauss); + base_input->GetInputValue(&bed,gauss); + base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss); + thickness_input->GetInputValue(&thickness,gauss); + gap_input->GetInputValue(&gap,gauss); + head_input->GetInputValue(&head,gauss); + head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); + englacial_input->GetInputValue(&ieb,gauss); + lr_input->GetInputValue(&lr,gauss); + br_input->GetInputValue(&br,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + headold_input->GetInputValue(&head_old,gauss); + + /*Get ice A parameter*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + A=pow(B,-n); + + /*Compute beta term*/ + if(gapGetAlpha2(&alpha2,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + frictionheat=alpha2*(vx*vx+vy*vy); + + /*Get water and ice pressures*/ + IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); + IssmDouble pressure_water = rho_water*g*(head-bed); + if(pressure_water>pressure_ice) pressure_water = pressure_ice; + + /*Get water pressure from previous time step to use in lagged creep term*/ + IssmDouble pressure_water_old = rho_water*g*(head_old-bed); + if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice; + + /*Compute change in sensible heat due to changes in pressure melting point*/ + dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); + dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); + PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]); + + meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); + _assert_(meltrate>0.); + + for(int i=0;ivalues[i]+=Jdet*gauss->weight* + ( + meltrate*(1/rho_water-1/rho_ice) + +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap + -beta*sqrt(vx*vx+vy*vy) + +ieb + +storage*head_old/dt + )*basis[i]; + } + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete friction; + delete gauss; + return pe; +}/*}}}*/ +void HydrologyShaktiAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum); +}/*}}}*/ +void HydrologyShaktiAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void HydrologyShaktiAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + /*Intermediary*/ + IssmDouble dh[3]; + int* doflist = NULL; + IssmDouble* xyz_list = NULL; + + /*Get gravity from parameters*/ + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + + /*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); + + /*Get thickness and base on nodes to apply cap on water head*/ + IssmDouble* eff_pressure = xNew(numnodes); + IssmDouble* thickness = xNew(numnodes); + IssmDouble* bed = xNew(numnodes); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + element->GetInputListOnNodes(&thickness[0],ThicknessEnum); + element->GetInputListOnNodes(&bed[0],BaseEnum); + + /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/ + IssmDouble* head_old = xNew(numnodes); + element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum); + IssmDouble relaxation; + element->FindParam(&relaxation,HydrologyRelaxationEnum); + + /*Use the dof list to index into the solution vector: */ + for(int i=0;i hrho_ice*thickness[i]/rho_water+bed[i]){ + values[i] = rho_ice*thickness[i]/rho_water+bed[i]; + } + + /*Make sure that negative pressure is not allowed*/ + // if(values[i](values[i])) _error_("NaN found in solution vector"); + if(xIsInf(values[i])) _error_("Inf found in solution vector"); + } + + /*Add input to the element: */ + element->AddInput(HydrologyHeadEnum,values,element->GetElementType()); + element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum); + + /*Update reynolds number according to new solution*/ + element->GetVerticesCoordinates(&xyz_list); + Input* head_input = element->GetInput(HydrologyHeadEnum);_assert_(head_input); + head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list); + IssmDouble conductivity = GetConductivity(element); + + IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU; + element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum); + + /*Free resources:*/ + xDelete(values); + xDelete(thickness); + xDelete(bed); + xDelete(xyz_list); + xDelete(doflist); + xDelete(eff_pressure); + xDelete(head_old); +}/*}}}*/ +void HydrologyShaktiAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ + +/*Additional methods*/ +IssmDouble HydrologyShaktiAnalysis::GetConductivity(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble gap,reynolds; + + /*Get gravity from parameters*/ + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + + /*Get Reynolds and gap average values*/ + Input* reynolds_input = element->GetInput(HydrologyReynoldsEnum); _assert_(reynolds_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + reynolds_input->GetInputAverage(&reynolds); + gap_input->GetInputAverage(&gap); + + /*Compute conductivity*/ + IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds)); + _assert_(conductivity>0); + + /*Clean up and return*/ + return conductivity; +}/*}}}*/ +void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/ + + + for(int j=0;jelements->Size();j++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + UpdateGapHeight(element); + } + +}/*}}}*/ +void HydrologyShaktiAnalysis::UpdateGapHeight(Element* element){/*{{{*/ + + /*Skip if water or ice shelf element*/ + if(element->IsFloating()) return; + + /*Intermediaries */ + IssmDouble newgap = 0.; + IssmDouble Jdet,meltrate,G,dh[2],B,A,n,dt; + IssmDouble gap,bed,thickness,head,ieb; + IssmDouble lr,br,vx,vy,beta; + IssmDouble alpha2,frictionheat; + IssmDouble* xyz_list = NULL; + IssmDouble dpressure_water[2],dbed[2],PMPheat; + IssmDouble q = 0.; + IssmDouble channelization = 0.; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); + Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); + Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); + + /*Get conductivity from inputs*/ + IssmDouble conductivity = GetConductivity(element); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,2); + + /*Keep track of weights*/ + IssmDouble totalweights=0.; + + /* 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); + + geothermalflux_input->GetInputValue(&G,gauss); + base_input->GetInputValue(&bed,gauss); + base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss); + thickness_input->GetInputValue(&thickness,gauss); + gap_input->GetInputValue(&gap,gauss); + head_input->GetInputValue(&head,gauss); + head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); + englacial_input->GetInputValue(&ieb,gauss); + lr_input->GetInputValue(&lr,gauss); + br_input->GetInputValue(&br,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + /*Get ice A parameter*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + A=pow(B,-n); + + /*Compute beta term*/ + if(gapGetAlpha2(&alpha2,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + frictionheat=alpha2*(vx*vx+vy*vy); + + /*Get water and ice pressures*/ + IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); + IssmDouble pressure_water = rho_water*g*(head-bed); + if(pressure_water>pressure_ice) pressure_water = pressure_ice; + + /* Compute change in sensible heat due to changes in pressure melting point*/ + dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); + dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); + PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]); + + meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); + _assert_(meltrate>0.); + + newgap += gauss->weight*Jdet*(gap+dt*( + meltrate/rho_ice + -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap + +beta*sqrt(vx*vx+vy*vy) + )); + totalweights +=gauss->weight*Jdet; + + /* Compute basal water flux */ + q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])); + + /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */ + channelization += gauss->weight*Jdet*(meltrate/rho_ice/(meltrate/rho_ice+beta*sqrt(vx*vx+vy*vy))); + } + + /*Divide by connectivity*/ + newgap = newgap/totalweights; + IssmDouble mingap = 1e-6; + if(newgapthickness) + newgap = thickness; + + /*Add new gap as an input*/ + element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum); + + /*Divide by connectivity, add basal flux as an input*/ + q = q/totalweights; + element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum); + + /* Divide by connectivity, add degree of channelization as an input */ + channelization = channelization/totalweights; + element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum); + + /*Clean up and return*/ + xDelete(xyz_list); + delete friction; + delete gauss; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.h (revision 23020) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.h (revision 23020) @@ -0,0 +1,38 @@ +/*! \file HydrologyShaktiAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _HydrologyShaktiAnalysis_ +#define _HydrologyShaktiAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class HydrologyShaktiAnalysis: 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); + + /*Intermediaries*/ + IssmDouble GetConductivity(Element* element); + void UpdateGapHeight(FemModel* femmodel); + void UpdateGapHeight(Element* element); +}; +#endif Index: sm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 23019) +++ (revision ) @@ -1,596 +1,0 @@ -#include "./HydrologySommersAnalysis.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" - -/*Define 2 hardcoded parameters*/ -#define OMEGA 0.001 // parameter controlling transition to nonlinear resistance in basal system (dimensionless) -#define NU 1.787e-6 //kinematic water viscosity m^2/s -#define CT 7.5e-8 // Clapeyron slope (K/Pa) -#define CW 4.22e3 // specific heat capacity of water (J/kg/K) - -/*Model processing*/ -void HydrologySommersAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - - /*retrieve some parameters: */ - int hydrology_model; - iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); - - if(hydrology_model!=HydrologysommersEnum) return; - - IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologySommersAnalysisEnum,P1Enum); - -}/*}}}*/ -void HydrologySommersAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - - /*Fetch parameters: */ - int hydrology_model; - iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); - - /*Now, do we really want Sommers?*/ - if(hydrology_model!=HydrologysommersEnum) return; - - /*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(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum)); - } - } - else if(reCast(iomodel->Data("md.mesh.vertexonbase")[i])){ - if(iomodel->my_vertices[i]){ - loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum)); - } - } - } - 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(iomodel->loadcounter+i+1,i,iomodel,segments,HydrologySommersAnalysisEnum)); - } - } - - xDelete(segments); - -}/*}}}*/ -void HydrologySommersAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - /*Fetch parameters: */ - int hydrology_model; - iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); - - /*Now, do we really want Sommers?*/ - if(hydrology_model!=HydrologysommersEnum) return; - - if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); - ::CreateNodes(nodes,iomodel,HydrologySommersAnalysisEnum,P1Enum); - iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); -}/*}}}*/ -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,frictionlaw; - iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); - - /*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,"md.geometry.thickness",ThicknessEnum); - iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); - if(iomodel->domaintype!=Domain2DhorizontalEnum){ - iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); - iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); - } - iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); - iomodel->FetchDataToInput(elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.head",HydrologyHeadEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.gap_height",HydrologyGapHeightEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.englacial_input",HydrologyEnglacialInputEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.bump_spacing",HydrologyBumpSpacingEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.reynolds",HydrologyReynoldsEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); - iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); - iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); - iomodel->FindConstant(&frictionlaw,"md.friction.law"); - - /*Friction law variables*/ - switch(frictionlaw){ - case 1: - iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); - iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); - iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - break; - case 8: - iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); - break; - default: - _error_("Friction law "<< frictionlaw <<" not supported"); - } -}/*}}}*/ -void HydrologySommersAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ - - /*retrieve some parameters: */ - int hydrology_model; - int numoutputs; - char** requestedoutputs = NULL; - iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); - - /*Now, do we really want Sommers?*/ - if(hydrology_model!=HydrologysommersEnum) return; - - parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); - parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum)); - - /*Requested outputs*/ - iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs"); - parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs)); - if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs)); - iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs"); -}/*}}}*/ - -/*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 Jdet; - IssmDouble* xyz_list = NULL; - - /*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); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - - /*Get conductivity from inputs*/ - IssmDouble conductivity = GetConductivity(element); - - /*Get englacial storage coefficient*/ - IssmDouble storage,dt; - element->FindParam(&storage,HydrologyStorageEnum); - element->FindParam(&dt,TimesteppingTimeStepEnum); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(1); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - element->NodalFunctions(basis,gauss); - - 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]) - + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]; - } - } - } - - /*Clean up and return*/ - xDelete(xyz_list); - xDelete(basis); - xDelete(dbasis); - delete gauss; - return Ke; -}/*}}}*/ -ElementVector* HydrologySommersAnalysis::CreatePVector(Element* element){/*{{{*/ - - /*Skip if water or ice shelf element*/ - if(element->IsFloating()) return NULL; - - /*Intermediaries */ - IssmDouble Jdet,meltrate,G,dh[2],B,A,n; - IssmDouble gap,bed,thickness,head,ieb,head_old; - IssmDouble lr,br,vx,vy,beta; - IssmDouble alpha2,frictionheat; - IssmDouble PMPheat,dpressure_water[2],dbed[2]; - 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); - IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); - IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); - IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); - Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); - Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); - Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); - Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); - - /*Get conductivity from inputs*/ - IssmDouble conductivity = GetConductivity(element); - - /*Get englacial storage coefficient*/ - IssmDouble storage,dt; - element->FindParam(&storage,HydrologyStorageEnum); - element->FindParam(&dt,TimesteppingTimeStepEnum); - - /*Build friction element, needed later: */ - Friction* friction=new Friction(element,2); - - /* 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); - geothermalflux_input->GetInputValue(&G,gauss); - base_input->GetInputValue(&bed,gauss); - base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss); - thickness_input->GetInputValue(&thickness,gauss); - gap_input->GetInputValue(&gap,gauss); - head_input->GetInputValue(&head,gauss); - head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); - englacial_input->GetInputValue(&ieb,gauss); - lr_input->GetInputValue(&lr,gauss); - br_input->GetInputValue(&br,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - headold_input->GetInputValue(&head_old,gauss); - - /*Get ice A parameter*/ - B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); - A=pow(B,-n); - - /*Compute beta term*/ - if(gapGetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - frictionheat=alpha2*(vx*vx+vy*vy); - - /*Get water and ice pressures*/ - IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); - IssmDouble pressure_water = rho_water*g*(head-bed); - if(pressure_water>pressure_ice) pressure_water = pressure_ice; - - /*Get water pressure from previous time step to use in lagged creep term*/ - IssmDouble pressure_water_old = rho_water*g*(head_old-bed); - if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice; - - /*Compute change in sensible heat due to changes in pressure melting point*/ - dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); - dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); - PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]); - - meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); - _assert_(meltrate>0.); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight* - ( - meltrate*(1/rho_water-1/rho_ice) - +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap - -beta*sqrt(vx*vx+vy*vy) - +ieb - +storage*head_old/dt - )*basis[i]; - } - /*Clean up and return*/ - xDelete(xyz_list); - xDelete(basis); - delete friction; - 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*/ - IssmDouble dh[3]; - int* doflist = NULL; - IssmDouble* xyz_list = NULL; - - /*Get gravity from parameters*/ - IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); - - /*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); - - /*Get thickness and base on nodes to apply cap on water head*/ - IssmDouble* eff_pressure = xNew(numnodes); - IssmDouble* thickness = xNew(numnodes); - IssmDouble* bed = xNew(numnodes); - IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - element->GetInputListOnNodes(&thickness[0],ThicknessEnum); - element->GetInputListOnNodes(&bed[0],BaseEnum); - - /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/ - IssmDouble* head_old = xNew(numnodes); - element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum); - IssmDouble relaxation; - element->FindParam(&relaxation,HydrologyRelaxationEnum); - - /*Use the dof list to index into the solution vector: */ - for(int i=0;i hrho_ice*thickness[i]/rho_water+bed[i]){ - values[i] = rho_ice*thickness[i]/rho_water+bed[i]; - } - - /*Make sure that negative pressure is not allowed*/ - // if(values[i](values[i])) _error_("NaN found in solution vector"); - if(xIsInf(values[i])) _error_("Inf found in solution vector"); - } - - /*Add input to the element: */ - element->AddInput(HydrologyHeadEnum,values,element->GetElementType()); - element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum); - - /*Update reynolds number according to new solution*/ - element->GetVerticesCoordinates(&xyz_list); - Input* head_input = element->GetInput(HydrologyHeadEnum);_assert_(head_input); - head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list); - IssmDouble conductivity = GetConductivity(element); - - IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU; - element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum); - - /*Free resources:*/ - xDelete(values); - xDelete(thickness); - xDelete(bed); - xDelete(xyz_list); - xDelete(doflist); - xDelete(eff_pressure); - xDelete(head_old); -}/*}}}*/ -void HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ - /*Default, do nothing*/ - return; -}/*}}}*/ - -/*Additional methods*/ -IssmDouble HydrologySommersAnalysis::GetConductivity(Element* element){/*{{{*/ - - /*Intermediaries */ - IssmDouble gap,reynolds; - - /*Get gravity from parameters*/ - IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); - - /*Get Reynolds and gap average values*/ - Input* reynolds_input = element->GetInput(HydrologyReynoldsEnum); _assert_(reynolds_input); - Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); - reynolds_input->GetInputAverage(&reynolds); - gap_input->GetInputAverage(&gap); - - /*Compute conductivity*/ - IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds)); - _assert_(conductivity>0); - - /*Clean up and return*/ - return conductivity; -}/*}}}*/ -void HydrologySommersAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/ - - - for(int j=0;jelements->Size();j++){ - Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); - UpdateGapHeight(element); - } - -}/*}}}*/ -void HydrologySommersAnalysis::UpdateGapHeight(Element* element){/*{{{*/ - - /*Skip if water or ice shelf element*/ - if(element->IsFloating()) return; - - /*Intermediaries */ - IssmDouble newgap = 0.; - IssmDouble Jdet,meltrate,G,dh[2],B,A,n,dt; - IssmDouble gap,bed,thickness,head,ieb; - IssmDouble lr,br,vx,vy,beta; - IssmDouble alpha2,frictionheat; - IssmDouble* xyz_list = NULL; - IssmDouble dpressure_water[2],dbed[2],PMPheat; - IssmDouble q = 0.; - IssmDouble channelization = 0.; - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&dt,TimesteppingTimeStepEnum); - IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); - IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); - IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); - Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); - Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); - Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); - - /*Get conductivity from inputs*/ - IssmDouble conductivity = GetConductivity(element); - - /*Build friction element, needed later: */ - Friction* friction=new Friction(element,2); - - /*Keep track of weights*/ - IssmDouble totalweights=0.; - - /* 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); - - geothermalflux_input->GetInputValue(&G,gauss); - base_input->GetInputValue(&bed,gauss); - base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss); - thickness_input->GetInputValue(&thickness,gauss); - gap_input->GetInputValue(&gap,gauss); - head_input->GetInputValue(&head,gauss); - head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); - englacial_input->GetInputValue(&ieb,gauss); - lr_input->GetInputValue(&lr,gauss); - br_input->GetInputValue(&br,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - - /*Get ice A parameter*/ - B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); - A=pow(B,-n); - - /*Compute beta term*/ - if(gapGetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - frictionheat=alpha2*(vx*vx+vy*vy); - - /*Get water and ice pressures*/ - IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); - IssmDouble pressure_water = rho_water*g*(head-bed); - if(pressure_water>pressure_ice) pressure_water = pressure_ice; - - /* Compute change in sensible heat due to changes in pressure melting point*/ - dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); - dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); - PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]); - - meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); - _assert_(meltrate>0.); - - newgap += gauss->weight*Jdet*(gap+dt*( - meltrate/rho_ice - -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap - +beta*sqrt(vx*vx+vy*vy) - )); - totalweights +=gauss->weight*Jdet; - - /* Compute basal water flux */ - q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])); - - /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */ - channelization += gauss->weight*Jdet*(meltrate/rho_ice/(meltrate/rho_ice+beta*sqrt(vx*vx+vy*vy))); - } - - /*Divide by connectivity*/ - newgap = newgap/totalweights; - IssmDouble mingap = 1e-6; - if(newgapthickness) - newgap = thickness; - - /*Add new gap as an input*/ - element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum); - - /*Divide by connectivity, add basal flux as an input*/ - q = q/totalweights; - element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum); - - /* Divide by connectivity, add degree of channelization as an input */ - channelization = channelization/totalweights; - element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum); - - /*Clean up and return*/ - xDelete(xyz_list); - delete friction; - delete gauss; -}/*}}}*/ Index: sm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h (revision 23019) +++ (revision ) @@ -1,38 +1,0 @@ -/*! \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); - - /*Intermediaries*/ - IssmDouble GetConductivity(Element* element); - void UpdateGapHeight(FemModel* femmodel); - void UpdateGapHeight(Element* element); -}; -#endif Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 23019) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 23020) @@ -30,5 +30,6 @@ #include "./HydrologyDCInefficientAnalysis.h" #include "./HydrologyShreveAnalysis.h" -#include "./HydrologySommersAnalysis.h" +#include "./HydrologyShaktiAnalysis.h" +#include "./HydrologyPismAnalysis.h" #include "./LevelsetAnalysis.h" #include "./MasstransportAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23020) @@ -771,5 +771,5 @@ if(ishydrology){ analyses_temp[numanalyses++]=HydrologyShreveAnalysisEnum; - analyses_temp[numanalyses++]=HydrologySommersAnalysisEnum; + analyses_temp[numanalyses++]=HydrologyShaktiAnalysisEnum; analyses_temp[numanalyses++]=HydrologyDCInefficientAnalysisEnum; analyses_temp[numanalyses++]=HydrologyDCEfficientAnalysisEnum; Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23020) @@ -229,5 +229,5 @@ break; case 8: - GetAlpha2Sommers(palpha2,gauss); + GetAlpha2Shakti(palpha2,gauss); break; case 9: @@ -382,5 +382,5 @@ *palpha2=alpha2; }/*}}}*/ -void Friction::GetAlpha2Sommers(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ +void Friction::GetAlpha2Shakti(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-base)*/ @@ -684,14 +684,16 @@ /*Compute effective pressure first */ IssmDouble N,delta,W,Wmax,e0,Cc,P0; - //element->parameters->FindParam(&delta,TimeEnum); - //element->parameters->FindParam(&e0,TimeEnum); - //element->parameters->FindParam(&P0,TimeEnum); - //element->GetInputValue(&Cc,gauss,FrictionPressureAdjustedTemperatureEnum); - //element->GetInputValue(&W,gauss,FrictionPressureAdjustedTemperatureEnum); + element->parameters->FindParam(&delta,FrictionDeltaEnum); + element->parameters->FindParam(&e0,FrictionVoidRatioEnum); + //element->parameters->FindParam(&P0,XXEnum); + _error_("don't know what P0 is"); + element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum); + element->GetInputValue(&W,gauss,WatercolumnEnum); + Wmax = 0.; _error_("Don't know what to do with Wmax yet"); N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax)); /*Compute yield stress following a Mohr-Colomb criterion*/ IssmDouble phi; - //element->GetInputValue(&phi,gauss,FrictionPressureAdjustedTemperatureEnum); + element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum); IssmDouble tau_c = N*tan(phi); @@ -702,6 +704,6 @@ /*now compute alpha^2*/ IssmDouble u0,q; - //element->parameters->FindParam(&u0,TimeEnum); - //element->parameters->FindParam(&q,TimeEnum); + element->parameters->FindParam(&u0,FrictionThresholdSpeedEnum); + element->parameters->FindParam(&q,FrictionPseudoplasticityExponentEnum); IssmDouble alpha2 = tau_c/(pow(ub,1.-q)*pow(u0,q)); Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 23020) @@ -37,5 +37,5 @@ void GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Josh(IssmDouble* palpha2,Gauss* gauss); - void GetAlpha2Sommers(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2Shakti(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss); Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23020) @@ -179,6 +179,6 @@ switch(analysis_type){ - case HydrologySommersAnalysisEnum: - pe = this->CreatePVectorHydrologySommers(); + case HydrologyShaktiAnalysisEnum: + pe = this->CreatePVectorHydrologyShakti(); break; case HydrologyDCInefficientAnalysisEnum: @@ -337,5 +337,5 @@ /*}}}*/ -ElementVector* Moulin::CreatePVectorHydrologySommers(void){/*{{{*/ +ElementVector* Moulin::CreatePVectorHydrologyShakti(void){/*{{{*/ /*If this node is not the master node (belongs to another partition of the Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 23020) @@ -80,5 +80,5 @@ /*}}}*/ - ElementVector* CreatePVectorHydrologySommers(void); + ElementVector* CreatePVectorHydrologyShakti(void); ElementVector* CreatePVectorHydrologyDCInefficient(void); ElementVector* CreatePVectorHydrologyDCEfficient(void); Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23020) @@ -190,5 +190,5 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ switch(analysis_type){ - case HydrologySommersAnalysisEnum: + case HydrologyShaktiAnalysisEnum: /*Nothing!*/ break; @@ -213,6 +213,6 @@ switch(analysis_type){ - case HydrologySommersAnalysisEnum: - pe=CreatePVectorHydrologySommers(); + case HydrologyShaktiAnalysisEnum: + pe=CreatePVectorHydrologyShakti(); break; default: @@ -346,5 +346,5 @@ /*Neumannflux management*/ -ElementVector* Neumannflux::CreatePVectorHydrologySommers(void){/*{{{*/ +ElementVector* Neumannflux::CreatePVectorHydrologyShakti(void){/*{{{*/ /* constants*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h (revision 23020) @@ -73,5 +73,5 @@ /*}}}*/ /*Neumannflux management:{{{*/ - ElementVector* CreatePVectorHydrologySommers(void); + ElementVector* CreatePVectorHydrologyShakti(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23020) @@ -149,5 +149,5 @@ /*Nothing to add*/ } - else if(hydrology_model==HydrologysommersEnum){ + else if(hydrology_model==HydrologyshaktiEnum){ /*Nothing to add*/ } @@ -250,5 +250,5 @@ /*Nothing to add*/ } - else if(hydrology_model==HydrologysommersEnum){ + else if(hydrology_model==HydrologyshaktiEnum){ /*Nothing to add*/ } Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23020) @@ -116,11 +116,18 @@ } } - else if (hydrology_model==HydrologysommersEnum){ - femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); + else if (hydrology_model==HydrologyshaktiEnum){ + femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum); InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); solutionsequence_shakti_nonlinear(femmodel); if(VerboseSolution()) _printf0_(" updating gap height\n"); - HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis(); + HydrologyShaktiAnalysis* analysis = new HydrologyShaktiAnalysis(); analysis->UpdateGapHeight(femmodel); + delete analysis; + } + else if (hydrology_model==HydrologypismEnum){ + femmodel->SetCurrentConfiguration(HydrologyPismAnalysisEnum); + if(VerboseSolution()) _printf0_(" updating water column\n"); + HydrologyPismAnalysis* analysis = new HydrologyPismAnalysis(); + analysis->UpdateWaterColumn(femmodel); delete analysis; } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23019) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23020) @@ -127,4 +127,8 @@ FrictionGammaEnum, FrictionLawEnum, + FrictionPseudoplasticityExponentEnum, + FrictionThresholdSpeedEnum, + FrictionDeltaEnum, + FrictionVoidRatioEnum, GiaCrossSectionShapeEnum, GroundinglineMigrationEnum, @@ -459,4 +463,6 @@ FrictionQEnum, FrictionWaterLayerEnum, + FrictionTillFrictionAngleEnum, + FrictionSedimentCompressibilityCoefficientEnum, GeometryHydrostaticRatioEnum, GiaLithosphereThicknessEnum, @@ -483,4 +489,5 @@ HydrologyWaterVxEnum, HydrologyWaterVyEnum, + HydrologyDrainageRateEnum, IceEnum, IceMaskNodeActivationEnum, @@ -806,6 +813,8 @@ HydrologyshreveEnum, HydrologySolutionEnum, - HydrologySommersAnalysisEnum, - HydrologysommersEnum, + HydrologyShaktiAnalysisEnum, + HydrologyPismAnalysisEnum, + HydrologyshaktiEnum, + HydrologypismEnum, IceMassEnum, IceMassScaledEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23020) @@ -135,4 +135,8 @@ case FrictionGammaEnum : return "FrictionGamma"; case FrictionLawEnum : return "FrictionLaw"; + case FrictionPseudoplasticityExponentEnum : return "FrictionPseudoplasticityExponent"; + case FrictionThresholdSpeedEnum : return "FrictionThresholdSpeed"; + case FrictionDeltaEnum : return "FrictionDelta"; + case FrictionVoidRatioEnum : return "FrictionVoidRatio"; case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape"; case GroundinglineMigrationEnum : return "GroundinglineMigration"; @@ -465,4 +469,6 @@ case FrictionQEnum : return "FrictionQ"; case FrictionWaterLayerEnum : return "FrictionWaterLayer"; + case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle"; + case FrictionSedimentCompressibilityCoefficientEnum : return "FrictionSedimentCompressibilityCoefficient"; case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio"; case GiaLithosphereThicknessEnum : return "GiaLithosphereThickness"; @@ -489,4 +495,5 @@ case HydrologyWaterVxEnum : return "HydrologyWaterVx"; case HydrologyWaterVyEnum : return "HydrologyWaterVy"; + case HydrologyDrainageRateEnum : return "HydrologyDrainageRate"; case IceEnum : return "Ice"; case IceMaskNodeActivationEnum : return "IceMaskNodeActivation"; @@ -810,6 +817,8 @@ case HydrologyshreveEnum : return "Hydrologyshreve"; case HydrologySolutionEnum : return "HydrologySolution"; - case HydrologySommersAnalysisEnum : return "HydrologySommersAnalysis"; - case HydrologysommersEnum : return "Hydrologysommers"; + case HydrologyShaktiAnalysisEnum : return "HydrologyShaktiAnalysis"; + case HydrologyPismAnalysisEnum : return "HydrologyPismAnalysis"; + case HydrologyshaktiEnum : return "Hydrologyshakti"; + case HydrologypismEnum : return "Hydrologypism"; case IceMassEnum : return "IceMass"; case IceMassScaledEnum : return "IceMassScaled"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23020) @@ -135,10 +135,14 @@ else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum; else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum; - else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; - else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; + else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum; + else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; + if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum; + else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum; + else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; + else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; + else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; else if (strcmp(name,"GroundinglineMeltInterpolation")==0) return GroundinglineMeltInterpolationEnum; else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; @@ -256,12 +260,12 @@ else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; else if (strcmp(name,"SealevelUGia")==0) return SealevelUGiaEnum; - else if (strcmp(name,"SealevelNGia")==0) return SealevelNGiaEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"SealevelNGia")==0) return SealevelNGiaEnum; else if (strcmp(name,"SealevelRSL")==0) return SealevelRSLEnum; else if (strcmp(name,"SealevelRSLEustatic")==0) return SealevelRSLEustaticEnum; else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; + else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; else if (strcmp(name,"SealevelNGiaRate")==0) return SealevelNGiaRateEnum; else if (strcmp(name,"SealevelUGiaRate")==0) return SealevelUGiaRateEnum; @@ -379,12 +383,12 @@ else if (strcmp(name,"TransientIscoupler")==0) return TransientIscouplerEnum; else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; - else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; + else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; @@ -474,4 +478,6 @@ else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; + else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; + else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum; else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; @@ -498,6 +504,10 @@ else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum; else if (strcmp(name,"Ice")==0) return IceEnum; - else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; @@ -506,8 +516,5 @@ else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; + else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; @@ -622,5 +629,8 @@ else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; - else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; @@ -629,8 +639,5 @@ else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; + else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; @@ -745,5 +752,8 @@ else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; - else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; else if (strcmp(name,"EffectivePressureStacked")==0) return EffectivePressureStackedEnum; else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; @@ -752,8 +762,5 @@ else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"EplHeadStacked")==0) return EplHeadStackedEnum; + else if (strcmp(name,"EplHeadStacked")==0) return EplHeadStackedEnum; else if (strcmp(name,"EplHead")==0) return EplHeadEnum; else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; @@ -828,6 +835,8 @@ else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum; else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; - else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum; - else if (strcmp(name,"Hydrologysommers")==0) return HydrologysommersEnum; + else if (strcmp(name,"HydrologyShaktiAnalysis")==0) return HydrologyShaktiAnalysisEnum; + else if (strcmp(name,"HydrologyPismAnalysis")==0) return HydrologyPismAnalysisEnum; + else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum; + else if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum; else if (strcmp(name,"IceMass")==0) return IceMassEnum; else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; @@ -866,5 +875,8 @@ else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; else if (strcmp(name,"LoveKi")==0) return LoveKiEnum; - else if (strcmp(name,"LoveKr")==0) return LoveKrEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"LoveKr")==0) return LoveKrEnum; else if (strcmp(name,"LoveLi")==0) return LoveLiEnum; else if (strcmp(name,"LoveLr")==0) return LoveLrEnum; @@ -875,8 +887,5 @@ else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"Masscon")==0) return MassconEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; + else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; @@ -989,5 +998,8 @@ else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; - else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; @@ -998,8 +1010,5 @@ else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; + else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; @@ -1112,5 +1121,8 @@ else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; - else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; @@ -1121,8 +1133,5 @@ else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; + else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23019) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23020) @@ -203,5 +203,5 @@ case 1: return HydrologydcEnum; case 2: return HydrologyshreveEnum; - case 3: return HydrologysommersEnum; + case 3: return HydrologyshaktiEnum; default: _error_("Marshalled hydrology code \""< - C++ 467180201848084161120661 + C++ 507213172178899398142503 - MATLAB 16068941184034219269536 + MATLAB 8759048133095761679973 - C/C++ Header 420353436621594023136 + Python 391608369733676249818 - Python 162252327511073516009 + C/C++ Header 432371838511752725096 - m4 81052149987711078 + m4 814871531008611726 - Bourne Shell 4113188524825 + Javascript 10414852228964313356 + + + Bourne Shell 613122910321392 + + + Fortran 90533135810291723 @@ -68,5 +74,5 @@ - SUM: 26743418743935163794241916 + SUM: 23354360449191233458326253 Index: /issm/trunk-jpl/src/m/classes/frictionpism.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionpism.m (revision 23020) +++ /issm/trunk-jpl/src/m/classes/frictionpism.m (revision 23020) @@ -0,0 +1,74 @@ +%FRICTIONPISM class definition +% +% Usage: +% frictionpism=frictionpism(); + +classdef frictionpism + properties (SetAccess=public) + pseudoplasticity_exponent = 0.; + threshold_speed = 0.; + delta = 0.; + void_ratio = 0.; + till_friction_angle = NaN; + sediment_compressibility_coefficient = NaN; + end + methods + function self = extrude(self,md) % {{{ + self.till_friction_angle=project3d(md,'vector',self.till_friction_angle,'type','node','layer',1); + self.sediment_compressibility_coefficient=project3d(md,'vector',self.sediment_compressibility_coefficient,'type','node','layer',1); + end % }}} + function self = frictionpism(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + self=structtoobj(frictionpism(),varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + + self.pseudoplasticity_exponent = 0.6; + self.threshold_speed = 100.; + self.delta = 0.02; + self.void_ratio = 0.; %unknown!! + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + %Early return + if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end + if (strcmp(solution,'TransientSolution') & md.transient.isstressbalance ==0 & md.transient.isthermal == 0), return; end + + md = checkfield(md,'fieldname','friction.pseudoplasticity_exponent','numel',[1],'>',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','friction.threshold_speed','numel',[1],'>',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','friction.delta','numel',[1],'>',0,'<',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','friction.void_ratio','numel',[1],'>',0,'<',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','friction.till_friction_angle','NaN',1,'Inf',1,'<',360.,'>',0.,'size',[md.mesh.vertices 1]); %TODO: check rad or deg + end % }}} + function disp(self) % {{{ + disp(sprintf('Basal shear stress parameters for the PISM friction law (See Aschwwanden et al. 2016 for more details)')); + fielddisplay(self,'pseudoplasticity_exponent','pseudoplasticity exponent [dimensionless]'); + fielddisplay(self,'threshold_speed','threshold speed [m/yr]'); + fielddisplay(self,'delta','lower limit of the effective pressure, expressed as a fraction of overburden pressure [dimensionless]'); + fielddisplay(self,'void_ratio','void ratio at a reference effective pressure [dimensionless]'); + fielddisplay(self,'till_friction_angle','till friction angle [deg or rad??]'); + fielddisplay(self,'sediment_compressibility_coefficient','coefficient of compressibility of the sediment [dimensionless?]'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.friction.law','data',10,'format','Integer'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','pseudoplasticity_exponent','format','Double'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','threshold_speed','format','Double','scale',1./yts); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','delta','format','Double'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','void_ratio','format','Double'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','till_friction_angle','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','sediment_compressibility_coefficient','format','DoubleMat','mattype',1); + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + error('not implemented yet!'); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/frictionshakti.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionshakti.m (revision 23020) +++ /issm/trunk-jpl/src/m/classes/frictionshakti.m (revision 23020) @@ -0,0 +1,45 @@ +%FRICTIONSHAKTI class definition +% +% Usage: +% friction=frictionshakti(); + +classdef frictionshakti + properties (SetAccess=public) + coefficient = NaN; + end + methods + function self = extrude(self,md) % {{{ + self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1); + end % }}} + function self = frictionshakti(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + self=structtoobj(frictionshakti(),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('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end + md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1); + end % }}} + function disp(self) % {{{ + disp(sprintf('Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))')); + fielddisplay(self,'coefficient','friction coefficient [SI]'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/frictionshakti.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionshakti.py (revision 23020) +++ /issm/trunk-jpl/src/m/classes/frictionshakti.py (revision 23020) @@ -0,0 +1,47 @@ +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +class frictionshakti(object): + """ + FRICTIONSHAKTI class definition + + Usage: + friction=frictionshakti() + """ + + def __init__(self,md): # {{{ + self.coefficient = md.friction.coefficient + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))" + + string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]")) + return string + #}}} + def extrude(self,md): # {{{ + self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) + return self + #}}} + def setdefaultparameters(self): # {{{ + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + #Early return + if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: + return md + + md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) + return md + + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + # }}} Index: sm/trunk-jpl/src/m/classes/frictionsommers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionsommers.m (revision 23019) +++ (revision ) @@ -1,45 +1,0 @@ -%FRICTIONSOMMERS class definition -% -% Usage: -% friction=frictionsommers(); - -classdef frictionsommers - properties (SetAccess=public) - coefficient = NaN; - end - methods - function self = extrude(self,md) % {{{ - self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1); - end % }}} - function self = frictionsommers(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - case 1 - self=structtoobj(frictionsommers(),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('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end - md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1); - end % }}} - function disp(self) % {{{ - disp(sprintf('Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))')); - fielddisplay(self,'coefficient','friction coefficient [SI]'); - end % }}} - function marshall(self,prefix,md,fid) % {{{ - yts=md.constants.yts; - - WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer'); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - - end % }}} - end -end Index: sm/trunk-jpl/src/m/classes/frictionsommers.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionsommers.py (revision 23019) +++ (revision ) @@ -1,47 +1,0 @@ -from fielddisplay import fielddisplay -from project3d import project3d -from checkfield import checkfield -from WriteData import WriteData - -class frictionsommers(object): - """ - FRICTIONSOMMERS class definition - - Usage: - friction=frictionsommers() - """ - - def __init__(self,md): # {{{ - self.coefficient = md.friction.coefficient - #set defaults - self.setdefaultparameters() - - #}}} - def __repr__(self): # {{{ - string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))" - - string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]")) - return string - #}}} - def extrude(self,md): # {{{ - self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) - return self - #}}} - def setdefaultparameters(self): # {{{ - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ - - #Early return - if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: - return md - - md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) - return md - - # }}} - def marshall(self,prefix,md,fid): # {{{ - yts=md.constants.yts - WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - # }}} Index: /issm/trunk-jpl/src/m/classes/hydrologypism.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologypism.m (revision 23019) +++ /issm/trunk-jpl/src/m/classes/hydrologypism.m (revision 23020) @@ -46,5 +46,5 @@ WriteData(fid,prefix,'name','md.hydrology.model','data',4,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','drainage_rate','format','DoubleMat','mattype',1,'scale',1./(1000.*yts)); %from mm/yr to m/s + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','drainage_rate','format','DoubleMat','mattype',1,'scale',1./(1000.*yts)); %from mm/yr to m/s end % }}} end Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.js =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.js (revision 23019) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.js (revision 23020) @@ -79,5 +79,5 @@ }else if(solutiontype ==='TransientSolution'){ - analyses=['StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis','ThermalAnalysis','MeltingAnalysis','EnthalpyAnalysis','MasstransportAnalysis','HydrologySommersAnalysis']; + analyses=['StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis','ThermalAnalysis','MeltingAnalysis','EnthalpyAnalysis','MasstransportAnalysis','HydrologyShaktiAnalysis']; }else if(solutiontype ==='SealevelriseSolution'){ Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 23019) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 23020) @@ -71,5 +71,5 @@ analyses={'EsaAnalysis'}; elseif strcmp(solutiontype,'TransientSolution') - analyses={'StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis','ThermalAnalysis','MeltingAnalysis','EnthalpyAnalysis','MasstransportAnalysis','HydrologySommersAnalysis'}; + analyses={'StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis','ThermalAnalysis','MeltingAnalysis','EnthalpyAnalysis','MasstransportAnalysis','HydrologyShaktiAnalysis'}; elseif strcmp(solutiontype,'SealevelriseSolution') analyses={'SealevelriseAnalysis'};