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

Last change on this file since 13325 was 13325, checked in by Eric.Larour, 13 years ago

NEW: large change. Replaced all our instances of ExternalResult objects with one
generic templated class. Much simpler, though it hit some of the code in the solutions
in modules.
File size: 4.5 KB
Line 
1/*!\file: solver_nonlinear.cpp
2 * \brief: core of a non-linear solution, using fixed-point method
3 */
4
5#include "../toolkits/toolkits.h"
6#include "../classes/objects/objects.h"
7#include "../io/io.h"
8#include "../EnumDefinitions/EnumDefinitions.h"
9#include "../modules/modules.h"
10#include "../solutions/solutions.h"
11#include "./solvers.h"
12
13void solver_newton(FemModel* femmodel){
14
15 /*intermediary: */
16 bool converged;
17 int num_unstable_constraints;
18 int count;
19 IssmDouble kmax;
20 Matrix<IssmDouble>* Kff = NULL;
21 Matrix<IssmDouble>* Kfs = NULL;
22 Matrix<IssmDouble>* Jff = NULL;
23 Vector<IssmDouble>* ug = NULL;
24 Vector<IssmDouble>* old_ug = NULL;
25 Vector<IssmDouble>* uf = NULL;
26 Vector<IssmDouble>* old_uf = NULL;
27 Vector<IssmDouble>* duf = NULL;
28 Vector<IssmDouble>* pf = NULL;
29 Vector<IssmDouble>* pJf = NULL;
30 Vector<IssmDouble>* df = NULL;
31 Vector<IssmDouble>* ys = NULL;
32
33 /*parameters:*/
34 int max_nonlinear_iterations;
35 int configuration_type;
36
37 /*Recover parameters: */
38 femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
39 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
40 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
41
42 count=1;
43 converged=false;
44
45 /*Start non-linear iteration using input velocity: */
46 GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
47 Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters);
48
49 //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
50 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
51 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
52
53 for(;;){
54
55 xdelete(&old_ug);old_ug=ug;
56 xdelete(&old_uf);old_uf=uf;
57
58 /*Solver forward model*/
59 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
60 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
61 Reduceloadx(pf,Kfs,ys);xdelete(&Kfs);
62 Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df);
63 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
64 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug);
65
66 /*Check convergence*/
67 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
68 xdelete(&Kff); xdelete(&pf);
69 if(converged==true){
70 bool max_iteration_state=false;
71 int tempStep=1;
72 IssmDouble tempTime=1.0;
73 femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
74 break;
75 }
76 if(count>=max_nonlinear_iterations){
77 _pprintLine_(" maximum number of Newton iterations (" << max_nonlinear_iterations << ") exceeded");
78 bool max_iteration_state=true;
79 int tempStep=1;
80 IssmDouble tempTime=1.0;
81 femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
82 break;
83 }
84
85 /*Prepare next iteration using Newton's method*/
86 SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
87 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
88 Reduceloadx(pf,Kfs,ys); xdelete(&Kfs);
89
90 pJf=pf->Duplicate(); Kff->MatMult(uf,pJf); xdelete(&Kff);
91 pJf->Scale(-1.0); pJf->AXPY(pf,+1.0); xdelete(&pf);
92
93 CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
94 Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf);
95 uf->AXPY(duf, 1.0); xdelete(&duf);
96 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
97 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
98
99 count++;
100 }
101
102 if(VerboseConvergence()) _pprintLine_("\n total number of iterations: " << count-1);
103
104 /*clean-up*/
105 xdelete(&uf);
106 xdelete(&ug);
107 xdelete(&old_ug);
108 xdelete(&old_uf);
109}
Note: See TracBrowser for help on using the repository browser.