Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8278) +++ /issm/trunk/src/c/Makefile.am (revision 8279) @@ -1,3 +1,3 @@ -INCLUDES = @DAKOTAINCL@ @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@ @METISINCL@ @CHACOINCL@ @SCOTCHINCL@ @PLAPACKINCL@ @BLASLAPACKINCL@ @MKLINCL@ @MUMPSINCL@ @TRIANGLEINCL@ @HYPREINCL@ @MLINCL@ +INCLUDES = @DAKOTAINCL@ @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@ @METISINCL@ @CHACOINCL@ @SCOTCHINCL@ @PLAPACKINCL@ @BLASLAPACKINCL@ @MKLINCL@ @MUMPSINCL@ @TRIANGLEINCL@ @HYPREINCL@ @MLINCL@ @TAOINCL@ #Compile serial library, and then try and compile parallel library @@ -375,4 +375,5 @@ ./toolkits/petsc/patches/MatMultPatch.cpp\ ./toolkits/petsc/petscincludes.h\ + ./toolkits/tao/taoincludes.h\ ./toolkits/mpi/mpiincludes.h\ ./toolkits/mpi/patches/mpipatches.h\ @@ -1015,4 +1016,5 @@ ./toolkits/petsc/patches/MatMultPatch.cpp\ ./toolkits/petsc/petscincludes.h\ + ./toolkits/tao/taoincludes.h\ ./toolkits/mpi/mpiincludes.h\ ./toolkits/mpi/patches/mpipatches.h\ @@ -1285,4 +1287,5 @@ ./solutions/WriteLockFile.cpp\ ./solutions/control_core.cpp\ + ./solutions/controltao_core.cpp\ ./solutions/controlrestart.cpp\ ./solutions/controlconvergence.cpp\ @@ -1325,5 +1328,5 @@ endif -LDADD = ./libpISSM.a $(PETSCLIB) $(FLIBS) $(PLAPACKLIB) $(MUMPSLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) libOverload.a $(MULTITHREADINGLIB) +LDADD = ./libpISSM.a $(PETSCLIB) $(TAOLIB) $(FLIBS) $(PLAPACKLIB) $(MUMPSLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) libOverload.a $(MULTITHREADINGLIB) issm_exe_SOURCES = solutions/issm.cpp Index: /issm/trunk/src/c/shared/Matlab/ModuleBoot.cpp =================================================================== --- /issm/trunk/src/c/shared/Matlab/ModuleBoot.cpp (revision 8278) +++ /issm/trunk/src/c/shared/Matlab/ModuleBoot.cpp (revision 8279) @@ -24,4 +24,9 @@ PetscInitializeNoArguments(); + /*Initialize Tao*/ + #ifdef _HAVE_TAO_ + TaoInitialize(0,0,0,0); + #endif + return 1; } Index: /issm/trunk/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/controltao_core.cpp (revision 8279) +++ /issm/trunk/src/c/solutions/controltao_core.cpp (revision 8279) @@ -0,0 +1,143 @@ +/*!\file: control_core.cpp + * \brief: core of the control solution + */ + +#include "../toolkits/toolkits.h" +#include "../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 controltao_core(FemModel* femmodel){ + + int i,n; + + /*parameters: */ + int num_controls; + int nsteps; + double eps_cm; + double tolx; + bool cm_gradient; + int dim; + int solution_type; + bool isstokes; + bool qmu_analysis=false; + + int* control_type = NULL; + double* responses=NULL; + double* maxiter=NULL; + double* cm_jump=NULL; + + + /*intermediary: */ + double search_scalar=0; + OptArgs optargs; + OptPars optpars; + + /*Solution and Adjoint core pointer*/ + void (*solutioncore)(FemModel*)=NULL; + void (*adjointcore)(FemModel*)=NULL; + + /*output: */ + double* J=NULL; + + printf("I am in TAO!!\n"); + + /*Recover parameters used throughout the solution:{{{1*/ + femmodel->parameters->FindParam(&nsteps,NStepsEnum); + femmodel->parameters->FindParam(&num_controls,NumControlsEnum); + femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum); + femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); + femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum); + femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum); + femmodel->parameters->FindParam(&eps_cm,EpsCmEnum); + femmodel->parameters->FindParam(&tolx,TolXEnum); + femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + /*}}}*/ + + /*out of solution_type, figure out solution core and adjoint function pointer*/ + CorePointerFromSolutionEnum(&solutioncore,solution_type); + AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); + + /*Launch once a complete solution to set up all inputs*/ + _printf_(VerboseControl(),"%s\n"," preparing initial solution"); + if (isstokes) solutioncore(femmodel); + + /*Initialize responses: */ + J=(double*)xmalloc(nsteps*sizeof(double)); + + /*Initialize some of the BrentSearch arguments: */ + optargs.femmodel=femmodel; + optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ; + + /*Start looping: */ + for(n=0;nelements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n],CmResponseEnum); + + /*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */ + if (solution_type==SteadystateSolutionEnum) solutioncore(femmodel); + + _printf_(VerboseControl(),"%s\n"," compute adjoint state:"); + adjointcore(femmodel); + + gradient_core(femmodel,n,search_scalar); + + /*Return gradient if asked: */ + if (cm_gradient){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum); + goto cleanup_and_return; + } + + _printf_(VerboseControl(),"%s\n"," optimizing along gradient direction"); + optargs.n=n; optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n]; + BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs); + //OptimalSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs); + + _printf_(VerboseControl(),"%s\n"," updating parameter using optimized search scalar"); //true means update save controls + InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true); + + if(controlconvergence(J,responses,eps_cm,n)) break; + + /*Temporary saving every 5 control steps: */ + /*if (((n+1)%5)==0){ + _printf_(VerboseControl(),"%s\n"," saving temporary results"); + controlrestart(femmodel,J); + }*/ + } + + _printf_(VerboseControl(),"%s\n"," preparing final solution"); + femmodel->parameters->SetParam(false,ControlAnalysisEnum); //needed to turn control result output in solutioncore + solutioncore(femmodel); + + /*some results not computed by steadystate_core or diagnostic_core: */ + if(!qmu_analysis){ //do not save this if we are running the control core from a qmu run! + for(i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]); + femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0)); + //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToStringx(control_type),1,0)); + } + + cleanup_and_return: + /*Free ressources: */ + xfree((void**)&control_type); + xfree((void**)&responses); + xfree((void**)&maxiter); + xfree((void**)&cm_jump); + xfree((void**)&J); + + /*control_core might be used in Qmu, so leave everything similar to where it started: */ + femmodel->parameters->SetParam(true,ControlAnalysisEnum); + + /* Finalize TAO */ + TaoFinalize(); +} Index: /issm/trunk/src/c/solutions/issm.cpp =================================================================== --- /issm/trunk/src/c/solutions/issm.cpp (revision 8278) +++ /issm/trunk/src/c/solutions/issm.cpp (revision 8279) @@ -79,5 +79,11 @@ /*if control is being run on top of a solution, change core: */ - if(control_analysis)solutioncore=&control_core; + if(control_analysis){ + #ifdef _HAVE_TAO_ + solutioncore=&controltao_core; + #else + solutioncore=&control_core; + #endif + } /*are we running the solution sequence, or a qmu wrapper around it? : */ Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 8278) +++ /issm/trunk/src/c/solutions/solutions.h (revision 8279) @@ -24,4 +24,5 @@ void bedslope_core(FemModel* femmodel); void control_core(FemModel* femmodel); +void controltao_core(FemModel* femmodel); void prognostic_core(FemModel* femmodel); void balancedthickness_core(FemModel* femmodel); Index: /issm/trunk/src/c/toolkits/tao/taoincludes.h =================================================================== --- /issm/trunk/src/c/toolkits/tao/taoincludes.h (revision 8279) +++ /issm/trunk/src/c/toolkits/tao/taoincludes.h (revision 8279) @@ -0,0 +1,17 @@ +/* \file toolkits.h + * \brief: this API allows use of external packages, provides patches, etc ... + */ +#ifdef HAVE_CONFIG_H +#include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#ifndef _TAOINCLUDES_H_ +#define _TOOINCLUDES_H_ + +#ifdef _HAVE_TAO_ +#include "tao.h" +#endif + +#endif Index: /issm/trunk/src/c/toolkits/toolkits.h =================================================================== --- /issm/trunk/src/c/toolkits/toolkits.h (revision 8278) +++ /issm/trunk/src/c/toolkits/toolkits.h (revision 8279) @@ -10,4 +10,5 @@ #include "./metis/metisincludes.h" #include "./triangle/triangleincludes.h" +#include "./tao/taoincludes.h" #endif