Index: /issm/trunk-jpl/examples/AMR/runme.m =================================================================== --- /issm/trunk-jpl/examples/AMR/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/AMR/runme.m (revision 25915) @@ -1,13 +1,14 @@ % Mismip3D experiment with AMR using BAMG -steps=[1:3]; +steps=[1:3]; if any(steps==1) disp(' Step 1: Coarse mesh'); - + %Generate an unstructured coarse mesh on the MISMIP domain with typical element edge length equal to 10,000 m md=bamg(model,'domain','./domain.exp','hmax',10000,'splitcorners',1); save AMRCoarseMesh md -end +end + if any(steps==2) disp(' Step 2: Parameterization'); @@ -16,5 +17,5 @@ md=setmask(md,'',''); - + % Run parameterization script to set up geometry, inital velocity, material properties, etc. md=parameterize(md,'./mismip.par'); @@ -23,5 +24,5 @@ % Here, we are refining around the grounding line % We impose the element resolution at grounding equal to 1000 m (1 km) - % The criterion used is the element distance to the grounding line. + % The criterion used is the element distance to the grounding line % The distance used here is 10000 m (10 km), used in both side around the grouding line (upstream and downstream) md.amr.groundingline_resolution=1000; @@ -35,7 +36,8 @@ save AMRParam md end -if any(steps==3); - disp(' Step 3: Solve!'); - + +if any(steps==3) + disp(' Step 3: Solve!'); + md=loadmodel('AMRParam'); @@ -43,10 +45,11 @@ md.timestepping.time_step=1; md.timestepping.final_time=500; % here, as example, only 500 yr. - md.settings.output_frequency=10;% here, save results every 10 yr + md.settings.output_frequency=10;% here, save results every 10 yr md.stressbalance.maxiter=30; md.stressbalance.abstol=NaN; md.stressbalance.restol=1; + md.settings.solver_residue_threshold=1e-2; % relaxing (the first stress balance solver iteration presents values higher than the original threshold. This probably happens because the initial velocity is set to one). md.verbose=verbose('convergence',false,'solution',true); - + % Specify that you want to run the model on your current (local host) computer % Change the number of processors according to your machine (here np=2) Index: /issm/trunk-jpl/examples/EsaGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/EsaGRACE/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/EsaGRACE/runme.m (revision 25915) @@ -1,58 +1,54 @@ - clear all; addpath('../Data','../Functions'); -steps=[1]; % [1:5] +steps=[1]; % [1:5] -if any(steps==1) +if any(steps==1) %{{{ disp(' Step 1: Global mesh creation'); - resolution=300; % [km] - radius = 6.371012*10^3; % [km] + resolution=300; % [km] + radius = 6.371012*10^3; % [km] - md=model; - md.mask=maskpsl(); + md=model; md.mesh=gmshplanet('radius',radius,'resolution',resolution); - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); save ./Models/EsaGRACE_Mesh md; - + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); +end %}}} -end - -if any(steps==2) +if any(steps==2) %{{{ disp(' Step 2: Define loads in meters of ice height equivalent'); md = loadmodel('./Models/EsaGRACE_Mesh'); year_month = 2007+15/365; - time_range = [year_month year_month]; - - water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - - md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent + time_range = [year_month year_month]; - save ./Models/EsaGRACE_Loads md; - + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent + + save ./Models/EsaGRACE_Loads md; + plotmodel (md,'data',md.esa.deltathickness,... 'view',[90 -90],'caxis',[-.1 .1],... 'title','Ice height equivalent [m]'); +end %}}} -end - -if any(steps==3) +if any(steps==3) %{{{ disp(' Step 3: Parameterization'); md = loadmodel('./Models/EsaGRACE_Loads'); - nlove=10001; - md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; - md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; + love_numbers = lovenumbers('maxdeg',10000,'referenceframe','CF'); + md.esa.love_h = love_numbers.h; + md.esa.love_l = love_numbers.l; - md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); - md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); pos=find(md.esa.deltathickness~=0); - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - md.mask.land_levelset = 1-md.mask.ocean_levelset; + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + %md.mask.land_levelset = 1-md.mask.ocean_levelset; di=md.materials.rho_ice/md.materials.rho_water; @@ -61,52 +57,50 @@ md.geometry.base=md.geometry.surface-md.geometry.thickness; md.geometry.bed=md.geometry.base; - + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - + md.miscellaneous.name='EsaGRACE'; - - save ./Models/EsaGRACE_Parameterization md; -end + save ./Models/EsaGRACE_Parameterization md; +end %}}} -if any(steps==4) +if any(steps==4) %{{{ disp(' Step 4: Solve Esa solver'); md = loadmodel('./Models/EsaGRACE_Parameterization'); md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'}; - - md.cluster=generic('name',oshostname(),'np',3); + + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); md=solve(md,'Esa'); - save ./Models/EsaGRACE_Solution md; + save ./Models/EsaGRACE_Solution md; +end %}}} -end - -if any(steps==5) +if any(steps==5) %{{{ disp(' Step 5: Plot solutions'); md = loadmodel('./Models/EsaGRACE_Solution'); - sol1 = md.esa.deltathickness*100; % [cm] - sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] - sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] - sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] + sol1 = md.esa.deltathickness*100; % [cm] + sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] + sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] + sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... - 'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; - fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'}; + 'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; + fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'}; - res = 1.0; % [degree] + res = 1.0; % [degree] [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); - sol_grid = zeros(size(lat_grid)); + sol_grid = zeros(size(lat_grid)); - for kk=1:4 + for kk=1:4 sol=eval(sprintf('sol%d',kk)); - - if length(sol)==md.mesh.numberofelements + + if length(sol)==md.mesh.numberofelements for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -114,37 +108,35 @@ end for jj=1:md.mesh.numberofvertices - pos=ceil(find(pp==jj)/3); - temp(jj)=mean(sol(pos)); + pos=ceil(find(pp==jj)/3); + temp(jj)=mean(sol(pos)); end - sol=temp'; - end + sol=temp'; + end - F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); F.Method = 'linear'; - F.ExtrapolationMethod = 'linear'; + F.ExtrapolationMethod = 'linear'; sol_grid = F(lat_grid, lon_grid); - sol_grid(isnan(sol_grid))=0; - sol_grid(lat_grid>85 & sol_grid==0)=NaN; + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0)=NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) - figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; - pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; if (kk==1) - geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); - end - plot(long,lat,'k'); hold off; + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + end + plot(long,lat,'k'); hold off; c1=colorbar; colormap('haxby'); - caxis([-min(abs(min(sol)),abs(max(sol))) min(abs(min(sol)),abs(max(sol)))]); - xlim([-180 180]); - ylim([-90 90]); - grid on; - title(sol_name(kk)); + caxis([-min(abs(min(sol)),abs(max(sol))) min(abs(min(sol)),abs(max(sol)))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); set(gcf,'color','w'); - %export_fig(fig_name{kk}); + %export_fig(fig_name{kk}); end - -end - +end %}}} Index: /issm/trunk-jpl/examples/EsaWahr/runme.m =================================================================== --- /issm/trunk-jpl/examples/EsaWahr/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/EsaWahr/runme.m (revision 25915) @@ -1,67 +1,63 @@ - clear all; addpath('../Functions'); -steps=[0]; % [0:6] +steps=[1]; % [1:7] -if any(steps==0) - disp(' Step 0: Mesh creation'); +if any(steps==1) %{{{ + disp(' Step 1: Mesh creation'); - md=roundmesh(model,100000,10000); % Domain radius and element size [m] + md=roundmesh(model,100000,10000); % Domain radius and element size [m] save ./Models/EsaWahr_Mesh md; - + plotmodel(md,'data','mesh'); +end %}}} -end +if any(steps==2) %{{{ + disp(' Step 2: Anisotropic mesh creation'); -if any(steps==1) - disp(' Step 1: Anisotropic mesh creation'); + md=roundmesh(model,100000,1000); - md=roundmesh(model,100000,1000); + disc_radius=20*1000; + rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); + field = abs(rad_dist-disc_radius); - disc_radius=20*1000; - rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); - field = abs(rad_dist-disc_radius); - - md = bamg(md,'field',field,'err',50,'hmax',10000); + md = bamg(md,'field',field,'err',50,'hmax',10000); save ./Models/EsaWahr_Mesh md; - + plotmodel (md,'data','mesh'); +end %}}} -end - -if any(steps==2) - disp(' Step 2: Define loads'); +if any(steps==3) %{{{ + disp(' Step 3: Define loads'); md = loadmodel('./Models/EsaWahr_Mesh'); - rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; + rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; index=md.mesh.elements; - x_cent=mean(md.mesh.x(index),2); - y_cent=mean(md.mesh.y(index),2); + x_cent=mean(md.mesh.x(index),2); + y_cent=mean(md.mesh.y(index),2); - md.esa.deltathickness = zeros(md.mesh.numberofelements,1); - disc_radius=20; % [km] + md.esa.deltathickness = zeros(md.mesh.numberofelements,1); + disc_radius=20; % [km] rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; - save ./Models/EsaWahr_Loads md; - + save ./Models/EsaWahr_Loads md; + plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]'); +end %}}} -end - -if any(steps==3) - disp(' Step 3: Parameterization'); +if any(steps==4) %{{{ + disp(' Step 4: Parameterization'); md = loadmodel('./Models/EsaWahr_Loads'); - nlove=10001; - md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; - md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; + love_numbers = lovenumbers('maxdeg',10000,'referenceframe','CF'); + md.esa.love_h = love_numbers.h; + md.esa.love_l = love_numbers.l; - md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); - md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); + md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); di=md.materials.rho_ice/md.materials.rho_water; @@ -70,40 +66,38 @@ md.geometry.base=md.geometry.surface-md.geometry.thickness; md.geometry.bed=md.geometry.base; - + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - + md.miscellaneous.name='EsaWahr'; - - save ./Models/EsaWahr_Parameterization md; -end + save ./Models/EsaWahr_Parameterization md; +end %}}} -if any(steps==4) - disp(' Step 4: Solve Esa solver'); +if any(steps==5) %{{{ + disp(' Step 5: Solve Esa solver'); md = loadmodel('./Models/EsaWahr_Parameterization'); md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'}; - - md.cluster=generic('name',oshostname(),'np',3); + + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); md=solve(md,'Esa'); - save ./Models/EsaWahr_Solution md; + save ./Models/EsaWahr_Solution md; +end %}}} -end - -if any(steps==5) - disp(' Step 5: Plot solutions'); +if any(steps==6) %{{{ + disp(' Step 6: Plot solutions'); md = loadmodel('./Models/EsaWahr_Solution'); - vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] - horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] - horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] - horz = sqrt(horz_n.^2+horz_e.^2); % [mm] + vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] + horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] + horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] + horz = sqrt(horz_n.^2+horz_e.^2); % [mm] - set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); + set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); figure('Position', [100, 100, 800, 600]); plotmodel(md,'data',vert,... @@ -125,49 +119,46 @@ 'caxis#3-4',[-0.5 0.5],... 'axispos',[0.505 0.02 0.47 0.47],... - 'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); - %export_fig('Fig5.pdf'); + 'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); + %export_fig('Fig5.pdf'); +end %}}} -end - -if any(steps==6) - disp(' Step 6: Compare results against Wahr semi-analytic solutions'); +if any(steps==7) %{{{ + disp(' Step 7: Compare results against Wahr semi-analytic solutions'); md = loadmodel('./Models/EsaWahr_Solution'); - vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] - horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] - horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] - horz = sqrt(horz_n.^2+horz_e.^2); % [mm] - - xi=[0:500:100000]; % grid points [m] - yi=zeros(1,length(xi)); - vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); - horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); + vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] + horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] + horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] + horz = sqrt(horz_n.^2+horz_e.^2); % [mm] - % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) - disc_radius = 20*1000; % [m] + xi=[0:500:100000]; % grid points [m] + yi=zeros(1,length(xi)); + vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); + horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); + + % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) + disc_radius = 20*1000; % [m] [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l); - set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); + set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); figure1=figure('Position', [100, 100, 700, 400]); - ylabel_1={'0',' ','1','','2','','3',''}; + ylabel_1={'0',' ','1','','2','','3',''}; axes1 = axes('Layer','top','Position',[0.1 0.15 0.8 0.8],... 'XTick',[0:10:100],'xlim',[0 100],... - 'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); - box(axes1,'on'); hold(axes1,'all'); grid on; - xlabel(axes1,'Radial distance [km]'); + 'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); + box(axes1,'on'); hold(axes1,'all'); grid on; + xlabel(axes1,'Radial distance [km]'); ylabel(axes1,'Displacement [mm]'); - plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1); - % analytic soln - h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); - h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); - % ISSM soln - h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); - h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); + plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1); + % analytic soln + h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); + h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); + % ISSM soln + h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); + h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); ag1 = gca; leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)'); - set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); + set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); set(gcf,'color','w'); - %export_fig('Fig6.pdf'); - -end - + %export_fig('Fig6.pdf'); +end %}}} Index: /issm/trunk-jpl/examples/Functions/grace.m =================================================================== --- /issm/trunk-jpl/examples/Functions/grace.m (revision 25914) +++ /issm/trunk-jpl/examples/Functions/grace.m (revision 25915) @@ -43,9 +43,15 @@ % Monthly GRACE data - filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc']; - time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC - long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75 - lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75 - rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm] + filename=['../Data/GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc']; + try + time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC + long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75 + lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75 + rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm] + catch e + disp('If dataset file exists at another location, modify the path in examples/Functions/grace.m') + rethrow(e) + end + time_yr = 2002+time_0/365; % [yr] Index: /issm/trunk-jpl/examples/Greenland/runme.m =================================================================== --- /issm/trunk-jpl/examples/Greenland/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/Greenland/runme.m (revision 25915) @@ -5,5 +5,5 @@ ncdata='../Data/Greenland_5km_dev1.2.nc'; -if any(steps==1) +if any(steps==1) %{{{ disp(' Step 1: Mesh creation'); @@ -22,5 +22,5 @@ %Mesh Greenland md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,'err',8); - + %convert x,y coordinates (Polar stereo) to lat/lon [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,+1,39,71); @@ -29,7 +29,7 @@ plotmodel (md,'data','mesh'); -end - -if any(steps==2) +end %}}} + +if any(steps==2) %{{{ disp(' Step 2: Parameterization'); md = loadmodel('./Models/Greenland.Mesh_generation'); @@ -39,8 +39,8 @@ md = setflowequation(md,'SSA','all'); - save ./Models/Greenland.Parameterization md; -end - -if any(steps==3) + save ./Models/Greenland.Parameterization md; +end %}}} + +if any(steps==3) %{{{ disp(' Step 3: Control method friction'); md = loadmodel('./Models/Greenland.Parameterization'); @@ -66,5 +66,5 @@ %Additional parameters - md.stressbalance.restol=0.01; md.stressbalance.reltol=0.1; + md.stressbalance.restol=0.01; md.stressbalance.reltol=0.1; md.stressbalance.abstol=NaN; @@ -78,8 +78,8 @@ md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient; - save ./Models/Greenland.Control_drag md; -end - -if any(steps==4) + save ./Models/Greenland.Control_drag md; +end %}}} + +if any(steps==4) %{{{ disp(' Step 4: Transient run'); md = loadmodel('./Models/Greenland.Control_drag'); @@ -101,6 +101,5 @@ %Additional options md.inversion.iscontrol=0; - md.transient.requested_outputs={'IceVolume','TotalSmb', ... - 'SmbMassBalance'}; + md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'}; md.verbose=verbose('solution',true,'module',true,'convergence',true); @@ -109,9 +108,9 @@ md=solve(md,'Transient'); - save ./Models/Greenland.Transient md; -end - -if any(steps==5) - disp(' Step 5: Plotting'); + save ./Models/Greenland.Transient md; +end %}}} + +if any(steps==5) %{{{ + disp(' Step 5: Plotting'); md = loadmodel('./Models/Greenland.Transient'); @@ -143,7 +142,7 @@ subplot(3,1,3); plot([0.2:0.2:20],volume); title('Ice Volume'); xlabel('years') -end - -if any(steps==6) +end %}}} + +if any(steps==6) %{{{ disp(' Step 6: Extract Box SMB'); md = loadmodel('./Models/Greenland.Transient'); @@ -172,8 +171,7 @@ clear smbbox - -end - -if any(steps==7) +end %}}} + +if any(steps==7) %{{{ disp(' Step 7: Historical Relaxation run'); md = loadmodel('./Models/Greenland.Control_drag'); @@ -198,6 +196,5 @@ %Additional options md.inversion.iscontrol=0; - md.transient.requested_outputs={'IceVolume','TotalSmb', ... - 'SmbMassBalance'}; + md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'}; md.verbose=verbose('solution',true,'module',true); @@ -207,45 +204,35 @@ save ./Models/Greenland.HistoricTransient md; -end - -if any(steps==8) +end % step 7 end + +if any(steps==8) %{{{ disp(' Step 8: Plotting exercise'); %Load historic transient model - - %Create Line Plots of relaxation run. Create a figure. - + %Create Line Plots of relaxation run. Create a figure. %Save surface mass balance, by looping through 200 years (1000 steps) - % Note, the first output will always contain output from time step 1 - + %Note, the first output will always contain output from time step 1 %Plot surface mass balance time series in first subplot - %Title this plot Mean surface mass balance - %Save velocity by looping through 200 years - %Plot velocity time series in second subplot - %Title this plot Mean Velocity - - %Save Ice Volume by looping through 200 years - + %Save Ice Volume by looping through 200 years %Plot volume time series in third subplot - %Title this plot Mean Velocity and add an x label of years -end - -if any(steps==9) +end % step 8 end + +if any(steps==9) %{{{ disp(' Step 9: Box Transient run'); md = loadmodel('./Models/Greenland.HistoricTransient_200yr'); @@ -292,6 +279,5 @@ %Additional options md.inversion.iscontrol=0; - md.transient.requested_outputs={'IceVolume','TotalSmb', ... - 'SmbMassBalance'}; + md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'}; md.verbose=verbose('solution',true,'module',true); @@ -301,7 +287,7 @@ save ./Models/Greenland.BoxTransient md; -end - -if any(steps==10) +end %}}} + +if any(steps==10) %{{{ disp(' Step 10: Plot Box Transient'); md = loadmodel('./Models/Greenland.BoxTransient'); @@ -327,3 +313,3 @@ subplot(3,1,3); plot(t,volume); title('Ice Volume'); xlabel('years') -end +end %}}} Index: sm/trunk-jpl/examples/ISMIP/CheatyIsmipA.par =================================================================== --- /issm/trunk-jpl/examples/ISMIP/CheatyIsmipA.par (revision 25914) +++ (revision ) @@ -1,52 +1,0 @@ -%Parameterization for ISMIP A experiment - -%Set the Simulation generic name #md.miscellaneous -%-> - -%Geometry -disp(' Constructing Geometry'); - -%Define the geometry of the simulation #md.geometry -%surface is [-x*tan(0.5*pi/180)] #md.mesh -%-> -md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.); -%base is [surface-1000+500*sin(x*2*pi/L).*sin(y*2*pi/L)] -%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x) -%-> -L=max(md.mesh.x)-min(md.mesh.x); -md.geometry.base=md.geometry.surface-1000.0+500.0*sin(md.mesh.x*2.0*pi/L).*sin(md.mesh.y*2.0*pi/L); -%thickness is the difference between surface and base #md.geometry -%-> -md.geometry.thickness=md.geometry.surface-md.geometry.base; -%plot the geometry to check it out -%-> -plotmodel(md,'data',md.geometry.thickness); - -disp(' Defining friction parameters'); - -%These parameters will not be used but need to be fixed #md.friction -%one friciton coefficient per node (md.mesh.numberofvertices,1) -%-> -md.friction.coefficient=200.0*ones(md.mesh.numberofvertices,1); -%one friciton exponent (p,q) per element -%-> -md.friction.p=ones(md.mesh.numberofelements,1); -%-> -md.friction.q=ones(md.mesh.numberofelements,1); - -disp(' Construct ice rheological properties'); - -%The rheology parameters sit in the material section #md.materials -%B has one value per vertex -%-> -md.materials.rheology_B=6.8067e7*ones(md.mesh.numberofvertices,1); -%n has one value per element -%-> -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - -disp(' Set boundary conditions'); - -%Set the default boundary conditions for an ice-sheet -% #help SetIceSheetBC -%-> -md=SetIceSheetBC(md); Index: sm/trunk-jpl/examples/ISMIP/CheatyIsmipF.par =================================================================== --- /issm/trunk-jpl/examples/ISMIP/CheatyIsmipF.par (revision 25914) +++ (revision ) @@ -1,66 +1,0 @@ -%Parameterization for ISMIP F experiment - -%Set the Simulation generic name #md.miscellaneous -%-> - -%Geometry -disp(' Constructing Geometry'); - -%Define the geometry of the simulation #md.geometry -%surface is [-x*tan(3.0*pi/180)] #md.mesh -%-> -md.geometry.surface=md.mesh.x*tan(3.0*pi/180.0); -%base is [surface-1000+100*exp(-((x-L/2).^2+(y-L/2).^2)/(10000.^2))] -%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x) -%-> -L=max(md.mesh.x)-min(md.mesh.x); -%-> -md.geometry.base=md.geometry.surface-1000.0+100.0*exp(-((md.mesh.x-L/2.0).^2.0+(md.mesh.y-L/2.0).^2.0)/(10000.^2.0)); -%thickness is the difference between surface and base #md.geometry -%-> -md.geometry.thickness=md.geometry.surface-md.geometry.base; -%plot the geometry to check it out -%-> -plotmodel(md,'data',md.geometry.thickness); - -disp(' Defining friction parameters'); - -%These parameters will not be used but need to be fixed #md.friction -%one friciton coefficient per node (md.mesh.numberofvertices,1) -%conversion form year to seconds with #md.constants.yts -%-> -md.friction.coefficient=sqrt(md.constants.yts/(1000*2.140373*10^-7))*ones(md.mesh.numberofvertices,1); -%one friciton exponent (p,q) per element -%-> -md.friction.p=ones(md.mesh.numberofelements,1); -%-> -md.friction.q=zeros(md.mesh.numberofelements,1); - -disp(' Construct ice rheological properties'); - -%The rheology parameters sit in the material section #md.materials -%B has one value per vertex -%-> -md.materials.rheology_B=(1/(2.140373*10^-7/md.constants.yts))*ones(md.mesh.numberofvertices,1); -%n has one value per element -%-> -md.materials.rheology_n=1*ones(md.mesh.numberofelements,1); - -disp(' Set boundary conditions'); - -%Set the default boundary conditions for an ice-sheet -% #help SetIceSheetBC -%-> -md=SetIceSheetBC(md); - -disp(' Initializing velocity and pressure'); - -%initialize the velocity and pressurefields of #md.initialization -%-> -md.initialization.vx=zeros(md.mesh.numberofvertices,1); -%-> -md.initialization.vy=zeros(md.mesh.numberofvertices,1); -%-> -md.initialization.vz=zeros(md.mesh.numberofvertices,1); -%-> -md.initialization.pressure=zeros(md.mesh.numberofvertices,1); Index: sm/trunk-jpl/examples/ISMIP/CheatyRunme.m =================================================================== --- /issm/trunk-jpl/examples/ISMIP/CheatyRunme.m (revision 25914) +++ (revision ) @@ -1,221 +1,0 @@ -%wich steps to perform steps are from 1 to 8 -%step 7 is specific to ISMIPA -%step 8 is specific to ISMIPF - -steps=[1]; - -% parameter file to be used, choose between CheatyIsmipA.par or CheatyIsmipF.par -ParamFile='CheatyIsmipF.par' - -%Run Steps - -%Mesh Generation #1 -if any(steps==1) - - %initialize md as a new model #help model - %-> - md=model(); - % generate a squaremesh #help squaremesh - % Side is 80 km long with 20 points - %-> - if(ParamFile=='CheatyIsmipA.par'), - md=squaremesh(md,80000,80000,20,20); - elseif(ParamFile=='CheatyIsmipF.par'), - md=squaremesh(md,100000,100000,30,30); - end - % plot the given mesh #plotdoc - %-> - plotmodel(md,'data','mesh') - % save the given model - %-> - save ./Models/ISMIP.Mesh_generation md; -end - -%Masks #2 -if any(steps==2) - - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.Mesh_generation'); - % set the mask #help setmask - % all MISMIP nodes are grounded - %-> - md=setmask(md,'',''); - % plot the given mask #md.mask to locate the field - %-> - plotmodel(md,'data',md.mask.ocean_levelset); - % save the given model - %-> - save ./Models/ISMIP.SetMask md; -end - -%Parameterization #3 -if any(steps==3) - - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.SetMask'); - % parametrize the model # help parameterize - % you will need to fil-up the parameter file defined by the - % ParamFile variable - %-> - md=parameterize(md,ParamFile); - % save the given model - %-> - save ./Models/ISMIP.Parameterization md; -end - -%Extrusion #4 -if any(steps==4) - - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.Parameterization'); - % vertically extrude the preceding mesh #help extrude - % only 5 layers exponent 1 - %-> - md=extrude(md,9,1); - % plot the 3D geometry #plotdoc - %-> - plotmodel(md,'data',md.geometry.base) - % save the given model - %-> - save ./Models/ISMIP.Extrusion md; -end - -%Set the flow computing method #5 -if any(steps==5) - - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.Extrusion'); - % set the approximation for the flow computation #help setflowequation - % We will be using the Higher Order Model (HO) - %-> - md=setflowequation(md,'HO','all'); - % save the given model - %-> - save ./Models/ISMIP.SetFlow md; -end - -%Set Boundary Conditions #6 -if any(steps==6) - - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.SetFlow'); - % dirichlet boundary condition are known as SPCs - % ice frozen to the base, no velocity #md.stressbalance - % SPCs are initialized at NaN one value per vertex - %-> - 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); - % extract the nodenumbers at the base #md.mesh.vertexonbase - %-> - basalnodes=find(md.mesh.vertexonbase); - % set the sliding to zero on the bed - %-> - md.stressbalance.spcvx(basalnodes)=0.0; - %-> - md.stressbalance.spcvy(basalnodes)=0.0; - % periodic boundaries have to be fixed on the sides - % create tabs with the side of the domain - % for x - % create maxX #help find - %-> - maxX=find(md.mesh.x==max(md.mesh.x)); - % create minX - %-> - minX=find(md.mesh.x==min(md.mesh.x)); - % for y, max X and minX should be excluded - % create maxY - %-> - maxY=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x)); - % create minY - %-> - minY=find(md.mesh.y==min(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x)); - % set the node that should be paired together - % #md.stressbalance.vertex_pairing - %-> - md.stressbalance.vertex_pairing=[minX,maxX;minY,maxY]; - if (ParamFile=='CheatyIsmipF.par') - % if we are dealing with IsmipF the solution is in - % masstransport - md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing; - end - % save the given model - %-> - save ./Models/ISMIP.BoundaryCondition md; -end - -%Solving #7 -if any(steps==7) - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.BoundaryCondition'); - % Set cluster #md.cluster - % generic parameters #help generic - % set only the name and number of process - %-> - md.cluster=generic('name',oshostname(),'np',2); - % Set which control message you want to see #help verbose - %-> - md.verbose=verbose('convergence',true); - % Solve #help solve - % we are solving a StressBalanc - %-> - md=solve(md,'Stressbalance'); - % save the given model - %-> - save ./Models/ISMIP.StressBalance md; - % plot the surface velocities #plotdoc - %-> - plotmodel(md,'data',md.results.StressbalanceSolution.Vel) -end - -%Solving #8 -if any(steps==8) - % load the preceding step #help loadmodel - % path is given by the organizer with the name of the given step - %-> - md = loadmodel('./Models/ISMIP.BoundaryCondition'); - % Set cluster #md.cluster - % generic parameters #help generic - % set only the name and number of process - %-> - md.cluster=generic('name',oshostname(),'np',2); - % Set which control message you want to see #help verbose - %-> - md.verbose=verbose('convergence',true); - % set the transient model to ignore the thermal model - % #md.transient - %-> - md.transient.isthermal=0; - % define the timestepping scheme - % everything here should be provided in years #md.timestepping - % give the length of the time_step (4 years) - %-> - md.timestepping.time_step=4; - % give final_time (20*4 years time_steps) - %-> - md.timestepping.final_time=4*20; - % Solve #help solve - % we are solving a TransientSolution - %-> - md=solve(md,'Transient'); - % save the given model - %-> - save ./Models/ISMIP.Transient md; - % plot the surface velocities #plotdoc - %-> - plotmodel(md,'data',md.results.TransientSolution(20).Vel) -end Index: /issm/trunk-jpl/examples/ISMIP/runme.m =================================================================== --- /issm/trunk-jpl/examples/ISMIP/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/ISMIP/runme.m (revision 25915) @@ -1,3 +1,3 @@ -%wich steps to perform steps are from 1 to 8 +%which steps to perform; steps are from 1 to 8 %step 7 is specific to ISMIPA %step 8 is specific to ISMIPF @@ -11,5 +11,5 @@ %Mesh Generation #1 -if any(steps==1) +if any(steps==1) %{{{ %initialize md as a new model #help model @@ -26,8 +26,8 @@ %-> -end +end %}}} %Masks #2 -if any(steps==2) +if any(steps==2) %{{{ % load the preceding step #help loadmodel @@ -45,8 +45,8 @@ %-> -end +end %}}} %Parameterization #3 -if any(steps==3) +if any(steps==3) %{{{ % load the preceding step #help loadmodel @@ -62,8 +62,8 @@ %-> -end +end %}}} %Extrusion #4 -if any(steps==4) +if any(steps==4) %{{{ % load the preceding step #help loadmodel @@ -81,8 +81,8 @@ %-> -end +end %}}} %Set the flow computing method #5 -if any(steps==5) +if any(steps==5) %{{{ % load the preceding step #help loadmodel @@ -97,8 +97,8 @@ %-> -end +end %}}} %Set Boundary Conditions #6 -if any(steps==6) +if any(steps==6) %{{{ % load the preceding step #help loadmodel @@ -118,5 +118,5 @@ %-> - % set the sliding to zero on the bed (Vx and Vy) + % set the sliding to zero on the bed (Vx and Vy) %-> @@ -146,12 +146,12 @@ % if we are dealing with IsmipF the solution is in masstransport md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing; - end - % save the given model - %-> - -end + end + % save the given model + %-> + +end %}}} %Solving #7 -if any(steps==7) +if any(steps==7) %{{{ % load the preceding step #help loadmodel % path is given by the organizer with the name of the given step @@ -175,8 +175,8 @@ % plot the surface velocities #plotdoc %-> -end +end %}}} %Solving #8 -if any(steps==8) +if any(steps==8) %{{{ % load the preceding step #help loadmodel % path is given by the organizer with the name of the given step @@ -212,4 +212,3 @@ % plot the surface velocities #plotdoc %-> - -end +end %}}} Index: /issm/trunk-jpl/examples/IceBridge/Greenland.par =================================================================== --- /issm/trunk-jpl/examples/IceBridge/Greenland.par (revision 25914) +++ /issm/trunk-jpl/examples/IceBridge/Greenland.par (revision 25915) @@ -30,5 +30,5 @@ %Reading IceBridge data for Jakobshavn disp(' reading IceBridge Jakobshavn bedrock'); -fid = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt'); +fid = fopen('../Data/Jakobshavn_2008_2011_Composite/grids/Jakobshavn_2008_2011_Composite_XYZGrid.txt'); titles = fgets(fid); data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])'; Index: /issm/trunk-jpl/examples/IceBridge/runme.m =================================================================== --- /issm/trunk-jpl/examples/IceBridge/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/IceBridge/runme.m (revision 25915) @@ -1,4 +1,3 @@ clear all; -close all; steps=[1]; @@ -6,5 +5,5 @@ ncdata='../Data/Greenland_5km_dev1.2.nc'; -if any(steps==1) +if any(steps==1) %{{{ disp(' Step 1: Mesh creation'); @@ -24,5 +23,5 @@ md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,'err',8); return; - + %Refine mesh in the region of Jakobshavn (resolution = 3000 m) hmaxVertices=NaN*ones(md.mesh.numberofvertices,1); @@ -37,18 +36,18 @@ save ./Models/Greenland.Mesh_generation md; -end +end %}}} -if any(steps==2) +if any(steps==2) %{{{ disp(' Step 2: Parameterization'); md = loadmodel('./Models/Greenland.Mesh_generation'); md = setmask(md,'',''); - md = parameterize(md,'Greenland.par'); + md = parameterize(md,'./Greenland.par'); md = setflowequation(md,'SSA','all'); save ./Models/Greenland.Parameterization2 md; -end +end %}}} -if any(steps==3) +if any(steps==3) %{{{ disp(' Step 3: Control method friction'); md = loadmodel('./Models/Greenland.Parameterization2'); @@ -93,7 +92,7 @@ save ./Models/Greenland.Control_drag md; -end +end %}}} -if any(steps==4) +if any(steps==4) %{{{ disp(' Step 4: Transient run'); md = loadmodel('./Models/Greenland.Control_drag'); @@ -122,7 +121,7 @@ save ./Models/Greenland.Transient md; -end +end %}}} -if any(steps==5) +if any(steps==5) %{{{ disp(' Step 5: Plotting'); md = loadmodel('./Models/Greenland.Transient'); @@ -155,3 +154,3 @@ subplot(3,1,3); plot(time_plot,volume); title('Ice Volume'); xlabel('years') -end +end %}}} Index: /issm/trunk-jpl/examples/Inversion/runme.m =================================================================== --- /issm/trunk-jpl/examples/Inversion/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/Inversion/runme.m (revision 25915) @@ -1,4 +1,5 @@ -step=1; -if any(step==1) +steps=[1]; + +if any(steps==1) %{{{ %Generate observations md = model; @@ -9,9 +10,10 @@ md.cluster = generic('np',2); md = solve(md,'Stressbalance'); - plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'title','"True" B',... + plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[1.3 1.9]*10^8,'title','"True" B',... 'data',md.results.StressbalanceSolution.Vel,'title','"observed velocities"') save model1 md -end -if any(step==2) +end %}}} + +if any(steps==2) %{{{ %Modify rheology, now constant loadmodel('model1.mat'); @@ -20,14 +22,15 @@ %results of previous run are taken as observations md.inversion=m1qn3inversion(); - md.inversion.vx_obs = md.results.StressbalanceSolution.Vx; - md.inversion.vy_obs = md.results.StressbalanceSolution.Vy; - md.inversion.vel_obs = md.results.StressbalanceSolution.Vel; + md.inversion.vx_obs = md.results.StressbalanceSolution.Vx; + md.inversion.vy_obs = md.results.StressbalanceSolution.Vy; + md.inversion.vel_obs = md.results.StressbalanceSolution.Vel; md = solve(md,'Stressbalance'); - plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'title','B first guess',... + plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[1.3 1.9]*10^8,'title','B first guess',... 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities') save model2 md -end -if any(step==3) +end %}}} + +if any(steps==3) %{{{ %invert for ice rigidity loadmodel('model2.mat'); @@ -40,6 +43,6 @@ md.inversion.cost_functions = 101; md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1); - md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); - md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); + md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); + md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); %Go solve! @@ -48,6 +51,7 @@ plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',... 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities') -end -if any(step==4) +end %}}} + +if any(steps==4) %{{{ %invert for ice rigidity loadmodel('model2.mat'); @@ -59,8 +63,8 @@ md.inversion.maxsteps = maxsteps; md.inversion.cost_functions = [101 502]; - md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1); - md.inversion.cost_functions_coefficients(:,2) = 10^-16*ones(md.mesh.numberofvertices,1); - md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); - md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); + md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1); + md.inversion.cost_functions_coefficients(:,2) = 10^-16*ones(md.mesh.numberofvertices,1); + md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); + md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); %Go solve! @@ -69,3 +73,3 @@ plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',... 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities') -end +end %}}} Index: /issm/trunk-jpl/examples/Jakobshavn/runme.m =================================================================== --- /issm/trunk-jpl/examples/Jakobshavn/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/Jakobshavn/runme.m (revision 25915) @@ -1,38 +1,40 @@ steps=[1]; -if any(steps==1) - disp(' Step 1: Mesh creation'); +if any(steps==1) %{{{ + disp(' Step 1: Mesh creation'); md=triangle(model,'Domain.exp',2000); %Get observed velocity field on mesh nodes ncdata='../Data/Greenland_5km_dev1.2.nc'; - if ~exist(ncdata,'file'), + if ~exist(ncdata,'file'), error('File Greenland_5km_dev1.2.nc not downloaded in Data Directory'); end - x1 = ncread(ncdata,'x1'); - y1 = ncread(ncdata,'y1'); - velx = ncread(ncdata,'surfvelx'); - vely = ncread(ncdata,'surfvely'); - vx = InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0); - vy = InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0); - vel = sqrt(vx.^2+vy.^2); + x1 = ncread(ncdata,'x1'); + y1 = ncread(ncdata,'y1'); + velx = ncread(ncdata,'surfvelx'); + vely = ncread(ncdata,'surfvely'); + vx = InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0); + vy = InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0); + vel = sqrt(vx.^2+vy.^2); %refine mesh using surface velocities as metric md=bamg(md,'hmin',1200,'hmax',15000,'field',vel,'err',5); [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,39,71); - + save JksMesh md -end -if any(steps==2) - disp(' Step 2: Parameterization'); +end %}}} + +if any(steps==2) %{{{ + disp(' Step 2: Parameterization'); md=loadmodel('JksMesh'); - + md=setmask(md,'',''); - md=parameterize(md,'Jks.par'); + md=parameterize(md,'Jks.par'); save JksPar md -end -if any(steps==3) - disp(' Step 3: Control method friction'); +end %}}} + +if any(steps==3) %{{{ + disp(' Step 3: Control method friction'); md=loadmodel('JksPar'); @@ -66,9 +68,10 @@ md.cluster=generic('name',oshostname,'np',4); md=solve(md,'Stressbalance'); - + save JksControl md -end -if any(steps==4) - disp(' Plotting') +end %}}} + +if any(steps==4) %{{{ + disp(' Plotting') md=loadmodel('JksControl'); @@ -82,3 +85,3 @@ 'title','Friction Coefficient',... 'colorbartitle#3','(m)'); -end +end %}}} Index: /issm/trunk-jpl/examples/LcurveAnalysis/runme.m =================================================================== --- /issm/trunk-jpl/examples/LcurveAnalysis/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/LcurveAnalysis/runme.m (revision 25915) @@ -1,4 +1,5 @@ -step=[1]; -if any(step==1) +steps=[1]; + +if any(steps==1) %{{{ % Generate observations md = model; @@ -12,6 +13,7 @@ 'data',md.results.StressbalanceSolution.Vel,'title','"observed velocities"') save model1 md -end -if any(step==2) +end %}}} + +if any(steps==2) %{{{ % Modify rheology, now constant loadmodel('model1.mat'); @@ -28,6 +30,7 @@ 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities') save model2 md -end -if any(step==3) +end %}}} + +if any(steps==3) %{{{ % Perform L-curve analysis for ice rigidity inversion loadmodel('model2.mat'); @@ -42,6 +45,6 @@ md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); - md.verbose = verbose('solution',false,'control',true); - + md.verbose = verbose('solution',false,'control',true); + % Starting L-curve analysis: % @@ -56,11 +59,11 @@ % respective value of Jo and R (R versus Jo). % - min_alpha = 1.e-20; - max_alpha = 1.e-11; - nstep_alpha = 30; - log_step = (log10(max_alpha)-log10(min_alpha))/nstep_alpha; - log_alphas = [log10(min_alpha):log_step:log10(max_alpha)]; - alphas = 10.^log_alphas; - J = zeros(length(alphas),length(md.inversion.cost_functions)+1); + min_alpha = 1.e-20; + max_alpha = 1.e-11; + nstep_alpha = 30; + log_step = (log10(max_alpha)-log10(min_alpha))/nstep_alpha; + log_alphas = [log10(min_alpha):log_step:log10(max_alpha)]; + alphas = 10.^log_alphas; + J = zeros(length(alphas),length(md.inversion.cost_functions)+1); % Loop over the alphas for i=1:length(alphas), @@ -78,6 +81,6 @@ Jo = Jo + J(:,i); % sum of the cost functions (no regularization term). In this example, only 101 end - R = J(:,end-1)./alphas(:); % only the regularization term - + R = J(:,end-1)./alphas(:); % only the regularization term + % Tip: % A rescale in the axes may be useful to visualize the L-curve. @@ -86,5 +89,5 @@ % % Apply a linear transformation on the original axis (Jo, R): - % + % % | 1 alpha | | Jo | | Jo + alpha*R | | J | % | | | | = | | = | | @@ -100,9 +103,10 @@ hoffset=0.1*Jo; text(Jo+hoffset,R+voffset,[repmat('\alpha = ',length(alphas),1) num2str(alphas(:),'%2.0e')],... - 'FontSize',10,'HorizontalAlignment','left','VerticalAlignment','Middle') + 'FontSize',10,'HorizontalAlignment','left','VerticalAlignment','Middle') xlabel('$\mathrm{log}(\mathcal{J}_0$)','Interpreter','latex') ylabel('$\mathrm{log}(\mathcal{R})$','Interpreter','latex') -end -if any(step==4) +end %}}} + +if any(steps==4) %{{{ %invert for ice rigidity loadmodel('model2.mat'); @@ -114,8 +118,8 @@ md.inversion.maxsteps = maxsteps; md.inversion.cost_functions = [101 502]; - md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1); - md.inversion.cost_functions_coefficients(:,2) = 4.e-17*ones(md.mesh.numberofvertices,1); % here you can use the best value found for alpha - md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); - md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); + md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1); + md.inversion.cost_functions_coefficients(:,2) = 4.e-17*ones(md.mesh.numberofvertices,1); % here you can use the best value found for alpha + md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1); + md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1); %Go solve! @@ -124,3 +128,3 @@ plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',... 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities') -end +end %}}} Index: /issm/trunk-jpl/examples/Pig/runme.m =================================================================== --- /issm/trunk-jpl/examples/Pig/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/Pig/runme.m (revision 25915) @@ -1,12 +1,12 @@ steps=[1]; -if any(steps==1) %Mesh Generation #1 +if any(steps==1) %Mesh Generation #1 %{{{ %Mesh parameters domain =['./DomainOutline.exp']; - hinit=10000; % element size for the initial mesh - hmax=40000; % maximum element size of the final mesh - hmin=5000; % minimum element size of the final mesh - gradation=1.7; % maximum size ratio between two neighboring elements - err=8; % maximum error between interpolated and control field + hinit=10000; % element size for the initial mesh + hmax=40000; % maximum element size of the final mesh + hmin=5000; % minimum element size of the final mesh + gradation=1.7; % maximum size ratio between two neighboring elements + err=8; % maximum error between interpolated and control field % Generate an initial uniform mesh (resolution = hinit m) @@ -14,27 +14,27 @@ % Load Velocities - nsidc_vel='../Data/Antarctica_ice_velocity.nc'; + nsidc_vel='../Data/Antarctica_ice_velocity.nc'; % Get necessary data to build up the velocity grid - xmin = ncreadatt(nsidc_vel,'/','xmin'); - ymax = ncreadatt(nsidc_vel,'/','ymax'); - spacing = ncreadatt(nsidc_vel,'/','spacing'); - nx = double(ncreadatt(nsidc_vel,'/','nx')); - ny = double(ncreadatt(nsidc_vel,'/','ny')); - vx = double(ncread(nsidc_vel,'vx')); - vy = double(ncread(nsidc_vel,'vy')); + xmin = ncreadatt(nsidc_vel,'/','xmin'); + ymax = ncreadatt(nsidc_vel,'/','ymax'); + spacing = ncreadatt(nsidc_vel,'/','spacing'); + nx = double(ncreadatt(nsidc_vel,'/','nx')); + ny = double(ncreadatt(nsidc_vel,'/','ny')); + vx = double(ncread(nsidc_vel,'vx')); + vy = double(ncread(nsidc_vel,'vy')); % Read coordinates - xmin = strtrim(xmin); - xmin = str2num(xmin(1:end-2)); - ymax = strtrim(ymax); - ymax = str2num(ymax(1:end-2)); + xmin = strtrim(xmin); + xmin = str2num(xmin(1:end-2)); + ymax = strtrim(ymax); + ymax = str2num(ymax(1:end-2)); spacing = strtrim(spacing); - spacing = str2num(spacing(1:end-2)); + spacing = str2num(spacing(1:end-2)); % Build the coordinates x=xmin+(0:1:nx)'*spacing; y=(ymax-ny*spacing)+(0:1:ny)'*spacing; - + % Interpolate velocities onto coarse mesh vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0); @@ -45,5 +45,5 @@ % Adapt the mesh to minimize error in velocity interpolation md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err); - + %ploting plotmodel(md,'data','mesh') @@ -51,19 +51,18 @@ % Save model save ./Models/PIG_Mesh_generation md; -end +end %}}} -if any(steps==2) %Masks #2 +if any(steps==2) %Masks #2 %{{{ + md = loadmodel('./Models/PIG_Mesh_generation'); - md = loadmodel('./Models/PIG_Mesh_generation'); - - % Load SeaRISe dataset for Antarctica + % Load SeaRISe dataset for Antarctica % http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica searise='../Data/Antarctica_5km_withshelves_v0.75.nc'; - + %read thickness mask from SeaRISE x1=double(ncread(searise,'x1')); y1=double(ncread(searise,'y1')); thkmask=double(ncread(searise,'thkmask')); - + %interpolate onto our mesh vertices groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0)); @@ -77,11 +76,10 @@ %ploting plotmodel(md,'data',md.mask.ocean_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice') - + % Save model save ./Models/PIG_SetMask md; -end +end %}}} -if any(steps==3) %Parameterization #3 - +if any(steps==3) %Parameterization #3 %{{{ md = loadmodel('./Models/PIG_SetMask'); md = parameterize(md,'./Pig.par'); @@ -89,11 +87,10 @@ % Use a MacAyeal flow model md = setflowequation(md,'SSA','all'); - + % Save model save ./Models/PIG_Parameterization md; -end +end %}}} -if any(steps==4) %Control Method #4 - +if any(steps==4) %Control Method #4 %{{{ md = loadmodel('./Models/PIG_Parameterization'); @@ -135,8 +132,7 @@ % Save model save ./Models/PIG_Control_drag md; -end +end %}}} -if any(steps==5) %Plot #5 - +if any(steps==5) %Plot #5 %{{{ md = loadmodel('./Models/PIG_Control_drag'); @@ -149,8 +145,7 @@ 'caxis#1-2',([1.5,4000]),... 'colorbartitle#3','(m)', 'log#1-2',10); -end +end %}}} -if any(steps==6) %Higher-Order #6 - +if any(steps==6) %Higher-Order #6 %{{{ % Load Model @@ -165,8 +160,7 @@ % Save Model -end +end % step 6 end -if any(steps==7) %Plot #7 - +if any(steps==7) %Plot #7 %{{{ mdHO = loadmodel('./Models/PIG_ModelHO'); mdSSA = loadmodel('./Models/PIG_Control_drag'); @@ -176,13 +170,13 @@ plotmodel(mdHO,'nlines',3,'ncols',2,'axis#all','equal',... - 'data',mdHO.initialization.vel,'title','Observed velocity',... - 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',... - 'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',... - 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',... - 'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',... - 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',... - 'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),... - 'colorbar#all','on','view#all',2,... - 'colorbartitle#all','(m/yr)',... - 'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10); -end + 'data',mdHO.initialization.vel,'title','Observed velocity',... + 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',... + 'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',... + 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',... + 'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',... + 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',... + 'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),... + 'colorbar#all','on','view#all',2,... + 'colorbartitle#all','(m/yr)',... + 'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10); +end %}}} Index: /issm/trunk-jpl/examples/Pig2/runme.m =================================================================== --- /issm/trunk-jpl/examples/Pig2/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/Pig2/runme.m (revision 25915) @@ -1,13 +1,13 @@ -steps=1; - -if any(steps==1) %Mesh Generation #1 +steps=[1]; + +if any(steps==1) %Mesh Generation #1 %{{{ %Mesh parameters domain =['./DomainOutline.exp']; - hinit=5000; % element size for the initial mesh - hmax=40000; % maximum element size of the final mesh - hmin=4000; % minimum element size of the final mesh - gradation=1.7; % maximum size ratio between two neighboring elements - err=8; % maximum error between interpolated and control field + hinit=5000; % element size for the initial mesh + hmax=40000; % maximum element size of the final mesh + hmin=4000; % minimum element size of the final mesh + gradation=1.7; % maximum size ratio between two neighboring elements + err=8; % maximum error between interpolated and control field % Generate an initial uniform mesh (resolution = hinit m) @@ -15,17 +15,17 @@ % Get necessary data to build up the velocity grid - nsidc_vel='../Data/Antarctica_ice_velocity.nc'; - xmin = strsplit(ncreadatt(nsidc_vel,'/','xmin')); xmin = str2num(xmin{2}); - ymax = strsplit(ncreadatt(nsidc_vel,'/','ymax')); ymax = str2num(ymax{2}); - spacing = strsplit(ncreadatt(nsidc_vel,'/','spacing')); spacing = str2num(spacing{2}); - nx = double(ncreadatt(nsidc_vel,'/','nx')); - ny = double(ncreadatt(nsidc_vel,'/','ny')); - vx = double(ncread(nsidc_vel,'vx')); - vy = double(ncread(nsidc_vel,'vy')); + nsidc_vel ='../Data/Antarctica_ice_velocity.nc'; + xmin = strsplit(ncreadatt(nsidc_vel,'/','xmin')); xmin = str2num(xmin{2}); + ymax = strsplit(ncreadatt(nsidc_vel,'/','ymax')); ymax = str2num(ymax{2}); + spacing = strsplit(ncreadatt(nsidc_vel,'/','spacing')); spacing = str2num(spacing{2}); + nx = double(ncreadatt(nsidc_vel,'/','nx')); + ny = double(ncreadatt(nsidc_vel,'/','ny')); + vx = double(ncread(nsidc_vel,'vx')); + vy = double(ncread(nsidc_vel,'vy')); % Build the coordinates x=xmin+(0:1:nx)'*spacing; y=(ymax-ny*spacing)+(0:1:ny)'*spacing; - + % Interpolate velocities onto coarse mesh vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0); @@ -40,18 +40,17 @@ plotmodel(md,'data','mesh') save ./Models/PIG_Mesh_generation md; -end - -if any(steps==2) %Masks #2 - - md = loadmodel('./Models/PIG_Mesh_generation'); +end %}}} + +if any(steps==2) %Masks #2 %{{{ + md = loadmodel('./Models/PIG_Mesh_generation'); % Load SeaRISe dataset for Antarctica http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica searise='../Data/Antarctica_5km_withshelves_v0.75.nc'; - + %read thickness mask from SeaRISE x1=double(ncread(searise,'x1')); y1=double(ncread(searise,'y1')); thkmask=double(ncread(searise,'thkmask')); - + %interpolate onto our mesh vertices groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0)); @@ -64,19 +63,17 @@ plotmodel(md,'data',md.mask.ocean_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice') - + save ./Models/PIG_SetMask md; -end - -if any(steps==3) %Parameterization #3 - +end %}}} + +if any(steps==3) %Parameterization #3 %{{{ md = loadmodel('./Models/PIG_SetMask'); md = setflowequation(md,'SSA','all'); md = parameterize(md,'./Pig.par'); - + save ./Models/PIG_Parameterization md; -end - -if any(steps==4) %Rheology B inversion - +end %}}} + +if any(steps==4) %Rheology B inversion %{{{ md = loadmodel('./Models/PIG_Parameterization'); @@ -120,8 +117,7 @@ % Save model save ./Models/PIG_Control_B md; -end - -if any(steps==5) %drag inversion - +end %}}} + +if any(steps==5) %drag inversion %{{{ md = loadmodel('./Models/PIG_Control_B'); @@ -153,8 +149,9 @@ 'caxis#1-2',([1.5,4000]),... 'colorbartitle#3','(m)', 'log#1-2',10); + save ./Models/PIG_Control_drag md; -end - -if any(steps==6) %Transient Run #1 +end %}}} + +if any(steps==6) %Transient Run #1 %{{{ md = loadmodel('./Models/PIG_Control_drag'); @@ -183,5 +180,5 @@ % Save model save ./Models/PIG_Transient md; -end +end %}}} if any(steps==7) %High Melt #2 @@ -196,7 +193,7 @@ save ./Models/PIG_HighMelt md; -end - -if any(steps==8) %Ice Front retreat +end %}}} + +if any(steps==8) %Ice Front retreat %{{{ md = loadmodel('./Models/PIG_Transient'); @@ -209,7 +206,7 @@ save ./Models/PIG_FrontRetreat md; -end - -if any(steps==9) %High surface mass balance #3 +end %}}} + +if any(steps==9) %High surface mass balance #3 %{{{ %Load model from PIG_Transient %... @@ -229,3 +226,3 @@ %Save model save ./Models/PIG_SMB md; -end +end %}}} Index: /issm/trunk-jpl/examples/PigSensitivity/runme.m =================================================================== --- /issm/trunk-jpl/examples/PigSensitivity/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/PigSensitivity/runme.m (revision 25915) @@ -1,7 +1,10 @@ -steps=[1]; +steps=[1:4]; -if any(steps==1) %Transient Run #1 +try - md = loadmodel('../Pig/Models/PIG_Control_drag'); + +if any(steps==1) %Transient Run #1 %{{{ + + md = loadmodel('../Pig/Models/PIG_Control_drag'); md.inversion.iscontrol=0; @@ -11,9 +14,9 @@ md.transient.ismovingfront=0; md.transient.isthermal=0; - + pos=find(md.mask.ocean_levelset<0); md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); md.basalforcings.floatingice_melting_rate=25*ones(md.mesh.numberofvertices,1); - + md.timestepping.time_step=0.1; md.timestepping.final_time=10; @@ -24,12 +27,12 @@ % Save model save ./Models/PIG_Transient md; -end +end %}}} -if any(steps==2) %High Melt #2 - md = loadmodel('./Models/PIG_Transient'); +if any(steps==2) %High Melt #2 %{{{ + md = loadmodel('./Models/PIG_Transient'); md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); md.basalforcings.floatingice_melting_rate=60*ones(md.mesh.numberofvertices,1); - + md.timestepping.time_step=0.1; md.timestepping.final_time=10; @@ -39,8 +42,8 @@ save ./Models/PIG_HighMelt md; -end +end %}}} -if any(steps==3) %Ice Front retreat - md = loadmodel('./Models/PIG_Transient'); +if any(steps==3) %Ice Front retreat %{{{ + md = loadmodel('./Models/PIG_Transient'); md2=extract(md,'~FrontRetreat.exp'); @@ -58,11 +61,12 @@ save ./Models/PIG_FrontRetreat md2; -end +end %}}} -if any(steps==4) %High surface mass balance #3 +if any(steps==4) + disp('Needs work!') %High surface mass balance #3 %{{{ %Load model %Change external forcing basal melting rate and surface mass balance) - + %Refine time steps and time span of the simulation @@ -73,3 +77,3 @@ %Save model -end +end %}}} Index: /issm/trunk-jpl/examples/SlrFarrell/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrFarrell/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/SlrFarrell/runme.m (revision 25915) @@ -1,38 +1,38 @@ - clear all; -steps=[1]; % [1:5] +steps=[1:5]; -if any(steps==1) +try + % [1:5] + +if any(steps==1) %{{{ disp(' Step 1: Global mesh creation'); numrefine=1; - resolution=150*1e3; % inital resolution [m] + resolution=150*1e3; % inital resolution [m] radius = 6.371012*10^6; % mean radius of Earth, m - mindistance_coast=150*1e3; % coastal resolution [m] + mindistance_coast=150*1e3; % coastal resolution [m] mindistance_land=300*1e3; % resolution on the continents [m] - maxdistance=600*1e3; % max element size (on mid-oceans) [m] + maxdistance=600*1e3; % max element size (on mid-oceans) [m] - md=model; - md.mask=maskpsl(); - md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + md=model; + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] for i=1:numrefine, - - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); distance=zeros(md.mesh.numberofvertices,1); - pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); - pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); for j=1:md.mesh.numberofvertices - phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; + phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; if md.mask.ocean_levelset(j), - phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; + phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); else - phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; + phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); @@ -41,51 +41,44 @@ end pos=find(distancemindistance_land); - distance(pos2)=mindistance_land; - dist=min(maxdistance,distance); + pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); + distance(pos2)=mindistance_land; + + dist=min(maxdistance,distance); md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); - end - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); save ./Models/SlrFarrell_Mesh md; - + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); +end %}}} -end - -if any(steps==2) +if any(steps==2) %{{{ disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); md = loadmodel('./Models/SlrFarrell_Mesh'); - md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere - + md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere md.slr.deltathickness=zeros(md.mesh.numberofelements,1); md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); - save ./Models/SlrFarrell_Loads md; - + save ./Models/SlrFarrell_Loads md; + plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]'); +end %}}} -end - -if any(steps==3) +if any(steps==3) %{{{ disp(' Step 3: Parameterization'); md = loadmodel('./Models/SlrFarrell_Loads'); - nlove=10001; - md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; - md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; - md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; + md.solidearth.lovenumbers=lovenumbers('maxdeg',10000); - md.mask.land_levelset = 1-md.mask.ocean_levelset; - md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); - md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1); + md.mask.land_levelset = 1-md.mask.ocean_levelset; + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1); pos=find(md.mesh.lat <-80); - md.mask.ice_levelset(pos(1))=-1; % ice yes! - md.mask.ocean_levelset(pos(1))=1; % ice grounded! + md.mask.ice_levelset(pos(1))=-1; % ice yes! + md.mask.ocean_levelset(pos(1))=1; % ice grounded! di=md.materials.rho_ice/md.materials.rho_water; @@ -94,67 +87,63 @@ md.geometry.base=md.geometry.surface-md.geometry.thickness; md.geometry.bed=md.geometry.base; - + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - + md.miscellaneous.name='SlrFarrell'; - - save ./Models/SlrFarrell_Parameterization md; -end + save ./Models/SlrFarrell_Parameterization md; +end %}}} -if any(steps==4) +if any(steps==4) %{{{ disp(' Step 4: Solve Slr solver'); md = loadmodel('./Models/SlrFarrell_Parameterization'); - md.cluster=generic('name',oshostname(),'np',3); + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); - md.slr.reltol = 0.1/100; % per cent change in solution + md.slr.reltol = 0.1/100; % percent change in solution md=solve(md,'Slr'); - save ./Models/SlrFarrell_Solution md; + save ./Models/SlrFarrell_Solution md; +end %}}} -end - -if any(steps==5) +if any(steps==5) %{{{ disp(' Step 5: Plot solutions'); md = loadmodel('./Models/SlrFarrell_Solution'); - sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) + sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) res = 1; % [degree] [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); - sol_grid = zeros(size(lat_grid)); + sol_grid = zeros(size(lat_grid)); - F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); - F.Method = 'linear'; - F.ExtrapolationMethod = 'linear'; - + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + sol_grid = F(lat_grid, lon_grid); - sol_grid(isnan(sol_grid))=0; - sol_grid(lat_grid>85 & sol_grid==0) =NaN; + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0)=NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) - figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; [C,h]=contour(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105],'-k','linewidth',2); - clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500); - geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); - plot(long,lat,'k'); hold off; + clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500); + geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); + plot(long,lat,'k'); hold off; c1=colorbar; - colormap(flipud(haxby)); - caxis([96 105]); - xlim([-170 170]); - ylim([-85 85]); - grid on; - title('Relative sea-level [% of GMSL]'); + colormap(flipud(haxby)); + caxis([96 105]); + xlim([-170 170]); + ylim([-85 85]); + grid on; + title('Relative sea-level [% of GMSL]'); set(gcf,'color','w'); - %export_fig('Fig5.pdf'); - -end - + %export_fig('Fig5.pdf'); +end %}}} Index: /issm/trunk-jpl/examples/SlrGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 25915) @@ -1,39 +1,40 @@ - clear all; addpath('../Data','../Functions'); -steps=[1]; % [1:7] - -if any(steps==1) +steps=[1:7]; + +try + % [1:7] + +if any(steps==1) %{{{ disp(' Step 1: Global mesh creation'); numrefine=1; - resolution=150*1e3; % inital resolution [m] + resolution=150*1e3; % inital resolution [m] radius = 6.371012*10^6; % mean radius of Earth, m - mindistance_coast=150*1e3; % coastal resolution [m] + mindistance_coast=150*1e3; % coastal resolution [m] mindistance_land=300*1e3; % resolution on the continents [m] - maxdistance=600*1e3; % max element size (on mid-oceans) [m] - - md=model; - md.mask=maskpsl(); - md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + maxdistance=600*1e3; % max element size (on mid-oceans) [m] + + md=model; + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] for i=1:numrefine, - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); distance=zeros(md.mesh.numberofvertices,1); - pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); - pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); for j=1:md.mesh.numberofvertices - phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; + phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; if md.mask.ocean_levelset(j), - phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; + phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); else - phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; + phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); @@ -43,51 +44,45 @@ pos=find(distancemindistance_land); - distance(pos2)=mindistance_land; - - dist=min(maxdistance,distance); + pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); + distance(pos2)=mindistance_land; + + dist=min(maxdistance,distance); md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); - - end - - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + end + + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); save ./Models/SlrGRACE_Mesh md; - + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); - -end - -if any(steps==2) +end %}}} + +if any(steps==2) %{{{ disp(' Step 2: Define loads in meters of ice height equivalent'); md = loadmodel('./Models/SlrGRACE_Mesh'); year_month = 2007+15/365; - time_range = [year_month year_month]; - - water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - - md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; - - save ./Models/SlrGRACE_Loads md; - + time_range = [year_month year_month]; + + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + + save ./Models/SlrGRACE_Loads md; + plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]'); - -end - -if any(steps==3) +end %}}} + +if any(steps==3) %{{{ disp(' Step 3: Parameterization'); md = loadmodel('./Models/SlrGRACE_Loads'); - nlove=10001; - md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; - md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; - md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; - - md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); - md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + md.solidearth.lovenumbers = lovenumbers('maxdeg',10000); + + md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); pos=find(md.slr.deltathickness~=0); - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - md.mask.land_levelset = 1-md.mask.ocean_levelset; + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.land_levelset = 1-md.mask.ocean_levelset; md.slr.sealevel=zeros(md.mesh.numberofvertices,1); @@ -100,51 +95,49 @@ md.geometry.base=md.geometry.surface-md.geometry.thickness; md.geometry.bed=md.geometry.base; - + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - + md.miscellaneous.name='SlrGRACE'; - - save ./Models/SlrGRACE_Parameterization md; - -end - -if any(steps==4) + + save ./Models/SlrGRACE_Parameterization md; +end %}}} + +if any(steps==4) %{{{ disp(' Step 4: Solve Slr solver'); md = loadmodel('./Models/SlrGRACE_Parameterization'); - md.slr.rigid=1; - md.slr.elastic=1; + md.slr.rigid=1; + md.slr.elastic=1; md.slr.rotation=1; - md.cluster=generic('name',oshostname(),'np',3); + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); md=solve(md,'Slr'); - save ./Models/SlrGRACE_Solution md; - -end - -if any(steps==5) + save ./Models/SlrGRACE_Solution md; +end %}}} + +if any(steps==5) %{{{ disp(' Step 5: Plot solutions'); md = loadmodel('./Models/SlrGRACE_Solution'); - sol1 = md.slr.deltathickness*100; % [cm] - sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] - - sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; - fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; - - res = 1.0; + sol1 = md.slr.deltathickness*100; % [cm] + sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] + + sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; + fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; + + res = 1.0; [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); - sol_grid = zeros(size(lat_grid)); + sol_grid = zeros(size(lat_grid)); for kk=1:2 sol=eval(sprintf('sol%d',kk)); - - if length(sol)==md.mesh.numberofelements + + if length(sol)==md.mesh.numberofelements for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -152,98 +145,96 @@ end for jj=1:md.mesh.numberofvertices - pos=ceil(find(pp==jj)/3); - temp(jj)=mean(sol(pos)); - end - sol=temp'; - end - - F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + pos=ceil(find(pp==jj)/3); + temp(jj)=mean(sol(pos)); + end + sol=temp'; + end + + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); F.Method = 'linear'; - F.ExtrapolationMethod = 'linear'; + F.ExtrapolationMethod = 'linear'; sol_grid = F(lat_grid, lon_grid); - sol_grid(isnan(sol_grid))=0; - sol_grid(lat_grid>85 & sol_grid==0) = NaN; + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) = NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) - figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; if (kk==1) - geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); else - geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); - end - plot(long,lat,'k'); hold off; + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); + end + plot(long,lat,'k'); hold off; c1=colorbar; - colormap('haxby'); - caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); - xlim([-180 180]); - ylim([-90 90]); - grid on; - title(sol_name(kk)); + colormap('haxby'); + caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); set(gcf,'color','w'); - %export_fig(fig_name{kk}); - end - -end - -if any(steps==6) + %export_fig(fig_name{kk}); + end +end %}}} + +if any(steps==6) %{{{ disp(' Step 6: Transient run'); md = loadmodel('./Models/SlrGRACE_Parameterization'); disp('Projecting loads onto the mesh...'); - time_range = 2007 + [15 350]/365; - water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - - [ne,nt]=size(water_load); - md.slr.deltathickness = zeros(ne+1,nt); - md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; - md.slr.deltathickness(ne+1,:)=[1:nt]; % times - + time_range = 2007 + [15 350]/365; + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + [ne,nt]=size(water_load); + md.slr.deltathickness = zeros(ne+1,nt); + md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + md.slr.deltathickness(ne+1,:)=[1:nt]; % times + md.transient.issmb=0; - md.transient.ismasstransport=0; - md.transient.isstressbalance=0; + md.transient.ismasstransport=0; + md.transient.isstressbalance=0; md.transient.isthermal=0; md.transient.isslr=1; - + md.timestepping.start_time=-1; - md.timestepping.final_time=nt; - md.timestepping.time_step=1; - md.timestepping.interp_forcings=0; + md.timestepping.final_time=nt; + md.timestepping.time_step=1; + md.timestepping.interp_forcings=0; md.settings.output_frequency=1; - md.cluster=generic('name',oshostname(),'np',3); + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); md=solve(md,'tr'); - save ./Models/SlrGRACE_Transient md; - -end - -if any(steps==7) + save ./Models/SlrGRACE_Transient md; +end %}}} + +if any(steps==7) %{{{ disp(' Step 7: Plot transient'); md = loadmodel('./Models/SlrGRACE_Transient'); - time = md.slr.deltathickness(end,:); + time = md.slr.deltathickness(end,:); for tt=1:length(time) gmsl(tt) = md.results.TransientSolution(tt).SealevelEustatic*1000; % GMSL rate mm/yr - sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100; % ice equivalent height [cm/yr] - sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000; % mm/yr - end - sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; - movie_name = {'Movie_dH.avi','Movie_slr.avi'}; - - res = 1.0; + sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100; % ice equivalent height [cm/yr] + sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000; % mm/yr + end + sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; + movie_name = {'Movie_dH.avi','Movie_slr.avi'}; + + res = 1.0; [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); - sol_grid = zeros(size(lat_grid)); + sol_grid = zeros(size(lat_grid)); for kk=1:2 sol=eval(sprintf('sol%d',kk)); - - if length(sol)==md.mesh.numberofelements + + if length(sol)==md.mesh.numberofelements for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -251,55 +242,53 @@ end for jj=1:md.mesh.numberofvertices - pos=ceil(find(pp==jj)/3); - temp2(jj,:)=mean(sol(pos,:)); - end - sol=temp2; - end - - vidObj = VideoWriter(movie_name{kk}); - vidObj.FrameRate=2; % frames per second + pos=ceil(find(pp==jj)/3); + temp2(jj,:)=mean(sol(pos,:)); + end + sol=temp2; + end + + vidObj = VideoWriter(movie_name{kk}); + vidObj.FrameRate=2; % frames per second open(vidObj); for jj=1:length(time) - F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); F.Method = 'linear'; - F.ExtrapolationMethod = 'linear'; + F.ExtrapolationMethod = 'linear'; sol_grid = F(lat_grid, lon_grid); - sol_grid(isnan(sol_grid))=0; - sol_grid(lat_grid>85 & sol_grid==0) = NaN; + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) = NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) - figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; if (kk==1) - geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); else - geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); - end - plot(long,lat,'k'); hold off; + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); + end + plot(long,lat,'k'); hold off; c1=colorbar; - colormap('haxby'); - caxis([-floor(min(abs(min(min(sol))),abs(max(max(sol))))) floor(min(abs(min(min(sol))),abs(max(max(sol)))))]); - xlim([-180 180]); - ylim([-90 90]); - grid on; - title(sol_name(kk)); + colormap('haxby'); + caxis([-floor(min(abs(min(min(sol))),abs(max(max(sol))))) floor(min(abs(min(min(sol))),abs(max(max(sol)))))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); set(gcf,'color','w'); - writeVideo(vidObj,getframe(gcf)); - close + writeVideo(vidObj,getframe(gcf)); + close end close(vidObj); end - !open *.avi; - - % plot GMSL time series - plot(time,gmsl,'-*','linewidth',3); grid on; - xlabel('# month'); + !open *.avi; + + % plot GMSL time series + plot(time,gmsl,'-*','linewidth',3); grid on; + xlabel('# month'); ylabel('GMSL [mm/yr]'); set(gcf,'color','w'); %export_fig('Fig7.pdf'); - -end - +end %}}} Index: /issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m (revision 25915) @@ -1,58 +1,57 @@ - clear all; addpath('../Data','../Functions'); -steps=[1]; % [1:8] - -if any(steps==1) +steps=[1:8]; + +try + % [1:8] + +if any(steps==1) %{{{ disp(' Step 1: Global mesh creation'); - resolution=300; % [km] - radius = 6.371012*10^3; % [km] - - md=model; - md.mask=maskpsl(); + resolution=300; % [km] + radius = 6.371012*10^3; % [km] + + md=model; md.mesh=gmshplanet('radius',radius,'resolution',resolution); - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); save ./Models/SlrGRACE_Mesh md; - + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); - -end - -if any(steps==2) +end %}}} + +if any(steps==2) %{{{ disp(' Step 2: Global mesh creation: refined coastlines'); numrefine=1; - resolution=150*1e3; % inital resolution [m] + resolution=150*1e3; % inital resolution [m] radius = 6.371012*10^6; % mean radius of Earth, m - mindistance_coast=150*1e3; % coastal resolution [m] - mindistance_source=75*1e3; % source resolution [m] + mindistance_coast=150*1e3; % coastal resolution [m] + mindistance_source=75*1e3; % source resolution [m] mindistance_land=300*1e3; % resolution on the continents [m] - maxdistance=600*1e3; % max element size (on mid-oceans) [m] - - md=model; - md.mask=maskpsl(); - md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + maxdistance=600*1e3; % max element size (on mid-oceans) [m] + + md=model; + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] for i=1:numrefine, - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); distance=zeros(md.mesh.numberofvertices,1); - pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); - pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); for j=1:md.mesh.numberofvertices - phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; + phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; if md.mask.ocean_levelset(j), - phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; + phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); else - phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; + phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); @@ -62,94 +61,88 @@ pos=find(distancemindistance_land); - distance(pos2)=mindistance_land; + pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); + distance(pos2)=mindistance_land; domain = 'Greenland'; %'Antarctica'; %'HMA'; %'Alaska'; %'Glaciers' - mask = domain_mask(md.mesh.lat,md.mesh.long,domain); + mask = domain_mask(md.mesh.lat,md.mesh.long,domain); distance(mask>0)=mindistance_source; - dist=min(maxdistance,distance); + dist=min(maxdistance,distance); md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); - end - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); save ./Models/SlrGRACE_Mesh_refined md; - + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); - -end - -if any(steps==3) +end %}}} + +if any(steps==3) %{{{ disp(' Step 3: Define loads in meters of ice height equivalent'); md = loadmodel('./Models/SlrGRACE_Mesh_refined'); year_month = 2007+15/365; - time_range = [year_month year_month]; - water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - - index = md.mesh.elements; - lat=90-md.mesh.lat; - lon=md.mesh.long; - lon(lon<0)=180+(180+lon(lon<0)); - - ax_0=lat(index(:,1)); ay_0=lon(index(:,1)); - bx_0=lat(index(:,2)); by_0=lon(index(:,2)); - cx_0=lat(index(:,3)); cy_0=lon(index(:,3)); - + time_range = [year_month year_month]; + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + index = md.mesh.elements; + lat=90-md.mesh.lat; + lon=md.mesh.long; + lon(lon<0)=180+(180+lon(lon<0)); + + ax_0=lat(index(:,1)); ay_0=lon(index(:,1)); + bx_0=lat(index(:,2)); by_0=lon(index(:,2)); + cx_0=lat(index(:,3)); cy_0=lon(index(:,3)); + for ii=1:md.mesh.numberofelements if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180) if ay_0(ii)==0 ay_0(ii)=360; - end + end if by_0(ii)==0 - by_0(ii)=360; - end - if cy_0(ii)==0 - cy_0(ii)=360; - end - end + by_0(ii)=360; + end + if cy_0(ii)==0 + cy_0(ii)=360; + end + end end - - ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; - by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; - cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; - - ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; - by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; - cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; - - lat_element=(ax_0+bx_0+cx_0)/3; + + ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; + by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; + cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; + + ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; + by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; + cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; + + lat_element=(ax_0+bx_0+cx_0)/3; lon_element=(ay_0+by_0+cy_0)/3; - lat_element_0 = 90-lat_element; lon_element_0 = lon_element; + lat_element_0 = 90-lat_element; lon_element_0 = lon_element; lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180; - - domain = 'Greenland'; %'HMA'; %'Alaska'; %'Antarctica'; %'Glaciers'; - mask = domain_mask(lat_element_0,lon_element_0,domain); - - md.slr.deltathickness = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice; - - save ./Models/SlrGRACE_Loads md; - - plotmodel(md,'data',md.slr.deltathickness,'view',[90 90],'title','Ice height equivalent [m]'); - -end - -if any(steps==4) + + domain = 'Greenland'; %'HMA'; %'Alaska'; %'Antarctica'; %'Glaciers'; + mask = domain_mask(lat_element_0,lon_element_0,domain); + + md.slr.deltathickness = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice; + + save ./Models/SlrGRACE_Loads md; + + plotmodel(md,'data',md.slr.deltathickness,'view',[90 90],'title','Ice height equivalent [m]'); +end %}}} + +if any(steps==4) %{{{ disp(' Step 4: Parameterization'); md = loadmodel('./Models/SlrGRACE_Loads'); - nlove=10001; - md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; - md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; - md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; - - md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); - md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + md.solidearth.lovenumbers = lovenumbers('maxdeg',10000); + + md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); pos=find(md.slr.deltathickness~=0); - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - md.mask.land_levelset = 1-md.mask.ocean_levelset; + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.land_levelset = 1-md.mask.ocean_levelset; md.slr.sealevel=zeros(md.mesh.numberofvertices,1); @@ -161,5 +154,5 @@ md.slr.Ugia=zeros(md.mesh.numberofvertices,1); - md.slr.geodetic=1; + md.slr.geodetic=1; di=md.materials.rho_ice/md.materials.rho_water; @@ -168,51 +161,49 @@ md.geometry.base=md.geometry.surface-md.geometry.thickness; md.geometry.bed=md.geometry.base; - + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - + md.miscellaneous.name='SlrGRACE'; - - save ./Models/SlrGRACE_Parameterization md; - -end - -if any(steps==5) + + save ./Models/SlrGRACE_Parameterization md; +end %}}} + +if any(steps==5) %{{{ disp(' Step 5: Solve Slr solver'); md = loadmodel('./Models/SlrGRACE_Parameterization'); - md.slr.rigid=1; - md.slr.elastic=1; + md.slr.rigid=1; + md.slr.elastic=1; md.slr.rotation=1; - md.cluster=generic('name',oshostname(),'np',3); + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); md=solve(md,'Slr'); - save ./Models/SlrGRACE_Solution md; - -end - -if any(steps==6) + save ./Models/SlrGRACE_Solution md; +end %}}} + +if any(steps==6) %{{{ disp(' Step 6: Plot solutions'); md = loadmodel('./Models/SlrGRACE_Solution'); - sol1 = md.slr.deltathickness*100; % [cm] - sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] - - sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; - fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; - - res = 1.0; + sol1 = md.slr.deltathickness*100; % [cm] + sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] + + sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; + fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; + + res = 1.0; [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); - sol_grid = zeros(size(lat_grid)); + sol_grid = zeros(size(lat_grid)); for kk=1:2 sol=eval(sprintf('sol%d',kk)); - - if length(sol)==md.mesh.numberofelements + + if length(sol)==md.mesh.numberofelements for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -220,97 +211,94 @@ end for jj=1:md.mesh.numberofvertices - pos=ceil(find(pp==jj)/3); - temp(jj)=mean(sol(pos)); - end - sol=temp'; - end - - F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + pos=ceil(find(pp==jj)/3); + temp(jj)=mean(sol(pos)); + end + sol=temp'; + end + + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); F.Method = 'linear'; - F.ExtrapolationMethod = 'linear'; + F.ExtrapolationMethod = 'linear'; sol_grid = F(lat_grid, lon_grid); - sol_grid(isnan(sol_grid))=0; - sol_grid(lat_grid>85 & sol_grid==0) = NaN; + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) = NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) - figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; if (kk==1) - geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); else - geoshow(lat,long,'DisplayType','polygon','FaceColor',[0.75 0.75 0.75]); - end - plot(long,lat,'k'); hold off; + geoshow(lat,long,'DisplayType','polygon','FaceColor',[0.75 0.75 0.75]); + end + plot(long,lat,'k'); hold off; c1=colorbar; - colormap('haxby'); - caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); - xlim([-180 180]); - ylim([-90 90]); - grid on; - title(sol_name(kk)); + colormap('haxby'); + caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); set(gcf,'color','w'); - %export_fig(fig_name{kk}); + %export_fig(fig_name{kk}); end -end - -if any(steps==7) +end %}}} + +if any(steps==7) %{{{ disp(' Step 7: Transient run'); md = loadmodel('./Models/SlrGRACE_Parameterization'); disp('Projecting loads onto the mesh...'); - time_range = [2005+15/365 2014+350/365]; - %time_range = 2007 + [15 350]/365; - water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - - [ne,nt]=size(water_load); - md.slr.deltathickness = zeros(ne+1,nt); - md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; - md.slr.deltathickness(ne+1,:)=[1:nt]; % times - + time_range = [2005+15/365 2014+350/365]; + %time_range = 2007 + [15 350]/365; + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + [ne,nt]=size(water_load); + md.slr.deltathickness = zeros(ne+1,nt); + md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + md.slr.deltathickness(ne+1,:)=[1:nt]; % times + md.transient.issmb=0; - md.transient.ismasstransport=0; - md.transient.isstressbalance=0; + md.transient.ismasstransport=0; + md.transient.isstressbalance=0; md.transient.isthermal=0; md.transient.isslr=1; - + md.timestepping.start_time=-1; - md.timestepping.final_time=nt; - md.timestepping.time_step=1; - md.timestepping.interp_forcings=0; + md.timestepping.final_time=nt; + md.timestepping.time_step=1; + md.timestepping.interp_forcings=0; md.settings.output_frequency=1; - md.cluster=generic('name',oshostname(),'np',3); + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); md=solve(md,'tr'); - save ./Models/SlrGRACE_Transient md; - -end - -if any(steps==8) + save ./Models/SlrGRACE_Transient md; +end %}}} + +if any(steps==8) %{{{ disp(' Step 8: Plot transient'); md = loadmodel('./Models/SlrGRACE_Transient'); - time = md.slr.deltathickness(end,:); - - tide_x = 37+29/60; % degree N - tide_y = 126+20/60; % degree E + time = md.slr.deltathickness(end,:); + + tide_x = 37+29/60; % degree N + tide_y = 126+20/60; % degree E for tt=1:length(time) - gmsl(tt) = md.results.TransientSolution(tt).SealevelRSLEustatic*1000; - sol = (md.results.TransientSolution(tt+1).Sealevel-md.results.TransientSolution(tt).Sealevel)*1000; + gmsl(tt) = md.results.TransientSolution(tt).SealevelRSLEustatic*1000; + sol = (md.results.TransientSolution(tt+1).Sealevel-md.results.TransientSolution(tt).Sealevel)*1000; slr_incheon(tt) = griddata(md.mesh.lat,md.mesh.long,sol,tide_x,tide_y); end - plot(2005+time/12,gmsl-gmsl(1),'-*','linewidth',3); grid on; hold on; - plot(2005+time/12,slr_incheon-slr_incheon(1),'-*r','linewidth',3); hold off; - xlabel('Year'); + plot(2005+time/12,gmsl-gmsl(1),'-*','linewidth',3); grid on; hold on; + plot(2005+time/12,slr_incheon-slr_incheon(1),'-*r','linewidth',3); hold off; + xlabel('Year'); ylabel('GMSL [mm]'); set(gcf,'color','w'); - legend('GMSL','Incheon'); - -end - + legend('GMSL','Incheon'); +end %}}} Index: /issm/trunk-jpl/examples/UncertaintyQuantification/runme.m =================================================================== --- /issm/trunk-jpl/examples/UncertaintyQuantification/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/UncertaintyQuantification/runme.m (revision 25915) @@ -1,8 +1,8 @@ %PIG Uncertainty Quantification Application -steps=[1]; - -if any(steps==1) - disp(' Step 1: plot flux gates'); - +steps=[1]; + +if any(steps==1) %{{{ + disp(' Step 1: plot flux gates'); + md = loadmodel('../Pig/Models/PIG_Control_drag'); @@ -28,14 +28,14 @@ '8','9','10','11','12','13'},... 'textposition',textpositions); - -end -if any(steps==2) - disp(' Step 2: compute cross overs from CRESIS'); - +end %}}} + +if any(steps==2) %{{{ + disp(' Step 2: compute cross overs from CRESIS'); + md = loadmodel('../Pig/Models/PIG_Control_drag'); %load cross overs: CRESIS McCord Antarctica, 2009 (courtesy of John Paden) load('../Data/CrossOvers2009.mat'); - + %interpolate cross over errors over our mesh vertices DeltaHH=InterpFromMeshToMesh2d(index,x,y,dhh,md.mesh.x,md.mesh.y); @@ -46,21 +46,21 @@ %filter out unrealistic error ranges flags=ContourToNodes(md.mesh.x,md.mesh.y,'ErrorContour.exp',1); - pos=find(~flags); DeltaHH(pos)=0; + pos=find(~flags); DeltaHH(pos)=0; %avoid large unrealistic values pos=find(DeltaHH>1); DeltaHH(pos)=1; pos=find(DeltaHH<-1); DeltaHH(pos)=-1; - + %transform into absolute errors and setup a minimum error everywhere. DeltaHH=abs(DeltaHH); - pos=find(DeltaHH==0); pos2=find(DeltaHH~=0); + pos=find(DeltaHH==0); pos2=find(DeltaHH~=0); DeltaHH(pos)=min(DeltaHH(pos2)); save Models/PIG.CrossOvers DeltaHH - -end -if any(steps==3) - disp(' Step 3: sampling analysis'); - +end %}}} + +if any(steps==3) %{{{ + disp(' Step 3: sampling analysis'); + %load model and cross over errors md = loadmodel('../Pig/Models/PIG_Control_drag'); @@ -68,16 +68,16 @@ %partition the mesh - md.qmu.numberofpartitions=50; - md=partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,... - 'weighting','on'); - md.qmu.partition=md.qmu.partition-1; %switch partition to c-indexing + npart=50; + partition=partitioner(md,'package','chaco','npart',npart,'weighting','on')-1; %make DeltaHH into our 3 sigma deviation DeltaHH=DeltaHH/6; %2 (to transform DeltaHH into a radius) x 3 (for 3 sigma) - DeltaHH_on_partition=AreaAverageOntoPartition(md,DeltaHH); - DeltaHH_on_grids=DeltaHH_on_partition(md.qmu.partition+1); %just to check in case - - md.qmu.variables.thickness=normal_uncertain('scaled_Thickness',1,1); - md.qmu.variables.thickness.stddev=DeltaHH_on_partition; + DeltaHH_on_partition=AreaAverageOntoPartition(md,DeltaHH,partition); + DeltaHH_on_grids=DeltaHH_on_partition(partition+1); %just to check in case + + md.qmu.variables.thickness=normal_uncertain('descriptor','scaled_Thickness',... + 'mean',ones(npart,1),... + 'stddev',DeltaHH_on_partition,... + 'partition',partition); %responses @@ -110,25 +110,29 @@ %mass flux profiles - md.qmu.mass_flux_profiles={'MassFlux1.exp',... - 'MassFlux2.exp',... - 'MassFlux3.exp',... - 'MassFlux4.exp',... - 'MassFlux5.exp',... - 'MassFlux6.exp',... - 'MassFlux7.exp',... - 'MassFlux8.exp',... - 'MassFlux9.exp',... - 'MassFlux10.exp',... - 'MassFlux11.exp',... - 'MassFlux12.exp',... - 'MassFlux13.exp'}; + md.qmu.mass_flux_profiles={... + 'MassFlux1.exp',... + 'MassFlux2.exp',... + 'MassFlux3.exp',... + 'MassFlux4.exp',... + 'MassFlux5.exp',... + 'MassFlux6.exp',... + 'MassFlux7.exp',... + 'MassFlux8.exp',... + 'MassFlux9.exp',... + 'MassFlux10.exp',... + 'MassFlux11.exp',... + 'MassFlux12.exp',... + 'MassFlux13.exp'... + }; md.qmu.mass_flux_profile_directory='../MassFluxes/'; %% sampling analysis - md.qmu.method =dakota_method('nond_samp'); - md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... - 'seed',1234,... - 'samples',30,... - 'sample_type','lhs'); %random or lhs + md.qmu.method =dakota_method('nond_samp'); + md.qmu.method(end) =dmeth_params_set(... + md.qmu.method(end),... + 'seed',1234,... + 'samples',30,... + 'sample_type','lhs'... + ); %random or lhs %% a variety of parameters @@ -143,5 +147,5 @@ %Turn off verbosity md.verbose=verbose(0); - + %Here, we choose to run with 4 processors, 3 for DAKOTA % while one serves as the master @@ -157,9 +161,10 @@ md=solve(md,'Stressbalance','overwrite','y'); - save ./Models/PIG.Sampling md; -end -if any(steps==4) - disp(' Step 4: sensitivity analysis'); - + save ./Models/PIG.Sampling md; +end %}}} + +if any(steps==4) %{{{ + disp(' Step 4: sensitivity analysis'); + %load model md = loadmodel('../Pig/Models/PIG_Control_drag'); @@ -195,17 +200,19 @@ %mass flux profiles - md.qmu.mass_flux_profiles={'MassFlux1.exp',... - 'MassFlux2.exp',... - 'MassFlux3.exp',... - 'MassFlux4.exp',... - 'MassFlux5.exp',... - 'MassFlux6.exp',... - 'MassFlux7.exp',... - 'MassFlux8.exp',... - 'MassFlux9.exp',... - 'MassFlux10.exp',... - 'MassFlux11.exp',... - 'MassFlux12.exp',... - 'MassFlux13.exp'}; + md.qmu.mass_flux_profiles={... + 'MassFlux1.exp',... + 'MassFlux2.exp',... + 'MassFlux3.exp',... + 'MassFlux4.exp',... + 'MassFlux5.exp',... + 'MassFlux6.exp',... + 'MassFlux7.exp',... + 'MassFlux8.exp',... + 'MassFlux9.exp',... + 'MassFlux10.exp',... + 'MassFlux11.exp',... + 'MassFlux12.exp',... + 'MassFlux13.exp'... + }; md.qmu.mass_flux_profile_directory='../MassFluxes/'; @@ -227,5 +234,5 @@ % while one serves as the master md.cluster=generic('name',oshostname,'np',2); - + %Dakota runs in parallel with a master/slave configuration. % At least 2 cpu's are needed to run the UQ @@ -238,21 +245,22 @@ md.verbose=verbose('qmu',true); md=solve(md,'Stressbalance','overwrite','y'); - - save ./Models/PIG.Sensitivity md; -end -if any(steps==5) + + save ./Models/PIG.Sensitivity md; +end %}}} + +if any(steps==5) %{{{ disp(' Step 5: plot partition'); - + %load model md = loadmodel('./Models/PIG.Sampling'); - + plotmodel(md,'data','mesh','partitionedges','on',... 'linewidth',2, 'axis#all','image','unit','km','colorbar','off',... 'title','Partition Edges on ISSM mesh','grid','on'); - -end -if any(steps==6) - disp(' Step 6: plot histogram'); - +end %}}} + +if any(steps==6) %{{{ + disp(' Step 6: plot histogram'); + %load model md = loadmodel('./Models/PIG.Sampling'); @@ -270,8 +278,9 @@ 'none','EdgeColor','red'); -end -if any(steps==7) - disp(' Step 7: plot sensitivity'); - +end %}}} + +if any(steps==7) %{{{ + disp(' Step 7: plot sensitivity'); + %load model md = loadmodel('./Models/PIG.Sensitivity'); @@ -280,5 +289,5 @@ %which profile are we looking at? - index=1; + index=1; %To plot sensitivities @@ -309,4 +318,3 @@ 'colorbartitle#3','If_{B}','unit#all','km','figure',2,... 'title','Importance Factors: H, \alpha, B'); - -end +end %}}} Index: /issm/trunk-jpl/examples/shakti/runme.m =================================================================== --- /issm/trunk-jpl/examples/shakti/runme.m (revision 25914) +++ /issm/trunk-jpl/examples/shakti/runme.m (revision 25915) @@ -1,5 +1,5 @@ steps=[1:3]; -if any(steps==1) +if any(steps==1) %{{{ disp(' Step 1: Mesh'); @@ -8,6 +8,7 @@ save MoulinMesh md -end -if any(steps==2) +end %}}} + +if any(steps==2) %{{{ disp(' Step 2: Parameterization'); md=loadmodel('MoulinMesh'); @@ -21,5 +22,5 @@ % Change hydrology class to Sommers' SHaKTI model md.hydrology=hydrologyshakti(); - + % Define initial water head such that water pressure is 50% of ice overburden pressure md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base; @@ -48,6 +49,7 @@ save MoulinParam md; -end -if any(steps==3); +end %}}} + +if any(steps==3) %{{{ disp(' Step 3: Solve!'); md=loadmodel('MoulinParam'); @@ -80,3 +82,3 @@ save MoulinTransient md -end +end %}}} Index: /issm/trunk-jpl/jenkins/examples_tests.sh =================================================================== --- /issm/trunk-jpl/jenkins/examples_tests.sh (revision 25915) +++ /issm/trunk-jpl/jenkins/examples_tests.sh (revision 25915) @@ -0,0 +1,623 @@ +#!/bin/bash +################################################################################ +# This script runs the examples tests. It is intended to be called from +# jenkins/jenkins.sh and helps organization. +# +# TODO: +# - Figure out how to remove \ and \n\ from multiline string variables while +# preserving formatting when value is printed to file. +################################################################################ + +## Constants +# +RUNME_FILE='runme.m' +RUN_EXAMPLE=0 +STATUS_HANDLING="\ + disp('SUCCESS');\n\ + catch me\n\ + message=getReport(me);\n\ + fprintf('%s',message);\n\ + disp('FAILURE');\n\ + end\n\ + exit\n\ +" + +cd $ISSM_DIR/examples + +for dir in ./* ; do + if [ -d "${dir}" ]; then + # Some of the examples are incomplete (on purpose). As such, we + # will have to populate the missing steps in order to make sure + # that everything is working. + + cd ${dir} + + if [ "${dir}" == "./AMR" ]; then + sed -i '.bak' '1 s|^.*$|try\n\n&|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Data" ]; then + echo 'Directory contains datasets only; no example to run.' + RUN_EXAMPLE=0 + elif [ "${dir}" == "./EsaGRACE" ]; then + sed -i '.bak' 's|steps=\[1\];|steps=\[1:5\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./EsaWahr" ]; then + sed -i '.bak' 's|steps=\[1\];|steps=\[1:7\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Functions" ]; then + echo "Directory contains functions only; no example to run." + RUN_EXAMPLE=0 + elif [ "${dir}" == "./Greenland" ]; then + # STEP_SEVEN #{{{ + STEP_SEVEN="\ + if any(steps==7) %{{{\n\ + disp(' Step 7: Historical Relaxation run');\n\ + md = loadmodel('./Models/Greenland.Control_drag');\n\ + \n\ + load smbbox\n\ + \n\ + %convert mesh x,y into the Box projection\n\ + [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,39,71);\n\ + [xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);\n\ + \n\ + %Interpolate and set surface mass balance\n\ + index = BamgTriangulate(x1(:),y1(:));\n\ + smb_mo = InterpFromMeshToMesh2d(index,x1(:),y1(:),smbmean(:),xi,yi);\n\ + smb = smb_mo*12/1000*md.materials.rho_freshwater/md.materials.rho_ice;\n\ + md.smb.mass_balance = [smb;1 ];\n\ + \n\ + %Set transient options, run for 20 years, saving every 5 timesteps\n\ + md.timestepping.time_step=0.2;\n\ + md.timestepping.final_time=200;\n\ + md.settings.output_frequency=5;\n\ + \n\ + %Additional options\n\ + md.inversion.iscontrol=0;\n\ + md.transient.requested_outputs={'IceVolume','TotalSmb', ...\n\ + 'SmbMassBalance'};\n\ + md.verbose=verbose('solution',true,'module',true);\n\ + \n\ + %Go solve\n\ + md.cluster=generic('name',oshostname,'np',2);\n\ + md=solve(md,'Transient');\n\ + \n\ + save ./Models/Greenland.HistoricTransient_200yr md;\n\ + end %}}}\n\ + " + #}}} + # STEP_EIGHT #{{{ + STEP_EIGHT="\ + if any(steps==8) %{{{\n\ + %Load historic transient model\n\ + md=loadmodel('./Models/Greenland.HistoricTransient_200yr');\n\ + \n\ + %Create Line Plots of relaxation run. Create a figure.\n\ + figure;\n\ + \n\ + %Save surface mass balance, by looping through 200 years (1000 steps)\n\ + %Note, the first output will always contain output from time step 1\n\ + surfmb=[];\n\ + for i=2:201;\n\ + surfmb=[surfmb md.results.TransientSolution(i).SmbMassBalance];\n\ + end\n\ + \n\ + %Plot surface mass balance time series in first subplot\n\ + subplot(3,1,1);\n\ + plot([1:200],mean(surfmb));\n\ + \n\ + %Title this plot Mean surface mass balance\n\ + title('Mean Surface mass balance');\n\ + \n\ + %Save velocity by looping through 200 years\n\ + vel=[];\n\ + for i=2:201;\n\ + vel=[vel md.results.TransientSolution(i).Vel];\n\ + end\n\ + \n\ + %Plot velocity time series in second subplot\n\ + subplot(3,1,2);\n\ + plot([1:200],mean(vel));\n\ + \n\ + %Title this plot Mean Velocity\n\ + title('Mean Velocity');\n\ + \n\ + %Save Ice Volume by looping through 200 years\n\ + volume=[];\n\ + for i=2:201;\n\ + volume=[volume md.results.TransientSolution(i).IceVolume];\n\ + end\n\ + \n\ + %Plot volume time series in third subplot\n\ + subplot(3,1,3);\n\ + plot([1:200],volume);\n\ + \n\ + %Title this plot Mean Velocity and add an x label of years\n\ + title('Ice Volume');\n\ + xlabel('years');\n\ + end %}}}\n\ + " + #}}} + sed -i '.bak' 's|steps=\[1\];|steps=\[1:8\];\n\ntry\n|' $RUNME_FILE + perl -0755 -p -i -e "s|if any\(steps==7\).*% step 7 end|${STEP_SEVEN}|s" $RUNME_FILE + perl -0755 -p -i -e "s|if any\(steps==8\).*% step 8 end|${STEP_EIGHT}|s" $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./IceBridge" ]; then + sed -i '.bak' 's|steps=\[1\];|steps=\[1:5\];\n\ntry\n|' $RUNME_FILE + perl -0755 -p -i -e "s|\n\t%Mesh greenland without.*return;\n||s" $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./IceflowModels" ]; then + sed -i '.bak' '1 s|^.*$|try\n\n&|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Inversion" ]; then + sed -i '.bak' 's|steps=\[1\];|steps=\[1:4\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./ISMIP" ]; then + # TODO: + # - Run test again with ISMIPF configuration (will likely need to + # add conditional after 'RUN_EXAMPLE -eq 1' block) + # + # RUNME #{{{ + RUNME="\ + try\n\ + %which steps to perform; steps are from 1 to 8\n\ + %step 7 is specific to ISMIPA\n\ + %step 8 is specific to ISMIPF\n\ + \n\ + steps=[1:7]; %ISMIPA\n\ + %steps=[1:6,8]; %ISMIPF\n\ + \n\ + % parameter file to be used, choose between IsmipA.par or IsmipF.par\n\ + ParamFile='IsmipA.par';\n\ + %ParamFile='IsmipF.par';\n\ + \n\ + %Run Steps\n\ + \n\ + %Mesh Generation #1\n\ + if any(steps==1) %{{{\n\ + %initialize md as a new model #help model\n\ + %->\n\ + md=model();\n\ + % generate a squaremesh #help squaremesh\n\ + % Side is 80 km long with 20 points\n\ + %->\n\ + if(ParamFile=='IsmipA.par'),\n\ + md=squaremesh(md,80000,80000,20,20);\n\ + elseif(ParamFile=='IsmipF.par'),\n\ + md=squaremesh(md,100000,100000,30,30);\n\ + end\n\ + % plot the given mesh #plotdoc\n\ + %->\n\ + plotmodel(md,'data','mesh')\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.Mesh_generation md;\n\ + end %}}}\n\ + \n\ + %Masks #2\n\ + if any(steps==2) %{{{\n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.Mesh_generation');\n\ + % set the mask #help setmask\n\ + % all MISMIP nodes are grounded\n\ + %->\n\ + md=setmask(md,'','');\n\ + % plot the given mask #md.mask to locate the field\n\ + %->\n\ + plotmodel(md,'data',md.mask.ocean_levelset);\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.SetMask md;\n\ + end %}}}\n\ + \n\ + %Parameterization #3\n\ + if any(steps==3) %{{{\n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.SetMask');\n\ + % parametrize the model # help parameterize\n\ + % you will need to fil-up the parameter file defined by the\n\ + % ParamFile variable\n\ + %->\n\ + md=parameterize(md,ParamFile);\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.Parameterization md;\n\ + end %}}}\n\ + \n\ + %Extrusion #4\n\ + if any(steps==4) %{{{\n\ + \n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.Parameterization');\n\ + % vertically extrude the preceding mesh #help extrude\n\ + % only 5 layers exponent 1\n\ + %->\n\ + md=extrude(md,9,1);\n\ + % plot the 3D geometry #plotdoc\n\ + %->\n\ + plotmodel(md,'data',md.geometry.base)\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.Extrusion md;\n\ + end %}}}\n\ + \n\ + %Set the flow computing method #5\n\ + if any(steps==5) %{{{\n\ + \n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.Extrusion');\n\ + % set the approximation for the flow computation #help setflowequation\n\ + % We will be using the Higher Order Model (HO)\n\ + %->\n\ + md=setflowequation(md,'HO','all');\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.SetFlow md;\n\ + end %}}}\n\ + \n\ + %Set Boundary Conditions #6\n\ + if any(steps==6) %{{{\n\ + \n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.SetFlow');\n\ + % dirichlet boundary condition are known as SPCs\n\ + % ice frozen to the base, no velocity #md.stressbalance\n\ + % SPCs are initialized at NaN one value per vertex\n\ + %->\n\ + md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);\n\ + %->\n\ + md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);\n\ + %->\n\ + md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);\n\ + % extract the nodenumbers at the base #md.mesh.vertexonbase\n\ + %->\n\ + basalnodes=find(md.mesh.vertexonbase);\n\ + % set the sliding to zero on the bed\n\ + %->\n\ + md.stressbalance.spcvx(basalnodes)=0.0;\n\ + %->\n\ + md.stressbalance.spcvy(basalnodes)=0.0;\n\ + % periodic boundaries have to be fixed on the sides\n\ + % create tabs with the side of the domain\n\ + % for x\n\ + % create maxX #help find\n\ + %->\n\ + maxX=find(md.mesh.x==max(md.mesh.x));\n\ + % create minX\n\ + %->\n\ + minX=find(md.mesh.x==min(md.mesh.x));\n\ + % for y, max X and minX should be excluded\n\ + % create maxY\n\ + %->\n\ + maxY=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));\n\ + % create minY\n\ + %->\n\ + minY=find(md.mesh.y==min(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));\n\ + % set the node that should be paired together\n\ + % #md.stressbalance.vertex_pairing\n\ + %->\n\ + md.stressbalance.vertex_pairing=[minX,maxX;minY,maxY];\n\ + if (ParamFile=='IsmipF.par')\n\ + % if we are dealing with IsmipF the solution is in\n\ + % masstransport\n\ + md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing;\n\ + end\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.BoundaryCondition md;\n\ + end %}}}\n\ + \n\ + %Solving #7\n\ + if any(steps==7) %{{{\n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.BoundaryCondition');\n\ + % Set cluster #md.cluster\n\ + % generic parameters #help generic\n\ + % set only the name and number of process\n\ + %->\n\ + md.cluster=generic('name',oshostname(),'np',2);\n\ + % Set which control message you want to see #help verbose\n\ + %->\n\ + md.verbose=verbose('convergence',true);\n\ + % Solve #help solve\n\ + % we are solving a StressBalanc\n\ + %->\n\ + md=solve(md,'Stressbalance');\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.StressBalance md;\n\ + % plot the surface velocities #plotdoc\n\ + %->\n\ + plotmodel(md,'data',md.results.StressbalanceSolution.Vel)\n\ + end %}}}\n\ + \n\ + %Solving #8\n\ + if any(steps==8) %{{{\n\ + % load the preceding step #help loadmodel\n\ + % path is given by the organizer with the name of the given step\n\ + %->\n\ + md = loadmodel('./Models/ISMIP.BoundaryCondition');\n\ + % Set cluster #md.cluster\n\ + % generic parameters #help generic\n\ + % set only the name and number of process\n\ + %->\n\ + md.cluster=generic('name',oshostname(),'np',2);\n\ + % Set which control message you want to see #help verbose\n\ + %->\n\ + md.verbose=verbose('convergence',true);\n\ + % set the transient model to ignore the thermal model\n\ + % #md.transient \n\ + %->\n\ + md.transient.isthermal=0;\n\ + % define the timestepping scheme\n\ + % everything here should be provided in years #md.timestepping\n\ + % give the length of the time_step (4 years)\n\ + %->\n\ + md.timestepping.time_step=4;\n\ + % give final_time (20*4 years time_steps)\n\ + %->\n\ + md.timestepping.final_time=4*20;\n\ + % Solve #help solve\n\ + % we are solving a TransientSolution\n\ + %->\n\ + md=solve(md,'Transient');\n\ + % save the given model\n\ + %->\n\ + save ./Models/ISMIP.Transient md;\n\ + % plot the surface velocities #plotdoc\n\ + %->\n\ + plotmodel(md,'data',md.results.TransientSolution(20).Vel)\n\ + end %}}}\n\ + " + #}}} + # PAR_A #{{{ + PAR_A="\ + %Parameterization for ISMIP A experiment\n\ + \n\ + %Set the Simulation generic name #md.miscellaneous\n\ + %->\n\ + \n\ + %Geometry\n\ + disp(' Constructing Geometry');\n\ + \n\ + %Define the geometry of the simulation #md.geometry\n\ + %surface is [-x*tan(0.5*pi/180)] #md.mesh\n\ + %->\n\ + md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.);\n\ + %base is [surface-1000+500*sin(x*2*pi/L).*sin(y*2*pi/L)]\n\ + %L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)\n\ + %->\n\ + L=max(md.mesh.x)-min(md.mesh.x);\n\ + md.geometry.base=md.geometry.surface-1000.0+500.0*sin(md.mesh.x*2.0*pi/L).*sin(md.mesh.y*2.0*pi/L);\n\ + %thickness is the difference between surface and base #md.geometry\n\ + %->\n\ + md.geometry.thickness=md.geometry.surface-md.geometry.base;\n\ + %plot the geometry to check it out\n\ + %->\n\ + plotmodel(md,'data',md.geometry.thickness);\n\ + \n\ + disp(' Defining friction parameters');\n\ + \n\ + %These parameters will not be used but need to be fixed #md.friction\n\ + %one friciton coefficient per node (md.mesh.numberofvertices,1)\n\ + %->\n\ + md.friction.coefficient=200.0*ones(md.mesh.numberofvertices,1);\n\ + %one friciton exponent (p,q) per element\n\ + %->\n\ + md.friction.p=ones(md.mesh.numberofelements,1);\n\ + %->\n\ + md.friction.q=ones(md.mesh.numberofelements,1);\n\ + \n\ + disp(' Construct ice rheological properties');\n\ + \n\ + %The rheology parameters sit in the material section #md.materials\n\ + %B has one value per vertex\n\ + %->\n\ + md.materials.rheology_B=6.8067e7*ones(md.mesh.numberofvertices,1);\n\ + %n has one value per element\n\ + %->\n\ + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);\n\ + \n\ + disp(' Set boundary conditions');\n\ + \n\ + %Set the default boundary conditions for an ice-sheet \n\ + % #help SetIceSheetBC\n\ + %->\n\ + md=SetIceSheetBC(md);\n\ + " + #}}} + # PAR_F #{{{ + PAR_F="\ + %Parameterization for ISMIP F experiment\n\ + \n\ + %Set the Simulation generic name #md.miscellaneous\n\ + %->\n\ + \n\ + %Geometry\n\ + disp(' Constructing Geometry');\n\ + \n\ + %Define the geometry of the simulation #md.geometry\n\ + %surface is [-x*tan(3.0*pi/180)] #md.mesh\n\ + %->\n\ + md.geometry.surface=md.mesh.x*tan(3.0*pi/180.0);\n\ + %base is [surface-1000+100*exp(-((x-L/2).^2+(y-L/2).^2)/(10000.^2))]\n\ + %L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)\n\ + %->\n\ + L=max(md.mesh.x)-min(md.mesh.x);\n\ + %->\n\ + md.geometry.base=md.geometry.surface-1000.0+100.0*exp(-((md.mesh.x-L/2.0).^2.0+(md.mesh.y-L/2.0).^2.0)/(10000.^2.0));\n\ + %thickness is the difference between surface and base #md.geometry\n\ + %->\n\ + md.geometry.thickness=md.geometry.surface-md.geometry.base;\n\ + %plot the geometry to check it out\n\ + %->\n\ + plotmodel(md,'data',md.geometry.thickness);\n\ + \n\ + disp(' Defining friction parameters');\n\ + \n\ + %These parameters will not be used but need to be fixed #md.friction\n\ + %one friciton coefficient per node (md.mesh.numberofvertices,1)\n\ + %conversion form year to seconds with #md.constants.yts\n\ + %->\n\ + md.friction.coefficient=sqrt(md.constants.yts/(1000*2.140373*10^-7))*ones(md.mesh.numberofvertices,1);\n\ + %one friciton exponent (p,q) per element\n\ + %->\n\ + md.friction.p=ones(md.mesh.numberofelements,1);\n\ + %->\n\ + md.friction.q=zeros(md.mesh.numberofelements,1);\n\ + \n\ + disp(' Construct ice rheological properties');\n\ + \n\ + %The rheology parameters sit in the material section #md.materials\n\ + %B has one value per vertex\n\ + %->\n\ + md.materials.rheology_B=(1/(2.140373*10^-7/md.constants.yts))*ones(md.mesh.numberofvertices,1);\n\ + %n has one value per element\n\ + %->\n\ + md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);\n\ + \n\ + disp(' Set boundary conditions');\n\ + \n\ + %Set the default boundary conditions for an ice-sheet \n\ + % #help SetIceSheetBC\n\ + %->\n\ + md=SetIceSheetBC(md);\n\ + \n\ + disp(' Initializing velocity and pressure');\n\ + \n\ + %initialize the velocity and pressurefields of #md.initialization\n\ + %->\n\ + md.initialization.vx=zeros(md.mesh.numberofvertices,1);\n\ + %->\n\ + md.initialization.vy=zeros(md.mesh.numberofvertices,1);\n\ + %->\n\ + md.initialization.vz=zeros(md.mesh.numberofvertices,1);\n\ + %->\n\ + md.initialization.pressure=zeros(md.mesh.numberofvertices,1);\n\ + " + #}}} + perl -0755 -p -i'.bak' -e "s|^.*$|${RUNME}|s" $RUNME_FILE + perl -0755 -p -i'.bak' -e "s|^.*$|${PAR_A}|s" IsmipA.par + perl -0755 -p -i'.bak' -e "s|^.*$|${PAR_F}|s" IsmipF.par + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Jakobshavn" ]; then + sed -i '.bak' 's|steps=\[1\];|steps=\[1:4\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./LcurveAnalysis" ]; then + sed -i '.bak' 's|steps=\[1\];|steps=\[1:4\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Mesh" ]; then + # NOTE: Cannot test exptool region selection without GUI + sed -i '.bak' 's|steps=\[1\];|steps=\[1:7\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Pig" ]; then + # STEP_SIX #{{{ + STEP_SIX="\ + if any(steps==6) %{{{\n\ + % Load Model\n\ + md = loadmodel('./Models/PIG_Control_drag');\n\ + % Disable inversion\n\ + md.inversion.iscontrol=0;\n\ + % Extrude Mesh\n\ + disp(' Extruding mesh');\n\ + number_of_layers=3;\n\ + md=extrude(md,number_of_layers,1);\n\ + % Set Flowequation\n\ + disp(' Using HO Ice Flow Model');\n\ + md=setflowequation(md,'HO','all');\n\ + % Solve\n\ + md=solve(md,'Stressbalance');\n\ + % Save Model\n\ + save ./Models/PIG_ModelHO md;\n\ + end %}}}\n\ + " + #}}} + mv ./DomainOutline.bkp ./DomainOutline.exp > /dev/null 2>&1 + sed -i '.bak' "s|steps=\[1\];|steps=\[1:7\];\ntry\n|" $RUNME_FILE + perl -0755 -p -i -e "s|if any\(steps==6\).*% step 6 end|${STEP_SIX}|s" $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./Pig2" ]; then + STEP_NINE="\n disp('Needs work!'); exit" + sed -i '.bak' 's|steps=\[1\];|steps=\[1:9\];\n\ntry\n|' $RUNME_FILE + perl -0755 -p -i -e "s|if any\(steps==9\).*% step 9 end|${STEP_NINE}|s" $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./PigSensitivity" ]; then + # TODO: Implement step 4 + + # STEP_FOUR # {{{ + STEP_FOUR="\n disp('Needs work!')" + #}}} + sed -i '.bak' 's|steps=\[1\];|steps=\[1:4\];\n\ntry\n|' $RUNME_FILE + sed -i '' "s|if any(steps==4)|&${STEP_FOUR}|" $RUNME_FILE + RUN_EXAMPLE=0 + elif [ "${dir}" == "./shakti" ]; then + sed -i '.bak' 's|steps=\[1:3\];|steps=\[1:3\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./SlrFarrell" ]; then + # TODO: Convert from md.slr + sed -i '.bak' 's|steps=\[1\];|steps=\[1:5\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=0 + elif [ "${dir}" == "./SlrGRACE" ]; then + # TODO: Convert from md.slr + sed -i '.bak' 's|steps=\[1\];|steps=\[1:7\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=0 + elif [ "${dir}" == "./SlrGRACE_NIMS" ]; then + # TODO: Convert from md.slr + sed -i '.bak' 's|steps=\[1\];|steps=\[1:8\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=0 + elif [ "${dir}" == "./SquareIceShelf" ]; then + sed -i '.bak' '1 s|^.*$|try\n\n&|' $RUNME_FILE + RUN_EXAMPLE=1 + elif [ "${dir}" == "./UncertaintyQuantification" ]; then + #Step 3: sampling analysis + #Index exceeds the number of array elements (1). + # + #Error in AreaAverageOntoPartition (line 50) + # partvector(i)=sum(weightedvector(pos))/sum(md.qmu.vertex_weight(pos)); + # + #Error in runme (line 78) + # DeltaHH_on_partition=AreaAverageOntoPartition(md,DeltaHH,partition); + sed -i '.bak' 's|steps=\[1\];|steps=\[1:7\];\n\ntry\n|' $RUNME_FILE + RUN_EXAMPLE=0 + else + echo "Not implemented yet!" + exit 1 + fi + + if [ $RUN_EXAMPLE -eq 1 ]; then + echo "Testing example: $(basename $dir)" + LOG_RUNME_FILE="matlab_log_$(basename $dir)_examples.log" + echo -e ${STATUS_HANDLING} >> ${RUNME_FILE} + $MATLAB_PATH/bin/matlab -nodisplay -nosplash -r "addpath $ISSM_DIR/src/m/dev; devpath; addpath $ISSM_DIR/nightlylog/; runme" -logfile $ISSM_DIR/nightlylog/$LOG_RUNME_FILE + echo "starting: $(basename $dir)" >> $ISSM_DIR/nightlylog/matlab_log_examples.log + cat $ISSM_DIR/nightlylog/$LOG_RUNME_FILE >> $ISSM_DIR/nightlylog/matlab_log_examples.log + echo "finished: $(basename $dir)" >> $ISSM_DIR/nightlylog/matlab_log_examples.log + mv -f ${RUNME_FILE}.bak ${RUNME_FILE} + else + echo "Skipping: $(basename $dir)" + fi + + # Extra clean up + if [ "${dir}" == "./ISMIP" ]; then + mv -f IsmipA.par.bak IsmipA.par + mv -f IsmipF.par.bak IsmipF.par + fi + + if [ "${dir}" == "./Pig" ]; then + mv -f DomainOutline.exp DomainOutline.bkp + fi + + cd .. + fi +done Index: /issm/trunk-jpl/jenkins/pine_island-mac-examples =================================================================== --- /issm/trunk-jpl/jenkins/pine_island-mac-examples (revision 25915) +++ /issm/trunk-jpl/jenkins/pine_island-mac-examples (revision 25915) @@ -0,0 +1,83 @@ +# NOTE: Same configuration as pine_island-mac-full + +#--------------------# +# ISSM Configuration # +#--------------------# + +MATLAB_PATH="/Applications/MATLAB_R2019b.app" + +ISSM_CONFIG='\ + --prefix=${ISSM_DIR} \ + --disable-static \ + --enable-development \ + --enable-debugging \ + --with-numthreads=4 \ + --with-matlab-dir=${MATLAB_PATH} \ + --with-python-dir=/System/Library/Frameworks/Python.framework/Versions/2.7 \ + --with-python-numpy-dir=/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy \ + --with-fortran-lib="-L/usr/local/Cellar/gcc/10.2.0/lib/gcc/10 -lgfortran" \ + --with-mpi-include=${ISSM_DIR}/externalpackages/petsc/install/include \ + --with-mpi-libflags="-L${ISSM_EXT_SHARED_DIR}/petsc/install/lib -lmpi -lmpicxx -lmpifort" \ + --with-blas-lapack-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \ + --with-metis-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \ + --with-scalapack-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \ + --with-mumps-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \ + --with-hdf5-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \ + --with-petsc-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \ + --with-gsl-dir=${ISSM_EXT_SHARED_DIR}/gsl/install \ + --with-boost-dir=${ISSM_EXT_SHARED_DIR}/boost/install \ + --with-dakota-dir=${ISSM_EXT_SHARED_DIR}/dakota/install \ + --with-triangle-dir=${ISSM_EXT_SHARED_DIR}/triangle/install \ + --with-chaco-dir=${ISSM_EXT_DIR}/chaco/install \ + --with-m1qn3-dir=${ISSM_EXT_DIR}/m1qn3/install \ + --with-semic-dir=${ISSM_EXT_DIR}/semic/install \ +' + +#-------------------# +# External Packages # +#-------------------# + +EXTERNALPACKAGES=" + autotools install.sh + cmake install.sh + petsc install-3.12-mac.sh + gsl install.sh + boost install-1.7-mac.sh + dakota install-6.2-mac.sh + curl install-7-mac.sh + netcdf install-4.7-parallel.sh + proj install-6.2.sh + gdal install-3-python.sh + gshhg install.sh + gmt install-6-mac.sh + gmsh install-4.sh + triangle install-mac.sh + chaco install.sh + m1qn3 install.sh + semic install.sh + shell2junit install.sh +" + +#---------# +# Testing # +#---------# + +# Test suites +MATLAB_TEST=0 +PYTHON_TEST=0 +JAVASCRIPT_TEST=0 +EXAMPLES_TEST=1 + +# Number of CPUs used in ISSM compilation +# +# NOTE: One is usually safer as some packages are very sensitive to parallel +# compilation +# +NUMCPUS_INSTALL=8 + +# Number of CPUs used in the nightly runs +NUMCPUS_RUN=1 + +# Nightly run options +MATLAB_NROPTIONS="" +PYTHON_NROPTIONS="" Index: /issm/trunk-jpl/test/NightlyRun/test418.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test418.m (revision 25914) +++ /issm/trunk-jpl/test/NightlyRun/test418.m (revision 25915) @@ -9,6 +9,5 @@ %partitioning npart=100; -[partition, md]=partitioner(md,'package','chaco','npart',npart); -partition=partition-1; +partition=partitioner(md,'package','chaco','npart',npart)-1; vector=(1:1:md.mesh.numberofvertices)';