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

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

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