Index: /issm/trunk/src/m/solutions/adjoint_core.m =================================================================== --- /issm/trunk/src/m/solutions/adjoint_core.m (revision 4535) +++ /issm/trunk/src/m/solutions/adjoint_core.m (revision 4535) @@ -0,0 +1,50 @@ +function femmodel=adjoint_core(femmodel), +%ADJOINT_CORE - compute inverse method adjoint state + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + isstokes=femmodel.parameters.IsStokes; + dim=femmodel.parameters.Dim; + solution_type=femmodel.parameters.SolutionType; + conserve_loads=true; + + %set analysis type to compute velocity: + if(isstokes), + femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum); + else + femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum); + end + + displaystring('\n%s',[' recover solution for this stiffness and right hand side:']); + [femmodel,ug,K_ff0,K_fs0]=solver_diagnostic_nonlinear(femmodel,conserve_loads); + + displaystring('\n%s',[' Build Du, difference between observed and modeled velocities:']); + du_g=Du(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters); + + displaystring('\n%s',[' reduce adjoint load from g_set to f_set:']); + du_f=Reduceloadfromgtof(du_g,femmodel.Gmn,K_fs0,femmodel.ys,femmodel.nodesets,true); %true means that ys0 flag is activated: all spcs show 0 displacement + + displaystring('\n%s',[' solve for adjoint vector:']); + adjoint_f=Solver(K_ff0,du_f,[],femmodel.parameters); + + displaystring('\n%s',[' merge back to g set:']); + adjoint_g=Mergesolutionfromftog(adjoint_f,femmodel.Gmn,femmodel.ys,femmodel.nodesets,true);%true means that ys0 flag is activated: all spc are 0 + + %Update inputs using adjoint solution, and same type of setup as diagnostic solution: + if(isstokes), + femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum); + else + femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum); + end + + femmodel.elements=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,adjoint_g); + + displaystring(verbose,'\n%s',[' saving results...']); + if(solution_type==AdjointSolutionEnum) + femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointxEnum); + femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointyEnum); + if(dim==3), + femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointzEnum); + end + end + Index: /issm/trunk/src/m/solutions/control_core.m =================================================================== --- /issm/trunk/src/m/solutions/control_core.m (revision 4534) +++ /issm/trunk/src/m/solutions/control_core.m (revision 4535) @@ -7,5 +7,5 @@ %Preprocess models - femmodel=ControlInitialization(femmodel); + femmodel=stokescontrolinit(femmodel); %recover parameters common to all solutions @@ -24,68 +24,64 @@ 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 Index: /issm/trunk/src/m/solutions/gradient_core.m =================================================================== --- /issm/trunk/src/m/solutions/gradient_core.m (revision 4535) +++ /issm/trunk/src/m/solutions/gradient_core.m (revision 4535) @@ -0,0 +1,55 @@ +function femmodel=gradient_core(femmodel,varargin), +%GRADIENT_CORE - Brief compute inverse method gradient direction +% +% Usage: +% femmodel=gradient_core(femmodel,varargin); +% +% Examples: +% femmodel=gradient_core(femmodel); +% femmodel=gradient_core(femmodel,step,search_scalar); + +if nargin==3, + step=varargin{1}; + search_scalar=varargin{2}; +elseif nargin==1 + step=0; + search_scalar=0;; +else + help gradient_core + error('bad usage'); +end + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + control_type=femmodel.parameters.ControlType; + control_steady=femmodel.parameters.ControlSteady; + + displaystring(verbose,'\n%s',[' compute adjoint state...']); + femmodel=adjoint_core(femmodel); + + displaystring(verbose,'\n%s',[' compute gradient...']); + grad=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters); + + displaystring(verbose,'\n%s',[' retrieve old gradient...']); + old_gradient=GetVectorFromInputs(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,OldGradientEnum,VertexEnum); + + if control_steady; + femmodel=diagnostic_core(femmodel); + end + + if (step>0 && search_scalar==0), + displaystring(verbose,'\n%s',[' orthogonalization...']); + new_gradient=Orth(grad,old_gradient); + else + displaystring(verbose,'\n%s',[' normalizing direction...']); + new_gradient=Orth(grad,[]); + end + displaystring(verbose,'\n%s',[' done...']); + + %point gradient and old_gradient to new_gradient: + grad=new_gradient; + old_gradient=new_gradient; + + %plug back into inputs: + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputUpdateFromVector(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,grad,GradientEnum,VertexEnum); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputUpdateFromVector(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,old_gradient,OldGradientEnum,VertexEnum); Index: /issm/trunk/src/m/solutions/stokescontrolinit.m =================================================================== --- /issm/trunk/src/m/solutions/stokescontrolinit.m (revision 4535) +++ /issm/trunk/src/m/solutions/stokescontrolinit.m (revision 4535) @@ -0,0 +1,35 @@ +function femmodel=stokescontrolinit(femmodel), +%STOKESCONTROLINIT - initialize femmodel if control method on Stokes + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + isstokes=femmodel.parameters.IsStokes; + stokesreconditioning=femmodel.parameters.StokesReconditioning; + + %if no Stokes analysis carried out, just return + if (isstokes==0), + return; + end + + % For Stokes inverse control method, we are going to carry out the inversion only on the Stokes part. So we need to solve the Hutter or MacAyeal/Pattyn femmodel here, and constrain the Stokes femmodel using the Hutter or MacAyeal/Pattyn at the boundary. We don't want to have to do that at every inversion step, as it needs be done only once: + + %Compute slopes: + femmodel=bedslope_core(femmodel); + + %Run a complete diagnostic to update the Stokes spcs: + femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum); + femmodel=solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads); + + %vertical velocity + femmodel=solver_linear(NULL,femmodel); + + %recondition pressure computed previously: + femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureEnum,PressureStokesEnum); + femmode.elements=InputScale(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureStokesEnum,1.0/stokesreconditioning); + + displaystring(verbose,'\n%s',[' update boundary conditions for stokes using velocities previously computed...']); + femmodel=ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum); + + displaystring(verbose,'\n%s',[' computing stokes velocity and pressure...']); + femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum); + femmode=solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);