Ice Sheet System Model  4.18
Code documentation
Friction.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 "../classes.h"
14 #include "shared/shared.h"
15 /*}}}*/
16 
17 /*Constructors/destructors*/
18 Friction::Friction(){/*{{{*/
19  this->element=NULL;
20  this->dim=0;
21  this->law=0;
22 
23 }
24 /*}}}*/
25 Friction::Friction(Element* element_in,int dim_in){/*{{{*/
26 
27  this->element=element_in;
28  this->dim=dim_in;
29  element_in->FindParam(&this->law,FrictionLawEnum);
30 }
31 /*}}}*/
32 Friction::~Friction(){/*{{{*/
33 }
34 /*}}}*/
35 
36 /*methods: */
37 void Friction::Echo(void){/*{{{*/
38  _printf_("Friction:\n");
39  _printf_(" dim: " << this->dim<< "\n");
40 }
41 /*}}}*/
42 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
43 
44  switch(this->law){
45  case 1:
46  GetAlphaViscousComplement(palpha_complement,gauss);
47  break;
48  case 3:
49  GetAlphaHydroComplement(palpha_complement,gauss);
50  break;
51  case 4:
52  GetAlphaTempComplement(palpha_complement,gauss);
53  break;
54  default:
55  _error_("not supported");
56  }
57 
58 
59  /*Checks*/
60  _assert_(!xIsNan<IssmDouble>(*palpha_complement));
61  _assert_(!xIsInf<IssmDouble>(*palpha_complement));
62 
63 }/*}}}*/
64 void Friction::GetAlphaHydroComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
65 
66  /*diverse: */
67  IssmDouble q_exp;
68  IssmDouble C_param;
69  IssmDouble As;
70  IssmDouble n;
71  IssmDouble alpha;
72  IssmDouble Chi,Gamma;
73  IssmDouble alpha_complement;
74 
75  /*Recover parameters: */
76  element->GetInputValue(&q_exp,gauss,FrictionQEnum);
77  element->GetInputValue(&C_param,gauss,FrictionCEnum);
78  element->GetInputValue(&As,gauss,FrictionAsEnum);
79  element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
80 
81  /*Get effective pressure and velocity magnitude*/
82  IssmDouble Neff = EffectivePressure(gauss);
83  IssmDouble vmag = VelMag(gauss);
84 
85  if (q_exp==1){
86  alpha=1;
87  }
88  else{
89  alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
90  }
91  Chi = vmag/(pow(C_param,n)*pow(Neff,n)*As);
92  Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
93  /*Check to prevent dividing by zero if vmag==0*/
94  if(vmag==0.) alpha_complement=0.;
95  else if(Neff==0.) alpha_complement=0.;
96  else alpha_complement=-(C_param*Neff/(n*vmag)) *
97  pow(Gamma,((1.-n)/n)) *
98  (Gamma/As - (alpha*q_exp*pow(Chi,q_exp-1.)* Gamma * Gamma/As));
99 
100  /*Assign output pointers:*/
101  *palpha_complement=alpha_complement;
102 }
103 /*}}}*/
104 void Friction::GetAlphaTempComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
105  /*Here, we want to parameterize the friction as a function of temperature
106  *
107  * alpha2 = alpha2_viscous * 1/f(T)
108  *
109  * where f(T) = exp((T-Tpmp)/gamma)
110  */
111 
112  /*Intermediaries: */
113  IssmDouble f,T,pressure,Tpmp,gamma;
114  IssmDouble alpha_complement;
115 
116  /*Get viscous part*/
117  this->GetAlphaViscousComplement(&alpha_complement,gauss);
118 
119  /*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
120  element->GetInputValue(&T,gauss,TemperatureEnum);
121  element->GetInputValue(&pressure,gauss,PressureEnum);
122  Tpmp = element->TMeltingPoint(pressure);
123 
124  /*Compute scaling parameter*/
125  element->parameters->FindParam(&gamma,FrictionGammaEnum);
126  alpha_complement = alpha_complement/ (exp((T-Tpmp)/gamma)+1e-3);
127 
128  /*Assign output pointers:*/
129  *palpha_complement=alpha_complement;
130 }/*}}}*/
131 void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
132 
133  /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
134  * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
135  * alpha_complement= Neff ^r * vel ^s*/
136 
137  /*diverse: */
138  IssmDouble r,s;
139  IssmDouble drag_p,drag_q;
140  IssmDouble drag_coefficient;
141  IssmDouble alpha_complement;
142 
143  /*Recover parameters: */
144  element->GetInputValue(&drag_p,gauss,FrictionPEnum);
145  element->GetInputValue(&drag_q,gauss,FrictionQEnum);
146  element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
147 
148  /*compute r and q coefficients: */
149  r=drag_q/drag_p;
150  s=1./drag_p;
151 
152  /*Get effective pressure*/
153  IssmDouble Neff = EffectivePressure(gauss);
154  IssmDouble vmag = VelMag(gauss);
155 
156  /*Check to prevent dividing by zero if vmag==0*/
157  if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
158  else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
159 
160  /*Assign output pointers:*/
161  *palpha_complement=alpha_complement;
162 }
163 /*}}}*/
164 void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
165 
166  switch(this->law){
167  case 1:
168  GetAlpha2Viscous(palpha2,gauss);
169  break;
170  case 2:
171  GetAlpha2Weertman(palpha2,gauss);
172  break;
173  case 3:
174  GetAlpha2Hydro(palpha2,gauss);
175  break;
176  case 4:
177  GetAlpha2Temp(palpha2,gauss);
178  break;
179  case 5:
180  GetAlpha2WaterLayer(palpha2,gauss);
181  break;
182  case 6:
183  GetAlpha2WeertmanTemp(palpha2,gauss);
184  break;
185  case 7:
186  GetAlpha2Coulomb(palpha2,gauss);
187  break;
188  case 8:
189  GetAlpha2Shakti(palpha2,gauss);
190  break;
191  case 9:
192  GetAlpha2Josh(palpha2,gauss);
193  break;
194  case 10:
195  GetAlpha2PISM(palpha2,gauss);
196  break;
197  case 11:
198  GetAlpha2Schoof(palpha2,gauss);
199  break;
200  case 12:
201  GetAlpha2Tsai(palpha2,gauss);
202  break;
203  default:
204  _error_("Friction law "<< this->law <<" not supported");
205  }
206 
207  /*Checks*/
208  _assert_(!xIsNan<IssmDouble>(*palpha2));
209  _assert_(!xIsInf<IssmDouble>(*palpha2));
210  _assert_(*palpha2>=0);
211 
212 }/*}}}*/
213 void Friction::GetAlpha2Coulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
214 
215  /*This routine calculates the basal friction coefficient
216  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
217 
218  /*diverse: */
219  IssmDouble r,s;
220  IssmDouble drag_p, drag_q;
221  IssmDouble thickness,base,floatation_thickness,sealevel;
222  IssmDouble drag_coefficient,drag_coefficient_coulomb;
223  IssmDouble alpha2,alpha2_coulomb;
224 
225  /*Recover parameters: */
226  element->GetInputValue(&drag_p,gauss,FrictionPEnum);
227  element->GetInputValue(&drag_q,gauss,FrictionQEnum);
228  element->GetInputValue(&thickness, gauss,ThicknessEnum);
229  element->GetInputValue(&base, gauss,BaseEnum);
230  element->GetInputValue(&sealevel, gauss,SealevelEnum);
231  element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
232  element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
233  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
234  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
235  IssmDouble gravity = element->FindParam(ConstantsGEnum);
236 
237  //compute r and q coefficients: */
238  r=drag_q/drag_p;
239  s=1./drag_p;
240 
241  /*Get effective pressure*/
242  IssmDouble Neff = EffectivePressure(gauss);
243  IssmDouble vmag = VelMag(gauss);
244 
245  /*Check to prevent dividing by zero if vmag==0*/
246  if(vmag==0. && (s-1.)<0.){
247  alpha2=0.;
248  }
249  else{
250  alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
251  }
252 
253  floatation_thickness=0;
254  if(base<0) floatation_thickness=-(rho_water/rho_ice)*base;
255  if(vmag==0.){
256  alpha2_coulomb=0.;
257  }
258  else{
259  alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
260  }
261 
262  if(alpha2_coulomb<alpha2) alpha2=alpha2_coulomb;
263 
264  /*Assign output pointers:*/
265  *palpha2=alpha2;
266 }/*}}}*/
267 void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
268 
269  /*This routine calculates the basal friction coefficient
270  Based on Gagliardini 2007, needs a good effective pressure computation
271  Not tested so far so use at your own risks
272  alpha2= NeffC[Chi/(1+alpha*Chi^q)]^(1/n)*|vel|^(-1) with
273  Chi=|vel|/(C^n*Neff^n*As)
274  alpha=(q-1)^(q-1)/q^q */
275 
276  /*diverse: */
277  IssmDouble q_exp;
278  IssmDouble C_param;
279  IssmDouble As;
280  IssmDouble n;
281  IssmDouble alpha;
282  IssmDouble Chi,Gamma;
283  IssmDouble alpha2;
284 
285  /*Recover parameters: */
286  element->GetInputValue(&q_exp,gauss,FrictionQEnum);
287  element->GetInputValue(&C_param,gauss,FrictionCEnum);
288  element->GetInputValue(&As,gauss,FrictionAsEnum);
289  element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
290 
291  /*Get effective pressure*/
292  IssmDouble Neff = EffectivePressure(gauss);
293  IssmDouble vmag = VelMag(gauss);
294 
295  //compute alpha and Chi coefficients: */
296  if (q_exp==1){
297  alpha=1;
298  }
299  else{
300  alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
301  }
302  Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
303  Gamma=(Chi/(1. + alpha * pow(Chi,q_exp)));
304  /*Check to prevent dividing by zero if vmag==0*/
305  if(vmag==0.) alpha2=0.;
306  else if (Neff==0) alpha2=0.0;
307  else alpha2=Neff * C_param * pow(Gamma,1./n) * pow(vmag,-1);
308 
309  /*Assign output pointers:*/
310  *palpha2=alpha2;
311 }/*}}}*/
312 void Friction::GetAlpha2Shakti(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
313 
314  /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-base)*/
315 
316  /*diverse: */
317  IssmDouble pressure_ice,pressure_water;
318  IssmDouble Neff;
319  IssmDouble drag_coefficient;
320  IssmDouble base,thickness,head,sealevel;
321  IssmDouble alpha2;
322 
323  /*Recover parameters: */
324  element->GetInputValue(&thickness, gauss,ThicknessEnum);
325  element->GetInputValue(&base, gauss,BaseEnum);
326  element->GetInputValue(&head, gauss,HydrologyHeadEnum);
327  element->GetInputValue(&sealevel, gauss,SealevelEnum);
328  element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
329  IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
330  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
331  IssmDouble gravity = element->FindParam(ConstantsGEnum);
332 
333  //From base and thickness, compute effective pressure when drag is viscous:
334  pressure_ice = rho_ice*gravity*thickness;
335  pressure_water = rho_water*gravity*(head-base+sealevel);
336  Neff=pressure_ice-pressure_water;
337  if(Neff<0.) Neff=0.;
338 
339  alpha2=drag_coefficient*drag_coefficient*Neff;
340 
341  /*Assign output pointers:*/
342  *palpha2=alpha2;
343 }
344 /*}}}*/
345 void Friction::GetAlpha2Temp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
346  /*Here, we want to parameterize the friction as a function of temperature
347  *
348  * alpha2 = alpha2_viscous * 1/f(T)
349  *
350  * where f(T) = exp((T-Tpmp)/gamma)
351  */
352 
353  /*Intermediaries: */
354  IssmDouble f,T,pressure,Tpmp,gamma;
355  IssmDouble alpha2;
356 
357  /*Get viscous part*/
358  this->GetAlpha2Viscous(&alpha2,gauss);
359 
360  /*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
361  element->GetInputValue(&T,gauss,TemperatureEnum);
362  element->GetInputValue(&pressure,gauss,PressureEnum);
363  Tpmp = element->TMeltingPoint(pressure);
364 
365  /*Compute scaling parameter*/
366  element->parameters->FindParam(&gamma,FrictionGammaEnum);
367  alpha2 = alpha2 / (exp((T-Tpmp)/gamma)+1e-3);
368 
369  /*Assign output pointers:*/
370  *palpha2=alpha2;
371 }/*}}}*/
372 void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
373  /*Here, we want to parameterize the friction as a function of temperature
374  *
375  * alpha2 = alpha2_viscous * 1/f(T)
376  *
377  * where f(T) = exp((T-Tpmp)/gamma)
378  */
379 
380  /*Intermediaries: */
381  IssmDouble T,Tpmp,deltaT,deltaTref,pressure,diff,drag_coefficient;
382  IssmDouble alpha2,time,gamma,ref,alp_new,alphascaled;
383  const IssmDouble yts = 365*24*3600.;
384 
385  /*Get viscous part*/
386  this->GetAlpha2Viscous(&alpha2,gauss);
387 
388  /*Get delta Refs*/
389  element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum);
390  element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
391  /*New*/
392  /*element->GetInputValue(&deltaTrefsfc,gauss,FrictionSurfaceTemperatureEnum);
393  * element->GetInputValue(&Tpdd,gauss,TemperaturePDDEnum);
394  * */
395 
396  /*Compute delta T*/
397  element->GetInputValue(&T,gauss,TemperatureEnum);
398  element->GetInputValue(&pressure,gauss,PressureEnum);
399  Tpmp = element->TMeltingPoint(pressure);
400  deltaT = T-Tpmp;
401 
402 
403  /*Compute gamma*/
404  element->parameters->FindParam(&time,TimeEnum);
405  element->parameters->FindParam(&gamma,FrictionGammaEnum);
406 
407  ref = exp(deltaTref/gamma);
408  alp_new = ref/exp(deltaT/gamma);
409 
410  alphascaled = sqrt(alp_new)*drag_coefficient;
411  if (alphascaled > 300) alp_new = (300/drag_coefficient)*(300/drag_coefficient);
412 
413  alp_new=alp_new*alpha2;
414 
415  /*Assign output pointers:*/
416  *palpha2=alp_new;
417 }/*}}}*/
418 void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
419 
420  /*This routine calculates the basal friction coefficient
421  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
422 
423  /*diverse: */
424  IssmDouble r,s;
425  IssmDouble drag_p, drag_q;
426  IssmDouble drag_coefficient;
427  IssmDouble alpha2;
428 
429  /*Recover parameters: */
430  element->GetInputValue(&drag_p,gauss,FrictionPEnum);
431  element->GetInputValue(&drag_q,gauss,FrictionQEnum);
432  element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
433 
434  /*compute r and q coefficients: */
435  r=drag_q/drag_p;
436  s=1./drag_p;
437 
438  /*Get effective pressure and basal velocity*/
439  IssmDouble Neff = EffectivePressure(gauss);
440  IssmDouble vmag = VelMag(gauss);
441 
442  /*Check to prevent dividing by zero if vmag==0*/
443  if(vmag==0. && (s-1.)<0.) alpha2=0.;
444  else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
445 
446  /*Assign output pointers:*/
447  *palpha2=alpha2;
448 }/*}}}*/
449 void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
450 
451  /*This routine calculates the basal friction coefficient
452  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
453 
454  /*diverse: */
455  IssmDouble r,s;
456  IssmDouble drag_p, drag_q;
457  IssmDouble Neff,F;
458  IssmDouble thickness,base,sealevel;
459  IssmDouble drag_coefficient,water_layer;
460  IssmDouble alpha2;
461 
462  /*Recover parameters: */
463  element->parameters->FindParam(&F,FrictionFEnum);
464  element->GetInputValue(&drag_p,gauss,FrictionPEnum);
465  element->GetInputValue(&drag_q,gauss,FrictionQEnum);
466  element->GetInputValue(&thickness, gauss,ThicknessEnum);
467  element->GetInputValue(&base, gauss,BaseEnum);
468  element->GetInputValue(&sealevel, gauss,SealevelEnum);
469  element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
470  element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
471  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
472  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
473  IssmDouble gravity = element->FindParam(ConstantsGEnum);
474 
475  //compute r and q coefficients: */
476  r=drag_q/drag_p;
477  s=1./drag_p;
478 
479  //From base and thickness, compute effective pressure when drag is viscous:
480  if(base>0) base=0;
481  if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(base-sealevel);
482  else if(water_layer>0) Neff=gravity*rho_ice*thickness*F;
483  else _error_("negative water layer thickness");
484  if(Neff<0) Neff=0;
485 
486  IssmDouble vmag = VelMag(gauss);
487 
488  /*Check to prevent dividing by zero if vmag==0*/
489  if(vmag==0. && (s-1.)<0.) alpha2=0.;
490  else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
491 
492  /*Assign output pointers:*/
493  *palpha2=alpha2;
494 }/*}}}*/
495 void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
496 
497  /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
498 
499  /*diverse: */
500  IssmDouble C,m;
501  IssmDouble alpha2;
502 
503  /*Recover parameters: */
504  element->GetInputValue(&C,gauss,FrictionCEnum);
505  element->GetInputValue(&m,gauss,FrictionMEnum);
506 
507  /*Get velocity magnitude*/
508  IssmDouble vmag = VelMag(gauss);
509 
510  /*Check to prevent dividing by zero if vmag==0*/
511  if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
512  else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
513 
514  /*Assign output pointers:*/
515  *palpha2=alpha2;
516 }/*}}}*/
517 void Friction::GetAlpha2WeertmanTemp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
518  /*Here, we want to parameterize the friction as a function of temperature
519  *
520  * alpha2 = alpha2_weertman * 1/f(T)
521  *
522  * where f(T) = exp((T-Tpmp)/gamma)
523  */
524 
525  /*Intermediaries: */
526  IssmDouble f,T,pressure,Tpmp,gamma;
527  IssmDouble alpha2;
528 
529  /*Get viscous part*/
530  this->GetAlpha2Weertman(&alpha2,gauss);
531 
532  /*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
533  element->GetInputValue(&T,gauss,TemperatureEnum);
534  element->GetInputValue(&pressure,gauss,PressureEnum);
535  Tpmp = element->TMeltingPoint(pressure);
536 
537  /*Compute scaling parameter*/
538  element->parameters->FindParam(&gamma,FrictionGammaEnum);
539  alpha2 = alpha2 / exp((T-Tpmp)/gamma);
540 
541  /*Assign output pointers:*/
542  *palpha2=alpha2;
543 }/*}}}*/
544 void Friction::GetAlpha2PISM(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
545  /*Here, we want to parameterize the friction using a pseudoplastic friction law,
546  * computing the basal shear stress as
547  *
548  * alpha2 = tau_c (u_b/(abs(u_b)^(1-q)*u_0^q))
549  *
550  * The yield stress tau_c is a function of the effective pressure N
551  * using a Mohr-Coloumb criterion, so that
552  * tau_c = tan(phi)*N,
553  * where phi is the till friction angle, representing sediment strength
554  *
555  * The effective pressure is given by Eq. (5) in Aschwanden et al. 2016:
556  *
557  * N = delta * P0 * 10^((e_0/Cc)(1-(W/Wmax)))
558  *
559  * W is calculated by a non-conserving hydrology model in HydrologyPismAnalysis.cpp
560  *
561  * see Aschwanden et al. 2016 and Bueler and Brown, 2009 for more details
562  */
563 
564  /*compute ice overburden pressure P0*/
565  IssmDouble thickness,base,P0;
566  element->GetInputValue(&thickness, gauss,ThicknessEnum);
567  element->GetInputValue(&base, gauss,BaseEnum);
568  //element->GetInputValue(&sealevel, gauss,SealevelEnum);
569  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
570  IssmDouble gravity = element->FindParam(ConstantsGEnum);
571  P0 = gravity*rho_ice*thickness;
572 
573  /*Compute effective pressure*/
574  IssmDouble N,delta,W,Wmax,e0,Cc;
575  element->parameters->FindParam(&delta,FrictionDeltaEnum);
576  element->parameters->FindParam(&e0,FrictionVoidRatioEnum);
577  element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum);
578  element->GetInputValue(&W,gauss,WatercolumnEnum);
579  element->GetInputValue(&Wmax,gauss,HydrologyWatercolumnMaxEnum);
580 
581  /*Check that water column height is within 0 and upper bound, correct if needed*/
582  if(W>Wmax) W=Wmax;
583  if(W<0) W=0.;
584 
585  N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));
586 
587  /*Get till friction angles, defined by user [deg]*/
588  IssmDouble phi;
589  element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);
590 
591  /*Convert till friction angle from user-defined deg to rad, which Matlab uses*/
592  phi = phi*PI/180.;
593 
594  /*Compute yield stress following a Mohr-Colomb criterion*/
595  IssmDouble tau_c = N*tan(phi);
596 
597  /*Compute basal speed*/
598  IssmDouble ub;
599  element->GetInputValue(&ub,gauss,VelEnum);
600 
601  /*now compute alpha^2*/
602  IssmDouble u0,q;
603  element->parameters->FindParam(&u0,FrictionThresholdSpeedEnum);
604  element->parameters->FindParam(&q,FrictionPseudoplasticityExponentEnum);
605  IssmDouble alpha2 = tau_c/(pow(ub+1.e-10,1.-q)*pow(u0,q));
606 
607  /*Assign output pointers:*/
608  *palpha2=alpha2;
609 }/*}}}*/
610 void Friction::GetAlpha2Schoof(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
611 
612  /*This routine calculates the basal friction coefficient
613  *
614  * C |u_b|^(m-1)
615  * alpha2= __________________________
616  * (1+(C/(Cmax N))^1/m |u_b| )^m
617  *
618  * */
619 
620  /*diverse: */
621  IssmDouble C,Cmax,m,alpha2;
622 
623  /*Recover parameters: */
624  element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum);
625  element->GetInputValue(&C,gauss,FrictionCEnum);
626  element->GetInputValue(&m,gauss,FrictionMEnum);
627 
628  /*Get effective pressure and velocity magnitude*/
629  IssmDouble N = EffectivePressure(gauss);
630  IssmDouble ub = VelMag(gauss);
631 
632  /*Compute alpha^2*/
633  if(ub<1e-10){
634  alpha2 = 0.;
635  }
636  else{
637  alpha2= (C*pow(ub,m-1.)) / pow(1.+ pow(C/(Cmax*N),1./m)*ub,m);
638  }
639 
640  /*Assign output pointers:*/
641  *palpha2=alpha2;
642 }/*}}}*/
643 void Friction::GetAlpha2Tsai(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
644 
645  /*This routine calculates the basal friction coefficient
646  *
647  * alpha2= min(C |ub|^m , f N ) / |ub|
648  *
649  * */
650 
651  /*diverse: */
652  IssmDouble C,f,m,alpha2;
653 
654  /*Recover parameters: */
655  element->GetInputValue(&f,gauss,FrictionfEnum);
656  element->GetInputValue(&C,gauss,FrictionCEnum);
657  element->GetInputValue(&m,gauss,FrictionMEnum);
658 
659  /*Get effective pressure and velocity magnitude*/
660  IssmDouble N = EffectivePressure(gauss);
661  IssmDouble ub = VelMag(gauss);
662 
663  /*Compute alpha^2*/
664  if(ub<1e-10){
665  alpha2 = 0.;
666  }
667  else{
668  alpha2= C*pow(ub,m);
669 
670  if(alpha2>f*N) alpha2 = f*N;
671 
672  alpha2 = alpha2/ub;
673  }
674 
675  /*Assign output pointers:*/
676  *palpha2=alpha2;
677 }/*}}}*/
678 
679 IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/
680  /*Get effective pressure as a function of flag */
681 
682  /*diverse: */
683  int coupled_flag;
684  IssmDouble thickness,base,sealevel;
685  IssmDouble p_ice,p_water;
686  IssmDouble Neff,Neff_limit;
687 
688  /*Recover parameters: */
689  element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
690  element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
691 
692  /*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
693  switch(coupled_flag){
694  case 0:{
695  element->GetInputValue(&thickness, gauss,ThicknessEnum);
696  element->GetInputValue(&base, gauss,BaseEnum);
697  element->GetInputValue(&sealevel, gauss,SealevelEnum);
698  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
699  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
700  IssmDouble gravity = element->FindParam(ConstantsGEnum);
701  p_ice = gravity*rho_ice*thickness;
702  p_water = rho_water*gravity*(sealevel-base);
703  Neff = p_ice - p_water;
704  }
705  break;
706  case 1:{
707  element->GetInputValue(&thickness, gauss,ThicknessEnum);
708  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
709  IssmDouble gravity = element->FindParam(ConstantsGEnum);
710  p_ice = gravity*rho_ice*thickness;
711  p_water = 0.;
712  Neff = p_ice - p_water;
713  }
714  break;
715  case 2:{
716  element->GetInputValue(&thickness, gauss,ThicknessEnum);
717  element->GetInputValue(&base, gauss,BaseEnum);
718  element->GetInputValue(&sealevel, gauss,SealevelEnum);
719  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
720  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
721  IssmDouble gravity = element->FindParam(ConstantsGEnum);
722  p_ice = gravity*rho_ice*thickness;
723  p_water = max(0.,rho_water*gravity*(sealevel-base));
724  Neff = p_ice - p_water;
725  }
726  break;
727  case 3:{
728  element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
729  element->GetInputValue(&thickness, gauss,ThicknessEnum);
730  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
731  IssmDouble gravity = element->FindParam(ConstantsGEnum);
732  p_ice = gravity*rho_ice*thickness;
733  }
734  break;
735  case 4:{
736  element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
737  element->GetInputValue(&thickness, gauss,ThicknessEnum);
738  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
739  IssmDouble gravity = element->FindParam(ConstantsGEnum);
740  p_ice = gravity*rho_ice*thickness;
741  }
742  break;
743  default:
744  _error_("not supported");
745  }
746 
747  /*Make sure Neff is positive*/
748  if(Neff<Neff_limit*p_ice) Neff=Neff_limit*p_ice;
749 
750  /*Return effective pressure*/
751  return Neff;
752 
753 }/*}}}*/
754 IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/
755  /*Get effective pressure as a function of flag */
756 
757 
758  /*diverse*/
759  IssmDouble vx,vy,vz,vmag;
760 
761  switch(dim){
762  case 1:
763  element->GetInputValue(&vx,gauss,VxEnum);
764  vmag=sqrt(vx*vx);
765  break;
766  case 2:
767  element->GetInputValue(&vx,gauss,VxEnum);
768  element->GetInputValue(&vy,gauss,VyEnum);
769  vmag=sqrt(vx*vx+vy*vy);
770  break;
771  case 3:
772  element->GetInputValue(&vx,gauss,VxEnum);
773  element->GetInputValue(&vy,gauss,VyEnum);
774  element->GetInputValue(&vz,gauss,VzEnum);
775  vmag=sqrt(vx*vx+vy*vy+vz*vz);
776  break;
777  default:
778  _error_("not supported");
779  }
780 
781  return vmag;
782 
783 }/*}}}*/
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Friction::GetAlphaTempComplement
void GetAlphaTempComplement(IssmDouble *alpha_complement, Gauss *gauss)
Definition: Friction.cpp:104
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
Friction::VelMag
IssmDouble VelMag(Gauss *gauss)
Definition: Friction.cpp:754
IssmDouble
double IssmDouble
Definition: types.h:37
Friction::GetAlphaViscousComplement
void GetAlphaViscousComplement(IssmDouble *alpha_complement, Gauss *gauss)
Definition: Friction.cpp:131
FrictionPEnum
@ FrictionPEnum
Definition: EnumDefinitions.h:578
FrictionFEnum
@ FrictionFEnum
Definition: EnumDefinitions.h:146
FrictionEffectivePressureEnum
@ FrictionEffectivePressureEnum
Definition: EnumDefinitions.h:576
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Friction::GetAlphaComplement
void GetAlphaComplement(IssmDouble *alpha_complement, Gauss *gauss)
Definition: Friction.cpp:42
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
FrictionGammaEnum
@ FrictionGammaEnum
Definition: EnumDefinitions.h:147
Friction::GetAlpha2Shakti
void GetAlpha2Shakti(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:312
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
FrictionQEnum
@ FrictionQEnum
Definition: EnumDefinitions.h:580
Friction::Echo
void Echo(void)
Definition: Friction.cpp:37
Friction::GetAlphaHydroComplement
void GetAlphaHydroComplement(IssmDouble *alpha_complement, Gauss *gauss)
Definition: Friction.cpp:64
Friction::GetAlpha2WaterLayer
void GetAlpha2WaterLayer(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:449
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
Friction::~Friction
~Friction()
Definition: Friction.cpp:32
FrictionCouplingEnum
@ FrictionCouplingEnum
Definition: EnumDefinitions.h:143
Friction::GetAlpha2PISM
void GetAlpha2PISM(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:544
Friction::GetAlpha2Temp
void GetAlpha2Temp(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:345
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
FrictionAsEnum
@ FrictionAsEnum
Definition: EnumDefinitions.h:571
Element
Definition: Element.h:41
FrictionThresholdSpeedEnum
@ FrictionThresholdSpeedEnum
Definition: EnumDefinitions.h:150
Friction::GetAlpha2Weertman
void GetAlpha2Weertman(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:495
FrictionLawEnum
@ FrictionLawEnum
Definition: EnumDefinitions.h:148
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
FrictionCmaxEnum
@ FrictionCmaxEnum
Definition: EnumDefinitions.h:573
FrictionPseudoplasticityExponentEnum
@ FrictionPseudoplasticityExponentEnum
Definition: EnumDefinitions.h:149
Friction::GetAlpha2Josh
void GetAlpha2Josh(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:372
FrictionPressureAdjustedTemperatureEnum
@ FrictionPressureAdjustedTemperatureEnum
Definition: EnumDefinitions.h:579
Friction::GetAlpha2Coulomb
void GetAlpha2Coulomb(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:213
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
Friction::element
Element * element
Definition: Friction.h:16
MaterialsRhoSeawaterEnum
@ MaterialsRhoSeawaterEnum
Definition: EnumDefinitions.h:265
PI
const double PI
Definition: constants.h:11
FrictionDeltaEnum
@ FrictionDeltaEnum
Definition: EnumDefinitions.h:144
Element::TMeltingPoint
IssmDouble TMeltingPoint(IssmDouble pressure)
Definition: Element.cpp:4583
TemperatureEnum
@ TemperatureEnum
Definition: EnumDefinitions.h:831
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
FrictionSedimentCompressibilityCoefficientEnum
@ FrictionSedimentCompressibilityCoefficientEnum
Definition: EnumDefinitions.h:581
Friction::Friction
Friction()
Definition: Friction.cpp:18
FrictionfEnum
@ FrictionfEnum
Definition: EnumDefinitions.h:584
Friction::GetAlpha2Hydro
void GetAlpha2Hydro(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:267
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
HydrologyWatercolumnMaxEnum
@ HydrologyWatercolumnMaxEnum
Definition: EnumDefinitions.h:623
Element::parameters
Parameters * parameters
Definition: Element.h:51
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
Friction::law
int law
Definition: Friction.h:18
WatercolumnEnum
@ WatercolumnEnum
Definition: EnumDefinitions.h:859
FrictionEffectivePressureLimitEnum
@ FrictionEffectivePressureLimitEnum
Definition: EnumDefinitions.h:145
FrictionTillFrictionAngleEnum
@ FrictionTillFrictionAngleEnum
Definition: EnumDefinitions.h:582
Element::GetInputValue
void GetInputValue(bool *pvalue, int enum_type)
Definition: Element.cpp:1177
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
Friction::GetAlpha2
void GetAlpha2(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:164
Friction::GetAlpha2Tsai
void GetAlpha2Tsai(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:643
FrictionCoefficientcoulombEnum
@ FrictionCoefficientcoulombEnum
Definition: EnumDefinitions.h:575
VelEnum
@ VelEnum
Definition: EnumDefinitions.h:844
Friction::dim
int dim
Definition: Friction.h:17
FrictionMEnum
@ FrictionMEnum
Definition: EnumDefinitions.h:577
Friction::GetAlpha2WeertmanTemp
void GetAlpha2WeertmanTemp(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:517
EffectivePressureEnum
@ EffectivePressureEnum
Definition: EnumDefinitions.h:548
HydrologyHeadEnum
@ HydrologyHeadEnum
Definition: EnumDefinitions.h:615
FrictionVoidRatioEnum
@ FrictionVoidRatioEnum
Definition: EnumDefinitions.h:151
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
FrictionWaterLayerEnum
@ FrictionWaterLayerEnum
Definition: EnumDefinitions.h:583
Friction::GetAlpha2Schoof
void GetAlpha2Schoof(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:610
shared.h
FrictionCEnum
@ FrictionCEnum
Definition: EnumDefinitions.h:572
Friction::GetAlpha2Viscous
void GetAlpha2Viscous(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:418
Friction::EffectivePressure
IssmDouble EffectivePressure(Gauss *gauss)
Definition: Friction.cpp:679
Gauss
Definition: Gauss.h:8
Generated on Thu Jul 2 2020 08:09:18 for Ice Sheet System Model by   doxygen 1.8.19