Changeset 980

Timestamp:

06/13/09 11:21:49 (16 years ago)

Author:

Mathieu Morlighem

Message:

extended processresults to control

Location:

issm/trunk/src/m/solutions/cielo

Files:

• control.m (modified) (4 diffs)
• controlfinalsol.m (modified) (2 diffs)
• processresults.m (modified) (2 diffs)

issm/trunk/src/m/solutions/cielo/control.m

@@ -5 +5 @@
         
         %Build all models requested for control simulation
+        models.analysis_type='control'; 
         md.analysis_type='control';
         md.sub_analysis_type='';
-        m=CreateFemModel(md); 
+        models.dh=CreateFemModel(md); 
 
         % figure out number of dof: just for information purposes.
-        md.dof=m.nodesets.fsize; %biggest dof number
+        md.dof=modelsize(models);
 
         %initialize control parameters, gradients and observations
-        u_g_obs=m.parameters.u_g_obs;
-        u_g=m.parameters.u_g;
-        param_g=m.parameters.param_g;
-        grad_g=zeros(m.nodesets.gsize,1);
+        u_g_obs=models.dh.parameters.u_g_obs;
+        u_g=models.dh.parameters.u_g;
+        param_g=models.dh.parameters.param_g;
+        grad_g=zeros(models.dh.nodesets.gsize,1);
                 
         %set optimization options.
-        options=ControlOptions(m.parameters);
+        options=ControlOptions(models.dh.parameters);
 
-        %initialize inputs, ie m.nparameters on which we invert.
+        %initialize inputs, ie models.dh.nparameters on which we invert.
         inputs=inputlist;
-        inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
-        inputs=add(inputs,'velocity',u_g,'doublevec',3,m.parameters.numberofnodes);
+        inputs=add(inputs,models.dh.parameters.control_type,param_g,'doublevec',2,models.dh.parameters.numberofnodes);
+        inputs=add(inputs,'velocity',u_g,'doublevec',3,models.dh.parameters.numberofnodes);
 
-        for n=1:m.parameters.nsteps,
+        for n=1:models.dh.parameters.nsteps,
                 
                 %set options
-                options=optimset(options,'MaxFunEvals',m.parameters.maxiter(n));
+                options=optimset(options,'MaxFunEvals',models.dh.parameters.maxiter(n));
 
-                disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m.parameters.nsteps)]));
+                disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(models.dh.parameters.nsteps)]));
 
                 %update inputs with new fit
-                inputs=add(inputs,'fit',m.parameters.fit(n),'double');
-                inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
+                inputs=add(inputs,'fit',models.dh.parameters.fit(n),'double');
+                inputs=add(inputs,models.dh.parameters.control_type,param_g,'doublevec',2,models.dh.parameters.numberofnodes);
 
                 %Update inputs in datasets
-                [m.elements,m.nodes,m.loads,m.materials]=UpdateFromInputs(m.elements,m.nodes,m.loads,m.materials,inputs);
+                [models.dh.elements,models.dh.nodes,models.dh.loads,models.dh.materials]=UpdateFromInputs(models.dh.elements,models.dh.nodes,models.dh.loads,models.dh.materials,inputs);
 
                 disp('      computing gradJ...');
-                [u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type);
-                inputs=add(inputs,'velocity',u_g,'doublevec',2,m.parameters.numberofnodes);
+                [u_g c(n).grad_g]=GradJCompute(models.dh,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type);
+                inputs=add(inputs,'velocity',u_g,'doublevec',2,models.dh.parameters.numberofnodes);
                 disp('      done.');
 
@@ -54 +55 @@
                 
                 %visualize direction.
-                if m.parameters.plot
+                if models.dh.parameters.plot
                         plot_direction;
                 end
                 
                 disp('      optimizing along gradient direction...'); 
-                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,param_g,u_g_obs,c(n).grad_g,n,md.analysis_type,md.sub_analysis_type);
+                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,models.dh,inputs,param_g,u_g_obs,c(n).grad_g,n,md.analysis_type,md.sub_analysis_type);
                 disp('      done.');
 
                 disp('      updating parameter using optimized search scalar...');
-                param_g=param_g+search_scalar*m.parameters.optscal(n)*c(n).grad_g;
+                param_g=param_g+search_scalar*models.dh.parameters.optscal(n)*c(n).grad_g;
                 disp('      done.');
 
                 disp('      constraining the new distribution...');    
-                param_g=ControlConstrain(param_g,m.parameters);
+                param_g=ControlConstrain(param_g,models.dh.parameters);
                 disp('      done.');
 
                 %visualize direction.
-                if m.parameters.plot,
+                if models.dh.parameters.plot,
                         plot_newdistribution;
                 end
@@ -77 +78 @@
                 %some temporary saving 
                 if(mod(n,5)==0),
-                        solution=controlfinalsol(c,m,param_g,inputs,md.analysis_type,md.sub_analysis_type);
+                        solution=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
                         save temporary_control_results solution
                 end
@@ -86 +87 @@
         %Create final solution
         disp('      preparing final velocity solution...');
-        solution=controlfinalsol(c,m,param_g,inputs,md.analysis_type,md.sub_analysis_type);
+        results=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
         disp('      done.');
         
         disp('      load results ontol model...');
-        md=loadcontrolfinalsol(md,solution);
+        if ~isstruct(md.results), md.results=struct(); end
+        md.results.control=processresults(models,results);
         disp('      done.');
         

issm/trunk/src/m/solutions/cielo/controlfinalsol.m

@@ -1 @@
-function solution=controlfinalsol(c,m,p_g,inputs,analysis_type,sub_analysis_type);
+function results=controlfinalsol(c,m,param_g,inputs,analysis_type,sub_analysis_type);
+
+%initialize results
+results.time=0;
+results.step=1;
 
 %From parameters, build inputs for icediagnostic_core, using the final parameters
-inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
-u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
-
-%Build partitioning vectors to recover solution
-indx=[1:2:m.nodesets.gsize];
-indy=[2:2:m.nodesets.gsize];
-
-%Recover velocity, and parameters, in the correct partitioning.
-vx=u_g(indx);
-vy=u_g(indy);
-vel=sqrt(vx.^2+vy.^2);
-parameter=p_g(indx);
+inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
+results.u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
 
 %Recover misfit at each iteration of the control method 
@@ -20 +14 @@
         J(i)=c(i).J;
 end
+results.J=J;
 
-%Store in solution
-solution.vx=vx;
-solution.vy=vy;
-solution.vel=vel;
-solution.J=J;
-solution.parameter=parameter;
+%store param_g
+results.param_g=param_g;

issm/trunk/src/m/solutions/cielo/processresults.m

@@ -26 +26 @@
         isstokes=m_ds.parameters.isstokes;
 end
+if strcmpi(analysis_type,'control'),
+        m_dh=models.dh;
+
+        %some flags needed
+        dim=m_dh.parameters.dim;
+        ishutter=m_dh.parameters.ishutter;
+        ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+        isstokes=m_dh.parameters.isstokes;
+end
 if (strcmpi(analysis_type,'thermal') | strcmpi(analysis_type,'transient')),
         m_m=models.m;
@@ -43 +52 @@
 
                         u_g=results(i).u_g;
-                        yts=m_ds.parameters.yts;
+                        yts=m_dh.parameters.yts;
 
                         if dim==2,