source: issm/oecreview/Archive/21337-21723/ISSM-21550-21551.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: 18.8 KB
RevLine 
[21726]1Index: ../trunk-jpl/src/m/classes/frictionjosh.m
2===================================================================
3--- ../trunk-jpl/src/m/classes/frictionjosh.m (revision 21550)
4+++ ../trunk-jpl/src/m/classes/frictionjosh.m (revision 21551)
5@@ -6,12 +6,12 @@
6 classdef frictionjosh
7 properties (SetAccess=public)
8 coefficient = NaN;
9- temperature = NaN;
10+ pressure_adjusted_temperature = NaN;
11 end
12 methods
13 function self = extrude(self,md) % {{{
14 self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1);
15- self.temperature=project3d(md,'vector',self.temperature,'type','node','layer',1);
16+ self.pressure_adjusted_temperature=project3d(md,'vector',self.pressure_adjusted_temperature,'type','node','layer',1);
17 end % }}}
18 function self = frictionjosh(varargin) % {{{
19 switch nargin
20@@ -31,8 +31,8 @@
21 %Early return
22 if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end
23
24- md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1);
25- md = checkfield(md,'fieldname','friction.temperature','timeseries',1,'NaN',1,'Inf',1);
26+ md = checkfield(md,'fieldname','friction.coefficient','NaN',1,'Inf',1);
27+ md = checkfield(md,'fieldname','friction.pressure_adjusted_temperature','NaN',1,'Inf',1);
28
29 %Check that temperature is provided
30 md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
31@@ -40,13 +40,13 @@
32 function disp(self) % {{{
33 disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)'));
34 fielddisplay(self,'coefficient','friction coefficient [SI]');
35- fielddisplay(self,'temperature','friction temperature [K]');
36+ fielddisplay(self,'pressure_adjusted_temperature','friction pressure_adjusted_temperature (T - Tpmp) [K]');
37 end % }}}
38 function marshall(self,prefix,md,fid) % {{{
39
40 WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer');
41 WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
42- WriteData(fid,prefix,'class','friction','object',self,'fieldname','temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
43+ WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
44 end % }}}
45 end
46 end
47Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
48===================================================================
49--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21550)
50+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21551)
51@@ -2032,13 +2032,17 @@
52 yearlytemperatures[2] = s[2];
53 this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
54
55- /*Convert that to enthalpy for the enthalpy model*/
56- IssmDouble enthalpy[6];
57- GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
58- ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
59- ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.);
60- ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.);
61- this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum));
62+ bool isenthalpy;
63+ this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
64+ if(isenthalpy){
65+ /*Convert that to enthalpy for the enthalpy model*/
66+ IssmDouble enthalpy[6];
67+ GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
68+ ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
69+ ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.);
70+ ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.);
71+ this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum));
72+ }
73 }
74 this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum));
75 this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
76Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp
77===================================================================
78--- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 21550)
79+++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 21551)
80@@ -246,6 +246,9 @@
81 case 8:
82 GetAlpha2Sommers(palpha2,gauss);
83 break;
84+ case 9:
85+ GetAlpha2Josh(palpha2,gauss);
86+ break;
87 default:
88 _error_("Friction law "<< this->law <<" not supported");
89 }
90@@ -462,6 +465,47 @@
91 /*Assign output pointers:*/
92 *palpha2=alpha2;
93 }/*}}}*/
94+void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
95+ /*Here, we want to parameterize the friction as a function of temperature
96+ *
97+ * alpha2 = alpha2_viscous * 1/f(T)
98+ *
99+ * where f(T) = exp((T-Tpmp)/gamma)
100+ */
101+
102+ /*Intermediaries: */
103+ IssmDouble T,Tpmp,deltaT,deltaTref,pressure;
104+ IssmDouble alpha2,time,gamma;
105+ const IssmDouble yts = 365*24*3600.;
106+
107+ /*Get viscous part*/
108+ this->GetAlpha2Viscous(&alpha2,gauss);
109+
110+ /*Get delta Refs*/
111+ element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum);
112+
113+ /*Compute delta T*/
114+ element->GetInputValue(&T,gauss,TemperatureEnum);
115+ element->GetInputValue(&pressure,gauss,PressureEnum);
116+ Tpmp = element->TMeltingPoint(pressure);
117+ deltaT = T-Tpmp;
118+
119+ /*Compute gamma*/
120+ element->parameters->FindParam(&time,TimeEnum);
121+ if(time<25e3*yts){
122+ gamma = 10.;
123+ }
124+ else{
125+ gamma = 5.;
126+ }
127+ gamma = 5.;
128+
129+ /*Compute scaling parameter*/
130+ alpha2 = alpha2 * exp((deltaTref - deltaT)/(2*gamma));
131+
132+ /*Assign output pointers:*/
133+ *palpha2=alpha2;
134+}/*}}}*/
135 void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
136
137 /*This routine calculates the basal friction coefficient
138Index: ../trunk-jpl/src/c/classes/Loads/Friction.h
139===================================================================
140--- ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 21550)
141+++ ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 21551)
142@@ -35,6 +35,7 @@
143 void GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
144 void GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss);
145 void GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
146+ void GetAlpha2Josh(IssmDouble* palpha2,Gauss* gauss);
147 void GetAlpha2Sommers(IssmDouble* palpha2,Gauss* gauss);
148 void GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
149 void GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
150Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
151===================================================================
152--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 21550)
153+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 21551)
154@@ -822,6 +822,10 @@
155 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
156 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
157 break;
158+ case 9:
159+ iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
160+ iomodel->FetchDataToInput(elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum);
161+ break;
162 default:
163 _error_("not supported");
164 }
165Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
166===================================================================
167--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21550)
168+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21551)
169@@ -108,6 +108,7 @@
170 FrictionAsEnum,
171 FrictionCoefficientEnum,
172 FrictionCoefficientcoulombEnum,
173+ FrictionPressureAdjustedTemperatureEnum,
174 FrictionPEnum,
175 FrictionQEnum,
176 FrictionMEnum,
177@@ -273,6 +274,9 @@
178 MaterialsMantleDensityEnum,
179 MaterialsEarthDensityEnum,
180 MeshAverageVertexConnectivityEnum,
181+ MeshXEnum,
182+ MeshYEnum,
183+ MeshZEnum,
184 MeshElementsEnum,
185 MeshNumberofelementsEnum,
186 MeshNumberoflayersEnum,
187Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
188===================================================================
189--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21550)
190+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21551)
191@@ -114,6 +114,7 @@
192 case FrictionAsEnum : return "FrictionAs";
193 case FrictionCoefficientEnum : return "FrictionCoefficient";
194 case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb";
195+ case FrictionPressureAdjustedTemperatureEnum : return "FrictionPressureAdjustedTemperature";
196 case FrictionPEnum : return "FrictionP";
197 case FrictionQEnum : return "FrictionQ";
198 case FrictionMEnum : return "FrictionM";
199@@ -279,6 +280,9 @@
200 case MaterialsMantleDensityEnum : return "MaterialsMantleDensity";
201 case MaterialsEarthDensityEnum : return "MaterialsEarthDensity";
202 case MeshAverageVertexConnectivityEnum : return "MeshAverageVertexConnectivity";
203+ case MeshXEnum : return "MeshX";
204+ case MeshYEnum : return "MeshY";
205+ case MeshZEnum : return "MeshZ";
206 case MeshElementsEnum : return "MeshElements";
207 case MeshNumberofelementsEnum : return "MeshNumberofelements";
208 case MeshNumberoflayersEnum : return "MeshNumberoflayers";
209Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
210===================================================================
211--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21550)
212+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21551)
213@@ -114,6 +114,7 @@
214 else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
215 else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
216 else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
217+ else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
218 else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
219 else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
220 else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
221@@ -135,11 +136,11 @@
222 else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
223 else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
224 else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
225- else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
226 else stage=2;
227 }
228 if(stage==2){
229- if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
230+ if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
231+ else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
232 else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
233 else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum;
234 else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
235@@ -258,11 +259,11 @@
236 else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
237 else if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum;
238 else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
239- else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
240 else stage=3;
241 }
242 if(stage==3){
243- if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
244+ if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
245+ else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
246 else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
247 else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
248 else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
249@@ -285,6 +286,9 @@
250 else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum;
251 else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
252 else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
253+ else if (strcmp(name,"MeshX")==0) return MeshXEnum;
254+ else if (strcmp(name,"MeshY")==0) return MeshYEnum;
255+ else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
256 else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
257 else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum;
258 else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum;
259@@ -378,14 +382,14 @@
260 else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
261 else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
262 else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
263- else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
264+ else stage=4;
265+ }
266+ if(stage==4){
267+ if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
268 else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
269 else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
270 else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
271- else stage=4;
272- }
273- if(stage==4){
274- if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
275+ else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
276 else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
277 else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
278 else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
279@@ -501,14 +505,14 @@
280 else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
281 else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
282 else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
283- else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
284+ else stage=5;
285+ }
286+ if(stage==5){
287+ if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
288 else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
289 else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
290 else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
291- else stage=5;
292- }
293- if(stage==5){
294- if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
295+ else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
296 else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
297 else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
298 else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
299@@ -624,14 +628,14 @@
300 else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
301 else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
302 else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
303- else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
304+ else stage=6;
305+ }
306+ if(stage==6){
307+ if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
308 else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
309 else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
310 else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum;
311- else stage=6;
312- }
313- if(stage==6){
314- if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
315+ else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
316 else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
317 else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum;
318 else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum;
319@@ -747,14 +751,14 @@
320 else if (strcmp(name,"Sset")==0) return SsetEnum;
321 else if (strcmp(name,"Dense")==0) return DenseEnum;
322 else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
323- else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
324+ else stage=7;
325+ }
326+ if(stage==7){
327+ if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
328 else if (strcmp(name,"Seq")==0) return SeqEnum;
329 else if (strcmp(name,"Mpi")==0) return MpiEnum;
330 else if (strcmp(name,"Mumps")==0) return MumpsEnum;
331- else stage=7;
332- }
333- if(stage==7){
334- if (strcmp(name,"Gsl")==0) return GslEnum;
335+ else if (strcmp(name,"Gsl")==0) return GslEnum;
336 else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
337 else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
338 else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
339@@ -870,14 +874,14 @@
340 else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
341 else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
342 else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
343- else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
344+ else stage=8;
345+ }
346+ if(stage==8){
347+ if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
348 else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
349 else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
350 else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
351- else stage=8;
352- }
353- if(stage==8){
354- if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
355+ else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
356 else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
357 else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
358 else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
359@@ -993,14 +997,14 @@
360 else if (strcmp(name,"Melange")==0) return MelangeEnum;
361 else if (strcmp(name,"Water")==0) return WaterEnum;
362 else if (strcmp(name,"DataSet")==0) return DataSetEnum;
363- else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
364+ else stage=9;
365+ }
366+ if(stage==9){
367+ if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
368 else if (strcmp(name,"Loads")==0) return LoadsEnum;
369 else if (strcmp(name,"Materials")==0) return MaterialsEnum;
370 else if (strcmp(name,"Nodes")==0) return NodesEnum;
371- else stage=9;
372- }
373- if(stage==9){
374- if (strcmp(name,"Contours")==0) return ContoursEnum;
375+ else if (strcmp(name,"Contours")==0) return ContoursEnum;
376 else if (strcmp(name,"Parameters")==0) return ParametersEnum;
377 else if (strcmp(name,"Vertices")==0) return VerticesEnum;
378 else if (strcmp(name,"Results")==0) return ResultsEnum;
Note: See TracBrowser for help on using the repository browser.