Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16637) +++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16638) @@ -56,6 +56,8 @@ iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); iomodel->FetchDataToInput(elements,EplHeadEnum); + iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum); + }/*}}}*/ void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16637) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16638) @@ -113,7 +113,8 @@ HydrologydcEplCompressibilityEnum, HydrologydcEplPorosityEnum, HydrologydcEplThicknessEnum, - HydrologydcEplTransmitivityEnum, + HydrologydcEplThicknessOldEnum, + HydrologydcEplConductivityEnum, HydrologydcIsefficientlayerEnum, HydrologydcSedimentlimitFlagEnum, HydrologydcSedimentlimitEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16637) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16638) @@ -121,7 +121,8 @@ case HydrologydcEplCompressibilityEnum : return "HydrologydcEplCompressibility"; case HydrologydcEplPorosityEnum : return "HydrologydcEplPorosity"; case HydrologydcEplThicknessEnum : return "HydrologydcEplThickness"; - case HydrologydcEplTransmitivityEnum : return "HydrologydcEplTransmitivity"; + case HydrologydcEplThicknessOldEnum : return "HydrologydcEplThicknessOld"; + case HydrologydcEplConductivityEnum : return "HydrologydcEplConductivity"; case HydrologydcIsefficientlayerEnum : return "HydrologydcIsefficientlayer"; case HydrologydcSedimentlimitFlagEnum : return "HydrologydcSedimentlimitFlag"; case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16637) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16638) @@ -121,7 +121,8 @@ else if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum; else if (strcmp(name,"HydrologydcEplPorosity")==0) return HydrologydcEplPorosityEnum; else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum; - else if (strcmp(name,"HydrologydcEplTransmitivity")==0) return HydrologydcEplTransmitivityEnum; + else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum; + else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum; else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum; else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum; else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum; @@ -135,11 +136,11 @@ else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum; else if (strcmp(name,"WaterTransfer")==0) return WaterTransferEnum; else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum; - else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; + if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; + else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum; else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum; @@ -258,11 +259,11 @@ else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; - else if (strcmp(name,"Surface")==0) return SurfaceEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; + if (strcmp(name,"Surface")==0) return SurfaceEnum; + else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; else if (strcmp(name,"Element")==0) return ElementEnum; - else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"FileParam")==0) return FileParamEnum; + if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; + else if (strcmp(name,"FileParam")==0) return FileParamEnum; else if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum; else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; - else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; + else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; Index: ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp (revision 16637) +++ ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp (revision 16638) @@ -8,7 +8,6 @@ #include "../../toolkits/toolkits.h" void InputDuplicatex(FemModel* femmodel,int original_enum, int new_enum){ - /*Go through elemnets, and ask to reinitialie the input: */ for(int i=0;ielements->Size();i++){ Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16637) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16638) @@ -128,9 +128,12 @@ femmodel->HydrologyEPLupdateDomainx(); femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); + /*Reset constraint on the ZigZag Lock*/ + ResetConstraintsx(femmodel); /*Iteration on the EPL layer*/ eplconverged = false; for(;;){ + femmodel->HydrologyEPLThicknessx(); femmodel->HydrologyTransferx(); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16637) +++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16638) @@ -83,6 +83,7 @@ femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); if (isefficientlayer){ InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum); + InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum); } /*Proceed now to heads computations*/ @@ -91,8 +92,8 @@ if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ if(VerboseSolution()) _printf0_(" saving results \n"); if(isefficientlayer){ - int outputs[6] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6); + int outputs[7] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],7); } else{ int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum}; Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16637) +++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16638) @@ -144,6 +144,20 @@ return pow(B,-n); } /*}}}*/ +/*FUNCTION Matice::GetAbar {{{*/ +IssmDouble Matice::GetAbar(){ + /* + * A = 1/B^n + */ + + IssmDouble B,n; + + inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum); + n=this->GetN(); + + return pow(B,-n); +} +/*}}}*/ /*FUNCTION Matice::GetB {{{*/ IssmDouble Matice::GetB(){ Index: ../trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16637) +++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16638) @@ -35,6 +35,7 @@ virtual void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; virtual IssmDouble GetA()=0; + virtual IssmDouble GetAbar()=0; virtual IssmDouble GetB()=0; virtual IssmDouble GetBbar()=0; virtual IssmDouble GetN()=0; Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16637) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16638) @@ -65,8 +65,7 @@ if(isefficientlayer){ iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum); iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum); - iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum); - iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum); + iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum); } } else{ @@ -353,9 +352,9 @@ ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity)); } /*}}}*/ -/*FUNCTION Matpar::GetEplStoring {{{*/ -IssmDouble Matpar::GetEplStoring(){ - return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness* +/*FUNCTION Matpar::GetEplSpecificStoring {{{*/ +IssmDouble Matpar::GetEplSpecificStoring(){ + return this->rho_freshwater* this->g* this->epl_porosity* ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity)); } /*}}}*/ @@ -368,10 +367,10 @@ IssmDouble Matpar::GetSedimentThickness(){ return sediment_thickness; } -/*}}}*/ -/*FUNCTION Matpar::GetEplTransitivity {{{*/ -IssmDouble Matpar::GetEplTransmitivity(){ - return epl_transmitivity; +/*}}}*/ +/*FUNCTION Matpar::GetEplConductivity {{{*/ +IssmDouble Matpar::GetEplConductivity(){ + return epl_conductivity; } /*}}}*/ /*FUNCTION Matpar::TMeltingPoint {{{*/ Index: ../trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16637) +++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16638) @@ -62,6 +62,7 @@ void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon); void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); IssmDouble GetA(); + IssmDouble GetAbar(); IssmDouble GetB(); IssmDouble GetBbar(); IssmDouble GetD(); Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16637) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16638) @@ -48,8 +48,7 @@ IssmDouble epl_compressibility; IssmDouble epl_porosity; - IssmDouble epl_thickness; - IssmDouble epl_transmitivity; + IssmDouble epl_conductivity; /*gia: */ IssmDouble lithosphere_shear_modulus; @@ -93,6 +92,7 @@ void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; IssmDouble GetA(){_error_("not supported");}; + IssmDouble GetAbar(){_error_("not supported");}; IssmDouble GetB(){_error_("not supported");}; IssmDouble GetBbar(){_error_("not supported");}; IssmDouble GetN(){_error_("not supported");}; @@ -120,10 +120,10 @@ IssmDouble GetHydrologyCR(); IssmDouble GetHydrologyN(); IssmDouble GetSedimentStoring(); - IssmDouble GetEplStoring(); + IssmDouble GetEplSpecificStoring(); IssmDouble GetSedimentTransmitivity(); IssmDouble GetSedimentThickness(); - IssmDouble GetEplTransmitivity(); + IssmDouble GetEplConductivity(); IssmDouble TMeltingPoint(IssmDouble pressure); IssmDouble PureIceEnthalpy(IssmDouble pressure); IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16637) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16638) @@ -1577,8 +1577,9 @@ void Tria::InputDuplicate(int original_enum,int new_enum){ /*Call inputs method*/ - if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); - + if (IsInput(original_enum)) { + inputs->DuplicateInput(original_enum,new_enum); + } } /*}}}*/ /*FUNCTION Tria::InputScale{{{*/ @@ -2096,6 +2097,8 @@ name==SedimentHeadEnum || name==EplHeadOldEnum || name==EplHeadEnum || + name==HydrologydcEplThicknessOldEnum || + name==HydrologydcEplThicknessEnum || name==HydrologydcMaskEplactiveEnum || name==MeshVertexonbedEnum || name==WaterTransferEnum || @@ -2454,10 +2457,6 @@ this->ComputeStressTensor(); input=this->inputs->GetInput(output_enum); break; - case SurfaceNormalVelocityEnum: - this->ComputeSurfaceNormalVelocity(); - input=this->inputs->GetInput(output_enum); - break; default: _error_("input "<(connectivity[i]); vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); - - /*FIXME: Inputs needs to be updated first in Tria::ResultInterpolation, this is a hack*/ - this->inputs->DeleteInput(SurfaceNormalVelocityEnum); break; } default: @@ -2650,53 +2646,6 @@ *(surface_normal+2) = normal[2]/normal_norm; } /*}}}*/ -/*FUNCTION Tria::ComputeSurfaceNormalVelocity{{{*/ -void Tria::ComputeSurfaceNormalVelocity(){ - - IssmDouble sum,tangential_vector[2],normal_vector[2],time,ux,uy; - IssmDouble normal_velocity[NUMVERTICES],xyz_list[NUMVERTICES][3]; - IssmDouble value[NUMVERTICES],verticesonsurface[NUMVERTICES]; - - for(int iv=0;ivGetInput(SurfaceSlopeXEnum); _assert_(slope_input); - // Input* slope_input= inputs->GetInput(SurfaceEnum); _assert_(slope_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - - /*Get list of nodes on surface*/ - GetInputListOnVertices(&verticesonsurface[0],MeshVertexonsurfaceEnum); - sum = verticesonsurface[0]+verticesonsurface[1]+verticesonsurface[2]; - _assert_(sum==0. || sum==1. || sum==2.); - - /*Compute normal velocity for surface nodes from L2 projected slope*/ - for(int iv=0;ivGaussNode(P1Enum,iv); - slope_input->GetInputValue(&value[iv],gauss); - // slope_input->GetInputDerivativeValue(&value[iv],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&ux,gauss); - vy_input->GetInputValue(&uy,gauss); - tangential_vector[0]=sqrt(1./(pow(value[iv],2.)+1.)); - tangential_vector[1]=value[iv]*tangential_vector[0]; - normal_vector[0]=-1.*tangential_vector[1]; - normal_vector[1]=tangential_vector[0]; - normal_velocity[iv]=ux*normal_vector[0]+uy*normal_vector[1]; - } - } - - delete gauss; - this->inputs->AddInput(new TriaInput(SurfaceNormalVelocityEnum,&normal_velocity[0],P1Enum)); - -} -/*}}}*/ /*FUNCTION Tria::TimeAdapt{{{*/ IssmDouble Tria::TimeAdapt(void){ @@ -3277,7 +3226,7 @@ /*Now get the average SMB over the element*/ Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input); smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1 - Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 + Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 /*Return: */ return Total_Smb; @@ -6733,9 +6682,10 @@ ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){ /* Intermediaries */ - IssmDouble D_scalar,Jdet; - IssmDouble epl_transmitivity,dt; - IssmDouble epl_storing; + IssmDouble D_scalar,Jdet,dt; + IssmDouble epl_thickness; + IssmDouble epl_conductivity; + IssmDouble epl_specificstoring; IssmDouble xyz_list[NUMVERTICES][3]; /*Check that all nodes are active, else return empty matrix*/ @@ -6754,19 +6704,23 @@ /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - epl_storing = matpar->GetEplStoring(); - epl_transmitivity = matpar->GetEplTransmitivity(); + epl_specificstoring = matpar->GetEplSpecificStoring(); + epl_conductivity = matpar->GetEplConductivity(); + /* Start looping on the number of gaussian points: */ GaussTria* gauss=new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ - + + gauss->GaussPoint(ig); GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + thickness_input->GetInputValue(&epl_thickness,gauss); /*Diffusivity*/ - D_scalar=epl_transmitivity*gauss->weight*Jdet; + D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet; if(reCast(dt)) D_scalar=D_scalar*dt; D[0][0]=D_scalar; D[0][1]=0.; D[1][0]=0.; D[1][1]=D_scalar; @@ -6779,7 +6733,7 @@ /*Transient*/ if(reCast(dt)){ GetNodalFunctions(basis,gauss); - D_scalar=epl_storing*gauss->weight*Jdet; + D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet; TripleMultiply(basis,numnodes,1,0, &D_scalar,1,1,0, @@ -6915,7 +6869,8 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble dt,scalar,water_head; IssmDouble transfer,residual; - IssmDouble epl_storing; + IssmDouble epl_thickness; + IssmDouble epl_specificstoring; GaussTria* gauss = NULL; /*Check that all nodes are active, else return empty matrix*/ @@ -6931,11 +6886,12 @@ IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ - epl_storing = matpar->GetEplStoring(); + epl_specificstoring = matpar->GetEplSpecificStoring(); GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); - Input* transfer_input=inputs->GetInput(WaterTransferEnum); _assert_(transfer_input); + Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); + Input* transfer_input=inputs->GetInput(WaterTransferEnum); _assert_(transfer_input); + Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); Input* old_wh_input=NULL; if(reCast(dt)){ @@ -6958,8 +6914,9 @@ /*Transient term*/ if(reCast(dt)){ + thickness_input->GetInputValue(&epl_thickness,gauss); old_wh_input->GetInputValue(&water_head,gauss); - scalar = Jdet*gauss->weight*water_head*epl_storing; + scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness; for(int i=0;ivalues[i]+=scalar*basis[i]; } } @@ -7141,8 +7098,8 @@ IssmDouble sed_trans,sed_thick; IssmDouble leakage,h_max; IssmDouble wh_trans; - IssmDouble activeEpl[numdof]; - IssmDouble eplstoring[numdof],sedstoring[numdof]; + IssmDouble activeEpl[numdof],epl_thickness[numdof]; + IssmDouble epl_specificstoring[numdof],sedstoring[numdof]; IssmDouble epl_head[numdof],sed_head[numdof]; GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -7164,6 +7121,7 @@ GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); GetInputListOnVertices(&epl_head[0],EplHeadEnum); + GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum); @@ -7176,12 +7134,12 @@ wh_trans=0.0; } else{ - eplstoring[i]=matpar->GetEplStoring(); + epl_specificstoring[i]=matpar->GetEplSpecificStoring(); sedstoring[i]=matpar->GetSedimentStoring(); /*EPL head higher than sediment head, transfer from the epl to the sediment*/ if(epl_head[i]>sed_head[i]){ - wh_trans=eplstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick); /*No transfer if the sediment head is allready at the maximum*/ this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); @@ -7245,7 +7203,7 @@ GetInputListOnVertices(&eplhead[0],EplHeadEnum); GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum); - /*Get minimum sediment head*/ + /*Get minimum sediment head of the element*/ sedheadmin=sedhead[0]; for(i=1;iparameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + if(isefficientlayer){ + /*For now, assuming just one way to compute EPL thickness*/ + rho_water = matpar->GetRhoWater(); + rho_ice = matpar->GetRhoIce(); + gravity = matpar->GetG(); + latentheat = matpar->GetLatentHeat(); + n = material->GetN(); + A = material->GetAbar(); + + GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + GetInputListOnVertices(&eplhead[0],EplHeadEnum); + GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); + GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); + GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); + GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); + GetInputListOnVertices(&bed[0],BedEnum); + + epl_conductivity = matpar->GetEplConductivity(); + + for(int i=0;iinputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); + } } /*}}}*/ Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16637) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16638) @@ -2644,6 +2644,8 @@ name==EplHeadEnum || name==SedimentHeadOldEnum || name==EplHeadOldEnum || + name==HydrologydcEplThicknessOldEnum || + name==HydrologydcEplThicknessEnum || name==HydrologydcMaskEplactiveEnum || name==WaterTransferEnum @@ -2870,7 +2872,7 @@ switch(input->GetResultInterpolation()){ case P0Enum: - _error_("not implemented..."); + _error_("P0 not implemented yet for input "<ComputeEPLThickness(); - delete tria->material; delete tria; + /*Get the flag to know if the efficient layer is present*/ + this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + if(isefficientlayer){ + /*For now, assuming just one way to compute EPL thickness*/ + rho_water = matpar->GetRhoWater(); + rho_ice = matpar->GetRhoIce(); + gravity = matpar->GetG(); + latentheat = matpar->GetLatentHeat(); + epl_conductivity = matpar->GetEplConductivity(); + n = material->GetN(); + A = material->GetA(); + + + GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + GetInputListOnVertices(&eplhead[0],EplHeadEnum); + GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); + GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); + GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); + GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); + GetInputListOnVertices(&bed[0],BedEnum); + + for(int i=0;iinputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); + + /*Stop if we have reached the surface*/ + if (penta->IsOnSurface()) break; + + /* get upper Penta*/ + penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); + } + } } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/ Index: ../trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16637) +++ ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16638) @@ -26,7 +26,7 @@ epl_compressibility = 0; epl_porosity = 0; epl_thickness = 0; - epl_transmitivity = 0; + epl_conductivity = 0; end methods % {{{ function obj = hydrologydc(varargin) @@ -60,7 +60,7 @@ obj.epl_compressibility = 1.0e-08; obj.epl_porosity = 0.4; obj.epl_thickness = 1.0; - obj.epl_transmitivity = 8.0e-02; + obj.epl_conductivity = 8.0e-02; end % }}} @@ -96,7 +96,7 @@ md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1); md = checkfield(md,'hydrology.epl_porosity','>',0,'numel',1); md = checkfield(md,'hydrology.epl_thickness','>',0,'numel',1); - md = checkfield(md,'hydrology.epl_transmitivity','>',0,'numel',1); + md = checkfield(md,'hydrology.epl_conductivity','>',0,'numel',1); end end % }}} @@ -136,8 +136,8 @@ fielddisplay(obj,'mask_eplactive','active (1) or not (0) EPL'); fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]'); fielddisplay(obj,'epl_porosity','epl [dimensionless]'); - fielddisplay(obj,'epl_thickness','epl thickness [m]'); - fielddisplay(obj,'epl_transmitivity','epl transmitivity [m^2/s]'); + fielddisplay(obj,'epl_thickness','epl initial thickness [m]'); + fielddisplay(obj,'epl_conductivity','epl conductivity [m^2/s]'); end end @@ -171,7 +171,8 @@ WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double'); WriteData(fid,'object',obj,'fieldname','epl_porosity','format','Double'); WriteData(fid,'object',obj,'fieldname','epl_thickness','format','Double'); - WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double'); + %WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double'); + WriteData(fid,'object',obj,'fieldname','epl_conductivity','format','Double'); end end % }}} Index: ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m =================================================================== --- ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m (revision 0) +++ ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m (revision 16638) @@ -0,0 +1,11 @@ +function macro=HydrologydcEplConductivityEnum() +%HYDROLOGYDCEPLCONDUCTIVITYENUM - Enum of HydrologydcEplConductivity +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologydcEplConductivityEnum() + +macro=StringToEnum('HydrologydcEplConductivity'); Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16637) +++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16638) @@ -113,7 +113,8 @@ def HydrologydcEplCompressibilityEnum(): return StringToEnum("HydrologydcEplCompressibility")[0] def HydrologydcEplPorosityEnum(): return StringToEnum("HydrologydcEplPorosity")[0] def HydrologydcEplThicknessEnum(): return StringToEnum("HydrologydcEplThickness")[0] -def HydrologydcEplTransmitivityEnum(): return StringToEnum("HydrologydcEplTransmitivity")[0] +def HydrologydcEplThicknessOldEnum(): return StringToEnum("HydrologydcEplThicknessOld")[0] +def HydrologydcEplConductivityEnum(): return StringToEnum("HydrologydcEplConductivity")[0] def HydrologydcIsefficientlayerEnum(): return StringToEnum("HydrologydcIsefficientlayer")[0] def HydrologydcSedimentlimitFlagEnum(): return StringToEnum("HydrologydcSedimentlimitFlag")[0] def HydrologydcSedimentlimitEnum(): return StringToEnum("HydrologydcSedimentlimit")[0] Index: ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m =================================================================== --- ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m (revision 0) +++ ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m (revision 16638) @@ -0,0 +1,11 @@ +function macro=HydrologydcEplThicknessOldEnum() +%HYDROLOGYDCEPLTHICKNESSOLDENUM - Enum of HydrologydcEplThicknessOld +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologydcEplThicknessOldEnum() + +macro=StringToEnum('HydrologydcEplThicknessOld');