source: issm/trunk/test/NightlyRun/test233.m@ 8824

Last change on this file since 8824 was 8824, checked in by Eric.Larour, 14 years ago

Split spcvelocity into spcvx, spcvy and spcvz

File size: 2.0 KB
Line 
1md=mesh(model,'../Exp/Square.exp',150000);
2md=geography(md,'all','');
3md=parameterize(md,'../Par/SquareShelf.par');
4md=setelementstype(md,'macayeal','all');
5md.cluster=none;
6
7%redo the parameter file for this special shelf.
8%constant thickness, constrained (vy=0) flow into an icefront,
9%from 0 m/yr at the grounding line.
10
11%tighten
12md.eps_res=10^-4;
13
14%needed later
15ymin=min(md.y);
16ymax=max(md.y);
17xmin=min(md.x);
18xmax=max(md.x);
19
20di=md.rho_ice/md.rho_water;
21
22h=1000;
23md.thickness=h*ones(md.numberofnodes,1);
24md.bed=-md.rho_ice/md.rho_water*md.thickness;
25md.surface=md.bed+md.thickness;
26
27%Initial velocity and pressure
28md.vx=zeros(md.numberofnodes,1);
29md.vy=zeros(md.numberofnodes,1);
30md.vz=zeros(md.numberofnodes,1);
31md.pressure=zeros(md.numberofnodes,1);
32
33%Materials
34md.observed_temperature=(273-20)*ones(md.numberofnodes,1);
35md.rheology_B=paterson(md.observed_temperature);
36md.rheology_n=3*ones(md.numberofelements,1);
37md.temperature=md.observed_temperature;
38
39%Boundary conditions:
40md.spcvx=NaN*ones(md.numberofnodes,1);
41md.spcvy=NaN*ones(md.numberofnodes,1);
42md.spcvz=NaN*ones(md.numberofnodes,1);
43
44%constrain flanks to 0 normal velocity
45pos=find(md.x==xmin | md.x==xmax);
46md.spcvx(pos)=0;
47md.spcvz(pos)=NaN;
48
49%constrain grounding line to 0 velocity
50pos=find(md.y==ymin);
51md.spcvx(pos)=0;
52md.spcvy(pos)=0;
53
54%icefront
55nodeonicefront=zeros(md.numberofnodes,1);
56pos=find(md.y==ymax); nodeonicefront(pos)=1;
57pos=find(nodeonicefront(md.segments(:,1)) | nodeonicefront(md.segments(:,2))); pressureload=md.segments(pos,:);
58pressureload=[pressureload WaterEnum*md.elementoniceshelf(pressureload(:,end))];
59md.pressureload=pressureload;
60
61md=solve(md,DiagnosticSolutionEnum);
62
63%create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
64%ey_c=(md.rho_ice*md.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
65%vy_c=ey_c.*md.y*md.yts;
66
67%Fields and tolerances to track changes
68field_names ={'Vy'};
69field_tolerances={1e-13};
70field_values={PatchToVec(md.results.DiagnosticSolution.Vy)};
Note: See TracBrowser for help on using the repository browser.