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

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

Added Newton's method for Pattyn, parallel only for now, to be improved
File size: 3.4 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 Mat Kff = NULL, Kfs = NULL, Jff = NULL;
17 Vec ug = NULL, old_ug = NULL;
18 Vec uf = NULL, old_uf = NULL, duf = NULL;
19 Vec pf = NULL, pJf = NULL;
20 Vec df = NULL;
21 Vec ys = NULL;
22
23 bool converged;
24 int num_unstable_constraints;
25 int count;
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 //save pointer to old velocity
50 VecFree(&old_ug);old_ug=ug;
51 VecFree(&old_uf);old_uf=uf;
52
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);
57
58 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
59
60 CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
61
62 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);
63 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
64 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
65 VecDuplicate(pf,&pJf);
66 MatMultPatch(Kff,uf,pJf);
67 VecAXPY(pJf,-1.,pf);
68 Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
69 MatFree(&Kff);VecFree(&pf); VecFree(&df);
70 VecAXPY(uf,-1.,duf);
71 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
72
73 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
74
75 /*Increase count: */
76 count++;
77 if(converged==true)break;
78 if(count>=max_nonlinear_iterations){
79 _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
80 break;
81 }
82 }
83
84 _printf_(VerboseConvergence(),"\n total number of iterations: %i\n",count-1);
85
86 /*clean-up*/
87 VecFree(&uf);
88 VecFree(&ug);
89 VecFree(&old_ug);
90 VecFree(&old_uf);
91}
Note: See TracBrowser for help on using the repository browser.