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

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

CHG: significant changes to the organization of slr -> new solidearth class
that encompasses all of the slr fields.

File size: 14.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*/
8void HydrologyShreveAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9
10 /*retrieve some parameters: */
11 int hydrology_model;
12 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
13
14 if(hydrology_model!=HydrologyshreveEnum) return;
15
16 IoModelToConstraintsx(constraints,iomodel,"md.hydrologyshreve.spcwatercolumn",HydrologyShreveAnalysisEnum,P1Enum);
17
18}/*}}}*/
19void HydrologyShreveAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
20 /*No loads*/
21}/*}}}*/
22void HydrologyShreveAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
23
24 /*Fetch parameters: */
25 int hydrology_model;
26 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
27
28 /*Now, do we really want Shreve?*/
29 if(hydrology_model!=HydrologyshreveEnum) return;
30
31 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
32 ::CreateNodes(nodes,iomodel,HydrologyShreveAnalysisEnum,P1Enum);
33 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
34}/*}}}*/
35int HydrologyShreveAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
36 return 1;
37}/*}}}*/
38void HydrologyShreveAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
39
40 /*Fetch data needed: */
41 int hydrology_model;
42 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
43
44 /*Now, do we really want Shreve?*/
45 if(hydrology_model!=HydrologyshreveEnum) return;
46
47 /*Update elements: */
48 int counter=0;
49 for(int i=0;i<iomodel->numberofelements;i++){
50 if(iomodel->my_elements[i]){
51 Element* element=(Element*)elements->GetObjectByOffset(counter);
52 element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
53 counter++;
54 }
55 }
56
57 iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
58 iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
59 iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
60 iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
61 if(iomodel->domaintype!=Domain2DhorizontalEnum){
62 iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
63 iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
64 }
65 iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
66 iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
67 iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
68 iomodel->FetchDataToInput(inputs2,elements,"md.initialization.watercolumn",WatercolumnEnum);
69
70 inputs2->DuplicateInput(WatercolumnEnum,WaterColumnOldEnum);
71}/*}}}*/
72void HydrologyShreveAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
73
74 /*retrieve some parameters: */
75 int hydrology_model;
76 int numoutputs;
77 char** requestedoutputs = NULL;
78 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
79
80 /*Now, do we really want Shreve?*/
81 if(hydrology_model!=HydrologyshreveEnum) return;
82
83 parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
84 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.stabilization",HydrologyshreveStabilizationEnum));
85 /*Requested outputs*/
86 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
87 parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
88 if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
89 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
90
91}/*}}}*/
92
93/*Finite Element Analysis*/
94void HydrologyShreveAnalysis::Core(FemModel* femmodel){/*{{{*/
95 _error_("not implemented");
96}/*}}}*/
97ElementVector* HydrologyShreveAnalysis::CreateDVector(Element* element){/*{{{*/
98 /*Default, return NULL*/
99 return NULL;
100}/*}}}*/
101void HydrologyShreveAnalysis::CreateHydrologyWaterVelocityInput(Element* element){/*{{{*/
102
103 /*Intermediaries*/
104 IssmDouble dsdx,dsdy,dbdx,dbdy,w;
105
106 /*Retrieve all inputs and parameters*/
107 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
108 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
109 IssmDouble g = element->FindParam(ConstantsGEnum);
110 IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum);
111 Input2* surfaceslopex_input = element->GetInput2(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
112 Input2* surfaceslopey_input = element->GetInput2(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
113 Input2* bedslopex_input = element->GetInput2(BedSlopeXEnum); _assert_(bedslopex_input);
114 Input2* bedslopey_input = element->GetInput2(BedSlopeYEnum); _assert_(bedslopey_input);
115 Input2* watercolumn_input = element->GetInput2(WatercolumnEnum); _assert_(watercolumn_input);
116
117 /*Fetch number of vertices and allocate output*/
118 int numvertices = element->GetNumberOfVertices();
119 IssmDouble* vx = xNew<IssmDouble>(numvertices);
120 IssmDouble* vy = xNew<IssmDouble>(numvertices);
121
122 Gauss* gauss=element->NewGauss();
123 for(int iv=0;iv<numvertices;iv++){
124 gauss->GaussVertex(iv);
125 surfaceslopex_input->GetInputValue(&dsdx,gauss);
126 surfaceslopey_input->GetInputValue(&dsdy,gauss);
127 bedslopex_input->GetInputValue(&dbdx,gauss);
128 bedslopey_input->GetInputValue(&dbdy,gauss);
129 watercolumn_input->GetInputValue(&w,gauss);
130
131 /* Water velocity x and y components */
132 vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
133 vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
134 }
135
136 /*clean-up*/
137 delete gauss;
138
139 /*Add to inputs*/
140 element->AddInput2(HydrologyWaterVxEnum,vx,P1Enum);
141 element->AddInput2(HydrologyWaterVyEnum,vy,P1Enum);
142 xDelete<IssmDouble>(vx);
143 xDelete<IssmDouble>(vy);
144}/*}}}*/
145ElementMatrix* HydrologyShreveAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
146_error_("Not implemented");
147}/*}}}*/
148ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/
149
150 /*Intermediaries */
151 IssmDouble diffusivity;
152 IssmDouble Jdet,D_scalar,dt,h;
153 IssmDouble vx,vy,vel,dvxdx,dvydy;
154 IssmDouble dvx[2],dvy[2];
155 IssmDouble* xyz_list = NULL;
156
157 /*Fetch number of nodes and dof for this finite element*/
158 int numnodes = element->GetNumberOfNodes();
159
160 /*Initialize Element vector and other vectors*/
161 ElementMatrix* Ke = element->NewElementMatrix();
162 IssmDouble* basis = xNew<IssmDouble>(numnodes);
163 IssmDouble* B = xNew<IssmDouble>(2*numnodes);
164 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);
165 IssmDouble D[2][2]={0.};
166
167 /*Create water velocity vx and vy from current inputs*/
168 CreateHydrologyWaterVelocityInput(element);
169
170 /*Retrieve all inputs and parameters*/
171 element->GetVerticesCoordinates(&xyz_list);
172 element->FindParam(&dt,TimesteppingTimeStepEnum);
173 element->FindParam(&diffusivity,HydrologyshreveStabilizationEnum);
174 Input2* vx_input=element->GetInput2(HydrologyWaterVxEnum); _assert_(vx_input);
175 Input2* vy_input=element->GetInput2(HydrologyWaterVyEnum); _assert_(vy_input);
176 h = element->CharacteristicLength();
177
178 /* Start looping on the number of gaussian points: */
179 Gauss* gauss=element->NewGauss(2);
180 for(int ig=gauss->begin();ig<gauss->end();ig++){
181 gauss->GaussPoint(ig);
182
183 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
184 element->NodalFunctions(basis,gauss);
185
186 vx_input->GetInputValue(&vx,gauss);
187 vy_input->GetInputValue(&vy,gauss);
188 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
189 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
190
191 D_scalar=gauss->weight*Jdet;
192
193 TripleMultiply(basis,1,numnodes,1,
194 &D_scalar,1,1,0,
195 basis,1,numnodes,0,
196 Ke->values,1);
197
198 GetB(B,element,xyz_list,gauss);
199 GetBprime(Bprime,element,xyz_list,gauss);
200
201 dvxdx=dvx[0];
202 dvydy=dvy[1];
203 D_scalar=dt*gauss->weight*Jdet;
204
205 D[0][0]=D_scalar*dvxdx;
206 D[1][1]=D_scalar*dvydy;
207 TripleMultiply(B,2,numnodes,1,
208 &D[0][0],2,2,0,
209 B,2,numnodes,0,
210 &Ke->values[0],1);
211
212 D[0][0]=D_scalar*vx;
213 D[1][1]=D_scalar*vy;
214 TripleMultiply(B,2,numnodes,1,
215 &D[0][0],2,2,0,
216 Bprime,2,numnodes,0,
217 &Ke->values[0],1);
218
219 /*Artificial diffusivity*/
220 vel=sqrt(vx*vx+vy*vy);
221 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
222 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
223 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
224 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
225 TripleMultiply(Bprime,2,numnodes,1,
226 &D[0][0],2,2,0,
227 Bprime,2,numnodes,0,
228 &Ke->values[0],1);
229 }
230
231 /*Clean up and return*/
232 xDelete<IssmDouble>(xyz_list);
233 xDelete<IssmDouble>(basis);
234 xDelete<IssmDouble>(B);
235 xDelete<IssmDouble>(Bprime);
236 delete gauss;
237 return Ke;
238}/*}}}*/
239ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
240
241 /*Skip if water or ice shelf element*/
242 if(element->IsFloating()) return NULL;
243
244 /*Intermediaries */
245 IssmDouble Jdet,dt;
246 IssmDouble mb,oldw;
247 IssmDouble* xyz_list = NULL;
248
249 /*Fetch number of nodes and dof for this finite element*/
250 int numnodes = element->GetNumberOfNodes();
251
252 /*Initialize Element vector and other vectors*/
253 ElementVector* pe = element->NewElementVector();
254 IssmDouble* basis = xNew<IssmDouble>(numnodes);
255
256 /*Retrieve all inputs and parameters*/
257 element->GetVerticesCoordinates(&xyz_list);
258 element->FindParam(&dt,TimesteppingTimeStepEnum);
259 Input2* mb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
260 Input2* oldw_input = element->GetInput2(WaterColumnOldEnum); _assert_(oldw_input);
261
262 /*Initialize mb_correction to 0, do not forget!:*/
263 /* Start looping on the number of gaussian points: */
264 Gauss* gauss=element->NewGauss(2);
265 for(int ig=gauss->begin();ig<gauss->end();ig++){
266 gauss->GaussPoint(ig);
267
268 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
269 element->NodalFunctions(basis,gauss);
270
271 mb_input->GetInputValue(&mb,gauss);
272 oldw_input->GetInputValue(&oldw,gauss);
273
274 if(dt!=0.){
275 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
276 }
277 else{
278 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
279 }
280 }
281
282 /*Clean up and return*/
283 xDelete<IssmDouble>(xyz_list);
284 xDelete<IssmDouble>(basis);
285 delete gauss;
286 return pe;
287}/*}}}*/
288void HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
289 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2.
290 * For node i, Bi can be expressed in the actual coordinate system
291 * by:
292 * Bi=[ N ]
293 * [ N ]
294 * where N is the finiteelement function for node i.
295 *
296 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
297 */
298
299 /*Fetch number of nodes for this finite element*/
300 int numnodes = element->GetNumberOfNodes();
301
302 /*Get nodal functions*/
303 IssmDouble* basis=xNew<IssmDouble>(numnodes);
304 element->NodalFunctions(basis,gauss);
305
306 /*Build B: */
307 for(int i=0;i<numnodes;i++){
308 B[numnodes*0+i] = basis[i];
309 B[numnodes*1+i] = basis[i];
310 }
311
312 /*Clean-up*/
313 xDelete<IssmDouble>(basis);
314}/*}}}*/
315void HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
316 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
317 * For node i, Bi' can be expressed in the actual coordinate system
318 * by:
319 * Bi_prime=[ dN/dx ]
320 * [ dN/dy ]
321 * where N is the finiteelement function for node i.
322 *
323 * We assume B' has been allocated already, of size: 3x(2*numnodes)
324 */
325
326 /*Fetch number of nodes for this finite element*/
327 int numnodes = element->GetNumberOfNodes();
328
329 /*Get nodal functions derivatives*/
330 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
331 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
332
333 /*Build B': */
334 for(int i=0;i<numnodes;i++){
335 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
336 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
337 }
338
339 /*Clean-up*/
340 xDelete<IssmDouble>(dbasis);
341
342}/*}}}*/
343void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
344 element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
345}/*}}}*/
346void HydrologyShreveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
347 _error_("Not implemented yet");
348}/*}}}*/
349void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
350
351 /*Intermediary*/
352 int* doflist = NULL;
353
354 /*Fetch number of nodes for this finite element*/
355 int numnodes = element->GetNumberOfNodes();
356
357 /*Fetch dof list and allocate solution vector*/
358 element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
359 IssmDouble* values = xNew<IssmDouble>(numnodes);
360
361 /*Use the dof list to index into the solution vector: */
362 for(int i=0;i<numnodes;i++){
363 values[i]=solution[doflist[i]];
364 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
365 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
366 if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
367 }
368
369 /*Add input to the element: */
370 element->AddInput2(WatercolumnEnum,values,element->GetElementType());
371
372 /*Free ressources:*/
373 xDelete<IssmDouble>(values);
374 xDelete<int>(doflist);
375}/*}}}*/
376void HydrologyShreveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
377 /*Default, do nothing*/
378 return;
379}/*}}}*/
380
381/*Needed changes to switch to the Johnson formulation*//*{{{*/
382/*All the changes are to be done in the velocity computation.
383 The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
384 'p' and 'q' which are the exponent of the Manning formula for laminar (p=2,q=1) or turbulent (p=2/3,q=1/2) flow
385 'R' the hydraulic radius
386 'n' the manning roughness coeficient
387
388 With these, the velocity reads ;
389
390 v= - (1/n)* pow(R,p)*pow((grad phi(rho_water*g)),q)
391
392 you should also redefine the water pressure potential 'phi' with respect to the effective pressure deffinition given in Johson:
393 phi=(rho_ice*g*( surface + ((rho_water/rho_ice)-1)*base - k_n*((thickness* grad(base))/omega) )
394
395 where
396 'omega' is the fractional area of the base occupied by the water film
397 'k_n' is a parameter
398 This last equation derives from the effective pressure definition developped in Alley 1989
399*/
400/*}}}*/
Note: See TracBrowser for help on using the repository browser.