Changeset 4537

Timestamp:

07/12/10 16:36:46 (15 years ago)

Author:

seroussi

Message:

keep working on control serial

Location:

issm/trunk/src/m/solutions

Files:

• control_core.m (modified) (2 diffs)
• controlconvergence.m (modified) (1 diff)
• objectivefunctionC.m (modified) (1 diff)

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

@@ -28 +28 @@
         options.MaxIter=femmodel.parameters.MaxIter;
 
+        %Initialize misfits with a vector of zeros
+        J=zeros(nsteps,1);
+
         for n=1:nsteps,
 
@@ -48 +51 @@
 
                 displaystring(verbose,'\n%s',['      optimizing along gradient direction...']);
-                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,param_g,c(n).grad_g,n,femmodel.parameters);
+                [search_scalar J(n)]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,n,femmodel.parameters);
 
-                displaystring(verbose,'\n%s',['      updating parameter using optimized search scalar...']);
-                param_g=param_g+search_scalar*femmodel.parameters.Optscal(n)*grad_g;
+                displaystring('\n%s',['      updating parameter using optimized search scalar:']);
+                scalar=search_scalar*optscal(n);
+                femmodel.elements=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,ControlParameterEnum);
 
-                displaystring(verbose,'\n%s',['      constraining the new distribution...']);
-                param_g=ControlConstrain(param_g,femmodel.parameters);
+                displaystring('\n%s',['      constraning the new distribution...']);
+                [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
 
-                disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
+                displaystring('\n%s',['      save new parameter...']);
+                femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,ControlParameterEnum);
+
+                disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]);
 
                 %Has convergence been reached?

issm/trunk/src/m/solutions/controlconvergence.m

@@ -1 @@
-function converged=controlconvergence(J,fit,eps_cm,n)
+function convergence=controlconvergence(J,fit,eps_cm,n)
 %CONTROLCONVERGENCE - determine the convergence of control_core solution
 %

issm/trunk/src/m/solutions/objectivefunctionC.m

@@ -1 @@
-function J =objectivefunctionC(search_scalar,models,p_g,grad_g,n,analysis_type,sub_analysis_type);
+function J =objectivefunctionC(search_scalar,femmodel,n);
+%OBJECTIVEFUNCTIONC - objective function that return a parameter for a certain function
 
-%recover active model.
-m=models.active;
+conserve_loads=true;
+%recover some parameters
+optscal=femmodel.parameters.OptScal(n);
+fit=femmodel.parameters.Fit(n);
+control_type=femmodel.parameters.ControlType;
+control_steady=femmodel.parameters.ControlSteady;
+analysis_type=femmodel.parameters.AnalysisType;
+isstokes=femmodel.parameters.IsStokes;
+cm_min=femmodel.parameters.CmMin;
+cm_max=femmodel.parameters.CmMax;
 
-%recover some parameters
-optscal=m.parameters.Optscal(n);
-fit=m.parameters.Fit(n);
-control_type=m.parameters.ControlType;
-analysis_type=m.parameters.AnalysisType;
+%set current configuration
+if isstokes,
+        femmodel=SetCurrentConfiguration(femmodel,DiagnsoticStokesAnalysisEnum);
+else
+        femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
+end
 
-%Update along gradient using scalar supplied by fmincon optimization routine
-parameter=p_g+search_scalar*optscal*grad_g;
+%Use ControlParameterEnum input to  reinitialize our input parameter:
+femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);
 
-%Plug parameter into inputs
-inputs=add(inputs,m.parameters.ControlType,parameter,'doublevec',1,m.parameters.NumberOfNodes);
+%Use search scalar to shoot parameter in the gradient direction:
+scalar=search_scalar*optscal;
+femmodel.elements=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,ControlParameterEnum);
 
-%Run diagnostic with updated parameters.
-u_g=diagnostic_core_nonlinear(m,analysis_type,sub_analysis_type);
+%Constrain:
+[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
 
-%add velocity to inputs.
-inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes);
+%Run diagnostic with updated inputs:
+if(control_steady==0),
+        femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);  %true means we conserve loads at each diagnostic run
+else
+        femmodel=diagnostic_core(femmodel);  %We need a 3D velocity!! (vz is required for the next thermal run)
+end
 
-%Compute misfit for this velocity field. 
-J=CostFunction(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
+%Compute misfit for this velocity field
+[femmodel.elements,femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices, femmodel.loads,femmodel.materials,femmodel.parameters,fit,FitEnum);
+J=CostFunction(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials, femmodel.parameters);