Ice Sheet System Model  4.18
Code documentation
HydrologyGlaDSAnalysis.cpp
Go to the documentation of this file.
1 #include "./HydrologyGlaDSAnalysis.h"
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 
7 #define AEPS 2.2204460492503131E-015
8 
9 /*Model processing*/
10 void HydrologyGlaDSAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
11 
12  /*retrieve some parameters: */
13  int hydrology_model;
14  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
15 
16  if(hydrology_model!=HydrologyGlaDSEnum) return;
17 
18  IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spcphi",HydrologyGlaDSAnalysisEnum,P1Enum);
19 
20 }/*}}}*/
21 void HydrologyGlaDSAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
22 
23  /*Fetch parameters: */
24  int hydrology_model;
25  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
26 
27  /*Now, do we really want GlaDS?*/
28  if(hydrology_model!=HydrologyGlaDSEnum) return;
29 
30  /*Add channels?*/
31  int K,L;
32  bool ischannels;
33  IssmDouble* channelarea;
34  iomodel->FindConstant(&ischannels,"md.hydrology.ischannels");
35  iomodel->FetchData(&channelarea,&K,&L,"md.initialization.channelarea");
36  if(ischannels){
37  /*Get faces (edges in 2d)*/
38  CreateFaces(iomodel);
39  for(int i=0;i<iomodel->numberoffaces;i++){
40 
41  /*Get left and right elements*/
42  int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
43 
44  /*Now, if this element is not in the partition*/
45  if(!iomodel->my_elements[element]) continue;
46 
47  /*Add channelarea from initialization if exists*/
48  if(K!=0 && K!=iomodel->numberoffaces){
49  _error_("Unknown dimension for channel area initialization.");
50  }
51  if(K==0){
52  loads->AddObject(new Channel(i+1,0.,i,iomodel));
53  }
54  else{
55  loads->AddObject(new Channel(i+1,channelarea[i],i,iomodel));
56  }
57  iomodel->DeleteData(1,"md.initialization.channelarea");
58  }
59  }
60 
61  /*Create discrete loads for Moulins*/
62  CreateSingleNodeToElementConnectivity(iomodel);
63  for(int i=0;i<iomodel->numberofvertices;i++){
64  if (iomodel->domaintype!=Domain3DEnum){
65  /*keep only this partition's nodes:*/
66  if(iomodel->my_vertices[i]){
67  loads->AddObject(new Moulin(i+1,i,iomodel));
68  }
69  }
70  else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
71  if(iomodel->my_vertices[i]){
72  loads->AddObject(new Moulin(i+1,i,iomodel));
73  }
74  }
75  }
76  iomodel->DeleteData(1,"md.mesh.vertexonbase");
77 
78  /*Deal with Neumann BC*/
79  int M,N;
80  int *segments = NULL;
81  iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
82 
83  /*Check that the size seem right*/
84  _assert_(N==3); _assert_(M>=3);
85  for(int i=0;i<M;i++){
86  if(iomodel->my_elements[segments[i*3+2]-1]){
87  loads->AddObject(new Neumannflux(i+1,i,iomodel,segments));
88  }
89  }
90  xDelete<int>(segments);
91 
92 }/*}}}*/
93 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
94 
95  /*Fetch parameters: */
96  int hydrology_model;
97  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
98 
99  /*Now, do we really want GlaDS?*/
100  if(hydrology_model!=HydrologyGlaDSEnum) return;
101 
102  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
103  ::CreateNodes(nodes,iomodel,HydrologyGlaDSAnalysisEnum,P1Enum);
104  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
105 }/*}}}*/
106 int HydrologyGlaDSAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
107  return 1;
108 }/*}}}*/
109 void HydrologyGlaDSAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
110 
111  /*Fetch data needed: */
112  int hydrology_model,frictionlaw;
113  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
114 
115  /*Now, do we really want GlaDS?*/
116  if(hydrology_model!=HydrologyGlaDSEnum) return;
117 
118  /*Update elements: */
119  int counter=0;
120  for(int i=0;i<iomodel->numberofelements;i++){
121  if(iomodel->my_elements[i]){
122  Element* element=(Element*)elements->GetObjectByOffset(counter);
123  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
124  counter++;
125  }
126  }
127 
128  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
129  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
130  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
131  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
132  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
133  if(iomodel->domaintype!=Domain2DhorizontalEnum){
134  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
135  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
136  }
137  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
138  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
139  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
140  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
141  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
142  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
143  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
144  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);
145  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
146  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
147  iomodel->FindConstant(&frictionlaw,"md.friction.law");
148 
149  /*Friction law variables*/
150  switch(frictionlaw){
151  case 1:
152  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
153  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
154  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
155  break;
156  case 8:
157  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
158  break;
159  default:
160  _error_("Friction law "<< frictionlaw <<" not supported");
161  }
162 }/*}}}*/
163 void HydrologyGlaDSAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
164 
165  /*retrieve some parameters: */
166  int hydrology_model;
167  int numoutputs;
168  char** requestedoutputs = NULL;
169  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
170 
171  /*Now, do we really want GlaDS?*/
172  if(hydrology_model!=HydrologyGlaDSEnum) return;
173 
174  parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
175  parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
176  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.pressure_melt_coefficient",HydrologyPressureMeltCoefficientEnum));
177  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.cavity_spacing",HydrologyCavitySpacingEnum));
178  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.ischannels",HydrologyIschannelsEnum));
179  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.melt_flag",HydrologyMeltFlagEnum));
180  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum));
181  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_sheet_width",HydrologyChannelSheetWidthEnum));
182  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
183 
184  /*Deal with friction parameters*/
185  int frictionlaw;
186  iomodel->FindConstant(&frictionlaw,"md.friction.law");
187  if(frictionlaw==6){
188  parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
189  }
190  if(frictionlaw==4){
191  parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
192  parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
193  parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
194  }
195  if(frictionlaw==1 || frictionlaw==3 || frictionlaw==7){
196  parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
197  parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
198  }
199  if(frictionlaw==9){
200  parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
201  parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
202  parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
203  }
204 
205  /*Requested outputs*/
206  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
207  parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
208  if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
209  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
210 }/*}}}*/
211 
212 /*Finite Element Analysis*/
213 void HydrologyGlaDSAnalysis::Core(FemModel* femmodel){/*{{{*/
214  _error_("not implemented");
215 }/*}}}*/
216 ElementVector* HydrologyGlaDSAnalysis::CreateDVector(Element* element){/*{{{*/
217  /*Default, return NULL*/
218  return NULL;
219 }/*}}}*/
220 ElementMatrix* HydrologyGlaDSAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
221  _error_("Not implemented");
222 }/*}}}*/
223 ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/
224 
225  /*Intermediaries */
226  IssmDouble Jdet,dphi[3],h,k;
227  IssmDouble A,B,n,phi_old,phi,phi_0,H,b,v1;
228  IssmDouble* xyz_list = NULL;
229 
230  /*Hard coded coefficients*/
231  const IssmPDouble alpha = 5./4.;
232  const IssmPDouble beta = 3./2.;
233 
234  /*Fetch number of nodes and dof for this finite element*/
235  int numnodes = element->GetNumberOfNodes();
236 
237  /*Initialize Element vector and other vectors*/
238  ElementMatrix* Ke = element->NewElementMatrix();
239  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
240  IssmDouble* basis = xNew<IssmDouble>(numnodes);
241 
242  /*Retrieve all inputs and parameters*/
243  element->GetVerticesCoordinates(&xyz_list);
244 
245  /*Get all inputs and parameters*/
246  IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
247  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
248  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
249  IssmDouble g = element->FindParam(ConstantsGEnum);
250  IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum);
251  Input2* k_input = element->GetInput2(HydrologySheetConductivityEnum);_assert_(k_input);
252  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
253  Input2* h_input = element->GetInput2(HydrologySheetThicknessEnum); _assert_(h_input);
254  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
255  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
256  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
257  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
258 
259  /* Start looping on the number of gaussian points: */
260  Gauss* gauss=element->NewGauss(2);
261  for(int ig=gauss->begin();ig<gauss->end();ig++){
262  gauss->GaussPoint(ig);
263 
264  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
265  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
266  element->NodalFunctions(basis,gauss);
267 
268  phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
269  phi_input->GetInputValue(&phi,gauss);
270  h_input->GetInputValue(&h,gauss);
271  k_input->GetInputValue(&k,gauss);
272  B_input->GetInputValue(&B,gauss);
273  n_input->GetInputValue(&n,gauss);
274  b_input->GetInputValue(&b,gauss);
275  H_input->GetInputValue(&H,gauss);
276 
277  /*Get norm of gradient of hydraulic potential and make sure it is >0*/
278  IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
279  if(normgradphi < AEPS) normgradphi = AEPS;
280 
281  IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);
282 
283  /*Diffusive term*/
284  for(int i=0;i<numnodes;i++){
285  for(int j=0;j<numnodes;j++){
286  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
287  coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
288  + coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
289  }
290  }
291 
292  /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
293  phi_0 = rho_water*g*b + rho_ice*g*H;
294  A=pow(B,-n);
295  v1 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*( - n));
296  for(int i=0;i<numnodes;i++){
297  for(int j=0;j<numnodes;j++){
298  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
299  }
300  }
301 
302  /*Transient term if dt>0*/
303  if(dt>0.){
304  /*Diffusive term*/
305  for(int i=0;i<numnodes;i++){
306  for(int j=0;j<numnodes;j++){
307  Ke->values[i*numnodes+j] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*basis[i]*basis[j];
308  }
309  }
310  }
311 
312  }
313 
314  /*Clean up and return*/
315  xDelete<IssmDouble>(xyz_list);
316  xDelete<IssmDouble>(basis);
317  xDelete<IssmDouble>(dbasis);
318  delete gauss;
319  return Ke;
320 }/*}}}*/
321 ElementVector* HydrologyGlaDSAnalysis::CreatePVector(Element* element){/*{{{*/
322 
323  /*Skip if water or ice shelf element*/
324  if(element->IsFloating()) return NULL;
325 
326  /*Intermediaries */
327  bool meltflag;
328  IssmDouble Jdet,w,v2,vx,vy,ub,h,h_r;
329  IssmDouble G,m,frictionheat,alpha2;
330  IssmDouble A,B,n,phi_old,phi,phi_0;
331  IssmDouble H,b;
332  IssmDouble* xyz_list = NULL;
333 
334  /*Fetch number of nodes and dof for this finite element*/
335  int numnodes = element->GetNumberOfNodes();
336 
337  /*Initialize Element vector and other vectors*/
338  ElementVector* pe = element->NewElementVector();
339  IssmDouble* basis = xNew<IssmDouble>(numnodes);
340 
341  /*Retrieve all inputs and parameters*/
342  element->GetVerticesCoordinates(&xyz_list);
343  element->FindParam(&meltflag,HydrologyMeltFlagEnum);
344  IssmDouble L = element->FindParam(MaterialsLatentheatEnum);
345  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
346  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
347  IssmDouble l_r = element->FindParam(HydrologyCavitySpacingEnum);
348  IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
349  IssmDouble g = element->FindParam(ConstantsGEnum);
350  IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum);
351  Input2* hr_input = element->GetInput2(HydrologyBumpHeightEnum);_assert_(hr_input);
352  Input2* vx_input = element->GetInput2(VxEnum);_assert_(vx_input);
353  Input2* vy_input = element->GetInput2(VyEnum);_assert_(vy_input);
354  Input2* h_input = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
355  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
356  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
357  Input2* G_input = element->GetInput2(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
358  Input2* m_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);_assert_(m_input);
359  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
360  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
361  Input2* phiold_input = element->GetInput2(HydraulicPotentialOldEnum); _assert_(phiold_input);
362  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
363 
364  /*Build friction element, needed later: */
365  Friction* friction=new Friction(element,2);
366 
367  /* Start looping on the number of gaussian points: */
368  Gauss* gauss=element->NewGauss(2);
369  for(int ig=gauss->begin();ig<gauss->end();ig++){
370  gauss->GaussPoint(ig);
371 
372  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
373  element->NodalFunctions(basis,gauss);
374 
375  /*Get input values at gauss points*/
376  vx_input->GetInputValue(&vx,gauss);
377  vy_input->GetInputValue(&vy,gauss);
378  h_input->GetInputValue(&h,gauss);
379  G_input->GetInputValue(&G,gauss);
380  B_input->GetInputValue(&B,gauss);
381  n_input->GetInputValue(&n,gauss);
382  hr_input->GetInputValue(&h_r,gauss);
383  phi_input->GetInputValue(&phi,gauss);
384  b_input->GetInputValue(&b,gauss);
385  H_input->GetInputValue(&H,gauss);
386 
387  /*Get basal velocity*/
388  ub = sqrt(vx*vx + vy*vy);
389 
390  /*Compute cavity opening w*/
391  w = 0.;
392  if(h<h_r) w = ub*(h_r-h)/l_r;
393 
394  /*Compute frictional heat flux*/
395  friction->GetAlpha2(&alpha2,gauss);
396  frictionheat=alpha2*ub*ub;
397 
398  /*Compute melt (if necessary)*/
399  if(!meltflag){
400  m = (G + frictionheat)/(rho_ice*L);
401  }
402  else{
403  m_input->GetInputValue(&m,gauss);
404  }
405 
406  /*Compute closing rate*/
407  phi_0 = rho_water*g*b + rho_ice*g*H;
408  A=pow(B,-n);
409  v2 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
410 
411  for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
412 
413  /*Transient term if dt>0*/
414  if(dt>0.){
415  phiold_input->GetInputValue(&phi_old,gauss);
416  for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
417  }
418  }
419 
420  /*Clean up and return*/
421  xDelete<IssmDouble>(xyz_list);
422  xDelete<IssmDouble>(basis);
423  delete gauss;
424  delete friction;
425  return pe;
426 }/*}}}*/
427 void HydrologyGlaDSAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
428 
429  element->GetSolutionFromInputsOneDof(solution,HydraulicPotentialEnum);
430 
431 }/*}}}*/
432 void HydrologyGlaDSAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
433  _error_("Not implemented yet");
434 }/*}}}*/
435 void HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
436  element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
437 }/*}}}*/
438 void HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
439  /*Default, do nothing*/
440  return;
441 }/*}}}*/
442 
443 /*GlaDS specifics*/
444 void HydrologyGlaDSAnalysis::SetChannelCrossSectionOld(FemModel* femmodel){/*{{{*/
445 
446  bool ischannels;
447  femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
448  if(!ischannels) return;
449 
450  for(int i=0;i<femmodel->loads->Size();i++){
451  if(femmodel->loads->GetEnum(i)==ChannelEnum){
452  Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
453  channel->SetChannelCrossSectionOld();
454  }
455  }
456 
457 }/*}}}*/
458 void HydrologyGlaDSAnalysis::UpdateSheetThickness(FemModel* femmodel){/*{{{*/
459 
460  for(int j=0;j<femmodel->elements->Size();j++){
461  Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
462  UpdateSheetThickness(element);
463  }
464 
465 }/*}}}*/
466 void HydrologyGlaDSAnalysis::UpdateSheetThickness(Element* element){/*{{{*/
467 
468  /*Skip if water or ice shelf element*/
469  if(element->IsFloating()) return;
470 
471  /*Intermediaries */
472  IssmDouble Jdet,vx,vy,ub,h_old,N,h_r,H,b;
473  IssmDouble A,B,n,phi,phi_0;
474  IssmDouble alpha,beta;
475 
476  /*Fetch number vertices for this element*/
477  int numvertices = element->GetNumberOfVertices();
478 
479  /*Initialize new sheet thickness*/
480  IssmDouble* h_new = xNew<IssmDouble>(numvertices);
481 
482  /*Retrieve all inputs and parameters*/
483  IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
484  IssmDouble l_r = element->FindParam(HydrologyCavitySpacingEnum);
485  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
486  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
487  IssmDouble g = element->FindParam(ConstantsGEnum);
488  Input2* hr_input = element->GetInput2(HydrologyBumpHeightEnum);_assert_(hr_input);
489  Input2* vx_input = element->GetInput2(VxEnum);_assert_(vx_input);
490  Input2* vy_input = element->GetInput2(VyEnum);_assert_(vy_input);
491  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
492  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
493  Input2* hold_input = element->GetInput2(HydrologySheetThicknessOldEnum);_assert_(hold_input);
494  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
495  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
496  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
497 
498  /* Start looping on the number of gaussian points: */
499  Gauss* gauss=element->NewGauss();
500  for(int iv=0;iv<numvertices;iv++){
501  gauss->GaussVertex(iv);
502 
503  /*Get input values at gauss points*/
504  phi_input->GetInputValue(&phi,gauss);
505  vx_input->GetInputValue(&vx,gauss);
506  vy_input->GetInputValue(&vy,gauss);
507  hold_input->GetInputValue(&h_old,gauss);
508  B_input->GetInputValue(&B,gauss);
509  n_input->GetInputValue(&n,gauss);
510  hr_input->GetInputValue(&h_r,gauss);
511  b_input->GetInputValue(&b,gauss);
512  H_input->GetInputValue(&H,gauss);
513 
514  /*Get values for a few potentials*/
515  phi_0 = rho_water*g*b + rho_ice*g*H;
516  N = phi_0 - phi;
517 
518  /*Get basal velocity*/
519  ub = sqrt(vx*vx + vy*vy);
520 
521  /*Get A from B and n*/
522  A=pow(B,-n);
523 
524  /*Define alpha and beta*/
525  if(h_old<h_r){
526  alpha = -ub/l_r - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
527  beta = ub*h_r/l_r;
528  }
529  else{
530  alpha = - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
531  beta = 0.;
532  }
533 
534  /*Get new sheet thickness*/
535  h_new[iv] = ODE1(alpha,beta,h_old,dt,1);
536 
537  /*Make sure it is positive*/
538  if(h_new[iv]<AEPS) h_new[iv] = AEPS;
539  }
540 
541  element->AddInput2(HydrologySheetThicknessEnum,h_new,P1Enum);
542 
543  /*Clean up and return*/
544  xDelete<IssmDouble>(h_new);
545  delete gauss;
546 }/*}}}*/
547 void HydrologyGlaDSAnalysis::UpdateEffectivePressure(FemModel* femmodel){/*{{{*/
548 
549  for(int j=0;j<femmodel->elements->Size();j++){
550  Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
551  UpdateEffectivePressure(element);
552  }
553 
554 }/*}}}*/
555 void HydrologyGlaDSAnalysis::UpdateEffectivePressure(Element* element){/*{{{*/
556 
557  /*Skip if water or ice shelf element*/
558  if(element->IsFloating()) return;
559 
560  /*Intermediary*/
561  IssmDouble phi_0, phi_m, p_i;
562  IssmDouble H,b,phi;
563 
564  int numnodes = element->GetNumberOfNodes();
565 
566  /*Get thickness and base on nodes to apply cap on water head*/
567  IssmDouble* N = xNew<IssmDouble>(numnodes);
568  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
569  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
570  IssmDouble g = element->FindParam(ConstantsGEnum);
571  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
572  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
573  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
574 
575  /* Start looping on the number of gaussian points: */
576  Gauss* gauss=element->NewGauss();
577  for(int iv=0;iv<numnodes;iv++){
578  gauss->GaussNode(element->FiniteElement(),iv);
579 
580  /*Get input values at gauss points*/
581  H_input->GetInputValue(&H,gauss);
582  b_input->GetInputValue(&b,gauss);
583  phi_input->GetInputValue(&phi,gauss);
584 
585  /*Elevation potential*/
586  phi_m = rho_water*g*b;
587 
588  /*Compute overburden pressure*/
589  p_i = rho_ice*g*H;
590 
591  /*Copmute overburden potential*/
592  phi_0 = phi_m + p_i;
593 
594  /*Calculate effective pressure*/
595  N[iv] = phi_0 - phi;
596 
597  if(xIsNan<IssmDouble>(N[iv])) _error_("NaN found in solution vector");
598  if(xIsInf<IssmDouble>(N[iv])) _error_("Inf found in solution vector");
599  }
600 
601  element->AddInput2(EffectivePressureEnum,N,element->FiniteElement());
602 
603  /*Clean up and return*/
604  delete gauss;
605  xDelete<IssmDouble>(N);
606 }/*}}}*/
607 void HydrologyGlaDSAnalysis::UpdateChannelCrossSection(FemModel* femmodel){/*{{{*/
608 
609  bool ischannels;
610  femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
611  if(!ischannels) return;
612 
613  for(int i=0;i<femmodel->loads->Size();i++){
614  if(femmodel->loads->GetEnum(i)==ChannelEnum){
615  Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
616  channel->UpdateChannelCrossSection();
617  }
618  }
619 
620 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
HydrologyGlaDSAnalysis::Core
void Core(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:213
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Channel
Definition: Channel.h:18
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
HydrologyGlaDSAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: HydrologyGlaDSAnalysis.cpp:163
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
FrictionPEnum
@ FrictionPEnum
Definition: EnumDefinitions.h:578
HydrologySheetThicknessOldEnum
@ HydrologySheetThicknessOldEnum
Definition: EnumDefinitions.h:622
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
HydrologyRequestedOutputsEnum
@ HydrologyRequestedOutputsEnum
Definition: EnumDefinitions.h:173
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
HydrologyGlaDSAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: HydrologyGlaDSAnalysis.cpp:321
Moulin
Definition: Moulin.h:21
HydraulicPotentialEnum
@ HydraulicPotentialEnum
Definition: EnumDefinitions.h:597
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
FrictionGammaEnum
@ FrictionGammaEnum
Definition: EnumDefinitions.h:147
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
HydrologyGlaDSAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: HydrologyGlaDSAnalysis.cpp:109
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
HydrologyPressureMeltCoefficientEnum
@ HydrologyPressureMeltCoefficientEnum
Definition: EnumDefinitions.h:171
HydrologyGlaDSAnalysis::UpdateSheetThickness
void UpdateSheetThickness(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:458
MaterialsLatentheatEnum
@ MaterialsLatentheatEnum
Definition: EnumDefinitions.h:254
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
IoModel::numberoffaces
int numberoffaces
Definition: IoModel.h:97
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
FrictionQEnum
@ FrictionQEnum
Definition: EnumDefinitions.h:580
Gauss::GaussNode
virtual void GaussNode(int finitelement, int iv)=0
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
Element::IsFloating
bool IsFloating()
Definition: Element.cpp:1987
FrictionCouplingEnum
@ FrictionCouplingEnum
Definition: EnumDefinitions.h:143
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
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
HydrologyGlaDSAnalysisEnum
@ HydrologyGlaDSAnalysisEnum
Definition: EnumDefinitions.h:1100
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
ChannelEnum
@ ChannelEnum
Definition: EnumDefinitions.h:1008
HydrologyGlaDSAnalysis::UpdateEffectivePressure
void UpdateEffectivePressure(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:547
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
HydrologySheetConductivityEnum
@ HydrologySheetConductivityEnum
Definition: EnumDefinitions.h:620
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
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
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
HydrologyGlaDSAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: HydrologyGlaDSAnalysis.cpp:10
Element::NewGauss
virtual Gauss * NewGauss(void)=0
HydrologySheetThicknessEnum
@ HydrologySheetThicknessEnum
Definition: EnumDefinitions.h:621
HydrologyGlaDSAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: HydrologyGlaDSAnalysis.cpp:216
Element::InputUpdateFromSolutionOneDof
virtual void InputUpdateFromSolutionOneDof(IssmDouble *solution, int inputenum)=0
HydrologyGlaDSAnalysis.h
: header file for generic external result object
HydrologyBumpHeightEnum
@ HydrologyBumpHeightEnum
Definition: EnumDefinitions.h:600
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
HydrologyGlaDSAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: HydrologyGlaDSAnalysis.cpp:427
StringArrayParam
Definition: StringArrayParam.h:20
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
HydrologyGlaDSAnalysis::UpdateChannelCrossSection
void UpdateChannelCrossSection(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:607
IntParam
Definition: IntParam.h:20
HydrologyCavitySpacingEnum
@ HydrologyCavitySpacingEnum
Definition: EnumDefinitions.h:163
Channel::SetChannelCrossSectionOld
void SetChannelCrossSectionOld(void)
Definition: Channel.cpp:596
DataSet::GetEnum
int GetEnum()
Definition: DataSet.cpp:324
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
FemModel::loads
Loads * loads
Definition: FemModel.h:54
HydrologyNumRequestedOutputsEnum
@ HydrologyNumRequestedOutputsEnum
Definition: EnumDefinitions.h:170
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
Friction
Definition: Friction.h:13
CreateSingleNodeToElementConnectivity
void CreateSingleNodeToElementConnectivity(IoModel *iomodel)
Definition: CreateSingleNodeToElementConnectivity.cpp:10
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
Gauss::GaussVertex
virtual void GaussVertex(int iv)=0
ODE1
IssmDouble ODE1(IssmDouble alpha, IssmDouble beta, IssmDouble Si, IssmDouble dt, int method)
Definition: ODE1.cpp:5
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
HydrologyChannelConductivityEnum
@ HydrologyChannelConductivityEnum
Definition: EnumDefinitions.h:164
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
IoModel::faces
int * faces
Definition: IoModel.h:88
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
HydrologyEnglacialVoidRatioEnum
@ HydrologyEnglacialVoidRatioEnum
Definition: EnumDefinitions.h:166
HydrologyChannelSheetWidthEnum
@ HydrologyChannelSheetWidthEnum
Definition: EnumDefinitions.h:165
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
FrictionEffectivePressureLimitEnum
@ FrictionEffectivePressureLimitEnum
Definition: EnumDefinitions.h:145
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
HydrologyGlaDSAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: HydrologyGlaDSAnalysis.cpp:106
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
Element::FiniteElement
virtual int FiniteElement(void)=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
HydrologyGlaDSAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: HydrologyGlaDSAnalysis.cpp:220
HydrologyGlaDSEnum
@ HydrologyGlaDSEnum
Definition: EnumDefinitions.h:1101
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
HydrologyGlaDSAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: HydrologyGlaDSAnalysis.cpp:21
BasalforcingsGeothermalfluxEnum
@ BasalforcingsGeothermalfluxEnum
Definition: EnumDefinitions.h:477
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
HydrologyMeltFlagEnum
@ HydrologyMeltFlagEnum
Definition: EnumDefinitions.h:168
EffectivePressureEnum
@ EffectivePressureEnum
Definition: EnumDefinitions.h:548
IoModel
Definition: IoModel.h:48
HydrologyGlaDSAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:438
CreateFaces
void CreateFaces(IoModel *iomodel)
Definition: CreateFaces.cpp:9
HydrologyGlaDSAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: HydrologyGlaDSAnalysis.cpp:223
HydrologyGlaDSAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: HydrologyGlaDSAnalysis.cpp:432
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
AEPS
#define AEPS
Definition: HydrologyGlaDSAnalysis.cpp:7
ElementMatrix
Definition: ElementMatrix.h:19
HydrologyGlaDSAnalysis::SetChannelCrossSectionOld
void SetChannelCrossSectionOld(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:444
Vector< IssmDouble >
HydrologyGlaDSAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: HydrologyGlaDSAnalysis.cpp:435
Channel::UpdateChannelCrossSection
void UpdateChannelCrossSection(void)
Definition: Channel.cpp:601
HydrologyIschannelsEnum
@ HydrologyIschannelsEnum
Definition: EnumDefinitions.h:167
HydrologyGlaDSAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: HydrologyGlaDSAnalysis.cpp:93
Gauss
Definition: Gauss.h:8
HydraulicPotentialOldEnum
@ HydraulicPotentialOldEnum
Definition: EnumDefinitions.h:598
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
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