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

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

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