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

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

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