Ice Sheet System Model  4.18
Code documentation
Channel.cpp
Go to the documentation of this file.
1 
5 /*Headers:*/
6 /*{{{*/
7 #ifdef HAVE_CONFIG_H
8 #include <config.h>
9 #else
10 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11 #endif
12 
13 #include "shared/shared.h"
14 #include "../classes.h"
15 /*}}}*/
16 
17 /*Macros*/
18 #define NUMNODES 2
19 #define NUMVERTICES 2
20 
21 #define C_W 4.22e3 /*specific heat capacity of water (J/kg/K)*/
22 #define ALPHA_C 5./4.
23 #define BETA_C 3./2.
24 /*Make sure these are the same as in HydrologyGlaDSAnalysis::CreateKMatrix*/
25 #define ALPHA_S 5./4.
26 #define BETA_S 3./2.
27 #define AEPS 2.2204460492503131E-015
28 
29 /*Channel constructors and destructor*/
30 Channel::Channel(){/*{{{*/
31  this->id = -1;
32  this->sid = -1;
33  this->parameters = NULL;
34  this->helement = NULL;
35  this->element = NULL;
36  this->hnodes = NULL;
37  this->hvertices = NULL;
38  this->nodes = NULL;
39 }
40 /*}}}*/
41 Channel::Channel(int channel_id,IssmDouble channelarea,int index,IoModel* iomodel){/*{{{*/
42 //Channel::Channel(int channel_id,int i,int index,IoModel* iomodel)
43 
44  this->id=channel_id;
45  this->sid=channel_id-1;
46  this->parameters = NULL;
47  this->element = NULL;
48  this->nodes = NULL;
49 
50  /*Set channel cross section to 0*/
51  //this->S = 0.;
52  //this->Sold = 0.;
53  this->S = channelarea;
54  this->Sold = channelarea;
55 
56  /*Get edge info*/
57  int i1 = iomodel->faces[4*index+0];
58  int i2 = iomodel->faces[4*index+1];
59  int e1 = iomodel->faces[4*index+2];
60  int e2 = iomodel->faces[4*index+3];
61 
62  if(e2==-1){
63  this->boundary = true;
64  }
65  else{
66  this->boundary = false;
67  }
68 
69  /*Set Element hook (4th column may be -1 for boundary edges)*/
70  this->helement = new Hook(&e1,1);
71 
72  /*Set Vertices hooks (4th column may be -1 for boundary edges)*/
73  int channel_vertex_ids[2];
74  channel_vertex_ids[0]=i1;
75  channel_vertex_ids[1]=i2;
76  this->hvertices =new Hook(&channel_vertex_ids[0],2);
77 
78  /*Set Nodes hooks (!! Assumes P1 CG)*/
79  int channel_node_ids[2];
80  channel_node_ids[0]=i1;
81  channel_node_ids[1]=i2;
82  this->hnodes=new Hook(&channel_node_ids[0],2);
83 
84 }/*}}}*/
85 Channel::~Channel(){/*{{{*/
86  this->parameters=NULL;
87  delete helement;
88  delete hnodes;
89  delete hvertices;
90 }/*}}}*/
91 
92 /*Object virtual functions definitions:*/
93 Object* Channel::copy() {/*{{{*/
94 
95  Channel* channel=NULL;
96 
97  channel=new Channel();
98 
99  /*copy fields: */
100  channel->id=this->id;
101  channel->S=this->S;
102 
103  /*point parameters: */
104  channel->parameters=this->parameters;
105 
106  /*now deal with hooks and objects: */
107  channel->hnodes = (Hook*)this->hnodes->copy();
108  channel->hvertices = (Hook*)this->hvertices->copy();
109  channel->helement = (Hook*)this->helement->copy();
110 
111  /*corresponding fields*/
112  channel->nodes = (Node**)channel->hnodes->deliverp();
113  channel->vertices = (Vertex**)channel->hvertices->deliverp();
114  channel->element = (Element*)channel->helement->delivers();
115 
116  return channel;
117 }
118 /*}}}*/
119 void Channel::DeepEcho(void){/*{{{*/
120 
121  _printf_("Channel:\n");
122  _printf_(" id: " << id << "\n");
123  _printf_(" S: " << S << "\n");
124  hnodes->DeepEcho();
125  hvertices->DeepEcho();
126  helement->DeepEcho();
127  _printf_(" parameters\n");
128  if(parameters)
129  parameters->DeepEcho();
130  else
131  _printf_(" NULL\n");
132 }
133 /*}}}*/
134 void Channel::Echo(void){/*{{{*/
135  _printf_("Channel:\n");
136  _printf_(" id: " << id << "\n");
137  _printf_(" S: " << S << "\n");
138  hnodes->Echo();
139  hvertices->Echo();
140  helement->Echo();
141  _printf_(" parameters: " << parameters << "\n");
142 }
143 /*}}}*/
144 int Channel::Id(void){/*{{{*/
145  return id;
146 }
147 /*}}}*/
148 void Channel::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
149 
150  _assert_(this);
151 
152  /*ok, marshall operations: */
153  MARSHALLING_ENUM(ChannelEnum);
154  MARSHALLING(id);
155  MARSHALLING(S);
156 
157  if(marshall_direction==MARSHALLING_BACKWARD){
158  this->hnodes = new Hook();
159  this->hvertices = new Hook();
160  this->helement = new Hook();
161  }
162 
163  this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
164  this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
165  this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
166 
167  /*corresponding fields*/
168  nodes =(Node**)this->hnodes->deliverp();
169  vertices =(Vertex**)this->hvertices->deliverp();
170  element =(Element*)this->helement->delivers();
171 
172 }
173 /*}}}*/
174 int Channel::ObjectEnum(void){/*{{{*/
175  return ChannelEnum;
176 }/*}}}*/
177 
178 /*Load virtual functions definitions:*/
179 void Channel::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
180 
181  /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
182  * datasets, using internal ids and offsets hidden in hooks: */
183  hnodes->configure((DataSet*)nodesin);
184  hvertices->configure((DataSet*)verticesin);
185  helement->configure((DataSet*)elementsin);
186 
187  /*Initialize hooked fields*/
188  this->nodes = (Node**)hnodes->deliverp();
189  this->vertices = (Vertex**)hvertices->deliverp();
190  this->element = (Element*)helement->delivers();
191 
192  /*point parameters to real dataset: */
193  this->parameters=parametersin;
194 }
195 /*}}}*/
196 void Channel::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
197 
198  /*recover some parameters*/
199  ElementMatrix* Ke=NULL;
200  int analysis_type;
201  this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
202 
203  switch(analysis_type){
204  case HydrologyGlaDSAnalysisEnum:
205  Ke = this->CreateKMatrixHydrologyGlaDS();
206  break;
207  default:
208  _error_("Don't know why we should be here");
209  }
210 
211  /*Add to global matrix*/
212  if(Ke){
213  Ke->AddToGlobal(Kff,Kfs);
214  delete Ke;
215  }
216 
217 }
218 /*}}}*/
219 void Channel::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
220 
221  /*recover some parameters*/
222  ElementVector* pe=NULL;
223  int analysis_type;
224  this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
225 
226  switch(analysis_type){
227  case HydrologyGlaDSAnalysisEnum:
228  pe = this->CreatePVectorHydrologyGlaDS();
229  break;
230  default:
231  _error_("Don't know why we should be here");
232  }
233 
234  /*Add to global matrix*/
235  if(pe){
236  pe->AddToGlobal(pf);
237  delete pe;
238  }
239 
240 }
241 /*}}}*/
242 void Channel::GetNodesLidList(int* lidlist){/*{{{*/
243 
244  _assert_(lidlist);
245  _assert_(nodes);
246 
247  for(int i=0;i<NUMNODES;i++) lidlist[i]=nodes[i]->Lid();
248 }
249 /*}}}*/
250 void Channel::GetNodesSidList(int* sidlist){/*{{{*/
251 
252  _assert_(sidlist);
253  _assert_(nodes);
254 
255  for(int i=0;i<NUMNODES;i++) sidlist[i]=nodes[i]->Sid();
256 }
257 /*}}}*/
258 int Channel::GetNumberOfNodes(void){/*{{{*/
259  return NUMNODES;
260 }
261 /*}}}*/
262 bool Channel::IsPenalty(void){/*{{{*/
263  return false;
264 }
265 /*}}}*/
266 void Channel::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
267 
268  /*No stiffness loads applied, do nothing: */
269  return;
270 
271 }
272 /*}}}*/
273 void Channel::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
274 
275  /*No penalty loads applied, do nothing: */
276  return;
277 
278 }
279 /*}}}*/
280 void Channel::ResetHooks(){/*{{{*/
281 
282  this->nodes=NULL;
283  this->vertices=NULL;
284  this->element=NULL;
285  this->parameters=NULL;
286 
287  /*Get Element type*/
288  this->hnodes->reset();
289  this->hvertices->reset();
290  this->helement->reset();
291 
292 }
293 /*}}}*/
294 void Channel::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
295 
296 }
297 /*}}}*/
298 void Channel::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
299 
300  /*Output */
301  int d_nz = 0;
302  int o_nz = 0;
303 
304  /*Loop over all nodes*/
305  for(int i=0;i<this->GetNumberOfNodes();i++){
306 
307  if(!flags[this->nodes[i]->Lid()]){
308 
309  /*flag current node so that no other element processes it*/
310  flags[this->nodes[i]->Lid()]=true;
311 
312  int counter=0;
313  while(flagsindices[counter]>=0) counter++;
314  flagsindices[counter]=this->nodes[i]->Lid();
315 
316  /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
317  switch(set2_enum){
318  case FsetEnum:
319  if(nodes[i]->fsize){
320  if(this->nodes[i]->IsClone())
321  o_nz += 1;
322  else
323  d_nz += 1;
324  }
325  break;
326  case GsetEnum:
327  if(nodes[i]->gsize){
328  if(this->nodes[i]->IsClone())
329  o_nz += 1;
330  else
331  d_nz += 1;
332  }
333  break;
334  case SsetEnum:
335  if(nodes[i]->ssize){
336  if(this->nodes[i]->IsClone())
337  o_nz += 1;
338  else
339  d_nz += 1;
340  }
341  break;
342  default: _error_("not supported");
343  }
344  }
345  }
346 
347  /*Assign output pointers: */
348  *pd_nz=d_nz;
349  *po_nz=o_nz;
350 }
351 /*}}}*/
352 
353 /*Channel specific functions*/
354 ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
355 
356  /*Initialize Element matrix and return if necessary*/
357  Tria* tria=(Tria*)element;
358  if(!tria->IsIceInElement()) return NULL;
359  _assert_(tria->FiniteElement()==P1Enum);
360  int index1=tria->GetVertexIndex(vertices[0]);
361  int index2=tria->GetVertexIndex(vertices[1]);
362 
363  /*Intermediaries */
364  IssmDouble Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
365  IssmDouble A,B,n,phi_old,phi,phi_0,dPw,ks,Ngrad;
366  IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
367  IssmDouble xyz_list[NUMVERTICES][3];
368  IssmDouble xyz_list_tria[3][3];
369  const int numnodes = NUMNODES;
370 
371  /*Initialize Element vector and other vectors*/
372  ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters);
373  IssmDouble basis[NUMNODES];
374  IssmDouble dbasisdx[2*NUMNODES];
375  IssmDouble dbasisds[NUMNODES];
376 
377  /*Retrieve all inputs and parameters*/
378  GetVerticesCoordinates(&xyz_list[0][0] ,this->vertices,NUMVERTICES);
379  GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
380 
381  IssmDouble L = element->FindParam(MaterialsLatentheatEnum);
382  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
383  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
384  IssmDouble g = element->FindParam(ConstantsGEnum);
385  IssmDouble kc = element->FindParam(HydrologyChannelConductivityEnum);
386  IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum);
387  IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum);
388 
389  Input2* h_input = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
390  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
391  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
392  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
393  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
394  Input2* ks_input = element->GetInput2(HydrologySheetConductivityEnum); _assert_(ks_input);
395  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
396 
397  /*Get tangent vector*/
398  IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
399  IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
400  IssmDouble Lt = sqrt(tx*tx+ty*ty);
401  tx = tx/Lt;
402  ty = ty/Lt;
403 
404  /* Start looping on the number of gaussian points: */
405  Gauss* gauss=new GaussTria(index1,index2,2);
406  for(int ig=gauss->begin();ig<gauss->end();ig++){
407  gauss->GaussPoint(ig);
408 
409  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
410  tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
411  tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement());
412  dbasisds[0] = dbasisdx[0*2+0]*tx + dbasisdx[0*2+1]*ty;
413  dbasisds[1] = dbasisdx[1*2+0]*tx + dbasisdx[1*2+1]*ty;
414 
415  /*Get input values at gauss points*/
416  phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
417  b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
418  phi_input->GetInputValue(&phi,gauss);
419  h_input->GetInputValue(&h,gauss);
420  ks_input->GetInputValue(&ks,gauss);
421  B_input->GetInputValue(&B,gauss);
422  n_input->GetInputValue(&n,gauss);
423  b_input->GetInputValue(&b,gauss);
424  H_input->GetInputValue(&H,gauss);
425 
426  /*Get values for a few potentials*/
427  phi_0 = rho_water*g*b + rho_ice*g*H;
428  dphids = dphi[0]*tx + dphi[1]*ty;
429  dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
430  Ngrad = fabs(dphids);
431  if(Ngrad<AEPS) Ngrad = AEPS;
432 
433  /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
434  IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
435  IssmDouble Ks = ks * pow(h ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
436 
437  /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
438  qc = - Ks * dphids;
439 
440  /*d(phi - phi_m)/ds*/
441  dPw = dphids - dphimds;
442 
443  /*Compute f factor*/
444  fFactor = 0.;
445  if(this->S>0. || qc*dPw>0.){
446  fFactor = lc * qc;
447  }
448 
449  /*Compute Afactor and Bfactor*/
450  Afactor = C_W*c_t*rho_water;
451  Bfactor = 1./L * (1./rho_ice - 1./rho_water);
452  if(dphids>0){
453  Xifactor = + Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
454  }
455  else{
456  Xifactor = - Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
457  }
458 
459  /*Diffusive term*/
460  for(int i=0;i<numnodes;i++){
461  for(int j=0;j<numnodes;j++){
462  /*GlaDSCoupledSolver.F90 line 1659*/
463  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
464  +Kc*dbasisds[i]*dbasisds[j] /*Diffusion term*/
465  - Afactor * Bfactor* Kc * dPw * basis[i] * dbasisds[j] /*First part of Pi*/
466  + Afactor * fFactor * Bfactor * basis[i] * dbasisds[j] /*Second part of Pi*/
467  + Xifactor* basis[i] * dbasisds[j] /*Xi term*/
468  );
469  }
470  }
471 
472  /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
473  A=pow(B,-n);
474  v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n));
475  for(int i=0;i<numnodes;i++){
476  for(int j=0;j<numnodes;j++){
477  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
478  }
479  }
480  }
481 
482  /*Clean up and return*/
483  delete gauss;
484  return Ke;
485 }
486 /*}}}*/
487 ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/
488 
489  /*Initialize Element matrix and return if necessary*/
490  Tria* tria=(Tria*)element;
491  if(!tria->IsIceInElement()) return NULL;
492  _assert_(tria->FiniteElement()==P1Enum);
493  int index1=tria->GetVertexIndex(vertices[0]);
494  int index2=tria->GetVertexIndex(vertices[1]);
495 
496  /*Intermediaries */
497  IssmDouble Jdet,v2,Afactor,Bfactor,fFactor;
498  IssmDouble A,B,n,phi_old,phi,phi_0,dphimds,dphi[2];
499  IssmDouble H,h,b,db[2],dphids,qc,dPw,ks,Ngrad;
500  IssmDouble xyz_list[NUMVERTICES][3];
501  IssmDouble xyz_list_tria[3][3];
502  const int numnodes = NUMNODES;
503 
504  /*Initialize Element vector and other vectors*/
505  ElementVector* pe = new ElementVector(this->nodes,NUMNODES,this->parameters);
506  IssmDouble basis[NUMNODES];
507 
508  /*Retrieve all inputs and parameters*/
509  GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
510  GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
511 
512  IssmDouble L = element->FindParam(MaterialsLatentheatEnum);
513  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
514  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
515  IssmDouble kc = element->FindParam(HydrologyChannelConductivityEnum);
516  IssmDouble g = element->FindParam(ConstantsGEnum);
517  IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum);
518  IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum);
519 
520  Input2* h_input = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
521  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
522  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
523  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
524  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
525  Input2* ks_input = element->GetInput2(HydrologySheetConductivityEnum); _assert_(ks_input);
526  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
527 
528  /*Get tangent vector*/
529  IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
530  IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
531  IssmDouble Lt = sqrt(tx*tx+ty*ty);
532  tx = tx/Lt;
533  ty = ty/Lt;
534 
535  /* Start looping on the number of gaussian points: */
536  Gauss* gauss=new GaussTria(index1,index2,2);
537  for(int ig=gauss->begin();ig<gauss->end();ig++){
538  gauss->GaussPoint(ig);
539 
540  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
541  tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
542 
543  /*Get input values at gauss points*/
544  b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
545  phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
546  h_input->GetInputValue(&h,gauss);
547  ks_input->GetInputValue(&ks,gauss);
548  B_input->GetInputValue(&B,gauss);
549  n_input->GetInputValue(&n,gauss);
550  phi_input->GetInputValue(&phi,gauss);
551  b_input->GetInputValue(&b,gauss);
552  H_input->GetInputValue(&H,gauss);
553 
554  /*Get values for a few potentials*/
555  phi_0 = rho_water*g*b + rho_ice*g*H;
556  dphids = dphi[0]*tx + dphi[1]*ty;
557  dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
558  Ngrad = fabs(dphids);
559  if(Ngrad<AEPS) Ngrad = AEPS;
560 
561  /*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
562  IssmDouble Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.);
563 
564  /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
565  qc = - Ks * dphids;
566 
567  /*d(phi - phi_m)/ds*/
568  dPw = dphids - dphimds;
569 
570  /*Compute f factor*/
571  fFactor = 0.;
572  if(this->S>0. || qc*dPw>0.){
573  fFactor = lc * qc;
574  }
575 
576  /*Compute Afactor and Bfactor*/
577  Afactor = C_W*c_t*rho_water;
578  Bfactor = 1./L * (1./rho_ice - 1./rho_water);
579 
580  /*Compute closing rate*/
581  /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
582  A=pow(B,-n);
583  v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
584 
585  for(int i=0;i<numnodes;i++){
586  pe->values[i]+= - Jdet*gauss->weight*(-v2)*basis[i];
587  pe->values[i]+= + Jdet*gauss->weight*Afactor*Bfactor*fFactor*dphimds*basis[i];
588  }
589  }
590 
591  /*Clean up and return*/
592  delete gauss;
593  return pe;
594 }
595 /*}}}*/
596 void Channel::SetChannelCrossSectionOld(void){/*{{{*/
597 
598  this->Sold = this->S;
599 
600 } /*}}}*/
601 void Channel::UpdateChannelCrossSection(void){/*{{{*/
602 
603  /*Initialize Element matrix and return if necessary*/
604  Tria* tria=(Tria*)element;
605  if(!tria->IsIceInElement() || this->boundary){
606  this->S = 0.;
607  return;
608  }
609  _assert_(tria->FiniteElement()==P1Enum);
610  int index1=tria->GetVertexIndex(vertices[0]);
611  int index2=tria->GetVertexIndex(vertices[1]);
612 
613  /*Intermediaries */
614  IssmDouble A,B,n,phi,phi_0,ks,Ngrad;
615  IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
616  IssmDouble xyz_list[NUMVERTICES][3];
617  IssmDouble xyz_list_tria[3][3];
618 
619  /*Retrieve all inputs and parameters*/
620  GetVerticesCoordinates(&xyz_list[0][0] ,this->vertices,NUMVERTICES);
621  GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
622 
623  IssmDouble L = element->FindParam(MaterialsLatentheatEnum);
624  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
625  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
626  IssmDouble g = element->FindParam(ConstantsGEnum);
627  IssmDouble kc = element->FindParam(HydrologyChannelConductivityEnum);
628  IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum);
629  IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum);
630  IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
631 
632  Input2* h_input = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
633  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
634  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
635  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
636  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
637  Input2* ks_input = element->GetInput2(HydrologySheetConductivityEnum); _assert_(ks_input);
638  Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
639 
640  /*Get tangent vector*/
641  IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
642  IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
643  IssmDouble Lt = sqrt(tx*tx+ty*ty);
644  tx = tx/Lt;
645  ty = ty/Lt;
646 
647  /*Evaluate fields on center of edge*/
648  GaussTria* gauss=new GaussTria();
649  gauss->GaussEdgeCenter(index1,index2);
650 
651  /*Get input values at gauss points*/
652  phi_input->GetInputValue(&phi,gauss);
653  phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
654  h_input->GetInputValue(&h,gauss);
655  ks_input->GetInputValue(&ks,gauss);
656  B_input->GetInputValue(&B,gauss);
657  n_input->GetInputValue(&n,gauss);
658  b_input->GetInputValue(&b,gauss);
659  b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
660  H_input->GetInputValue(&H,gauss);
661 
662  /*Get values for a few potentials*/
663  phi_0 = rho_water*g*b + rho_ice*g*H;
664  dphids = dphi[0]*tx + dphi[1]*ty;
665  dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
666  Ngrad = fabs(dphids);
667  if(Ngrad<AEPS) Ngrad = AEPS;
668 
669  /*d(phi - phi_m)/ds*/
670  IssmDouble dPw = dphids - dphimds;
671 
672  /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
673  IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
674 
675  /*Ice rate factor*/
676  A=pow(B,-n);
677 
678  IssmDouble C = C_W*c_t*rho_water;
679  IssmDouble Qprime = -kc * pow(Ngrad,BETA_C-2.)*dphids;
680  IssmDouble N = phi_0 - phi;
681 
682  bool converged = false;
683  int count = 0;
684 
685  while(!converged){
686 
687  IssmDouble Snew = this->S;
688 
689  /*Compute f factor*/
690  IssmDouble fFactor = 0.;
691  if(this->S>0. || qc*dPw>0.){
692  fFactor = lc * qc;
693  }
694 
695  IssmDouble alpha = 1./(rho_ice*L)*(
696  fabs(Qprime*pow(Snew,ALPHA_C-1.)*dphids)
697  + C*Qprime*pow(Snew,ALPHA_C-1.)*dPw
698  ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
699 
700  IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
701 
702  /*Solve ODE*/
703  this->S = ODE1(alpha,beta,this->Sold,dt,2);
704  if(this->S<0.) this->S = 0.;
705  _assert_(!xIsNan<IssmDouble>(this->S));
706 
707  count++;
708 
709  if(fabs((this->S - Snew)/(Snew+AEPS))<1e-8 || count>=10) converged = true;
710  }
711 
712  /*Clean up and return*/
713  delete gauss;
714 }
715 /*}}}*/
716 void Channel::WriteChannelCrossSection(IssmPDouble* values){/*{{{*/
717 
718  _assert_(values);
719  values[this->sid] = reCast<IssmPDouble>(this->S);
720 }
721 /*}}}*/
Matrix< IssmDouble >
Channel::GetNumberOfNodes
int GetNumberOfNodes(void)
Definition: Channel.cpp:258
Channel::GetNodesSidList
void GetNodesSidList(int *sidlist)
Definition: Channel.cpp:250
Vertices
Declaration of Vertices class.
Definition: Vertices.h:15
BETA_C
#define BETA_C
Definition: Channel.cpp:23
Channel::ResetHooks
void ResetHooks()
Definition: Channel.cpp:280
Channel
Definition: Channel.h:18
_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::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
Channel::CreatePVectorHydrologyGlaDS
ElementVector * CreatePVectorHydrologyGlaDS(void)
Definition: Channel.cpp:487
TriaRef::GetSegmentJacobianDeterminant
void GetSegmentJacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)
Definition: TriaRef.cpp:229
Channel::hvertices
Hook * hvertices
Definition: Channel.h:32
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Channel::hnodes
Hook * hnodes
Definition: Channel.h:31
HydraulicPotentialEnum
@ HydraulicPotentialEnum
Definition: EnumDefinitions.h:597
Channel::CreateKMatrixHydrologyGlaDS
ElementMatrix * CreateKMatrixHydrologyGlaDS(void)
Definition: Channel.cpp:354
Hook::deliverp
Object ** deliverp(void)
Definition: Hook.cpp:187
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
Channel::S
IssmDouble S
Definition: Channel.h:21
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
MARSHALLING_ENUM
#define MARSHALLING_ENUM(EN)
Definition: Marshalling.h:14
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
Hook::DeepEcho
void DeepEcho(void)
Definition: Hook.cpp:77
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
Channel::id
int id
Definition: Channel.h:27
Channel::Sold
IssmDouble Sold
Definition: Channel.h:22
Elements
Declaration of Elements class.
Definition: Elements.h:17
HydrologyPressureMeltCoefficientEnum
@ HydrologyPressureMeltCoefficientEnum
Definition: EnumDefinitions.h:171
SsetEnum
@ SsetEnum
Definition: EnumDefinitions.h:1282
ALPHA_S
#define ALPHA_S
Definition: Channel.cpp:25
MaterialsLatentheatEnum
@ MaterialsLatentheatEnum
Definition: EnumDefinitions.h:254
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
Element::vertices
Vertex ** vertices
Definition: Element.h:49
Channel::Channel
Channel()
Definition: Channel.cpp:30
Channel::GetNodesLidList
void GetNodesLidList(int *lidlist)
Definition: Channel.cpp:242
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
BETA_S
#define BETA_S
Definition: Channel.cpp:26
TriaRef::GetSegmentNodalFunctions
void GetSegmentNodalFunctions(IssmDouble *basis, Gauss *gauss, int index1, int index2, int finiteelement)
Definition: TriaRef.cpp:243
Hook::reset
void reset(void)
Definition: Hook.cpp:211
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
HydrologyGlaDSAnalysisEnum
@ HydrologyGlaDSAnalysisEnum
Definition: EnumDefinitions.h:1100
Element
Definition: Element.h:41
Channel::ObjectEnum
int ObjectEnum()
Definition: Channel.cpp:174
Channel::element
Element * element
Definition: Channel.h:35
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Channel::vertices
Vertex ** vertices
Definition: Channel.h:36
Channel::PenaltyCreatePVector
void PenaltyCreatePVector(Vector< IssmDouble > *pf, IssmDouble kmax)
Definition: Channel.cpp:273
GaussTria
Definition: GaussTria.h:12
ChannelEnum
@ ChannelEnum
Definition: EnumDefinitions.h:1008
HydrologySheetConductivityEnum
@ HydrologySheetConductivityEnum
Definition: EnumDefinitions.h:620
Object
Definition: Object.h:13
Channel::SetwiseNodeConnectivity
void SetwiseNodeConnectivity(int *d_nz, int *o_nz, Node *node, bool *flags, int *flagsindices, int set1_enum, int set2_enum)
Definition: Channel.cpp:298
Parameters::DeepEcho
void DeepEcho()
Definition: Parameters.cpp:99
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
Hook::delivers
Object * delivers(void)
Definition: Hook.cpp:191
Materials
Declaration of Materials class.
Definition: Materials.h:16
Channel::Configure
void Configure(Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
Definition: Channel.cpp:179
Channel::IsPenalty
bool IsPenalty(void)
Definition: Channel.cpp:262
Channel::Echo
void Echo()
Definition: Channel.cpp:134
HydrologySheetThicknessEnum
@ HydrologySheetThicknessEnum
Definition: EnumDefinitions.h:621
Channel::copy
Object * copy()
Definition: Channel.cpp:93
Hook
Definition: Hook.h:16
Hook::configure
void configure(DataSet *dataset)
Definition: Hook.cpp:145
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
Channel::sid
int sid
Definition: Channel.h:26
Channel::boundary
bool boundary
Definition: Channel.h:23
Channel::CreateKMatrix
void CreateKMatrix(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: Channel.cpp:196
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
MARSHALLING
#define MARSHALLING(FIELD)
Definition: Marshalling.h:29
Channel::CreatePVector
void CreatePVector(Vector< IssmDouble > *pf)
Definition: Channel.cpp:219
Channel::SetChannelCrossSectionOld
void SetChannelCrossSectionOld(void)
Definition: Channel.cpp:596
GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble *xyz, Vertex **vertices, int numvertices, bool spherical)
Definition: Vertex.cpp:225
Channel::SetCurrentConfiguration
void SetCurrentConfiguration(Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
Definition: Channel.cpp:294
C_W
#define C_W
Definition: Channel.cpp:21
Tria::FiniteElement
int FiniteElement(void)
Definition: Tria.cpp:1318
Hook::copy
Object * copy(void)
Definition: Hook.cpp:61
Input2
Definition: Input2.h:18
MARSHALLING_BACKWARD
@ MARSHALLING_BACKWARD
Definition: Marshalling.h:10
NUMVERTICES
#define NUMVERTICES
Definition: Channel.cpp:19
Channel::~Channel
~Channel()
Definition: Channel.cpp:85
Channel::DeepEcho
void DeepEcho()
Definition: Channel.cpp:119
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
ODE1
IssmDouble ODE1(IssmDouble alpha, IssmDouble beta, IssmDouble Si, IssmDouble dt, int method)
Definition: ODE1.cpp:5
Loads
Declaration of Loads class.
Definition: Loads.h:16
Node
Definition: Node.h:23
Node::Lid
int Lid(void)
Definition: Node.cpp:618
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Tria
Definition: Tria.h:29
NUMNODES
#define NUMNODES
Definition: Channel.cpp:18
HydrologyChannelConductivityEnum
@ HydrologyChannelConductivityEnum
Definition: EnumDefinitions.h:164
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
IoModel::faces
int * faces
Definition: IoModel.h:88
Channel::parameters
Parameters * parameters
Definition: Channel.h:38
Gauss::begin
virtual int begin(void)=0
ALPHA_C
#define ALPHA_C
Definition: Channel.cpp:22
HydrologyChannelSheetWidthEnum
@ HydrologyChannelSheetWidthEnum
Definition: EnumDefinitions.h:165
TriaRef::GetSegmentNodalFunctionsDerivatives
void GetSegmentNodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list_tria, Gauss *gauss, int index1, int index2, int finiteelement)
Definition: TriaRef.cpp:286
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
GaussTria::GaussEdgeCenter
void GaussEdgeCenter(int index1, int index2)
Definition: GaussTria.cpp:423
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
ElementVector::AddToGlobal
void AddToGlobal(Vector< IssmDouble > *pf)
Definition: ElementVector.cpp:155
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
Channel::helement
Hook * helement
Definition: Channel.h:30
ElementVector
Definition: ElementVector.h:20
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
Hook::Echo
void Echo(void)
Definition: Hook.cpp:104
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Vertex
Definition: Vertex.h:19
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
IoModel
Definition: IoModel.h:48
FsetEnum
@ FsetEnum
Definition: EnumDefinitions.h:1075
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
DataSet
Declaration of DataSet class.
Definition: DataSet.h:14
Channel::WriteChannelCrossSection
void WriteChannelCrossSection(IssmPDouble *values)
Definition: Channel.cpp:716
shared.h
ElementMatrix
Definition: ElementMatrix.h:19
Channel::PenaltyCreateKMatrix
void PenaltyCreateKMatrix(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *kfs, IssmDouble kmax)
Definition: Channel.cpp:266
Channel::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Channel.cpp:148
Channel::nodes
Node ** nodes
Definition: Channel.h:37
Vector< IssmDouble >
Channel::Id
int Id()
Definition: Channel.cpp:144
Channel::UpdateChannelCrossSection
void UpdateChannelCrossSection(void)
Definition: Channel.cpp:601
Gauss
Definition: Gauss.h:8
Tria::GetVertexIndex
int GetVertexIndex(Vertex *vertex)
Definition: Tria.cpp:2368
AEPS
#define AEPS
Definition: Channel.cpp:27
Hook::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Hook.cpp:122
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Generated on Thu Jul 2 2020 08:09:18 for Ice Sheet System Model by   doxygen 1.8.19