source: issm/oecreview/Archive/21724-22754/ISSM-21825-21826.diff@ 22755

Last change on this file since 22755 was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 21.8 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/c/shared/Elements/elements.h
2===================================================================
3--- ../trunk-jpl/src/c/shared/Elements/elements.h (revision 21825)
4+++ ../trunk-jpl/src/c/shared/Elements/elements.h (revision 21826)
5@@ -16,7 +16,7 @@
6 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
7 IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm);
8 void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag);
9-void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
10+void EstarStrainrateQuantities(IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
11 IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
12 IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm,
13 IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t,
14Index: ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp
15===================================================================
16--- ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21825)
17+++ ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21826)
18@@ -2,12 +2,12 @@
19 #include "../Numerics/types.h"
20 #include "../Exceptions/exceptions.h"
21 #include "./elements.h"
22-void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
23+void EstarStrainrateQuantities(IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
24
25 /*Intermediaries*/
26 IssmDouble omega[3];
27 IssmDouble nrsp[3],nrsp_norm;
28- IssmDouble eps[3][3],epseff;
29+ IssmDouble eps[3][3];
30 IssmDouble epsprime[3],epsprime_norm;
31
32 /*Get omega, correction for rigid body rotation*/
33@@ -36,11 +36,6 @@
34 eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]);
35 eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];
36
37- /*effective strain rate*/
38- epseff = 0.;
39- /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
40- epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] + eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]);
41-
42 /*Compute the shear strain rate on the non rotating shear plane*/
43 epsprime[0]=0.;
44 epsprime[1]=0.;
45@@ -72,7 +67,6 @@
46 epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
47
48 /*Assign output pointers*/
49- *pepseff=epseff;
50 *pepsprime_norm=epsprime_norm;
51 }/*}}}*/
52 void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag){/*{{{*/
53Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
54===================================================================
55--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21825)
56+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21826)
57@@ -47,7 +47,7 @@
58 /*Intermediaries*/
59 IssmDouble vx,vy,vz,vmag;
60 IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3];
61- IssmDouble epseff,epsprime;
62+ IssmDouble eps[3][3],epseff,epsprime;
63 int dim;
64 IssmDouble *xyz_list = NULL;
65
66@@ -100,7 +100,18 @@
67 dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
68 dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
69 }
70- EstarStrainrateQuantities(&epseff,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
71+
72+ /*Build strain rate tensor*/
73+ eps[0][0] = dvx[0]; eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]);
74+ eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]);
75+ eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];
76+
77+ /*effective strain rate*/
78+ epseff = 0.;
79+ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
80+ epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] + eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]);
81+
82+ EstarStrainrateQuantities(&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
83 lambdas[iv]=EstarLambdaS(epseff,epsprime);
84 }
85
86@@ -676,7 +687,7 @@
87 material->GetViscosity_B(&dmudB,eps_eff);
88 break;
89 case MatestarEnum:
90- material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
91+ material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input,eps_eff);
92 break;
93 default: _error_("not supported");
94 }
95@@ -713,7 +724,7 @@
96 material->GetViscosity_B(&dmudB,eps_eff);
97 break;
98 case MatestarEnum:
99- material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
100+ material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff);
101 break;
102 default: _error_("not supported");
103 }
104@@ -750,7 +761,7 @@
105 material->GetViscosity_B(&dmudB,eps_eff);
106 break;
107 case MatestarEnum:
108- material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
109+ material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff);
110 break;
111 default: _error_("not supported");
112 }
113Index: ../trunk-jpl/src/c/classes/Materials/Matice.h
114===================================================================
115--- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21825)
116+++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21826)
117@@ -87,9 +87,9 @@
118 void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
119 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
120 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
121- void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
122- void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
123- void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
124+ void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){_error_("not supported");};
125+ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");};
126+ void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");};
127 /*}}}*/
128 };
129
130Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
131===================================================================
132--- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21825)
133+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21826)
134@@ -124,9 +124,9 @@
135 void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");};
136 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
137 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
138- void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
139- void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
140- void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
141+ void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){_error_("not supported");};
142+ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");};
143+ void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");};
144 /*}}}*/
145 /*Numerics: {{{*/
146 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
147Index: ../trunk-jpl/src/c/classes/Materials/Matestar.cpp
148===================================================================
149--- ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21825)
150+++ ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21826)
151@@ -248,13 +248,13 @@
152 _error_("not implemented yet");
153 }
154 /*}}}*/
155-IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
156+IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged){/*{{{*/
157
158 /*output: */
159 IssmDouble viscosity;
160
161 /*Intermediaries*/
162- IssmDouble epseff,epsprime_norm;
163+ IssmDouble epsprime_norm;
164 IssmDouble lambdas;
165 IssmDouble vmag,dvmag[3];
166 IssmDouble B,Ec,Es,E,n;
167@@ -272,8 +272,8 @@
168 dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
169 }
170
171- EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
172- lambdas=EstarLambdaS(epseff,epsprime_norm);
173+ EstarStrainrateQuantities(&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
174+ lambdas=EstarLambdaS(eps_eff,epsprime_norm);
175
176 /*Get B and enhancement*/
177 n=GetN(); _assert_(n>0.);
178@@ -293,24 +293,27 @@
179
180 /*Compute viscosity*/
181 /*if no strain rate, return maximum viscosity*/
182- if(epseff==0.){
183+ if(eps_eff==0.){
184 viscosity = 1.e+14/2.;
185 //viscosity = B;
186 //viscosity=2.5*pow(10.,17);
187 }
188 else{
189- viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n));
190+ viscosity = B/(2.*pow(E,1./n)*pow(eps_eff,2./n));
191 }
192
193+ /*Checks in debugging mode*/
194+ if(viscosity<=0) _error_("Negative viscosity");
195+
196 /*Assign output pointer*/
197 return viscosity;
198 }
199 /*}}}*/
200-IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
201+IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged){/*{{{*/
202
203 /*Intermediaries*/
204 IssmDouble dmudB;
205- IssmDouble epseff,epsprime_norm;
206+ IssmDouble epsprime_norm;
207 IssmDouble lambdas;
208 IssmDouble vmag,dvmag[3];
209 IssmDouble Ec,Es,E;
210@@ -328,8 +331,8 @@
211 dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
212 }
213
214- EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
215- lambdas=EstarLambdaS(epseff,epsprime_norm);
216+ EstarStrainrateQuantities(&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
217+ lambdas=EstarLambdaS(eps_eff,epsprime_norm);
218
219 /*Get enhancement*/
220 if (isdepthaveraged==0.){
221@@ -345,8 +348,8 @@
222 E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
223
224 /*Compute dmudB*/
225- if(epseff==0.) dmudB = 0.;
226- else dmudB = 1./(2.*pow(E,1./3.)*pow(epseff,2./3.));
227+ if(eps_eff==0.) dmudB = 0.;
228+ else dmudB = 1./(2.*pow(E,1./3.)*pow(eps_eff,2./3.));
229
230 /*Assign output*/
231 return dmudB;
232@@ -407,7 +410,7 @@
233
234 }
235 /*}}}*/
236-void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
237+void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){/*{{{*/
238
239 /*Intermediaries*/
240 IssmDouble vx,vy,vz;
241@@ -433,10 +436,10 @@
242 }
243
244 /*Compute dmudB*/
245- *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
246+ *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged);
247 }
248 /*}}}*/
249-void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
250+void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){/*{{{*/
251
252 /*Intermediaries*/
253 IssmDouble vx,vy,vz;
254@@ -462,9 +465,9 @@
255 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
256
257 /*Compute viscosity*/
258- *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
259+ *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged);
260 }/*}}}*/
261-void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
262+void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){/*{{{*/
263 /*Intermediaries*/
264 IssmDouble vx,vy,vz;
265 IssmDouble dvx[3],dvy[3],dvz[3];
266@@ -492,7 +495,7 @@
267 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
268
269 /*Compute viscosity*/
270- *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
271+ *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged);
272 }/*}}}*/
273 void Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
274
275@@ -499,6 +502,9 @@
276 /*Intermediaries*/
277 IssmDouble vx,vy,vz;
278 IssmDouble dvx[3],dvy[3],dvz[3];
279+ IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
280+ IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];*/
281+ IssmDouble eps_eff,eps0=1.e-27;
282 bool isdepthaveraged=0.;
283
284 /*Get velocity derivatives in all directions*/
285@@ -519,8 +525,19 @@
286 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
287 }
288
289+ if(dim==3){
290+ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
291+ element->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
292+ eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
293+ }
294+ else{
295+ /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
296+ element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
297+ eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
298+ }
299+
300 /*Compute viscosity*/
301- *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
302+ *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged);
303 }
304 /*}}}*/
305 void Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
306@@ -531,8 +548,22 @@
307 /*Intermediaries*/
308 IssmDouble vx,vy,vz;
309 IssmDouble dvx[3],dvy[3],dvz[3];
310+ IssmDouble epsilon3d[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
311+ IssmDouble epsilon2d[5]; /* epsilon=[exx,exy];*/
312+ IssmDouble eps_eff;
313 bool isdepthaveraged=0.;
314
315+ if(dim==3){
316+ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
317+ element->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
318+ eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
319+ }
320+ else{
321+ /* eps_eff^2 = 1/2 (2*exx^2 + 2*exy^2 ) (since eps_zz = - eps_xx)*/
322+ element->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
323+ eps_eff = 1./sqrt(2.)*sqrt(2*epsilon2d[0]*epsilon2d[0] + 2*epsilon2d[1]*epsilon2d[1]);
324+ }
325+
326 /*Get velocity derivatives in all directions*/
327 _assert_(dim==2 || dim==3);
328 _assert_(vx_input);
329@@ -552,7 +583,7 @@
330 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
331
332 /*Compute viscosity*/
333- *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
334+ *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged);
335 }/*}}}*/
336 void Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
337 _error_("not implemented yet");
338@@ -565,6 +596,9 @@
339 /*Intermediaries*/
340 IssmDouble vx,vy,vz;
341 IssmDouble dvx[3],dvy[3],dvz[3];
342+ IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */
343+ IssmDouble epsilon1d; /* epsilon=[exx]; */
344+ IssmDouble eps_eff;
345 bool isdepthaveraged=1.;
346
347 /*Get velocity derivatives in all directions*/
348@@ -588,8 +622,19 @@
349 vz = 0.;
350 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
351
352+ if(dim==2){
353+ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
354+ element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
355+ eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
356+ }
357+ else{
358+ /* eps_eff^2 = exx^2*/
359+ element->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
360+ eps_eff = fabs(epsilon1d);
361+ }
362+
363 /*Compute viscosity*/
364- *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
365+ *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged);
366 }/*}}}*/
367 void Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
368 _error_("not implemented yet");
369Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
370===================================================================
371--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21825)
372+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21826)
373@@ -723,7 +723,7 @@
374 IssmDouble viscosity;
375 IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
376 IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */
377- IssmDouble eps_eff,E=1.0;
378+ IssmDouble eps_eff;
379
380 if(dim==3){
381 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
382Index: ../trunk-jpl/src/c/classes/Materials/Material.h
383===================================================================
384--- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21825)
385+++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21826)
386@@ -52,9 +52,9 @@
387 virtual void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
388 virtual void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
389 virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
390- virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
391- virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
392- virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
393+ virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble epseff)=0;
394+ virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble epseff)=0;
395+ virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble epseff)=0;
396
397 };
398 #endif
399Index: ../trunk-jpl/src/c/classes/Materials/Matestar.h
400===================================================================
401--- ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21825)
402+++ ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21826)
403@@ -85,12 +85,12 @@
404 void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
405 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
406 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
407- void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
408- void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
409- void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
410+ void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff);
411+ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff);
412+ void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff);
413 /*}}}*/
414- IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
415- IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
416+ IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged);
417+ IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged);
418 };
419
420 #endif /* _MATESTAR_H_ */
Note: See TracBrowser for help on using the repository browser.