Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23242) @@ -3133,5 +3133,8 @@ } - #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable. + #if defined(_HAVE_AD_) + /*we want to avoid the round operation at all cost. Not differentiable.*/ + _error_("not implemented yet"); + #else dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd; dMass = round(dMass * 100.0)/100.0; Index: /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 23241) +++ /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 23242) @@ -622,5 +622,5 @@ /*Specific instantiations for IssmDouble*: */ -#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. +#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. template <> inline GenericExternalResult::~GenericExternalResult(){ /*{{{*/ xDelete(result_name); @@ -709,5 +709,5 @@ return new GenericExternalResult*>(this->id,StringToEnumx(this->result_name),this->value,this->step,this->time); } /*}}}*/ -#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. +#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. template <> inline void GenericExternalResult*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23242) @@ -881,5 +881,5 @@ results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops)); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ solution_time = profiler->TotalTime(ADCORE); solution_flops = profiler->TotalFlops(ADCORE); @@ -1845,5 +1845,5 @@ if(isautodiff){ - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); Index: /issm/trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 23242) @@ -97,5 +97,5 @@ femmodel->OutputControlsx(&femmodel->results); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ IssmPDouble* J_passive=xNew(nsteps); for(int i=0;i(J[i]); Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 23242) @@ -88,5 +88,5 @@ /*output*/ - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ IssmPDouble* J_passive=xNew(2*num); for(int i=0;i<2*num;i++) J_passive[i]=reCast(output[i]); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23242) @@ -73,6 +73,34 @@ xDelete(options); /*}}}*/ + #elif _HAVE_CODIPACK_ + //fprintf(stderr, "*** Codipack CreateParametersAutodiff()\n"); + /*initialize a placeholder to store solver pointers: {{{*/ + /*Solver pointers depend on what type of solver we are implementing: */ + options=OptionsFromAnalysis(parameters,DefaultAnalysisEnum); //options database is not filled in yet, use default. + ToolkitOptions::Init(options); + + switch(IssmSolverTypeFromToolkitOptions()){ + case MumpsEnum:{ + #ifndef _HAVE_MUMPS_ + _error_("CoDiPack: requesting mumps solver without MUMPS being compiled in!"); + #endif + break; + } + case GslEnum: { + #ifndef _HAVE_GSL_ + _error_("CoDiPack: requesting GSL solver without GSL being compiled in!"); + #endif + break; + } + default: + _error_("solver type not supported yet!"); + } + /*Free ressources: */ + xDelete(options); + #endif + #if defined(_HAVE_AD_) if(isautodiff){ +#if _HAVE_ADOLC_ /*Copy some parameters from IoModel to parameters dataset: {{{*/ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum)); @@ -83,4 +111,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum)); /*}}}*/ +#endif /*retrieve driver: {{{*/ iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver"); @@ -88,4 +117,8 @@ if(strcmp(autodiff_driver,"fos_forward")==0){ +#if _HAVE_CODIPACK_ + // FIXME codi support Foward Mode (scalar) + _error_("Foward Mode (scalar) not supported yet!"); +#endif parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_forward_index",AutodiffFosForwardIndexEnum)); } @@ -94,4 +127,8 @@ } else if(strcmp(autodiff_driver,"fov_forward")==0){ +#if _HAVE_CODIPACK_ + // FIXME codi support Foward Mode (vector) + _error_("Foward Mode (vector) not supported yet!"); +#endif /*Retrieve list of indices: */ iomodel->FetchData(&indices,&num_indices,&dummy,"md.autodiff.fov_forward_indices"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23242) @@ -374,5 +374,5 @@ ParseToolkitsOptionsx(parameters,toolkitsoptionsfid); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ if(VerboseMProcessor()) _printf0_(" starting autodiff parameters \n"); CreateParametersAutodiff(parameters,iomodel); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23242) @@ -62,5 +62,5 @@ //into specified top structure depth (zTop). Also make sure top grid cell //structure length (dzTop) is greater than 5 cm - #ifndef _HAVE_ADOLC_ //avoid the round operation check! + #ifndef _HAVE_AD_ //avoid the round operation check! if (dgpTop != round(dgpTop)){ _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop"); @@ -425,9 +425,9 @@ // spectral range: // 0.3 - 0.8um - IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5)); + IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5)); // 0.8 - 1.5um - IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5)); + IssmDouble a1 = max(0, 0.95 - 15.4 *pow(gsz,0.5)); // 1.5 - 2.8um - IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5)); + IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5)); // broadband surface albedo @@ -667,5 +667,5 @@ // NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability dt=1e12; - for(int i=0;i 1/2 @@ -752,5 +752,5 @@ // less when Ts is taken as the mean of the x top grid cells. Ts = (T[0] + T[1])/2.0; - Ts = fmin(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC) + Ts = min(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC) //TURBULENT HEAT FLUX @@ -767,5 +767,5 @@ // do not allow Ri to exceed 0.19 - Ri = fmin(Ri, 0.19); + Ri = min(Ri, 0.19); // calculate momentum 'coefM' stability factor @@ -947,9 +947,9 @@ // spectral albedos: // 0.3 - 0.8um - IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5)); + IssmDouble a0 = min(0.98, 1.0 - 1.58 *pow(gsz[0],0.5)); // 0.8 - 1.5um - IssmDouble a1 = fmax(0.0, 0.95 - 15.4 *pow(gsz[0],0.5)); + IssmDouble a1 = max(0.0, 0.95 - 15.4 *pow(gsz[0],0.5)); // 1.5 - 2.8um - IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5)); + IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5)); // separate net shortwave radiative flux into spectral ranges @@ -1204,5 +1204,5 @@ mass_diff = mass - massinit - P; - #ifndef _HAVE_ADOLC_ //avoid round operation. only check in forward mode. + #ifndef _HAVE_AD_ //avoid round operation. only check in forward mode. mass_diff = round(mass_diff * 100.0)/100.0; if (mass_diff > 0) _error_("mass not conserved in accumulation function"); @@ -1323,9 +1323,9 @@ // calculate temperature excess above 0 deg C exsT=xNewZeroInit(n); - for(int i=0;i(n); for(int i=0;i(n); for(int i=0;i 0.0 + Ttol ){ @@ -1380,8 +1380,8 @@ T[i+1] = surpE[i]/m[i+1]/CI + T[i+1]; - exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1]; - T[i+1] = fmin(CtoK, T[i+1]); - - surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI); + exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1]; + T[i+1] = min(CtoK, T[i+1]); + + surpT[i+1] = max(0.0, exsT[i+1] - LF/CI); surpE[i+1] = surpT[i+1] * CI * m[i+1]; @@ -1394,9 +1394,9 @@ // convert temperature excess to melt [kg] - for(int i=0;iW[i]-Wtol ){ dW[i]= -1*W[i]; @@ -1474,6 +1474,6 @@ if (dW[i] < 0.0-Wtol){ // excess pore water dMax = (dIce - d[i])*dz_0; // maximum refreeze - IssmDouble maxF2 = fmin(dMax, maxF[i]-F1); // maximum refreeze - F2 = fmin(-1*dW[i], maxF2); // pore water refreeze + IssmDouble maxF2 = min(dMax, maxF[i]-F1); // maximum refreeze + F2 = min(-1*dW[i], maxF2); // pore water refreeze m[i] = m[i] + F2; // mass after refreeze d[i] = m[i]/dz_0; @@ -1753,5 +1753,5 @@ /*only in forward mode! avoid round in AD mode as it is not differentiable: */ - #ifndef _HAVE_ADOLC_ + #ifndef _HAVE_AD_ dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0; dE = round(sumE0 - sumE1 - sumER + addE); @@ -1915,6 +1915,6 @@ c0arth = 0.07 * H; c1arth = 0.03 * H; - M0 = fmax(1.435 - (0.151 * log(C)),0.25); - M1 = fmax(2.366 - (0.293 * log(C)),0.25); + M0 = max(1.435 - (0.151 * log(C)),0.25); + M1 = max(2.366 - (0.293 * log(C)),0.25); c0 = M0*c0arth; c1 = M1*c1arth; @@ -1926,6 +1926,6 @@ c0arth = 0.07 * H; c1arth = 0.03 * H; - M0 = fmax(1.042 - (0.0916 * log(C)),0.25); - M1 = fmax(1.734 - (0.2039 * log(C)),0.25); + M0 = max(1.042 - (0.0916 * log(C)),0.25); + M1 = max(1.734 - (0.2039 * log(C)),0.25); c0 = M0*c0arth; c1 = M1*c1arth; @@ -2021,5 +2021,5 @@ // do not allow Ri to exceed 0.19 - Ri = fmin(Ri, 0.19); + Ri = min(Ri, 0.19); // calculate momentum 'coefM' stability factor Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23241) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23242) @@ -53,4 +53,5 @@ AutodiffNumIndependentsEnum, AutodiffObufsizeEnum, + AutodiffTapeAllocEnum, AutodiffTbufsizeEnum, AutodiffXpEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23242) @@ -61,4 +61,5 @@ case AutodiffNumIndependentsEnum : return "AutodiffNumIndependents"; case AutodiffObufsizeEnum : return "AutodiffObufsize"; + case AutodiffTapeAllocEnum : return "AutodiffTapeAlloc"; case AutodiffTbufsizeEnum : return "AutodiffTbufsize"; case AutodiffXpEnum : return "AutodiffXp"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23241) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23242) @@ -61,4 +61,5 @@ else if (strcmp(name,"AutodiffNumIndependents")==0) return AutodiffNumIndependentsEnum; else if (strcmp(name,"AutodiffObufsize")==0) return AutodiffObufsizeEnum; + else if (strcmp(name,"AutodiffTapeAlloc")==0) return AutodiffTapeAllocEnum; else if (strcmp(name,"AutodiffTbufsize")==0) return AutodiffTbufsizeEnum; else if (strcmp(name,"AutodiffXp")==0) return AutodiffXpEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum; else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum; - else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum; + if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; + else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum; else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum; else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; - else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum; + if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum; + else if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum; else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum; else if (strcmp(name,"SealevelriseHoriz")==0) return SealevelriseHorizEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; - else if (strcmp(name,"WorldComm")==0) return WorldCommEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; + if (strcmp(name,"WorldComm")==0) return WorldCommEnum; + else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum; else if (strcmp(name,"Adjoint")==0) return AdjointEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; - else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; + if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; + else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; - else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; + if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; + else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Divergence")==0) return DivergenceEnum; else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum; - else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; + if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; + else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"LoveHr")==0) return LoveHrEnum; else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; - else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"LoveKi")==0) return LoveKiEnum; + if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; + else if (strcmp(name,"LoveKi")==0) return LoveKiEnum; else if (strcmp(name,"LoveKr")==0) return LoveKrEnum; else if (strcmp(name,"LoveLi")==0) return LoveLiEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; - else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; + if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; + else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; - else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; + if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; + else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;