Changeset 11275

Timestamp:

01/31/12 16:58:36 (13 years ago)

Author:

Mathieu Morlighem

Message:

TAO is now optional depending on a model flag (md.inversion.tao), to avoid having to recompile all the time

Location:

issm/trunk-jpl/src/c

Files:

• EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• modules/StringToEnumx/StringToEnumx.cpp (modified) (1 diff)
• solutions/controltao_core.cpp (modified) (6 diffs)
• solutions/issm.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

@@ -71 +71 @@
         InversionGradientScalingEnum,
         InversionIscontrolEnum,
+        InversionTaoEnum,
         InversionMaxParametersEnum,
         InversionMaxiterPerStepEnum,

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -75 +75 @@
                 case InversionGradientScalingEnum : return "InversionGradientScaling";
                 case InversionIscontrolEnum : return "InversionIscontrol";
+                case InversionTaoEnum : return "InversionTao";
                 case InversionMaxParametersEnum : return "InversionMaxParameters";
                 case InversionMaxiterPerStepEnum : return "InversionMaxiterPerStep";

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -84 +84 @@
         parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
         parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum));
 
         /*some parameters that did not come with the iomodel: */

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -73 +73 @@
         else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
         else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
+        else if (strcmp(name,"InversionTao")==0) return InversionTaoEnum;
         else if (strcmp(name,"InversionMaxParameters")==0) return InversionMaxParametersEnum;
         else if (strcmp(name,"InversionMaxiterPerStep")==0) return InversionMaxiterPerStepEnum;

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*);
 typedef struct {
@@ -39 +40 @@
 
         /*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","8"); if(info) _error_("STOP");
+        //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","40"); if(info) _error_("STOP");
-
-
+        info = PetscOptionsSetValue("-tao_max_funcs","30"); if(info) _error_("STOP");
 
         /*Additional options*/
@@ -58 +57 @@
         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_blmvm"); if(info) _error_("STOP");
+        info = TaoSetType(tao,"tao_cg"); if(info) _error_("STOP");
         info = TaoSetInitialVector(tao,initial_solution);  if(info) _error_("STOP");
+        info = TaoSetObjectiveRoutine(tao,FormFunction,(void*)&user);  if(info) _error_("STOP");
         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);
-        info = TaoSolve(tao); //if(info) _error_("STOP");
+        info = TaoSetTolerances(tao,10e-18,10e-18,10e-18,10e-18,10e-18); if(info) _error_("STOP");
+        info = TaoSolve(tao); if(info) _error_("STOP");
 
         /*Get solution status*/
@@ -96 +97 @@
 }
 
+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){
 
@@ -102 +115 @@
         FemModel *femmodel = user->femmodel;
         Vec       gradient = NULL;
-
-        /*Temp*/
-//      double*   Xserial=NULL;
-//      VecToMPISerial(&Xserial,X);
-//      printf("X= [%20.20g %20.20g %20.20g]\n",Xserial[0],Xserial[1],Xserial[2]);
-        /*End Temp*/
-
 
         int costfunction=SurfaceAbsVelMisfitEnum;
@@ -120 +126 @@
         CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 
-        //printf("X\n");
-        //VecView(X,PETSC_VIEWER_STDOUT_SELF);
-
-        //printf("Gradient\n");
-        //VecView(G,PETSC_VIEWER_STDOUT_SELF);
-
-        printf("f(x) = %g\n",*fcn);
+        printf("f2(x) = %g\n",*fcn);
         return 0;
 }

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

@@ -16 +16 @@
         FILE *petscoptionsfid  = NULL;
         bool  waitonlock       = false;
-        bool  dakota_analysis;
-        bool  control_analysis; 
+        bool  dakota_analysis,control_analysis,tao_analysis;
 
         /*FemModel: */
@@ -81 +80 @@
         femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
         femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
+        femmodel->parameters->FindParam(&tao_analysis,InversionTaoEnum);
         MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
-
-        /*are we running the solution sequence, or a qmu wrapper around it? : */
 
         _printf_(true,"call computational core:\n");
@@ -96 +94 @@
         else if(control_analysis){
                 #ifdef _HAVE_CONTROL_
-                #ifdef _HAVE_TAO_
-                controltao_core(femmodel);
-                #else
-                control_core(femmodel);
-                #endif
+                if(tao_analysis)
+                 controltao_core(femmodel);
+                else
+                 control_core(femmodel);
                 #else
                 _error_("ISSM was not compiled with control support, cannot carry out dakota analysis!");