Changeset 8390
- Timestamp:
- 05/22/11 18:58:03 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/solutions/controltao_core.cpp
r8385 r8390 3 3 */ 4 4 #include "../issm.h" 5 5 6 #ifdef _HAVE_TAO_ 6 7 #include "tao.h" 7 #endif 8 9 /*Local prototype*/ 10 int FormFunctionGradient(TAO_APPLICATION,Vec,double *,Vec,void*); 11 typedef struct { 12 FemModel* femmodel; 13 } AppCtx; 8 14 9 15 void controltao_core(FemModel* femmodel){ 10 #ifdef _HAVE_TAO_ //only works if tao library has been compiled in.11 16 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 */ 42 28 43 29 #ifdef _HAVE_TAO_ 44 int argc; 45 char **args; 46 PetscGetArgs(&argc,&args); 30 int argc; char **args; PetscGetArgs(&argc,&args); 47 31 int ierr=TaoInitialize(&argc,&args,(char*)0,""); 48 32 if(ierr) _error_("Could not initialize Tao"); 49 33 #endif 50 34 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; 66 37 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); 70 50 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"); 74 76 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); 141 79 142 80 /* Finalize TAO */ 143 81 TaoFinalize(); 82 } 144 83 84 int 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 104 void controltao_core(FemModel* femmodel){ 105 _error_("TAO not installed"); 106 } 145 107 #endif //_HAVE_TAO_ 146 }
Note:
See TracChangeset
for help on using the changeset viewer.