From fa33a8e95368b07c5dc849d66d8cd68ca9515570 Mon Sep 17 00:00:00 2001 From: aleahsommers Date: Mon, 6 Jan 2025 10:35:48 -0500 Subject: [PATCH] ADD:adding Helheim SHAKTI tutorials --- examples/Helheim/Greenland.par | 94 ++++++++++++++++++++++++++++++++ examples/Helheim/Helheim_tut.exp | 20 +++++++ examples/Helheim/runme.m | 86 +++++++++++++++++++++++++++++ examples/HelheimSHAKTI/runme.m | 72 ++++++++++++++++++++++++ 4 files changed, 272 insertions(+) create mode 100644 examples/Helheim/Greenland.par create mode 100644 examples/Helheim/Helheim_tut.exp create mode 100644 examples/Helheim/runme.m create mode 100644 examples/HelheimSHAKTI/runme.m diff --git a/examples/Helheim/Greenland.par b/examples/Helheim/Greenland.par new file mode 100644 index 000000000..156a1226f --- /dev/null +++ b/examples/Helheim/Greenland.par @@ -0,0 +1,94 @@ +%required Data +accpath = ['/home/ModelData/Greenland/RACMO2Accumulation/SMBGreenland/smb_1961-1990.mat']; +velpath = ['/home/ModelData/Greenland/VelMouginot/RignotGreenland2012Vel.mat']; + +%WARNING: we do have a levelset here, so we must set it to one even for stressbalances... +md.transient.ismovingfront=1; + +disp(' Interpolating mask'); +mask = int8(interpBedmachineGreenland(md.mesh.x,md.mesh.y,'mask')); +md.mask.ice_levelset= -1*ones(md.mesh.numberofvertices,1); +pos = find(mask<1); +md.mask.ice_levelset(pos)=1; + +disp(' reading MC bed (assumes no floating ice)'); +md.geometry.bed = interpBedmachineGreenland(md.mesh.x,md.mesh.y,'bed'); +md.geometry.base = md.geometry.bed; + +disp(' reading Howat surface'); +md.geometry.surface=interpBedmachineGreenland(md.mesh.x,md.mesh.y,'surface'); +pos = find(md.mask.ice_levelset>0); +md.geometry.surface(pos) = md.geometry.base(pos)+10; %Minimum thickness + +md.geometry.thickness = md.geometry.surface - md.geometry.bed; +pos=find(md.geometry.thickness<=10); +md.geometry.surface(pos) = md.geometry.base(pos)+10; %Minimum thickness +md.geometry.thickness = md.geometry.surface - md.geometry.bed; + +md.masstransport.min_thickness = 10; + +disp(' Adjusting ice mask'); +%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 +pos = find(max(md.mask.ice_levelset(md.mesh.elements),[],2)>0); +md.mask.ice_levelset(md.mesh.elements(pos,:)) = 1; +% For the region where surface is NaN, set thickness to small value (consistency requires >0) +pos=find((md.mask.ice_levelset<0).*(md.geometry.surface<0)); +md.mask.ice_levelset(pos)=1; +pos=find((md.mask.ice_levelset<0).*(isnan(md.geometry.surface))); +md.mask.ice_levelset(pos)=1; + +disp(' -- reconstruct thickness'); +md.geometry.thickness=md.geometry.surface-md.geometry.base; + +disp(' reading velocities '); +[md.inversion.vx_obs md.inversion.vy_obs]=interpJoughinCompositeGreenland(md.mesh.x,md.mesh.y); +pos=find(isnan(md.inversion.vx_obs) | isnan(md.inversion.vy_obs)); +md.inversion.vx_obs(pos)=0; +md.inversion.vy_obs(pos)=0; +md.inversion.vel_obs = sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2); +md.initialization.vx = md.inversion.vx_obs; +md.initialization.vy = md.inversion.vy_obs; +md.initialization.vz = zeros(md.mesh.numberofvertices,1); +md.initialization.vel = md.inversion.vel_obs; + +%drag md.drag or stress +md.friction.coefficient = 20*ones(md.mesh.numberofvertices,1); %q = 1. +md.friction.p = ones(md.mesh.numberofelements,1); +md.friction.q = ones(md.mesh.numberofelements,1); + +%No friction on PURELY ocean element +pos_e = find(min(md.mask.ice_levelset(md.mesh.elements),[],2)<0); +flags=ones(md.mesh.numberofvertices,1); +flags(md.mesh.elements(pos_e,:))=0; +md.friction.coefficient(find(flags))=0; + +%flow law +disp(' Creating flow law parameters (assume ice is at 0°C for now)'); +md.materials.rheology_n = 3*ones(md.mesh.numberofelements,1); +md.materials.rheology_B = cuffey(273.15 - 10)*ones(md.mesh.numberofvertices,1); + +md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1); + +disp(' Geothermal flux from Shapiro et al.'); +md.basalforcings.geothermalflux=interpSeaRISE(md.mesh.x,md.mesh.y,'bheatflx'); + +disp(' Setting up thermal model'); +md.initialization.temperature=min(0,interpSeaRISE(md.mesh.x,md.mesh.y,'surftemp'))+273.15; +md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1); +md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1); +md.thermal.spctemperature=md.initialization.temperature; +md.thermal.isenthalpy=1; +md.thermal.isdynamicbasalspc=1; + +%Deal with boundary conditions: +disp(' Set Boundary conditions'); +md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1); +md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1); +md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1); +md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6); +md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3); +pos=find((md.mask.ice_levelset<0).*(md.mesh.vertexonboundary)); +md.stressbalance.spcvx(pos)=md.initialization.vx(pos); +md.stressbalance.spcvy(pos)=md.initialization.vy(pos); +md.stressbalance.spcvz(pos)=0; diff --git a/examples/Helheim/Helheim_tut.exp b/examples/Helheim/Helheim_tut.exp new file mode 100644 index 000000000..20ba131cd --- /dev/null +++ b/examples/Helheim/Helheim_tut.exp @@ -0,0 +1,20 @@ +## Name:./Exp/Domain_TEMP.exp +## Icon:0 +# Points Count Value +14 1.000000 +# X pos Y pos +251808.9599315871 -2519299.1448521991 +232017.0235357451 -2526258.4544728124 +232815.0354161858 -2578466.8389138547 +256982.9213724715 -2601667.2788287275 +302263.2150846492 -2594953.0557058132 +309967.1091796352 -2586747.6107035289 +315046.5997917140 -2581507.9891960463 +313861.3853155622 -2571621.9108800408 +310644.3745945790 -2567964.0619031191 +309120.5274109554 -2557188.2365386733 +306496.1239280481 -2552245.1973806708 +305818.8585131043 -2537910.3838224630 +285125.3861979210 -2517350.5381584275 +251808.9599315871 -2519299.1448521991 + diff --git a/examples/Helheim/runme.m b/examples/Helheim/runme.m new file mode 100644 index 000000000..a28cc8b45 --- /dev/null +++ b/examples/Helheim/runme.m @@ -0,0 +1,86 @@ +% runme script to perform a stress balance inversion for basal friction +% on Helheim Glacier (SE Greenland) + +% Step 1: Set up model domain outline and mesh +% Step 2: Set parameters +% Step 3: Inversion for basal friction + +% Note that Steps 1 and 2 require previously downloading data for velocities, +% surface elevation and bed elevation. + +steps = [1:3]; % Choose which steps to run here + +cluster=generic('name',oshostname(),'np',8); + +org=organizer('repository',['./Models'],'prefix',['Model_Helheim_'],'steps',steps); clear steps; + +if perform(org,'Mesh'),% {{{ + + % Initialize uniform mesh + md=triangle(model,['Helheim_tut.exp'],500); + + [velx vely]=interpJoughinCompositeGreenland(md.mesh.x,md.mesh.y); + vel = sqrt(velx.^2+vely.^2); + + % Refine mesh based on surface velocities + md=bamg(md,'hmin',100,'hmax',1500,'field',vel,'err',5); + [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,45,70); + md.mesh.epsg=3413; + + savemodel(org,md); +end %}}} +if perform(org,'Param'),% {{{ + + md=loadmodel(org,'Mesh'); + md=setflowequation(md,'SSA','all'); % Set flow law to SSA + + md=setmask(md,'',''); + md=parameterize(md,'Greenland.par'); + md.miscellaneous.name = 'Helheim'; + + savemodel(org,md); +end%}}} +if perform(org,'Inversion_drag'),% {{{ + + md=loadmodel(org,'Param'); + + %Control general + md.inversion=m1qn3inversion(md.inversion); + md.inversion.iscontrol=1; + md.verbose=verbose('solution',false,'control',true); + md.transient.amr_frequency = 0; + + %Cost functions + md.inversion.cost_functions=[101 103 501]; + md.inversion.cost_functions_coefficients=zeros(md.mesh.numberofvertices,numel(md.inversion.cost_functions)); + md.inversion.cost_functions_coefficients(:,1)=5000; + md.inversion.cost_functions_coefficients(:,2)=10; + md.inversion.cost_functions_coefficients(:,3)=2*50^-3; + pos=find(md.mask.ice_levelset>0); + md.inversion.cost_functions_coefficients(pos,1:2)=0; + + %Controls + md.inversion.control_parameters={'FrictionCoefficient'}; + md.inversion.maxsteps=100; + md.inversion.maxiter =100; + md.inversion.min_parameters=0.05*ones(md.mesh.numberofvertices,1); + md.inversion.max_parameters=4000*ones(md.mesh.numberofvertices,1); + md.inversion.control_scaling_factors=1; + + %Additional parameters + md.stressbalance.restol=0.01; + md.stressbalance.reltol=0.1; + md.stressbalance.abstol=NaN; + + %Go solve + md.cluster=cluster; + + md=solve(md,'sb'); + + %Put results back into the model + md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient; + md.initialization.vx=md.results.StressbalanceSolution.Vx; + md.initialization.vy=md.results.StressbalanceSolution.Vy; + + savemodel(org,md); +end%}}} \ No newline at end of file diff --git a/examples/HelheimSHAKTI/runme.m b/examples/HelheimSHAKTI/runme.m new file mode 100644 index 000000000..6c11dcaea --- /dev/null +++ b/examples/HelheimSHAKTI/runme.m @@ -0,0 +1,72 @@ +% Coupled SHAKTI-ISSM simulation of Helheim, starting from previous inversion +clear all;close all + +% Load Helheim Glacier model +load Models/Model_Helheim_Inversion_drag.mat + +% Turn off inversion +md.inversion.iscontrol=0; + +% HYDROLOGY SPECIFIC PARAMETERIZATION: +% Change hydrology class to SHAKTI model +md.hydrology=hydrologyshakti(); + +% Define distributed englacial input to the subglacial system (m/yr) +md.hydrology.englacial_input = .0*ones(md.mesh.numberofvertices,1); + +% Define initial water head such that water pressure is 80% of ice overburden pressure +md.hydrology.head = 0.8*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base; + +% Initial subglacial gap height of 0.01m everywhere +md.hydrology.gap_height = 0.01*ones(md.mesh.numberofelements,1); + +% Typical bed bump bump spacing -- must be non-zero +md.hydrology.bump_spacing = 1.0*ones(md.mesh.numberofelements,1); + +% Typical bed bump height -- set to zero to turn off opening by sliding +md.hydrology.bump_height = 0.0*ones(md.mesh.numberofelements,1); + +% Initial Reynolds number (start at Re=1000 everywhere) +md.hydrology.reynolds= 1000*ones(md.mesh.numberofelements,1); + +% Boundary conditions +md.hydrology.spchead = NaN(md.mesh.numberofvertices,1); + +% Set head=0 for thin ice, everywhere <=10m ice thickness +pos=find(md.geometry.thickness==10); +md.hydrology.spchead(pos)=0; + +% Set pressure in fjord equal to hydrostatic pressure of fjord water +pos=find(md.mask.ice_levelset>0); +md.hydrology.spchead(pos)=0; + +% % Deal with terminus b.c. in SHAKTI +% terminus=ContourToNodes(md.mesh.x,md.mesh.y,'Exp/terminus_big.exp',1); +% pos=find(terminus & md.mask.ice_levelset<=0); +% md.hydrology.spchead(pos)=0; + +md.hydrology.moulin_input = zeros(md.mesh.numberofvertices,1); % No moulin inputs +md.hydrology.neumannflux=zeros(md.mesh.numberofelements,1); % No-flux b.c. on domain boundary + +% Coupling and friction +md.transient=deactivateall(md.transient); +md.transient.isstressbalance=1; % 1=Solve for ice velocity, 0=turn off stress balance +md.transient.ishydrology=1; + +% Friction +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 +md.friction.effective_pressure=Neff; +md.friction.coupling = 4; % Change this to couple (4) /uncouple with prescribed N (3) *** + +% Specify that you want to run the model on your current computer +% Change the number of processors according to your machine +md.cluster=generic('np',8); + +% Define the time stepping scheme +md.timestepping.time_step=1*3600/md.constants.yts; % Time step (in years) +md.timestepping.final_time=30/365; % Final time (in years) +md.settings.output_frequency=24; % To save model output less frequently than every time step + +md.verbose.solution=1; + +md=solve(md,'Transient');