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

  • ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h

     
    621621} /*}}}*/
    622622
    623623        /*Specific instantiations for IssmDouble*: */
    624 #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.
     624#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.
    625625template <> inline GenericExternalResult<IssmDouble*>::~GenericExternalResult(){ /*{{{*/
    626626        xDelete<char>(result_name);
    627627        xDelete<IssmDouble>(value);
    … …  
    708708        template <> inline Object* GenericExternalResult<Vector<IssmPDouble>*>::copy(void){ /*{{{*/
    709709                return new GenericExternalResult<Vector<IssmPDouble>*>(this->id,StringToEnumx(this->result_name),this->value,this->step,this->time);
    710710        } /*}}}*/
    711 #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.
     711#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.
    712712        template <> inline void GenericExternalResult<Vector<IssmPDouble>*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/
    713713
    714714                char *name   = NULL;

  • ../trunk-jpl/src/c/classes/FemModel.cpp

     
    880880                results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentMemEnum,  solution_memory));
    881881                results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops));
    882882
    883                 #ifdef _HAVE_ADOLC_
     883                #ifdef _HAVE_AD_
    884884                solution_time   = profiler->TotalTime(ADCORE);
    885885                solution_flops  = profiler->TotalFlops(ADCORE);
    886886                solution_memory = profiler->Memory(ADCORE);
    … …  
    18441844        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    18451845
    18461846        if(isautodiff){
    1847                 #ifdef _HAVE_ADOLC_
     1847                #ifdef _HAVE_AD_
    18481848                parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    18491849                parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
    18501850                if(num_dependents){

  • ../trunk-jpl/src/c/classes/Elements/Element.cpp

     
    31323132                        fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
    31333133                }
    31343134
    3135                 #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
     3135                #if defined(_HAVE_AD_)
     3136                /*we want to avoid the round operation at all cost. Not differentiable.*/
     3137                _error_("not implemented yet");
     3138                #else
    31363139                dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
    31373140                dMass = round(dMass * 100.0)/100.0;
    31383141

  • ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

     
    5252        AutodiffNumDependentsEnum,
    5353        AutodiffNumIndependentsEnum,
    5454        AutodiffObufsizeEnum,
     55        AutodiffTapeAllocEnum,
    5556        AutodiffTbufsizeEnum,
    5657        AutodiffXpEnum,
    5758        BalancethicknessStabilizationEnum,

  • ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

     
    6060                case AutodiffNumDependentsEnum : return "AutodiffNumDependents";
    6161                case AutodiffNumIndependentsEnum : return "AutodiffNumIndependents";
    6262                case AutodiffObufsizeEnum : return "AutodiffObufsize";
     63                case AutodiffTapeAllocEnum : return "AutodiffTapeAlloc";
    6364                case AutodiffTbufsizeEnum : return "AutodiffTbufsize";
    6465                case AutodiffXpEnum : return "AutodiffXp";
    6566                case BalancethicknessStabilizationEnum : return "BalancethicknessStabilization";

  • ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

     
    6060              else if (strcmp(name,"AutodiffNumDependents")==0) return AutodiffNumDependentsEnum;
    6161              else if (strcmp(name,"AutodiffNumIndependents")==0) return AutodiffNumIndependentsEnum;
    6262              else if (strcmp(name,"AutodiffObufsize")==0) return AutodiffObufsizeEnum;
     63              else if (strcmp(name,"AutodiffTapeAlloc")==0) return AutodiffTapeAllocEnum;
    6364              else if (strcmp(name,"AutodiffTbufsize")==0) return AutodiffTbufsizeEnum;
    6465              else if (strcmp(name,"AutodiffXp")==0) return AutodiffXpEnum;
    6566              else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum;
    … …  
    135136              else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum;
    136137              else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum;
    137138              else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum;
    138               else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum;
     142              if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum;
     143              else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum;
    143144              else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum;
    144145              else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
    145146              else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
    … …  
    258259              else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum;
    259260              else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum;
    260261              else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
    261               else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum;
     265              if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum;
     266              else if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum;
    266267              else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum;
    267268              else if (strcmp(name,"SealevelriseHoriz")==0) return SealevelriseHorizEnum;
    268269              else if (strcmp(name,"SealevelriseLoopIncrement")==0) return SealevelriseLoopIncrementEnum;
    … …  
    381382              else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
    382383              else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
    383384              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
    384               else if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
     388              if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
     389              else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
    389390              else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum;
    390391              else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
    391392              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    … …  
    504505              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    505506              else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
    506507              else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
    507               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
     511              if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
     512              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    512513              else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
    513514              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
    514515              else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
    … …  
    627628              else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
    628629              else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
    629630              else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
    630               else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
     634              if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
     635              else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
    635636              else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
    636637              else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
    637638              else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
    … …  
    750751              else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
    751752              else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
    752753              else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
    753               else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
     757              if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
     758              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
    758759              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    759760              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    760761              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    … …  
    873874              else if (strcmp(name,"LoveHi")==0) return LoveHiEnum;
    874875              else if (strcmp(name,"LoveHr")==0) return LoveHrEnum;
    875876              else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    876               else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
     880              if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
     881              else if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
    881882              else if (strcmp(name,"LoveKr")==0) return LoveKrEnum;
    882883              else if (strcmp(name,"LoveLi")==0) return LoveLiEnum;
    883884              else if (strcmp(name,"LoveLr")==0) return LoveLrEnum;
    … …  
    996997              else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum;
    997998              else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
    998999              else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
    999               else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
     1003              if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
     1004              else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
    10041005              else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
    10051006              else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
    10061007              else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
    … …  
    11191120              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
    11201121              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
    11211122              else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
    1122               else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
     1126              if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
     1127              else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
    11271128              else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    11281129              else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;
    11291130              else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;

  • ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

     
    373373        /*Now, deal with toolkits options, which need to be put into the parameters dataset: */
    374374        ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);
    375375
    376         #ifdef _HAVE_ADOLC_
     376        #ifdef _HAVE_AD_
    377377        if(VerboseMProcessor()) _printf0_("   starting autodiff parameters \n");
    378378        CreateParametersAutodiff(parameters,iomodel);
    379379        if(VerboseMProcessor()) _printf0_("   ending autodiff parameters \n");

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp

     
    7272                /*Free ressources: */
    7373                xDelete<char>(options);
    7474                /*}}}*/
     75                #elif _HAVE_CODIPACK_
     76                //fprintf(stderr, "*** Codipack CreateParametersAutodiff()\n");
     77                /*initialize a placeholder to store solver pointers: {{{*/
     78                /*Solver pointers depend on what type of solver we are implementing: */
     79                options=OptionsFromAnalysis(parameters,DefaultAnalysisEnum); //options database is not filled in yet, use default.
     80                ToolkitOptions::Init(options);
    7581
     82                switch(IssmSolverTypeFromToolkitOptions()){
     83                        case MumpsEnum:{
     84                                #ifndef _HAVE_MUMPS_
     85                                _error_("CoDiPack: requesting mumps solver without MUMPS being compiled in!");
     86                                #endif
     87                                break;
     88                                }
     89                        case GslEnum: {
     90                                #ifndef _HAVE_GSL_
     91                                _error_("CoDiPack: requesting GSL solver without GSL being compiled in!");
     92                                #endif
     93                                break;
     94                                }
     95                        default:
     96                                                        _error_("solver type not supported yet!");
     97                }
     98                /*Free ressources: */
     99                xDelete<char>(options);
     100                #endif
     101                #if defined(_HAVE_AD_)
     102
    76103        if(isautodiff){
     104#if _HAVE_ADOLC_
    77105                /*Copy some parameters from IoModel to parameters dataset: {{{*/
    78106                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
    79107                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
    … …  
    82110                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
    83111                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
    84112                /*}}}*/
     113#endif
    85114                /*retrieve driver: {{{*/
    86115                iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");
    87116                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum));
    88117
    89118                if(strcmp(autodiff_driver,"fos_forward")==0){
     119#if _HAVE_CODIPACK_
     120                        // FIXME codi support Foward Mode (scalar)
     121                        _error_("Foward Mode (scalar) not supported yet!");
     122#endif
    90123                        parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_forward_index",AutodiffFosForwardIndexEnum));
    91124                }
    92125                else if(strcmp(autodiff_driver,"fos_reverse")==0){
    … …  
    93126                        parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_reverse_index",AutodiffFosReverseIndexEnum));
    94127                }
    95128                else if(strcmp(autodiff_driver,"fov_forward")==0){
     129#if _HAVE_CODIPACK_
     130                        // FIXME codi support Foward Mode (vector)
     131                        _error_("Foward Mode (vector) not supported yet!");
     132#endif
    96133                        /*Retrieve list of indices: */
    97134                        iomodel->FetchData(&indices,&num_indices,&dummy,"md.autodiff.fov_forward_indices");
    98135                        parameters->AddObject(new IntMatParam(AutodiffFovForwardIndicesEnum,indices,num_indices,1));

  • ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

     
    6161        //check to see if the top grid cell structure length (dzTop) goes evenly
    6262        //into specified top structure depth (zTop). Also make sure top grid cell
    6363        //structure length (dzTop) is greater than 5 cm
    64         #ifndef _HAVE_ADOLC_  //avoid the round operation check!
     64        #ifndef _HAVE_AD_  //avoid the round operation check!
    6565        if (dgpTop != round(dgpTop)){
    6666                _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop");
    6767        }
    … …  
    424424
    425425                        // spectral range:
    426426                        // 0.3 - 0.8um
    427                         IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5));
     427                        IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5));
    428428                        // 0.8 - 1.5um
    429                         IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5));
     429                        IssmDouble a1 = max(0, 0.95 - 15.4 *pow(gsz,0.5));
    430430                        // 1.5 - 2.8um
    431                         IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
     431                        IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
    432432
    433433                        // broadband surface albedo
    434434                        a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
    … …  
    666666        // determine minimum acceptable delta t (diffusion number > 1/2) [s]
    667667        // NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability
    668668        dt=1e12;
    669         for(int i=0;i<m;i++)dt = fmin(dt,CI * pow(dz[i],2) * d[i]  / (3 * K[i]) * thermo_scaling);
     669        for(int i=0;i<m;i++)dt = min(dt,CI * pow(dz[i],2) * d[i]  / (3 * K[i]) * thermo_scaling);
    670670
    671671        // smallest possible even integer of 60 min where diffusion number > 1/2
    672672        // must go evenly into one hour or the data frequency if it is smaller
    … …  
    751751                // calculated.  The estimated enegy balance & melt are significanly
    752752                // less when Ts is taken as the mean of the x top grid cells.
    753753                Ts = (T[0] + T[1])/2.0;
    754                 Ts = fmin(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
     754                Ts = min(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
    755755
    756756                //TURBULENT HEAT FLUX
    757757
    … …  
    766766                // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
    767767
    768768                // do not allow Ri to exceed 0.19
    769                 Ri = fmin(Ri, 0.19);
     769                Ri = min(Ri, 0.19);
    770770
    771771                // calculate momentum 'coefM' stability factor
    772772                if (Ri > 0.0+Ttol){
    … …  
    946946
    947947                        // spectral albedos:
    948948                        // 0.3 - 0.8um
    949                         IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
     949                        IssmDouble a0 = min(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
    950950                        // 0.8 - 1.5um
    951                         IssmDouble a1 = fmax(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
     951                        IssmDouble a1 = max(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
    952952                        // 1.5 - 2.8um
    953                         IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
     953                        IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
    954954
    955955                        // separate net shortwave radiative flux into spectral ranges
    956956                        IssmDouble swfS[3];
    … …  
    12031203
    12041204                mass_diff = mass - massinit - P;
    12051205
    1206                 #ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
     1206                #ifndef _HAVE_AD_  //avoid round operation. only check in forward mode.
    12071207                mass_diff = round(mass_diff * 100.0)/100.0;
    12081208                if (mass_diff > 0) _error_("mass not conserved in accumulation function");
    12091209                #endif
    … …  
    13221322
    13231323        // calculate temperature excess above 0 deg C
    13241324        exsT=xNewZeroInit<IssmDouble>(n);
    1325         for(int i=0;i<n;i++) exsT[i]= fmax(0.0, T[i] - CtoK);        // [K] to [degC]
     1325        for(int i=0;i<n;i++) exsT[i]= max(0.0, T[i] - CtoK);        // [K] to [degC]
    13261326
    13271327        // new grid point center temperature, T [K]
    13281328        // for(int i=0;i<n;i++) T[i]-=exsT[i];
    1329         for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
     1329        for(int i=0;i<n;i++) T[i]=min(T[i],CtoK);
    13301330
    13311331        // specify irreducible water content saturation [fraction]
    13321332        const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
    … …  
    13361336        if (cellsum(W,n) >0.0+Wtol){
    13371337                if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("      pore water refreeze\n");
    13381338                // calculate maximum freeze amount, maxF [kg]
    1339                 for(int i=0;i<n;i++) maxF[i] = fmax(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
     1339                for(int i=0;i<n;i++) maxF[i] = max(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
    13401340
    13411341                // freeze pore water and change snow/ice properties
    1342                 for(int i=0;i<n;i++) dW[i] = fmin(maxF[i], W[i]);    // freeze mass [kg]   
     1342                for(int i=0;i<n;i++) dW[i] = min(maxF[i], W[i]);    // freeze mass [kg]   
    13431343                for(int i=0;i<n;i++) W[i] -= dW[i];                                            // pore water mass [kg]
    13441344                for(int i=0;i<n;i++) m[i] += dW[i];                                            // new mass [kg]
    13451345                for(int i=0;i<n;i++) d[i] = m[i] / dz[i];                                    // density [kg m-3]   
    … …  
    13551355        exsW=xNew<IssmDouble>(n);
    13561356        for(int i=0;i<n;i++){
    13571357                Wi= (dIce - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
    1358                 exsW[i] = fmax(0.0, W[i] - Wi);                  // water "squeezed" from snow [kg]
     1358                exsW[i] = max(0.0, W[i] - Wi);                  // water "squeezed" from snow [kg]
    13591359        }
    13601360
    13611361        //// MELT, PERCOLATION AND REFREEZE
    … …  
    13671367                // if so redistribute temperature to lower cells (temperature surplus)
    13681368                // (maximum T of snow before entire grid cell melts is a constant
    13691369                // LF/CI = 159.1342)
    1370                 surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- LF/CI);
     1370                surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = max(0.0, exsT[i]- LF/CI);
    13711371
    13721372                if (cellsum(surpT,n) > 0.0 + Ttol ){
    13731373                        // _printf_("T Surplus");
    … …  
    13791379                                // use surplus energy to increase the temperature of lower cell
    13801380                                T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
    13811381
    1382                                 exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1];
    1383                                 T[i+1] = fmin(CtoK, T[i+1]);
     1382                                exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
     1383                                T[i+1] = min(CtoK, T[i+1]);
    13841384
    1385                                 surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI);
     1385                                surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
    13861386                                surpE[i+1] = surpT[i+1] * CI * m[i+1];
    13871387
    13881388                                // adjust current cell properties (again 159.1342 is the max T)
    … …  
    13931393                }
    13941394
    13951395                // convert temperature excess to melt [kg]
    1396                 for(int i=0;i<n;i++) M[i] = fmin(exsT[i] * d[i] * dz[i] * CI / LF, m[i]);  // melt
     1396                for(int i=0;i<n;i++) M[i] = min(exsT[i] * d[i] * dz[i] * CI / LF, m[i]);  // melt
    13971397                sumM = cellsum(M,n);                                                       // total melt [kg]
    13981398
    13991399                // calculate maximum refreeze amount, maxF [kg]
    1400                 for(int i=0;i<n;i++)maxF[i] = fmax(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
     1400                for(int i=0;i<n;i++)maxF[i] = max(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
    14011401
    14021402                // initialize refreeze, runoff, flxDn and dW vectors [kg]
    14031403                IssmDouble* F = xNewZeroInit<IssmDouble>(n);
    … …  
    14341434
    14351435                                m[i] = m[i] - M[i];                     // mass after melt
    14361436                                Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water
    1437                                 dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]);            // change in pore water
    1438                                 R[i] = fmax(0.0, inM - dW[i]);             // runoff
     1437                                dW[i] = max(min(inM, Wi - W[i]),-1*W[i]);            // change in pore water
     1438                                R[i] = max(0.0, inM - dW[i]);             // runoff
    14391439                                F[i] = 0.0;
    14401440                        }
    14411441                        // check if no energy to refreeze meltwater
    … …  
    14461446
    14471447                                m[i] = m[i] - M[i];                     // mass after melt
    14481448                                Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water
    1449                                 dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]);              // change in pore water
    1450                                 flxDn[i+1] = fmax(0.0, inM - dW[i]);         // meltwater out
     1449                                dW[i] = max(min(inM, Wi - W[i]),-1*W[i]);              // change in pore water
     1450                                flxDn[i+1] = max(0.0, inM - dW[i]);         // meltwater out
    14511451                                R[i] = 0.0;
    14521452                                F[i] = 0.0;                               // no freeze
    14531453                        }
    … …  
    14591459                                m[i] = m[i] - M[i];
    14601460                                IssmDouble dz_0 = m[i]/d[i];         
    14611461                                IssmDouble dMax = (dIce - d[i])*dz_0;              // d max = dIce
    1462                                 IssmDouble F1 = fmin(fmin(inM,dMax),maxF[i]);         // maximum refreeze               
     1462                                IssmDouble F1 = min(min(inM,dMax),maxF[i]);         // maximum refreeze               
    14631463                                m[i] = m[i] + F1;                       // mass after refreeze
    14641464                                d[i] = m[i]/dz_0;
    14651465
    14661466                                //-----------------------pore water-----------------------------
    14671467                                Wi = (dIce-d[i])* Swi * dz_0;            // irreducible water
    1468                                 dW[i] = fmin(inM - F1, Wi-W[i]);         // change in pore water
     1468                                dW[i] = min(inM - F1, Wi-W[i]);         // change in pore water
    14691469                                if (dW[i] < 0.0-Wtol && -1*dW[i]>W[i]-Wtol ){
    14701470                                        dW[i]= -1*W[i];
    14711471                                }
    … …  
    14731473
    14741474                                if (dW[i] < 0.0-Wtol){                         // excess pore water
    14751475                                        dMax = (dIce - d[i])*dz_0;          // maximum refreeze                                             
    1476                                         IssmDouble maxF2 = fmin(dMax, maxF[i]-F1);      // maximum refreeze
    1477                                         F2 = fmin(-1*dW[i], maxF2);            // pore water refreeze
     1476                                        IssmDouble maxF2 = min(dMax, maxF[i]-F1);      // maximum refreeze
     1477                                        F2 = min(-1*dW[i], maxF2);            // pore water refreeze
    14781478                                        m[i] = m[i] + F2;                   // mass after refreeze
    14791479                                        d[i] = m[i]/dz_0;
    14801480                                        dW[i] = dW[i] - F2;
    … …  
    17521752        for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
    17531753
    17541754        /*only in forward mode! avoid round in AD mode as it is not differentiable: */
    1755         #ifndef _HAVE_ADOLC_
     1755        #ifndef _HAVE_AD_
    17561756        dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
    17571757        dE = round(sumE0 - sumE1 - sumER +  addE);
    17581758        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
    … …  
    19141914                                H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
    19151915                                c0arth = 0.07 * H;
    19161916                                c1arth = 0.03 * H;
    1917                                 M0 = fmax(1.435 - (0.151 * log(C)),0.25);
    1918                                 M1 = fmax(2.366 - (0.293 * log(C)),0.25);
     1917                                M0 = max(1.435 - (0.151 * log(C)),0.25);
     1918                                M1 = max(2.366 - (0.293 * log(C)),0.25);
    19191919                                c0 = M0*c0arth;
    19201920                                c1 = M1*c1arth;
    19211921                                break;
    … …  
    19251925                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
    19261926                                c0arth = 0.07 * H;
    19271927                                c1arth = 0.03 * H;
    1928                                 M0 = fmax(1.042 - (0.0916 * log(C)),0.25);
    1929                                 M1 = fmax(1.734 - (0.2039 * log(C)),0.25);
     1928                                M0 = max(1.042 - (0.0916 * log(C)),0.25);
     1929                                M1 = max(1.734 - (0.2039 * log(C)),0.25);
    19301930                                c0 = M0*c0arth;
    19311931                                c1 = M1*c1arth;
    19321932                                break;
    … …  
    20202020        // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
    20212021
    20222022        // do not allow Ri to exceed 0.19
    2023         Ri = fmin(Ri, 0.19);
     2023        Ri = min(Ri, 0.19);
    20242024
    20252025        // calculate momentum 'coefM' stability factor
    20262026        if (Ri > 0.0+Ttol){

  • ../trunk-jpl/src/c/cores/controlvalidation_core.cpp

     
    8787        }
    8888
    8989        /*output*/
    90         #ifdef _HAVE_ADOLC_
     90        #ifdef _HAVE_AD_
    9191        IssmPDouble* J_passive=xNew<IssmPDouble>(2*num);
    9292        for(int i=0;i<2*num;i++) J_passive[i]=reCast<IssmPDouble>(output[i]);
    9393        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0));

  • ../trunk-jpl/src/c/cores/control_core.cpp

     
    9696        if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run!
    9797                femmodel->OutputControlsx(&femmodel->results);
    9898
    99                 #ifdef _HAVE_ADOLC_
     99                #ifdef _HAVE_AD_
    100100                IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
    101101                for(int i=0;i<nsteps;i++) J_passive[i]=reCast<IssmPDouble>(J[i]);
    102102                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.