icediagnostic_core_nonlinear

PURPOSE ^

ICEDIAGNOSTIC_CORE_NONLINEAR - core solution of non linear problems

SYNOPSIS ^

function [u_g varargout]=icediagnostic_core_nonlinear(m,analysis_type, varargin);

DESCRIPTION ^

ICEDIAGNOSTIC_CORE_NONLINEAR - core solution of non linear problems

   Core of the diagnostic solution, deals with the Direct Shooting Method
   to linearize the non-linear material equations.
   varargin is for extra inputs to solution sequence.
   
   Usage:
      [u_g varargout]=icediagnostic_core_nonlinear(m,analysis_type,varargin);

   See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_LINEAR

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [u_g varargout]=icediagnostic_core_nonlinear(m,analysis_type, varargin);
0002 %ICEDIAGNOSTIC_CORE_NONLINEAR - core solution of non linear problems
0003 %
0004 %   Core of the diagnostic solution, deals with the Direct Shooting Method
0005 %   to linearize the non-linear material equations.
0006 %   varargin is for extra inputs to solution sequence.
0007 %
0008 %   Usage:
0009 %      [u_g varargout]=icediagnostic_core_nonlinear(m,analysis_type,varargin);
0010 %
0011 %   See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_LINEAR
0012 
0013 %global variables
0014 global cluster gridset element_debug element_debugid
0015 
0016 %recover fem model fields
0017 elements=m.elements;
0018 grids=m.grids;
0019 materials=m.materials;
0020 loads=m.loads;
0021 ys=m.ys;
0022 gridset=m.gridset;
0023 params=m.params;
0024 G_mn=m.Gmn;
0025 
0026 %recover parameters
0027 sparsity=params.sparsity;
0028 solver_type=params.solver_type;
0029 eps_abs=params.eps_abs;
0030 eps_rel=params.eps_rel;
0031 yts=params.yts;
0032 debug=params.debug;
0033 element_debug=params.element_debug;
0034 element_debugid=params.element_debugid;
0035 
0036 %recover existing velocity if given in input and initialize solution
0037 if nargin==3,
0038     inputs=varargin{1};
0039 
0040     [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0041     if velocity_is_present
0042         soln(1).u_g=velocity_param;
0043         soln(1).u_g(4:6:gridset.gsize)=0; %no pressure in the velocity
0044         soln(2).u_g=velocity_param;
0045         soln(2).u_g(4:6:gridset.gsize)=0; %no pressure in the velocity
0046     else    
0047         soln(1).u_g={};
0048         soln(2).u_g={};
0049     end
0050 else
0051     velocity_is_present=0;
0052     inputs=struct();
0053     soln(1).u_g={};
0054     soln(2).u_g={};
0055 end
0056 
0057 
0058 %Initialization
0059 count=2;
0060 converged=0;
0061     
0062 if debug,
0063     disp(sprintf('%s','   starting direct shooting method'));
0064 end
0065 
0066 while(~converged),
0067 
0068     % Generate system matrices (stiffness and load)
0069     %compute stiffness matrix flag
0070     kflag=1;
0071     %compute loads ?
0072     if strcmp(analysis_type, 'diagnostic_stokes'),
0073         pflag=1;
0074     else
0075         if count==2, pflag=1; else pflag=0; end
0076     end
0077 
0078     %add velocity to inputs
0079     inputs.velocity=soln(count).u_g;
0080     inputs.oldvelocity=soln(count-1).u_g;
0081     
0082     %generate stiffness and loads
0083     [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type);
0084     [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g);
0085 
0086     %Save loads
0087     if count==2, 
0088         p_g_old=p_g;
0089     else
0090         p_g=p_g_old;
0091     end
0092     
0093     if cluster, 
0094         K_gg=distributed(gplus(K_gg),'convert');
0095         p_g=gplus(p_g);
0096     end
0097     
0098     % Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints),
0099     % and compute modifications to loads from single point constraints.
0100     
0101     [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); 
0102 
0103     % Reduce load from g set to f set
0104     p_f=IceReducerightside(p_g,G_mn,K_fs,ys);
0105 
0106     % Solve
0107     u_f=Solver(K_ff,p_f,solver_type);
0108     if debug,
0109         disp(sprintf('%s%g','      condition number of stiffness matrix: ',condest(K_ff)));
0110     end
0111    
0112     %increment index
0113     count=count+1;
0114 
0115     % Add single point constraints back, ie increase f-set by s-set into the global g-set.
0116     soln(count).u_g= full(Mergesolution_g(u_f,G_mn, ys)); %make solution full
0117 
0118     %Figure out if convergence is reached.
0119     if((count>=4)| velocity_is_present), %make sure we already iterated at least once.
0120 
0121         %compute relative velocity difference for this step.
0122         relative_change=norm(soln(count).u_g-soln(count-1).u_g,2)/norm(soln(count-1).u_g,2);
0123 
0124         if relative_change<eps_rel, 
0125             if debug, disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',relative_change,' < ',eps_rel)); end
0126             converged=1;
0127         else
0128             if debug, disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',relative_change,' > ',eps_rel)); end
0129             converged=0;
0130         end
0131 
0132         if ~isnan(eps_abs)
0133             %compute velocity difference for this step.
0134             change=max(abs(soln(count).u_g-soln(count-1).u_g))*yts;
0135 
0136             if change<eps_abs, 
0137                 if debug, disp(sprintf('%s %g %s %g %s','      convergence criterion: max(du)=',change,' < ',eps_abs,'m/yr')); end
0138             else
0139                 if debug, disp(sprintf('%s %g %s %g %s','      convergence criterion: max(du)=',change,' > ',eps_abs,'m/yr')); end
0140                 converged=0;
0141             end
0142         end
0143     end
0144 end
0145 
0146 %prepare output
0147 u_g=soln(end).u_g;
0148 
0149 %more output might be needed, for ex when running in cielocontrol.m
0150 nout=max(nargout,1)-1;
0151 if nout==2,
0152     %K_ff and K_fs are requested, for the case where we have no loads (ie basal drag)
0153     inputs.drag=zeros(gridset.gsize,1); kflag=1;pflag=1; 
0154     inputs.velocity=soln(count).u_g;
0155 
0156     [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type);
0157     [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g);
0158     [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); 
0159     varargout(1)={K_ff};
0160     varargout(2)={K_fs};
0161 end

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003