Changeset 4048 for issm/trunk/src/c/solutions/control_core.cpp

Timestamp:

06/16/10 18:20:52 (15 years ago)

Author:

Eric.Larour

Message:

Serious cleanup of the control solutoin

File:

• issm/trunk/src/c/solutions/control_core.cpp (modified) (4 diffs)

issm/trunk/src/c/solutions/control_core.cpp

@@ -14 +14 @@
 void control_core(FemModel* femmodel){
 
-        Vec     u_g=NULL;
-        Vec     t_g=NULL;
-        Vec     m_g=NULL;
-        double  search_scalar;
+        int     i,n;
+        
+        /*parameters: */
+        int     verbose=0;
         int     control_type;
+        int     control_steady;
+        int     nsteps;
+        double  eps_cm;
+        double  tolx;
+        double  cm_min;
+        double  cm_max;
+        int     cm_gradient;
+
         double* fit=NULL;
         double* optscal=NULL;
         double* maxiter=NULL;
         double* cm_jump=NULL;
-        double  eps_cm;
-        double  tolx;
-        double* param_g=NULL;
-        Vec     grad_g=NULL;
-        Vec     new_grad_g=NULL;
-        Vec     grad_g_old=NULL;
-        double* grad_g_double=NULL;
-        double  cm_min;
-        double  cm_max;
-        int     cm_gradient;
-        int     nsteps,n,i;
-        double* J=NULL;
+
+                
+        /*intermediary: */
+        double  search_scalar;
         OptArgs optargs;
         OptPars optpars;
-        Param*  param=NULL;
 
-        /*flags: */
-        int analysis_type;
-        int sub_analysis_type;
-        int     control_steady;
-        int verbose=0;
-        int converged=0;
-        int numberofnodes;
+        /*output: */
+        double* J=NULL;
 
         /*Process models*/
@@ -62 +56 @@
         femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
         femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
-        femmodel->parameters->FindParam(&param_g,NULL,NULL,ControlParameterEnum);
-        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
-        femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
         femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
 
         /*Initialize misfit: */
         J=(double*)xmalloc(nsteps*sizeof(double));
-
+                
+        /*Initialize some of the BrentSearch arguments: */
+        optargs.femmodel=femmodel;
+        optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ;
+        
         /*Start looping: */
         for(n=0;n<nsteps;n++){
 
                 _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
-                inputs->Add(control_type,param_g,1,numberofnodes);
                 femmodel->UpdateInputsFromConstant(fit[n],FitEnum);
                 
-                /*In case we are running a steady state control method, compute new temperature field using new parameter 
-                 * distribution: */
+                /*In case we are running a steady state control method, compute new temperature field using new parameter * distribution: */
                 if (control_steady) steadystate_core(model);
         
@@ -87 +79 @@
                 /*Return gradient if asked: */
                 if (cm_gradient){
-                        /*Transfer gradient from input to results: */
                         InputToResultx(femodel->elements,femodel->nodes,femodel->vertices,femodel->loads,femodel->materials,femodel->parameters,GradientEnum);
-                        return;
+                        goto cleanup_and_return;
                 }
 
-                /*Normalize if last gradient not satisfying (search_scalar==0)*/
-                if (n>0 && search_scalar==0){
-                        _printf_("%s","      orthogonalization...");
-                        Orthx(&new_grad_g,grad_g,grad_g_old);
-                        _printf_("%s\n"," done.");
-                }
-                else{ 
-                        _printf_("%s","      normalizing directions...");
-                        Orthx(&new_grad_g,grad_g,NULL);
-                        _printf_("%s\n"," done.");
-                }
-                VecFree(&grad_g); VecFree(&grad_g_old); 
-                grad_g_old=new_grad_g;
-                VecToMPISerial(&grad_g_double,new_grad_g);
-
                 _printf_("%s\n","      optimizing along gradient direction");
-                optargs.model=model;
-                optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.n=n;
-                optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];optpars.cm_jump=cm_jump[n];
+                optargs.n=n; optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
                 BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
 
                 _printf_("%s","      updating parameter using optimized search scalar...");
-                for(i=0;i<numberofnodes;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
-                _printf_("%s\n"," done.");
+                AXPYInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],ControlParameterEnum);
 
                 _printf_("%s","      constraining the new distribution...");    
-                ControlConstrainx(param_g,numberofnodes,cm_min,cm_max,control_type);
-                _printf_("%s\n"," done.");
+                ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
+                
+                _printf_("%s","      save new parameter...");
+                DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum);
                 
                 _printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
+                
+                if(controlconvergence)goto cleanup_and_return;
 
-                /*some freeing:*/
-                xfree((void**)&grad_g_double);
-
-                /*Has convergence been reached?*/
-                if (!isnan(eps_cm)){
-                        i=n-2;
-                        //go through the previous misfits(starting from n-2)
-                        while(i>=0){
-                                if (fit[i]==fit[n]){
-                                        //convergence test only if we have the same misfits
-                                        if ((J[i]-J[n])/J[n] <= eps_cm){
-                                                //convergence if convergence criteria fullfilled
-                                                converged=1;
-                                                _printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],"<",eps_cm);
-                                        }
-                                        else{
-                                                _printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],">",eps_cm);
-                                        }
-                                        break;
-                                }
-                                i=i-1;
-                        }
-                }
-                //stop if convergence has been reached
-                if(converged) break;
-
-                //some temporary saving
+                /*Temporary saving every 5 control steps: */
                 if (((n+1)%5)==0){
                         _printf_("%s","      saving temporary results...");
-                        ControlRestart(model,param_g);
-                        _printf_("%s\n"," done.");
+                        ControlRestart(femmodel);
                 }
         }
 
-        /*Write results to disk: */
         _printf_("%s","      preparing final velocity solution");
-        /*Launch diagnostic with the last parameter distribution*/
-        if (control_steady){
-                model->UpdateInputsFromVector(param_g,control_type,VertexEnum);
-                steadystate_results=steadystate_core(model);
+        if (control_steady) steadystate_core(femmodel);
+        else diagnostic_core(femmodel);
 
-                //extract u_g ,t_g and m_g from steadystate results, and erase diagnostic_results;
-                steadystate_results->FindResult(&u_g,"u_g");
-                steadystate_results->FindResult(&m_g,"m_g");
-                steadystate_results->FindResult(&t_g,"t_g");
-                delete steadystate_results;
-        }
-        else{
-                model->UpdateInputsFromVector(param_g,control_type,VertexEnum);
-                diagnostic_results=diagnostic_core(model);
+        /*some results not computed by steadystate_core or diagnostic_core: */
+        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
+        femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
+        femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
 
-                //extract u_g from diagnostic_results, and erase diagnostic_results;
-                diagnostic_results->FindResult(&u_g,"u_g");
-                delete diagnostic_results;
-        }
-
-        /*Plug results into output dataset: */
-        results->AddObject(new Result(results->Size()+1,0,1,"u_g",u_g));
-        results->AddObject(new Result(results->Size()+1,0,1,"param_g",param_g,numberofnodes));
-        results->AddObject(new Result(results->Size()+1,0,1,"J",J,nsteps));
-        if (control_steady){
-                results->AddObject(new Result(results->Size()+1,0,1,"t_g",t_g));
-                results->AddObject(new Result(results->Size()+1,0,1,"m_g",m_g));
-        }
+        cleanup_and_return:
         
-        /*Add analysis_type and control_type to results: */
-        results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(analysis_type)));
-        results->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
-
         /*Free ressources: */
         xfree((void**)&control_type);
@@ -198 +124 @@
         xfree((void**)&maxiter);
         xfree((void**)&cm_jump);
-        VecFree(&new_grad_g); //do not VecFree grad_g and grad_g_old, they point to new_grad_g
-        xfree((void**)&grad_g_double);
-        xfree((void**)&param_g);
-        VecFree(&u_g);
-        VecFree(&t_g);
-        VecFree(&m_g);
         xfree((void**)&J);
-
-        //return: 
-        return results;
 }