source: issm/oecreview/Archive/21337-21723/ISSM-21699-21700.diff@ 21726

Last change on this file since 21726 was 21726, checked in by Mathieu Morlighem, 8 years ago

CHG added Archive/21337-21723
File size: 15.4 KB
RevLine 
[21726]1Index: ../trunk-jpl/src/c/classes/Materials/Matestar.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21699)
4+++ ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21700)
5@@ -159,6 +159,7 @@
6 }
7 /*}}}*/
8 IssmDouble Matestar::GetB(){/*{{{*/
9+
10 /*Output*/
11 IssmDouble B;
12
13@@ -167,6 +168,7 @@
14 }
15 /*}}}*/
16 IssmDouble Matestar::GetBbar(){/*{{{*/
17+
18 /*Output*/
19 IssmDouble Bbar;
20
21@@ -192,6 +194,15 @@
22 return Ec;
23 }
24 /*}}}*/
25+IssmDouble Matestar::GetEcbar(){/*{{{*/
26+
27+ /*Output*/
28+ IssmDouble Ecbar;
29+
30+ element->inputs->GetInputAverage(&Ecbar,MaterialsRheologyEcbarEnum);
31+ return Ecbar;
32+}
33+/*}}}*/
34 IssmDouble Matestar::GetEs(){/*{{{*/
35
36 /*Output*/
37@@ -201,6 +212,15 @@
38 return Es;
39 }
40 /*}}}*/
41+IssmDouble Matestar::GetEsbar(){/*{{{*/
42+
43+ /*Output*/
44+ IssmDouble Esbar;
45+
46+ element->inputs->GetInputAverage(&Esbar,MaterialsRheologyEsbarEnum);
47+ return Esbar;
48+}
49+/*}}}*/
50 IssmDouble Matestar::GetN(){/*{{{*/
51
52 /*Output*/
53@@ -210,51 +230,6 @@
54 /*}}}*/
55 void Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
56 _error_("not implemented yet");
57- /*From a string tensor and a material object, return viscosity, using Glen's flow law.
58- (1-D) B
59- viscosity= -------------------------
60- 2 eps_eff ^[(n-1)/n]
61-
62- where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
63- and n the flow law exponent.
64-
65- If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
66- return 10^14, initial viscosity.
67- */
68-
69- /*output: */
70- IssmDouble viscosity;
71-
72- /*Intermediary: */
73- IssmDouble B,D=0.,n;
74-
75- /*Get B and n*/
76- B=GetB(); _assert_(B>0.);
77- n=GetN(); _assert_(n>0.);
78-
79- if (n==1.){
80- /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
81- viscosity=(1.-D)*B/2.;
82- }
83- else{
84-
85- /*if no strain rate, return maximum viscosity*/
86- if(eps_eff==0.){
87- viscosity = 1.e+14/2.;
88- //viscosity = B;
89- //viscosity=2.5*pow(10.,17);
90- }
91-
92- else{
93- viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
94- }
95- }
96-
97- /*Checks in debugging mode*/
98- if(viscosity<=0) _error_("Negative viscosity");
99-
100- /*Return: */
101- *pviscosity=viscosity;
102 }
103 /*}}}*/
104 void Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
105@@ -273,14 +248,16 @@
106 _error_("not implemented yet");
107 }
108 /*}}}*/
109-IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
110+IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
111
112+ /*output: */
113+ IssmDouble viscosity;
114+
115 /*Intermediaries*/
116- IssmDouble viscosity;
117 IssmDouble epseff,epsprime_norm;
118 IssmDouble lambdas;
119 IssmDouble vmag,dvmag[3];
120- IssmDouble B,Ec,Es,E;
121+ IssmDouble B,Ec,Es,E,n;
122
123 /*Calculate velocity magnitude and its derivative*/
124 vmag = sqrt(vx*vx+vy*vy+vz*vz);
125@@ -299,9 +276,17 @@
126 lambdas=EstarLambdaS(epseff,epsprime_norm);
127
128 /*Get B and enhancement*/
129- B=GetB(); _assert_(B>0.);
130- Ec=GetEc(); _assert_(Ec>=0.);
131- Es=GetEs(); _assert_(Es>=0.);
132+ n=GetN(); _assert_(n>0.);
133+ if (isdepthaveraged==0.){
134+ B=GetB(); _assert_(B>0.);
135+ Ec=GetEc(); _assert_(Ec>=0.);
136+ Es=GetEs(); _assert_(Es>=0.);
137+ }
138+ else{
139+ B=GetBbar(); _assert_(B>0.);
140+ Ec=GetEcbar(); _assert_(Ec>=0.);
141+ Es=GetEsbar(); _assert_(Es>=0.);
142+ }
143
144 /*Get total enhancement factor E(lambdas)*/
145 E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
146@@ -314,34 +299,64 @@
147 //viscosity=2.5*pow(10.,17);
148 }
149 else{
150- viscosity = B/(2.*pow(E*epseff*epseff,1./3.));
151+ viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n));
152 }
153
154 /*Assign output pointer*/
155 return viscosity;
156 }
157 /*}}}*/
158-void Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
159- /*output: */
160+IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
161+
162+ /*Intermediaries*/
163 IssmDouble dmudB;
164+ IssmDouble epseff,epsprime_norm;
165+ IssmDouble lambdas;
166+ IssmDouble vmag,dvmag[3];
167+ IssmDouble Ec,Es,E;
168
169- /*Intermediary: */
170- IssmDouble E=1.,n;
171+ /*Calculate velocity magnitude and its derivative*/
172+ vmag = sqrt(vx*vx+vy*vy+vz*vz);
173+ if(vmag<1e-12){
174+ dvmag[0]=0;
175+ dvmag[1]=0;
176+ dvmag[2]=0;
177+ }
178+ else{
179+ dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
180+ dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
181+ dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
182+ }
183
184- n=GetN(); _assert_(n>0.);
185- if(n==1.){
186- /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */
187- dmudB=1./(2.*E);
188+ EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
189+ lambdas=EstarLambdaS(epseff,epsprime_norm);
190+
191+ /*Get enhancement*/
192+ if (isdepthaveraged==0.){
193+ Ec=GetEc(); _assert_(Ec>=0.);
194+ Es=GetEs(); _assert_(Es>=0.);
195 }
196 else{
197- if(eps_eff==0.) dmudB = 0.;
198- else dmudB = 1./(2.*pow(E*eps_eff*eps_eff,1./3.));
199+ Ec=GetEcbar(); _assert_(Ec>=0.);
200+ Es=GetEsbar(); _assert_(Es>=0.);
201 }
202
203- /*Return: */
204- *pdmudB=dmudB;
205+ /*Get total enhancement factor E(lambdas)*/
206+ E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
207+
208+ /*Compute dmudB*/
209+ if(epseff==0.) dmudB = 0.;
210+ else dmudB = 1./(2.*pow(E,1./3.)*pow(epseff,2./3.));
211+
212+ /*Assign output*/
213+ return dmudB;
214+
215 }
216 /*}}}*/
217+void Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
218+ _error_("not implemented yet");
219+}
220+/*}}}*/
221 void Matestar::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/
222 _error_("not implemented yet");
223 }
224@@ -392,12 +407,99 @@
225
226 }
227 /*}}}*/
228+void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
229+
230+ /*Intermediaries*/
231+ IssmDouble vx,vy,vz;
232+ IssmDouble dvx[3],dvy[3],dvz[3];
233+ bool isdepthaveraged=0.;
234+
235+ /*Get velocity derivatives in all directions*/
236+ _assert_(dim>1);
237+ _assert_(vx_input);
238+ vx_input->GetInputValue(&vx,gauss);
239+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
240+ _assert_(vy_input);
241+ vy_input->GetInputValue(&vy,gauss);
242+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
243+ if(dim==3){
244+ _assert_(vz_input);
245+ vz_input->GetInputValue(&vz,gauss);
246+ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
247+ }
248+ else{
249+ vz = 0.;
250+ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
251+ }
252+
253+ /*Compute dmudB*/
254+ *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
255+}
256+/*}}}*/
257+void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
258+
259+ /*Intermediaries*/
260+ IssmDouble vx,vy,vz;
261+ IssmDouble dvx[3],dvy[3],dvz[3];
262+ bool isdepthaveraged=0.;
263+
264+ /*Get velocity derivatives in all directions*/
265+ _assert_(dim==2 || dim==3);
266+ _assert_(vx_input);
267+ vx_input->GetInputValue(&vx,gauss);
268+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
269+ if(dim==3){
270+ _assert_(vy_input);
271+ vy_input->GetInputValue(&vy,gauss);
272+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
273+ }
274+ else{
275+ dvx[2] = 0.;
276+ vy = 0.;
277+ dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
278+ }
279+ vz = 0.;
280+ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
281+
282+ /*Compute viscosity*/
283+ *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
284+}/*}}}*/
285+void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
286+ /*Intermediaries*/
287+ IssmDouble vx,vy,vz;
288+ IssmDouble dvx[3],dvy[3],dvz[3];
289+ bool isdepthaveraged=1.;
290+
291+ /*Get velocity derivatives in all directions*/
292+ _assert_(dim==1 || dim==2);
293+ _assert_(vx_input);
294+ vx_input->GetInputValue(&vx,gauss);
295+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
296+ if(dim==2){
297+ _assert_(vy_input);
298+ vy_input->GetInputValue(&vy,gauss);
299+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
300+ }
301+ else{
302+ dvx[1] = 0.;
303+ dvx[2] = 0.;
304+ vy = 0.;
305+ dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
306+ }
307+ dvx[2] = 0.;
308+ dvy[2] = 0.;
309+ vz = 0.;
310+ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
311+
312+ /*Compute viscosity*/
313+ *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
314+}/*}}}*/
315 void Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
316
317 /*Intermediaries*/
318 IssmDouble vx,vy,vz;
319 IssmDouble dvx[3],dvy[3],dvz[3];
320- IssmDouble B,Ec,Es;
321+ bool isdepthaveraged=0.;
322
323 /*Get velocity derivatives in all directions*/
324 _assert_(dim>1);
325@@ -418,7 +520,7 @@
326 }
327
328 /*Compute viscosity*/
329- *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
330+ *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
331 }
332 /*}}}*/
333 void Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
334@@ -429,6 +531,7 @@
335 /*Intermediaries*/
336 IssmDouble vx,vy,vz;
337 IssmDouble dvx[3],dvy[3],dvz[3];
338+ bool isdepthaveraged=0.;
339
340 /*Get velocity derivatives in all directions*/
341 _assert_(dim==2 || dim==3);
342@@ -449,7 +552,7 @@
343 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
344
345 /*Compute viscosity*/
346- *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
347+ *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
348 }/*}}}*/
349 void Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
350 _error_("not implemented yet");
351@@ -458,10 +561,11 @@
352 _error_("not implemented yet");
353 }/*}}}*/
354 void Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
355+
356 /*Intermediaries*/
357 IssmDouble vx,vy,vz;
358 IssmDouble dvx[3],dvy[3],dvz[3];
359- IssmDouble B,Ec,Es;
360+ bool isdepthaveraged=1.;
361
362 /*Get velocity derivatives in all directions*/
363 _assert_(dim==1 || dim==2);
364@@ -485,7 +589,7 @@
365 dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
366
367 /*Compute viscosity*/
368- *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
369+ *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
370 }/*}}}*/
371 void Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
372 _error_("not implemented yet");
373Index: ../trunk-jpl/src/c/classes/Materials/Material.h
374===================================================================
375--- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21699)
376+++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21700)
377@@ -52,6 +52,9 @@
378 virtual void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
379 virtual void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
380 virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
381+ virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
382+ virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
383+ virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
384
385 };
386 #endif
387Index: ../trunk-jpl/src/c/classes/Materials/Matestar.h
388===================================================================
389--- ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21699)
390+++ ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21700)
391@@ -69,7 +69,9 @@
392 IssmDouble GetD();
393 IssmDouble GetDbar();
394 IssmDouble GetEc();
395+ IssmDouble GetEcbar();
396 IssmDouble GetEs();
397+ IssmDouble GetEsbar();
398 IssmDouble GetN();
399 bool IsDamage();
400 bool IsEnhanced(){_error_("not supported");};
401@@ -83,8 +85,12 @@
402 void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
403 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
404 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
405+ void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
406+ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
407+ void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
408 /*}}}*/
409- IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
410+ IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
411+ IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
412 };
413
414 #endif /* _MATESTAR_H_ */
415Index: ../trunk-jpl/src/c/classes/Materials/Matice.h
416===================================================================
417--- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21699)
418+++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21700)
419@@ -87,6 +87,9 @@
420 void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
421 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
422 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
423+ void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
424+ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
425+ void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
426 /*}}}*/
427 };
428
429Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
430===================================================================
431--- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21699)
432+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21700)
433@@ -124,6 +124,9 @@
434 void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");};
435 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
436 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
437+ void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
438+ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
439+ void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
440 /*}}}*/
441 /*Numerics: {{{*/
442 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
Note: See TracBrowser for help on using the repository browser.