Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/EISMINT.par =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/EISMINT.par (revision 14124) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/EISMINT.par (revision 14124) @@ -0,0 +1,57 @@ +disp(' creating thickness'); +hmin=0.01; +hmax=2756.7; +radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); +radiusmax=max(radius); +md.geometry.thickness=hmin*ones(size(md.mesh.x,1),1)+hmax*(4.*((1./2.)^(4./3.)*ones(size(md.mesh.x,1),1)-((radius)./(2.*radiusmax)).^(4./3.))).^(3./8.); +md.geometry.bed=0.*md.geometry.thickness; +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + +disp(' creating drag'); +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); %q=1. no drag is specified in the analytical solution +%Take care of iceshelves: no basal drag +pos=find(md.mask.elementonfloatingice); +md.friction.coefficient(md.mesh.elements(pos,:))=0.; +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +disp(' creating temperatures'); +tmin=238.15; %K +st=1.67*10^-2/1000.; %k/m +md.initialization.temperature=tmin+st*radius; +md.basalforcings.geothermalflux=4.2*10^-2*ones(md.mesh.numberofvertices,1); + +disp(' creating flow law parameter'); +md.materials.rheology_B=6.81*10^7*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); + +disp(' creating surface mass balance'); +smb_max=0.5; %m/yr +sb=10^-2/1000.; %m/yr/m +rel=450.*1000.; %m +md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius)); + +disp(' creating velocities'); +constant=0.3; +md.inversion.vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1; +md.inversion.vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1; +md.inversion.vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2); +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); + +%Deal with boundary conditions: +disp(' boundary conditions for diagnostic model:'); +md=SetMarineIceSheetBC(md,'RoundFront.exp'); + +radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); +pos=find(radius==min(radius)); +md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center + +md.diagnostic.spcvx(pos)=0.; +md.diagnostic.spcvy(pos)=0.; +md.diagnostic.spcvz(pos)=0.; + +%Convergence criterion +md.diagnostic.restol=0.01; Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/RoundFront.exp =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/RoundFront.exp (revision 14124) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/RoundFront.exp (revision 14124) @@ -0,0 +1,47 @@ +## Name: +## Icon:0 +# Points Count Value +41 1 +# X pos Y pos +1000000.000000 0.000000 +987688.340595 156434.465040 +951056.516295 309016.994375 +891006.524188 453990.499740 +809016.994375 587785.252292 +707106.781187 707106.781187 +587785.252292 809016.994375 +453990.499740 891006.524188 +309016.994375 951056.516295 +156434.465040 987688.340595 +0.000000 1000000.000000 +-156434.465040 987688.340595 +-309016.994375 951056.516295 +-453990.499740 891006.524188 +-587785.252292 809016.994375 +-707106.781187 707106.781187 +-809016.994375 587785.252292 +-891006.524188 453990.499740 +-951056.516295 309016.994375 +-987688.340595 156434.465040 +-1000000.000000 0.000000 +-987688.340595 -156434.465040 +-951056.516295 -309016.994375 +-891006.524188 -453990.499740 +-809016.994375 -587785.252292 +-707106.781187 -707106.781187 +-587785.252292 -809016.994375 +-453990.499740 -891006.524188 +-309016.994375 -951056.516295 +-156434.465040 -987688.340595 +-0.000000 -1000000.000000 +156434.465040 -987688.340595 +309016.994375 -951056.516295 +453990.499740 -891006.524188 +587785.252292 -809016.994375 +707106.781187 -707106.781187 +809016.994375 -587785.252292 +891006.524188 -453990.499740 +951056.516295 -309016.994375 +987688.340595 -156434.465040 +1000000.000000 0.000000 + Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/Square.exp =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/Square.exp (revision 14124) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/Square.exp (revision 14124) @@ -0,0 +1,10 @@ +## Name:domainoutline +## Icon:0 +# Points Count Value +5 1. +# X pos Y pos +0 0 +1000000 0 +1000000 1000000 +0 1000000 +0 0 Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SquareFront.exp =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SquareFront.exp (revision 14124) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SquareFront.exp (revision 14124) @@ -0,0 +1,10 @@ +## Name:icefront +## Icon:0 +# Points Count Value +5 1. +# X pos Y pos +-1000 900000 +-1000 1100000 +1100000 1100000 +1100000 900000 +-1000 900000 Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SquareShelf.par =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SquareShelf.par (revision 14124) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SquareShelf.par (revision 14124) @@ -0,0 +1,54 @@ +%Start defining model parameters here + +%Geometry +hmin=300; +hmax=1000; +ymin=min(md.mesh.y); +ymax=max(md.mesh.y); +md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); +md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + +%Initial velocity and pressure +x = transpose(ncread('SquareShelf.nc','x')); +y = transpose(ncread('SquareShelf.nc','y')); +vx = transpose(ncread('SquareShelf.nc','vx')); +vy = transpose(ncread('SquareShelf.nc','vy')); +index = transpose(ncread('SquareShelf.nc','index')); +md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y); +md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y); +clear vx vy x y index; +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); + +%Materials +md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); +md.materials.rheology_B=paterson(md.initialization.temperature); +md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + +%Friction +pos=find(md.mask.elementonfloatingice); +md.friction.coefficient=20*ones(md.mesh.numberofvertices,1); +md.friction.coefficient(md.mesh.elements(pos,:))=0; +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +%Numerical parameters +md.diagnostic.viscosity_overshoot=0.3; +md.prognostic.stabilization=1; +md.thermal.stabilization=1; +md.settings.waitonlock=30; +md.verbose=verbose(0); +md.diagnostic.restol=0.10; +md.steadystate.reltol=0.02; +md.diagnostic.reltol=0.02; +md.diagnostic.abstol=NaN; +md.timestepping.time_step=1; +md.timestepping.final_time=3; + +%Boundary conditions: +... + +%Change name so that no test have the same name +A=dbstack; +if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/eismint.m =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/eismint.m (revision 14124) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/eismint.m (revision 14124) @@ -0,0 +1,45 @@ +md=model(); +%Create mesh with roundmesh +... +%Set mask +... +%Parameterize model +... + +%We extrude the model to have a 3d model +... + +%Set ice flow approximation +... + +%Create boundary conditions: zero velocity on the bed +... + +%Go Solve +md.cluster=generic('np',2); +md.verbose.convergence=1; +md=solve(md,DiagnosticSolutionEnum()); +vel=DepthAverage(md,sqrt(md.results.DiagnosticSolution.Vx.^2+md.results.DiagnosticSolution.Vy.^2)); + +%Calculate analytical velocity +constant=0.3; +vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; +vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; +vel_obs=sqrt(vx_obs.^2+vy_obs.^2); +vel_obs=project2d(md,vel_obs,1); + +plotmodel(md,... + 'data',vel ,'view',2,'caxis',[0 200],'title','Modelled velocity',... + 'data',vel_obs,'view',2,'caxis',[0 200],'title','Analytical velocity',... + 'data',abs(vel-vel_obs)./(vel_obs+eps)*100,'caxis',[0 30],'view',2,'title','Relative misfit (%)'); + +subplot(2,2,4) +hold on; +plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel,'r.'); +plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.'); +title('Analytical vs calculated velocity'); +xlabel('distance to the center of the icesheet [m]'); +ylabel('velocity (m/yr)'); +legend('calculated velocity','exact velocity'); +axis([0 750000 0 200]); +hold off; Index: /issm/workshop/2012/Talks/07_Ice_flow_models/Code/runme.m =================================================================== --- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/runme.m (revision 14123) +++ /issm/workshop/2012/Talks/07_Ice_flow_models/Code/runme.m (revision 14124) @@ -1,69 +1,7 @@ ->> 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') +md=triangle(model(),'Square.exp',50000.); +md=setmask(md,'all',''); +md=parameterize(md,'SquareShelf.par'); +md=extrude(md,3,1); +%Set flow equation +... +%Solve