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

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

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