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

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

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