Skip to content

Commit fa33a8e

Browse files
aleahsommersaleahsommers
committed
ADD:adding Helheim SHAKTI tutorials
1 parent a43ae18 commit fa33a8e

File tree

4 files changed

+272
-0
lines changed

4 files changed

+272
-0
lines changed

‎examples/Helheim/Greenland.par‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
%required Data
2+
accpath = ['/home/ModelData/Greenland/RACMO2Accumulation/SMBGreenland/smb_1961-1990.mat'];
3+
velpath = ['/home/ModelData/Greenland/VelMouginot/RignotGreenland2012Vel.mat'];
4+
5+
%WARNING: we do have a levelset here, so we must set it to one even for stressbalances...
6+
md.transient.ismovingfront=1;
7+
8+
disp(' Interpolating mask');
9+
mask = int8(interpBedmachineGreenland(md.mesh.x,md.mesh.y,'mask'));
10+
md.mask.ice_levelset= -1*ones(md.mesh.numberofvertices,1);
11+
pos = find(mask<1);
12+
md.mask.ice_levelset(pos)=1;
13+
14+
disp(' reading MC bed (assumes no floating ice)');
15+
md.geometry.bed = interpBedmachineGreenland(md.mesh.x,md.mesh.y,'bed');
16+
md.geometry.base = md.geometry.bed;
17+
18+
disp(' reading Howat surface');
19+
md.geometry.surface=interpBedmachineGreenland(md.mesh.x,md.mesh.y,'surface');
20+
pos = find(md.mask.ice_levelset>0);
21+
md.geometry.surface(pos) = md.geometry.base(pos)+10; %Minimum thickness
22+
23+
md.geometry.thickness = md.geometry.surface - md.geometry.bed;
24+
pos=find(md.geometry.thickness<=10);
25+
md.geometry.surface(pos) = md.geometry.base(pos)+10; %Minimum thickness
26+
md.geometry.thickness = md.geometry.surface - md.geometry.bed;
27+
28+
md.masstransport.min_thickness = 10;
29+
30+
disp(' Adjusting ice mask');
31+
%Tricky part here: we want to offset the mask by one element so that we don't end up with a cliff at the transition
32+
pos = find(max(md.mask.ice_levelset(md.mesh.elements),[],2)>0);
33+
md.mask.ice_levelset(md.mesh.elements(pos,:)) = 1;
34+
% For the region where surface is NaN, set thickness to small value (consistency requires >0)
35+
pos=find((md.mask.ice_levelset<0).*(md.geometry.surface<0));
36+
md.mask.ice_levelset(pos)=1;
37+
pos=find((md.mask.ice_levelset<0).*(isnan(md.geometry.surface)));
38+
md.mask.ice_levelset(pos)=1;
39+
40+
disp(' -- reconstruct thickness');
41+
md.geometry.thickness=md.geometry.surface-md.geometry.base;
42+
43+
disp(' reading velocities ');
44+
[md.inversion.vx_obs md.inversion.vy_obs]=interpJoughinCompositeGreenland(md.mesh.x,md.mesh.y);
45+
pos=find(isnan(md.inversion.vx_obs) | isnan(md.inversion.vy_obs));
46+
md.inversion.vx_obs(pos)=0;
47+
md.inversion.vy_obs(pos)=0;
48+
md.inversion.vel_obs = sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
49+
md.initialization.vx = md.inversion.vx_obs;
50+
md.initialization.vy = md.inversion.vy_obs;
51+
md.initialization.vz = zeros(md.mesh.numberofvertices,1);
52+
md.initialization.vel = md.inversion.vel_obs;
53+
54+
%drag md.drag or stress
55+
md.friction.coefficient = 20*ones(md.mesh.numberofvertices,1); %q = 1.
56+
md.friction.p = ones(md.mesh.numberofelements,1);
57+
md.friction.q = ones(md.mesh.numberofelements,1);
58+
59+
%No friction on PURELY ocean element
60+
pos_e = find(min(md.mask.ice_levelset(md.mesh.elements),[],2)<0);
61+
flags=ones(md.mesh.numberofvertices,1);
62+
flags(md.mesh.elements(pos_e,:))=0;
63+
md.friction.coefficient(find(flags))=0;
64+
65+
%flow law
66+
disp(' Creating flow law parameters (assume ice is at 0°C for now)');
67+
md.materials.rheology_n = 3*ones(md.mesh.numberofelements,1);
68+
md.materials.rheology_B = cuffey(273.15 - 10)*ones(md.mesh.numberofvertices,1);
69+
70+
md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
71+
md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
72+
73+
disp(' Geothermal flux from Shapiro et al.');
74+
md.basalforcings.geothermalflux=interpSeaRISE(md.mesh.x,md.mesh.y,'bheatflx');
75+
76+
disp(' Setting up thermal model');
77+
md.initialization.temperature=min(0,interpSeaRISE(md.mesh.x,md.mesh.y,'surftemp'))+273.15;
78+
md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);
79+
md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
80+
md.thermal.spctemperature=md.initialization.temperature;
81+
md.thermal.isenthalpy=1;
82+
md.thermal.isdynamicbasalspc=1;
83+
84+
%Deal with boundary conditions:
85+
disp(' Set Boundary conditions');
86+
md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
87+
md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
88+
md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
89+
md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
90+
md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
91+
pos=find((md.mask.ice_levelset<0).*(md.mesh.vertexonboundary));
92+
md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
93+
md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
94+
md.stressbalance.spcvz(pos)=0;

‎examples/Helheim/Helheim_tut.exp‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
## Name:./Exp/Domain_TEMP.exp
2+
## Icon:0
3+
# Points Count Value
4+
14 1.000000
5+
# X pos Y pos
6+
251808.9599315871 -2519299.1448521991
7+
232017.0235357451 -2526258.4544728124
8+
232815.0354161858 -2578466.8389138547
9+
256982.9213724715 -2601667.2788287275
10+
302263.2150846492 -2594953.0557058132
11+
309967.1091796352 -2586747.6107035289
12+
315046.5997917140 -2581507.9891960463
13+
313861.3853155622 -2571621.9108800408
14+
310644.3745945790 -2567964.0619031191
15+
309120.5274109554 -2557188.2365386733
16+
306496.1239280481 -2552245.1973806708
17+
305818.8585131043 -2537910.3838224630
18+
285125.3861979210 -2517350.5381584275
19+
251808.9599315871 -2519299.1448521991
20+

‎examples/Helheim/runme.m‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
% runme script to perform a stress balance inversion for basal friction
2+
% on Helheim Glacier (SE Greenland)
3+
4+
% Step 1: Set up model domain outline and mesh
5+
% Step 2: Set parameters
6+
% Step 3: Inversion for basal friction
7+
8+
% Note that Steps 1 and 2 require previously downloading data for velocities,
9+
% surface elevation and bed elevation.
10+
11+
steps = [1:3]; % Choose which steps to run here
12+
13+
cluster=generic('name',oshostname(),'np',8);
14+
15+
org=organizer('repository',['./Models'],'prefix',['Model_Helheim_'],'steps',steps); clear steps;
16+
17+
if perform(org,'Mesh'),% {{{
18+
19+
% Initialize uniform mesh
20+
md=triangle(model,['Helheim_tut.exp'],500);
21+
22+
[velx vely]=interpJoughinCompositeGreenland(md.mesh.x,md.mesh.y);
23+
vel = sqrt(velx.^2+vely.^2);
24+
25+
% Refine mesh based on surface velocities
26+
md=bamg(md,'hmin',100,'hmax',1500,'field',vel,'err',5);
27+
[md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,45,70);
28+
md.mesh.epsg=3413;
29+
30+
savemodel(org,md);
31+
end %}}}
32+
if perform(org,'Param'),% {{{
33+
34+
md=loadmodel(org,'Mesh');
35+
md=setflowequation(md,'SSA','all'); % Set flow law to SSA
36+
37+
md=setmask(md,'','');
38+
md=parameterize(md,'Greenland.par');
39+
md.miscellaneous.name = 'Helheim';
40+
41+
savemodel(org,md);
42+
end%}}}
43+
if perform(org,'Inversion_drag'),% {{{
44+
45+
md=loadmodel(org,'Param');
46+
47+
%Control general
48+
md.inversion=m1qn3inversion(md.inversion);
49+
md.inversion.iscontrol=1;
50+
md.verbose=verbose('solution',false,'control',true);
51+
md.transient.amr_frequency = 0;
52+
53+
%Cost functions
54+
md.inversion.cost_functions=[101 103 501];
55+
md.inversion.cost_functions_coefficients=zeros(md.mesh.numberofvertices,numel(md.inversion.cost_functions));
56+
md.inversion.cost_functions_coefficients(:,1)=5000;
57+
md.inversion.cost_functions_coefficients(:,2)=10;
58+
md.inversion.cost_functions_coefficients(:,3)=2*50^-3;
59+
pos=find(md.mask.ice_levelset>0);
60+
md.inversion.cost_functions_coefficients(pos,1:2)=0;
61+
62+
%Controls
63+
md.inversion.control_parameters={'FrictionCoefficient'};
64+
md.inversion.maxsteps=100;
65+
md.inversion.maxiter =100;
66+
md.inversion.min_parameters=0.05*ones(md.mesh.numberofvertices,1);
67+
md.inversion.max_parameters=4000*ones(md.mesh.numberofvertices,1);
68+
md.inversion.control_scaling_factors=1;
69+
70+
%Additional parameters
71+
md.stressbalance.restol=0.01;
72+
md.stressbalance.reltol=0.1;
73+
md.stressbalance.abstol=NaN;
74+
75+
%Go solve
76+
md.cluster=cluster;
77+
78+
md=solve(md,'sb');
79+
80+
%Put results back into the model
81+
md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
82+
md.initialization.vx=md.results.StressbalanceSolution.Vx;
83+
md.initialization.vy=md.results.StressbalanceSolution.Vy;
84+
85+
savemodel(org,md);
86+
end%}}}

‎examples/HelheimSHAKTI/runme.m‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
% Coupled SHAKTI-ISSM simulation of Helheim, starting from previous inversion
2+
clear all;close all
3+
4+
% Load Helheim Glacier model
5+
load Models/Model_Helheim_Inversion_drag.mat
6+
7+
% Turn off inversion
8+
md.inversion.iscontrol=0;
9+
10+
% HYDROLOGY SPECIFIC PARAMETERIZATION:
11+
% Change hydrology class to SHAKTI model
12+
md.hydrology=hydrologyshakti();
13+
14+
% Define distributed englacial input to the subglacial system (m/yr)
15+
md.hydrology.englacial_input = .0*ones(md.mesh.numberofvertices,1);
16+
17+
% Define initial water head such that water pressure is 80% of ice overburden pressure
18+
md.hydrology.head = 0.8*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base;
19+
20+
% Initial subglacial gap height of 0.01m everywhere
21+
md.hydrology.gap_height = 0.01*ones(md.mesh.numberofelements,1);
22+
23+
% Typical bed bump bump spacing -- must be non-zero
24+
md.hydrology.bump_spacing = 1.0*ones(md.mesh.numberofelements,1);
25+
26+
% Typical bed bump height -- set to zero to turn off opening by sliding
27+
md.hydrology.bump_height = 0.0*ones(md.mesh.numberofelements,1);
28+
29+
% Initial Reynolds number (start at Re=1000 everywhere)
30+
md.hydrology.reynolds= 1000*ones(md.mesh.numberofelements,1);
31+
32+
% Boundary conditions
33+
md.hydrology.spchead = NaN(md.mesh.numberofvertices,1);
34+
35+
% Set head=0 for thin ice, everywhere <=10m ice thickness
36+
pos=find(md.geometry.thickness==10);
37+
md.hydrology.spchead(pos)=0;
38+
39+
% Set pressure in fjord equal to hydrostatic pressure of fjord water
40+
pos=find(md.mask.ice_levelset>0);
41+
md.hydrology.spchead(pos)=0;
42+
43+
% % Deal with terminus b.c. in SHAKTI
44+
% terminus=ContourToNodes(md.mesh.x,md.mesh.y,'Exp/terminus_big.exp',1);
45+
% pos=find(terminus & md.mask.ice_levelset<=0);
46+
% md.hydrology.spchead(pos)=0;
47+
48+
md.hydrology.moulin_input = zeros(md.mesh.numberofvertices,1); % No moulin inputs
49+
md.hydrology.neumannflux=zeros(md.mesh.numberofelements,1); % No-flux b.c. on domain boundary
50+
51+
% Coupling and friction
52+
md.transient=deactivateall(md.transient);
53+
md.transient.isstressbalance=1; % 1=Solve for ice velocity, 0=turn off stress balance
54+
md.transient.ishydrology=1;
55+
56+
% Friction
57+
Neff = md.materials.rho_ice*md.constants.g*md.geometry.thickness-md.materials.rho_water*md.constants.g*(md.hydrology.head - md.geometry.base); %Initial effective pressure
58+
md.friction.effective_pressure=Neff;
59+
md.friction.coupling = 4; % Change this to couple (4) /uncouple with prescribed N (3) ***
60+
61+
% Specify that you want to run the model on your current computer
62+
% Change the number of processors according to your machine
63+
md.cluster=generic('np',8);
64+
65+
% Define the time stepping scheme
66+
md.timestepping.time_step=1*3600/md.constants.yts; % Time step (in years)
67+
md.timestepping.final_time=30/365; % Final time (in years)
68+
md.settings.output_frequency=24; % To save model output less frequently than every time step
69+
70+
md.verbose.solution=1;
71+
72+
md=solve(md,'Transient');

0 commit comments

Comments
 (0)