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

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

NEW: large change to the code, to adapt to ADOLC requirements.

This change relates to the introduction of template classes and functions for the
Option.h abstract class. This is needed, because we want to make the Matlab
API independent from the libCore objects, which are dependent on the IssmDouble*
ADOLC type (adouble), when the Matlab API is dependent on the IssmPDouble* type (double).

To make them independent, we need to be able to specify at run time Options, Matrix and
Vector objects that hold either IssmDouble or IssmPDouble objects. The only way to do
that is through the use of templated classes for Option.h, Matrix and Vector.

The change gets rid of a lot of useless code (especially in the classes/objects/Options
directory), by introducing template versions of the same code.

The bulk of the changes to src/modules and src/mex modules is to adapt to this
new runtime declaration of templated Matrix, Vector and Option objects.
File size: 4.5 KB
RevLine 
[11322]1/*!\file: solver_nonlinear.cpp
2 * \brief: core of a non-linear solution, using fixed-point method
3 */
4
5#include "../toolkits/toolkits.h"
[12832]6#include "../classes/objects/objects.h"
[11322]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: */
[11327]16 bool converged;
17 int num_unstable_constraints;
18 int count;
[12477]19 IssmDouble kmax;
[13216]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;
[11322]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
[11679]55 xdelete(&old_ug);old_ug=ug;
56 xdelete(&old_uf);old_uf=uf;
[11322]57
[11324]58 /*Solver forward model*/
[11322]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);
[11679]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);
[11322]65
[11324]66 /*Check convergence*/
67 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
[11679]68 xdelete(&Kff); xdelete(&pf);
[12271]69 if(converged==true){
70 bool max_iteration_state=false;
71 int tempStep=1;
[12477]72 IssmDouble tempTime=1.0;
[12271]73 femmodel->results->AddObject(new BoolExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
74 break;
75 }
[11322]76 if(count>=max_nonlinear_iterations){
[12519]77 _pprintLine_(" maximum number of Newton iterations (" << max_nonlinear_iterations << ") exceeded");
[12271]78 bool max_iteration_state=true;
79 int tempStep=1;
[12477]80 IssmDouble tempTime=1.0;
[12271]81 femmodel->results->AddObject(new BoolExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
[11322]82 break;
83 }
[11324]84
85 /*Prepare next iteration using Newton's method*/
[11337]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);
[11679]88 Reduceloadx(pf,Kfs,ys); xdelete(&Kfs);
[11337]89
[11679]90 pJf=pf->Duplicate(); Kff->MatMult(uf,pJf); xdelete(&Kff);
91 pJf->Scale(-1.0); pJf->AXPY(pf,+1.0); xdelete(&pf);
[11337]92
93 CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
[11679]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);
[11324]97 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
[11332]98
99 count++;
[11322]100 }
101
[12515]102 if(VerboseConvergence()) _pprintLine_("\n total number of iterations: " << count-1);
[11322]103
104 /*clean-up*/
[11679]105 xdelete(&uf);
106 xdelete(&ug);
107 xdelete(&old_ug);
108 xdelete(&old_uf);
[11322]109}
Note: See TracBrowser for help on using the repository browser.