Changeset 11283

Timestamp:

02/01/12 11:04:08 (13 years ago)

Author:

Mathieu Morlighem

Message:

Fixed tao, now working smoothly

File:

• issm/trunk-jpl/src/c/solutions/controltao_core.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/solutions/controltao_core.cpp

@@ -17 +17 @@
 
 /*Local prototype*/
-int FormFunction(TaoSolver tao,Vec,double *,void*);
-int FormFunctionGradient(TaoSolver tao,Vec,double *,Vec,void*);
+int FormFunctionGradient(TaoSolver tao,Vec,double*,Vec,void*);
 typedef struct {
         FemModel* femmodel;
@@ -26 +25 @@
 
         /*TAO*/
-        int                        i,n,info;
-        TaoSolverTerminationReason reason;
-        TaoSolver                  tao;
-        Vec                        initial_solution = NULL;
-        AppCtx                     user;
-        PetscInt                   iter;
-        double                     ff                 ,gnorm;
+        int       i,n,info;
+        TaoSolver tao;
+        Vec       initial_solution = NULL;
+        AppCtx    user;
 
         /*Initialize TAO*/
@@ -40 +36 @@
 
         /*Line search options*/
-        //info = PetscOptionsSetValue("-tao_ls_stepmax","10e11"); if(info) _error_("STOP"); //does not work
-        //info = PetscOptionsSetValue("-tao_ls_stepmin","10e5"); if(info) _error_("STOP");    //does not work
-        info = PetscOptionsSetValue("-tao_ls_maxfev","3"); if(info) _error_("STOP");
         /*TAO options: http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoSetFromOptions.html*/
         info = PetscOptionsSetValue("-tao_monitor",""); if(info) _error_("STOP");
-        info = PetscOptionsSetValue("-tao_max_its","10"); if(info) _error_("STOP");
-        info = PetscOptionsSetValue("-tao_max_funcs","30"); if(info) _error_("STOP");
-
-        /*Additional options*/
-        //info = PetscOptionsSetValue("-info","/u/astrid-r1b/morlighe/svn/issm/trunk/test/NightlyRun/taolog.txt"); if(info) _error_("STOP");
+        info = PetscOptionsSetValue("-tao_max_its","20"); if(info) _error_("STOP");
+        info = PetscOptionsSetValue("-tao_max_funcs","50"); if(info) _error_("STOP");
 
         /*Initialize argument*/
@@ -55 +45 @@
 
         /*Set up and solve TAO*/
-        GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,MaterialsRheologyBbarEnum,VertexEnum);
         info = TaoCreate(PETSC_COMM_WORLD,&tao); if(info) _error_("STOP");
-        //info = TaoSetType(tao,"tao_blmvm"); if(info) _error_("STOP");
-        info = TaoSetType(tao,"tao_cg"); if(info) _error_("STOP");
+        info = TaoSetType(tao,"tao_blmvm"); if(info) _error_("STOP");
+        //info = TaoSetType(tao,"tao_cg"); if(info) _error_("STOP");
+
+        int size=femmodel->vertices->NumberOfVertices();
+        Vec XL=NULL;
+        Vec XU=NULL;
+        VecCreate(PETSC_COMM_WORLD,&XL);                               VecCreate(PETSC_COMM_WORLD,&XU);
+        VecSetSizes(XL,PetscDetermineLocalSize(size),PETSC_DECIDE);    VecSetSizes(XU,PetscDetermineLocalSize(size),PETSC_DECIDE);
+        VecSetFromOptions(XL);                                         VecSetFromOptions(XU);
+        VecSet(XL,1);                                                  VecSet(XU,200.);
+        TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);
+
+        GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,FrictionCoefficientEnum,VertexEnum);
         info = TaoSetInitialVector(tao,initial_solution);  if(info) _error_("STOP");
-        info = TaoSetObjectiveRoutine(tao,FormFunction,(void*)&user);  if(info) _error_("STOP");
+        VecFree(&initial_solution);
+
         info = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);  if(info) _error_("STOP");
         info = TaoSetFromOptions(tao);  if(info) _error_("STOP");
         /* http://www.mcs.anl.gov/research/projects/tao/docs/manpages/taosolver/TaoSetTolerances.html*/
-        /*                          fatol, frtol, gatol, grtol, gttol*/
-        info = TaoSetTolerances(tao,10e-18,10e-18,10e-18,10e-18,10e-18); if(info) _error_("STOP");
+        /*                          fatol ,frtol ,gatol ,grtol ,gttol*/
+        info = TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28); if(info) _error_("STOP");
         info = TaoSolve(tao); if(info) _error_("STOP");
 
         /*Get solution status*/
-        info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); //CHKERRQ(info);
-        switch(reason){
-                /*http://www.mcs.anl.gov/research/projects/tao/docs/manpages/taosolver/TaoGetTerminationReason.html*/
-                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,PETSC_VIEWER_STDOUT_SELF);  if(info) _error_("STOP");
+        info = TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);  if(info) _error_("STOP");
 
         /*Clean up*/
         info = TaoDestroy(&tao);  if(info) _error_("STOP");
-        VecFree(&initial_solution);
+        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
 
         /* Finalize TAO */
         TaoFinalize();
-}
-
-int FormFunction(TaoSolver tao, Vec X, double *fcn,void *userCtx){
-        AppCtx   *user     = (AppCtx*)userCtx;
-        FemModel *femmodel = user->femmodel;
-
-        int costfunction=SurfaceAbsVelMisfitEnum;
-        femmodel->parameters->SetParam(&costfunction,1,1,StepResponsesEnum);
-        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum);
-        diagnostic_core(user->femmodel);
-        CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
-        printf("f(x) = %g\n",*fcn);
-        return 0;
 }
 int FormFunctionGradient(TaoSolver tao, Vec X, double *fcn,Vec G,void *userCtx){
@@ -118 +88 @@
         int costfunction=SurfaceAbsVelMisfitEnum;
         femmodel->parameters->SetParam(&costfunction,1,1,StepResponsesEnum);
-        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum);
+        //InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum);
+        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
         adjointdiagnostic_core(user->femmodel);
-        Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, MaterialsRheologyBbarEnum);
-        VecScale(gradient,-10e7);
-        //VecScale(gradient,-1.);
+        //Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, MaterialsRheologyBbarEnum);
+        Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
+        VecScale(gradient,-1.);
         VecCopy(gradient,G);
+        VecFree(&gradient);
         CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
-
-        printf("f2(x) = %g\n",*fcn);
         return 0;
 }
+
 #else
 void controltao_core(FemModel* femmodel){