Changeset 11317 for issm/trunk-jpl/src/c/solutions/controltao_core.cpp

Timestamp:

02/03/12 13:50:06 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added 2 modules to make TAO work

File:

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

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

@@ -25 +25 @@
 
         /*TAO*/
-        int       ierr,numberofvertices;
-        AppCtx    user;
-        TaoSolver tao;
-        Vec       X  = NULL;
-        Vec       XL = NULL;
-        Vec       XU = NULL;
+        int        ierr,numberofvertices;
+        int        num_controls;
+        AppCtx     user;
+        TaoSolver  tao;
+        int       *control_list   = NULL;
+        Vec        X              = NULL;
+        Vec        XL             = NULL;
+        Vec        XU             = NULL;
 
         /*Initialize TAO*/
@@ -37 +39 @@
         ierr = TaoInitialize(&argc,&args,(char*)0,"");
         if(ierr) _error_("Could not initialize Tao");
+
+        /*Recover some parameters*/
+        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+        femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum);
 
         /*Initialize TAO*/
@@ -52 +58 @@
         XL=NewVec(numberofvertices); VecSet(XL,1.);
         XU=NewVec(numberofvertices); VecSet(XU,200.);
-        TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);
+        //TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);
 
-        GetVectorFromInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum,VertexEnum);
+        GetVectorFromControlInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
         TaoSetInitialVector(tao,X);
 
@@ -64 +70 @@
         TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
         TaoGetSolutionVector(tao,&X);
-        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
-        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
+        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+        for(int i=0;i<num_controls;i++){
+                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]);
+        }
 
         /*Clean up and return*/
+        xfree((void**)&control_list);
         VecFree(&X);
         TaoDestroy(&tao);
@@ -75 +84 @@
 
         /*Retreive arguments*/
-        AppCtx   *user     = (AppCtx*)userCtx;
-        FemModel *femmodel = user->femmodel;
-        Vec       gradient = NULL;
+        int       solution_type,num_cost_functions;
+        AppCtx   *user           = (AppCtx *)userCtx;
+        FemModel *femmodel       = user->femmodel;
+        int      *cost_functions = NULL;
+        double   *cost_functionsd= NULL;
+        Vec       gradient       = NULL;
 
-        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,FrictionCoefficientEnum,VertexEnum);
-        adjointdiagnostic_core(user->femmodel);
-        //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum);
+        /*Set new variable*/
+        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+
+        /*Recover some parameters*/
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+        femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
+        femmodel->parameters->FindParam(&cost_functionsd,NULL,NULL,InversionCostFunctionsEnum);
+
+        /*Prepare objective function*/
+        cost_functions=(int*)xmalloc(num_cost_functions*sizeof(int));
+        for(int i=0;i<num_cost_functions;i++) cost_functions[i]=(int)cost_functionsd[i]; //FIXME
+        femmodel->parameters->SetParam(cost_functions,1,num_cost_functions,StepResponsesEnum);
+
+        /*Compute solution and adjoint*/
+        void (*solutioncore)(FemModel*)=NULL;
+        void (*adjointcore)(FemModel*)=NULL;
+        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+        solutioncore(femmodel);
+        adjointcore(femmodel);
+
+        /*Compute objective function*/
+        CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+        /*Compute gradient*/
         Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-        VecScale(gradient,-1.);
-        VecCopy(gradient,G);
-        VecFree(&gradient);
-        CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+        VecCopy(gradient,G); VecFree(&gradient);
+        VecScale(G,-1.);
+
+        /*Clean-up and return*/
+        xfree((void**)&cost_functions);
+        xfree((void**)&cost_functionsd);
         return 0;
 }