Changeset 11306

Timestamp:

02/02/12 11:31:06 (13 years ago)

Author:

Mathieu Morlighem

Message:

improving CM algorithm

Location:

issm/trunk-jpl/src

Files:

• c/solutions/control_core.cpp (modified) (4 diffs)
• c/solutions/gradient_core.cpp (modified) (3 diffs)
• c/solutions/solutions.h (modified) (1 diff)
• m/solutions/control_core.m (modified) (3 diffs)
• m/solutions/gradient_core.m (modified) (3 diffs)

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

@@ -34 +34 @@
                 
         /*intermediary: */
-        double  search_scalar=0;
+        double  search_scalar=1;
         OptArgs optargs;
         OptPars optpars;
@@ -84 +84 @@
                 femmodel->parameters->SetParam(step_responses,1,num_responses,StepResponsesEnum);
                 
-                /*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */
+                /*In steady state inversion, compute new temperature field now*/
                 if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
 
                 _printf_(VerboseControl(),"%s\n","   compute adjoint state:");
                 adjointcore(femmodel);
-                gradient_core(femmodel,n,search_scalar);
+                gradient_core(femmodel,n,search_scalar==0);
 
                 /*Return gradient if asked: */
@@ -100 +100 @@
                 optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
                 BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
-                //OptimalSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
 
                 _printf_(VerboseControl(),"%s\n","   updating parameter using optimized search scalar"); //true means update save controls
@@ -116 +115 @@
                 for(i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);
                 femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
-                //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,InversionControlParametersEnum,EnumToStringx(control_type),1,0));
         }
 

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

@@ -13 +13 @@
 #include "../solvers/solvers.h"
 
-
-void gradient_core(FemModel* femmodel,int step, double search_scalar){ //step defaults to 0, search_scalar defaults to 0
+void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ 
         
         /*parameters: */
@@ -40 +39 @@
                 Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type[i]);
 
-                if (step>0 && search_scalar==0){
+                if (orthogonalize){
                         _printf_(VerboseControl(),"   orthogonalization\n");
                         ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]);
@@ -53 +52 @@
                 if(norm_grad<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm of J with respect to %s is zero",EnumToStringx(control_type[i]));
                 if(isnan(norm_grad))_error_("||∂J/∂α||∞ = NaN  gradient norm of J with respect to %s is NaN" ,EnumToStringx(control_type[i]));
-                if(i==0 || (optscal_list[num_controls*step+i]/norm_grad)<optscal) optscal=optscal_list[num_controls*step+i]/norm_grad;
+                optscal=optscal_list[num_controls*step+i]/norm_grad;
 
                 /*plug back into inputs: */

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

@@ -15 +15 @@
 void adjointdiagnostic_core(FemModel* femmodel);
 void adjointbalancethickness_core(FemModel* femmodel);
-void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0);
+void gradient_core(FemModel* femmodel,int n=0,bool orthogonalize=false);
 void diagnostic_core(FemModel* femmodel);
 void hydrology_core(FemModel* femmodel);

issm/trunk-jpl/src/m/solutions/control_core.m

@@ -25 +25 @@
         %Initialize misfits with a vector of zeros
         J=zeros(nsteps,1);
-        search_scalar=0;
+        search_scalar=1;
 
         %Get core from solution type
@@ -42 +42 @@
                 femmodel.parameters.StepResponses=responses(n,:);
 
-                %In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
+                %In steady state inversion, compute new temperature field now
                 if (solution_type==SteadystateSolutionEnum)
                         femmodel=steadystate_core(femmodel);
@@ -50 +50 @@
                 eval(['femmodel=' adjointcore '(femmodel);']);
 
-                femmodel=gradient_core(femmodel,n,search_scalar);
+                femmodel=gradient_core(femmodel,n,search_scalar==0);
 
                 %Return gradient if asked

issm/trunk-jpl/src/m/solutions/gradient_core.m

@@ -11 +11 @@
 if nargin==3,
         step=varargin{1};
-        search_scalar=varargin{2};
+        orthogonalize=varargin{2};
 elseif nargin==1
         step=0;
-        search_scalar=0;;
+        orthogonalize=false;
 else
         help gradient_core
@@ -31 +31 @@
                 grad=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type(i));
 
-                if (step>1 && search_scalar==0),
+                if orthogonalize,
                         issmprintf(VerboseControl,'%s',['   orthogonalization']);
                         old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i));
@@ -43 +43 @@
                  if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm of J with respect to ' EnumToString(control_type(i))  ' is zero']); end
                  if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i))  ' is NaN' ]); end
-                 if(i==1 | (gradient_scaling_list(step,i)/norm_grad)<gradient_scaling) gradient_scaling=gradient_scaling_list(step,i)/norm_grad; end
+                 gradient_scaling=gradient_scaling_list(step,i)/norm_grad;
 
                 %plug back into inputs: