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

Last change on this file since 27794 was 27794, checked in by Mathieu Morlighem, 22 months ago

BUG: was not pulling the right field
File size: 13.6 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.hydrology.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,Inputs* inputs,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(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
53 counter++;
54 }
55 }
56
57 iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
58 iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
59 iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
60 iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
61 if(iomodel->domaintype!=Domain2DhorizontalEnum){
62 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
63 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
64 }
65 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
66 iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
67 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
68 iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum);
69
70 inputs->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}/*}}}*/
97void HydrologyShreveAnalysis::PreCore(FemModel* femmodel){/*{{{*/
98 _error_("not implemented");
99}/*}}}*/
100ElementVector* HydrologyShreveAnalysis::CreateDVector(Element* element){/*{{{*/
101 /*Default, return NULL*/
102 return NULL;
103}/*}}}*/
104void HydrologyShreveAnalysis::CreateHydrologyWaterVelocityInput(Element* element){/*{{{*/
105
106 /*Intermediaries*/
107 IssmDouble dsdx,dsdy,dbdx,dbdy,w;
108
109 /*Retrieve all inputs and parameters*/
110 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
111 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
112 IssmDouble g = element->FindParam(ConstantsGEnum);
113 IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum);
114 Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
115 Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
116 Input* bedslopex_input = element->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
117 Input* bedslopey_input = element->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
118 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
119
120 /*Fetch number of vertices and allocate output*/
121 int numvertices = element->GetNumberOfVertices();
122 IssmDouble* vx = xNew<IssmDouble>(numvertices);
123 IssmDouble* vy = xNew<IssmDouble>(numvertices);
124
125 Gauss* gauss=element->NewGauss();
126 for(int iv=0;iv<numvertices;iv++){
127 gauss->GaussVertex(iv);
128 surfaceslopex_input->GetInputValue(&dsdx,gauss);
129 surfaceslopey_input->GetInputValue(&dsdy,gauss);
130 bedslopex_input->GetInputValue(&dbdx,gauss);
131 bedslopey_input->GetInputValue(&dbdy,gauss);
132 watercolumn_input->GetInputValue(&w,gauss);
133
134 /* Water velocity x and y components */
135 vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
136 vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
137 }
138
139 /*clean-up*/
140 delete gauss;
141
142 /*Add to inputs*/
143 element->AddInput(HydrologyWaterVxEnum,vx,P1Enum);
144 element->AddInput(HydrologyWaterVyEnum,vy,P1Enum);
145 xDelete<IssmDouble>(vx);
146 xDelete<IssmDouble>(vy);
147}/*}}}*/
148ElementMatrix* HydrologyShreveAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
149_error_("Not implemented");
150}/*}}}*/
151ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/
152
153 /*Intermediaries */
154 IssmDouble diffusivity;
155 IssmDouble Jdet,D_scalar,dt,h;
156 IssmDouble vx,vy,vel,dvxdx,dvydy;
157 IssmDouble dvx[2],dvy[2];
158 IssmDouble* xyz_list = NULL;
159
160 /*Fetch number of nodes and dof for this finite element*/
161 int numnodes = element->GetNumberOfNodes();
162
163 /*Initialize Element vector and other vectors*/
164 ElementMatrix* Ke = element->NewElementMatrix();
165 IssmDouble* basis = xNew<IssmDouble>(numnodes);
166 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
167 IssmDouble D[2][2]={0.};
168
169 /*Create water velocity vx and vy from current inputs*/
170 CreateHydrologyWaterVelocityInput(element);
171
172 /*Retrieve all inputs and parameters*/
173 element->GetVerticesCoordinates(&xyz_list);
174 element->FindParam(&dt,TimesteppingTimeStepEnum);
175 element->FindParam(&diffusivity,HydrologyshreveStabilizationEnum);
176 Input* vx_input=element->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
177 Input* vy_input=element->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
178 h = element->CharacteristicLength();
179
180 /* Start looping on the number of gaussian points: */
181 Gauss* gauss=element->NewGauss(2);
182 while(gauss->next()){
183
184 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
185 element->NodalFunctions(basis,gauss);
186 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
187
188 vx_input->GetInputValue(&vx,gauss);
189 vy_input->GetInputValue(&vy,gauss);
190 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
191 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
192
193 /*Transient term*/
194 D_scalar=gauss->weight*Jdet;
195 for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
196
197 /*Advection terms*/
198 dvxdx=dvx[0];
199 dvydy=dvy[1];
200 D_scalar=dt*gauss->weight*Jdet;
201 for(int i=0;i<numnodes;i++){
202 for(int j=0;j<numnodes;j++){
203 /*\phi_i \phi_j \nabla\cdot v*/
204 Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
205 /*\phi_i v\cdot\nabla\phi_j*/
206 Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
207 }
208 }
209
210 /*Artificial diffusivity*/
211 vel=sqrt(vx*vx+vy*vy);
212 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
213 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
214 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
215 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
216 for(int i=0;i<numnodes;i++){
217 for(int j=0;j<numnodes;j++){
218 Ke->values[i*numnodes+j] += (
219 dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) +
220 dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j])
221 );
222 }
223 }
224 }
225
226 /*Clean up and return*/
227 xDelete<IssmDouble>(xyz_list);
228 xDelete<IssmDouble>(basis);
229 xDelete<IssmDouble>(dbasis);
230 delete gauss;
231 return Ke;
232}/*}}}*/
233ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
234
235 /*Skip if water or ice shelf element*/
236 if(element->IsAllFloating()) return NULL;
237
238 /*Intermediaries */
239 IssmDouble Jdet,dt;
240 IssmDouble mb,oldw;
241 IssmDouble* xyz_list = NULL;
242
243 /*Fetch number of nodes and dof for this finite element*/
244 int numnodes = element->GetNumberOfNodes();
245
246 /*Initialize Element vector and other vectors*/
247 ElementVector* pe = element->NewElementVector();
248 IssmDouble* basis = xNew<IssmDouble>(numnodes);
249
250 /*Retrieve all inputs and parameters*/
251 element->GetVerticesCoordinates(&xyz_list);
252 element->FindParam(&dt,TimesteppingTimeStepEnum);
253 Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
254 Input* oldw_input = element->GetInput(WaterColumnOldEnum); _assert_(oldw_input);
255
256 /*Initialize mb_correction to 0, do not forget!:*/
257 /* Start looping on the number of gaussian points: */
258 Gauss* gauss=element->NewGauss(2);
259 while(gauss->next()){
260
261 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
262 element->NodalFunctions(basis,gauss);
263
264 mb_input->GetInputValue(&mb,gauss);
265 oldw_input->GetInputValue(&oldw,gauss);
266
267 if(dt!=0.){
268 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
269 }
270 else{
271 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
272 }
273 }
274
275 /*Clean up and return*/
276 xDelete<IssmDouble>(xyz_list);
277 xDelete<IssmDouble>(basis);
278 delete gauss;
279 return pe;
280}/*}}}*/
281void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
282 element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
283}/*}}}*/
284void HydrologyShreveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
285 _error_("Not implemented yet");
286}/*}}}*/
287void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
288
289 /*Intermediary*/
290 int* doflist = NULL;
291
292 /*Fetch number of nodes for this finite element*/
293 int numnodes = element->GetNumberOfNodes();
294
295 /*Fetch dof list and allocate solution vector*/
296 element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
297 IssmDouble* watercolumn = xNew<IssmDouble>(numnodes);
298
299 /*Use the dof list to index into the solution vector: */
300 for(int i=0;i<numnodes;i++){
301 watercolumn[i]=solution[doflist[i]];
302 if(xIsNan<IssmDouble>(watercolumn[i])) _error_("NaN found in solution vector");
303 if(xIsInf<IssmDouble>(watercolumn[i])) _error_("Inf found in solution vector");
304 if (watercolumn[i]<10e-10) watercolumn[i]=10e-10; //correcting the water column to positive watercolumn
305 }
306
307 /*Add input to the element: */
308 element->AddInput(WatercolumnEnum,watercolumn,element->GetElementType());
309
310 /*Also update the hydrological loads for the sealevel core: */
311 IssmDouble* oldwatercolumn = xNew<IssmDouble>(numnodes);
312 IssmDouble* deltawatercolumn = xNew<IssmDouble>(numnodes);
313
314 element->GetInputListOnVertices(&watercolumn[0],WatercolumnEnum);
315 element->GetInputListOnVertices(&oldwatercolumn[0],WaterColumnOldEnum);
316 element->GetInputListOnVertices(&deltawatercolumn[0],AccumulatedDeltaTwsEnum);
317 for(int i=0;i<numnodes;i++){
318 deltawatercolumn[i] += watercolumn[i]-oldwatercolumn[i];
319 }
320 element->AddInput(AccumulatedDeltaTwsEnum,deltawatercolumn,P1Enum);
321
322 /*Free resources:*/
323 xDelete<IssmDouble>(oldwatercolumn);
324 xDelete<IssmDouble>(deltawatercolumn);
325 xDelete<IssmDouble>(watercolumn);
326 xDelete<int>(doflist);
327}/*}}}*/
328void HydrologyShreveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
329 /*Default, do nothing*/
330 return;
331}/*}}}*/
332
333/*Needed changes to switch to the Johnson formulation*//*{{{*/
334/*All the changes are to be done in the velocity computation.
335 The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
336 '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
337 'R' the hydraulic radius
338 'n' the manning roughness coeficient
339
340 With these, the velocity reads ;
341
342 v= - (1/n)* pow(R,p)*pow((grad phi(rho_water*g)),q)
343
344 you should also redefine the water pressure potential 'phi' with respect to the effective pressure deffinition given in Johson:
345 phi=(rho_ice*g*( surface + ((rho_water/rho_ice)-1)*base - k_n*((thickness* grad(base))/omega) )
346
347 where
348 'omega' is the fractional area of the base occupied by the water film
349 'k_n' is a parameter
350 This last equation derives from the effective pressure definition developped in Alley 1989
351*/
352/*}}}*/
Note: See TracBrowser for help on using the repository browser.