Index: sm/workshop/2012/Talks/06_Mesh_generation/Code/circles.m =================================================================== --- /issm/workshop/2012/Talks/06_Mesh_generation/Code/circles.m (revision 14071) +++ (revision ) @@ -1,7 +1,0 @@ -function vel=circles(x,y), - - u=4*x-2; v=4*y-2; - - vel=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; Index: sm/workshop/2012/Talks/06_Mesh_generation/Code/ex1.m =================================================================== --- /issm/workshop/2012/Talks/06_Mesh_generation/Code/ex1.m (revision 14071) +++ (revision ) @@ -1,11 +1,0 @@ -md=triangle(model,'Square.exp',.2); -plotmodel(md,'data','mesh','axis','off equal','title','');export_fig('../Images/Mesh1.pdf'); -md=triangle(model,'Square.exp',.05); -plotmodel(md,'data','mesh','axis','off equal','title','');export_fig('../Images/Mesh2.pdf'); -md=bamg(model,'domain','Square.exp','hmax',.05); -plotmodel(md,'data','mesh','axis','off equal','title','');export_fig('../Images/Mesh_bamg1.pdf'); -md=squaremesh(model,100,200,15,25); -plotmodel(md,'data','mesh','axis','on equal','title','','xlim',[0 100],'ylim',[0 120]);export_fig('../Images/Mesh_sq1.pdf'); -hvertices=[0.2;0.2;0.005;0.2]; -md=bamg(model,'domain','Square.exp','hVertices',hvertices); -plotmodel(md,'data','mesh','axis','on equal','title','','xlim',[0 1],'ylim',[0 1]);export_fig('../Images/Mesh_bamg2.pdf'); Index: sm/workshop/2012/Talks/06_Mesh_generation/Code/ex2.m =================================================================== --- /issm/workshop/2012/Talks/06_Mesh_generation/Code/ex2.m (revision 14071) +++ (revision ) @@ -1,50 +1,0 @@ -%Initial mesh to display field -steps=6; -if any(steps==1), -md=bamg(model,'domain','Square.exp','hmax',.002); -vel=shock(md.mesh.x,md.mesh.y); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0); export_fig('../Images/Bamgshock0.pdf'); -end -if any(steps==2), -md=bamg(model,'domain','Square.exp','hmax',.05); -vel=shock(md.mesh.x,md.mesh.y); -plotmodel(md,'data',vel,'edgecolor','w'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0); export_fig('../Images/Bamgshock1.pdf'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0,'edgecolor',[.99 .99 .99]); export_fig('../Images/Bamgshock1b.pdf'); -end -if any(steps==3), -md=bamg(model,'domain','Square.exp','hmax',.005); -vel=shock(md.mesh.x,md.mesh.y); -md=bamg(md,'field',vel,'err',0.05,'hmin',0.005,'hmax',0.3); -vel=shock(md.mesh.x,md.mesh.y); -plotmodel(md,'data',vel,'edgecolor','w'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0); export_fig('../Images/Bamgshock2.pdf'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0,'edgecolor',[.99 .99 .99]); export_fig('../Images/Bamgshock2b.pdf'); -end -if any(steps==4), -md=bamg(model,'domain','Square.exp','hmax',.005); -vel=shock(md.mesh.x,md.mesh.y); -md=bamg(md,'field',vel,'err',0.03,'hmin',0.005,'hmax',0.3); -vel=shock(md.mesh.x,md.mesh.y); -plotmodel(md,'data',vel,'edgecolor','w'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0); export_fig('../Images/Bamgshock3.pdf'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0,'edgecolor',[.99 .99 .99]); export_fig('../Images/Bamgshock3b.pdf'); -end -if any(steps==5), -md=bamg(model,'domain','Square.exp','hmax',.005); -vel=shock(md.mesh.x,md.mesh.y); -md=bamg(md,'field',vel,'err',0.03,'hmin',0.005,'hmax',0.3,'gradation',3); -vel=shock(md.mesh.x,md.mesh.y); -plotmodel(md,'data',vel,'edgecolor','w'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0); export_fig('../Images/Bamgshock4.pdf'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0,'edgecolor',[.99 .99 .99]); export_fig('../Images/Bamgshock4b.pdf'); -end -if any(steps==6), -md=bamg(model,'domain','Square.exp','hmax',.005); -vel=shock(md.mesh.x,md.mesh.y); -md=bamg(md,'field',vel,'err',0.03,'hmin',0.005,'hmax',0.3,'gradation',1.3,'anisomax',1); -vel=shock(md.mesh.x,md.mesh.y); -plotmodel(md,'data',vel,'edgecolor','w'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0); export_fig('../Images/Bamgshock5.pdf'); -plotmodel(md,'data',vel,'xlim',[0 1],'ylim',[0 1],'gridded',1,'whitepos',0,'edgecolor',[.99 .99 .99]); export_fig('../Images/Bamgshock5b.pdf'); -end Index: sm/workshop/2012/Talks/06_Mesh_generation/Code/figs.m =================================================================== --- /issm/workshop/2012/Talks/06_Mesh_generation/Code/figs.m (revision 14071) +++ (revision ) @@ -1,48 +1,0 @@ -steps=[2]; -if any(steps==1), - md=squaremesh(model,1,1,2,2); - vel=zeros(md.mesh.numberofvertices,1); - plotmodel(md,'data',vel,'title','Constant field','colorbar',0,'xlim',[0 1],'ylim',[0 1]);export_fig('../Images/Interp1.pdf'); - plotmodel(md,'data',vel,'title','Constant field','colorbar',0,'xlim',[0 1],'ylim',[0 1],'edgecolor','k');export_fig('../Images/Interp2.pdf'); -end -if any(steps==2), - md=squaremesh(model,1,1,2,2); - vel=md.mesh.x+.5*md.mesh.y; - md.mesh.z=vel; - plotmodel(md,'data',vel,'title','Linear field','colorbar',0,'xlim',[0 1],'ylim',[0 1],'view',[-30 40]);export_fig('../Images/Interp3.pdf'); - plotmodel(md,'data',vel,'title','Linear field','colorbar',0,'xlim',[0 1],'ylim',[0 1],'view',[-30 40],'edgecolor','k');export_fig('../Images/Interp4.pdf'); -end -if any(steps==3), - md=squaremesh(model,1,1,40,2); - vel=-5*abs(md.mesh.x-.5).^3 +.5^3; md.mesh.z=vel; - plotmodel(md,'data',vel,'title','Non linear field','colorbar',0,'xlim',[0 1],'ylim',[0 1],'view',[40 10]);export_fig('../Images/Interp5.pdf'); -end -if any(steps==4), - md=squaremesh(model,1,1,5,5); - vel=-5*abs(md.mesh.x-.5).^3 +.5^3; md.mesh.z=vel; - plotmodel(md,'data',vel,'title','32 elements','colorbar',0,'edgecolor','k','xlim',[0 1],'ylim',[0 1],'view',[40 10]);export_fig('../Images/Interp5b.pdf'); -end -if any(steps==5), - %1500 elements - md=squaremesh(model,1,1,30,30); - vel=-5*abs(md.mesh.x-.5).^3 +.5^3; md.mesh.z=vel; - plotmodel(md,'data',vel,'title','1,500 elements','colorbar',0,'xlim',[0 1],'ylim',[0 1],'edgecolor','k','view',[40 10]);export_fig('../Images/Interp6.pdf'); -end -if any(steps==6), - %1000 elements - md=squaremesh(model,1,1,20,30); - md.mesh.x=sin((md.mesh.x-0.5)*pi); - md.mesh.x=md.mesh.x-min(md.mesh.x); - md.mesh.x=md.mesh.x/(max(md.mesh.x)); - vel=-5*abs(md.mesh.x-.5).^3 +.5^3; md.mesh.z=vel; - plotmodel(md,'data',vel,'title','1,000 elements','colorbar',0,'xlim',[0 1],'ylim',[0 1],'edgecolor','k','view',[40 10]);export_fig('../Images/Interp7.pdf'); -end -if any(steps==7), - %40 elements - md=squaremesh(model,1,1,20,2); - md.mesh.x=sin((md.mesh.x-0.5)*pi); - md.mesh.x=md.mesh.x-min(md.mesh.x); - md.mesh.x=md.mesh.x/(max(md.mesh.x)); - vel=-5*abs(md.mesh.x-.5).^3 +.5^3; md.mesh.z=vel; - plotmodel(md,'data',vel,'title','40 elements','colorbar',0,'xlim',[0 1],'ylim',[0 1],'edgecolor','k','view',[40 10]);export_fig('../Images/Interp8.pdf'); -end Index: sm/workshop/2012/Talks/07_Ice_flow_models/Code/ex1.m =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/ex1.m (revision 14071) +++ (revision ) @@ -1,69 +1,0 @@ ->> md.flowequation - -ans = - - flow equation parameters: - ismacayealpattyn : 0 -- is the macayeal or pattyn approximation used ? - ishutter : 0 -- is the shallow ice approximation used ? - isstokes : 0 -- are the Full-Stokes equations used ? - vertex_equation : N/A -- flow equation for each vertex - element_equation : N/A -- flow equation for each element - bordermacayeal : N/A -- vertices on MacAyeal's border (for tiling) - borderpattyn : N/A -- vertices on Pattyn's border (for tiling) - borderstokes : N/A -- vertices on Stokes' border (for tiling) - -md=setflowequation(md,'hutter','all') -md=setflowequation(md,'stokes','all') -md=setflowequation(md,'macayeal','all') -md=setflowequation(md,'pattyn','all') - ->> plotmodel(md,'data','elements_type') - ->> md.diagnostic - -ans = - - Diagnostic solution parameters: - - Convergence criteria: - restol : 0.0001 -- mechanical equilibrium residue convergence criterion - reltol : 0.01 -- velocity relative convergence criterion, NaN -> not applied - abstol : 10 -- velocity absolute convergence criterion, NaN -> not applied - maxiter : 100 -- maximum number of nonlinear iterations - viscosity_overshoot : 0 -- over-shooting constant new=new+C*(new-old) - - boundary conditions: - spcvx : N/A -- x-axis velocity constraint (NaN means no constraint) - spcvy : N/A -- y-axis velocity constraint (NaN means no constraint) - spcvz : N/A -- z-axis velocity constraint (NaN means no constraint) - icefront : N/A -- segments on ice front list (last column 0-> Air, 1-> Water, 2->Ice - - Rift options: - rift_penalty_threshold : 0 -- threshold for instability of mechanical constraints - rift_penalty_lock : 10 -- number of iterations before rift penalties are locked - - Penalty options: - penalty_factor : 3 -- offset used by penalties: penalty = Kmax*10^offset - vertex_pairing : N/A -- pairs of vertices that are penalized - - Other: - shelf_dampening : 0 -- use dampening for floating ice ? Only for Stokes model - stokesreconditioning : 10000000000000 -- multiplier for incompressibility equation. Only for Stokes model - referential : N/A -- local referential - requested_outputs : N/A -- additional outputs requested - ->> md=SetIceSheetBC(md) ->> md=SetIceShelfBC(md,'Front.exp') ->> md=SetMarineIceSheefBC(md,'Front.exp') - ->> plotmodel(md,'data','BC') - -md=setflowequation(md,'pattyn',md.elementongroundedice,'fill','macayeal','coupling','penalties') -md=setflowequation(md,'pattyn',md.elementongroundedice,'fill','macayeal','coupling','tiling') -md=setflowequation(md,'stokes','Contour.exp','fill','pattyn') - ->> exptool('Contour.exp') - ->> plotmodel(md,'data','elements_type','edgecolor','k','expdisp','Contour.exp') - -md=setflowequation(md,'macayeal','FloatingIce.exp','fill','pattyn','coupling','penalties') Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/runme.m =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/runme.m (revision 14072) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/runme.m (revision 14072) @@ -0,0 +1,69 @@ +>> md.flowequation + +ans = + + flow equation parameters: + ismacayealpattyn : 0 -- is the macayeal or pattyn approximation used ? + ishutter : 0 -- is the shallow ice approximation used ? + isstokes : 0 -- are the Full-Stokes equations used ? + vertex_equation : N/A -- flow equation for each vertex + element_equation : N/A -- flow equation for each element + bordermacayeal : N/A -- vertices on MacAyeal's border (for tiling) + borderpattyn : N/A -- vertices on Pattyn's border (for tiling) + borderstokes : N/A -- vertices on Stokes' border (for tiling) + +md=setflowequation(md,'hutter','all') +md=setflowequation(md,'stokes','all') +md=setflowequation(md,'macayeal','all') +md=setflowequation(md,'pattyn','all') + +>> plotmodel(md,'data','elements_type') + +>> md.diagnostic + +ans = + + Diagnostic solution parameters: + + Convergence criteria: + restol : 0.0001 -- mechanical equilibrium residue convergence criterion + reltol : 0.01 -- velocity relative convergence criterion, NaN -> not applied + abstol : 10 -- velocity absolute convergence criterion, NaN -> not applied + maxiter : 100 -- maximum number of nonlinear iterations + viscosity_overshoot : 0 -- over-shooting constant new=new+C*(new-old) + + boundary conditions: + spcvx : N/A -- x-axis velocity constraint (NaN means no constraint) + spcvy : N/A -- y-axis velocity constraint (NaN means no constraint) + spcvz : N/A -- z-axis velocity constraint (NaN means no constraint) + icefront : N/A -- segments on ice front list (last column 0-> Air, 1-> Water, 2->Ice + + Rift options: + rift_penalty_threshold : 0 -- threshold for instability of mechanical constraints + rift_penalty_lock : 10 -- number of iterations before rift penalties are locked + + Penalty options: + penalty_factor : 3 -- offset used by penalties: penalty = Kmax*10^offset + vertex_pairing : N/A -- pairs of vertices that are penalized + + Other: + shelf_dampening : 0 -- use dampening for floating ice ? Only for Stokes model + stokesreconditioning : 10000000000000 -- multiplier for incompressibility equation. Only for Stokes model + referential : N/A -- local referential + requested_outputs : N/A -- additional outputs requested + +>> md=SetIceSheetBC(md) +>> md=SetIceShelfBC(md,'Front.exp') +>> md=SetMarineIceSheefBC(md,'Front.exp') + +>> plotmodel(md,'data','BC') + +md=setflowequation(md,'pattyn',md.elementongroundedice,'fill','macayeal','coupling','penalties') +md=setflowequation(md,'pattyn',md.elementongroundedice,'fill','macayeal','coupling','tiling') +md=setflowequation(md,'stokes','Contour.exp','fill','pattyn') + +>> exptool('Contour.exp') + +>> plotmodel(md,'data','elements_type','edgecolor','k','expdisp','Contour.exp') + +md=setflowequation(md,'macayeal','FloatingIce.exp','fill','pattyn','coupling','penalties')