Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17005) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17006) @@ -4659,5 +4659,5 @@ bool isefficientlayer; bool active_element; - int transfermethod; + int transfermethod,step; IssmDouble sed_trans,sed_thick; IssmDouble leakage,h_max; @@ -4666,11 +4666,15 @@ IssmDouble epl_specificstoring[numdof],sedstoring[numdof]; IssmDouble epl_head[numdof],sed_head[numdof]; + IssmDouble old_transfer[numdof]; Input* active_element_input=NULL; + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get the flag to know if the efficient layer is present*/ this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + + this->parameters->FindParam(&step,StepEnum); if(isefficientlayer){ @@ -4685,4 +4689,5 @@ case 1: + active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); active_element_input->GetInputValue(&active_element); @@ -4698,37 +4703,39 @@ if(!active_element){ - /*No transfer if the EPL is not active*/ - for(int i=0;iSetValue(doflist[i],wh_trans,INS_VAL); - } + /* for(int i=0;iSetValue(doflist[i],wh_trans,INS_VAL); */ + /* } */ } else{ - + + //GetInputListOnVertices(&old_transfer[0],WaterTransferEnum); + for(int i=0;iGetEplSpecificStoring(); sedstoring[i]=matpar->GetSedimentStoring(); - + this->GetHydrologyDCInefficientHmax(&h_max,i); + /*EPL head higher than sediment head, transfer from the epl to the sediment*/ if(epl_head[i]>sed_head[i]){ 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,i); - if(sed_head[i]>=h_max)wh_trans=0.0; + if(sed_head[i]>=h_max){ + wh_trans=0.0; + } } /*EPL head lower than sediment head, transfer from the sediment to the epl*/ else if(epl_head[i]<=sed_head[i]){ - wh_trans=sedstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=sedstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick); } + + /*Introduce relaxation*/ + //wh_trans=old_transfer[i]+0.9*(wh_trans-old_transfer[i]); + /*Assign output pointer*/ transfer->SetValue(doflist[i],wh_trans,INS_VAL); - /* if(nodes[i]->id>=54){ */ - /* printf("%i %e %e %e \n",nodes[i]->id-54,wh_trans,sed_head[i],epl_head[i]); */ - /* } */ - /* else{*/ - /* printf("%i %e %e %e \n",nodes[i]->id,wh_trans,sed_head[i],epl_head[i]); */ } } @@ -4799,7 +4806,7 @@ } /*If epl thickness gets under 0, close the layer*/ - /* else if(epl_thickness[i]<0.0){ */ - /* vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); */ - /* } */ + else if(epl_thickness[i]<0.0){ + vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); + } /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ this->GetHydrologyDCInefficientHmax(&h_max,i); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 17005) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 17006) @@ -382,12 +382,12 @@ /*FUNCTION Matpar::GetSedimentStoring {{{*/ IssmDouble Matpar::GetSedimentStoring(){ - return this->rho_freshwater* this->g* this->sediment_porosity* this->sediment_thickness* - ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity)); + return this->rho_freshwater * this->g * this->sediment_porosity * this->sediment_thickness * + (this->water_compressibility + (this->sediment_compressibility / this->sediment_porosity)); } /*}}}*/ /*FUNCTION Matpar::GetEplSpecificStoring {{{*/ IssmDouble Matpar::GetEplSpecificStoring(){ - return this->rho_freshwater* this->g* this->epl_porosity* - ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity)); + return this->rho_freshwater * this->g * this->epl_porosity * + (this->water_compressibility + (this->epl_compressibility / this->epl_porosity)); } /*}}}*/ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17005) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17006) @@ -80,14 +80,7 @@ ug_sed->Copy(ug_sed_main_iter); - //test - /* uf_sed_sub_iter=uf_sed->Duplicate(); */ - /* uf_sed->Copy(uf_sed_sub_iter); */ - if(isefficientlayer){ ug_epl_main_iter=ug_epl->Duplicate(); ug_epl->Copy(ug_epl_main_iter); - //test - ug_epl_sub_iter=ug_epl->Duplicate(); - ug_epl->Copy(ug_epl_sub_iter); } @@ -104,4 +97,5 @@ sedconverged=false; for(;;){ + uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter); uf_sed->Copy(uf_sed_sub_iter); @@ -121,4 +115,5 @@ if(num_unstable_constraints==0) sedconverged = true; if (sedcount>=hydro_maxiter){ + //sedconverged=true; _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); } @@ -128,5 +123,13 @@ if(sedconverged){ - sedconverged=false; + if(isefficientlayer){ + /*Updating Elemental Mask*/ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + femmodel->HydrologyTransferx(); + } + sedconverged=false; + /*Checking convegence on the value of the sediment head*/ duf=uf_sed_sub_iter->Duplicate();_assert_(duf); uf_sed_sub_iter->Copy(duf); @@ -137,4 +140,5 @@ if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("convergence criterion is NaN!"); if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the EPL is used but empty*/ + if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); if((ndu_sed/nu_sed)ElementizeEplMask(femmodel); - delete analysis; - femmodel->HydrologyTransferx(); - } femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum); @@ -166,5 +163,4 @@ InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); InputUpdateFromConstantx(femmodel,false,ConvergedEnum); - femmodel->HydrologyEPLupdateDomainx(); femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); @@ -188,7 +184,7 @@ //updating mask after the computation of the epl thickness (Allow to close too thin EPL) femmodel->HydrologyEPLupdateDomainx(); + ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); - ug_epl->Copy(ug_epl_sub_iter);_assert_(ug_epl); - + ug_epl->Copy(ug_epl_sub_iter); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); @@ -204,5 +200,5 @@ InputUpdateFromSolutionx(femmodel,ug_epl); ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - femmodel->HydrologyEPLupdateDomainx(); + if (!eplconverged){ @@ -210,6 +206,5 @@ if(num_unstable_constraints==0) eplconverged = true; if (eplcount>=hydro_maxiter){ - /*Hacking to get the results of non converged runs*/ - //eplconverged = true; + //eplconverged =true; _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); } @@ -220,4 +215,5 @@ if(eplconverged){ eplconverged=false; + dug=ug_epl_sub_iter->Duplicate();_assert_(dug); ug_epl_sub_iter->Copy(dug); @@ -229,4 +225,5 @@ if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("convergence criterion is NaN!"); if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ + if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); if((ndu_epl/nu_epl)ElementizeEplMask(femmodel); delete analysis; femmodel->HydrologyTransferx(); - + InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); InputUpdateFromConstantx(femmodel,sediment_kmax,MeltingOffsetEnum); @@ -248,5 +244,4 @@ } } - /*System convergence check*/ @@ -303,5 +298,5 @@ hydrocount++; if(hydroconverged)break; - } +} InputUpdateFromSolutionx(femmodel,ug_sed);