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

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

Large commit, to strip out Petsc from the code. Created new Matrix and Vector
wrapper objects to the Petsc Mat and Vec objects. These wrappers make it possible
to hook up a different package to handle matrix and vector operations.

File size: 3.8 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 Matrix* Kff = NULL;
21 Matrix* Kfs = NULL;
22 Matrix* Jff = NULL;
23 Vector* ug = NULL;
24 Vector* old_ug = NULL;
25 Vector* uf = NULL;
26 Vector* old_uf = NULL;
27 Vector* duf = NULL;
28 Vector* pf = NULL;
29 Vector* pJf = NULL;
30 Vector* df = NULL;
31 Vector* 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) break;
70 if(count>=max_nonlinear_iterations){
71 _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
72 break;
73 }
74
75 /*Prepare next iteration using Newton's method*/
76 SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
77 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
78 Reduceloadx(pf,Kfs,ys); xdelete(&Kfs);
79
80 pJf=pf->Duplicate(); Kff->MatMult(uf,pJf); xdelete(&Kff);
81 pJf->Scale(-1.0); pJf->AXPY(pf,+1.0); xdelete(&pf);
82
83 CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
84 Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf);
85 uf->AXPY(duf, 1.0); xdelete(&duf);
86 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
87 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
88
89 count++;
90 }
91
92 _printf_(VerboseConvergence(),"\n total number of iterations: %i\n",count-1);
93
94 /*clean-up*/
95 xdelete(&uf);
96 xdelete(&ug);
97 xdelete(&old_ug);
98 xdelete(&old_uf);
99}
Note: See TracBrowser for help on using the repository browser.