Index: /issm/trunk-jpl/externalpackages/python/install-2.7.3-macosx-lion.sh =================================================================== --- /issm/trunk-jpl/externalpackages/python/install-2.7.3-macosx-lion.sh (revision 13539) +++ /issm/trunk-jpl/externalpackages/python/install-2.7.3-macosx-lion.sh (revision 13540) @@ -8,5 +8,5 @@ #Download from ISSM server -$ISSM_DIR/scripts/DownloadExternalPackage.py http://issm.jpl.nasa.gov/files/externalpackages/Python-2.7.3.tgz Python-2.7.3.tgz +#$ISSM_DIR/scripts/DownloadExternalPackage.py http://issm.jpl.nasa.gov/files/externalpackages/Python-2.7.3.tgz Python-2.7.3.tgz #Untar and move python into install directory @@ -17,6 +17,5 @@ #Configure and compile cd src -./configure \ - --enable-framework="$ISSM_DIR/externalpackages/python/install/Library/Frameworks" +./configure --enable-framework="$ISSM_DIR/externalpackages/python/install/Library/Frameworks" if [ $# -eq 0 ]; then make Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13539) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13540) @@ -327,6 +327,4 @@ ./modules/RequestedDependentsx/RequestedDependentsx.h\ ./modules/RequestedDependentsx/RequestedDependentsx.cpp\ - ./modules/AutodiffDriversx/AutodiffDriversx.h\ - ./modules/AutodiffDriversx/AutodiffDriversx.cpp\ ./modules/ResetConstraintsx/ResetConstraintsx.h\ ./modules/ResetConstraintsx/ResetConstraintsx.cpp\ @@ -350,6 +348,8 @@ ./solutions/AnalysisConfiguration.cpp\ ./solutions/CorePointerFromSolutionEnum.cpp\ + ./solutions/PureCorePointerFromSolutionEnum.cpp\ ./solutions/EnvironmentInit.cpp\ ./solutions/EnvironmentFinalize.cpp\ + ./solutions/ad_core.cpp\ ./solvers/solver_linear.cpp\ ./solvers/solver_nonlinear.cpp\ Index: /issm/trunk-jpl/src/c/classes/objects/Profiler.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Profiler.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/classes/objects/Profiler.cpp (revision 13540) @@ -17,4 +17,6 @@ Profiler::Profiler(){ this->time=new Parameters(); + this->flops=new Parameters(); + this->memory=new Parameters(); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Profiler.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Profiler.h (revision 13539) +++ /issm/trunk-jpl/src/c/classes/objects/Profiler.h (revision 13540) @@ -21,4 +21,6 @@ StartCore, FinishCore, + StartAdCore, + FinishAdCore, Finish }; Index: sm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp (revision 13539) +++ (revision ) @@ -1,286 +1,0 @@ -/*!\file AutodiffDriversx - * \brief: compute outputs from the AD mode, using our dependents and independents, and drivers available in Adolc. - */ -#include -#include "../../modules/modules.h" -#include "../../shared/shared.h" -#include "../../include/include.h" -#include "../../toolkits/toolkits.h" -#include "../../EnumDefinitions/EnumDefinitions.h" - - -void AutodiffDriversx(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,Results* results){ - - - /*diverse: */ - int i; - int dummy; - - bool isautodiff = false; - int num_dependents; - int num_independents; - char* driver=NULL; - - /*state variables: */ - IssmDouble *axp = NULL; - double *xp = NULL; - - /*AD mode on?: */ - parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - - if(isautodiff){ - - #ifdef _HAVE_ADOLC_ - - if(VerboseAutodiff())_pprintLine_(" start AD driver"); - - /*preliminary checks: */ - parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); - parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); - if(!(num_dependents*num_independents)) return; - - /*retrieve state variable: */ - parameters->FindParam(&axp,&dummy,AutodiffXpEnum); - - /* driver argument */ - xp=xNew(num_independents); - for(i=0;i(axp[i]); - } - - /*get the EDF pointer:*/ - ext_diff_fct *anEDF_for_solverx_p=dynamic_cast * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; - - /*Branch according to AD driver: */ - parameters->FindParam(&driver,AutodiffDriverEnum); - - /* these are always needed regardless of the interpreter */ - anEDF_for_solverx_p->dp_x=xNew(anEDF_for_solverx_p->max_n); - anEDF_for_solverx_p->dp_y=xNew(anEDF_for_solverx_p->max_m); - - if (strcmp(driver,"fos_forward")==0){ - - int anIndepIndex; - double *tangentDir=NULL; - double *jacTimesTangentDir=NULL; - double *theOutput=NULL; - - /*retrieve direction index: */ - parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum); - - if (anIndepIndex<0 || anIndepIndex>=num_independents) _error_("index value for AutodiffFosForwardIndexEnum should be in [0,num_independents-1]"); - - tangentDir=xNewZeroInit(num_independents); - tangentDir[anIndepIndex]=1.0; - - jacTimesTangentDir=xNew(num_dependents); - theOutput=xNew(num_dependents); - - /*set the forward method function pointer: */ - anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; - - /*allocate the space for the parameters to invoke the EDF fos_forward:*/ - anEDF_for_solverx_p->dp_X=xNew(anEDF_for_solverx_p->max_n); - anEDF_for_solverx_p->dp_Y=xNew(anEDF_for_solverx_p->max_m); - - /*call driver: */ - fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir ); - - /*add to results*/ - results->AddObject(new GenericExternalResult(results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,1,0.0)); - - /*free resources :*/ - xDelete(theOutput); - xDelete(jacTimesTangentDir); - xDelete(tangentDir); - } - else if ((strcmp(driver,"fov_forward")==0) || (strcmp(driver,"fov_forward_all")==0)){ - - int* indepIndices=NULL; - int tangentDirNum; - int dummy; - double **jacTimesSeed=NULL; - double **seed=NULL; - double *theOutput=NULL; - std::set anIndexSet; - - - /*retrieve directions:*/ - if (strcmp(driver,"fov_forward_all")==0){ - tangentDirNum=num_independents; - indepIndices=xNewZeroInit(tangentDirNum); - for(i=0;iFindParam(&indepIndices,&tangentDirNum,&dummy,AutodiffFovForwardIndicesEnum); - } - - /*Some checks: */ - if (tangentDirNum<1 || tangentDirNum>num_independents) _error_("tangentDirNum should be in [1,num_independents]"); - - /* full Jacobian or Jacobian projection:*/ - jacTimesSeed=xNew(num_dependents,tangentDirNum); - - /*set the forward method function pointers: */ - anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx; - // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; - - /*allocate the space for the parameters to invoke EDF fov_forward:*/ - anEDF_for_solverx_p->dpp_X=xNew(anEDF_for_solverx_p->max_n, tangentDirNum); - anEDF_for_solverx_p->dpp_Y=xNew(anEDF_for_solverx_p->max_m, tangentDirNum); - - /*seed matrix: */ - seed=xNewZeroInit(num_independents,tangentDirNum); - - /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/ - for (int i=0; inum_independents) { - _error_("indepIndices values must be in [0,num_independents-1]"); - } - if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) { - _error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression"); - } - anIndexSet.insert(indepIndices[i]); - /* now populate the seed matrix from the set of independent indices; - * simple setup with a single 1.0 per column and at most a single 1.0 per row*/ - seed[indepIndices[i]][i]=1.0; - } - - /*allocate output: */ - theOutput=xNew(num_dependents); - - /*call driver: */ - fov_forward(1,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed ); - /*Free resources: */ - xDelete(theOutput); - xDelete(indepIndices); - xDelete(seed); - - /*add to results: */ - results->AddObject(new GenericExternalResult(results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,1,0.0)); - - /*Free resources: */ - xDelete(jacTimesSeed); - xDelete(indepIndices); - } - else if (strcmp(driver,"fos_reverse")==0) { - - int aDepIndex=0; - double *aWeightVector=NULL; - double *weightVectorTimesJac=NULL; - - /*retrieve direction index: */ - parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum); - - if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]"); - - aWeightVector=xNewZeroInit(num_dependents); - aWeightVector[aDepIndex]=1.0; - - weightVectorTimesJac=xNew(num_independents); - - /*set the forward method function pointer: */ - anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; - - /*allocate the space for the parameters to invoke the EDF fos_reverse :*/ - anEDF_for_solverx_p->dp_U=xNew(anEDF_for_solverx_p->max_m); - anEDF_for_solverx_p->dp_Z=xNew(anEDF_for_solverx_p->max_n); - - /*call driver: */ - fos_reverse(1,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); - - /*add to results*/ - results->AddObject(new GenericExternalResult(results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,1,0.0)); - - /*free resources :*/ - xDelete(weightVectorTimesJac); - xDelete(aWeightVector); - } - else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){ - - int* depIndices=NULL; - int weightNum; - int dummy; - double **weightsTimesJac=NULL; - double **weights=NULL; - std::set anIndexSet; - - - /*retrieve directions:*/ - if (strcmp(driver,"fov_reverse_all")==0){ - weightNum=num_dependents; - depIndices=xNewZeroInit(weightNum); - for(i=0;iFindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum); - } - - /*Some checks: */ - if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]"); - - /* full Jacobian or Jacobian projection:*/ - weightsTimesJac=xNew(weightNum,num_independents); - - /*set the forward method function pointers: */ - anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; - - /*allocate the space for the parameters to invoke the EDF fos_reverse :*/ - anEDF_for_solverx_p->dpp_U=xNew(weightNum,anEDF_for_solverx_p->max_m); - anEDF_for_solverx_p->dpp_Z=xNew(weightNum,anEDF_for_solverx_p->max_n); - - /*seed matrix: */ - weights=xNewZeroInit(weightNum,num_dependents); - - /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/ - for (int i=0; inum_dependents) { - _error_("depIndices values must be in [0,num_dependents-1]"); - } - if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) { - _error_("duplicate depIndices values are not allowed until we implement Jacobian decompression"); - } - anIndexSet.insert(depIndices[i]); - /* now populate the seed matrix from the set of independent indices; - * simple setup with a single 1.0 per column and at most a single 1.0 per row*/ - weights[depIndices[i]][i]=1.0; - } - - /*call driver: */ - fov_reverse(1,num_dependents,num_independents, weightNum, weights, weightsTimesJac ); - - /*add to results: */ - results->AddObject(new GenericExternalResult(results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,1,0.0)); - - /*Free resources: */ - xDelete(weights); - xDelete(weightsTimesJac); - xDelete(depIndices); - } - else _error_("driver: " << driver << " not yet supported!"); - - /* delete the allocated space for the parameters:*/ - xDelete(anEDF_for_solverx_p->dp_x); - xDelete(anEDF_for_solverx_p->dp_X); - xDelete(anEDF_for_solverx_p->dpp_X); - xDelete(anEDF_for_solverx_p->dp_y); - xDelete(anEDF_for_solverx_p->dp_Y); - xDelete(anEDF_for_solverx_p->dpp_Y); - xDelete(anEDF_for_solverx_p->dp_U); - xDelete(anEDF_for_solverx_p->dpp_U); - xDelete(anEDF_for_solverx_p->dp_Z); - xDelete(anEDF_for_solverx_p->dpp_Z); - - /*Free resources: */ - xDelete(xp); - xDelete(axp); - - if(VerboseAutodiff())_pprintLine_(" end AD driver"); - - #else - _error_("Should not be requesting AD drivers when an AD library is not available!"); - #endif - } -} Index: sm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.h (revision 13539) +++ (revision ) @@ -1,13 +1,0 @@ -/*!\file: AutodiffDriversx.h - * \brief header file for requesting AD results (gradient, jacobian, etc ...) - */ - -#ifndef _AUTODIFF_DRIVERSX_H_ -#define _AUTODIFF_DRIVERSX_H_ - -#include "../../Container/Container.h" - -/* local prototypes: */ -void AutodiffDriversx(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,Results* results); - -#endif /* _AUTODIFF_DRIVERSX_H_*/ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 13539) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 13540) @@ -8,5 +8,4 @@ /*Modules: */ #include "./AddExternalResultx/AddExternalResultx.h" -#include "./AutodiffDriversx/AutodiffDriversx.h" #include "./AverageFilterx/AverageFilterx.h" #include "./AverageOntoPartitionx/AverageOntoPartitionx.h" Index: /issm/trunk-jpl/src/c/solutions/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 13540) @@ -18,86 +18,40 @@ #include "../solvers/solvers.h" -void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){ +void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore){ /*output: */ void (*solutioncore)(FemModel*)=NULL; - switch(solutiontype){ + /*parameters: */ + bool control_analysis=false; + bool tao_analysis=false; + bool dakota_analysis=false; - case DiagnosticSolutionEnum: - #ifdef _HAVE_DIAGNOSTIC_ - solutioncore=&diagnostic_core; - #else - _error_("ISSM was not compiled with diagnostic capabilities. Exiting"); - #endif - break; - case SteadystateSolutionEnum: - #ifdef _HAVE_STEADYSTATE_ - solutioncore=&steadystate_core; - #else - _error_("ISSM was not compiled with steady state capabilities. Exiting"); - #endif - break; - case ThermalSolutionEnum: - #ifdef _HAVE_THERMAL_ - solutioncore=&thermal_core; - #else - _error_("ISSM was not compiled with thermal capabilities. Exiting"); - #endif - break; - case EnthalpySolutionEnum: - #ifdef _HAVE_THERMAL_ - solutioncore=&enthalpy_core; - #else - _error_("ISSM was not compiled with thermal capabilities. Exiting"); - #endif - break; - case BalancethicknessSolutionEnum: - #ifdef _HAVE_BALANCED_ - solutioncore=&balancethickness_core; - #else - _error_("ISSM was not compiled with balanced capabilities. Exiting"); - #endif - break; - case HydrologySolutionEnum: - #ifdef _HAVE_HYDROLOGY_ - solutioncore=&hydrology_core; - #else - _error_("ISSM was not compiled with hydrology capabilities. Exiting"); - #endif - break; - case SurfaceSlopeSolutionEnum: - #ifdef _HAVE_SLOPE_ - solutioncore=&surfaceslope_core; - #else - _error_("ISSM was not compiled with slope capabilities. Exiting"); - #endif - break; - case BedSlopeSolutionEnum: - #ifdef _HAVE_SLOPE_ - solutioncore=&bedslope_core; - #else - _error_("ISSM was not compiled with slope capabilities. Exiting"); - #endif - break; - case TransientSolutionEnum: - #ifdef _HAVE_TRANSIENT_ - solutioncore=&transient_core; - #else - _error_("ISSM was not compiled with transient capabilities. Exiting"); - #endif - break; - case PrognosticSolutionEnum: - #ifdef _HAVE_PROGNOSTIC_ - solutioncore=&prognostic_core; - #else - _error_("ISSM was not compiled with prognostic capabilities. Exiting"); - #endif - break; - default: - _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!"); - break; + /* retrieve some parameters that tell us whether wrappers are allowed, or whether we return + * a pure core. Wrappers can be dakota_core (which samples many solution_cores) or control_core (which + * carries out adjoint based inversion on a certain core: */ + parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + parameters->FindParam(&control_analysis,InversionIscontrolEnum); + parameters->FindParam(&tao_analysis,InversionTaoEnum); + + if(nodakotacore)dakota_analysis=false; + + if(dakota_analysis){ + #ifdef _HAVE_DAKOTA_ + solutioncore=dakota_core; + #else + _error_("ISSM was not compiled with dakota support, cannot carry out dakota analysis!"); + #endif } - + else if(control_analysis){ + #ifdef _HAVE_CONTROL_ + if(tao_analysis) solutioncore=controltao_core; + else solutioncore=control_core; + #else + _error_("ISSM was not compiled with control support, cannot carry out control analysis!"); + #endif + } + else PureCorePointerFromSolutionEnum(&solutioncore,parameters,solutiontype); /*This means we retrieve a core solution that is not a wrapper*/ + /*Assign output pointer:*/ _assert_(psolutioncore); Index: /issm/trunk-jpl/src/c/solutions/DakotaSpawnCore.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/DakotaSpawnCore.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/DakotaSpawnCore.cpp (revision 13540) @@ -55,4 +55,5 @@ void (*solutioncore)(FemModel*) = NULL; FemModel* femmodel = NULL; + bool nodakotacore = true; @@ -81,10 +82,5 @@ /*Determine solution sequence: */ if(VerboseQmu()) _pprintLine_("Starting " << EnumToStringx(solution_type) << " core:"); - CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); - #ifdef _HAVE_CONTROL_ - if(control_analysis)solutioncore=&control_core; - #else - _error_("ISSM was not compiled with control capabilities, exiting!"); - #endif + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type,nodakotacore); /*Run the core solution sequence: */ @@ -97,4 +93,6 @@ /*Free ressources:*/ DakotaFree(&d_variables,&d_variables_descriptors,&responses_descriptors, d_numvariables, numresponsedescriptors); + + return 1; //this is critical! do not return 0, otherwise, dakota_core will stop running! } Index: /issm/trunk-jpl/src/c/solutions/PureCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/PureCorePointerFromSolutionEnum.cpp (revision 13540) +++ /issm/trunk-jpl/src/c/solutions/PureCorePointerFromSolutionEnum.cpp (revision 13540) @@ -0,0 +1,106 @@ +/*!\file: PureCorePointerFromSolutionEnum.cpp + * \brief: return type of analyses, number of analyses and core solution function. + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" + +void PureCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){ + + /*output: */ + void (*solutioncore)(FemModel*)=NULL; + + switch(solutiontype){ + + case DiagnosticSolutionEnum: + #ifdef _HAVE_DIAGNOSTIC_ + solutioncore=&diagnostic_core; + #else + _error_("ISSM was not compiled with diagnostic capabilities. Exiting"); + #endif + break; + case SteadystateSolutionEnum: + #ifdef _HAVE_STEADYSTATE_ + solutioncore=&steadystate_core; + #else + _error_("ISSM was not compiled with steady state capabilities. Exiting"); + #endif + break; + case ThermalSolutionEnum: + #ifdef _HAVE_THERMAL_ + solutioncore=&thermal_core; + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + break; + case EnthalpySolutionEnum: + #ifdef _HAVE_THERMAL_ + solutioncore=&enthalpy_core; + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + break; + case BalancethicknessSolutionEnum: + #ifdef _HAVE_BALANCED_ + solutioncore=&balancethickness_core; + #else + _error_("ISSM was not compiled with balanced capabilities. Exiting"); + #endif + break; + case HydrologySolutionEnum: + #ifdef _HAVE_HYDROLOGY_ + solutioncore=&hydrology_core; + #else + _error_("ISSM was not compiled with hydrology capabilities. Exiting"); + #endif + break; + case SurfaceSlopeSolutionEnum: + #ifdef _HAVE_SLOPE_ + solutioncore=&surfaceslope_core; + #else + _error_("ISSM was not compiled with slope capabilities. Exiting"); + #endif + break; + case BedSlopeSolutionEnum: + #ifdef _HAVE_SLOPE_ + solutioncore=&bedslope_core; + #else + _error_("ISSM was not compiled with slope capabilities. Exiting"); + #endif + break; + case TransientSolutionEnum: + #ifdef _HAVE_TRANSIENT_ + solutioncore=&transient_core; + #else + _error_("ISSM was not compiled with transient capabilities. Exiting"); + #endif + break; + case PrognosticSolutionEnum: + #ifdef _HAVE_PROGNOSTIC_ + solutioncore=&prognostic_core; + #else + _error_("ISSM was not compiled with prognostic capabilities. Exiting"); + #endif + break; + default: + _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!"); + break; + } + + /*Assign output pointer:*/ + _assert_(psolutioncore); + *psolutioncore=solutioncore; + +} Index: /issm/trunk-jpl/src/c/solutions/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/ad_core.cpp (revision 13540) +++ /issm/trunk-jpl/src/c/solutions/ad_core.cpp (revision 13540) @@ -0,0 +1,301 @@ +/*!\file ad_core + * \brief: compute outputs from the AD mode, using our dependents and independents, and drivers available in Adolc. + */ + +/*Includes: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +#include +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/shared.h" +#include "../io/io.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" +/*}}}*/ + +void ad_core(FemModel* femmodel){ + + /*diverse: */ + int i; + int dummy; + + bool isautodiff = false; + int num_dependents; + int num_independents; + char* driver=NULL; + + /*state variables: */ + IssmDouble *axp = NULL; + double *xp = NULL; + + /*AD mode on?: */ + femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + + if(isautodiff){ + + #ifdef _HAVE_ADOLC_ + + if(VerboseAutodiff())_pprintLine_(" start ad core"); + + /*First, stop tracing: */ + trace_off(); + + /*preliminary checks: */ + femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); + femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); + if(!(num_dependents*num_independents)) return; + + /*retrieve state variable: */ + femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); + + /* driver argument */ + xp=xNew(num_independents); + for(i=0;i(axp[i]); + } + + /*get the EDF pointer:*/ + ext_diff_fct *anEDF_for_solverx_p=dynamic_cast * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; + + /*Branch according to AD driver: */ + femmodel->parameters->FindParam(&driver,AutodiffDriverEnum); + + /* these are always needed regardless of the interpreter */ + anEDF_for_solverx_p->dp_x=xNew(anEDF_for_solverx_p->max_n); + anEDF_for_solverx_p->dp_y=xNew(anEDF_for_solverx_p->max_m); + + if (strcmp(driver,"fos_forward")==0){ + + int anIndepIndex; + double *tangentDir=NULL; + double *jacTimesTangentDir=NULL; + double *theOutput=NULL; + + /*retrieve direction index: */ + femmodel->parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum); + + if (anIndepIndex<0 || anIndepIndex>=num_independents) _error_("index value for AutodiffFosForwardIndexEnum should be in [0,num_independents-1]"); + + tangentDir=xNewZeroInit(num_independents); + tangentDir[anIndepIndex]=1.0; + + jacTimesTangentDir=xNew(num_dependents); + theOutput=xNew(num_dependents); + + /*set the forward method function pointer: */ + anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; + + /*allocate the space for the parameters to invoke the EDF fos_forward:*/ + anEDF_for_solverx_p->dp_X=xNew(anEDF_for_solverx_p->max_n); + anEDF_for_solverx_p->dp_Y=xNew(anEDF_for_solverx_p->max_m); + + /*call driver: */ + fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir ); + + /*add to results*/ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,1,0.0)); + + /*free resources :*/ + xDelete(theOutput); + xDelete(jacTimesTangentDir); + xDelete(tangentDir); + } + else if ((strcmp(driver,"fov_forward")==0) || (strcmp(driver,"fov_forward_all")==0)){ + + int* indepIndices=NULL; + int tangentDirNum; + int dummy; + double **jacTimesSeed=NULL; + double **seed=NULL; + double *theOutput=NULL; + std::set anIndexSet; + + + /*retrieve directions:*/ + if (strcmp(driver,"fov_forward_all")==0){ + tangentDirNum=num_independents; + indepIndices=xNewZeroInit(tangentDirNum); + for(i=0;iparameters->FindParam(&indepIndices,&tangentDirNum,&dummy,AutodiffFovForwardIndicesEnum); + } + + /*Some checks: */ + if (tangentDirNum<1 || tangentDirNum>num_independents) _error_("tangentDirNum should be in [1,num_independents]"); + + /* full Jacobian or Jacobian projection:*/ + jacTimesSeed=xNew(num_dependents,tangentDirNum); + + /*set the forward method function pointers: */ + anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx; + // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; + + /*allocate the space for the parameters to invoke EDF fov_forward:*/ + anEDF_for_solverx_p->dpp_X=xNew(anEDF_for_solverx_p->max_n, tangentDirNum); + anEDF_for_solverx_p->dpp_Y=xNew(anEDF_for_solverx_p->max_m, tangentDirNum); + + /*seed matrix: */ + seed=xNewZeroInit(num_independents,tangentDirNum); + + /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/ + for (int i=0; inum_independents) { + _error_("indepIndices values must be in [0,num_independents-1]"); + } + if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) { + _error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression"); + } + anIndexSet.insert(indepIndices[i]); + /* now populate the seed matrix from the set of independent indices; + * simple setup with a single 1.0 per column and at most a single 1.0 per row*/ + seed[indepIndices[i]][i]=1.0; + } + + /*allocate output: */ + theOutput=xNew(num_dependents); + + /*call driver: */ + fov_forward(1,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed ); + /*Free resources: */ + xDelete(theOutput); + xDelete(indepIndices); + xDelete(seed); + + /*add to results: */ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,1,0.0)); + + /*Free resources: */ + xDelete(jacTimesSeed); + xDelete(indepIndices); + } + else if (strcmp(driver,"fos_reverse")==0) { + + int aDepIndex=0; + double *aWeightVector=NULL; + double *weightVectorTimesJac=NULL; + + /*retrieve direction index: */ + femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum); + + if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]"); + + aWeightVector=xNewZeroInit(num_dependents); + aWeightVector[aDepIndex]=1.0; + + weightVectorTimesJac=xNew(num_independents); + + /*set the forward method function pointer: */ + anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; + + /*allocate the space for the parameters to invoke the EDF fos_reverse :*/ + anEDF_for_solverx_p->dp_U=xNew(anEDF_for_solverx_p->max_m); + anEDF_for_solverx_p->dp_Z=xNew(anEDF_for_solverx_p->max_n); + + /*call driver: */ + fos_reverse(1,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); + + /*add to results*/ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,1,0.0)); + + /*free resources :*/ + xDelete(weightVectorTimesJac); + xDelete(aWeightVector); + } + else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){ + + int* depIndices=NULL; + int weightNum; + int dummy; + double **weightsTimesJac=NULL; + double **weights=NULL; + std::set anIndexSet; + + + /*retrieve directions:*/ + if (strcmp(driver,"fov_reverse_all")==0){ + weightNum=num_dependents; + depIndices=xNewZeroInit(weightNum); + for(i=0;iparameters->FindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum); + } + + /*Some checks: */ + if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]"); + + /* full Jacobian or Jacobian projection:*/ + weightsTimesJac=xNew(weightNum,num_independents); + + /*set the forward method function pointers: */ + anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; + + /*allocate the space for the parameters to invoke the EDF fos_reverse :*/ + anEDF_for_solverx_p->dpp_U=xNew(weightNum,anEDF_for_solverx_p->max_m); + anEDF_for_solverx_p->dpp_Z=xNew(weightNum,anEDF_for_solverx_p->max_n); + + /*seed matrix: */ + weights=xNewZeroInit(weightNum,num_dependents); + + /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/ + for (int i=0; inum_dependents) { + _error_("depIndices values must be in [0,num_dependents-1]"); + } + if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) { + _error_("duplicate depIndices values are not allowed until we implement Jacobian decompression"); + } + anIndexSet.insert(depIndices[i]); + /* now populate the seed matrix from the set of independent indices; + * simple setup with a single 1.0 per column and at most a single 1.0 per row*/ + weights[depIndices[i]][i]=1.0; + } + + /*call driver: */ + fov_reverse(1,num_dependents,num_independents, weightNum, weights, weightsTimesJac ); + + /*add to results: */ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,1,0.0)); + + /*Free resources: */ + xDelete(weights); + xDelete(weightsTimesJac); + xDelete(depIndices); + } + else _error_("driver: " << driver << " not yet supported!"); + + /* delete the allocated space for the parameters:*/ + xDelete(anEDF_for_solverx_p->dp_x); + xDelete(anEDF_for_solverx_p->dp_X); + xDelete(anEDF_for_solverx_p->dpp_X); + xDelete(anEDF_for_solverx_p->dp_y); + xDelete(anEDF_for_solverx_p->dp_Y); + xDelete(anEDF_for_solverx_p->dpp_Y); + xDelete(anEDF_for_solverx_p->dp_U); + xDelete(anEDF_for_solverx_p->dpp_U); + xDelete(anEDF_for_solverx_p->dp_Z); + xDelete(anEDF_for_solverx_p->dpp_Z); + + /*Free resources: */ + xDelete(xp); + xDelete(axp); + + if(VerboseAutodiff())_pprintLine_(" end AD core"); + + #else + _error_("Should not be requesting AD drivers when an AD library is not available!"); + #endif + } +} Index: /issm/trunk-jpl/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/control_core.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/control_core.cpp (revision 13540) @@ -62,5 +62,5 @@ /*out of solution_type, figure out solution core and adjoint function pointer*/ - CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + PureCorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 13540) @@ -94,5 +94,5 @@ femmodel->parameters->SetParam(true,SaveResultsEnum); void (*solutioncore)(FemModel*)=NULL; - CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + PureCorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); solutioncore(femmodel); Index: /issm/trunk-jpl/src/c/solutions/dakota_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/dakota_core.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/dakota_core.cpp (revision 13540) @@ -65,11 +65,8 @@ Parameters* parameters = NULL; - /*Retrieve parameters: */ - parameters=femmodel->parameters; - /*Recover dakota_input_file, dakota_output_file and dakota_error_file, in the parameters dataset in parallel */ - parameters->FindParam(&dakota_input_file,QmuInNameEnum); - parameters->FindParam(&dakota_output_file,QmuOutNameEnum); - parameters->FindParam(&dakota_error_file,QmuErrNameEnum); + femmodel->parameters->FindParam(&dakota_input_file,QmuInNameEnum); + femmodel->parameters->FindParam(&dakota_output_file,QmuOutNameEnum); + femmodel->parameters->FindParam(&dakota_error_file,QmuErrNameEnum); if(my_rank==0){ Index: /issm/trunk-jpl/src/c/solutions/issm.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/issm.cpp (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/issm.cpp (revision 13540) @@ -15,5 +15,4 @@ FILE *petscoptionsfid = NULL; bool waitonlock = false; - bool dakota_analysis,control_analysis,tao_analysis; /*AD: */ @@ -79,7 +78,4 @@ /*get parameters: */ femmodel->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum); - femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); - femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum); - femmodel->parameters->FindParam(&tao_analysis,InversionTaoEnum); femmodel->parameters->FindParam(&autodiff,AutodiffIsautodiffEnum); @@ -104,34 +100,7 @@ _pprintLine_("call computational core:"); - - - profiler->Tag(StartCore); - if(dakota_analysis){ - #ifdef _HAVE_DAKOTA_ - dakota_core(femmodel); - #else - _error_("ISSM was not compiled with dakota support, cannot carry out dakota analysis!"); - #endif - } - else if(control_analysis){ - #ifdef _HAVE_CONTROL_ - if(tao_analysis) controltao_core(femmodel); - else control_core(femmodel); - #else - _error_("ISSM was not compiled with control support, cannot carry out control analysis!"); - #endif - } - else{ - solutioncore(femmodel); - } - - #ifdef _HAVE_ADOLC_ - if(autodiff){ - trace_off(); - AutodiffDriversx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); - } - #endif - profiler->Tag(FinishCore); - + profiler->Tag(StartCore); solutioncore(femmodel); profiler->Tag(FinishCore); + profiler->Tag(StartAdCore); ad_core(femmodel); profiler->Tag(FinishAdCore); + ProfilerEnd(profiler,femmodel->results,femmodel->parameters); @@ -220,7 +189,7 @@ /*Add to results: */ - results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingSolutionTimeEnum, solution_time, 1, 0)); - results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentMemEnum, solution_memory, 1, 0)); - results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops, 1, 0)); + results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingSolutionTimeEnum, solution_time, 1, 0)); + results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentMemEnum, solution_memory, 1, 0)); + results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops, 1, 0)); } Index: /issm/trunk-jpl/src/c/solutions/solutions.h =================================================================== --- /issm/trunk-jpl/src/c/solutions/solutions.h (revision 13539) +++ /issm/trunk-jpl/src/c/solutions/solutions.h (revision 13540) @@ -31,4 +31,5 @@ void transient_core(FemModel* femmodel); void dakota_core(FemModel* femmodel); +void ad_core(FemModel* femmodel); IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs); @@ -52,5 +53,6 @@ //solution configuration void AnalysisConfiguration(int** panalyses,int* pnumanalyses, int solutiontype); -void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype); +void PureCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype); +void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore=false); void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype);