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

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

CHG: cosmetics, removing all deboule blank lines and indent single white lines correctly
File size: 3.5 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 bool lowmem=0;
36 int configuration_type;
37
38 /*Recover parameters: */
39 kflag=1; pflag=1;
40 femmodel->parameters->FindParam(&lowmem,SettingsLowmemEnum);
41 femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
42 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
43 femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
44
45 count=1;
46 converged=false;
47
48 if(VerboseSolution()) _pprintLine_("starting direct shooting method");
49 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum);
50 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum);
51 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
52
53 for(;;){
54
55 SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
56 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
57 Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf);
58 Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
59 xdelete(&tf_old); tf_old=tf->Duplicate();
60 xdelete(&Kff);xdelete(&pf);xdelete(&tg); xdelete(&df);
61 Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys);
62 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,tg);
63
64 ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
65
66 if (!converged){
67 if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints);
68 if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
69 if (count>=thermal_maxiter){
70 converged=true;
71 _pprintLine_(" maximum number of iterations (" << thermal_maxiter << ") exceeded");
72 }
73 }
74 count++;
75
76 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
77
78 if(converged)break;
79 }
80
81 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,tg);
82 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,melting_offset,MeltingOffsetEnum);
83
84 /*Free ressources: */
85 xdelete(&tg);
86 xdelete(&tf);
87 xdelete(&tf_old);
88 xdelete(&ys);
89}
Note: See TracBrowser for help on using the repository browser.