Index: /issm/trunk-jpl/examples/SlrFarrell/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrFarrell/runme.m (revision 22805) +++ /issm/trunk-jpl/examples/SlrFarrell/runme.m (revision 22805) @@ -0,0 +1,200 @@ + +clear all; +steps=[2]; % + +if any(steps==1) % Global mesh creation {{{ + disp(' Step 1: Global mesh creation'); + + numrefine=1; + resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. + radius = 6.371012*10^6; % mean radius of Earth, 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] + + %mesh earth: + md=model; + md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. + + for i=1:numrefine, + + %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + + %figure out distance to the coastline, in lat,long (not x,y,z): + 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); + + for j=1:md.mesh.numberofvertices + %figure out nearest coastline (using the great circle distance) + 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; + 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; + 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)); + end + distance(j)=min(d); + end + pos=find(distancemindistance_land); + distance(pos2)=mindistance_land; + + dist=min(maxdistance,distance); % max size 1000 km + %use distance to the coastline to refine mesh: + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); + end + + %figure out mask: + 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'); + %export_fig('Fig1.pdf'); + +end % }}} + +if any(steps==2) % Define source {{{ + disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); + md = loadmodel('./Models/SlrFarrell.Mesh'); + + % initial sea-level: 1 m RSL everywhere. + md.slr.sealevel=md.mask.ocean_levelset; + + save ./Models/SlrFarrell.Loads md; + + plotmodel (md,'data',md.slr.sealevel,'view',[90 90],... + 'title#all','Initial sea-level [m]'); + %export_fig('Fig2.pdf'); + +end % }}} + +if any(steps==3) % Parameterization {{{ + disp(' Step 3: Parameterization'); + md = loadmodel('./Models/SlrFarrell.Loads'); + + % Love numbers and reference frame: CF or CM (choose one!) + nlove=10001; % up to 10,000 degree + md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; + md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; + md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + + % Mask: for computational efficiency only those elements that have loads are convolved! + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + pos=find(md.slr.deltathickness~=0); + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads + md.mask.land_levelset = 1-md.mask.ocean_levelset; + + %% IGNORE BUT DO NOT DELETE %% {{{ + % Geometry: Important only when you want to couple with Ice Flow Model + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.thickness=ones(md.mesh.numberofvertices,1); + md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.surface-md.geometry.thickness; + md.geometry.bed=md.geometry.base; + % Materials: + 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); + % Miscellaneous: + md.miscellaneous.name='SlrFarrell'; + %% IGNORE BUT DO NOT DELETE %% }}} + + save ./Models/SlrFarrell.Parameterization md; + +end % }}} + +if any(steps==4) % Solve {{{ + disp(' Step 4: Solve Slr solver'); + md = loadmodel('./Models/SlrFarrell.Parameterization'); + + % Request outputs + md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; + + % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + + % Solve + md=solve(md,'Slr'); + + save ./Models/SlrFarrell.Solution md; + +end % }}} + +if any(steps==5) % Plot solutions {{{ + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/SlrFarrell.Solution'); + + % loads and solutions. + sol1 = md.slr.deltathickness*100; % WEH cm + sol2 = md.results.SlrSolution.SlrUmotion*1000; % [mm] + sol3 = md.results.SlrSolution.SlrNmotion*1000; % [mm] + sol4 = md.results.SlrSolution.SlrEmotion*1000; % [mm] + sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... + 'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; + + res = 1.0; % degree + + % Make a grid of lats and lons, based on the min and max of the original vectors + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + + for kk=1:4 + sol=eval(sprintf('sol%d',kk)); + + % if data are on elements, map those on to the vertices {{{ + if length(sol)==md.mesh.numberofelements + % map on to the vertices + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); + end + for jj=1:md.mesh.numberofvertices + pos=ceil(find(pp==jj)/3); + temp(jj)=mean(sol(pos)); + end + sol=temp'; + end % }}} + + % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + + % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; + plot(long,lat,'k'); hold off; + c1=colorbar; + colormap(jet); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); + set(gcf,'color','w'); + + %export_fig('Fig5.pdf'); + end + +end % }}} + +