Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15061) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15062) @@ -1485,4 +1485,10 @@ case HydrologyShreveAnalysisEnum: GetSolutionFromInputsHydrologyShreve(solution); + break; + case HydrologyDCInefficientAnalysisEnum: + GetSolutionFromInputsHydrologyDCInefficient(solution); + break; + case HydrologyDCEfficientAnalysisEnum: + GetSolutionFromInputsHydrologyDCEfficient(solution); break; #endif @@ -6381,4 +6387,69 @@ } /*}}}*/ +/*FUNCTION Tria::GetSolutionFromInputsHydrologyDCInefficient{{{*/ +void Tria::GetSolutionFromInputsHydrologyDCInefficient(Vector* solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int *doflist = NULL; + IssmDouble sedimenthead; + IssmDouble values[numdof]; + GaussTria *gauss = NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Get inputs*/ + Input* sedimenthead_input=inputs->GetInput(SedimentHeadEnum); _assert_(sedimenthead_input); + + /*P1 element only for now*/ + gauss=new GaussTria(); + for(i=0;iGaussVertex(i); + sedimenthead_input->GetInputValue(&sedimenthead,gauss); + values[i]=sedimenthead; + } + + solution->SetValues(numdof,doflist,values,INS_VAL); + + /*Free ressources:*/ + delete gauss; + xDelete(doflist); +} +/*}}}*/ +/*FUNCTION Tria::GetSolutionFromInputsHydrologyDCEfficient{{{*/ +void Tria::GetSolutionFromInputsHydrologyDCEfficient(Vector* solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int *doflist = NULL; + IssmDouble eplhead; + IssmDouble values[numdof]; + GaussTria *gauss = NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Get inputs*/ + Input* eplhead_input=inputs->GetInput(EplHeadEnum); _assert_(eplhead_input); + + /*P1 element only for now*/ + gauss=new GaussTria(); + for(i=0;iGaussVertex(i); + eplhead_input->GetInputValue(&eplhead,gauss); + values[i]=eplhead; + } + + solution->SetValues(numdof,doflist,values,INS_VAL); + + /*Free ressources:*/ + delete gauss; + xDelete(doflist); +} +/*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionHydrologyShreve{{{*/ void Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15061) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15062) @@ -249,4 +249,6 @@ ElementVector* CreatePVectorHydrologyDCEfficient(void); void GetSolutionFromInputsHydrologyShreve(Vector* solution); + void GetSolutionFromInputsHydrologyDCInefficient(Vector* solution); + void GetSolutionFromInputsHydrologyDCEfficient(Vector* solution); void CreateHydrologyWaterVelocityInput(void); void InputUpdateFromSolutionHydrology(IssmDouble* solution); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp (revision 15061) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp (revision 15062) @@ -18,4 +18,5 @@ IssmDouble penalty_factor; IssmDouble leakagefactor; + IssmDouble rel_tol; /*Get parameters: */ @@ -35,4 +36,5 @@ iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum); iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum); + iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum); if(sedimentlimit_flag==1){ @@ -51,4 +53,5 @@ parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag)); parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag)); + parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); /*Assign output pointer: */ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15061) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15062) @@ -90,4 +90,5 @@ SedimentHeadResidualEnum, EplHeadEnum, + HydrologydcRelTolEnum, HydrologydcSpcsedimentHeadEnum, HydrologydcSedimentCompressibilityEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15061) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15062) @@ -98,4 +98,5 @@ case SedimentHeadResidualEnum : return "SedimentHeadResidual"; case EplHeadEnum : return "EplHead"; + case HydrologydcRelTolEnum : return "HydrologydcRelTol"; case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead"; case HydrologydcSedimentCompressibilityEnum : return "HydrologydcSedimentCompressibility"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15061) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15062) @@ -98,4 +98,5 @@ else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum; else if (strcmp(name,"EplHead")==0) return EplHeadEnum; + else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum; else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum; else if (strcmp(name,"HydrologydcSedimentCompressibility")==0) return HydrologydcSedimentCompressibilityEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum; else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum; - else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum; + if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum; + else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum; else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum; else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum; - else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum; + if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum; + else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum; else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum; else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; - else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; + if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum; + else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; else if (strcmp(name,"StringParam")==0) return StringParamEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum; else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum; - else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; + if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum; + else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; else if (strcmp(name,"Step")==0) return StepEnum; else if (strcmp(name,"Time")==0) return TimeEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15061) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15062) @@ -9,35 +9,62 @@ void solutionsequence_hydro_nonlinear(FemModel* femmodel){ - + /*solution : */ Vector* ug=NULL; - Vector* uf=NULL; - Vector* uf_old=NULL; + Vector* uf_sed=NULL; + Vector* uf_epl=NULL; + Vector* old_uf=NULL; + Vector* aged_uf_sed=NULL; + Vector* aged_uf_epl=NULL; Vector* ys=NULL; - IssmDouble sediment_kmax,time; - - /*intermediary: */ + Vector* duf=NULL; + Matrix* Kff=NULL; Matrix* Kfs=NULL; Vector* pf=NULL; Vector* df=NULL; - - bool converged,isefficientlayer; - int constraints_converged; - int num_unstable_constraints; - int count; - int hydro_maxiter; + + bool sedconverged,eplconverged,hydroconverged; + bool isefficientlayer; + int constraints_converged; + int num_unstable_constraints; + int count; + int hydro_maxiter; + IssmDouble sediment_kmax,time; + IssmDouble eps_hyd; + IssmDouble ndu_sed,nu_sed; + IssmDouble ndu_epl,nu_epl; /*Recover parameters: */ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME + femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); + femmodel->parameters->FindParam(&time,TimeEnum); femmodel->BasisIntegralsx(); hydro_maxiter=100; count=1; - converged=false; + sedconverged=false; + eplconverged=false; + hydroconverged=false; - femmodel->parameters->FindParam(&time,TimeEnum); + /*Iteration on the two layers + transfer*/ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf_sed, ug, femmodel->nodes,femmodel->parameters); _assert_(uf_sed); + if(isefficientlayer) { + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf_epl, ug, femmodel->nodes,femmodel->parameters); + } + + for(;;){ - /*First layer*/ + //save pointer to old velocity + delete aged_uf_sed;aged_uf_sed=uf_sed; + if(isefficientlayer){ + delete aged_uf_epl; + aged_uf_epl=uf_epl; + } + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum); @@ -46,20 +73,21 @@ femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); femmodel->HydrologyTransferx(); + + /*Iteration on the sediment layer*/ for(;;){ - femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); - Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf; - Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters); - delete uf_old; uf_old=uf->Duplicate(); + Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf_sed; + Solverx(&uf_sed, Kff, pf,old_uf, df, femmodel->parameters); + delete old_uf; old_uf=uf_sed->Duplicate(); delete Kff; delete pf; delete ug; delete df; - Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; + + Mergesolutionfromftogx(&ug,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); - ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - if (!converged){ + if (!sedconverged){ if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); - if(num_unstable_constraints==0) converged = true; + if(num_unstable_constraints==0) sedconverged = true; if (count>=hydro_maxiter){ _error_(" maximum number of iterations (" << hydro_maxiter << ") exceeded"); @@ -68,7 +96,7 @@ count++; - if(converged){ + if(sedconverged){ femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); - InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sedconverged,ConvergedEnum); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum); @@ -84,20 +112,21 @@ femmodel->UpdateConstraintsx(); femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); - converged = false; + eplconverged = false; + /*Iteration on the EPL layer*/ for(;;){ femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); - Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf; - Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters); - delete uf_old; uf_old=uf->Duplicate(); + Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf_epl; + Solverx(&uf_epl, Kff, pf,old_uf, df, femmodel->parameters); + delete old_uf; old_uf=uf_epl->Duplicate(); delete Kff;delete pf; delete ug; delete df; - Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; + Mergesolutionfromftogx(&ug,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - if (!converged){ + if (!eplconverged){ if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); - if(num_unstable_constraints==0) converged = true; + if(num_unstable_constraints==0) eplconverged = true; if (count>=hydro_maxiter){ _error_(" maximum number of iterations (" << hydro_maxiter << ") exceeded"); @@ -106,6 +135,6 @@ count++; - if(converged){ - InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + if(eplconverged){ + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum); InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); @@ -115,12 +144,53 @@ /*System convergence check*/ - break; + if(!hydroconverged){ + //compute norm(du)/norm(u) + _assert_(aged_uf_sed); _assert_(uf_sed); + duf=uf_sed->Duplicate(); _assert_(duf); + aged_uf_sed->Copy(duf); duf->AYPX(uf_sed,-1.0); + ndu_sed=duf->Norm(NORM_TWO); nu_sed=aged_uf_sed->Norm(NORM_TWO); + if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("convergence criterion is NaN!"); + + if(isefficientlayer){ + duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0); + ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO); + if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("convergence criterion is NaN!"); + } + //print + if (!xIsNan(eps_hyd)){ + if((ndu_sed/nu_sed) " << eps_hyd*100 << " %"); + hydroconverged=false; + } + if(isefficientlayer){ + if((ndu_epl/nu_epl) " << eps_hyd*100 << " %"); + hydroconverged=false; + } + } + } + else _pprintLine_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu_sed/nu_sed*100 << " %"); + } + + if(hydroconverged)break; } InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - + /*Free ressources: */ delete ug; - delete uf; - delete uf_old; + delete uf_sed; + delete uf_epl; + delete old_uf; + delete aged_uf_sed; + delete aged_uf_epl; + delete duf; + } Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 15061) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 15062) @@ -6,4 +6,13 @@ classdef hydrologydc properties (SetAccess=public) + water_compressibility = 0; + isefficientlayer = 0; + penalty_factor = 0; + rel_tol = 0; + sedimentlimit_flag = 0; + sedimentlimit = 0; + transfer_flag = 0; + leakage_factor = 0; + spcsediment_head = NaN; sediment_compressibility = 0; @@ -11,17 +20,10 @@ sediment_thickness = 0; sediment_transmitivity = 0; - water_compressibility = 0; - isefficientlayer = 0; - penalty_factor = 0; - sedimentlimit_flag = 0; - sedimentlimit = 0; - transfer_flag = 0; - leakage_factor = 0; - spcepl_head = NaN; - epl_compressibility = 0; - epl_porosity = 0; - epl_thickness = 0; - epl_transmitivity = 0; + spcepl_head = NaN; + epl_compressibility = 0; + epl_porosity = 0; + epl_thickness = 0; + epl_transmitivity = 0; end @@ -38,11 +40,8 @@ %Parameters from de Fleurian 2013 - obj.sediment_compressibility = 1.0e-08; - obj.sediment_porosity = .4; - obj.sediment_thickness = 20.0; - obj.sediment_transmitivity = 8.0e-04; obj.water_compressibility = 5.04e-10; obj.isefficientlayer = 1; obj.penalty_factor = 3; + obj.rel_tol = 1.0e-06; obj.sedimentlimit_flag = 0; obj.sedimentlimit = 0; @@ -50,8 +49,13 @@ obj.leakage_factor = 10.0; - obj.epl_compressibility = 1.0e-08; - obj.epl_porosity = 0.4; - obj.epl_thickness = 1.0; - obj.epl_transmitivity = 8.0e-02; + obj.sediment_compressibility = 1.0e-08; + obj.sediment_porosity = .4; + obj.sediment_thickness = 20.0; + obj.sediment_transmitivity = 8.0e-04; + + obj.epl_compressibility = 1.0e-08; + obj.epl_porosity = 0.4; + obj.epl_thickness = 1.0; + obj.epl_transmitivity = 8.0e-02; end % }}} @@ -63,4 +67,17 @@ end + md = checkfield(md,'hydrology.water_compressibility','>',0,'numel',1); + md = checkfield(md,'hydrology.isefficientlayer','numel',[1],'values',[0 1]); + md = checkfield(md,'hydrology.penalty_factor','>',0,'numel',1); + md = checkfield(md,'hydrology.rel_tol','>',0,'numel',1); + md = checkfield(md,'hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]); + md = checkfield(md,'hydrology.transfer_flag','numel',[1],'values',[0 1]); + if obj.sedimentlimit_flag==1, + md = checkfield(md,'hydrology.sedimentlimit','>',0,'numel',1); + end + if obj.transfer_flag==1, + md = checkfield(md,'hydrology.leakage_factor','>',0,'numel',1); + end + md = checkfield(md,'hydrology.spcsediment_head','forcing',1); md = checkfield(md,'hydrology.sediment_compressibility','>',0,'numel',1); @@ -68,20 +85,6 @@ md = checkfield(md,'hydrology.sediment_thickness','>',0,'numel',1); md = checkfield(md,'hydrology.sediment_transmitivity','>',0,'numel',1); - md = checkfield(md,'hydrology.water_compressibility','>',0,'numel',1); - md = checkfield(md,'hydrology.isefficientlayer','numel',[1],'values',[0 1]); - md = checkfield(md,'hydrology.penalty_factor','>',0,'numel',1); - md = checkfield(md,'hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]); - md = checkfield(md,'hydrology.transfer_flag','numel',[1],'values',[0 1]); - - if obj.sedimentlimit_flag==1, - md = checkfield(md,'hydrology.sedimentlimit','>',0,'numel',1); - end - - if obj.transfer_flag==1, - md = checkfield(md,'hydrology.leakage_factor','>',0,'numel',1); - end if obj.isefficientlayer==1, - md = checkfield(md,'hydrology.spcepl_head','forcing',1); md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1); @@ -93,14 +96,9 @@ function disp(obj) % {{{ disp(sprintf(' hydrology Dual Porous Continuum Equivalent parameters:')); - disp(sprintf(' - for the sediment layer')); - - fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]'); - fielddisplay(obj,'sediment_compressibility','sediment compressibility [Pa^-1]'); - fielddisplay(obj,'sediment_porosity','sediment [dimensionless]'); - fielddisplay(obj,'sediment_thickness','sediment thickness [m]'); - fielddisplay(obj,'sediment_transmitivity','sediment transmitivity [m^2/s]'); + disp(sprintf(' - general parameters')); fielddisplay(obj,'water_compressibility','compressibility of water [Pa^-1]'); fielddisplay(obj,'isefficientlayer','do we use an efficient drainage system [1: true; 0: false]'); - fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method'); + fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]'); + fielddisplay(obj,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]'); fielddisplay(obj,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer'); disp(sprintf('%55s 0: no limit',' ')); @@ -117,8 +115,13 @@ fielddisplay(obj,'leakage_factor','user defined leakage factor [m]'); end - + disp(sprintf(' - for the sediment layer')); + fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]'); + fielddisplay(obj,'sediment_compressibility','sediment compressibility [Pa^-1]'); + fielddisplay(obj,'sediment_porosity','sediment [dimensionless]'); + fielddisplay(obj,'sediment_thickness','sediment thickness [m]'); + fielddisplay(obj,'sediment_transmitivity','sediment transmitivity [m^2/s]'); + if obj.isefficientlayer==1, disp(sprintf(' - for the epl layer')); - fielddisplay(obj,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]'); fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]'); @@ -130,26 +133,25 @@ end % }}} function marshall(obj,fid) % {{{ - WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer'); + WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer'); + WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double'); + WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean'); + WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'fieldname','rel_tol','format','Double'); + WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer'); + WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer'); + if obj.sedimentlimit_flag==1, + WriteData(fid,'object',obj,'fieldname','sedimentlimit','format','Double'); + end + if obj.transfer_flag==1, + WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double'); + end + WriteData(fid,'object',obj,'fieldname','spcsediment_head','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','sediment_compressibility','format','Double'); WriteData(fid,'object',obj,'fieldname','sediment_porosity','format','Double'); WriteData(fid,'object',obj,'fieldname','sediment_thickness','format','Double'); - WriteData(fid,'object',obj,'fieldname','sediment_transmitivity','format','Double'); - WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double'); - WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean'); - WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); - WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer'); - WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer'); - - if obj.sedimentlimit_flag==1, - WriteData(fid,'object',obj,'fieldname','sedimentlimit','format','Double'); - end - - if obj.transfer_flag==1, - WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double'); - end - - if obj.isefficientlayer==1, - + WriteData(fid,'object',obj,'fieldname','sediment_transmitivity','format','Double'); + + if obj.isefficientlayer==1, WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double'); Index: /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m (revision 15061) +++ /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m (revision 15062) @@ -1,4 +1,5 @@ function writeVTKcell(filename,model,Solution) % vtk export +% function writeVTKcell(filename,model,Solution) % creates a vtk-file filename.vtk containing simplicial mesh data % (only work for triangle now) Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 15061) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 15062) @@ -1143,4 +1143,18 @@ return StringToEnum('EplHead')[0] +def HydrologydcRelTolEnum(): + """ + HYDROLOGYDCRELTOLENUM - Enum of HydrologydcRelTol + + 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=HydrologydcRelTolEnum() + """ + + return StringToEnum('HydrologydcRelTol')[0] + def HydrologydcSpcsedimentHeadEnum(): """ @@ -7777,4 +7791,4 @@ """ - return 554 - + return 555 + Index: /issm/trunk-jpl/src/m/enum/HydrologydcRelTolEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologydcRelTolEnum.m (revision 15062) +++ /issm/trunk-jpl/src/m/enum/HydrologydcRelTolEnum.m (revision 15062) @@ -0,0 +1,11 @@ +function macro=HydrologydcRelTolEnum() +%HYDROLOGYDCRELTOLENUM - Enum of HydrologydcRelTol +% +% 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=HydrologydcRelTolEnum() + +macro=StringToEnum('HydrologydcRelTol'); Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 15061) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 15062) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=554; +macro=555;