Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_singleloop.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_singleloop.cpp (revision 27166) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_singleloop.cpp (revision 27166) @@ -0,0 +1,271 @@ +/*!\file: solutionsequence_hydro_singleloop.cpp + * \brief: core of the hydro solution + trying a single loop set-up to see the effect on convergence + */ + +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +/*FIXME, dirty hack to get the solutionsequence linear needed to compute the slopes*/ +#include "../solutionsequences/solutionsequences.h" + +void solutionsequence_hydro_nonlinear(FemModel* femmodel, bool* pconv_fail){ + /*solution : */ + Vector* ug_sed=NULL; + Vector* uf_sed=NULL; + Vector* uf_sed_sub_iter=NULL; + Vector* ug_sed_sub_iter=NULL; + Vector* ug_sed_main_iter=NULL; + Vector* ug_sed_init=NULL; + + Vector* ug_epl=NULL; + Vector* uf_epl=NULL; + Vector* uf_epl_sub_iter=NULL; + Vector* ug_epl_sub_iter=NULL; + Vector* ug_epl_main_iter=NULL; + Vector* ug_epl_init=NULL; + + Vector* ys=NULL; + Vector* dug=NULL; + Vector* duf=NULL; + + Matrix* Kff=NULL; + Matrix* Kfs=NULL; + + Vector* pf=NULL; + Vector* df=NULL; + + HydrologyDCInefficientAnalysis* inefanalysis = NULL; + HydrologyDCEfficientAnalysis* effanalysis = NULL; + + bool sedconverged,eplconverged,hydroconverged, headconverged, newthick; + bool isefficientlayer; + bool sliceadapt; + int constraints_converged; + int num_unstable_constraints; + int sedcount,eplcount,hydrocount; + int hydro_maxiter; + IssmDouble sediment_kmax; + IssmDouble eps_hyd; + IssmDouble ndu_sed,nu_sed; + IssmDouble ndu_epl,nu_epl; + IssmDouble ThickCount,L2Count; + /*Recover parameters: */ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + femmodel->parameters->FindParam(&sliceadapt,HydrologyStepAdaptEnum); + femmodel->parameters->FindParam(&hydro_maxiter,HydrologydcMaxIterEnum); + femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); + hydrocount=1; + hydroconverged=false; + + /*{{{ Retrieve inputs as the initial state for the non linear iteration*/ + GetBasalSolutionFromInputsx(&ug_sed,femmodel); + ug_sed_init=ug_sed->Duplicate(); + ug_sed->Copy(ug_sed_init); + Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters); + inefanalysis = new HydrologyDCInefficientAnalysis(); + /*Initialize the IDS element mask to exclude frozen nodes*/ + inefanalysis->ElementizeIdsMask(femmodel); + + if(isefficientlayer) { + effanalysis = new HydrologyDCEfficientAnalysis(); + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + GetBasalSolutionFromInputsx(&ug_epl,femmodel); + ug_epl_init=ug_epl->Duplicate(); + ug_epl->Copy(ug_epl_init); + /*Initialize the EPL element mask*/ + inefanalysis->ElementizeEplMask(femmodel); + effanalysis->InitZigZagCounter(femmodel); + } + /*}}}*/ + /*The real computation starts here,single loop on the two layer system*/ + for(;;){ + sedcount=1; + eplcount=1; + + /*keep the initial vlaue to compute convergence (on ug as activation can change for both layers)*/ + ug_sed_main_iter=ug_sed->Duplicate(); + ug_sed->Copy(ug_sed_main_iter); + if(isefficientlayer){ + 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); + ResetZigzagCounterx(femmodel); + InputUpdateFromConstantx(femmodel,false,ConvergedEnum); + femmodel->UpdateConstraintsx(); + + /*Reset constraint on the ZigZag Lock*/ + ResetConstraintsx(femmodel); + + /*{{{ Loop for the sediment penalization*/ + femmodel->profiler->Start(SEDLOOP); + + for(;;){ + /*{{{ Core of the computation*/ + if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n"); + femmodel->profiler->Start(SEDMatrix); + sedconverged=false; + SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes); + Reduceloadx(pf,Kfs,ys); delete Kfs; + delete uf_sed; + femmodel->profiler->Stop(SEDMatrix); + + femmodel->profiler->Start(SOLVER); + Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); + delete Kff; delete pf; delete df; delete uf_sed_sub_iter; + uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter); + uf_sed->Copy(uf_sed_sub_iter); + femmodel->profiler->Stop(SOLVER); + + delete ug_sed; + femmodel->profiler->Start(SEDUpdate); + Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel,ug_sed); + ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); + femmodel->profiler->Stop(SEDUpdate); + + /*}}}*/ + if (!sedconverged){ + /*First check that all the penalizations are applied*/ + 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; + } + + if(isefficientlayer){ + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + /*updating mask*/ + if(VerboseSolution()) _printf0_("==updating mask...\n"); + femmodel->HydrologyEPLupdateDomainx(&ThickCount); + ResetZigzagCounterx(femmodel); + InputUpdateFromConstantx(femmodel,false,ConvergedEnum); + + /*{{{ Retrieve the EPL head slopes and compute EPL Thickness*/ + if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); + + femmodel->profiler->Start(EPLMasking); + femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); + femmodel->UpdateConstraintsL2ProjectionEPLx(&L2Count); + femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); + solutionsequence_linear(femmodel); + femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + effanalysis->ComputeEPLThickness(femmodel); + //updating mask after the computation of the epl thickness (Allow to close too thin EPL) + femmodel->HydrologyEPLupdateDomainx(&ThickCount); + /*}}}*/ + femmodel->profiler->Stop(EPLMasking); + + femmodel->profiler->Start(EPLMatrices); + if(VerboseSolution()) _printf0_("Building EPL Matrix...\n"); + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes); + Reduceloadx(pf,Kfs,ys); delete Kfs; + delete uf_epl; + femmodel->profiler->Stop(EPLMatrices); + + femmodel->profiler->Start(SOLVER); + Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); + femmodel->profiler->Stop(SOLVER); + + delete Kff; delete pf; delete df; + delete uf_epl_sub_iter; + uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter); + uf_epl->Copy(uf_epl_sub_iter); + delete ug_epl; + + femmodel->profiler->Start(EPLUpdate); + Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel,ug_epl); + ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); + femmodel->profiler->Stop(EPLUpdate); + } + if(!hydroconverged){ + //compute norm(du)/norm(u) + dug=ug_sed_main_iter->Duplicate(); _assert_(dug); + ug_sed_main_iter->Copy(dug); + dug->AYPX(ug_sed,-1.0); + ndu_sed=dug->Norm(NORM_TWO); + delete dug; + nu_sed=ug_sed_main_iter->Norm(NORM_TWO); + delete ug_sed_main_iter; + if (xIsNan(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*/ + + if(isefficientlayer){ + 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 ((ndu_epl/nu_epl)=hydro_maxiter*9/10 && sliceadapt && !hydroconverged) { + if(VerboseSolution()) _printf0_("hydrology main loop did not converged after "<SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + InputUpdateFromSolutionx(femmodel,ug_sed_init); + delete ug_sed_init; + break; + } + if (hydrocount>=hydro_maxiter){ + _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); + } + } + hydrocount++; + if(hydroconverged)break; + + if(*pconv_fail)break; + } + /*To deal with adaptative stepping we only save results if we are actually converged*/ + if(hydroconverged){ + if(isefficientlayer){ + InputUpdateFromSolutionx(femmodel,ug_epl); + delete ug_epl_init; + } + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + InputUpdateFromSolutionx(femmodel,ug_sed); + delete ug_sed_init; + } + /*Free resources: */ + delete ug_epl; + delete ug_sed; + delete uf_sed; + delete uf_epl; + delete uf_epl_sub_iter; + delete uf_sed_sub_iter; + delete inefanalysis; + delete effanalysis; +}