Index: sm/trunk/src/m/solutions/ice/IceNormalizeconstraints.m =================================================================== --- /issm/trunk/src/m/solutions/ice/IceNormalizeconstraints.m (revision 33) +++ (revision ) @@ -1,20 +1,0 @@ -function G_mn=NormalizeConstraints(R_mg) -%NORMALIZECONSTRAINTS - create G_mn -% -% see NastranRgbDoc documentation in theoryguide directory. -% -% Usage: -% G_mn=NormalizeConstraints(R_mg) - -global gridset - -if(gridset.msize>0) - - R_mm = R_mg(:, gridset.pv_m); - R_mn = R_mg(:, gridset.pv_n); - - G_mn= - inv(R_mm) * R_mn; - -else - G_mn= []; -end Index: sm/trunk/src/m/solutions/ice/IcePenaltyConstraints.m =================================================================== --- /issm/trunk/src/m/solutions/ice/IcePenaltyConstraints.m (revision 33) +++ (revision ) @@ -1,25 +1,0 @@ -function [loads,converged]=IcePenaltyConstraints(loads,grids,materials,params,inputs) -%ICEPENALTYCONSTRAINTS - test if the thermal solution has converged -% -% This routine check if each element have converged in term of penalty activation -% once all the penalties remain the same between 2 iterations, return converged=1 -% -% Usage: -% [loads,converged]=IcePenaltyConstraints(loads,grids,materials,params,inputs) - -min_thermal_constraints=params.min_thermal_constraints; -unstable_count=0; - -%Go through constraints and set active flags -for n=1:length(loads), - [loads(n).load unstable]=PenaltyConstrain(loads(n).load,grids,materials,inputs); - unstable_count=unstable_count+unstable; -end - -disp(sprintf('%s %g',' number of unstable melting constraints: ',unstable_count)); - -if unstable_count<=min_thermal_constraints, - converged=1; -else - converged=0; -end Index: sm/trunk/src/m/solutions/ice/IceReducematrix.m =================================================================== --- /issm/trunk/src/m/solutions/ice/IceReducematrix.m (revision 33) +++ (revision ) @@ -1,22 +1,0 @@ -function Knn= IceReducematrix( Kgg, Gmn, pv_m,pv_n,msize,flag) -%ICEREDUCEMATRIX - reduce stiffness matrix from g to n -% -% Usage: -% Knn= IceReducematrix( Kgg, Gmn, pv_m,pv_n,msize,flag) - -if msize>0, - Kmm=Kgg(pv_m,pv_m); - Kmn=Kgg(pv_m,pv_n); - Knm=Kgg(pv_n,pv_m); - Knn=Kgg(pv_n,pv_n); - - Knn= Knn + Knm * Gmn; - - if flag~=1, - Kmn = Kmn + Kmm * Gmn; - Knn = Knn + Gmn' * Kmn; - end - -else - Knn=Kgg; -end Index: sm/trunk/src/m/solutions/ice/IceReducematrixfromgtof.m =================================================================== --- /issm/trunk/src/m/solutions/ice/IceReducematrixfromgtof.m (revision 33) +++ (revision ) @@ -1,20 +1,0 @@ -function [Kff Kfs]=IceReducematrixfromgtof(Kgg,Gmn) -%ICEREDUCEMATRIXFROMGTOF - reduce stiffness matrix from gset to fset -% -% Usage: -% [Kff Kfs]=IceReducematrixfromgtof(Kgg,Gmn) - -global gridset - -%Reduce matrix from g-size to n-size -Knn= IceReducematrix( Kgg, Gmn, gridset.pv_m, gridset.pv_n, gridset.msize, 2 ); - -%Reduce matrix from n-size to f-size - -if gridset.ssize>0 - Kff= Knn(gridset.pv_f, gridset.pv_f); - Kfs= Knn(gridset.pv_f,gridset.pv_s); -else - Kff= Knn; - Kfs=[]; -end Index: sm/trunk/src/m/solutions/ice/IceReducerightside.m =================================================================== --- /issm/trunk/src/m/solutions/ice/IceReducerightside.m (revision 33) +++ (revision ) @@ -1,32 +1,0 @@ -function p_f= IceReducerightside( p_g, G_mn, K_fs, y_s) -%ICEREDUCERIGHTSIDE - reduce load vector from g to f -% -% Usage: -% p_f= IceReducerightside( p_g, G_mn, K_fs, y_s) - -global gridset - -% Reduce p_g to p_f -% Reduce p_g to p_n -if ( gridset.msize > 0 ) - p_n = p_g( gridset.pv_n, :); - p_m = p_g( gridset.pv_m, :); - p_n = p_n + G_mn' * p_m; -else - p_n = p_g; -end - -% Reduce p_n to p_f -if ( gridset.ssize > 0 ) - p_f = p_n( gridset.pv_f, :); -else - p_f = p_n; -end - -% for nonzero boundary conditions, subtract coupling forces, -% this operation is only executed if flag_y_s > 0, -% the flag is set in the calling script for linear analysis -% and the first iteration step in nonlinear analysis -if ~isempty(y_s) - p_f = p_f - K_fs * y_s; -end Index: /issm/trunk/src/m/solutions/ice/Normalizeconstraints.m =================================================================== --- /issm/trunk/src/m/solutions/ice/Normalizeconstraints.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/Normalizeconstraints.m (revision 34) @@ -0,0 +1,20 @@ +function G_mn=NormalizeConstraints(R_mg) +%NORMALIZECONSTRAINTS - create G_mn +% +% see NastranRgbDoc documentation in theoryguide directory. +% +% Usage: +% G_mn=NormalizeConstraints(R_mg) + +global gridset + +if(gridset.msize>0) + + R_mm = R_mg(:, gridset.pv_m); + R_mn = R_mg(:, gridset.pv_n); + + G_mn= - inv(R_mm) * R_mn; + +else + G_mn= []; +end Index: /issm/trunk/src/m/solutions/ice/PenaltyConstraints.m =================================================================== --- /issm/trunk/src/m/solutions/ice/PenaltyConstraints.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/PenaltyConstraints.m (revision 34) @@ -0,0 +1,25 @@ +function [loads,converged]=PenaltyConstraints(loads,grids,materials,params,inputs) +%PENALTYCONSTRAINTS - test if the thermal solution has converged +% +% This routine check if each element have converged in term of penalty activation +% once all the penalties remain the same between 2 iterations, return converged=1 +% +% Usage: +% [loads,converged]=PenaltyConstraints(loads,grids,materials,params,inputs) + +min_thermal_constraints=params.min_thermal_constraints; +unstable_count=0; + +%Go through constraints and set active flags +for n=1:length(loads), + [loads(n).load unstable]=PenaltyConstrain(loads(n).load,grids,materials,inputs); + unstable_count=unstable_count+unstable; +end + +disp(sprintf('%s %g',' number of unstable melting constraints: ',unstable_count)); + +if unstable_count<=min_thermal_constraints, + converged=1; +else + converged=0; +end Index: /issm/trunk/src/m/solutions/ice/Reducematrix.m =================================================================== --- /issm/trunk/src/m/solutions/ice/Reducematrix.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/Reducematrix.m (revision 34) @@ -0,0 +1,22 @@ +function Knn= Reducematrix( Kgg, Gmn, pv_m,pv_n,msize,flag) +%REDUCEMATRIX - reduce stiffness matrix from g to n +% +% Usage: +% Knn= Reducematrix( Kgg, Gmn, pv_m,pv_n,msize,flag) + +if msize>0, + Kmm=Kgg(pv_m,pv_m); + Kmn=Kgg(pv_m,pv_n); + Knm=Kgg(pv_n,pv_m); + Knn=Kgg(pv_n,pv_n); + + Knn= Knn + Knm * Gmn; + + if flag~=1, + Kmn = Kmn + Kmm * Gmn; + Knn = Knn + Gmn' * Kmn; + end + +else + Knn=Kgg; +end Index: /issm/trunk/src/m/solutions/ice/Reducematrixfromgtof.m =================================================================== --- /issm/trunk/src/m/solutions/ice/Reducematrixfromgtof.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/Reducematrixfromgtof.m (revision 34) @@ -0,0 +1,20 @@ +function [Kff Kfs]=Reducematrixfromgtof(Kgg,Gmn) +%REDUCEMATRIXFROMGTOF - reduce stiffness matrix from gset to fset +% +% Usage: +% [Kff Kfs]=Reducematrixfromgtof(Kgg,Gmn) + +global gridset + +%Reduce matrix from g-size to n-size +Knn= Reducematrix( Kgg, Gmn, gridset.pv_m, gridset.pv_n, gridset.msize, 2 ); + +%Reduce matrix from n-size to f-size + +if gridset.ssize>0 + Kff= Knn(gridset.pv_f, gridset.pv_f); + Kfs= Knn(gridset.pv_f,gridset.pv_s); +else + Kff= Knn; + Kfs=[]; +end Index: /issm/trunk/src/m/solutions/ice/Reducerightside.m =================================================================== --- /issm/trunk/src/m/solutions/ice/Reducerightside.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/Reducerightside.m (revision 34) @@ -0,0 +1,32 @@ +function p_f= Reducerightside( p_g, G_mn, K_fs, y_s) +%REDUCERIGHTSIDE - reduce load vector from g to f +% +% Usage: +% p_f= Reducerightside( p_g, G_mn, K_fs, y_s) + +global gridset + +% Reduce p_g to p_f +% Reduce p_g to p_n +if ( gridset.msize > 0 ) + p_n = p_g( gridset.pv_n, :); + p_m = p_g( gridset.pv_m, :); + p_n = p_n + G_mn' * p_m; +else + p_n = p_g; +end + +% Reduce p_n to p_f +if ( gridset.ssize > 0 ) + p_f = p_n( gridset.pv_f, :); +else + p_f = p_n; +end + +% for nonzero boundary conditions, subtract coupling forces, +% this operation is only executed if flag_y_s > 0, +% the flag is set in the calling script for linear analysis +% and the first iteration step in nonlinear analysis +if ~isempty(y_s) + p_f = p_f - K_fs * y_s; +end Index: /issm/trunk/src/m/solutions/ice/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic.m (revision 33) +++ /issm/trunk/src/m/solutions/ice/diagnostic.m (revision 34) @@ -1,24 +1,19 @@ -function md=icediagnostic(md); -%ICEDIAGNOSTIC - compute the velocity field of a model +function md=diagnostic(md); +%DIAGNOSTIC - compute the velocity field of a model % % this routine is a Matlab drivn solution % % Usage: -% u_g=icediagnostic_core_linear(m,analysis_type,varargin); +% u_g=diagnostic_core_linear(m,analysis_type,varargin); % % See also: ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR + +%start timing +t1=clock; %define global variables iceglobal -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('icediagnostic error message: parallel support not implemented yet'); -end - -%Create fem structure (input of icediagnostic3d) +%Create fem structure (input of diagnostic3d) fem=struct(); %Figure out which type of elements are present @@ -48,5 +43,5 @@ %compute solution - u_g=icediagnostic2d(md,fem); + u_g=diagnostic2d(md,fem); %Load results onto model @@ -83,7 +78,11 @@ %compute solution - u_g=icediagnostic3d(md,fem); + u_g=diagnostic3d(md,fem); %Load results onto model md=Loadresults(md,fem,u_g); end + +%stop timing +t2=clock; +disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: /issm/trunk/src/m/solutions/ice/diagnostic2d.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic2d.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/diagnostic2d.m (revision 34) @@ -0,0 +1,55 @@ +function u_g=diagnostic2d(md,fem), +%DIAGNOSTIC2D - diagnostic of a 2d model +% +% this routine computes the velocity field of a glacier in 2d +% +% Usage: +% u_g=diagnostic2d(md,fem) +% +% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR + +%recover inputs +inputs=fem.inputs; + +%Compute surface slope if there are hutter elements +if fem.ishutter, + + %Get field of fem + m_ss=fem.m_ss; + m_dhu=fem.m_dhu; + + %Compute slope of the surface + if fem.m_dhu.params.debug + disp(sprintf('\n%s',['computing surface slope (x and y derivatives)...'])); + end + slopex=diagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs); + slopey=diagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs); + + %initialize slope surface input for hutter solution sequence + slopesurface=zeros(m_ss.gridset.gsize,1); + slopesurface(1:6:m_ss.gridset.gsize,1)=slopex(1:6:end); + slopesurface(2:6:m_ss.gridset.gsize,1)=slopey(1:6:end); + + %Compute Hutter solution + disp(sprintf('\n%s',['computing hutter velocities...'])); + inputs.slopesurface=slopesurface; + u_g=diagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs); + + %Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa) + u_g(4:6:m_dhu.gridset.gsize)=1/2*md.rho_ice*md.g*md.thickness; +end + +if fem.ismacayealpattyn, + + %Get field of fem + m_dh=fem.m_dh; + + %Compute horizontal solution. + disp(sprintf('\n%s',['computing horizontal velocities...'])); + + %Run core solution + u_g=diagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs); + + %Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa) + u_g(4:6:m_dh.gridset.gsize)=1/2*md.rho_ice*md.g*md.thickness; +end Index: /issm/trunk/src/m/solutions/ice/diagnostic3d.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic3d.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/diagnostic3d.m (revision 34) @@ -0,0 +1,110 @@ +function u_g=diagnostic3d(md,fem), +%DIAGNOSTIC3D - diagnostic of a 3d model +% +% this routine computes the velocity field of a glacier in 3d +% +% Usage: +% u_g=diagnostic3d(md,fem) +% +% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR + +global gridset + +%recover inputs +inputs=fem.inputs; + +%Compute surface slope if there are hutter elements +if fem.ishutter, + m_ss=fem.m_ss; + m_dhu=fem.m_dhu; + + %Compute slope of the surface + if m_ss.params.debug + disp(sprintf('\n%s',['computing surface slope (x and y derivatives)...'])); + end + slopex=diagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs); + slopey=diagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs); + + %Project slopes computed at surface, onto the 3d volume + slopex=project3d(md,project2d(md,slopex(1:6:end),md.numlayers),'node'); + slopey=project3d(md,project2d(md,slopey(1:6:end),md.numlayers),'node'); + + %Feed hutter solution with the slope input + slopesurface=zeros(m_ss.gridset.gsize,1); + slopesurface(1:6:m_ss.gridset.gsize,1)=slopex; + slopesurface(2:6:m_ss.gridset.gsize,1)=slopey; + inputs.slopesurface=slopesurface; + + %Compute Hutter solution + disp(sprintf('\n%s',['computing hutter velocities...'])); + u_g=diagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs); +end + +if fem.ismacayealpattyn, + m_dh=fem.m_dh; + + %Compute horizontal solution. + disp(sprintf('\n%s',['computing horizontal velocities...'])); + + %Run core solution + u_g=diagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs); +end + +%Extrude velocities for collapsed penta elements +u_g=VelocityExtrude(md,u_g); + +%Compute depth averaged velocity and add it to inputs +gridset=fem.m_dv.gridset; +velocity_average=HorizontalVelocityDepthAverage(md,u_g); +inputs.velocity=u_g; +inputs.velocity_average=velocity_average; + +%Compute basal vertical velocity, in 3d +disp(sprintf('\n%s',['computing basal vertical velocity...'])); +u_g_basevert=diagnostic_core_linear(fem.m_dbv,'diagnostic_basevert',inputs); + +%update dirichlet boundary conditions with these new base vertical velocities +m_dv=fem.m_dv; +gridset=m_dv.gridset; +m_dv.y_g=u_g_basevert; m_dv.ys=Reducevector_g(m_dv.y_g); + +%compute vertical velocities +disp(sprintf('\n%s',['computing vertical velocity...'])); +u_g_vert=diagnostic_core_linear(m_dv,'diagnostic_vert',inputs); + +%add contribution to vertical velocity +u_g=u_g+u_g_vert; + +%Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa) +u_g(4:6:m_dv.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z)/md.stokesreconditioning; + +if fem.isstokes, + m_bs=fem.m_bs; + m_ds=fem.m_ds; + + %Compute slope of the bed + if m_bs.params.debug + disp(sprintf('\n%s',['computing bed slope (x and y derivatives)...'])); + end + slopex=diagnostic_core_linear(m_bs,'bed_slope_compute_x',inputs); + slopey=diagnostic_core_linear(m_bs,'bed_slope_compute_y',inputs); + + slopex=project3d(md,project2d(md,slopex(1:6:end),1),'node'); + slopey=project3d(md,project2d(md,slopey(1:6:end),1),'node'); + + slopebed(1:6:m_bs.gridset.gsize,1)=slopex; + slopebed(2:6:m_bs.gridset.gsize,1)=slopey; + inputs.slopebed=slopebed; + inputs.velocity=u_g; + + %update dirichlet boundary conditions with the velocities computed previously + gridset=m_ds.gridset; + m_ds.y_g=u_g; m_ds.ys=Reducevector_g(m_ds.y_g); + + %compute stokes velocities and pressure + disp(sprintf('\n%s',['computing stokes velocities and pressure ...'])); + u_g=diagnostic_core_nonlinear(m_ds,'diagnostic_stokes',inputs); +end + +%Reconditioning of the pressure +u_g(4:6:m_dv.gridset.gsize)=u_g(4:6:m_dv.gridset.gsize)*md.stokesreconditioning; Index: /issm/trunk/src/m/solutions/ice/diagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic_core_linear.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/diagnostic_core_linear.m (revision 34) @@ -0,0 +1,66 @@ +function u_g=diagnostic_core_linear(m,analysis_type,varargin); +%DIAGNOSTIC_CORE_LINEAR - core solution of linear problems +% +% Core of the vertical diagnostic solution. No need for a non-linear iteration. +% Otherwise similar to ICEDIAGNOSTIC_CORE_NONLINEAR +% +% Usage: +% u_g=diagnostic_core_linear(m,analysis_type,varargin); +% +% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR + +%global variables +global cluster gridset + +if nargin==3, + inputs=varargin{1}; +else + inputs=struct(); +end + +%recover fem model fields +elements=m.elements; +grids=m.grids; +materials=m.materials; +loads=m.loads; +ys=m.ys; +gridset=m.gridset; +G_mn=m.Gmn; +params=m.params; + + +%recover parameters +sparsity=params.sparsity; +solver_type=params.solver_type; +debug=params.debug; + +disp(sprintf('%s',' solving')); + +% Generate system matrices (stiffness and load) +kflag=1; pflag=1; + +%generate stiffness and loads +[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); +[K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); + +if cluster, + K_gg=distributed(gplus(K_gg),'convert'); + p_g=gplus(p_g); +end + +% Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), +% and compute modifications to loads from single point constraints. +[K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); + +% Reduce load from g set to f set +p_f=IceReducerightside(p_g,G_mn,K_fs,ys); + +% Solve +uf=Solver(K_ff,p_f,solver_type); +if debug, + disp(sprintf('%s%g',' condition number of stiffness matrix: ',condest(K_ff))); +end + + +% Add single point constraints back, ie increase f-set by s-set into the global g-set. +u_g= full(Mergesolution_g(uf,G_mn,ys)); Index: /issm/trunk/src/m/solutions/ice/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic_core_nonlinear.m (revision 34) +++ /issm/trunk/src/m/solutions/ice/diagnostic_core_nonlinear.m (revision 34) @@ -0,0 +1,161 @@ +function [u_g varargout]=diagnostic_core_nonlinear(m,analysis_type, varargin); +%DIAGNOSTIC_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]=diagnostic_core_nonlinear(m,analysis_type,varargin); +% +% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_LINEAR + +%global variables +global cluster gridset element_debug element_debugid + +%recover fem model fields +elements=m.elements; +grids=m.grids; +materials=m.materials; +loads=m.loads; +ys=m.ys; +gridset=m.gridset; +params=m.params; +G_mn=m.Gmn; + +%recover parameters +sparsity=params.sparsity; +solver_type=params.solver_type; +eps_abs=params.eps_abs; +eps_rel=params.eps_rel; +yts=params.yts; +debug=params.debug; +element_debug=params.element_debug; +element_debugid=params.element_debugid; + +%recover existing velocity if given in input and initialize solution +if nargin==3, + inputs=varargin{1}; + + [velocity_param velocity_is_present]=recover_input(inputs,'velocity'); + if velocity_is_present + soln(1).u_g=velocity_param; + soln(1).u_g(4:6:gridset.gsize)=0; %no pressure in the velocity + soln(2).u_g=velocity_param; + soln(2).u_g(4:6:gridset.gsize)=0; %no pressure in the velocity + else + soln(1).u_g={}; + soln(2).u_g={}; + end +else + velocity_is_present=0; + inputs=struct(); + soln(1).u_g={}; + soln(2).u_g={}; +end + + +%Initialization +count=2; +converged=0; + +if debug, + disp(sprintf('%s',' starting direct shooting method')); +end + +while(~converged), + + % Generate system matrices (stiffness and load) + %compute stiffness matrix flag + kflag=1; + %compute loads ? + if strcmp(analysis_type, 'diagnostic_stokes'), + pflag=1; + else + if count==2, pflag=1; else pflag=0; end + end + + %add velocity to inputs + inputs.velocity=soln(count).u_g; + inputs.oldvelocity=soln(count-1).u_g; + + %generate stiffness and loads + [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); + [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); + + %Save loads + if count==2, + p_g_old=p_g; + else + p_g=p_g_old; + end + + if cluster, + K_gg=distributed(gplus(K_gg),'convert'); + p_g=gplus(p_g); + end + + % Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), + % and compute modifications to loads from single point constraints. + + [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); + + % Reduce load from g set to f set + p_f=IceReducerightside(p_g,G_mn,K_fs,ys); + + % Solve + u_f=Solver(K_ff,p_f,solver_type); + if debug, + disp(sprintf('%s%g',' condition number of stiffness matrix: ',condest(K_ff))); + end + + %increment index + count=count+1; + + % Add single point constraints back, ie increase f-set by s-set into the global g-set. + soln(count).u_g= full(Mergesolution_g(u_f,G_mn, ys)); %make solution full + + %Figure out if convergence is reached. + if((count>=4)| velocity_is_present), %make sure we already iterated at least once. + + %compute relative velocity difference for this step. + relative_change=norm(soln(count).u_g-soln(count-1).u_g,2)/norm(soln(count-1).u_g,2); + + if relative_change ',eps_rel)); end + converged=0; + end + + if ~isnan(eps_abs) + %compute velocity difference for this step. + change=max(abs(soln(count).u_g-soln(count-1).u_g))*yts; + + if change ',eps_abs,'m/yr')); end + converged=0; + end + end + end +end + +%prepare output +u_g=soln(end).u_g; + +%more output might be needed, for ex when running in cielocontrol.m +nout=max(nargout,1)-1; +if nout==2, + %K_ff and K_fs are requested, for the case where we have no loads (ie basal drag) + inputs.drag=zeros(gridset.gsize,1); kflag=1;pflag=1; + inputs.velocity=soln(count).u_g; + + [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); + [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); + [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); + varargout(1)={K_ff}; + varargout(2)={K_fs}; +end Index: sm/trunk/src/m/solutions/ice/diagnostic_wrapper.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic_wrapper.m (revision 33) +++ (revision ) @@ -1,36 +1,0 @@ -function md=icediagnostic_wrapper(md); -%ICEDIAGNOSTIC_WRAPPER - wrapper of diagnostic solution -% -% diagnostic wrapper for diagnostic matlab driven solution. -% we need a wrapper so that we can launch the diagnostic solution using Matlab parallel, -% -% Usage: -% md=icediagnostic_wrapper(md) -% -% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR - -%start timing -t1=clock; - -%If running in parallel, launch core diagnostic routine within a parallel job. -if strcmpi(md.cluster,'yes'), - sched = findResource('scheduler', 'configuration', md.scheduler_configuration); - pjob = createParallelJob(sched); - set(pjob, 'FileDependencies', {'icediagnostic.m'}); - set(pjob, 'MaximumNumberOfWorkers', md.np); set(pjob, 'MinimumNumberOfWorkers', md.np); - t = createTask(pjob, @icediagnostic, 1, {md}); - set(t, 'CaptureCommandWindowOutput', true) - disp(' Submitting job to cluster'); - submit(pjob); - waitForState(pjob); - get(t, 'CommandWindowOutput') - get(t,'ErrorMessage') - disp(' Done running on cluster'); -else - md=icediagnostic(md); -end - -%stop timing -t2=clock; - -disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: /issm/trunk/src/m/solutions/ice/huttercompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/huttercompute.m (revision 33) +++ /issm/trunk/src/m/solutions/ice/huttercompute.m (revision 34) @@ -7,12 +7,4 @@ %define global variables iceglobal - -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('icediagnostic error message: parallel support not implemented yet'); -end %setup some parameters to be passed to the core solution @@ -45,5 +37,5 @@ slopesurface(2:6:gridset.gsize,1)=md.surface_slopey; disp(sprintf('\n%s',['computing Hutter''s velocities...'])); -u_g=icediagnostic_core_linear(elements,grids,materials,loads,params,ys,'diagnostic_hutter',struct('slopesurface',slopesurface)); +u_g=diagnostic_core_linear(elements,grids,materials,loads,params,ys,'diagnostic_hutter',struct('slopesurface',slopesurface)); md.vx=u_g(1:6:end)*md.yts; Index: sm/trunk/src/m/solutions/ice/icediagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic.m (revision 33) +++ (revision ) @@ -1,89 +1,0 @@ -function md=icediagnostic(md); -%ICEDIAGNOSTIC - compute the velocity field of a model -% -% this routine is a Matlab drivn solution -% -% Usage: -% u_g=icediagnostic_core_linear(m,analysis_type,varargin); -% -% See also: ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR - -%define global variables -iceglobal - -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('icediagnostic error message: parallel support not implemented yet'); -end - -%Create fem structure (input of icediagnostic3d) -fem=struct(); -%Figure out which type of elements are present -[fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type); - -if strcmpi(md.type,'2d'), - - %First, build elements,grids,loads, etc ... for horizontal model - if fem.ishutter, - fem.m_ss=CreateFemModel(md,'surface_slope_compute'); - fem.m_dhu=CreateFemModel(md,'diagnostic_hutter'); - end - if fem.ismacayealpattyn, - fem.m_dh=CreateFemModel(md,'diagnostic_horiz'); - end - - %plug inputs in fem (existing velocity) - if ~isnan(md.vx) & ~isnan(md.vy) & fem.ismacayealpattyn, - %an input velocity is present in the model, use it to bootstrap the diagnostic core nonlinear solution. - m_dh=fem.m_dh; - fem.inputs.velocity=zeros(m_dh.gridset.gsize,1); - fem.inputs.velocity(1:6:m_dh.gridset.gsize)=md.vx/md.yts; - fem.inputs.velocity(2:6:m_dh.gridset.gsize)=md.vy/md.yts; - else - fem.inputs=struct(); - end - - %compute solution - u_g=icediagnostic2d(md,fem); - - %Load results onto model - md=Loadresults(md,fem,u_g); - -else - - %First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model - fem.m_dbv=CreateFemModel(md,'diagnostic_basevert'); - fem.m_dv=CreateFemModel(md,'diagnostic_vert'); - - if fem.ismacayealpattyn, - fem.m_dh=CreateFemModel(md,'diagnostic_horiz'); - end - if fem.ishutter, - fem.m_ss=CreateFemModel(md,'surface_slope_compute'); - fem.m_dhu=CreateFemModel(md,'diagnostic_hutter'); - end - if fem.isstokes, - fem.m_bs=CreateFemModel(md,'bed_slope_compute'); - fem.m_ds=CreateFemModel(md,'diagnostic_stokes'); - end - - %plug inputs in fem (existing velocity) - if ~isnan(md.vx) & ~isnan(md.vy) & fem.ismacayealpattyn, - %an input velocity is present in the model, use it to bootstrap the diagnostic core nonlinear solution. - m_dh=fem.m_dh; - fem.inputs.velocity=zeros(m_dh.gridset.gsize,1); - fem.inputs.velocity(1:6:m_dh.gridset.gsize)=md.vx/md.yts; - fem.inputs.velocity(2:6:m_dh.gridset.gsize)=md.vy/md.yts; - else - fem.inputs=struct(); - end - - %compute solution - u_g=icediagnostic3d(md,fem); - - %Load results onto model - md=Loadresults(md,fem,u_g); -end Index: sm/trunk/src/m/solutions/ice/icediagnostic2d.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic2d.m (revision 33) +++ (revision ) @@ -1,55 +1,0 @@ -function u_g=icediagnostic2d(md,fem), -%ICEDIAGNOSTIC2D - diagnostic of a 2d model -% -% this routine computes the velocity field of a glacier in 2d -% -% Usage: -% u_g=icediagnostic2d(md,fem) -% -% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR - -%recover inputs -inputs=fem.inputs; - -%Compute surface slope if there are hutter elements -if fem.ishutter, - - %Get field of fem - m_ss=fem.m_ss; - m_dhu=fem.m_dhu; - - %Compute slope of the surface - if fem.m_dhu.params.debug - disp(sprintf('\n%s',['computing surface slope (x and y derivatives)...'])); - end - slopex=icediagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs); - slopey=icediagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs); - - %initialize slope surface input for hutter solution sequence - slopesurface=zeros(m_ss.gridset.gsize,1); - slopesurface(1:6:m_ss.gridset.gsize,1)=slopex(1:6:end); - slopesurface(2:6:m_ss.gridset.gsize,1)=slopey(1:6:end); - - %Compute Hutter solution - disp(sprintf('\n%s',['computing hutter velocities...'])); - inputs.slopesurface=slopesurface; - u_g=icediagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs); - - %Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa) - u_g(4:6:m_dhu.gridset.gsize)=1/2*md.rho_ice*md.g*md.thickness; -end - -if fem.ismacayealpattyn, - - %Get field of fem - m_dh=fem.m_dh; - - %Compute horizontal solution. - disp(sprintf('\n%s',['computing horizontal velocities...'])); - - %Run core solution - u_g=icediagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs); - - %Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa) - u_g(4:6:m_dh.gridset.gsize)=1/2*md.rho_ice*md.g*md.thickness; -end Index: sm/trunk/src/m/solutions/ice/icediagnostic3d.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic3d.m (revision 33) +++ (revision ) @@ -1,110 +1,0 @@ -function u_g=icediagnostic3d(md,fem), -%ICEDIAGNOSTIC3D - diagnostic of a 3d model -% -% this routine computes the velocity field of a glacier in 3d -% -% Usage: -% u_g=icediagnostic3d(md,fem) -% -% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR - -global gridset - -%recover inputs -inputs=fem.inputs; - -%Compute surface slope if there are hutter elements -if fem.ishutter, - m_ss=fem.m_ss; - m_dhu=fem.m_dhu; - - %Compute slope of the surface - if m_ss.params.debug - disp(sprintf('\n%s',['computing surface slope (x and y derivatives)...'])); - end - slopex=icediagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs); - slopey=icediagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs); - - %Project slopes computed at surface, onto the 3d volume - slopex=project3d(md,project2d(md,slopex(1:6:end),md.numlayers),'node'); - slopey=project3d(md,project2d(md,slopey(1:6:end),md.numlayers),'node'); - - %Feed hutter solution with the slope input - slopesurface=zeros(m_ss.gridset.gsize,1); - slopesurface(1:6:m_ss.gridset.gsize,1)=slopex; - slopesurface(2:6:m_ss.gridset.gsize,1)=slopey; - inputs.slopesurface=slopesurface; - - %Compute Hutter solution - disp(sprintf('\n%s',['computing hutter velocities...'])); - u_g=icediagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs); -end - -if fem.ismacayealpattyn, - m_dh=fem.m_dh; - - %Compute horizontal solution. - disp(sprintf('\n%s',['computing horizontal velocities...'])); - - %Run core solution - u_g=icediagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs); -end - -%Extrude velocities for collapsed penta elements -u_g=VelocityExtrude(md,u_g); - -%Compute depth averaged velocity and add it to inputs -gridset=fem.m_dv.gridset; -velocity_average=HorizontalVelocityDepthAverage(md,u_g); -inputs.velocity=u_g; -inputs.velocity_average=velocity_average; - -%Compute basal vertical velocity, in 3d -disp(sprintf('\n%s',['computing basal vertical velocity...'])); -u_g_basevert=icediagnostic_core_linear(fem.m_dbv,'diagnostic_basevert',inputs); - -%update dirichlet boundary conditions with these new base vertical velocities -m_dv=fem.m_dv; -gridset=m_dv.gridset; -m_dv.y_g=u_g_basevert; m_dv.ys=Reducevector_g(m_dv.y_g); - -%compute vertical velocities -disp(sprintf('\n%s',['computing vertical velocity...'])); -u_g_vert=icediagnostic_core_linear(m_dv,'diagnostic_vert',inputs); - -%add contribution to vertical velocity -u_g=u_g+u_g_vert; - -%Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa) -u_g(4:6:m_dv.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z)/md.stokesreconditioning; - -if fem.isstokes, - m_bs=fem.m_bs; - m_ds=fem.m_ds; - - %Compute slope of the bed - if m_bs.params.debug - disp(sprintf('\n%s',['computing bed slope (x and y derivatives)...'])); - end - slopex=icediagnostic_core_linear(m_bs,'bed_slope_compute_x',inputs); - slopey=icediagnostic_core_linear(m_bs,'bed_slope_compute_y',inputs); - - slopex=project3d(md,project2d(md,slopex(1:6:end),1),'node'); - slopey=project3d(md,project2d(md,slopey(1:6:end),1),'node'); - - slopebed(1:6:m_bs.gridset.gsize,1)=slopex; - slopebed(2:6:m_bs.gridset.gsize,1)=slopey; - inputs.slopebed=slopebed; - inputs.velocity=u_g; - - %update dirichlet boundary conditions with the velocities computed previously - gridset=m_ds.gridset; - m_ds.y_g=u_g; m_ds.ys=Reducevector_g(m_ds.y_g); - - %compute stokes velocities and pressure - disp(sprintf('\n%s',['computing stokes velocities and pressure ...'])); - u_g=icediagnostic_core_nonlinear(m_ds,'diagnostic_stokes',inputs); -end - -%Reconditioning of the pressure -u_g(4:6:m_dv.gridset.gsize)=u_g(4:6:m_dv.gridset.gsize)*md.stokesreconditioning; Index: sm/trunk/src/m/solutions/ice/icediagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic_core_linear.m (revision 33) +++ (revision ) @@ -1,66 +1,0 @@ -function u_g=icediagnostic_core_linear(m,analysis_type,varargin); -%ICEDIAGNOSTIC_CORE_LINEAR - core solution of linear problems -% -% Core of the vertical diagnostic solution. No need for a non-linear iteration. -% Otherwise similar to ICEDIAGNOSTIC_CORE_NONLINEAR -% -% Usage: -% u_g=icediagnostic_core_linear(m,analysis_type,varargin); -% -% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR - -%global variables -global cluster gridset - -if nargin==3, - inputs=varargin{1}; -else - inputs=struct(); -end - -%recover fem model fields -elements=m.elements; -grids=m.grids; -materials=m.materials; -loads=m.loads; -ys=m.ys; -gridset=m.gridset; -G_mn=m.Gmn; -params=m.params; - - -%recover parameters -sparsity=params.sparsity; -solver_type=params.solver_type; -debug=params.debug; - -disp(sprintf('%s',' solving')); - -% Generate system matrices (stiffness and load) -kflag=1; pflag=1; - -%generate stiffness and loads -[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); -[K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); - -if cluster, - K_gg=distributed(gplus(K_gg),'convert'); - p_g=gplus(p_g); -end - -% Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), -% and compute modifications to loads from single point constraints. -[K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); - -% Reduce load from g set to f set -p_f=IceReducerightside(p_g,G_mn,K_fs,ys); - -% Solve -uf=Solver(K_ff,p_f,solver_type); -if debug, - disp(sprintf('%s%g',' condition number of stiffness matrix: ',condest(K_ff))); -end - - -% Add single point constraints back, ie increase f-set by s-set into the global g-set. -u_g= full(Mergesolution_g(uf,G_mn,ys)); Index: sm/trunk/src/m/solutions/ice/icediagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic_core_nonlinear.m (revision 33) +++ (revision ) @@ -1,161 +1,0 @@ -function [u_g varargout]=icediagnostic_core_nonlinear(m,analysis_type, varargin); -%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 - -%global variables -global cluster gridset element_debug element_debugid - -%recover fem model fields -elements=m.elements; -grids=m.grids; -materials=m.materials; -loads=m.loads; -ys=m.ys; -gridset=m.gridset; -params=m.params; -G_mn=m.Gmn; - -%recover parameters -sparsity=params.sparsity; -solver_type=params.solver_type; -eps_abs=params.eps_abs; -eps_rel=params.eps_rel; -yts=params.yts; -debug=params.debug; -element_debug=params.element_debug; -element_debugid=params.element_debugid; - -%recover existing velocity if given in input and initialize solution -if nargin==3, - inputs=varargin{1}; - - [velocity_param velocity_is_present]=recover_input(inputs,'velocity'); - if velocity_is_present - soln(1).u_g=velocity_param; - soln(1).u_g(4:6:gridset.gsize)=0; %no pressure in the velocity - soln(2).u_g=velocity_param; - soln(2).u_g(4:6:gridset.gsize)=0; %no pressure in the velocity - else - soln(1).u_g={}; - soln(2).u_g={}; - end -else - velocity_is_present=0; - inputs=struct(); - soln(1).u_g={}; - soln(2).u_g={}; -end - - -%Initialization -count=2; -converged=0; - -if debug, - disp(sprintf('%s',' starting direct shooting method')); -end - -while(~converged), - - % Generate system matrices (stiffness and load) - %compute stiffness matrix flag - kflag=1; - %compute loads ? - if strcmp(analysis_type, 'diagnostic_stokes'), - pflag=1; - else - if count==2, pflag=1; else pflag=0; end - end - - %add velocity to inputs - inputs.velocity=soln(count).u_g; - inputs.oldvelocity=soln(count-1).u_g; - - %generate stiffness and loads - [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); - [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); - - %Save loads - if count==2, - p_g_old=p_g; - else - p_g=p_g_old; - end - - if cluster, - K_gg=distributed(gplus(K_gg),'convert'); - p_g=gplus(p_g); - end - - % Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), - % and compute modifications to loads from single point constraints. - - [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); - - % Reduce load from g set to f set - p_f=IceReducerightside(p_g,G_mn,K_fs,ys); - - % Solve - u_f=Solver(K_ff,p_f,solver_type); - if debug, - disp(sprintf('%s%g',' condition number of stiffness matrix: ',condest(K_ff))); - end - - %increment index - count=count+1; - - % Add single point constraints back, ie increase f-set by s-set into the global g-set. - soln(count).u_g= full(Mergesolution_g(u_f,G_mn, ys)); %make solution full - - %Figure out if convergence is reached. - if((count>=4)| velocity_is_present), %make sure we already iterated at least once. - - %compute relative velocity difference for this step. - relative_change=norm(soln(count).u_g-soln(count-1).u_g,2)/norm(soln(count-1).u_g,2); - - if relative_change ',eps_rel)); end - converged=0; - end - - if ~isnan(eps_abs) - %compute velocity difference for this step. - change=max(abs(soln(count).u_g-soln(count-1).u_g))*yts; - - if change ',eps_abs,'m/yr')); end - converged=0; - end - end - end -end - -%prepare output -u_g=soln(end).u_g; - -%more output might be needed, for ex when running in cielocontrol.m -nout=max(nargout,1)-1; -if nout==2, - %K_ff and K_fs are requested, for the case where we have no loads (ie basal drag) - inputs.drag=zeros(gridset.gsize,1); kflag=1;pflag=1; - inputs.velocity=soln(count).u_g; - - [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); - [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); - [K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); - varargout(1)={K_ff}; - varargout(2)={K_fs}; -end Index: sm/trunk/src/m/solutions/ice/icediagnostic_wrapper.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic_wrapper.m (revision 33) +++ (revision ) @@ -1,36 +1,0 @@ -function md=icediagnostic_wrapper(md); -%ICEDIAGNOSTIC_WRAPPER - wrapper of diagnostic solution -% -% diagnostic wrapper for diagnostic matlab driven solution. -% we need a wrapper so that we can launch the diagnostic solution using Matlab parallel, -% -% Usage: -% md=icediagnostic_wrapper(md) -% -% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR - -%start timing -t1=clock; - -%If running in parallel, launch core diagnostic routine within a parallel job. -if strcmpi(md.cluster,'yes'), - sched = findResource('scheduler', 'configuration', md.scheduler_configuration); - pjob = createParallelJob(sched); - set(pjob, 'FileDependencies', {'icediagnostic.m'}); - set(pjob, 'MaximumNumberOfWorkers', md.np); set(pjob, 'MinimumNumberOfWorkers', md.np); - t = createTask(pjob, @icediagnostic, 1, {md}); - set(t, 'CaptureCommandWindowOutput', true) - disp(' Submitting job to cluster'); - submit(pjob); - waitForState(pjob); - get(t, 'CommandWindowOutput') - get(t,'ErrorMessage') - disp(' Done running on cluster'); -else - md=icediagnostic(md); -end - -%stop timing -t2=clock; - -disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: sm/trunk/src/m/solutions/ice/icemelting_core.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icemelting_core.m (revision 33) +++ (revision ) @@ -1,55 +1,0 @@ -function melting_g=icemelting_core(m,analysis_type,varargin); -%ICEMELTING_CORE - core solution of melting computation -% -% Usage: -% melting_g=icemelting_core(m,analysis_type,varargin) - -%global variables -global cluster gridset - -%recover fem model fields -elements=m.elements; -grids=m.grids; -materials=m.materials; -loads=m.loads; -ys=m.ys; -gridset=m.gridset; -G_mn=m.Gmn; -params=m.params; - -%recover parameters -sparsity=params.sparsity; -solver_type=params.solver_type; - -%recover extra inputs and initialize solution -if nargin==3, - inputs=varargin{1}; -else - inputs=struct(); -end - -%initialize solution vector -melting_g=[]; - -disp(sprintf('%s',' solving')); - -% Generate system matrices (stiffness and load) -%stiffness is computed at each step, not the loads. -kflag=1;pflag=1; - -%generate stiffness and loads -[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, params.sparsity,inputs,'melting'); -[K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, params.sparsity,inputs,'melting',K_gg,p_g); - -% Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), -% and compute modifications to loads from single point constraints. -[K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); - -% Reduce load from g set to f set -p_f=IceReducerightside(p_g,G_mn,K_fs,ys); - -% Solve -melting_f=Solver(K_ff,p_f,params.solver_type); - -% Add single point constraints back, ie increase f-set by s-set into the global g-set. -melting_g= Mergesolution_g(melting_f,G_mn,ys); Index: sm/trunk/src/m/solutions/ice/iceparameters.m =================================================================== --- /issm/trunk/src/m/solutions/ice/iceparameters.m (revision 33) +++ (revision ) @@ -1,104 +1,0 @@ -function md=iceparameters(md); -%ICEPARAMETERS - compute parameters -% -% This routine is caled by the solver to compute the parameters -% that are given in md.outputparameters (stress, strain rate,...) -% -% Usage: -% md=iceparameters(md) - -%define global variables -iceglobal - -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('iceparameters error message: parallel support not implemented yet'); -end - -if strcmpi(md.type,'2d'), - %First, build elements,grids,loads, etc ... for thermal and melting - m_dh=CreateFemModel(md,'diagnostic_horiz'); - - %Are any more computations requested? - if ~isempty(md.parameteroutput), - disp(sprintf('\n%s',['computing outputs...'])); - %velocity field - u_g=zeros(gridset.gsize,1); - u_g(1:6:gridset.gsize)=md.vx/md.yts; %from m/a to m/s - u_g(2:6:gridset.gsize)=md.vy/md.yts; - - %Go through fields of parameteroutput and compute - for i=1:length(md.parameteroutput), - if strcmpi(md.parameteroutput{i},'viscousheating'), - disp(sprintf('%s',[' viscousheating'])); - md.viscousheating=ViscousHeatingCompute(m_dh,struct('velocity',u_g),md.type); - elseif strcmpi(md.parameteroutput{i},'pressure_elem'), - disp(sprintf('%s',[' pressure_elem'])); - md.pressure_elem=PressureElemCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'strainrate'), - disp(sprintf('%s',[' strainrate'])); - md.strainrate=StrainRateCompute(m_dh,struct('velocity',u_g),md.type); - elseif strcmpi(md.parameteroutput{i},'stress'), - disp(sprintf('%s',[' stress'])); - md.stress=StressCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'stress_bed'), - disp(sprintf('%s',[' stress_bed'])); - md.stress_bed=StressBedCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'stress_surface'), - disp(sprintf('%s',[' stress_surface'])); - md.stress_surface=StressSurfaceCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'deviatoricstress'), - disp(sprintf('%s',[' deviatoricstress'])); - md.deviatoricstress=DeviatoricStressCompute(m_dh,struct('velocity',u_g),md.type); - end - end - end - -else - - %First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model - m_dv=CreateFemModel(md,'diagnostic_vert'); - - %Are any more computations requested? - if ~isempty(md.parameteroutput), - disp(sprintf('\n%s',['computing outputs...'])); - %velocity field - u_g=zeros(gridset.gsize,1); - u_g(1:6:gridset.gsize)=md.vx/md.yts; %from m/a to m/s - u_g(2:6:gridset.gsize)=md.vy/md.yts; - u_g(3:6:gridset.gsize)=md.vz/md.yts; - if ~isnan(md.pressure), - u_g(4:6:gridset.gsize)=md.pressure; - end - - %Go through fields of parameteroutput and compute - for i=1:length(md.parameteroutput), - if strcmpi(md.parameteroutput{i},'viscousheating'), - disp(sprintf('%s',[' viscousheating'])); - md.viscousheating=ViscousHeatingCompute(m_dv,struct('velocity',u_g),md.type); - elseif strcmpi(md.parameteroutput{i},'pressure_elem'), - disp(sprintf('%s',[' pressure_elem'])); - md.pressure_elem=PressureElemCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'strainrate'), - disp(sprintf('%s',[' strainrate'])); - md.strainrate=StrainRateCompute(m_dv,struct('velocity',u_g),md.type); - elseif strcmpi(md.parameteroutput{i},'stress'), - disp(sprintf('%s',[' stress'])); - md.stress=StressCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'stress_bed'), - disp(sprintf('%s',[' stress_bed'])); - md.stress_bed=StressBedCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'stress_surface'), - disp(sprintf('%s',[' stress_surface'])); - md.stress_surface=StressSurfaceCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type); - elseif strcmpi(md.parameteroutput{i},'deviatoricstress'), - disp(sprintf('%s',[' deviatoricstress'])); - md.deviatoricstress=DeviatoricStressCompute(m_dv,struct('velocity',u_g),md.type); - end - end - end - -end Index: sm/trunk/src/m/solutions/ice/iceparameters_wrapper.m =================================================================== --- /issm/trunk/src/m/solutions/ice/iceparameters_wrapper.m (revision 33) +++ (revision ) @@ -1,34 +1,0 @@ -function md=iceparameters_wrapper(md); -%ICEPARAMETERS_WRAPPER - wrapper of the parameters solution sequence -% -% parameters wrapper for parameters matlab driven solution. -% we need a wrapper so that we can launch the parameters solution using Matlab parallel, -% -% Usage: -% md=iceparameters_wrapper(md) - -%start timing -t1=clock; - -%If running in parallel, launch core parameters routine within a parallel job. -if strcmpi(md.cluster,'yes'), - sched = findResource('scheduler', 'configuration', md.scheduler_configuration); - pjob = createParallelJob(sched); - set(pjob, 'FileDependencies', {'iceparameters.m'}); - set(pjob, 'MaximumNumberOfWorkers', md.np); set(pjob, 'MinimumNumberOfWorkers', md.np); - t = createTask(pjob, @iceparameters, 1, {md}); - set(t, 'CaptureCommandWindowOutput', true) - disp(' Submitting job to cluster'); - submit(pjob); - waitForState(pjob); - get(t, 'CommandWindowOutput') - get(t,'ErrorMessage') - disp(' Done running on cluster'); -else - md=iceparameters(md); -end - -%stop timing -t2=clock; - -disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: sm/trunk/src/m/solutions/ice/iceprognostic.m =================================================================== --- /issm/trunk/src/m/solutions/ice/iceprognostic.m (revision 33) +++ (revision ) @@ -1,70 +1,0 @@ -function md=iceprognostic(md) -%ICEPROGNOSTIC - prognostic solution -% -% This routine is used to compute the evolution of the geometry -% of a model -% -% Usage: -% md=iceprognostic(md) - -%define global variables -iceglobal - -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('iceprognostic error message: parallel support not implemented yet'); -end - -%setup some parameters to be passed to the core solution -params.sparsity=md.sparsity; -params.solver_type=md.solver_type; -params.eps_rel=md.eps_rel; -params.eps_abs=md.eps_abs; -params.debug=md.debug; - -%First, build elements,grids,loads, etc ... for prognostic -m=CreateFemModel(md,'prognostic'); - -%compute depth averaged horizontal velocity -u_g=zeros(gridset.gsize,1); -u_g(1:6:end)=md.vx/md.yts; %from m/yr to m/s -u_g(2:6:end)=md.vy/md.yts; %from m/yr to m/s -u_g(3:6:end)=md.vz/md.yts; %from m/yr to m/s - -if strcmpi(md.type,'3d'), - velocity_average=HorizontalVelocityDepthAverage(md,u_g); -else - %the velocity is already depthaverdged - velocity_average=u_g; -end - -%move surface and bed velocities to first lower layer. This will be needed by the prognostic model -melting=zeros(gridset.gsize,1); -accumulation=zeros(gridset.gsize,1); -if strcmpi(md.type,'3d'), - melting(1:6:end)=ShiftLayers(md,md.melting,1,1); - accumulation(1:6:end)=ShiftLayers(md,md.accumulation,md.numlayers,1); -else - melting(1:6:end)=md.melting/md.yts; %from m/yr to m/s - accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s -end - -%prepare inputs -thickness=zeros(gridset.gsize,1);thickness(1:6:end)=md.thickness; -%inputs=struct('thickness',thickness,'melting',melting,'accumulation',accumulation,... -%'surface_vertical_velocity',ws,'basal_vertical_velocity',wb,'dt',md.dt,'velocity_average',velocity_average); -inputs=struct('thickness',thickness,'melting',melting,'accumulation',accumulation,'dt',md.dt,'velocity_average',velocity_average); - -%Run core solution -h_g=iceprognostic_core(m,'prognostic',inputs); - -%project collapsed thickness onto 3d mesh -indx=1:6:m.gridset.gsize; indx=indx(m.tpart); -if strcmpi(md.type,'3d'), - md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node'); -else - md.new_thickness=h_g(indx); -end Index: sm/trunk/src/m/solutions/ice/iceprognostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/ice/iceprognostic_core.m (revision 33) +++ (revision ) @@ -1,58 +1,0 @@ -function h_g=iceprognostic_core(m,analysis_type,varargin) -%ICEPROGNOSTIC_CORE - core of prognostic solution -% -% This routine is used to compute the new thickness of a model after a time step ∆t -% -% Usage: -% h_g=iceprognostic_core(m,analysis_type,varargin) - -%global variables -global cluster gridset - -%recover fem model fields -elements=m.elements; -grids=m.grids; -materials=m.materials; -loads=m.loads; -ys=m.ys; -gridset=m.gridset; -G_mn=m.Gmn; -params=m.params; - -%recover parameters -sparsity=params.sparsity; -solver_type=params.solver_type; - -%recover extra inputs and initialize solution -if nargin==3, - inputs=varargin{1}; -else - inputs=struct(); -end - -disp(sprintf('%s',' solving')); - -% Generate system matrices (stiffness and load) -kflag=1; pflag=1; - -%generate stiffness and loads -[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type); -[K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); - -if cluster, - K_gg=distributed(gplus(K_gg),'convert'); - p_g=gplus(p_g); -end - -% Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), -% and compute modifications to loads from single point constraints. -[K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); - -% Reduce load from g set to f set -p_f=IceReducerightside(p_g,G_mn,K_fs,ys); - -% Solve -uf=Solver(K_ff,p_f,solver_type); - -% Add single point constraints back, ie increase f-set by s-set into the global g-set. -h_g= full(Mergesolution_g(uf,G_mn,ys)); Index: sm/trunk/src/m/solutions/ice/iceprognostic_wrapper.m =================================================================== --- /issm/trunk/src/m/solutions/ice/iceprognostic_wrapper.m (revision 33) +++ (revision ) @@ -1,34 +1,0 @@ -function md=iceprognostic_wrapper(md); -%ICEPROGNOSTIC_WRAPPER - wrapper of the prognostic solution -% -% prognostic wrapper for prognostic matlab driven solution. -% we need a wrapper so that we can launch the prognostic solution using Matlab parallel, -% -% Usage: -% md=iceprognostic_wrapper(md); - -%start timing -t1=clock; - -%If running in parallel, launch core prognostic routine within a parallel job. -if strcmpi(md.cluster,'yes'), - sched = findResource('scheduler', 'configuration', md.scheduler_configuration); - pjob = createParallelJob(sched); - set(pjob, 'FileDependencies', {'iceprognostic.m'}); - set(pjob, 'MaximumNumberOfWorkers', md.np); set(pjob, 'MinimumNumberOfWorkers', md.np); - t = createTask(pjob, @iceprognostic, 1, {md}); - set(t, 'CaptureCommandWindowOutput', true) - disp(' Submitting job to cluster'); - submit(pjob); - waitForState(pjob); - get(t, 'CommandWindowOutput') - get(t,'ErrorMessage') - disp(' Done running on cluster'); -else - md=iceprognostic(md); -end - -%stop timing -t2=clock; - -disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: sm/trunk/src/m/solutions/ice/icethermal.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icethermal.m (revision 33) +++ (revision ) @@ -1,97 +1,0 @@ -function md=icethermal(md,analysis_type) -%ICETHERMAL - thermal matlab driven solution -% -% Usage: -% md=icethermal(md,analysis_type) - -%define global variables -iceglobal - -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('icethermal error message: parallel support not implemented yet'); -end - -%First, build elements,grids,loads, etc ... for thermal and melting -m_t=CreateFemModel(md,analysis_type); -if strcmpi(analysis_type,'thermalsteady'), - m_m=CreateFemModel(md,'meltingsteady'); -else - m_m=CreateFemModel(md,'meltingtransient'); -end -%initialize velocity,pressure and material properties -velocity=zeros(m_t.gridset.gsize,1); -velocity(1:6:m_t.gridset.gsize)=md.vx(m_t.part)/md.yts; -velocity(2:6:m_t.gridset.gsize)=md.vy(m_t.part)/md.yts; -velocity(3:6:m_t.gridset.gsize)=md.vz(m_t.part)/md.yts; - -pressure=zeros(m_t.gridset.gsize,1); -pressure(1:6:m_t.gridset.gsize)=md.pressure; - -if strcmpi(analysis_type,'thermalsteady'), - - %Call core thermal computation - disp(' computing temperature...'); - gridset=m_t.gridset; - inputs=struct('pressure',pressure,'velocity',velocity,'dt',md.dt); - [t_g loads_t melting_offset]=icethermal_core(m_t,analysis_type,inputs); - - %Call core melting computation - disp(' computing melting...'); - gridset=m_m.gridset; - inputs=struct('pressure',pressure,'temperature',t_g,'melting_offset',melting_offset); - melting_g=icemelting_core(m_m,analysis_type,inputs); - - %Wrap up - indx=1:6:gridset.gsize; - indx=indx(m_t.tpart); - md.temperature=t_g(indx); - md.melting=melting_g(indx)*md.yts; %from m/s to m/a -else - - nsteps=md.ndt/md.dt; - - %initialize temperature and melting - temperature=zeros(gridset.gsize,1); - melting=zeros(gridset.gsize,1); - temperature(1:6:gridset.gsize)=md.temperature; - melting(1:6:gridset.gsize)=md.melting/md.yts; %in m/s - soln.t_g=temperature; - soln.melting_g=melting; - soln.time=0; - - - for n=1:nsteps, - - soln(n+1).time=(n)*md.dt; - disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt)); - - %Call core thermal computation - disp(' computing temperature...'); - gridset=m_t.gridset; - inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'velocity',velocity,'dt',md.dt); - [soln(n+1).t_g loads_t melting_offset]=icethermal_core(m_t,analysis_type,inputs); - - %Call core melting computation - disp(' computing melting...'); - gridset=m_m.gridset; - inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'melting_offset',melting_offset,'dt',md.dt); - soln(n+1).melting_g=icemelting_core(m_m,analysis_type,inputs); - - end - - %Wrap up - indx=1:6:m_t.gridset.gsize; - indx=indx(m_t.tpart); - solution=struct('time',{},'temperature',{},'melting',{}); - - for n=1:nsteps+1, - solution(n).temperature=soln(n).t_g(indx); - solution(n).melting=soln(n).melting_g(indx)*md.yts; %in m/a - solution(n).time=soln(n).time/md.yts; %in year - end - md.thermaltransient_results=solution; -end Index: sm/trunk/src/m/solutions/ice/icethermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icethermal_core.m (revision 33) +++ (revision ) @@ -1,73 +1,0 @@ -function [t_g loads melting_offset]=icethermal_core(m,analysis_type,varargin); -%ICETHERMAL_CORE - core of the thermal solution -% -% Usage: -% [t_g loads melting_offset]=icethermal_core(m,analysis_type,varargin) - -%global variables -global cluster gridset - -%recover fem model fields -elements=m.elements; -grids=m.grids; -materials=m.materials; -loads=m.loads; -ys=m.ys; -gridset=m.gridset; -G_mn=m.Gmn; -params=m.params; - -%recover parameters -sparsity=params.sparsity; -solver_type=params.solver_type; - -%recover extra inputs and initialize solution -if nargin==3, - inputs=varargin{1}; -else - inputs=struct(); -end - -%initialize solution vector -t_g=[]; - -converged=0; -firstiteration=1; - -disp(sprintf('%s',' starting direct shooting method')); - -dt=inputs.dt; - -%Recover pressure from inputs -pressure=inputs.pressure; - -while(~converged), - - % Generate system matrices (stiffness and load) - - %loads and stiffness are computed only once, not the penalties - kflag=1;pflag=1; - - %generate stiffness and loads - if firstiteration - [K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, params.sparsity,inputs,analysis_type); - firstiteration=0; - end - [K_gg_penalized,p_g_penalized melting_offset]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, params.sparsity,inputs,analysis_type,K_gg,p_g); - - % Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), - % and compute modifications to loads from single point constraints. - [K_ff K_fs]=IceReducematrixfromgtof(K_gg_penalized,G_mn); - - % Reduce load from g set to f set - p_f=IceReducerightside(p_g_penalized,G_mn,K_fs,ys); - - % Solve - t_f=Solver(K_ff,p_f,params.solver_type); - - % Add single point constraints back, ie increase f-set by s-set into the global g-set. - t_g= Mergesolution_g(t_f,G_mn,ys); - - [loads,converged]=IcePenaltyConstraints(loads,grids,materials,params,struct('pressure',pressure, 'dt',dt, 'temperature',t_g)); - -end Index: sm/trunk/src/m/solutions/ice/icethermal_wrapper.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icethermal_wrapper.m (revision 33) +++ (revision ) @@ -1,34 +1,0 @@ -function md=icethermal_wrapper(md,analysis_type); -%ICETHERMAL_WRAPPER - wrapper of the thermal solution -% -% thermal wrapper for thermal matlab driven solution. -% we need a wrapper so that we can launch the thermal solution using Matlab parallel, -% -% Usage: -% md=icethermal_wrapper(md,analysis_type); - -%start timing -t1=clock; - -%If running in parallel, launch core thermal routine within a parallel job. -if strcmpi(md.cluster,'yes'), - sched = findResource('scheduler', 'configuration', md.scheduler_configuration); - pjob = createParallelJob(sched); - set(pjob, 'FileDependencies', {'icethermal.m'}); - set(pjob, 'MaximumNumberOfWorkers', md.np); set(pjob, 'MinimumNumberOfWorkers', md.np); - t = createTask(pjob, @icethermal, 1, {md}); - set(t, 'CaptureCommandWindowOutput', true) - disp(' Submitting job to cluster'); - submit(pjob); - waitForState(pjob); - get(t, 'CommandWindowOutput') - get(t,'ErrorMessage') - disp(' Done running on cluster'); -else - md=icethermal(md,analysis_type); -end - -%stop timing -t2=clock; - -disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: sm/trunk/src/m/solutions/ice/icetransient.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icetransient.m (revision 33) +++ (revision ) @@ -1,22 +1,0 @@ -function md=icetransient(md); -%ICETRANSIENT - transient solution -% -% Usage: -% md=icetransient(md) - -%define global variables -iceglobal - -%determine if run is parallel -if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; - -%for now, only serial support is in -if cluster, - error('icetransient error message: parallel support not implemented yet'); -end - -if strcmpi(md.type,'2d'), - md=icetransient2d(md); -else - md=icetransient3d(md); -end Index: sm/trunk/src/m/solutions/ice/icetransient2d.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icetransient2d.m (revision 33) +++ (revision ) @@ -1,111 +1,0 @@ -function md=icetransient2d(md); -%ICETRANSIENT2D - transient solution of 2d models -% -% Usage: -% md=icetransient2d(md) - -global gridset - -%Create fem structure (input of icediagnostic2d) -fem=struct(); - -%Figure out which type of elements are present -[fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type); - -%First, build elements,grids,loads, etc ... for horizontal model -if fem.ishutter, - fem.m_ss=CreateFemModel(md,'surface_slope_compute'); - fem.m_dhu=CreateFemModel(md,'diagnostic_hutter'); -end -if fem.ismacayealpattyn, - fem.m_dh=CreateFemModel(md,'diagnostic_horiz'); -end -m_p=CreateFemModel(md,'prognostic'); - -%initialize (velocity,pressure) -solution.u_g=zeros(m_p.gridset.gsize,1); -if ~isnan(md.vx), solution.u_g(1:6:m_p.gridset.gsize)=md.vx/md.yts; end -if ~isnan(md.vy), solution.u_g(2:6:m_p.gridset.gsize)=md.vy/md.yts; end - -%initialize geometry and time -solution.thickness=zeros(m_p.gridset.gsize,1);solution.thickness(1:6:m_p.gridset.gsize)=md.thickness; -solution.surface=zeros(m_p.gridset.gsize,1);solution.surface(1:6:m_p.gridset.gsize)=md.surface; -solution.bed=zeros(m_p.gridset.gsize,1);solution.bed(1:6:m_p.gridset.gsize)=md.bed; -solution.time=0; - -%melting and accumulation are taken constant for all iterations -melting=zeros(m_p.gridset.gsize,1); -accumulation=zeros(m_p.gridset.gsize,1); -melting(1:6:end)=md.melting/md.yts; %from m/yr to m/s -accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s - -%build constant inputs -fem.inputs=struct(); -fem.inputs.accumulation=accumulation; -fem.inputs.melting=melting; -fem.inputs.dt=md.dt; - -%first time step is given by model. -dt=md.dt; -finaltime=md.ndt; -time=dt; -n=1; %counter - -while time