0001 function [u_g varargout]=icediagnostic_core_nonlinear(m,analysis_type, varargin);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 global cluster gridset element_debug element_debugid
0015
0016
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
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
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;
0044 soln(2).u_g=velocity_param;
0045 soln(2).u_g(4:6:gridset.gsize)=0;
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
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
0069
0070 kflag=1;
0071
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
0079 inputs.velocity=soln(count).u_g;
0080 inputs.oldvelocity=soln(count-1).u_g;
0081
0082
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
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
0099
0100
0101 [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn);
0102
0103
0104 p_f=IceReducerightside(p_g,G_mn,K_fs,ys);
0105
0106
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
0113 count=count+1;
0114
0115
0116 soln(count).u_g= full(Mergesolution_g(u_f,G_mn, ys));
0117
0118
0119 if((count>=4)| velocity_is_present),
0120
0121
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
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
0147 u_g=soln(end).u_g;
0148
0149
0150 nout=max(nargout,1)-1;
0151 if nout==2,
0152
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