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

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

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