[5337] | 1 | md=mesh(model,'../Exp/Square.exp',150000);
|
---|
| 2 | md=geography(md,'all','');
|
---|
| 3 | md=parameterize(md,'../Par/SquareShelf.par');
|
---|
| 4 | md=setelementstype(md,'macayeal','all');
|
---|
[5955] | 5 | md.cluster=none;
|
---|
[5337] | 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
|
---|
| 12 | md.eps_res=10^-4;
|
---|
| 13 |
|
---|
| 14 | %needed later
|
---|
| 15 | ymin=min(md.y);
|
---|
| 16 | ymax=max(md.y);
|
---|
| 17 | xmin=min(md.x);
|
---|
| 18 | xmax=max(md.x);
|
---|
| 19 |
|
---|
| 20 | di=md.rho_ice/md.rho_water;
|
---|
| 21 |
|
---|
| 22 | h=1000;
|
---|
[8304] | 23 | md.thickness=h*ones(md.numberofnodes,1);
|
---|
[5337] | 24 | md.bed=-md.rho_ice/md.rho_water*md.thickness;
|
---|
| 25 | md.surface=md.bed+md.thickness;
|
---|
| 26 |
|
---|
| 27 | %Initial velocity and pressure
|
---|
[8304] | 28 | md.vx=zeros(md.numberofnodes,1);
|
---|
| 29 | md.vy=zeros(md.numberofnodes,1);
|
---|
| 30 | md.vz=zeros(md.numberofnodes,1);
|
---|
| 31 | md.pressure=zeros(md.numberofnodes,1);
|
---|
[5337] | 32 |
|
---|
| 33 | %Materials
|
---|
[8304] | 34 | md.observed_temperature=(273-20)*ones(md.numberofnodes,1);
|
---|
[5337] | 35 | md.rheology_B=paterson(md.observed_temperature);
|
---|
| 36 | md.rheology_n=3*ones(md.numberofelements,1);
|
---|
| 37 | md.temperature=md.observed_temperature;
|
---|
| 38 |
|
---|
| 39 | %Boundary conditions:
|
---|
[8304] | 40 | md.spcvelocity=zeros(md.numberofnodes,6);
|
---|
[5337] | 41 |
|
---|
| 42 | %constrain flanks to 0 normal velocity
|
---|
| 43 | pos=find(md.x==xmin | md.x==xmax);
|
---|
| 44 | md.spcvelocity(pos,1)=1;
|
---|
| 45 | md.spcvelocity(pos,3)=0;
|
---|
| 46 |
|
---|
| 47 | %constrain grounding line to 0 velocity
|
---|
| 48 | pos=find(md.y==ymin);
|
---|
| 49 | md.spcvelocity(pos,1:2)=1;
|
---|
| 50 | md.spcvelocity(pos,3:4)=0;
|
---|
| 51 |
|
---|
| 52 | %icefront
|
---|
[8304] | 53 | nodeonicefront=zeros(md.numberofnodes,1);
|
---|
| 54 | pos=find(md.y==ymax); nodeonicefront(pos)=1;
|
---|
| 55 | pos=find(nodeonicefront(md.segments(:,1)) | nodeonicefront(md.segments(:,2))); pressureload=md.segments(pos,:);
|
---|
[5337] | 56 | pressureload=[pressureload WaterEnum*md.elementoniceshelf(pressureload(:,end))];
|
---|
| 57 | md.pressureload=pressureload;
|
---|
| 58 |
|
---|
[8295] | 59 | md=solve(md,DiagnosticSolutionEnum);
|
---|
[5337] | 60 |
|
---|
| 61 | %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
|
---|
| 62 | %ey_c=(md.rho_ice*md.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
|
---|
| 63 | %vy_c=ey_c.*md.y*md.yts;
|
---|
| 64 |
|
---|
| 65 | %Fields and tolerances to track changes
|
---|
| 66 | field_names ={'Vy'};
|
---|
| 67 | field_tolerances={1e-13};
|
---|
| 68 | field_values={PatchToVec(md.results.DiagnosticSolution.Vy)};
|
---|