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

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

CHG: fixing issues with elastic vs rigid vs rotation.
File size: 12.8 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,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.solidearth.initialsealevel",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}/*}}}*/
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 Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
112 Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
113 Input* bedslopex_input = element->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
114 Input* bedslopey_input = element->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
115 Input* watercolumn_input = element->GetInput(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->AddInput(HydrologyWaterVxEnum,vx,P1Enum);
141 element->AddInput(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* dbasis = xNew<IssmDouble>(2*numnodes);
164 IssmDouble D[2][2]={0.};
165
166 /*Create water velocity vx and vy from current inputs*/
167 CreateHydrologyWaterVelocityInput(element);
168
169 /*Retrieve all inputs and parameters*/
170 element->GetVerticesCoordinates(&xyz_list);
171 element->FindParam(&dt,TimesteppingTimeStepEnum);
172 element->FindParam(&diffusivity,HydrologyshreveStabilizationEnum);
173 Input* vx_input=element->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
174 Input* vy_input=element->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
175 h = element->CharacteristicLength();
176
177 /* Start looping on the number of gaussian points: */
178 Gauss* gauss=element->NewGauss(2);
179 while(gauss->next()){
180
181 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
182 element->NodalFunctions(basis,gauss);
183 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
184
185 vx_input->GetInputValue(&vx,gauss);
186 vy_input->GetInputValue(&vy,gauss);
187 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
188 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
189
190 /*Transient term*/
191 D_scalar=gauss->weight*Jdet;
192 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];
193
194 /*Advection terms*/
195 dvxdx=dvx[0];
196 dvydy=dvy[1];
197 D_scalar=dt*gauss->weight*Jdet;
198 for(int i=0;i<numnodes;i++){
199 for(int j=0;j<numnodes;j++){
200 /*\phi_i \phi_j \nabla\cdot v*/
201 Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
202 /*\phi_i v\cdot\nabla\phi_j*/
203 Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
204 }
205 }
206
207 /*Artificial diffusivity*/
208 vel=sqrt(vx*vx+vy*vy);
209 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
210 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
211 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
212 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
213 for(int i=0;i<numnodes;i++){
214 for(int j=0;j<numnodes;j++){
215 Ke->values[i*numnodes+j] += (
216 dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) +
217 dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j])
218 );
219 }
220 }
221 }
222
223 /*Clean up and return*/
224 xDelete<IssmDouble>(xyz_list);
225 xDelete<IssmDouble>(basis);
226 xDelete<IssmDouble>(dbasis);
227 delete gauss;
228 return Ke;
229}/*}}}*/
230ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
231
232 /*Skip if water or ice shelf element*/
233 if(element->IsFloating()) return NULL;
234
235 /*Intermediaries */
236 IssmDouble Jdet,dt;
237 IssmDouble mb,oldw;
238 IssmDouble* xyz_list = NULL;
239
240 /*Fetch number of nodes and dof for this finite element*/
241 int numnodes = element->GetNumberOfNodes();
242
243 /*Initialize Element vector and other vectors*/
244 ElementVector* pe = element->NewElementVector();
245 IssmDouble* basis = xNew<IssmDouble>(numnodes);
246
247 /*Retrieve all inputs and parameters*/
248 element->GetVerticesCoordinates(&xyz_list);
249 element->FindParam(&dt,TimesteppingTimeStepEnum);
250 Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
251 Input* oldw_input = element->GetInput(WaterColumnOldEnum); _assert_(oldw_input);
252
253 /*Initialize mb_correction to 0, do not forget!:*/
254 /* Start looping on the number of gaussian points: */
255 Gauss* gauss=element->NewGauss(2);
256 while(gauss->next()){
257
258 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
259 element->NodalFunctions(basis,gauss);
260
261 mb_input->GetInputValue(&mb,gauss);
262 oldw_input->GetInputValue(&oldw,gauss);
263
264 if(dt!=0.){
265 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
266 }
267 else{
268 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
269 }
270 }
271
272 /*Clean up and return*/
273 xDelete<IssmDouble>(xyz_list);
274 xDelete<IssmDouble>(basis);
275 delete gauss;
276 return pe;
277}/*}}}*/
278void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
279 element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
280}/*}}}*/
281void HydrologyShreveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
282 _error_("Not implemented yet");
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->GetDofListLocal(&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(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
301 if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
302 }
303
304 /*Add input to the element: */
305 element->AddInput(WatercolumnEnum,values,element->GetElementType());
306
307 /*Free ressources:*/
308 xDelete<IssmDouble>(values);
309 xDelete<int>(doflist);
310}/*}}}*/
311void HydrologyShreveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
312 /*Default, do nothing*/
313 return;
314}/*}}}*/
315
316/*Needed changes to switch to the Johnson formulation*//*{{{*/
317/*All the changes are to be done in the velocity computation.
318 The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
319 '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
320 'R' the hydraulic radius
321 'n' the manning roughness coeficient
322
323 With these, the velocity reads ;
324
325 v= - (1/n)* pow(R,p)*pow((grad phi(rho_water*g)),q)
326
327 you should also redefine the water pressure potential 'phi' with respect to the effective pressure deffinition given in Johson:
328 phi=(rho_ice*g*( surface + ((rho_water/rho_ice)-1)*base - k_n*((thickness* grad(base))/omega) )
329
330 where
331 'omega' is the fractional area of the base occupied by the water film
332 'k_n' is a parameter
333 This last equation derives from the effective pressure definition developped in Alley 1989
334*/
335/*}}}*/
Note: See TracBrowser for help on using the repository browser.