source: issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp@ 14263

Last change on this file since 14263 was 14263, checked in by Mathieu Morlighem, 12 years ago

CHG: lowmem not needed
File size: 3.2 KB
Line 
1/*!\file: solver_thermal_nonlinear.cpp
2 * \brief: core of the thermal solution
3 */
4
5#include "../toolkits/toolkits.h"
6#include "../classes/objects/objects.h"
7#include "../io/io.h"
8#include "../EnumDefinitions/EnumDefinitions.h"
9#include "../modules/modules.h"
10
11void solver_thermal_nonlinear(FemModel* femmodel){
12
13 /*solution : */
14 Vector<IssmDouble>* tg=NULL;
15 Vector<IssmDouble>* tf=NULL;
16 Vector<IssmDouble>* tf_old=NULL;
17 Vector<IssmDouble>* ys=NULL;
18 IssmDouble melting_offset;
19
20 /*intermediary: */
21 Matrix<IssmDouble>* Kff=NULL;
22 Matrix<IssmDouble>* Kfs=NULL;
23 Vector<IssmDouble>* pf=NULL;
24 Vector<IssmDouble>* df=NULL;
25
26 bool converged;
27 int constraints_converged;
28 int num_unstable_constraints;
29 int count;
30 int thermal_penalty_threshold;
31 int thermal_maxiter;
32
33 /*parameters:*/
34 int kflag,pflag;
35 int configuration_type;
36
37 /*Recover parameters: */
38 kflag=1; pflag=1;
39 femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
40 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
41 femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
42
43 count=1;
44 converged=false;
45
46 if(VerboseSolution()) _pprintLine_("starting direct shooting method");
47 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum);
48 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum);
49 femmodel->UpdateConstraintsx();
50
51 for(;;){
52
53 femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset);
54 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
55 Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf);
56 Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
57 xdelete(&tf_old); tf_old=tf->Duplicate();
58 xdelete(&Kff);xdelete(&pf);xdelete(&tg); xdelete(&df);
59 Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys);
60 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,tg);
61
62 ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
63
64 if (!converged){
65 if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints);
66 if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
67 if (count>=thermal_maxiter){
68 converged=true;
69 _pprintLine_(" maximum number of iterations (" << thermal_maxiter << ") exceeded");
70 }
71 }
72 count++;
73
74 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
75
76 if(converged)break;
77 }
78
79 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,tg);
80 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,melting_offset,MeltingOffsetEnum);
81
82 /*Free ressources: */
83 xdelete(&tg);
84 xdelete(&tf);
85 xdelete(&tf_old);
86 xdelete(&ys);
87}
Note: See TracBrowser for help on using the repository browser.