source: issm/oecreview/Archive/14312-15392/ISSM-15351-15352.diff@ 15393

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

NEW: adding Archive/14312-15392 for oecreview
File size: 8.6 KB
RevLine 
[15393]1Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
2===================================================================
3--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15351)
4+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15352)
5@@ -8,17 +8,17 @@
6 #include "../modules/modules.h"
7
8 void solutionsequence_hydro_nonlinear(FemModel* femmodel){
9-
10 /*solution : */
11 Vector<IssmDouble>* ug_sed=NULL;
12 Vector<IssmDouble>* ug_epl=NULL;
13 Vector<IssmDouble>* uf=NULL;
14- Vector<IssmDouble>* old_uf=NULL;
15- Vector<IssmDouble>* aged_ug_sed=NULL;
16- Vector<IssmDouble>* aged_ug_epl=NULL;
17+ Vector<IssmDouble>* uf_int_iter=NULL;
18+ Vector<IssmDouble>* ug_sed_main_iter=NULL;
19+ Vector<IssmDouble>* ug_epl_main_iter=NULL;
20 Vector<IssmDouble>* ys=NULL;
21 Vector<IssmDouble>* dug=NULL;
22-
23+ Vector<IssmDouble>* old_ug=NULL;
24+
25 Matrix<IssmDouble>* Kff=NULL;
26 Matrix<IssmDouble>* Kfs=NULL;
27 Vector<IssmDouble>* pf=NULL;
28@@ -33,14 +33,14 @@
29 IssmDouble sediment_kmax,time;
30 IssmDouble eps_hyd;
31 IssmDouble ndu_sed,nu_sed;
32- IssmDouble ndu_epl,nu_epl;
33+ IssmDouble ndu_epl,nu_epl;
34
35 /*Recover parameters: */
36 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME
37 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
38 femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
39 femmodel->parameters->FindParam(&time,TimeEnum);
40- hydro_maxiter=100;
41+ hydro_maxiter=150;
42 hydrocount=1;
43
44 /*Iteration on the two layers*/
45@@ -56,11 +56,11 @@
46 sedcount=1;
47 eplcount=1;
48 //save pointer to old velocity
49- delete aged_ug_sed;
50- aged_ug_sed=ug_sed;
51+ delete ug_sed_main_iter;
52+ ug_sed_main_iter=ug_sed;
53 if(isefficientlayer){
54- delete aged_ug_epl;
55- aged_ug_epl=ug_epl;
56+ delete ug_epl_main_iter;
57+ ug_epl_main_iter=ug_epl;
58 }
59
60 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
61@@ -77,9 +77,9 @@
62 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
63 Reduceloadx(pf,Kfs,ys); delete Kfs;
64 delete uf;
65- Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
66- delete old_uf; old_uf=uf->Duplicate();
67- if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting aged_ug_sed*/
68+ Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
69+ delete uf_int_iter; uf_int_iter=uf->Duplicate();
70+ if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting ug_sed_main_iter*/
71 delete Kff; delete pf; delete df;
72
73 Mergesolutionfromftogx(&ug_sed,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
74@@ -119,15 +119,15 @@
75 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
76 Reduceloadx(pf,Kfs,ys); delete Kfs;
77 delete uf;
78- Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
79- delete old_uf; old_uf=uf->Duplicate();
80+ Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
81+ delete uf_int_iter; uf_int_iter=uf->Duplicate();
82 if(eplcount>1) delete ug_epl;
83 delete Kff;delete pf;
84 delete df;
85 Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
86 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
87- //femmodel->HydrologyEPLupdateDomainx();
88 ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
89+ femmodel->HydrologyEPLupdateDomainx();
90
91 if (!eplconverged){
92 if(VerboseConvergence()) _printf0_(" # EPL unstable constraints = " << num_unstable_constraints << "\n");
93@@ -137,12 +137,11 @@
94 }
95 }
96 eplcount++;
97-
98+
99 if(eplconverged){
100 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum);
101 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
102 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
103- //femmodel->HydrologyEPLupdateDomainx();
104 break;
105 }
106 }
107@@ -151,14 +150,14 @@
108 if(!hydroconverged){
109 //compute norm(du)/norm(u)
110 dug=ug_sed->Duplicate(); _assert_(dug);
111- aged_ug_sed->Copy(dug);
112+ ug_sed_main_iter->Copy(dug);
113 dug->AYPX(ug_sed,-1.0);
114- ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO);
115+ ndu_sed=dug->Norm(NORM_TWO); nu_sed=ug_sed_main_iter->Norm(NORM_TWO);
116 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
117 if (!xIsNan<IssmDouble>(eps_hyd)){
118 if (!isefficientlayer){
119 if ((ndu_sed/nu_sed)<eps_hyd){
120- if(VerboseConvergence()) _printf0_(setw(50) << left << " Converged \n");
121+ if(VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n");
122 hydroconverged=true;
123 }
124 else{
125@@ -168,15 +167,15 @@
126 }
127 else{
128 dug=ug_epl->Duplicate();_assert_(dug);
129- aged_ug_epl->Copy(dug);_assert_(aged_ug_epl);
130+ ug_epl_main_iter->Copy(dug);_assert_(ug_epl_main_iter);
131 dug->AYPX(ug_epl,-1.0);
132 ndu_epl=dug->Norm(NORM_TWO);
133- nu_epl=aged_ug_epl->Norm(NORM_TWO);
134+ nu_epl=ug_epl_main_iter->Norm(NORM_TWO);
135
136 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!");
137 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
138- if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<eps_hyd){
139- if (VerboseConvergence()) _printf0_(setw(50) << left << " Converged \n");
140+ if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<(eps_hyd*10)){
141+ if (VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n");
142 hydroconverged=true;
143 }
144 else{
145@@ -201,9 +200,8 @@
146 delete ug_epl;
147 delete ug_sed;
148 delete uf;
149- delete old_uf;
150- delete aged_ug_sed;
151- delete aged_ug_epl;
152+ delete uf_int_iter;
153+ delete ug_sed_main_iter;
154+ delete ug_epl_main_iter;
155 delete dug;
156-
157 }
158Index: ../trunk-jpl/src/c/classes/FemModel.cpp
159===================================================================
160--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 15351)
161+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 15352)
162@@ -852,6 +852,7 @@
163 element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
164 element->CreateKMatrix(Kff_temp,NULL);
165 }
166+
167 for (i=0;i<this->loads->Size();i++){
168 load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
169 if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
170@@ -871,6 +872,7 @@
171 element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
172 element->CreateKMatrix(Kff,Kfs);
173 }
174+
175 for (i=0;i<this->loads->Size();i++){
176 load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
177 if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
178Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
179===================================================================
180--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15351)
181+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15352)
182@@ -6131,7 +6131,6 @@
183 for(int ig=gauss->begin();ig<gauss->end();ig++){
184
185 gauss->GaussPoint(ig);
186-
187 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
188
189 /*Diffusivity*/
190@@ -6156,7 +6155,6 @@
191 &Ke->values[0],1);
192 }
193 }
194-
195 /*Clean up and return*/
196 delete gauss;
197 return Ke;
198@@ -6243,11 +6241,9 @@
199
200 /* Start looping on the number of gaussian points: */
201 gauss=new GaussTria(2);
202-
203 for(int ig=gauss->begin();ig<gauss->end();ig++){
204
205 gauss->GaussPoint(ig);
206-
207 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
208 GetNodalFunctions(basis, gauss);
209
210@@ -6310,7 +6306,6 @@
211 for(int ig=gauss->begin();ig<gauss->end();ig++){
212
213 gauss->GaussPoint(ig);
214-
215 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
216 GetNodalFunctions(basis, gauss);
217
218@@ -6690,7 +6685,7 @@
219
220 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
221 this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
222- if(eplhead[i]>=h_max && this->AnyActive()){
223+ if(eplhead[i]>=h_max && this->AnyActive()){
224 for(j=0;j<numdof;j++){
225 if(old_active[j]>0.){
226 vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
Note: See TracBrowser for help on using the repository browser.