source: issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp@ 20669

Last change on this file since 20669 was 20669, checked in by Mathieu Morlighem, 9 years ago

NEW: added check for isinf in solution vector
File size: 13.6 KB
Line 
1#include "./HydrologyShreveAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
6
7/*Model processing*/
8void HydrologyShreveAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9
10 /*retrieve some parameters: */
11 int hydrology_model;
12 iomodel->FindConstant(&hydrology_model,HydrologyModelEnum);
13
14 if(hydrology_model!=HydrologyshreveEnum) return;
15
16 IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum);
17
18}/*}}}*/
19void HydrologyShreveAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
20 /*No loads*/
21}/*}}}*/
22void HydrologyShreveAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
23
24 /*Fetch parameters: */
25 int hydrology_model;
26 iomodel->FindConstant(&hydrology_model,HydrologyModelEnum);
27
28 /*Now, do we really want Shreve?*/
29 if(hydrology_model!=HydrologyshreveEnum) return;
30
31 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
32 ::CreateNodes(nodes,iomodel,HydrologyShreveAnalysisEnum,P1Enum);
33 iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
34}/*}}}*/
35int HydrologyShreveAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
36 return 1;
37}/*}}}*/
38void HydrologyShreveAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
39
40 /*Fetch data needed: */
41 int hydrology_model;
42 iomodel->FindConstant(&hydrology_model,HydrologyModelEnum);
43
44 /*Now, do we really want Shreve?*/
45 if(hydrology_model!=HydrologyshreveEnum) return;
46
47 /*Update elements: */
48 int counter=0;
49 for(int i=0;i<iomodel->numberofelements;i++){
50 if(iomodel->my_elements[i]){
51 Element* element=(Element*)elements->GetObjectByOffset(counter);
52 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
53 counter++;
54 }
55 }
56
57 iomodel->FetchDataToInput(elements,ThicknessEnum);
58 iomodel->FetchDataToInput(elements,SurfaceEnum);
59 iomodel->FetchDataToInput(elements,BaseEnum);
60 iomodel->FetchDataToInput(elements,SealevelEnum,0);
61 if(iomodel->domaintype!=Domain2DhorizontalEnum){
62 iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
63 iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
64 }
65 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
66 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
67 iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
68 iomodel->FetchDataToInput(elements,WatercolumnEnum);
69
70 elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum);
71}/*}}}*/
72void HydrologyShreveAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
73
74 /*retrieve some parameters: */
75 int hydrology_model;
76 iomodel->FindConstant(&hydrology_model,HydrologyModelEnum);
77
78 /*Now, do we really want Shreve?*/
79 if(hydrology_model!=HydrologyshreveEnum) return;
80
81 parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
82 parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum));
83
84}/*}}}*/
85
86/*Finite Element Analysis*/
87void HydrologyShreveAnalysis::Core(FemModel* femmodel){/*{{{*/
88 _error_("not implemented");
89}/*}}}*/
90ElementVector* HydrologyShreveAnalysis::CreateDVector(Element* element){/*{{{*/
91 /*Default, return NULL*/
92 return NULL;
93}/*}}}*/
94void HydrologyShreveAnalysis::CreateHydrologyWaterVelocityInput(Element* element){/*{{{*/
95
96 /*Intermediaries*/
97 IssmDouble dsdx,dsdy,dbdx,dbdy,w;
98
99 /*Retrieve all inputs and parameters*/
100 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
101 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
102 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum);
103 IssmDouble mu_water = element->GetMaterialParameter(MaterialsMuWaterEnum);
104 Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
105 Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
106 Input* bedslopex_input = element->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
107 Input* bedslopey_input = element->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
108 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
109
110 /*Fetch number of vertices and allocate output*/
111 int numvertices = element->GetNumberOfVertices();
112 IssmDouble* vx = xNew<IssmDouble>(numvertices);
113 IssmDouble* vy = xNew<IssmDouble>(numvertices);
114
115 Gauss* gauss=element->NewGauss();
116 for(int iv=0;iv<numvertices;iv++){
117 gauss->GaussVertex(iv);
118 surfaceslopex_input->GetInputValue(&dsdx,gauss);
119 surfaceslopey_input->GetInputValue(&dsdy,gauss);
120 bedslopex_input->GetInputValue(&dbdx,gauss);
121 bedslopey_input->GetInputValue(&dbdy,gauss);
122 watercolumn_input->GetInputValue(&w,gauss);
123
124 /* Water velocity x and y components */
125 vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
126 vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
127 }
128
129 /*clean-up*/
130 delete gauss;
131
132 /*Add to inputs*/
133 element->AddInput(HydrologyWaterVxEnum,vx,P1Enum);
134 element->AddInput(HydrologyWaterVyEnum,vy,P1Enum);
135 xDelete<IssmDouble>(vx);
136 xDelete<IssmDouble>(vy);
137}/*}}}*/
138ElementMatrix* HydrologyShreveAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
139_error_("Not implemented");
140}/*}}}*/
141ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/
142
143 /*Intermediaries */
144 IssmDouble diffusivity;
145 IssmDouble Jdet,D_scalar,dt,h;
146 IssmDouble vx,vy,vel,dvxdx,dvydy;
147 IssmDouble dvx[2],dvy[2];
148 IssmDouble* xyz_list = NULL;
149
150 /*Fetch number of nodes and dof for this finite element*/
151 int numnodes = element->GetNumberOfNodes();
152
153 /*Initialize Element vector and other vectors*/
154 ElementMatrix* Ke = element->NewElementMatrix();
155 IssmDouble* basis = xNew<IssmDouble>(numnodes);
156 IssmDouble* B = xNew<IssmDouble>(2*numnodes);
157 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);
158 IssmDouble D[2][2]={0.};
159
160 /*Create water velocity vx and vy from current inputs*/
161 CreateHydrologyWaterVelocityInput(element);
162
163 /*Retrieve all inputs and parameters*/
164 element->GetVerticesCoordinates(&xyz_list);
165 element->FindParam(&dt,TimesteppingTimeStepEnum);
166 element->FindParam(&diffusivity,HydrologyshreveStabilizationEnum);
167 Input* vx_input=element->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
168 Input* vy_input=element->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
169 h = element->CharacteristicLength();
170
171 /* Start looping on the number of gaussian points: */
172 Gauss* gauss=element->NewGauss(2);
173 for(int ig=gauss->begin();ig<gauss->end();ig++){
174 gauss->GaussPoint(ig);
175
176 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
177 element->NodalFunctions(basis,gauss);
178
179 vx_input->GetInputValue(&vx,gauss);
180 vy_input->GetInputValue(&vy,gauss);
181 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
182 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
183
184 D_scalar=gauss->weight*Jdet;
185
186 TripleMultiply(basis,1,numnodes,1,
187 &D_scalar,1,1,0,
188 basis,1,numnodes,0,
189 Ke->values,1);
190
191 GetB(B,element,xyz_list,gauss);
192 GetBprime(Bprime,element,xyz_list,gauss);
193
194 dvxdx=dvx[0];
195 dvydy=dvy[1];
196 D_scalar=dt*gauss->weight*Jdet;
197
198 D[0][0]=D_scalar*dvxdx;
199 D[1][1]=D_scalar*dvydy;
200 TripleMultiply(B,2,numnodes,1,
201 &D[0][0],2,2,0,
202 B,2,numnodes,0,
203 &Ke->values[0],1);
204
205 D[0][0]=D_scalar*vx;
206 D[1][1]=D_scalar*vy;
207 TripleMultiply(B,2,numnodes,1,
208 &D[0][0],2,2,0,
209 Bprime,2,numnodes,0,
210 &Ke->values[0],1);
211
212 /*Artificial diffusivity*/
213 vel=sqrt(vx*vx+vy*vy);
214 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
215 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
216 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
217 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
218 TripleMultiply(Bprime,2,numnodes,1,
219 &D[0][0],2,2,0,
220 Bprime,2,numnodes,0,
221 &Ke->values[0],1);
222 }
223
224 /*Clean up and return*/
225 xDelete<IssmDouble>(xyz_list);
226 xDelete<IssmDouble>(basis);
227 xDelete<IssmDouble>(B);
228 xDelete<IssmDouble>(Bprime);
229 delete gauss;
230 return Ke;
231}/*}}}*/
232ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
233
234 /*Skip if water or ice shelf element*/
235 if(element->IsFloating()) return NULL;
236
237 /*Intermediaries */
238 IssmDouble Jdet,dt;
239 IssmDouble mb,oldw;
240 IssmDouble* xyz_list = NULL;
241
242 /*Fetch number of nodes and dof for this finite element*/
243 int numnodes = element->GetNumberOfNodes();
244
245 /*Initialize Element vector and other vectors*/
246 ElementVector* pe = element->NewElementVector();
247 IssmDouble* basis = xNew<IssmDouble>(numnodes);
248
249 /*Retrieve all inputs and parameters*/
250 element->GetVerticesCoordinates(&xyz_list);
251 element->FindParam(&dt,TimesteppingTimeStepEnum);
252 Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
253 Input* oldw_input = element->GetInput(WaterColumnOldEnum); _assert_(oldw_input);
254
255 /*Initialize mb_correction to 0, do not forget!:*/
256 /* Start looping on the number of gaussian points: */
257 Gauss* gauss=element->NewGauss(2);
258 for(int ig=gauss->begin();ig<gauss->end();ig++){
259 gauss->GaussPoint(ig);
260
261 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
262 element->NodalFunctions(basis,gauss);
263
264 mb_input->GetInputValue(&mb,gauss);
265 oldw_input->GetInputValue(&oldw,gauss);
266
267 if(dt!=0.){
268 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
269 }
270 else{
271 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
272 }
273 }
274
275 /*Clean up and return*/
276 xDelete<IssmDouble>(xyz_list);
277 xDelete<IssmDouble>(basis);
278 delete gauss;
279 return pe;
280}/*}}}*/
281void HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
282 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
283 * For node i, Bi can be expressed in the actual coordinate system
284 * by:
285 * Bi=[ N ]
286 * [ N ]
287 * where N is the finiteelement function for node i.
288 *
289 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
290 */
291
292 /*Fetch number of nodes for this finite element*/
293 int numnodes = element->GetNumberOfNodes();
294
295 /*Get nodal functions*/
296 IssmDouble* basis=xNew<IssmDouble>(numnodes);
297 element->NodalFunctions(basis,gauss);
298
299 /*Build B: */
300 for(int i=0;i<numnodes;i++){
301 B[numnodes*0+i] = basis[i];
302 B[numnodes*1+i] = basis[i];
303 }
304
305 /*Clean-up*/
306 xDelete<IssmDouble>(basis);
307}/*}}}*/
308void HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
309 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
310 * For node i, Bi' can be expressed in the actual coordinate system
311 * by:
312 * Bi_prime=[ dN/dx ]
313 * [ dN/dy ]
314 * where N is the finiteelement function for node i.
315 *
316 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
317 */
318
319 /*Fetch number of nodes for this finite element*/
320 int numnodes = element->GetNumberOfNodes();
321
322 /*Get nodal functions derivatives*/
323 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
324 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
325
326 /*Build B': */
327 for(int i=0;i<numnodes;i++){
328 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
329 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
330 }
331
332 /*Clean-up*/
333 xDelete<IssmDouble>(dbasis);
334
335}/*}}}*/
336void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
337 element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
338}/*}}}*/
339void HydrologyShreveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
340 _error_("Not implemented yet");
341}/*}}}*/
342void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
343
344 /*Intermediary*/
345 int* doflist = NULL;
346
347 /*Fetch number of nodes for this finite element*/
348 int numnodes = element->GetNumberOfNodes();
349
350 /*Fetch dof list and allocate solution vector*/
351 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
352 IssmDouble* values = xNew<IssmDouble>(numnodes);
353
354 /*Use the dof list to index into the solution vector: */
355 for(int i=0;i<numnodes;i++){
356 values[i]=solution[doflist[i]];
357 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
358 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
359 if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
360 }
361
362 /*Add input to the element: */
363 element->AddInput(WatercolumnEnum,values,element->GetElementType());
364
365 /*Free ressources:*/
366 xDelete<IssmDouble>(values);
367 xDelete<int>(doflist);
368}/*}}}*/
369void HydrologyShreveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
370 /*Default, do nothing*/
371 return;
372}/*}}}*/
373
374
375/*Needed changes to switch to the Johnson formulation*//*{{{*/
376/*All the changes are to be done in the velocity computation.
377 The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
378 'p' and 'q' which are the exponent of the Manning formula for laminar (p=2,q=1) or turbulent (p=2/3,q=1/2) flow
379 'R' the hydraulic radius
380 'n' the manning roughness coeficient
381
382 With these, the velocity reads ;
383
384 v= - (1/n)* pow(R,p)*pow((grad phi(rho_water*g)),q)
385
386 you should also redefine the water pressure potential 'phi' with respect to the effective pressure deffinition given in Johson:
387 phi=(rho_ice*g*( surface + ((rho_water/rho_ice)-1)*base - k_n*((thickness* grad(base))/omega) )
388
389 where
390 'omega' is the fractional area of the base occupied by the water film
391 'k_n' is a parameter
392 This last equation derives from the effective pressure definition developped in Alley 1989
393*/
394/*}}}*/
Note: See TracBrowser for help on using the repository browser.