Changeset 4535

Timestamp:

07/12/10 12:37:02 (15 years ago)

Author:

seroussi

Message:

started to work on control methods serial

Location:

issm/trunk/src/m/solutions

Files:

• adjoint_core.m (added)
• control_core.m (modified) (2 diffs)
• gradient_core.m (added)
• stokescontrolinit.m (added)

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

@@ -7 +7 @@
 
         %Preprocess models
-        femmodel=ControlInitialization(femmodel);
+        femmodel=stokescontrolinit(femmodel);
 
         %recover parameters common to all solutions
@@ -24 +24 @@
         control_steady=femmodel.parameters.ControlSteady;
 
-%initialize control parameters
-param_g=femmodel.parameters.param_g;
+        %Initialise options with tolerance and maxiter
+        options.TolX=femmodel.parameters.TolX;
+        options.MaxIter=femmodel.parameters.MaxIter;
 
-%set optimization options.
-options=ControlOptions(femmodel.parameters);
+        for n=1:nsteps,
 
-for n=1: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);
 
-        %set options
-        options=optimset(options,'MaxFunEvals',femmodel.parameters.MaxIter(n));
+                %In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
+                if control_steady;
+                        femmodel=steadystate_core(femmodel);
+                end
 
-        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);
+                displaystring(verbose,'\n%s',['      computing gradJ...']);
+                femmodel=gradient_core(femmodel);
 
-        %In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
-        if control_steady;
-                femmodel=steadystate_core(femmodel);
+                %Return gradient if asked
+                if cm_gradient,
+                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,GradientEnum);
+                        return;
+                end
+
+                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);
+
+                displaystring(verbose,'\n%s',['      updating parameter using optimized search scalar...']);
+                param_g=param_g+search_scalar*femmodel.parameters.Optscal(n)*grad_g;
+
+                displaystring(verbose,'\n%s',['      constraining the new distribution...']);
+                param_g=ControlConstrain(param_g,femmodel.parameters);
+
+                disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
+
+                %Has convergence been reached?
+                converged=controlconvergence(J,fit,eps_cm,n);
+                if converged,
+                        break;
+                end
+
         end
 
-        %Update inputs in datasets
-        [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...']);
-        femmodel=gradient_core(femmodel);
-
-        %Return gradient if asked
-        if cm_gradient,
-                femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,GradientEnum);
-                %%%%%%%%%%%%%%%%%%%%%% PROBLEM MUST GO TO THE END (SEE PARALLEL)
-                error
-                break;
+        %generate output
+        displaystring(verbose,'\n%s',['      preparing final velocity solution...']);
+        if control_steady,
+                femmodel=steadystate_core(femmodel);
+        else
+                femmodel=diagnostic_core(femmodel);
         end
 
-        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);
-
-        displaystring(verbose,'\n%s',['      updating parameter using optimized search scalar...']);
-        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,femmodel.parameters);
-        
-        disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
-
-        %Has convergence been reached?
-        converged=controlconvergence(J,fit,eps_cm,n);
-        if converged,
-                break;
-        end
-
-end
-
-%generate output
-displaystring(verbose,'\n%s',['      preparing final velocity solution...']);
-if control_steady,
-        femmodel=steadystate_core(femmodel);
-else
-        femmodel=diagnostic_core(femmodel);
-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);
+        %Some results not computed by diagnostic or steadystate
+        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum);
+        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum);
+        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VelEnum);
+        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,GradientEnum);
+        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parametersn,AdjointxEnum);
+        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parametersn,AdjointyEnum);
+        if (dim==3) femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum);
+        femmodel.results.JEnum=J;
+        femmodel.results.ControlTypeEnum=EnumAsString(control_type);
 
 end %end function