source: issm/oecreview/Archive/16554-17801/ISSM-16760-16761.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 9.9 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16760)
4+++ ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16761)
5@@ -85,5 +85,28 @@
6 element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
7 }/*}}}*/
8 void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
9- _error_("not implemented yet");
10+
11+ /*Intermediary*/
12+ int* doflist = NULL;
13+
14+ /*Fetch number of nodes for this finite element*/
15+ int numnodes = element->GetNumberOfNodes();
16+
17+ /*Fetch dof list and allocate solution vector*/
18+ element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
19+ IssmDouble* values = xNew<IssmDouble>(numnodes);
20+
21+ /*Use the dof list to index into the solution vector: */
22+ for(int i=0;i<numnodes;i++){
23+ values[i]=solution[doflist[i]];
24+ if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
25+ if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
26+ }
27+
28+ /*Add input to the element: */
29+ element->AddInput(WatercolumnEnum,values,P1Enum);
30+
31+ /*Free ressources:*/
32+ xDelete<IssmDouble>(values);
33+ xDelete<int>(doflist);
34 }/*}}}*/
35Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
36===================================================================
37--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16760)
38+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16761)
39@@ -99,5 +99,15 @@
40 element->GetSolutionFromInputsOneDof(solution,EplHeadEnum);
41 }/*}}}*/
42 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
43- _error_("not implemented yet");
44+ int meshtype;
45+ element->FindParam(&meshtype,MeshTypeEnum);
46+ switch(meshtype){
47+ case Mesh2DhorizontalEnum:
48+ element->InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
49+ break;
50+ case Mesh3DEnum:
51+ element->InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum);
52+ break;
53+ default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
54+ }
55 }/*}}}*/
56Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
57===================================================================
58--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16760)
59+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16761)
60@@ -140,5 +140,60 @@
61 element->GetSolutionFromInputsOneDof(solution,SedimentHeadEnum);
62 }/*}}}*/
63 void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
64- _error_("not implemented yet");
65+
66+ int meshtype;
67+ bool converged;
68+ int* doflist=NULL;
69+ Element* basalelement=NULL;
70+
71+ element->FindParam(&meshtype,MeshTypeEnum);
72+ if(meshtype!=Mesh2DhorizontalEnum){
73+ if(!element->IsOnBed()) return;
74+ basalelement=element->SpawnBasalElement();
75+ }
76+ else{
77+ basalelement = element;
78+ }
79+
80+ /*Fetch number of nodes for this finite element*/
81+ int numnodes = basalelement->GetNumberOfNodes();
82+
83+ /*Fetch dof list and allocate solution vector*/
84+ basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
85+ IssmDouble* values = xNew<IssmDouble>(numnodes);
86+ IssmDouble* residual = xNew<IssmDouble>(numnodes);
87+
88+ /*Use the dof list to index into the solution vector: */
89+ for(int i=0;i<numnodes;i++){
90+ values[i] =solution[doflist[i]];
91+ if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
92+ }
93+
94+ /*If converged keep the residual in mind*/
95+ element->GetInputValue(&converged,ConvergedEnum);
96+
97+ /*Get inputs*/
98+ if(converged){
99+ IssmDouble penalty_factor,kmax,kappa,h_max;
100+ element->FindParam(&kmax,HydrologySedimentKmaxEnum);
101+ element->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
102+
103+ kappa=kmax*pow(10.,penalty_factor);
104+
105+ for(int i=0;i<numnodes;i++){
106+ basalelement->GetHydrologyDCInefficientHmax(&h_max,i);
107+ if(values[i]>h_max) residual[i] = kappa*(values[i]-h_max);
108+ else residual[i] = 0.;
109+ }
110+ }
111+
112+ /*Add input to the element: */
113+ element->AddBasalInput(SedimentHeadEnum,values,P1Enum);
114+ element->AddBasalInput(SedimentHeadResidualEnum,residual,P1Enum);
115+
116+ /*Free ressources:*/
117+ xDelete<IssmDouble>(values);
118+ xDelete<IssmDouble>(residual);
119+ xDelete<int>(doflist);
120+ if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
121 }/*}}}*/
122Index: ../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp
123===================================================================
124--- ../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 16760)
125+++ ../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 16761)
126@@ -73,6 +73,21 @@
127 _error_("not implemented yet");
128 }/*}}}*/
129 void L2ProjectionEPLAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
130- _error_("not implemented yet");
131+ int inputenum,meshtype;
132+
133+ element->FindParam(&inputenum,InputToL2ProjectEnum);
134+ element->FindParam(&meshtype,MeshTypeEnum);
135+ switch(meshtype){
136+ case Mesh2DhorizontalEnum:
137+ element->InputUpdateFromSolutionOneDof(solution,inputenum);
138+ break;
139+ case Mesh2DverticalEnum:
140+ element->InputUpdateFromSolutionOneDof(solution,inputenum);
141+ break;
142+ case Mesh3DEnum:
143+ element->InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
144+ break;
145+ default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
146+ }
147 }/*}}}*/
148
149Index: ../trunk-jpl/src/c/classes/Elements/Element.h
150===================================================================
151--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 16760)
152+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 16761)
153@@ -168,6 +168,7 @@
154
155 #ifdef _HAVE_HYDROLOGY_
156 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
157+ virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, int index)=0;
158 virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
159 virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
160 virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
161Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
162===================================================================
163--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16760)
164+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16761)
165@@ -7226,6 +7226,41 @@
166 *ph_max=h_max;
167 }
168 /*}}}*/
169+/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/
170+void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){
171+
172+ int hmax_flag;
173+ IssmDouble h_max;
174+ IssmDouble rho_ice,rho_water;
175+ IssmDouble thickness,bed;
176+ /*Get the flag to the limitation method*/
177+ this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
178+
179+ /*Switch between the different cases*/
180+ switch(hmax_flag){
181+ case 0:
182+ h_max=1.0e+10;
183+ break;
184+ case 1:
185+ parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum);
186+ break;
187+ case 2:
188+ rho_ice=matpar->GetRhoIce();
189+ rho_water=matpar->GetRhoFreshwater();
190+ this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum);
191+ this->GetInputValue(&bed,this->nodes[index],BedEnum);
192+ h_max=((rho_ice*thickness)/rho_water)+bed;
193+ break;
194+ case 3:
195+ _error_("Using normal stress not supported yet");
196+ break;
197+ default:
198+ _error_("no case higher than 3 for SedimentlimitFlag");
199+ }
200+ /*Assign output pointer*/
201+ *ph_max=h_max;
202+}
203+/*}}}*/
204 /*FUNCTION Tria::GetHydrologyTransfer{{{*/
205 void Tria::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
206
207Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
208===================================================================
209--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16760)
210+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16761)
211@@ -327,6 +327,7 @@
212 void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
213 void InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
214 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
215+ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index);
216 void GetHydrologyTransfer(Vector<IssmDouble>* transfer);
217 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
218 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
219Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
220===================================================================
221--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16760)
222+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16761)
223@@ -343,6 +343,7 @@
224 ElementVector* CreatePVectorHydrologyDCEfficient(void);
225 ElementVector* CreatePVectorL2ProjectionEPL(void);
226 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
227+ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");};
228 void GetHydrologyTransfer(Vector<IssmDouble>* transfer);
229 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
230 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
231Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
232===================================================================
233--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16760)
234+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16761)
235@@ -127,6 +127,7 @@
236 #endif
237 #ifdef _HAVE_HYDROLOGY_
238 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");};
239+ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");};
240 void GetHydrologyTransfer(Vector<IssmDouble>* transfer){_error_("not implemented yet");};
241 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
242 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
Note: See TracBrowser for help on using the repository browser.