Changeset 1033

Timestamp:

06/19/09 16:58:46 (16 years ago)

Author:

Mathieu Morlighem

Message:

added core solution in control

Location:

issm/trunk/src/m/solutions/cielo

Files:

• control.m (modified) (2 diffs)
• control_core.m (added)
• diagnostic_core.m (modified) (1 diff)

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

@@ -1 +1 @@
 function md=control(md)
+%CONTROL - launch control solution sequence
+%
+%   Usage:
+%      md=control(md)
+%
 
         %timing
@@ -5 +10 @@
         
         %Build all models requested for control simulation
-        models.analysis_type='control'; 
-        md.analysis_type='control';
-        md.sub_analysis_type='';
-        models.dh=CreateFemModel(md); 
+        models.analysis_type='control'; %needed for processresults
+
+        displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
+        md.analysis_type='control'; md.sub_analysis_type=''; models.dh=CreateFemModel(md); 
 
         % figure out number of dof: just for information purposes.
         md.dof=modelsize(models);
 
-        %initialize control parameters, gradients and observations
-        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(models.dh.parameters);
+        %initialize inputs
+        inputs=inputlist;
+        inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
 
-        %initialize inputs, ie models.dh.nparameters on which we invert.
-        inputs=inputlist;
-        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);
+        %compute solution
+        if ~models.dh.parameters.qmu_analysis,
+                %launch core of control solution.
+                results=control_core(models,inputs);
 
-        for n=1:models.dh.parameters.nsteps,
-                
-                %set options
-                options=optimset(options,'MaxFunEvals',models.dh.parameters.maxiter(n));
-
-                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',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
-                [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(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.');
-
-                disp('      normalizing directions...');
-                if n>=2,
-                        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
-                disp('      done.');
-                
-                %visualize direction.
-                if models.dh.parameters.plot
-                        plot_direction;
-                end
-                
-                disp('      optimizing along gradient direction...'); 
-                [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*models.dh.parameters.optscal(n)*c(n).grad_g;
-                disp('      done.');
-
-                disp('      constraining the new distribution...');    
-                param_g=ControlConstrain(param_g,models.dh.parameters);
-                disp('      done.');
-
-                %visualize direction.
-                if models.dh.parameters.plot,
-                        plot_newdistribution;
-                end
-
-                %some temporary saving 
-                if(mod(n,5)==0),
-                        solution=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
-                        save temporary_control_results solution
-                end
-                disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
-
+                %process results
+                if ~isstruct(md.results), md.results=struct(); end
+                md.results.control=processresults(models,results);
+        else
+                %launch dakota driver for diagnostic core solution
+                Qmu(models,inputs,models.dh.parameters);
         end
 
-        %Create final solution
-        disp('      preparing final velocity solution...');
-        results=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
-        disp('      done.');
-        
-        disp('      load results ontol model...');
-        if ~isstruct(md.results), md.results=struct(); end
-        md.results.control=processresults(models,results);
-        disp('      done.');
-        
+        %stop timing
         t2=clock;
-        disp(['      overall time spent on control code ' num2str(etime(t2,t1)) ' seconds'])
+        displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);

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

@@ -3 +3 @@
 %
 %   Usage:
-%      results=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
+%      results=diagnostic_core(models,inputs);
 %