Changeset 2892

Timestamp:

01/22/10 07:26:54 (15 years ago)

Author:

Eric.Larour

Message:

Added new solution type: control gradient. Can now retrieve control method gradient
and return, without computing whole control method solution.

Location:

issm/trunk/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/ModelProcessorx/Control/CreateParametersControl.cpp (modified) (2 diffs)
• c/ModelProcessorx/IoModel.cpp (modified) (2 diffs)
• c/ModelProcessorx/IoModel.h (modified) (1 diff)
• c/issm.h (modified) (1 diff)
• c/parallel/ProcessResults.cpp (modified) (2 diffs)
• c/parallel/control_core.cpp (modified) (3 diffs)
• m/classes/@model/model.m (modified) (1 diff)
• m/classes/@model/setdefaultparameters.m (modified) (1 diff)
• mex/Makefile.am (modified) (2 diffs)

issm/trunk/src/c/Makefile.am

@@ -308 +308 @@
                                         ./ProcessParamsx/ProcessParamsx.cpp\
                                         ./ProcessParamsx/ProcessParamsx.h\
+                                        ./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp\
+                                        ./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h\
                                         ./Qmux/Qmux.h\
                                         ./Qmux/Qmux.cpp\

issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp

@@ -20 +20 @@
         Param*   param = NULL;
         int      count;
-        int      analysis_type;
-        int      numberofdofspernode;
 
         double* fit=NULL;
@@ -140 +138 @@
                 param= new Param(count,"cm_max",DOUBLE);
                 param->SetDouble(iomodel->cm_max);
+                parameters->AddObject(param);
+
+                /*cm_gradient: */
+                count++;
+                param= new Param(count,"cm_gradient",DOUBLE);
+                param->SetDouble(iomodel->cm_gradient);
                 parameters->AddObject(param);
 

issm/trunk/src/c/ModelProcessorx/IoModel.cpp

@@ -125 +125 @@
         iomodel->cm_min=0;
         iomodel->cm_max=0;
+        iomodel->cm_gradient=0;
         iomodel->verbose=0;
         iomodel->plot=0;
@@ -353 +354 @@
         IoModelFetchData(&iomodel->cm_min,iomodel_handle,"cm_min");
         IoModelFetchData(&iomodel->cm_max,iomodel_handle,"cm_max");
+        IoModelFetchData(&iomodel->cm_gradient,iomodel_handle,"cm_gradient");
         IoModelFetchData(&iomodel->eps_res,iomodel_handle,"eps_res");
         IoModelFetchData(&iomodel->eps_rel,iomodel_handle,"eps_rel");

issm/trunk/src/c/ModelProcessorx/IoModel.h

@@ -117 +117 @@
         double  cm_min;
         double  cm_max;
+        int     cm_gradient;;
 
         double  cm_noisedampening;

issm/trunk/src/c/issm.h

@@ -49 +49 @@
 #include "./PenaltyConstraintsx/PenaltyConstraintsx.h"
 #include "./ProcessParamsx/ProcessParamsx.h"
+#include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
 #include "./Gradjx/Gradjx.h"
 #include "./Orthx/Orthx.h"

issm/trunk/src/c/parallel/ProcessResults.cpp

@@ -90 +90 @@
         double* m_g_serial=NULL;
         double* melting=NULL;
+
+        Vec     grad_g=NULL;
+        double* grad_g_serial=NULL;
+        double* gradient=NULL;
 
         Vec     v_g=NULL;
@@ -286 +290 @@
                         VecFree(&t_g);
                 }
+                else if(strcmp(result->GetFieldName(),"grad_g")==0){
+                        
+                        /*easy, grad_g is of size numberofnodes, on 1 dof, just repartition: */
+                        result->GetField(&grad_g);
+                        VecToMPISerial(&grad_g_serial,grad_g);
+                        fem_c->parameters->FindParam(&numberofnodes,"numberofnodes");
+                        VecToMPISerial(&partition,fem_c->partition->vector);
+
+                        gradient=(double*)xmalloc(numberofnodes*sizeof(double));
+
+                        for(i=0;i<numberofnodes;i++){
+                                gradient[i]=grad_g_serial[(int)partition[i]];
+                        }
+
+                        /*Ok, add gradient to newresults: */
+                        newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"gradient",gradient,numberofnodes);
+                        newresults->AddObject(newresult);
+
+                        /*do some cleanup: */
+                        xfree((void**)&grad_g_serial);
+                        xfree((void**)&partition);
+                        xfree((void**)&gradient);
+                        VecFree(&grad_g);
+                }
                 else if(strcmp(result->GetFieldName(),"m_g")==0){
                         /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */

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

@@ -42 +42 @@
         double  cm_min;
         double  cm_max;
+        int     cm_gradient;
         int     nsteps,n,i;
         double* J=NULL;
@@ -71 +72 @@
         model->FindParam(&cm_min,"cm_min");
         model->FindParam(&cm_max,"cm_max");
+        model->FindParam(&cm_gradient,"cm_gradient");
         model->FindParam(&param_g,NULL,NULL,"param_g");
         model->FindParam(&analysis_type,"analysis_type");
@@ -113 +115 @@
                 gradjcompute_results->FindResult(&grad_g,"grad_g");
                 delete gradjcompute_results;
+
+                /*Return gradient if asked: */
+                if (cm_gradient){
+                        
+                        /*Plug results into output dataset: */
+                        result=new Result(results->Size()+1,0,1,"grad_g",grad_g);
+                        results->AddObject(result);
+                
+                        /*Free ressources: */
+                        xfree((void**)&control_type);
+                        xfree((void**)&fit);
+                        xfree((void**)&optscal);
+                        xfree((void**)&maxiter);
+                        xfree((void**)&cm_jump);
+                        xfree((void**)&grad_g_double);
+                        xfree((void**)&param_g);
+                        VecFree(&u_g);
+                        VecFree(&t_g);
+                        VecFree(&m_g);
+                        xfree((void**)&J);
+                        return;
+                }
 
                 /*Normalize if last gradient not satisfying (search_scalar==0)*/

issm/trunk/src/m/classes/@model/model.m

@@ -198 +198 @@
         md.fit=[];
         md.cm_jump=[];
+        md.cm_gradient=0;
         md.epsvel=0;
         md.meanvel=0;

issm/trunk/src/m/classes/@model/setdefaultparameters.m

@@ -204 +204 @@
 md.cm_max=200;
 
+%stop control solution at the gradient computation and return it? 
+md.cm_gradient=0;
+
 %eps_cm is a criteria to stop the control methods.
 %if J[n]-J[n-1]/J[n] < criteria, the control run stops

issm/trunk/src/mex/Makefile.am

@@ -37 +37 @@
                                 PenaltySystemMatrices\
                                 ProcessParams\
+                                PropagateFlagsFromConnectivity\
                                 Qmu\
                                 Reduceloadfromgtof\
@@ -178 +179 @@
                           ProcessParams/ProcessParams.h
 
+PropagateFlagsFromConnectivity_SOURCES = PropagateFlagsFromConnectivity/PropagateFlagsFromConnectivity.cpp\
+                          PropagateFlagsFromConnectivity/PropagateFlagsFromConnectivity.h
+
 Qmu_SOURCES = Qmu/Qmu.cpp\
                           Qmu/Qmu.h