source: issm/oecreview/Archive/23185-23389/ISSM-23241-23242.diff@ 23390

Last change on this file since 23390 was 23390, checked in by Mathieu Morlighem, 6 years ago

CHG: added Archive/23185-23389
File size: 28.5 KB
RevLine 
[23390]1Index: ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h
2===================================================================
3--- ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 23241)
4+++ ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 23242)
5@@ -621,7 +621,7 @@
6 } /*}}}*/
7
8 /*Specific instantiations for IssmDouble*: */
9-#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization.
10+#if defined(_HAVE_AD_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization.
11 template <> inline GenericExternalResult<IssmDouble*>::~GenericExternalResult(){ /*{{{*/
12 xDelete<char>(result_name);
13 xDelete<IssmDouble>(value);
14@@ -708,7 +708,7 @@
15 template <> inline Object* GenericExternalResult<Vector<IssmPDouble>*>::copy(void){ /*{{{*/
16 return new GenericExternalResult<Vector<IssmPDouble>*>(this->id,StringToEnumx(this->result_name),this->value,this->step,this->time);
17 } /*}}}*/
18-#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization.
19+#if defined(_HAVE_AD_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization.
20 template <> inline void GenericExternalResult<Vector<IssmPDouble>*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/
21
22 char *name = NULL;
23Index: ../trunk-jpl/src/c/classes/FemModel.cpp
24===================================================================
25--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23241)
26+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23242)
27@@ -880,7 +880,7 @@
28 results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentMemEnum, solution_memory));
29 results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops));
30
31- #ifdef _HAVE_ADOLC_
32+ #ifdef _HAVE_AD_
33 solution_time = profiler->TotalTime(ADCORE);
34 solution_flops = profiler->TotalFlops(ADCORE);
35 solution_memory = profiler->Memory(ADCORE);
36@@ -1844,7 +1844,7 @@
37 parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
38
39 if(isautodiff){
40- #ifdef _HAVE_ADOLC_
41+ #ifdef _HAVE_AD_
42 parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
43 parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
44 if(num_dependents){
45Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
46===================================================================
47--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23241)
48+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23242)
49@@ -3132,7 +3132,10 @@
50 fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
51 }
52
53- #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
54+ #if defined(_HAVE_AD_)
55+ /*we want to avoid the round operation at all cost. Not differentiable.*/
56+ _error_("not implemented yet");
57+ #else
58 dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
59 dMass = round(dMass * 100.0)/100.0;
60
61Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
62===================================================================
63--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23241)
64+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23242)
65@@ -52,6 +52,7 @@
66 AutodiffNumDependentsEnum,
67 AutodiffNumIndependentsEnum,
68 AutodiffObufsizeEnum,
69+ AutodiffTapeAllocEnum,
70 AutodiffTbufsizeEnum,
71 AutodiffXpEnum,
72 BalancethicknessStabilizationEnum,
73Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
74===================================================================
75--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23241)
76+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23242)
77@@ -60,6 +60,7 @@
78 case AutodiffNumDependentsEnum : return "AutodiffNumDependents";
79 case AutodiffNumIndependentsEnum : return "AutodiffNumIndependents";
80 case AutodiffObufsizeEnum : return "AutodiffObufsize";
81+ case AutodiffTapeAllocEnum : return "AutodiffTapeAlloc";
82 case AutodiffTbufsizeEnum : return "AutodiffTbufsize";
83 case AutodiffXpEnum : return "AutodiffXp";
84 case BalancethicknessStabilizationEnum : return "BalancethicknessStabilization";
85Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
86===================================================================
87--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23241)
88+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23242)
89@@ -60,6 +60,7 @@
90 else if (strcmp(name,"AutodiffNumDependents")==0) return AutodiffNumDependentsEnum;
91 else if (strcmp(name,"AutodiffNumIndependents")==0) return AutodiffNumIndependentsEnum;
92 else if (strcmp(name,"AutodiffObufsize")==0) return AutodiffObufsizeEnum;
93+ else if (strcmp(name,"AutodiffTapeAlloc")==0) return AutodiffTapeAllocEnum;
94 else if (strcmp(name,"AutodiffTbufsize")==0) return AutodiffTbufsizeEnum;
95 else if (strcmp(name,"AutodiffXp")==0) return AutodiffXpEnum;
96 else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum;
97@@ -135,11 +136,11 @@
98 else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum;
99 else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum;
100 else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum;
101- else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum;
102 else stage=2;
103 }
104 if(stage==2){
105- if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum;
106+ if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum;
107+ else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum;
108 else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum;
109 else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
110 else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
111@@ -258,11 +259,11 @@
112 else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum;
113 else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum;
114 else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
115- else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum;
116 else stage=3;
117 }
118 if(stage==3){
119- if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum;
120+ if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum;
121+ else if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum;
122 else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum;
123 else if (strcmp(name,"SealevelriseHoriz")==0) return SealevelriseHorizEnum;
124 else if (strcmp(name,"SealevelriseLoopIncrement")==0) return SealevelriseLoopIncrementEnum;
125@@ -381,11 +382,11 @@
126 else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
127 else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
128 else if (strcmp(name,"Velocity")==0) return VelocityEnum;
129- else if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
130 else stage=4;
131 }
132 if(stage==4){
133- if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
134+ if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
135+ else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
136 else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum;
137 else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
138 else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
139@@ -504,11 +505,11 @@
140 else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
141 else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
142 else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
143- else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
144 else stage=5;
145 }
146 if(stage==5){
147- if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
148+ if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
149+ else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
150 else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
151 else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
152 else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
153@@ -627,11 +628,11 @@
154 else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
155 else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
156 else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
157- else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
158 else stage=6;
159 }
160 if(stage==6){
161- if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
162+ if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
163+ else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
164 else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
165 else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
166 else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
167@@ -750,11 +751,11 @@
168 else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
169 else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
170 else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
171- else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
172 else stage=7;
173 }
174 if(stage==7){
175- if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
176+ if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
177+ else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
178 else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
179 else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
180 else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
181@@ -873,11 +874,11 @@
182 else if (strcmp(name,"LoveHi")==0) return LoveHiEnum;
183 else if (strcmp(name,"LoveHr")==0) return LoveHrEnum;
184 else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
185- else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
186 else stage=8;
187 }
188 if(stage==8){
189- if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
190+ if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
191+ else if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
192 else if (strcmp(name,"LoveKr")==0) return LoveKrEnum;
193 else if (strcmp(name,"LoveLi")==0) return LoveLiEnum;
194 else if (strcmp(name,"LoveLr")==0) return LoveLrEnum;
195@@ -996,11 +997,11 @@
196 else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum;
197 else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
198 else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
199- else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
200 else stage=9;
201 }
202 if(stage==9){
203- if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
204+ if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
205+ else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
206 else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
207 else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
208 else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
209@@ -1119,11 +1120,11 @@
210 else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
211 else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
212 else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
213- else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
214 else stage=10;
215 }
216 if(stage==10){
217- if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
218+ if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
219+ else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
220 else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
221 else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;
222 else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
223Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
224===================================================================
225--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23241)
226+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23242)
227@@ -373,7 +373,7 @@
228 /*Now, deal with toolkits options, which need to be put into the parameters dataset: */
229 ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);
230
231- #ifdef _HAVE_ADOLC_
232+ #ifdef _HAVE_AD_
233 if(VerboseMProcessor()) _printf0_(" starting autodiff parameters \n");
234 CreateParametersAutodiff(parameters,iomodel);
235 if(VerboseMProcessor()) _printf0_(" ending autodiff parameters \n");
236Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp
237===================================================================
238--- ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23241)
239+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23242)
240@@ -72,8 +72,36 @@
241 /*Free ressources: */
242 xDelete<char>(options);
243 /*}}}*/
244+ #elif _HAVE_CODIPACK_
245+ //fprintf(stderr, "*** Codipack CreateParametersAutodiff()\n");
246+ /*initialize a placeholder to store solver pointers: {{{*/
247+ /*Solver pointers depend on what type of solver we are implementing: */
248+ options=OptionsFromAnalysis(parameters,DefaultAnalysisEnum); //options database is not filled in yet, use default.
249+ ToolkitOptions::Init(options);
250
251+ switch(IssmSolverTypeFromToolkitOptions()){
252+ case MumpsEnum:{
253+ #ifndef _HAVE_MUMPS_
254+ _error_("CoDiPack: requesting mumps solver without MUMPS being compiled in!");
255+ #endif
256+ break;
257+ }
258+ case GslEnum: {
259+ #ifndef _HAVE_GSL_
260+ _error_("CoDiPack: requesting GSL solver without GSL being compiled in!");
261+ #endif
262+ break;
263+ }
264+ default:
265+ _error_("solver type not supported yet!");
266+ }
267+ /*Free ressources: */
268+ xDelete<char>(options);
269+ #endif
270+ #if defined(_HAVE_AD_)
271+
272 if(isautodiff){
273+#if _HAVE_ADOLC_
274 /*Copy some parameters from IoModel to parameters dataset: {{{*/
275 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
276 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
277@@ -82,11 +110,16 @@
278 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
279 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
280 /*}}}*/
281+#endif
282 /*retrieve driver: {{{*/
283 iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");
284 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum));
285
286 if(strcmp(autodiff_driver,"fos_forward")==0){
287+#if _HAVE_CODIPACK_
288+ // FIXME codi support Foward Mode (scalar)
289+ _error_("Foward Mode (scalar) not supported yet!");
290+#endif
291 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_forward_index",AutodiffFosForwardIndexEnum));
292 }
293 else if(strcmp(autodiff_driver,"fos_reverse")==0){
294@@ -93,6 +126,10 @@
295 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_reverse_index",AutodiffFosReverseIndexEnum));
296 }
297 else if(strcmp(autodiff_driver,"fov_forward")==0){
298+#if _HAVE_CODIPACK_
299+ // FIXME codi support Foward Mode (vector)
300+ _error_("Foward Mode (vector) not supported yet!");
301+#endif
302 /*Retrieve list of indices: */
303 iomodel->FetchData(&indices,&num_indices,&dummy,"md.autodiff.fov_forward_indices");
304 parameters->AddObject(new IntMatParam(AutodiffFovForwardIndicesEnum,indices,num_indices,1));
305Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
306===================================================================
307--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23241)
308+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23242)
309@@ -61,7 +61,7 @@
310 //check to see if the top grid cell structure length (dzTop) goes evenly
311 //into specified top structure depth (zTop). Also make sure top grid cell
312 //structure length (dzTop) is greater than 5 cm
313- #ifndef _HAVE_ADOLC_ //avoid the round operation check!
314+ #ifndef _HAVE_AD_ //avoid the round operation check!
315 if (dgpTop != round(dgpTop)){
316 _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop");
317 }
318@@ -424,11 +424,11 @@
319
320 // spectral range:
321 // 0.3 - 0.8um
322- IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5));
323+ IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5));
324 // 0.8 - 1.5um
325- IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5));
326+ IssmDouble a1 = max(0, 0.95 - 15.4 *pow(gsz,0.5));
327 // 1.5 - 2.8um
328- IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
329+ IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
330
331 // broadband surface albedo
332 a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
333@@ -666,7 +666,7 @@
334 // determine minimum acceptable delta t (diffusion number > 1/2) [s]
335 // NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability
336 dt=1e12;
337- for(int i=0;i<m;i++)dt = fmin(dt,CI * pow(dz[i],2) * d[i] / (3 * K[i]) * thermo_scaling);
338+ for(int i=0;i<m;i++)dt = min(dt,CI * pow(dz[i],2) * d[i] / (3 * K[i]) * thermo_scaling);
339
340 // smallest possible even integer of 60 min where diffusion number > 1/2
341 // must go evenly into one hour or the data frequency if it is smaller
342@@ -751,7 +751,7 @@
343 // calculated. The estimated enegy balance & melt are significanly
344 // less when Ts is taken as the mean of the x top grid cells.
345 Ts = (T[0] + T[1])/2.0;
346- Ts = fmin(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC)
347+ Ts = min(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC)
348
349 //TURBULENT HEAT FLUX
350
351@@ -766,7 +766,7 @@
352 // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
353
354 // do not allow Ri to exceed 0.19
355- Ri = fmin(Ri, 0.19);
356+ Ri = min(Ri, 0.19);
357
358 // calculate momentum 'coefM' stability factor
359 if (Ri > 0.0+Ttol){
360@@ -946,11 +946,11 @@
361
362 // spectral albedos:
363 // 0.3 - 0.8um
364- IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
365+ IssmDouble a0 = min(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
366 // 0.8 - 1.5um
367- IssmDouble a1 = fmax(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
368+ IssmDouble a1 = max(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
369 // 1.5 - 2.8um
370- IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
371+ IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
372
373 // separate net shortwave radiative flux into spectral ranges
374 IssmDouble swfS[3];
375@@ -1203,7 +1203,7 @@
376
377 mass_diff = mass - massinit - P;
378
379- #ifndef _HAVE_ADOLC_ //avoid round operation. only check in forward mode.
380+ #ifndef _HAVE_AD_ //avoid round operation. only check in forward mode.
381 mass_diff = round(mass_diff * 100.0)/100.0;
382 if (mass_diff > 0) _error_("mass not conserved in accumulation function");
383 #endif
384@@ -1322,11 +1322,11 @@
385
386 // calculate temperature excess above 0 deg C
387 exsT=xNewZeroInit<IssmDouble>(n);
388- for(int i=0;i<n;i++) exsT[i]= fmax(0.0, T[i] - CtoK); // [K] to [degC]
389+ for(int i=0;i<n;i++) exsT[i]= max(0.0, T[i] - CtoK); // [K] to [degC]
390
391 // new grid point center temperature, T [K]
392 // for(int i=0;i<n;i++) T[i]-=exsT[i];
393- for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
394+ for(int i=0;i<n;i++) T[i]=min(T[i],CtoK);
395
396 // specify irreducible water content saturation [fraction]
397 const IssmDouble Swi = 0.07; // assumed constant after Colbeck, 1974
398@@ -1336,10 +1336,10 @@
399 if (cellsum(W,n) >0.0+Wtol){
400 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" pore water refreeze\n");
401 // calculate maximum freeze amount, maxF [kg]
402- for(int i=0;i<n;i++) maxF[i] = fmax(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
403+ for(int i=0;i<n;i++) maxF[i] = max(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
404
405 // freeze pore water and change snow/ice properties
406- for(int i=0;i<n;i++) dW[i] = fmin(maxF[i], W[i]); // freeze mass [kg]
407+ for(int i=0;i<n;i++) dW[i] = min(maxF[i], W[i]); // freeze mass [kg]
408 for(int i=0;i<n;i++) W[i] -= dW[i]; // pore water mass [kg]
409 for(int i=0;i<n;i++) m[i] += dW[i]; // new mass [kg]
410 for(int i=0;i<n;i++) d[i] = m[i] / dz[i]; // density [kg m-3]
411@@ -1355,7 +1355,7 @@
412 exsW=xNew<IssmDouble>(n);
413 for(int i=0;i<n;i++){
414 Wi= (dIce - d[i]) * Swi * (m[i] / d[i]); // irreducible water content [kg]
415- exsW[i] = fmax(0.0, W[i] - Wi); // water "squeezed" from snow [kg]
416+ exsW[i] = max(0.0, W[i] - Wi); // water "squeezed" from snow [kg]
417 }
418
419 //// MELT, PERCOLATION AND REFREEZE
420@@ -1367,7 +1367,7 @@
421 // if so redistribute temperature to lower cells (temperature surplus)
422 // (maximum T of snow before entire grid cell melts is a constant
423 // LF/CI = 159.1342)
424- surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- LF/CI);
425+ surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = max(0.0, exsT[i]- LF/CI);
426
427 if (cellsum(surpT,n) > 0.0 + Ttol ){
428 // _printf_("T Surplus");
429@@ -1379,10 +1379,10 @@
430 // use surplus energy to increase the temperature of lower cell
431 T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
432
433- exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1];
434- T[i+1] = fmin(CtoK, T[i+1]);
435+ exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
436+ T[i+1] = min(CtoK, T[i+1]);
437
438- surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI);
439+ surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
440 surpE[i+1] = surpT[i+1] * CI * m[i+1];
441
442 // adjust current cell properties (again 159.1342 is the max T)
443@@ -1393,11 +1393,11 @@
444 }
445
446 // convert temperature excess to melt [kg]
447- for(int i=0;i<n;i++) M[i] = fmin(exsT[i] * d[i] * dz[i] * CI / LF, m[i]); // melt
448+ for(int i=0;i<n;i++) M[i] = min(exsT[i] * d[i] * dz[i] * CI / LF, m[i]); // melt
449 sumM = cellsum(M,n); // total melt [kg]
450
451 // calculate maximum refreeze amount, maxF [kg]
452- for(int i=0;i<n;i++)maxF[i] = fmax(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
453+ for(int i=0;i<n;i++)maxF[i] = max(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
454
455 // initialize refreeze, runoff, flxDn and dW vectors [kg]
456 IssmDouble* F = xNewZeroInit<IssmDouble>(n);
457@@ -1434,8 +1434,8 @@
458
459 m[i] = m[i] - M[i]; // mass after melt
460 Wi = (dIce-d[i]) * Swi * (m[i]/d[i]); // irreducible water
461- dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]); // change in pore water
462- R[i] = fmax(0.0, inM - dW[i]); // runoff
463+ dW[i] = max(min(inM, Wi - W[i]),-1*W[i]); // change in pore water
464+ R[i] = max(0.0, inM - dW[i]); // runoff
465 F[i] = 0.0;
466 }
467 // check if no energy to refreeze meltwater
468@@ -1446,8 +1446,8 @@
469
470 m[i] = m[i] - M[i]; // mass after melt
471 Wi = (dIce-d[i]) * Swi * (m[i]/d[i]); // irreducible water
472- dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]); // change in pore water
473- flxDn[i+1] = fmax(0.0, inM - dW[i]); // meltwater out
474+ dW[i] = max(min(inM, Wi - W[i]),-1*W[i]); // change in pore water
475+ flxDn[i+1] = max(0.0, inM - dW[i]); // meltwater out
476 R[i] = 0.0;
477 F[i] = 0.0; // no freeze
478 }
479@@ -1459,13 +1459,13 @@
480 m[i] = m[i] - M[i];
481 IssmDouble dz_0 = m[i]/d[i];
482 IssmDouble dMax = (dIce - d[i])*dz_0; // d max = dIce
483- IssmDouble F1 = fmin(fmin(inM,dMax),maxF[i]); // maximum refreeze
484+ IssmDouble F1 = min(min(inM,dMax),maxF[i]); // maximum refreeze
485 m[i] = m[i] + F1; // mass after refreeze
486 d[i] = m[i]/dz_0;
487
488 //-----------------------pore water-----------------------------
489 Wi = (dIce-d[i])* Swi * dz_0; // irreducible water
490- dW[i] = fmin(inM - F1, Wi-W[i]); // change in pore water
491+ dW[i] = min(inM - F1, Wi-W[i]); // change in pore water
492 if (dW[i] < 0.0-Wtol && -1*dW[i]>W[i]-Wtol ){
493 dW[i]= -1*W[i];
494 }
495@@ -1473,8 +1473,8 @@
496
497 if (dW[i] < 0.0-Wtol){ // excess pore water
498 dMax = (dIce - d[i])*dz_0; // maximum refreeze
499- IssmDouble maxF2 = fmin(dMax, maxF[i]-F1); // maximum refreeze
500- F2 = fmin(-1*dW[i], maxF2); // pore water refreeze
501+ IssmDouble maxF2 = min(dMax, maxF[i]-F1); // maximum refreeze
502+ F2 = min(-1*dW[i], maxF2); // pore water refreeze
503 m[i] = m[i] + F2; // mass after refreeze
504 d[i] = m[i]/dz_0;
505 dW[i] = dW[i] - F2;
506@@ -1752,7 +1752,7 @@
507 for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
508
509 /*only in forward mode! avoid round in AD mode as it is not differentiable: */
510- #ifndef _HAVE_ADOLC_
511+ #ifndef _HAVE_AD_
512 dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
513 dE = round(sumE0 - sumE1 - sumER + addE);
514 if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
515@@ -1914,8 +1914,8 @@
516 H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
517 c0arth = 0.07 * H;
518 c1arth = 0.03 * H;
519- M0 = fmax(1.435 - (0.151 * log(C)),0.25);
520- M1 = fmax(2.366 - (0.293 * log(C)),0.25);
521+ M0 = max(1.435 - (0.151 * log(C)),0.25);
522+ M1 = max(2.366 - (0.293 * log(C)),0.25);
523 c0 = M0*c0arth;
524 c1 = M1*c1arth;
525 break;
526@@ -1925,8 +1925,8 @@
527 H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
528 c0arth = 0.07 * H;
529 c1arth = 0.03 * H;
530- M0 = fmax(1.042 - (0.0916 * log(C)),0.25);
531- M1 = fmax(1.734 - (0.2039 * log(C)),0.25);
532+ M0 = max(1.042 - (0.0916 * log(C)),0.25);
533+ M1 = max(1.734 - (0.2039 * log(C)),0.25);
534 c0 = M0*c0arth;
535 c1 = M1*c1arth;
536 break;
537@@ -2020,7 +2020,7 @@
538 // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
539
540 // do not allow Ri to exceed 0.19
541- Ri = fmin(Ri, 0.19);
542+ Ri = min(Ri, 0.19);
543
544 // calculate momentum 'coefM' stability factor
545 if (Ri > 0.0+Ttol){
546Index: ../trunk-jpl/src/c/cores/controlvalidation_core.cpp
547===================================================================
548--- ../trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 23241)
549+++ ../trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 23242)
550@@ -87,7 +87,7 @@
551 }
552
553 /*output*/
554- #ifdef _HAVE_ADOLC_
555+ #ifdef _HAVE_AD_
556 IssmPDouble* J_passive=xNew<IssmPDouble>(2*num);
557 for(int i=0;i<2*num;i++) J_passive[i]=reCast<IssmPDouble>(output[i]);
558 femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0));
559Index: ../trunk-jpl/src/c/cores/control_core.cpp
560===================================================================
561--- ../trunk-jpl/src/c/cores/control_core.cpp (revision 23241)
562+++ ../trunk-jpl/src/c/cores/control_core.cpp (revision 23242)
563@@ -96,7 +96,7 @@
564 if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run!
565 femmodel->OutputControlsx(&femmodel->results);
566
567- #ifdef _HAVE_ADOLC_
568+ #ifdef _HAVE_AD_
569 IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
570 for(int i=0;i<nsteps;i++) J_passive[i]=reCast<IssmPDouble>(J[i]);
571 femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,0,0));
Note: See TracBrowser for help on using the repository browser.