Index: /issm/trunk/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/controltao_core.cpp (revision 8389) +++ /issm/trunk/src/c/solutions/controltao_core.cpp (revision 8390) @@ -3,144 +3,105 @@ */ #include "../issm.h" + #ifdef _HAVE_TAO_ #include "tao.h" -#endif + +/*Local prototype*/ +int FormFunctionGradient(TAO_APPLICATION,Vec,double *,Vec,void*); +typedef struct { + FemModel* femmodel; +} AppCtx; void controltao_core(FemModel* femmodel){ - #ifdef _HAVE_TAO_ //only works if tao library has been compiled in. - 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; + /*TAO*/ + int i,n,info; + TaoMethod method = "tao_blmvm"; + //TaoMethod method = "tao_lmvm"; + //TaoMethod method = "tao_cg"; + //TaoMethod method = "tao_gpcg"; -> Hessian + TaoTerminateReason reason; + TAO_SOLVER tao; + TAO_APPLICATION controlapp; + Vec initial_solution=NULL; + AppCtx user; /* user-defined work context */ #ifdef _HAVE_TAO_ - int argc; - char **args; - PetscGetArgs(&argc,&args); + int argc; char **args; PetscGetArgs(&argc,&args); int ierr=TaoInitialize(&argc,&args,(char*)0,""); if(ierr) _error_("Could not initialize Tao"); #endif - /*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); - /*}}}*/ + /*Initialize argument*/ + user.femmodel=femmodel; - /*out of solution_type, figure out solution core and adjoint function pointer*/ - CorePointerFromSolutionEnum(&solutioncore,solution_type); - AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); + info = TaoCreate(PETSC_COMM_WORLD,method,&tao); if(info) _error_("STOP"); + info = TaoApplicationCreate(PETSC_COMM_WORLD,&controlapp); if(info) _error_("STOP"); + GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,RheologyBbarEnum,VertexEnum); + info = TaoAppSetInitialSolutionVec(controlapp,initial_solution); if(info) _error_("STOP"); + info = TaoAppSetObjectiveAndGradientRoutine(controlapp,FormFunctionGradient,(void*)&user); if(info) _error_("STOP"); + info = TaoSetOptions(controlapp,tao); if(info) _error_("STOP"); + info = TaoSetTolerances(tao,1e-24,1e-24,1e-24,1e-24); if(info) _error_("STOP"); + //info = TaoSetFunctionLowerBound(tao,0.77); if(info) _error_("STOP"); + info = TaoSolveApplication(controlapp,tao); //if(info) _error_("STOP"); + PetscInt iter; + double ff,gnorm; + info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); //CHKERRQ(info); - /*Launch once a complete solution to set up all inputs*/ - _printf_(VerboseControl(),"%s\n"," preparing initial solution"); - if (isstokes) solutioncore(femmodel); + switch(reason){ + /*http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoGetTerminationReason.html*/ + case TAO_CONVERGED_ATOL: _printf_(true,"TAO_CONVERGED_ATOL (res <= atol)\n"); break; + case TAO_CONVERGED_RTOL: _printf_(true,"TAO_CONVERGED_RTOL (res/res0 <= rtol)\n"); break; + case TAO_CONVERGED_TRTOL: _printf_(true,"TAO_CONVERGED_TRTOL (xdiff <= trtol) \n"); break; + case TAO_CONVERGED_MINF: _printf_(true,"TAO_CONVERGED_MINF (f <= fmin)\n"); break; + case TAO_CONVERGED_USER: _printf_(true,"TAO_CONVERGED_USER (user defined)\n"); break; + case TAO_DIVERGED_MAXITS: _printf_(true,"TAO_DIVERGED_MAXITS (its>maxits)\n"); break; + case TAO_DIVERGED_NAN: _printf_(true,"TAO_DIVERGED_NAN (Numerical problems)\n"); break; + case TAO_DIVERGED_MAXFCN: _printf_(true,"TAO_DIVERGED_MAXFCN (nfunc > maxnfuncts)\n"); break; + case TAO_DIVERGED_LS_FAILURE: _printf_(true,"TAO_DIVERGED_LS_FAILURE (line search failure)\n"); break; + case TAO_DIVERGED_TR_REDUCTION:_printf_(true,"TAO_DIVERGED_TR_REDUCTION \n"); break; + case TAO_DIVERGED_USER: _printf_(true,"TAO_DIVERGED_USER (user defined)\n"); break; + default: _printf_(true,"unknown reason"); + } + if (reason<=0){ + _printf_(true,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n"); + _printf_(true," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm); + } + else{ + _printf_(true,"TAO found a solution"); + } + info = TaoView(tao); if(info) _error_("STOP"); + info = TaoDestroy(tao); if(info) _error_("STOP"); + info = TaoAppDestroy(controlapp); if(info) _error_("STOP"); - /*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); + /*Clean up*/ + VecFree(&initial_solution); /* Finalize TAO */ TaoFinalize(); +} +int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn,Vec G,void *userCtx){ + + /*Retreive arguments*/ + AppCtx *user = (AppCtx*)userCtx; + FemModel *femmodel = user->femmodel; + Vec gradient = NULL; + + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum,CmResponseEnum); + InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,RheologyBbarEnum,VertexEnum); + adjointdiagnostic_core(user->femmodel); + Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, RheologyBbarEnum); + //VecView(gradient,PETSC_VIEWER_STDOUT_SELF); + //VecScale(gradient,-1.0); + VecCopy(gradient,G); + //VecView(G,PETSC_VIEWER_STDOUT_SELF); + CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum); + printf("f(x) = %g\n",*fcn); + return 0; +} +#else +void controltao_core(FemModel* femmodel){ + _error_("TAO not installed"); +} #endif //_HAVE_TAO_ -}