source: issm/trunk-jpl/src/c/solvers/solver_newton.cpp@ 14891

/*!\file: solver_nonlinear.cpp
 * \brief: core of a non-linear solution, using fixed-point method 
 */ 

#include "../toolkits/toolkits.h"
#include "../classes/objects/objects.h"
#include "../io/io.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../modules/modules.h"
#include "../solutions/solutions.h"
#include "./solvers.h"

void solver_newton(FemModel* femmodel){

        /*intermediary: */
        bool   converged;
        int    count,newton;
        IssmDouble kmax;
        Matrix<IssmDouble>* Kff = NULL;
        Matrix<IssmDouble>* Kfs    = NULL;
        Matrix<IssmDouble>* Jff = NULL;
        Vector<IssmDouble>* ug  = NULL;
        Vector<IssmDouble>* old_ug = NULL;
        Vector<IssmDouble>* uf  = NULL;
        Vector<IssmDouble>* old_uf = NULL;
        Vector<IssmDouble>* duf = NULL;
        Vector<IssmDouble>* pf  = NULL;
        Vector<IssmDouble>* pJf    = NULL;
        Vector<IssmDouble>* df  = NULL;
        Vector<IssmDouble>* ys  = NULL;

        /*parameters:*/
        int max_nonlinear_iterations;
        int  configuration_type;

        /*Recover parameters: */
        femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
        femmodel->parameters->FindParam(&newton,DiagnosticIsnewtonEnum);
        femmodel->UpdateConstraintsx();

        count=1;
        converged=false;

        /*Start non-linear iteration using input velocity: */
        GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
        Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters);

        //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
        InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);

        for(;;){

                delete old_ug;old_ug=ug;
                delete old_uf;old_uf=uf;

                /*Solver forward model*/
                if(count==1 || newton==2){
                        femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
                        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
                        Reduceloadx(pf,Kfs,ys);delete Kfs;
                        Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);delete df;
                        Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);delete ug;
                        delete old_ug;old_ug=ug;
                        delete old_uf;old_uf=uf;
                }
                uf=old_uf->Duplicate(); old_uf->Copy(uf);

                /*Prepare next iteration using Newton's method*/
                femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);delete df;
                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
                Reduceloadx(pf,Kfs,ys);delete Kfs;

                pJf=pf->Duplicate();
                Kff->MatMult(uf,pJf);// delete Kff);
                pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);     //delete pf);

                femmodel->CreateJacobianMatrixx(&Jff,kmax);
                Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); delete Jff; delete pJf;
                uf->AXPY(duf, 1.0); delete duf;
                Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);

                /*Check convergence*/
                convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); 
                delete Kff; delete pf;
                if(converged==true){    
                        bool max_iteration_state=false;
                        int tempStep=1;
                        IssmDouble tempTime=1.0;
                        femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
                        break;
                }
                if(count>=max_nonlinear_iterations){
                        _pprintLine_("   maximum number of Newton iterations (" << max_nonlinear_iterations << ") exceeded"); 
                        bool max_iteration_state=true;
                        int tempStep=1;
                        IssmDouble tempTime=1.0;
                        femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
                        break;
                }

                count++;
        }

        if(VerboseConvergence()) _pprintLine_("\n   total number of iterations: " << count-1);

        /*clean-up*/
        delete uf;
        delete ug;
        delete old_ug;
        delete old_uf;
}