Changeset 4134

Timestamp:

06/22/10 18:15:17 (15 years ago)

Author:

seroussi

Message:

start changing control_core

Location:

issm/trunk/src/m/solutions

Files:

• control_core.m (modified) (2 diffs)
• controlconvergence.m (added)
• controlrestart.m (added)

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

@@ -25 +25 @@
 
 %initialize control parameters
-param_g=model.parameters.param_g;
+param_g=femmodel.parameters.param_g;
 
 %set optimization options.
-options=ControlOptions(model.parameters);
+options=ControlOptions(femmodel.parameters);
 
-for n=1:model.parameters.NSteps,
+for n=1:nsteps,
 
         %set options
-        options=optimset(options,'MaxFunEvals',model.parameters.MaxIter(n));
+        options=optimset(options,'MaxFunEvals',femmodel.parameters.MaxIter(n));
 
-        disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(model.parameters.NSteps)]));
+        disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(femmodel.parameters.NSteps)]));
+        [femmodel.elements,femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,fit(n),FitEnum);
 
         %In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
-        if model.parameters.ControlSteady;
-                steadystate_results=steadystate_core(models); t_g=steadystate_results.t_g; 
-                inputs=add(inputs,'temperature',t_g,'doublevec',1,model.parameters.NumberOfNodes);
+        if control_steady;
+                femmodel=steadystate_core(femmodel);
         end
 
-        %update inputs with new fit
-        inputs=add(inputs,'fit',model.parameters.Fit(n),'double');
-        inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes);
-
         %Update inputs in datasets
-        [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters);
+        [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=UpdateFromInputs(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 
         displaystring(verbose,'\n%s',['      computing gradJ...']);
-        results_grad=gradjcompute_core(models);
-        u_g=results_grad.u_g; c(n).grad_g=results_grad.grad_g;
-        if dim==3,
-                if isstokes,
-                        inputs=add(inputs,'velocity',u_g,'doublevec',4,model.parameters.NumberOfNodes);
-                else
-                        if model.parameters.ControlSteady;
-                                inputs=add(inputs,'velocity',u_g,'doublevec',3,model.parameters.NumberOfNodes);
-                        else
-                                inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes);
-                        end
-                end
-        else
-                inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes);
-        end
+        femmodel=gradient_core(femmodel);
 
-        if n>=2 & search_scalar==0,
-                displaystring(verbose,'\n%s',['      normalizing directions...']);
-                c(n).grad_g=Orth(c(n).grad_g,c(n-1).grad_g);
-        else
-                c(n).grad_g=Orth(c(n).grad_g,{});
-        end
-
-        %visualize direction.
-        if model.parameters.Plot
-                plot_direction;
+        %Return gradient if asked
+        if cm_gradient,
+                InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel.parameters,GradientEnum);
+                %%%%%%%%%%%%%%%%%%%%%% PROBLEM MUST GO TO THE END (SEE PARALLEL)
+                error
+                break;
         end
 
         displaystring(verbose,'\n%s',['      optimizing along gradient direction...']);
-        [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,models,param_g,c(n).grad_g,n,model.parameters);
+        [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,param_g,c(n).grad_g,n,femmodel.parameters);
 
         displaystring(verbose,'\n%s',['      updating parameter using optimized search scalar...']);
-        param_g=param_g+search_scalar*model.parameters.Optscal(n)*c(n).grad_g;
+        param_g=param_g+search_scalar*femmodel.parameters.Optscal(n)*c(n).grad_g;
 
         displaystring(verbose,'\n%s',['      constraining the new distribution...']);
-        param_g=ControlConstrain(param_g,model.parameters);
+        param_g=ControlConstrain(param_g,femmodel.parameters);
         
         disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
 
         %Has convergence been reached?
-        convergence=0;
-        if ~isnan(eps_cm),
-                i=n-2;
-                %go through the previous misfits(starting from n-2)
-                while (i>=1),
-                        if (fit(i)==fit(n)),
-                                %convergence test only if we have the same misfits
-                                if ((c(i).J-c(n).J)/c(n).J <= eps_cm),
-                                        %convergence if convergence criteria fullfilled
-                                        convergence=1;
-                                        displaystring(verbose,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'<',eps_cm);
-                                else
-                                        displaystring(verbose,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'>',eps_cm);
-                                end
-                                break;
-                        end
-                        i=i-1;                                                                                                                                         
-                end
+        converged=controlconvergence(J,fit,eps_cm,n);
+        if converged,
+                break;
         end
-        %stop if convergence has been reached                                                                                                               
-        if (convergence), break; end
+
+        %Temporary saving every five iterations
+        if mod(n+1,5)==0;
+                displaystring(verbose,'\n%s',['      saving temporary results...']);
+                controlrestart(femmmodel,J);
+        end
 
 end
@@ -117 +84 @@
 %generate output
 displaystring(verbose,'\n%s',['      preparing final velocity solution...']);
-
-%compute final velocity from diagnostic_core (horiz+vertical)
-if model.parameters.ControlSteady;
-        inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes);
-        steadystate_results=steadystate_core(models); t_g=steadystate_results.t_g; 
-        u_g=steadystate_results.u_g;
-        t_g=steadystate_results.t_g;
-        m_g=steadystate_results.m_g;
+if control_steady,
+        femmodel=steadystate_core(femmodel);
 else
-        inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes);
-        results_diag=diagnostic_core(models);
-        u_g=results_diag.u_g;
+        femmodel=diagnostic_core(femmodel);
 end
 
-%Recover misfit at each iteration of the control method 
-J=zeros(length(c),1);
-for i=1:length(c),
-        J(i)=c(i).J;
-end
+%Some results not computed by diagnostic or steadystate
+InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type);
+femmodel.results.JEnum=J;
+femmodel.results.ControlTypeEnum=EnumAsString(control_type);
 
-%build results
-results.time=0;
-results.step=1;
-results.J=J;
-results.param_g=param_g;
-results.u_g=u_g;
-if model.parameters.ControlSteady,
-        results.t_g=t_g;
-        results.m_g=m_g;
-end
+end %end function