Ice Sheet System Model  4.18
Code documentation
HydrologyShaktiAnalysis.cpp
Go to the documentation of this file.
1 #include "./HydrologyShaktiAnalysis.h"
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 
7 /*Define 2 hardcoded parameters*/
8 #define OMEGA 0.001 // parameter controlling transition to nonlinear resistance in basal system (dimensionless)
9 #define NU 1.787e-6 //kinematic water viscosity m^2/s
10 #define CT 7.5e-8 // Clapeyron slope (K/Pa)
11 #define CW 4.22e3 // specific heat capacity of water (J/kg/K)
12 
13 /*Model processing*/
14 void HydrologyShaktiAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
15 
16  /*retrieve some parameters: */
17  int hydrology_model;
18  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
19 
20  if(hydrology_model!=HydrologyshaktiEnum) return;
21 
22  IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologyShaktiAnalysisEnum,P1Enum);
23 
24 }/*}}}*/
25 void HydrologyShaktiAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
26 
27  /*Fetch parameters: */
28  int hydrology_model;
29  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
30 
31  /*Now, do we really want Shakti?*/
32  if(hydrology_model!=HydrologyshaktiEnum) return;
33 
34  /*Create discrete loads for Moulins*/
35  CreateSingleNodeToElementConnectivity(iomodel);
36  for(int i=0;i<iomodel->numberofvertices;i++){
37  if (iomodel->domaintype!=Domain3DEnum){
38  /*keep only this partition's nodes:*/
39  if(iomodel->my_vertices[i]){
40  loads->AddObject(new Moulin(i+1,i,iomodel));
41  }
42  }
43  else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
44  if(iomodel->my_vertices[i]){
45  loads->AddObject(new Moulin(i+1,i,iomodel));
46  }
47  }
48  }
49  iomodel->DeleteData(1,"md.mesh.vertexonbase");
50 
51  /*Deal with Neumann BC*/
52  int M,N;
53  int *segments = NULL;
54  iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
55 
56  /*Check that the size seem right*/
57  _assert_(N==3); _assert_(M>=3);
58  for(int i=0;i<M;i++){
59  if(iomodel->my_elements[segments[i*3+2]-1]){
60  loads->AddObject(new Neumannflux(i+1,i,iomodel,segments));
61  }
62  }
63  xDelete<int>(segments);
64 
65 }/*}}}*/
66 void HydrologyShaktiAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
67 
68  /*Fetch parameters: */
69  int hydrology_model;
70  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
71 
72  /*Now, do we really want Shakti?*/
73  if(hydrology_model!=HydrologyshaktiEnum) return;
74 
75  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
76  ::CreateNodes(nodes,iomodel,HydrologyShaktiAnalysisEnum,P1Enum);
77  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
78 }/*}}}*/
79 int HydrologyShaktiAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
80  return 1;
81 }/*}}}*/
82 void HydrologyShaktiAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
83 
84  /*Fetch data needed: */
85  int hydrology_model,frictionlaw;
86  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
87 
88  /*Now, do we really want Shakti?*/
89  if(hydrology_model!=HydrologyshaktiEnum) return;
90 
91  /*Update elements: */
92  int counter=0;
93  for(int i=0;i<iomodel->numberofelements;i++){
94  if(iomodel->my_elements[i]){
95  Element* element=(Element*)elements->GetObjectByOffset(counter);
96  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
97  counter++;
98  }
99  }
100 
101  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
102  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
103  if(iomodel->domaintype!=Domain2DhorizontalEnum){
104  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
105  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
106  }
107  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
108  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
109  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
110  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
111  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.head",HydrologyHeadEnum);
112  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.gap_height",HydrologyGapHeightEnum);
113  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.englacial_input",HydrologyEnglacialInputEnum);
114  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
115  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.bump_spacing",HydrologyBumpSpacingEnum);
116  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
117  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.reynolds",HydrologyReynoldsEnum);
118  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
119  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
120  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
121  iomodel->FindConstant(&frictionlaw,"md.friction.law");
122 
123  /*Friction law variables*/
124  switch(frictionlaw){
125  case 1:
126  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
127  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
128  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
129  break;
130  case 8:
131  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
132  break;
133  default:
134  _error_("Friction law "<< frictionlaw <<" not supported");
135  }
136 }/*}}}*/
137 void HydrologyShaktiAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
138 
139  /*retrieve some parameters: */
140  int hydrology_model;
141  int numoutputs;
142  char** requestedoutputs = NULL;
143  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
144 
145  /*Now, do we really want Shakti?*/
146  if(hydrology_model!=HydrologyshaktiEnum) return;
147 
148  parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
149  parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
150  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
151  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
152 
153  /*Requested outputs*/
154  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
155  parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
156  if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
157  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
158 }/*}}}*/
159 
160 /*Finite Element Analysis*/
161 void HydrologyShaktiAnalysis::Core(FemModel* femmodel){/*{{{*/
162  _error_("not implemented");
163 }/*}}}*/
164 ElementVector* HydrologyShaktiAnalysis::CreateDVector(Element* element){/*{{{*/
165  /*Default, return NULL*/
166  return NULL;
167 }/*}}}*/
168 ElementMatrix* HydrologyShaktiAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
169 _error_("Not implemented");
170 }/*}}}*/
171 ElementMatrix* HydrologyShaktiAnalysis::CreateKMatrix(Element* element){/*{{{*/
172 
173  /*Intermediaries */
174  IssmDouble Jdet;
175  IssmDouble* xyz_list = NULL;
176 
177  /*Fetch number of nodes and dof for this finite element*/
178  int numnodes = element->GetNumberOfNodes();
179 
180  /*Initialize Element vector and other vectors*/
181  ElementMatrix* Ke = element->NewElementMatrix();
182  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
183  IssmDouble* basis = xNew<IssmDouble>(numnodes);
184 
185  /*Retrieve all inputs and parameters*/
186  element->GetVerticesCoordinates(&xyz_list);
187 
188  /*Get conductivity from inputs*/
189  IssmDouble conductivity = GetConductivity(element);
190 
191  /*Get englacial storage coefficient*/
192  IssmDouble storage,dt;
193  element->FindParam(&storage,HydrologyStorageEnum);
194  element->FindParam(&dt,TimesteppingTimeStepEnum);
195 
196  /* Start looping on the number of gaussian points: */
197  Gauss* gauss=element->NewGauss(1);
198  for(int ig=gauss->begin();ig<gauss->end();ig++){
199  gauss->GaussPoint(ig);
200 
201  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
202  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
203  element->NodalFunctions(basis,gauss);
204 
205  for(int i=0;i<numnodes;i++){
206  for(int j=0;j<numnodes;j++){
207  Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
208  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j];
209  }
210  }
211  }
212 
213  /*Clean up and return*/
214  xDelete<IssmDouble>(xyz_list);
215  xDelete<IssmDouble>(basis);
216  xDelete<IssmDouble>(dbasis);
217  delete gauss;
218  return Ke;
219 }/*}}}*/
220 ElementVector* HydrologyShaktiAnalysis::CreatePVector(Element* element){/*{{{*/
221 
222  /*Skip if water or ice shelf element*/
223  if(element->IsFloating()) return NULL;
224 
225  /*Intermediaries */
226  IssmDouble Jdet,meltrate,G,dh[2],B,A,n;
227  IssmDouble gap,bed,thickness,head,ieb,head_old;
228  IssmDouble lr,br,vx,vy,beta;
229  IssmDouble alpha2,frictionheat;
230  IssmDouble PMPheat,dpressure_water[2],dbed[2];
231  IssmDouble* xyz_list = NULL;
232 
233  /*Fetch number of nodes and dof for this finite element*/
234  int numnodes = element->GetNumberOfNodes();
235 
236  /*Initialize Element vector and other vectors*/
237  ElementVector* pe = element->NewElementVector();
238  IssmDouble* basis = xNew<IssmDouble>(numnodes);
239 
240  /*Retrieve all inputs and parameters*/
241  element->GetVerticesCoordinates(&xyz_list);
242  IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum);
243  IssmDouble g = element->FindParam(ConstantsGEnum);
244  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
245  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
246  Input2* geothermalflux_input = element->GetInput2(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
247  Input2* head_input = element->GetInput2(HydrologyHeadEnum); _assert_(head_input);
248  Input2* gap_input = element->GetInput2(HydrologyGapHeightEnum); _assert_(gap_input);
249  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
250  Input2* base_input = element->GetInput2(BaseEnum); _assert_(base_input);
251  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
252  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
253  Input2* englacial_input = element->GetInput2(HydrologyEnglacialInputEnum); _assert_(englacial_input);
254  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
255  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
256  Input2* lr_input = element->GetInput2(HydrologyBumpSpacingEnum); _assert_(lr_input);
257  Input2* br_input = element->GetInput2(HydrologyBumpHeightEnum); _assert_(br_input);
258  Input2* headold_input = element->GetInput2(HydrologyHeadOldEnum); _assert_(headold_input);
259 
260  /*Get conductivity from inputs*/
261  IssmDouble conductivity = GetConductivity(element);
262 
263  /*Get englacial storage coefficient*/
264  IssmDouble storage,dt;
265  element->FindParam(&storage,HydrologyStorageEnum);
266  element->FindParam(&dt,TimesteppingTimeStepEnum);
267 
268  /*Build friction element, needed later: */
269  Friction* friction=new Friction(element,2);
270 
271  /* Start looping on the number of gaussian points: */
272  Gauss* gauss=element->NewGauss(2);
273  for(int ig=gauss->begin();ig<gauss->end();ig++){
274  gauss->GaussPoint(ig);
275 
276  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
277  element->NodalFunctions(basis,gauss);
278  geothermalflux_input->GetInputValue(&G,gauss);
279  base_input->GetInputValue(&bed,gauss);
280  base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
281  thickness_input->GetInputValue(&thickness,gauss);
282  gap_input->GetInputValue(&gap,gauss);
283  head_input->GetInputValue(&head,gauss);
284  head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
285  englacial_input->GetInputValue(&ieb,gauss);
286  lr_input->GetInputValue(&lr,gauss);
287  br_input->GetInputValue(&br,gauss);
288  vx_input->GetInputValue(&vx,gauss);
289  vy_input->GetInputValue(&vy,gauss);
290  headold_input->GetInputValue(&head_old,gauss);
291 
292  /*Get ice A parameter*/
293  B_input->GetInputValue(&B,gauss);
294  n_input->GetInputValue(&n,gauss);
295  A=pow(B,-n);
296 
297  /*Compute beta term*/
298  if(gap<br)
299  beta = (br-gap)/lr;
300  else
301  beta = 0.;
302 
303  /*Compute frictional heat flux*/
304  friction->GetAlpha2(&alpha2,gauss);
305  vx_input->GetInputValue(&vx,gauss);
306  vy_input->GetInputValue(&vy,gauss);
307  frictionheat=alpha2*(vx*vx+vy*vy);
308 
309  /*Get water and ice pressures*/
310  IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.);
311  IssmDouble pressure_water = rho_water*g*(head-bed);
312  if(pressure_water>pressure_ice) pressure_water = pressure_ice;
313 
314  /*Get water pressure from previous time step to use in lagged creep term*/
315  IssmDouble pressure_water_old = rho_water*g*(head_old-bed);
316  if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice;
317 
318  /*Compute change in sensible heat due to changes in pressure melting point*/
319  dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
320  dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
321  PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
322 
323  meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
324  _assert_(meltrate>0.);
325 
326  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
327  (
328  meltrate*(1/rho_water-1/rho_ice)
329  +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
330  -beta*sqrt(vx*vx+vy*vy)
331  +ieb
332  +storage*head_old/dt
333  )*basis[i];
334  }
335  /*Clean up and return*/
336  xDelete<IssmDouble>(xyz_list);
337  xDelete<IssmDouble>(basis);
338  delete friction;
339  delete gauss;
340  return pe;
341 }/*}}}*/
342 void HydrologyShaktiAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
343  element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum);
344 }/*}}}*/
345 void HydrologyShaktiAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
346  _error_("Not implemented yet");
347 }/*}}}*/
348 void HydrologyShaktiAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
349 
350  /*Intermediary*/
351  IssmDouble dh[3];
352  int* doflist = NULL;
353  IssmDouble* xyz_list = NULL;
354 
355  /*Get gravity from parameters*/
356  IssmDouble g = element->FindParam(ConstantsGEnum);
357 
358  /*Fetch number of nodes for this finite element*/
359  int numnodes = element->GetNumberOfNodes();
360 
361  /*Fetch dof list and allocate solution vector*/
362  element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
363  IssmDouble* values = xNew<IssmDouble>(numnodes);
364 
365  /*Get thickness and base on nodes to apply cap on water head*/
366  IssmDouble* eff_pressure = xNew<IssmDouble>(numnodes);
367  IssmDouble* thickness = xNew<IssmDouble>(numnodes);
368  IssmDouble* bed = xNew<IssmDouble>(numnodes);
369  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
370  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
371  element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
372  element->GetInputListOnNodes(&bed[0],BaseEnum);
373 
374  /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/
375  IssmDouble* head_old = xNew<IssmDouble>(numnodes);
376  element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum);
377  IssmDouble relaxation;
378  element->FindParam(&relaxation,HydrologyRelaxationEnum);
379 
380  /*Use the dof list to index into the solution vector: */
381  for(int i=0;i<numnodes;i++){
382  values[i]=solution[doflist[i]];
383 
384  /*make sure that p_water<p_ice -> h<rho_i H/rho_w + zb*/
385  if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
386  values[i] = rho_ice*thickness[i]/rho_water+bed[i];
387  }
388 
389  /*Make sure that negative pressure is not allowed*/
390  // if(values[i]<bed[i]){
391  // values[i] = bed[i];
392  // }
393 
394  /*Under-relaxation*/
395  values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
396 
397  /*Calculate effective pressure*/
398  eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
399 
400  if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
401  if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
402  }
403 
404  /*Add input to the element: */
405  element->AddInput2(HydrologyHeadEnum,values,element->GetElementType());
406  element->AddInput2(EffectivePressureEnum,eff_pressure,P1Enum);
407 
408  /*Update reynolds number according to new solution*/
409  element->GetVerticesCoordinates(&xyz_list);
410  Input2* head_input = element->GetInput2(HydrologyHeadEnum);_assert_(head_input);
411  IssmDouble conductivity = GetConductivity(element);
412 
413  /*Get gap height derivatives at the center of the element*/
414  Gauss* gauss=element->NewGauss(1);
415  head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
416  delete gauss;
417 
418  IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
419  element->AddInput2(HydrologyReynoldsEnum,&reynolds,P0Enum);
420 
421  /*Free resources:*/
422  xDelete<IssmDouble>(values);
423  xDelete<IssmDouble>(thickness);
424  xDelete<IssmDouble>(bed);
425  xDelete<IssmDouble>(xyz_list);
426  xDelete<int>(doflist);
427  xDelete<IssmDouble>(eff_pressure);
428  xDelete<IssmDouble>(head_old);
429 }/*}}}*/
430 void HydrologyShaktiAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
431  /*Default, do nothing*/
432  return;
433 }/*}}}*/
434 
435 /*Additional methods*/
436 IssmDouble HydrologyShaktiAnalysis::GetConductivity(Element* element){/*{{{*/
437 
438  /*Intermediaries */
439  IssmDouble gap,reynolds;
440 
441  /*Get gravity from parameters*/
442  IssmDouble g = element->FindParam(ConstantsGEnum);
443 
444  /*Get Reynolds and gap average values*/
445  Input2* reynolds_input = element->GetInput2(HydrologyReynoldsEnum); _assert_(reynolds_input);
446  Input2* gap_input = element->GetInput2(HydrologyGapHeightEnum); _assert_(gap_input);
447  reynolds_input->GetInputAverage(&reynolds);
448  gap_input->GetInputAverage(&gap);
449 
450  /*Compute conductivity*/
451  IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
452  _assert_(conductivity>0);
453 
454  /*Clean up and return*/
455  return conductivity;
456 }/*}}}*/
457 void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/
458 
459  for(int j=0;j<femmodel->elements->Size();j++){
460  Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
461  UpdateGapHeight(element);
462  }
463 
464 }/*}}}*/
465 void HydrologyShaktiAnalysis::UpdateGapHeight(Element* element){/*{{{*/
466 
467  /*Skip if water or ice shelf element*/
468  if(element->IsFloating()) return;
469 
470  /*Intermediaries */
471  IssmDouble newgap = 0.;
472  IssmDouble Jdet,meltrate,G,dh[2],B,A,n,dt;
473  IssmDouble gap,bed,thickness,head,ieb;
474  IssmDouble lr,br,vx,vy,beta;
475  IssmDouble alpha2,frictionheat;
476  IssmDouble* xyz_list = NULL;
477  IssmDouble dpressure_water[2],dbed[2],PMPheat;
478  IssmDouble q = 0.;
479  IssmDouble channelization = 0.;
480 
481  /*Retrieve all inputs and parameters*/
482  element->GetVerticesCoordinates(&xyz_list);
483  element->FindParam(&dt,TimesteppingTimeStepEnum);
484  IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum);
485  IssmDouble g = element->FindParam(ConstantsGEnum);
486  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
487  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
488  Input2* geothermalflux_input = element->GetInput2(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
489  Input2* head_input = element->GetInput2(HydrologyHeadEnum); _assert_(head_input);
490  Input2* gap_input = element->GetInput2(HydrologyGapHeightEnum); _assert_(gap_input);
491  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
492  Input2* base_input = element->GetInput2(BaseEnum); _assert_(base_input);
493  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
494  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
495  Input2* englacial_input = element->GetInput2(HydrologyEnglacialInputEnum); _assert_(englacial_input);
496  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
497  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
498  Input2* lr_input = element->GetInput2(HydrologyBumpSpacingEnum); _assert_(lr_input);
499  Input2* br_input = element->GetInput2(HydrologyBumpHeightEnum); _assert_(br_input);
500 
501  /*Get conductivity from inputs*/
502  IssmDouble conductivity = GetConductivity(element);
503 
504  /*Build friction element, needed later: */
505  Friction* friction=new Friction(element,2);
506 
507  /*Keep track of weights*/
508  IssmDouble totalweights=0.;
509 
510  /* Start looping on the number of gaussian points: */
511  Gauss* gauss=element->NewGauss(2);
512  for(int ig=gauss->begin();ig<gauss->end();ig++){
513  gauss->GaussPoint(ig);
514 
515  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
516 
517  geothermalflux_input->GetInputValue(&G,gauss);
518  base_input->GetInputValue(&bed,gauss);
519  base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
520  thickness_input->GetInputValue(&thickness,gauss);
521  gap_input->GetInputValue(&gap,gauss);
522  head_input->GetInputValue(&head,gauss);
523  head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
524  englacial_input->GetInputValue(&ieb,gauss);
525  lr_input->GetInputValue(&lr,gauss);
526  br_input->GetInputValue(&br,gauss);
527  vx_input->GetInputValue(&vx,gauss);
528  vy_input->GetInputValue(&vy,gauss);
529 
530  /*Get ice A parameter*/
531  B_input->GetInputValue(&B,gauss);
532  n_input->GetInputValue(&n,gauss);
533  A=pow(B,-n);
534 
535  /*Compute beta term*/
536  if(gap<br)
537  beta = (br-gap)/lr;
538  else
539  beta = 0.;
540 
541  /*Compute frictional heat flux*/
542  friction->GetAlpha2(&alpha2,gauss);
543  vx_input->GetInputValue(&vx,gauss);
544  vy_input->GetInputValue(&vy,gauss);
545  frictionheat=alpha2*(vx*vx+vy*vy);
546 
547  /*Get water and ice pressures*/
548  IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.);
549  IssmDouble pressure_water = rho_water*g*(head-bed);
550  if(pressure_water>pressure_ice) pressure_water = pressure_ice;
551 
552  /* Compute change in sensible heat due to changes in pressure melting point*/
553  dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
554  dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
555  PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
556 
557  meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
558  _assert_(meltrate>0.);
559 
560  newgap += gauss->weight*Jdet*(gap+dt*(
561  meltrate/rho_ice
562  -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
563  +beta*sqrt(vx*vx+vy*vy)
564  ));
565  totalweights +=gauss->weight*Jdet;
566 
567  /* Compute basal water flux */
568  q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]));
569 
570  /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */
571  channelization += gauss->weight*Jdet*(meltrate/rho_ice/(meltrate/rho_ice+beta*sqrt(vx*vx+vy*vy)));
572  }
573 
574  /*Divide by connectivity*/
575  newgap = newgap/totalweights;
576  IssmDouble mingap = 1e-6;
577  if(newgap<mingap) newgap=mingap;
578 
579  /*Limit gap height to grow to surface*/
580  if(newgap>thickness)
581  newgap = thickness;
582 
583  /*Add new gap as an input*/
584  element->AddInput2(HydrologyGapHeightEnum,&newgap,P0Enum);
585 
586  /*Divide by connectivity, add basal flux as an input*/
587  q = q/totalweights;
588  element->AddInput2(HydrologyBasalFluxEnum,&q,P0Enum);
589 
590  /* Divide by connectivity, add degree of channelization as an input */
591  channelization = channelization/totalweights;
592  element->AddInput2(DegreeOfChannelizationEnum,&channelization,P0Enum);
593 
594  /*Clean up and return*/
595  xDelete<IssmDouble>(xyz_list);
596  delete friction;
597  delete gauss;
598 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Element::GetElementType
virtual int GetElementType(void)=0
HydrologyShaktiAnalysis::UpdateGapHeight
void UpdateGapHeight(FemModel *femmodel)
Definition: HydrologyShaktiAnalysis.cpp:457
HydrologyShaktiAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: HydrologyShaktiAnalysis.cpp:14
Element::GetInputListOnNodes
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1106
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
FrictionPEnum
@ FrictionPEnum
Definition: EnumDefinitions.h:578
HydrologyBumpSpacingEnum
@ HydrologyBumpSpacingEnum
Definition: EnumDefinitions.h:601
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
HydrologyGapHeightEnum
@ HydrologyGapHeightEnum
Definition: EnumDefinitions.h:614
HydrologyRequestedOutputsEnum
@ HydrologyRequestedOutputsEnum
Definition: EnumDefinitions.h:173
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
HydrologyShaktiAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: HydrologyShaktiAnalysis.cpp:220
Moulin
Definition: Moulin.h:21
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
HydrologyShaktiAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: HydrologyShaktiAnalysis.cpp:345
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
P0Enum
@ P0Enum
Definition: EnumDefinitions.h:661
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
HydrologyShaktiAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: HydrologyShaktiAnalysis.cpp:79
Elements
Declaration of Elements class.
Definition: Elements.h:17
HydrologyShaktiAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: HydrologyShaktiAnalysis.cpp:25
HydrologyHeadOldEnum
@ HydrologyHeadOldEnum
Definition: EnumDefinitions.h:616
MaterialsLatentheatEnum
@ MaterialsLatentheatEnum
Definition: EnumDefinitions.h:254
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
HydrologyRelaxationEnum
@ HydrologyRelaxationEnum
Definition: EnumDefinitions.h:172
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
OMEGA
#define OMEGA
Definition: HydrologyShaktiAnalysis.cpp:8
HydrologyShaktiAnalysisEnum
@ HydrologyShaktiAnalysisEnum
Definition: EnumDefinitions.h:1103
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
HydrologyShaktiAnalysis::Core
void Core(FemModel *femmodel)
Definition: HydrologyShaktiAnalysis.cpp:161
FrictionQEnum
@ FrictionQEnum
Definition: EnumDefinitions.h:580
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
Element::IsFloating
bool IsFloating()
Definition: Element.cpp:1987
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
HydrologyNeumannfluxEnum
@ HydrologyNeumannfluxEnum
Definition: EnumDefinitions.h:618
HydrologyShaktiAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: HydrologyShaktiAnalysis.cpp:164
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
HydrologyMoulinInputEnum
@ HydrologyMoulinInputEnum
Definition: EnumDefinitions.h:617
HydrologyShaktiAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: HydrologyShaktiAnalysis.cpp:82
HydrologyShaktiAnalysis.h
: header file for generic external result object
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
HydrologyShaktiAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: HydrologyShaktiAnalysis.cpp:342
FrictionLawEnum
@ FrictionLawEnum
Definition: EnumDefinitions.h:148
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
Neumannflux
Definition: Neumannflux.h:18
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
Element::NewGauss
virtual Gauss * NewGauss(void)=0
HydrologyEnglacialInputEnum
@ HydrologyEnglacialInputEnum
Definition: EnumDefinitions.h:613
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
HydrologyBumpHeightEnum
@ HydrologyBumpHeightEnum
Definition: EnumDefinitions.h:600
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
HydrologyShaktiAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: HydrologyShaktiAnalysis.cpp:168
StringArrayParam
Definition: StringArrayParam.h:20
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
HydrologyShaktiAnalysis::GetConductivity
IssmDouble GetConductivity(Element *element)
Definition: HydrologyShaktiAnalysis.cpp:436
NU
#define NU
Definition: HydrologyShaktiAnalysis.cpp:9
IntParam
Definition: IntParam.h:20
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
HydrologyNumRequestedOutputsEnum
@ HydrologyNumRequestedOutputsEnum
Definition: EnumDefinitions.h:170
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
Friction
Definition: Friction.h:13
HydrologyShaktiAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: HydrologyShaktiAnalysis.cpp:430
CreateSingleNodeToElementConnectivity
void CreateSingleNodeToElementConnectivity(IoModel *iomodel)
Definition: CreateSingleNodeToElementConnectivity.cpp:10
FemModel::elements
Elements * elements
Definition: FemModel.h:44
DegreeOfChannelizationEnum
@ DegreeOfChannelizationEnum
Definition: EnumDefinitions.h:521
Input2
Definition: Input2.h:18
Element::GetDofListLocal
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
Definition: Element.cpp:984
FemModel
Definition: FemModel.h:31
HydrologyshaktiEnum
@ HydrologyshaktiEnum
Definition: EnumDefinitions.h:1108
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
Loads
Declaration of Loads class.
Definition: Loads.h:16
HydrologyModelEnum
@ HydrologyModelEnum
Definition: EnumDefinitions.h:169
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
HydrologyReynoldsEnum
@ HydrologyReynoldsEnum
Definition: EnumDefinitions.h:619
HydrologyStorageEnum
@ HydrologyStorageEnum
Definition: EnumDefinitions.h:176
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
Friction::GetAlpha2
void GetAlpha2(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:164
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
ElementVector
Definition: ElementVector.h:20
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
BasalforcingsGeothermalfluxEnum
@ BasalforcingsGeothermalfluxEnum
Definition: EnumDefinitions.h:477
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
EffectivePressureEnum
@ EffectivePressureEnum
Definition: EnumDefinitions.h:548
IoModel
Definition: IoModel.h:48
HydrologyHeadEnum
@ HydrologyHeadEnum
Definition: EnumDefinitions.h:615
HydrologyShaktiAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: HydrologyShaktiAnalysis.cpp:171
CW
#define CW
Definition: HydrologyShaktiAnalysis.cpp:11
CT
#define CT
Definition: HydrologyShaktiAnalysis.cpp:10
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
HydrologyShaktiAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: HydrologyShaktiAnalysis.cpp:66
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Input2::GetInputAverage
virtual void GetInputAverage(IssmDouble *pvalue)
Definition: Input2.h:33
HydrologyShaktiAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: HydrologyShaktiAnalysis.cpp:137
Gauss
Definition: Gauss.h:8
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
HydrologyBasalFluxEnum
@ HydrologyBasalFluxEnum
Definition: EnumDefinitions.h:599
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
HydrologyShaktiAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: HydrologyShaktiAnalysis.cpp:348
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:15 for Ice Sheet System Model by   doxygen 1.8.19