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

Last change on this file since 11337 was 11337, checked in by Mathieu Morlighem, 13 years ago

Fixed Newton's method: do not change ys
File size: 3.7 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 "../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 double kmax;
20 Mat Kff = NULL, Kfs = NULL, Jff = NULL;
21 Vec ug = NULL, old_ug = NULL;
22 Vec uf = NULL, old_uf = NULL, duf = NULL;
23 Vec pf = NULL, pJf = NULL;
24 Vec df = NULL;
25 Vec ys = NULL;
26
27 /*parameters:*/
28 int max_nonlinear_iterations;
29 int configuration_type;
30
31 /*Recover parameters: */
32 femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
33 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
34 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
35
36 count=1;
37 converged=false;
38
39 /*Start non-linear iteration using input velocity: */
40 GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
41 Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters);
42
43 //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
44 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
45 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
46
47 for(;;){
48
49 VecFree(&old_ug);old_ug=ug;
50 VecFree(&old_uf);old_uf=uf;
51
52 /*Solver forward model*/
53 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
54 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
55 Reduceloadx(pf,Kfs,ys);MatFree(&Kfs);
56 Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);VecFree(&df);
57 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
58 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);VecFree(&ug);
59
60 /*Check convergence*/
61 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
62 MatFree(&Kff);VecFree(&pf);
63 if(converged==true) break;
64 if(count>=max_nonlinear_iterations){
65 _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
66 break;
67 }
68
69 /*Prepare next iteration using Newton's method*/
70 SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
71 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
72 Reduceloadx(pf,Kfs,ys); MatFree(&Kfs);
73
74 VecDuplicate(pf,&pJf);
75 MatMultPatch(Kff,uf,pJf); MatFree(&Kff);
76 VecScale(pJf,-1.);
77 VecAXPY(pJf,+1.,pf); VecFree(&pf);
78
79 CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
80 Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);//VecFree(&pJf);
81 VecAXPY(uf,1.,duf); VecFree(&duf);
82 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
83 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
84
85 count++;
86 }
87
88 _printf_(VerboseConvergence(),"\n total number of iterations: %i\n",count-1);
89
90 /*clean-up*/
91 VecFree(&uf);
92 VecFree(&ug);
93 VecFree(&old_ug);
94 VecFree(&old_uf);
95}
Note: See TracBrowser for help on using the repository browser.