Changeset 1652

Timestamp:

08/11/09 16:51:13 (16 years ago)

Author:

Mathieu Morlighem

Message:

simplified CM

Location:

issm/trunk/src/m/solutions/cielo

Files:

• ControlInitialization.m (added)
• ControlPrepareStokes.m (deleted)
• control_core.m (modified) (3 diffs)

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

@@ -6 +6 @@
 %
 
-%recover models
-m_dh=models.dh;
-m_dv=models.dv;
-m_ds=models.ds;
-m_sl=models.sl;
+%Preprocess models
+[inputs model]=ControlInitialization(models,inputs);
 
 %recover parameters common to all solutions
-debug=m_dh.parameters.debug;
-dim=m_dh.parameters.dim;
-ishutter=m_dh.parameters.ishutter;
-ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
-isstokes=m_dh.parameters.isstokes;
+debug=model.parameters.debug;
 
 %initialize control parameters
-param_g=models.dh.parameters.param_g;
+param_g=model.parameters.param_g;
 
 %set optimization options.
-options=ControlOptions(m_dh.parameters);
+options=ControlOptions(model.parameters);
 
-%Take car of Stokes : compute slope and get spc once for all
-if isstokes,
-        [inputs m_ds]=ControlPrepareStokes(models,inputs);
-end
-
-for n=1:m_dh.parameters.nsteps,
+for n=1:model.parameters.nsteps,
 
         %set options
-        options=optimset(options,'MaxFunEvals',m_dh.parameters.maxiter(n));
+        options=optimset(options,'MaxFunEvals',model.parameters.maxiter(n));
 
-        disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m_dh.parameters.nsteps)]));
+        disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(model.parameters.nsteps)]));
 
         %update inputs with new fit
-        inputs=add(inputs,'fit',m_dh.parameters.fit(n),'double');
-        inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',1,m_dh.parameters.numberofnodes);
+        inputs=add(inputs,'fit',model.parameters.fit(n),'double');
+        inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes);
 
         %Update inputs in datasets
-        if isstokes,
-                [m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials]=UpdateFromInputs(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,inputs);
-        else
-                [m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials]=UpdateFromInputs(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,inputs);
-        end
+        [model.elements,model.nodes,model.loads,model.materials]=UpdateFromInputs(model.elements,model.nodes,model.loads,model.materials,inputs);
 
         displaystring(debug,'\n%s',['      computing gradJ...']);
-        if isstokes,
-                [u_g c(n).grad_g]=GradJCompute(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
-                inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes);
-        else
-                [u_g c(n).grad_g]=GradJCompute(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-                inputs=add(inputs,'velocity',u_g,'doublevec',2,m_dh.parameters.numberofnodes);
-        end
+        [u_g c(n).grad_g]=GradJCompute(model,inputs,model.parameters.analysis_type,model.parameters.sub_analysis_type);
+        inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.numberofnodes);
 
         displaystring(debug,'\n%s',['      normalizing directions...']);
@@ -65 +44 @@
 
         %visualize direction.
-        if m_dh.parameters.plot
+        if model.parameters.plot
                 plot_direction;
         end
 
         displaystring(debug,'\n%s',['      optimizing along gradient direction...']);
-        if isstokes,
-                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_ds,inputs,param_g,c(n).grad_g,n,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
-        else
-                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_dh,inputs,param_g,c(n).grad_g,n,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-        end
+        [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,model,inputs,param_g,c(n).grad_g,n,model.parameters.analysis_type,model.parameters.sub_analysis_type);
 
         displaystring(debug,'\n%s',['      updating parameter using optimized search scalar...']);
-        param_g=param_g+search_scalar*m_dh.parameters.optscal(n)*c(n).grad_g;
+        param_g=param_g+search_scalar*model.parameters.optscal(n)*c(n).grad_g;
 
         displaystring(debug,'\n%s',['      constraining the new distribution...']);
-        param_g=ControlConstrain(param_g,m_dh.parameters);
+        param_g=ControlConstrain(param_g,model.parameters);
 
         %visualize direction.
-        if m_dh.parameters.plot,
+        if model.parameters.plot,
                 plot_newdistribution;
         end
@@ -94 +69 @@
 
 %compute final velocity from diagnostic_core (horiz+vertical)
-inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',1,m_dh.parameters.numberofnodes);
+inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes);
 results_diag=diagnostic_core(models,inputs);