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

Last change on this file since 16939 was 16939, checked in by Eric.Larour, 11 years ago

CHG: fixed assert.
File size: 12.1 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 meshtype,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,BedEnum);
46 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
47 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
48 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
49 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
50 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
51 iomodel->FetchDataToInput(elements,WatercolumnEnum);
52
53 elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum);
54}/*}}}*/
55void HydrologyShreveAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
56
57 /*Fetch parameters: */
58 int hydrology_model;
59 iomodel->Constant(&hydrology_model,HydrologyModelEnum);
60
61 /*Now, do we really want Shreve?*/
62 if(hydrology_model!=HydrologyshreveEnum) return;
63
64 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
65 ::CreateNodes(nodes,iomodel,HydrologyShreveAnalysisEnum,P1Enum);
66 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
67}/*}}}*/
68void HydrologyShreveAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
69
70 /*retrieve some parameters: */
71 int hydrology_model;
72 iomodel->Constant(&hydrology_model,HydrologyModelEnum);
73
74 if(hydrology_model!=HydrologyshreveEnum) return;
75
76 IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum);
77
78}/*}}}*/
79void HydrologyShreveAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
80 /*No loads*/
81}/*}}}*/
82
83/*Finite Element Analysis*/
84ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/
85
86 /*Intermediaries */
87 IssmDouble diffusivity;
88 IssmDouble Jdet,D_scalar,dt,h;
89 IssmDouble vx,vy,vel,dvxdx,dvydy;
90 IssmDouble dvx[2],dvy[2];
91 IssmDouble* xyz_list = NULL;
92
93 /*Fetch number of nodes and dof for this finite element*/
94 int numnodes = element->GetNumberOfNodes();
95
96 /*Initialize Element vector and other vectors*/
97 ElementMatrix* Ke = element->NewElementMatrix();
98 IssmDouble* basis = xNew<IssmDouble>(numnodes);
99 IssmDouble* B = xNew<IssmDouble>(2*numnodes);
100 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);
101 IssmDouble D[2][2]={0.};
102
103 /*Create water velocity vx and vy from current inputs*/
104 CreateHydrologyWaterVelocityInput(element);
105
106 /*Retrieve all inputs and parameters*/
107 element->GetVerticesCoordinates(&xyz_list);
108 element->FindParam(&dt,TimesteppingTimeStepEnum);
109 Input* vx_input=element->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
110 Input* vy_input=element->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
111 h = element->CharacteristicLength();
112
113 /* Start looping on the number of gaussian points: */
114 Gauss* gauss=element->NewGauss(2);
115 for(int ig=gauss->begin();ig<gauss->end();ig++){
116 gauss->GaussPoint(ig);
117
118 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
119 element->NodalFunctions(basis,gauss);
120
121 vx_input->GetInputValue(&vx,gauss);
122 vy_input->GetInputValue(&vy,gauss);
123 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
124 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
125
126 D_scalar=gauss->weight*Jdet;
127
128 TripleMultiply(basis,1,numnodes,1,
129 &D_scalar,1,1,0,
130 basis,1,numnodes,0,
131 Ke->values,1);
132
133 GetB(B,element,xyz_list,gauss);
134 GetBprime(Bprime,element,xyz_list,gauss);
135
136 dvxdx=dvx[0];
137 dvydy=dvy[1];
138 D_scalar=dt*gauss->weight*Jdet;
139
140 D[0][0]=D_scalar*dvxdx;
141 D[1][1]=D_scalar*dvydy;
142 TripleMultiply(B,2,numnodes,1,
143 &D[0][0],2,2,0,
144 B,2,numnodes,0,
145 &Ke->values[0],1);
146
147 D[0][0]=D_scalar*vx;
148 D[1][1]=D_scalar*vy;
149 TripleMultiply(B,2,numnodes,1,
150 &D[0][0],2,2,0,
151 Bprime,2,numnodes,0,
152 &Ke->values[0],1);
153
154 /*Artificial diffusivity*/
155 vel=sqrt(vx*vx+vy*vy);
156 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
157 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
158 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
159 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
160 TripleMultiply(Bprime,2,numnodes,1,
161 &D[0][0],2,2,0,
162 Bprime,2,numnodes,0,
163 &Ke->values[0],1);
164 }
165
166 /*Clean up and return*/
167 xDelete<IssmDouble>(xyz_list);
168 xDelete<IssmDouble>(basis);
169 xDelete<IssmDouble>(B);
170 xDelete<IssmDouble>(Bprime);
171 delete gauss;
172 return Ke;
173}/*}}}*/
174ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
175
176 /*Skip if water or ice shelf element*/
177 if(element->IsFloating()) return NULL;
178
179 /*Intermediaries */
180 IssmDouble Jdet,dt;
181 IssmDouble mb,oldw;
182 IssmDouble* xyz_list = NULL;
183
184 /*Fetch number of nodes and dof for this finite element*/
185 int numnodes = element->GetNumberOfNodes();
186
187 /*Initialize Element vector and other vectors*/
188 ElementVector* pe = element->NewElementVector();
189 IssmDouble* basis = xNew<IssmDouble>(numnodes);
190
191 /*Retrieve all inputs and parameters*/
192 element->GetVerticesCoordinates(&xyz_list);
193 element->FindParam(&dt,TimesteppingTimeStepEnum);
194 Input* mb_input = element->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input);
195 Input* oldw_input = element->GetInput(WaterColumnOldEnum); _assert_(oldw_input);
196
197 /*Initialize mb_correction to 0, do not forget!:*/
198 /* Start looping on the number of gaussian points: */
199 Gauss* gauss=element->NewGauss(2);
200 for(int ig=gauss->begin();ig<gauss->end();ig++){
201 gauss->GaussPoint(ig);
202
203 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
204 element->NodalFunctions(basis,gauss);
205
206 mb_input->GetInputValue(&mb,gauss);
207 oldw_input->GetInputValue(&oldw,gauss);
208
209 if(dt!=0.){
210 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
211 }
212 else{
213 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
214 }
215 }
216
217 /*Clean up and return*/
218 xDelete<IssmDouble>(xyz_list);
219 xDelete<IssmDouble>(basis);
220 delete gauss;
221 return pe;
222}/*}}}*/
223void HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
224 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
225 * For node i, Bi can be expressed in the actual coordinate system
226 * by:
227 * Bi=[ N ]
228 * [ N ]
229 * where N is the finiteelement function for node i.
230 *
231 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
232 */
233
234 /*Fetch number of nodes for this finite element*/
235 int numnodes = element->GetNumberOfNodes();
236
237 /*Get nodal functions*/
238 IssmDouble* basis=xNew<IssmDouble>(numnodes);
239 element->NodalFunctions(basis,gauss);
240
241 /*Build B: */
242 for(int i=0;i<numnodes;i++){
243 B[numnodes*0+i] = basis[i];
244 B[numnodes*1+i] = basis[i];
245 }
246
247 /*Clean-up*/
248 xDelete<IssmDouble>(basis);
249}/*}}}*/
250void HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
251 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
252 * For node i, Bi' can be expressed in the actual coordinate system
253 * by:
254 * Bi_prime=[ dN/dx ]
255 * [ dN/dy ]
256 * where N is the finiteelement function for node i.
257 *
258 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
259 */
260
261 /*Fetch number of nodes for this finite element*/
262 int numnodes = element->GetNumberOfNodes();
263
264 /*Get nodal functions derivatives*/
265 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
266 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
267
268 /*Build B': */
269 for(int i=0;i<numnodes;i++){
270 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
271 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
272 }
273
274 /*Clean-up*/
275 xDelete<IssmDouble>(dbasis);
276
277}/*}}}*/
278void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
279 element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
280}/*}}}*/
281void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
282
283 /*Intermediary*/
284 int* doflist = NULL;
285
286 /*Fetch number of nodes for this finite element*/
287 int numnodes = element->GetNumberOfNodes();
288
289 /*Fetch dof list and allocate solution vector*/
290 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
291 IssmDouble* values = xNew<IssmDouble>(numnodes);
292
293 /*Use the dof list to index into the solution vector: */
294 for(int i=0;i<numnodes;i++){
295 values[i]=solution[doflist[i]];
296 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
297 if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
298 }
299
300 /*Add input to the element: */
301 element->AddInput(WatercolumnEnum,values,P1Enum);
302
303 /*Free ressources:*/
304 xDelete<IssmDouble>(values);
305 xDelete<int>(doflist);
306}/*}}}*/
307
308/*Intermediaries*/
309void HydrologyShreveAnalysis::CreateHydrologyWaterVelocityInput(Element* element){/*{{{*/
310
311 /*Intermediaries*/
312 IssmDouble dsdx,dsdy,dbdx,dbdy,w;
313
314 /*Retrieve all inputs and parameters*/
315 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
316 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
317 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum);
318 IssmDouble CR = element->GetMaterialParameter(HydrologyshreveCREnum);
319 IssmDouble n_man = element->GetMaterialParameter(HydrologyshreveNEnum);
320 IssmDouble mu_water = element->GetMaterialParameter(MaterialsMuWaterEnum);
321 Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
322 Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
323 Input* bedslopex_input = element->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
324 Input* bedslopey_input = element->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
325 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
326
327 /* compute VelocityFactor */
328 IssmDouble VelocityFactor = n_man*CR*CR*rho_water*g/mu_water;
329
330 /*Fetch number of vertices and allocate output*/
331 int numvertices = element->GetNumberOfVertices();
332 IssmDouble* vx = xNew<IssmDouble>(numvertices);
333 IssmDouble* vy = xNew<IssmDouble>(numvertices);
334
335 Gauss* gauss=element->NewGauss();
336 for(int iv=0;iv<numvertices;iv++){
337 gauss->GaussVertex(iv);
338 surfaceslopex_input->GetInputValue(&dsdx,gauss);
339 surfaceslopey_input->GetInputValue(&dsdy,gauss);
340 bedslopex_input->GetInputValue(&dbdx,gauss);
341 bedslopey_input->GetInputValue(&dbdy,gauss);
342 watercolumn_input->GetInputValue(&w,gauss);
343
344 /* Water velocity x and y components */
345 // vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
346 // vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
347 vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
348 vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
349 }
350
351 /*clean-up*/
352 delete gauss;
353
354 /*Add to inputs*/
355 element->AddInput(HydrologyWaterVxEnum,vx,P1Enum);
356 element->AddInput(HydrologyWaterVyEnum,vy,P1Enum);
357 xDelete<IssmDouble>(vx);
358 xDelete<IssmDouble>(vy);
359}/*}}}*/
Note: See TracBrowser for help on using the repository browser.