1 | %Start defining model parameters here
|
---|
2 |
|
---|
3 | %Geometry and observation
|
---|
4 | x = archread('../Data/Pig.arch','x');
|
---|
5 | y = archread('../Data/Pig.arch','y');
|
---|
6 | vx_obs = archread('../Data/Pig.arch','vx_obs');
|
---|
7 | vy_obs = archread('../Data/Pig.arch','vy_obs');
|
---|
8 | index = archread('../Data/Pig.arch','index');
|
---|
9 | surface = archread('../Data/Pig.arch','surface');
|
---|
10 | thickness = archread('../Data/Pig.arch','thickness');
|
---|
11 | bed = archread('../Data/Pig.arch','bed');
|
---|
12 |
|
---|
13 | x = x{1};
|
---|
14 | y = y{1};
|
---|
15 | vx_obs = vx_obs{1};
|
---|
16 | vy_obs = vy_obs{1};
|
---|
17 | index = index{1};
|
---|
18 | surface = surface{1};
|
---|
19 | thickness = thickness{1};
|
---|
20 | bed = bed{1};
|
---|
21 |
|
---|
22 | md.inversion.vx_obs =InterpFromMeshToMesh2d(index,x,y,vx_obs,md.mesh.x,md.mesh.y);
|
---|
23 | md.inversion.vy_obs =InterpFromMeshToMesh2d(index,x,y,vy_obs,md.mesh.x,md.mesh.y);
|
---|
24 | md.geometry.surface =InterpFromMeshToMesh2d(index,x,y,surface,md.mesh.x,md.mesh.y);
|
---|
25 | md.geometry.thickness=InterpFromMeshToMesh2d(index,x,y,thickness,md.mesh.x,md.mesh.y);
|
---|
26 | md.geometry.base=md.geometry.surface-md.geometry.thickness;
|
---|
27 | md.geometry.bed =md.geometry.base;
|
---|
28 | pos = find(md.mask.ocean_levelset<0);
|
---|
29 | md.geometry.bed(pos) =InterpFromMeshToMesh2d(index,x,y,bed,md.mesh.x(pos),md.mesh.y(pos));
|
---|
30 | clear surface thickness vx_obs vy_obs x y index;
|
---|
31 | md.initialization.vx=md.inversion.vx_obs;
|
---|
32 | md.initialization.vy=md.inversion.vy_obs;
|
---|
33 | md.initialization.vz=zeros(md.mesh.numberofvertices,1);
|
---|
34 | md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
|
---|
35 |
|
---|
36 | %Materials
|
---|
37 | md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
|
---|
38 | md.materials.rheology_B=paterson(md.initialization.temperature);
|
---|
39 | md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
|
---|
40 | md.initialization.temperature=md.initialization.temperature;
|
---|
41 |
|
---|
42 | %Friction
|
---|
43 | md.friction.coefficient=50*ones(md.mesh.numberofvertices,1);
|
---|
44 | md.friction.coefficient(find(md.mask.ocean_levelset<0.))=0.;
|
---|
45 | md.friction.p=ones(md.mesh.numberofelements,1);
|
---|
46 | md.friction.q=ones(md.mesh.numberofelements,1);
|
---|
47 |
|
---|
48 | %Numerical parameters
|
---|
49 | md.masstransport.stabilization=1;
|
---|
50 | md.verbose=verbose(0);
|
---|
51 | md.settings.waitonlock=30;
|
---|
52 | md.timestepping.time_step=1;
|
---|
53 | md.timestepping.final_time=2;
|
---|
54 | md.stressbalance.restol=0.05;
|
---|
55 | md.stressbalance.reltol=1;
|
---|
56 | md.steadystate.reltol=1;
|
---|
57 | md.stressbalance.abstol=NaN;
|
---|
58 | md.groundingline.migration='None';
|
---|
59 |
|
---|
60 | %Boundary conditions:
|
---|
61 | md=SetMarineIceSheetBC(md);
|
---|
62 |
|
---|
63 | %Change name so that no test have the same name
|
---|
64 | A=dbstack;
|
---|
65 | if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end
|
---|