Changeset 8390


Ignore:
Timestamp:
05/22/11 18:58:03 (14 years ago)
Author:
Mathieu Morlighem
Message:

trunk: improved TAO interface but still not working correctly

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/solutions/controltao_core.cpp

    r8385 r8390  
    33 */
    44#include "../issm.h"
     5
    56#ifdef _HAVE_TAO_
    67#include "tao.h"
    7 #endif
     8
     9/*Local prototype*/
     10int FormFunctionGradient(TAO_APPLICATION,Vec,double *,Vec,void*);
     11typedef struct {
     12        FemModel* femmodel;
     13} AppCtx;
    814
    915void controltao_core(FemModel* femmodel){
    10         #ifdef _HAVE_TAO_ //only works if tao library has been compiled in.
    1116
    12         int     i,n;
    13        
    14         /*parameters: */
    15         int     num_controls;
    16         int     nsteps;
    17         double  eps_cm;
    18         double  tolx;
    19         bool    cm_gradient;
    20         int     dim;
    21         int     solution_type;
    22         bool    isstokes;
    23         bool    qmu_analysis=false;
    24 
    25         int*    control_type = NULL;
    26         double* responses=NULL;
    27         double* maxiter=NULL;
    28         double* cm_jump=NULL;
    29 
    30                
    31         /*intermediary: */
    32         double  search_scalar=0;
    33         OptArgs optargs;
    34         OptPars optpars;
    35 
    36         /*Solution and Adjoint core pointer*/
    37         void (*solutioncore)(FemModel*)=NULL;
    38         void (*adjointcore)(FemModel*)=NULL;
    39 
    40         /*output: */
    41         double* J=NULL;
     17        /*TAO*/
     18        int                i,n,info;
     19        TaoMethod          method     = "tao_blmvm";
     20        //TaoMethod          method     = "tao_lmvm";
     21        //TaoMethod          method     = "tao_cg";
     22        //TaoMethod          method     = "tao_gpcg"; -> Hessian
     23        TaoTerminateReason reason;
     24        TAO_SOLVER         tao;
     25        TAO_APPLICATION    controlapp;
     26        Vec                initial_solution=NULL;
     27        AppCtx          user;                /* user-defined work context */
    4228
    4329#ifdef _HAVE_TAO_
    44         int argc;
    45         char **args;
    46         PetscGetArgs(&argc,&args);
     30        int argc; char **args; PetscGetArgs(&argc,&args);
    4731        int ierr=TaoInitialize(&argc,&args,(char*)0,"");
    4832        if(ierr) _error_("Could not initialize Tao");
    4933#endif
    5034
    51         /*Recover parameters used throughout the solution:{{{1*/
    52         femmodel->parameters->FindParam(&nsteps,NStepsEnum);
    53         femmodel->parameters->FindParam(&num_controls,NumControlsEnum);
    54         femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    55         femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    56         femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum);
    57         femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum);
    58         femmodel->parameters->FindParam(&eps_cm,EpsCmEnum);
    59         femmodel->parameters->FindParam(&tolx,TolXEnum);
    60         femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
    61         femmodel->parameters->FindParam(&dim,DimEnum);
    62         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    63         femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    64         femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
    65         /*}}}*/
     35        /*Initialize argument*/
     36        user.femmodel=femmodel;
    6637
    67         /*out of solution_type, figure out solution core and adjoint function pointer*/
    68         CorePointerFromSolutionEnum(&solutioncore,solution_type);
    69         AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
     38        info = TaoCreate(PETSC_COMM_WORLD,method,&tao); if(info) _error_("STOP");
     39        info = TaoApplicationCreate(PETSC_COMM_WORLD,&controlapp); if(info) _error_("STOP");
     40        GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,RheologyBbarEnum,VertexEnum);
     41        info = TaoAppSetInitialSolutionVec(controlapp,initial_solution);  if(info) _error_("STOP");
     42        info = TaoAppSetObjectiveAndGradientRoutine(controlapp,FormFunctionGradient,(void*)&user);  if(info) _error_("STOP");
     43        info = TaoSetOptions(controlapp,tao);  if(info) _error_("STOP");
     44        info = TaoSetTolerances(tao,1e-24,1e-24,1e-24,1e-24); if(info) _error_("STOP");
     45        //info = TaoSetFunctionLowerBound(tao,0.77); if(info) _error_("STOP");
     46        info = TaoSolveApplication(controlapp,tao); //if(info) _error_("STOP");
     47        PetscInt          iter;
     48        double          ff,gnorm;
     49        info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); //CHKERRQ(info);
    7050
    71         /*Launch once a complete solution to set up all inputs*/
    72         _printf_(VerboseControl(),"%s\n","   preparing initial solution");
    73         if (isstokes) solutioncore(femmodel);
     51        switch(reason){
     52                /*http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoGetTerminationReason.html*/
     53                case TAO_CONVERGED_ATOL:       _printf_(true,"TAO_CONVERGED_ATOL (res <= atol)\n"); break;
     54                case TAO_CONVERGED_RTOL:       _printf_(true,"TAO_CONVERGED_RTOL (res/res0 <= rtol)\n"); break;
     55                case TAO_CONVERGED_TRTOL:      _printf_(true,"TAO_CONVERGED_TRTOL (xdiff <= trtol) \n"); break;
     56                case TAO_CONVERGED_MINF:       _printf_(true,"TAO_CONVERGED_MINF  (f <= fmin)\n"); break;
     57                case TAO_CONVERGED_USER:       _printf_(true,"TAO_CONVERGED_USER (user defined)\n"); break;
     58                case TAO_DIVERGED_MAXITS:      _printf_(true,"TAO_DIVERGED_MAXITS (its>maxits)\n"); break;
     59                case TAO_DIVERGED_NAN:         _printf_(true,"TAO_DIVERGED_NAN (Numerical problems)\n"); break;
     60                case TAO_DIVERGED_MAXFCN:      _printf_(true,"TAO_DIVERGED_MAXFCN (nfunc > maxnfuncts)\n"); break;
     61                case TAO_DIVERGED_LS_FAILURE:  _printf_(true,"TAO_DIVERGED_LS_FAILURE (line search failure)\n"); break;
     62                case TAO_DIVERGED_TR_REDUCTION:_printf_(true,"TAO_DIVERGED_TR_REDUCTION \n"); break;
     63                case TAO_DIVERGED_USER:        _printf_(true,"TAO_DIVERGED_USER (user defined)\n"); break;
     64                default: _printf_(true,"unknown reason");
     65        }
     66        if (reason<=0){
     67                _printf_(true,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
     68                _printf_(true," Iterations: %d,  Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
     69        }
     70        else{
     71         _printf_(true,"TAO found a solution");
     72        }
     73        info = TaoView(tao);  if(info) _error_("STOP");
     74        info = TaoDestroy(tao);  if(info) _error_("STOP");
     75        info = TaoAppDestroy(controlapp);  if(info) _error_("STOP");
    7476
    75         /*Initialize responses: */
    76         J=(double*)xmalloc(nsteps*sizeof(double));
    77                
    78         /*Initialize some of the BrentSearch arguments: */
    79         optargs.femmodel=femmodel;
    80         optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ;
    81        
    82         /*Start looping: */
    83         for(n=0;n<nsteps;n++){
    84 
    85                 /*Display info*/
    86                 _printf_(VerboseControl(),"\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    87                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n],CmResponseEnum);
    88                
    89                 /*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */
    90                 if (solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
    91 
    92                 _printf_(VerboseControl(),"%s\n","   compute adjoint state:");
    93                 adjointcore(femmodel);
    94        
    95                 gradient_core(femmodel,n,search_scalar);
    96 
    97                 /*Return gradient if asked: */
    98                 if (cm_gradient){
    99                         InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum);
    100                         goto cleanup_and_return;
    101                 }
    102 
    103                 _printf_(VerboseControl(),"%s\n","   optimizing along gradient direction");
    104                 optargs.n=n; optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
    105                 BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
    106                 //OptimalSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
    107 
    108                 _printf_(VerboseControl(),"%s\n","   updating parameter using optimized search scalar"); //true means update save controls
    109                 InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
    110                
    111                 if(controlconvergence(J,responses,eps_cm,n)) break;
    112 
    113                 /*Temporary saving every 5 control steps: */
    114                 /*if (((n+1)%5)==0){
    115                         _printf_(VerboseControl(),"%s\n","   saving temporary results");
    116                         controlrestart(femmodel,J);
    117                 }*/
    118         }
    119 
    120         _printf_(VerboseControl(),"%s\n","   preparing final solution");
    121         femmodel->parameters->SetParam(false,ControlAnalysisEnum); //needed to turn control result output in solutioncore
    122         solutioncore(femmodel);
    123 
    124         /*some results not computed by steadystate_core or diagnostic_core: */
    125         if(!qmu_analysis){ //do not save this if we are running the control core from a qmu run!
    126                 for(i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);
    127                 femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
    128                 //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToStringx(control_type),1,0));
    129         }
    130 
    131         cleanup_and_return:
    132         /*Free ressources: */
    133         xfree((void**)&control_type);
    134         xfree((void**)&responses);
    135         xfree((void**)&maxiter);
    136         xfree((void**)&cm_jump);
    137         xfree((void**)&J);
    138        
    139         /*control_core might be used in Qmu, so leave everything similar to where it started: */
    140         femmodel->parameters->SetParam(true,ControlAnalysisEnum);
     77        /*Clean up*/
     78        VecFree(&initial_solution);
    14179
    14280        /* Finalize TAO */
    14381        TaoFinalize();
     82}
    14483
     84int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn,Vec G,void *userCtx){
     85
     86        /*Retreive arguments*/
     87        AppCtx   *user     = (AppCtx*)userCtx;
     88        FemModel *femmodel = user->femmodel;
     89        Vec       gradient = NULL;
     90
     91        InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum,CmResponseEnum);
     92        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,RheologyBbarEnum,VertexEnum);
     93        adjointdiagnostic_core(user->femmodel);
     94        Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, RheologyBbarEnum);
     95        //VecView(gradient,PETSC_VIEWER_STDOUT_SELF);
     96        //VecScale(gradient,-1.0);
     97        VecCopy(gradient,G);
     98        //VecView(G,PETSC_VIEWER_STDOUT_SELF);
     99        CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum);
     100        printf("f(x) = %g\n",*fcn);
     101        return 0;
     102}
     103#else
     104void controltao_core(FemModel* femmodel){
     105        _error_("TAO not installed");
     106}
    145107#endif //_HAVE_TAO_
    146 }
Note: See TracChangeset for help on using the changeset viewer.