source: issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp@ 25549

Last change on this file since 25549 was 25549, checked in by Mathieu Morlighem, 5 years ago

CHG: minor esthetics
File size: 4.5 KB
Line 
1/*!\file: solutionsequence_nonlinear.cpp
2 * \brief: core of a non-linear solution, using fixed-point method
3 */
4
5#include "./solutionsequences.h"
6#include "../toolkits/toolkits.h"
7#include "../classes/classes.h"
8#include "../shared/shared.h"
9#include "../modules/modules.h"
10
11void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads){
12
13 /*intermediary: */
14 Matrix<IssmDouble>* Kff = NULL;
15 Matrix<IssmDouble>* Kfs = NULL;
16 Vector<IssmDouble>* ug = NULL;
17 Vector<IssmDouble>* uf = NULL;
18 Vector<IssmDouble>* old_uf = NULL;
19 Vector<IssmDouble>* pf = NULL;
20 Vector<IssmDouble>* df = NULL;
21 Vector<IssmDouble>* ys = NULL;
22
23 int constraints_converged;
24 int num_unstable_constraints;
25
26 /*parameters:*/
27 int min_mechanical_constraints;
28 int max_nonlinear_iterations;
29 int configuration_type;
30 IssmDouble eps_res,eps_rel,eps_abs;
31
32 /*Recover parameters: */
33 femmodel->parameters->FindParam(&min_mechanical_constraints,StressbalanceRiftPenaltyThresholdEnum);
34 femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
35 femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
36 femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
37 femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
38 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
39 femmodel->UpdateConstraintsx();
40
41 /*Were loads requested as output? : */
42 Loads* savedloads=NULL;
43 if(conserve_loads){
44 savedloads = static_cast<Loads*>(femmodel->loads->Copy());
45 }
46
47 int count=0;
48 bool converged=false;
49
50 /*Start non-linear iteration using input velocity: */
51 GetSolutionFromInputsx(&ug,femmodel);
52 Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
53
54 /*Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)*/
55 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
56 InputUpdateFromSolutionx(femmodel,ug);
57
58 /*allocate the matrices once and reuse them per iteration*/
59 if(femmodel->loads->numrifts == 0){
60 AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
61 }
62
63 for(;;){
64
65 //save pointer to old velocity
66 delete old_uf;old_uf=uf;
67 delete ug;
68
69 if(femmodel->loads->numrifts){
70 AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
71 }
72 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel, true);
73 CreateNodalConstraintsx(&ys,femmodel->nodes);
74 Reduceloadx(pf, Kfs, ys);
75 femmodel->profiler->Start(SOLVER);
76 Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
77 femmodel->profiler->Stop(SOLVER);
78
79 Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
80
81 convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs);
82 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
83 InputUpdateFromSolutionx(femmodel,ug);
84
85 /*Clean up if rifts*/
86 if(femmodel->loads->numrifts){
87 delete Kfs; delete Kff; delete pf; delete df;
88 }
89
90 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
91 if(VerboseConvergence()) _printf0_(" number of unstable constraints: " << num_unstable_constraints << "\n");
92
93 //rift convergence
94 if (!constraints_converged) {
95 if (converged){
96 if (num_unstable_constraints <= min_mechanical_constraints) converged=true;
97 else converged=false;
98 }
99 }
100
101 /*Increase count: */
102 count++;
103 if(converged==true){
104 femmodel->results->AddResult(new GenericExternalResult<int>(femmodel->results->Size()+1,StressbalanceConvergenceNumStepsEnum,count));
105 break;
106 }
107 if(count>=max_nonlinear_iterations){
108 _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
109 converged=true;
110 femmodel->results->AddResult(new GenericExternalResult<int>(femmodel->results->Size()+1,StressbalanceConvergenceNumStepsEnum,max_nonlinear_iterations));
111 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
112 InputUpdateFromSolutionx(femmodel,ug);
113 break;
114 }
115
116 /*Set the matrix entries to zero if we do an other iteration*/
117 if(femmodel->loads->numrifts==0){
118 Kff->SetZero();
119 Kfs->SetZero();
120 df->Set(0);
121 pf->Set(0);
122 }
123 }
124
125 /*delete matrices after the iteration loop*/
126 if(femmodel->loads->numrifts==0){
127 delete Kff; delete pf; delete df; delete Kfs;
128 }
129
130 if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count << "\n");
131
132 /*clean-up*/
133 if(conserve_loads){
134 delete femmodel->loads;
135 int index=femmodel->AnalysisIndex(configuration_type);
136 femmodel->loads_list[index]=savedloads;
137 femmodel->loads=savedloads;
138 }
139 delete uf;
140 delete ug;
141 delete old_uf;
142}
Note: See TracBrowser for help on using the repository browser.