Index: sm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_singleloop.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_singleloop.cpp (revision 27189) +++ (revision ) @@ -1,271 +1,0 @@ -/*!\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; -} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 27189) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 27190) @@ -14,5 +14,4 @@ void solutionsequence_thermal_nonlinear(FemModel* femmodel); void solutionsequence_hydro_nonlinear(FemModel* femmodel, bool* conv_fail); -void solutionsequence_hydro_singleloop(FemModel* femmodel, bool* conv_fail); void solutionsequence_shakti_nonlinear(FemModel* femmodel); void solutionsequence_glads_nonlinear(FemModel* femmodel);