1 | %Start defining model parameters here
|
---|
2 |
|
---|
3 | %Geometry
|
---|
4 | hmin=300;
|
---|
5 | hmax=1000;
|
---|
6 | ymin=min(md.mesh.y);
|
---|
7 | ymax=max(md.mesh.y);
|
---|
8 | xmin=min(md.mesh.x);
|
---|
9 | xmax=max(md.mesh.x);
|
---|
10 | md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin);
|
---|
11 | md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
|
---|
12 | md.geometry.surface=md.geometry.base+md.geometry.thickness;
|
---|
13 |
|
---|
14 | %Initial velocity and pressure
|
---|
15 | x = transpose(ncread('../Data/SquareShelf.nc','x'));
|
---|
16 | y = transpose(ncread('../Data/SquareShelf.nc','y'));
|
---|
17 | vx = transpose(ncread('../Data/SquareShelf.nc','vx'));
|
---|
18 | vy = transpose(ncread('../Data/SquareShelf.nc','vy'));
|
---|
19 | index = transpose(ncread('../Data/SquareShelf.nc','index'));
|
---|
20 | md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y);
|
---|
21 | md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y);
|
---|
22 | clear vx vy x y index;
|
---|
23 | md.initialization.vz=zeros(md.mesh.numberofvertices,1);
|
---|
24 | md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
|
---|
25 |
|
---|
26 | %Materials
|
---|
27 | md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
|
---|
28 | md.materials.rheology_B=paterson(md.initialization.temperature);
|
---|
29 | md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
|
---|
30 |
|
---|
31 | %Friction
|
---|
32 | md.friction.coefficient=20*ones(md.mesh.numberofvertices,1);
|
---|
33 | md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.;
|
---|
34 | md.friction.p=ones(md.mesh.numberofelements,1);
|
---|
35 | md.friction.q=ones(md.mesh.numberofelements,1);
|
---|
36 |
|
---|
37 | %Numerical parameters
|
---|
38 | md.masstransport.stabilization=1;
|
---|
39 | md.thermal.stabilization=1;
|
---|
40 | md.settings.waitonlock=30;
|
---|
41 | md.verbose=verbose(0);
|
---|
42 | md.stressbalance.restol=0.10;
|
---|
43 | md.steadystate.reltol=0.02;
|
---|
44 | md.stressbalance.reltol=0.02;
|
---|
45 | md.stressbalance.abstol=NaN;
|
---|
46 | md.timestepping.time_step=1;
|
---|
47 | md.timestepping.final_time=3;
|
---|
48 |
|
---|
49 | %Boundary conditions:
|
---|
50 | md=SetIceShelfBC(md,'./SquareFront.exp');
|
---|
51 |
|
---|
52 | %Change name so that no test have the same name
|
---|
53 | A=dbstack;
|
---|
54 | if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end
|
---|