Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17023) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17024) @@ -128,6 +128,11 @@ /*Check that all nodes are active, else return empty matrix*/ - if(!basalelement->AllActive()) return NULL; - + if(!basalelement->AllActive()) { + if(meshtype!=Mesh2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + return NULL; + } /* Intermediaries */ IssmDouble D_scalar,Jdet,dt; @@ -186,4 +191,5 @@ xDelete(B); delete gauss; + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; @@ -209,6 +215,11 @@ /*Check that all nodes are active, else return empty matrix*/ - if(!basalelement->AllActive()) return NULL; - + if(!basalelement->AllActive()){ + if(meshtype!=Mesh2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + return NULL; + } /*Intermediaries */ IssmDouble dt,scalar,water_head,connectivity; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17023) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17024) @@ -23,13 +23,11 @@ Vector* ug_epl_main_iter=NULL; - Vector* ys=NULL; Vector* dug=NULL; - - //testing stuff Vector* duf=NULL; Matrix* Kff=NULL; Matrix* Kfs=NULL; + Vector* pf=NULL; Vector* df=NULL; @@ -51,8 +49,12 @@ femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); femmodel->parameters->FindParam(&time,TimeEnum); + /*FIXME, hardcoded, put on an enum*/ hydro_maxiter=150; + hydrocount=1; hydroconverged=false; + /*We don't need the outer loop if only one layer is used*/ + if(!isefficientlayer) hydroconverged=true; /*Retrieve inputs as the initial state for the non linear iteration*/ @@ -63,6 +65,6 @@ femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); GetSolutionFromInputsx(&ug_epl,femmodel); - - /*Initialize the transfer input*/ + + /*Initialize the element mask*/ HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); analysis->ElementizeEplMask(femmodel); @@ -72,17 +74,17 @@ femmodel->HydrologyTransferx(); - /*Iteration on the two layers*/ + /*The real computation starts here, outermost loop is on the two layer system*/ for(;;){ sedcount=1; eplcount=1; - //save pointer to old velocity - ug_sed_main_iter=ug_sed->Duplicate(); - ug_sed->Copy(ug_sed_main_iter); - + + /*If there is two layers we need an outer loop value to compute convergence*/ if(isefficientlayer){ + ug_sed_main_iter=ug_sed->Duplicate(); + ug_sed->Copy(ug_sed_main_iter); ug_epl_main_iter=ug_epl->Duplicate(); ug_epl->Copy(ug_epl_main_iter); } - + /*Loop on sediment layer to deal with transfer and head value*/ femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); @@ -90,69 +92,62 @@ femmodel->UpdateConstraintsx(); femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); - - /*Reset constraint on the ZigZag Lock*/ + + /*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/ ResetConstraintsx(femmodel); - - /*Iteration on the sediment layer*/ sedconverged=false; + /* {{{ *//*Treating the sediment*/ for(;;){ - /* if(isefficientlayer){ */ - /* /\*Updating Elemental Mask*\/ */ - /* HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); */ - /* analysis->ElementizeEplMask(femmodel); */ - /* delete analysis; */ - /* femmodel->HydrologyTransferx(); */ - /* } */ uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter); uf_sed->Copy(uf_sed_sub_iter); - SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); - CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); - Reduceloadx(pf,Kfs,ys); delete Kfs; - delete uf_sed; - Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); - delete Kff; delete pf; delete df; - delete ug_sed; - Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; - InputUpdateFromSolutionx(femmodel,ug_sed); - ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - - if (!sedconverged){ - if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n"); - if(num_unstable_constraints==0) sedconverged = true; - if (sedcount>=hydro_maxiter){ - //sedconverged=true; - _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); - } - } - - sedcount++; - + /* {{{ *//*Loop on the sediment layer to deal with the penalization*/ + for(;;){ + /* {{{ *//*Core of the computation*/ + SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); delete Kfs; + delete uf_sed; + Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); + delete Kff; delete pf; delete df; + delete ug_sed; + Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel,ug_sed); + ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); + /* }}} */ + if (!sedconverged){ + if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n"); + if(num_unstable_constraints==0) sedconverged = true; + if (sedcount>=hydro_maxiter){ + _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); + } + } + /*Add an iteration and get out of the loop if the penalisation is converged*/ + sedcount++; + if(sedconverged)break; + } + /* }}} *//*End of the sediment penalization loop*/ + sedconverged=false; + /*Update Elemental Mask and compute transfer*/ + if(isefficientlayer){ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + femmodel->HydrologyTransferx(); + } + /*Checking convegence on the value of the sediment head*/ + duf=uf_sed_sub_iter->Duplicate();_assert_(duf); + uf_sed_sub_iter->Copy(duf); + duf->AYPX(uf_sed,-1.0); + ndu_sed=duf->Norm(NORM_TWO); + delete duf; + nu_sed=uf_sed_sub_iter->Norm(NORM_TWO); + 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 layer is 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)Duplicate();_assert_(duf); - uf_sed_sub_iter->Copy(duf); - duf->AYPX(uf_sed,-1.0); - ndu_sed=duf->Norm(NORM_TWO); - delete duf; - nu_sed=uf_sed_sub_iter->Norm(NORM_TWO); - 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); @@ -162,22 +157,23 @@ } } - - /*Second layer*/ + /* }}} *//*End of the global sediment loop*/ + + /* {{{ *//*Now dealing with the EPL in the same way*/ if(isefficientlayer){ - + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + /*updating mask*/ + femmodel->HydrologyEPLupdateDomainx(); InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); InputUpdateFromConstantx(femmodel,false,ConvergedEnum); femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); + + eplconverged = false; - - - /*Iteration on the EPL layer*/ - eplconverged = false; for(;;){ - - //updating mask after the computation of the epl thickness (Allow to close too thin EPL) - femmodel->HydrologyEPLupdateDomainx(); - /*Start by retrieving the EPL head slopes and compute EPL Thickness*/ + ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); + ug_epl->Copy(ug_epl_sub_iter); + + /* {{{ *//*Start by retrieving the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); @@ -188,60 +184,58 @@ solutionsequence_linear(femmodel); femmodel->HydrologyEPLThicknessx(); + //updating mask after the computation of the epl thickness (Allow to close too thin EPL) + femmodel->HydrologyEPLupdateDomainx(); + /* }}} */ + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - - /*Updating Elemental Mask*/ + /* {{{ *//*Loop on the epl layer to deal with the penalization*/ + for(;;){ + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); delete Kfs; + delete uf_epl; + Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); + delete Kff; delete pf; delete df; + delete uf_epl_sub_iter; + uf_epl_sub_iter=uf_epl->Duplicate(); + uf_epl->Copy(uf_epl_sub_iter); + delete ug_epl; + Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel,ug_epl); + ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); + + if (!eplconverged){ + if(VerboseConvergence()) _printf0_(" # EPL unstable constraints = " << num_unstable_constraints << "\n"); + if(num_unstable_constraints==0) eplconverged = true; + if (eplcount>=hydro_maxiter){ + _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); + } + } + /*Add an iteration and get out of the loop if the penalisation is converged*/ + eplcount++; + if(eplconverged) break; + } + /* }}} */ /*End of the EPL penalization loop*/ + + eplconverged=false; + /*Update Elemental Mask and compute transfer*/ HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); analysis->ElementizeEplMask(femmodel); delete analysis; femmodel->HydrologyTransferx(); - - ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); - ug_epl->Copy(ug_epl_sub_iter); - - SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); - CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); - Reduceloadx(pf,Kfs,ys); delete Kfs; - delete uf_epl; - Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); - delete Kff; delete pf; delete df; - delete uf_epl_sub_iter; - uf_epl_sub_iter=uf_epl->Duplicate(); - uf_epl->Copy(uf_epl_sub_iter); - delete ug_epl; - Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; - InputUpdateFromSolutionx(femmodel,ug_epl); - ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - - - if (!eplconverged){ - if(VerboseConvergence()) _printf0_(" # EPL unstable constraints = " << num_unstable_constraints << "\n"); - if(num_unstable_constraints==0) eplconverged = true; - if (eplcount>=hydro_maxiter){ - //eplconverged =true; - _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); - } - } - eplcount++; - + + dug=ug_epl_sub_iter->Duplicate();_assert_(dug); + ug_epl_sub_iter->Copy(dug); + dug->AYPX(ug_epl,-1.0); + ndu_epl=dug->Norm(NORM_TWO); + delete dug; + nu_epl=ug_epl_sub_iter->Norm(NORM_TWO); + 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)Duplicate();_assert_(dug); - ug_epl_sub_iter->Copy(dug); - dug->AYPX(ug_epl,-1.0); - ndu_epl=dug->Norm(NORM_TWO); - delete dug; - nu_epl=ug_epl_sub_iter->Norm(NORM_TWO); - - 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)(ndu_sed) || xIsNan(nu_sed)) _error_("Sed convergence criterion is NaN!"); if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the Sediment is used but empty*/ + dug=ug_epl_main_iter->Duplicate();_assert_(dug); + ug_epl_main_iter->Copy(dug); + dug->AYPX(ug_epl,-1.0); + ndu_epl=dug->Norm(NORM_TWO); + delete dug; + nu_epl=ug_epl_main_iter->Norm(NORM_TWO); + delete ug_epl_main_iter; + if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("EPL 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 (!xIsNan(eps_hyd)){ - if (!isefficientlayer){ - if ((ndu_sed/nu_sed) " << eps_hyd*100 << " %\n"); - hydroconverged=false; - } - } - else{ - dug=ug_epl_main_iter->Duplicate();_assert_(dug); - ug_epl_main_iter->Copy(dug); - dug->AYPX(ug_epl,-1.0); - ndu_epl=dug->Norm(NORM_TWO); - delete dug; - nu_epl=ug_epl_main_iter->Norm(NORM_TWO); - delete ug_epl_main_iter; - if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("EPL 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 ((ndu_epl/nu_epl)=hydro_maxiter){ - /*Hacking to get the results of non converged runs*/ - //hydroconverged = true; - _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); + _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); } } hydrocount++; if(hydroconverged)break; -} - + } + /* }}} */ InputUpdateFromSolutionx(femmodel,ug_sed); if(isefficientlayer)InputUpdateFromSolutionx(femmodel,ug_epl); - /*Free ressources: */ delete ug_epl; @@ -316,4 +296,4 @@ delete uf_epl; delete uf_epl_sub_iter; - + }