Index: sm/trunk-jpl/examples/MeltingSensitivity/DomainOutline.exp =================================================================== --- /issm/trunk-jpl/examples/MeltingSensitivity/DomainOutline.exp (revision 20736) +++ (revision ) @@ -1,101 +1,0 @@ -## Name: -## Icon:0 -# Points Count Value -95 25000.000000 -# X pos Y pos --1712113.0179281300 -349656.0205056490 --1711509.1917759699 -349354.1074295690 --1709240.6586108401 -349354.1074295690 --1707080.1508345299 -348922.0058743060 --1705459.7700022999 -347193.5996532570 --1703083.2114483600 -345249.1426545770 --1701138.7544496800 -344168.8887664210 --1698438.1197292900 -343952.8379887900 --1698006.0181740201 -341900.3556012940 --1693793.0280102200 -340279.9747690600 --1691092.3932898301 -339631.8224361670 --1688391.7585694401 -339415.7716585360 --1686447.3015707601 -340063.9239914290 --1683530.6160727399 -341036.1524907690 --1682558.3875734000 -343304.6856558960 --1682134.7209454600 -344182.1487658350 --1682126.2860181299 -345681.2442098390 --1681072.0444963200 -345932.4393879640 --1680634.4718407800 -346682.5639403050 --1680571.9614614199 -348557.8753211580 --1677446.4424933300 -348370.3441830730 --1675508.6207331200 -349495.5310115840 --1676616.9911885399 -353567.0975933760 --1672991.8601872099 -353799.6175753630 --1668839.1631938200 -353999.1991486380 --1667132.2298986400 -355121.4651541420 --1666069.5534494901 -357496.8595698890 --1662944.0344814099 -358309.4945015920 --1660756.1712037399 -358747.0671571240 --1657922.9836098100 -361682.8074238540 --1656579.2228470000 -364003.8487410000 --1653004.8841628900 -364748.0635758530 --1653442.4568184200 -361622.5446077650 --1651129.5727820301 -355684.0585683980 --1648929.2928810001 -350845.9692010000 --1646260.7698359799 -347802.2804785470 --1644571.8968503401 -341785.6704672340 --1643503.3064999001 -340869.0986596620 --1643190.7546030900 -336680.9032424240 --1642690.6715682000 -335243.1645171040 --1641913.0916588283 -332341.0814454452 --1640683.9739821283 -330436.6112348017 --1637816.0327364956 -328532.1410241582 --1635241.2508976001 -328652.8685499910 --1631765.4107996500 -327017.1790921330 --1625934.5618617306 -322342.6128395668 --1620608.3852626982 -320438.1426289233 --1616101.6204481323 -320438.1426289233 --1611594.8556335662 -320914.2601815842 --1607907.5026034669 -321866.4952869060 --1602991.0318966676 -324723.2006028712 --1599713.3847588014 -327579.9059188365 --1593977.5022675355 -328532.1410241582 --1590699.8551296694 -334721.6692087495 --1589259.8391499999 -338300.0840570000 --1588627.0898328500 -343572.9950327750 --1588056.6000060199 -345432.3563530700 --1586993.9235568701 -347432.6884926460 --1589259.8391499999 -349315.9832070000 --1585493.6744521901 -352996.1122558430 --1585118.6121760199 -354496.3613605250 --1585861.9157163899 -356240.3342513760 --1586368.8197632500 -358747.0671571240 --1587889.0972700799 -360996.6006787950 --1586681.3716600600 -363935.4286441500 --1585806.2263489999 -365748.2296456410 --1586743.8820394201 -368373.6655788350 --1586056.2678664399 -370123.9562009640 --1585493.6744521901 -371124.1222707520 --1586243.7990045301 -373562.0270658610 --1414257.6683934701 -327460.5482541260 --1412220.4351378100 -325972.8942023180 --1404515.9240538401 -296165.1737279960 --1387984.9332285200 -267002.0106682320 --1360067.9373810000 -107884.1934940000 --1423274.3497420000 -81217.9259650000 --1402446.7987670000 -45916.9921070000 --1478343.8065599999 60338.8188030000 --1503760.4789370000 36334.1837800000 --1528509.5619775900 63667.6809196406 --1623055.4918765700 62841.1506965373 --1725734.5897359999 -61372.6193030000 --1719614.6457630000 -69022.5492690000 --1759088.2843859999 -114004.1374660000 --1755416.3180020000 -117064.1094530000 --1761230.2647760001 -126856.0198080000 --1758170.2927900001 -130527.9861920000 --1762148.2563720001 -135423.9413700000 --1764902.2311590000 -131445.9777880000 --1769798.1863370000 -134199.9525750000 --1775000.1387139999 -131445.9777880000 --1794277.9622269999 -153477.7760890000 --1712576.7101950001 -225081.1205660000 --1729939.5891459100 -335652.4567604720 --1712113.0179281300 -349656.0205056490 - Index: sm/trunk-jpl/examples/MeltingSensitivity/Pig.par =================================================================== --- /issm/trunk-jpl/examples/MeltingSensitivity/Pig.par (revision 20736) +++ (revision ) @@ -1,129 +1,0 @@ -% Parameters to change/Try -friction_coefficient = 10; % default [10] -Temp_change = 0; % default [0 K] - -%Name and Coordinate system -md.miscellaneous.name='PIG'; -md.mesh.epsg=3031; - -%NetCdf Loading -disp(' Loading SeaRISE data from NetCDF'); -ncdata='../Data/Antarctica_5km_withshelves_v0.75.nc'; -x1 = ncread(ncdata,'x1'); -y1 = ncread(ncdata,'y1'); -usrf = ncread(ncdata,'usrf')'; -topg = ncread(ncdata,'topg')'; -temp = ncread(ncdata,'presartm')'; -smb = ncread(ncdata,'presprcp')'; -gflux = ncread(ncdata,'bheatflx_fox')'; - -disp(' Loading velocities data from NetCDF'); -nsidc_vel='../Data/Antarctica_ice_velocity.nc'; - -xmin = ncreadatt(nsidc_vel,'/','xmin'); -xmin = strtrim(xmin); %this is a string, and we need to recover the double value -xmin = xmin(1:end-2); %get rid of the unit -xmin = str2num(xmin); %convert to double - -ymax = ncreadatt(nsidc_vel,'/','ymax'); -ymax = strtrim(ymax); -ymax = ymax(1:end-2); -ymax = str2num(ymax); - -nx = ncreadatt(nsidc_vel,'/','nx'); -ny = ncreadatt(nsidc_vel,'/','ny'); - -spacing = ncreadatt(nsidc_vel,'/','spacing'); -spacing = strtrim(spacing); -spacing = spacing(1:end-2); -spacing = str2num(spacing); - -velx = double(ncread(nsidc_vel,'vx')); -vely = double(ncread(nsidc_vel,'vy')); - -x2=xmin+(0:1:nx)'*spacing; -x2=double(x2); - -y2=(ymax-ny*spacing)+(0:1:ny)'*spacing; -y2=double(y2); - -%Geometry -disp(' Interpolating surface and ice base'); -md.geometry.base = InterpFromGridToMesh(x1,y1,topg,md.mesh.x,md.mesh.y,0); -md.geometry.surface = InterpFromGridToMesh(x1,y1,usrf,md.mesh.x,md.mesh.y,0); -clear usrf, topg; - -disp(' Constructing thickness'); -md.geometry.thickness=md.geometry.surface-md.geometry.base; - -%ensure hydrostatic equilibrium on ice shelf: -di=md.materials.rho_ice/md.materials.rho_water; - -%Get the node numbers of floating nodes -pos=find(md.mask.groundedice_levelset<0); - -%apply a flotation criterion on the precedingly defined nodes and -%redefine base and thickness accordingly -md.geometry.thickness(pos)=1/(1-di)*md.geometry.surface(pos); -md.geometry.base(pos)=md.geometry.surface(pos)-md.geometry.thickness(pos); -md.geometry.hydrostatic_ratio=ones(md.mesh.numberofvertices,1); - -%Set min thickness to 1 meter -pos0=find(md.geometry.thickness<=0); -md.geometry.thickness(pos0)=1; -md.geometry.surface=md.geometry.thickness+md.geometry.base; - -%Initialization parameters -disp(' Interpolating temperatures'); -md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15+Temp_change; -clear temp; - -disp(' Set observed velocities') -vx_obs=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0); -vy_obs=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0); -clear velx vely; - -vel_obs=sqrt(vx_obs.^2+vy_obs.^2); -md.initialization.vx=vx_obs; -md.initialization.vy=vy_obs; -md.initialization.vz=zeros(md.mesh.numberofvertices,1); -md.initialization.vel=vel_obs; - -disp(' Set Pressure'); -md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness; - -disp(' Construct ice rheological properties'); -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); -md.materials.rheology_B=paterson(md.initialization.temperature); - -%Forcings -disp(' Interpolating surface mass balance'); -md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0); -md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice; -clear smb; - -disp(' Set geothermal heat flux'); -md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux,md.mesh.x,md.mesh.y,0); -clear gflux; - -%Friction and inversion set up -disp(' Construct basal friction parameters'); -md.friction.coefficient=friction_coefficient*ones(md.mesh.numberofvertices,1); -md.friction.p=ones(md.mesh.numberofelements,1); -md.friction.q=ones(md.mesh.numberofelements,1); - -%no friction applied on floating ice -pos=find(md.mask.groundedice_levelset<0); -md.friction.coefficient(pos)=0; - -md.inversion=m1qn3inversion(); -md.inversion.vx_obs=vx_obs; -md.inversion.vy_obs=vy_obs; -md.inversion.vel_obs=vel_obs; - -disp(' Set boundary conditions'); -md=SetMarineIceSheetBC(md); -md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1); -md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1); -md.thermal.spctemperature = [md.initialization.temperature;1]; %impose observed temperature on surface -md.masstransport.spcthickness = NaN*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/examples/MeltingSensitivity/runme.m =================================================================== --- /issm/trunk-jpl/examples/MeltingSensitivity/runme.m (revision 20736) +++ /issm/trunk-jpl/examples/MeltingSensitivity/runme.m (revision 20737) @@ -1,221 +1,73 @@ -%Which steps to be performed -steps=[1] ; +step=[1]; -%Run Steps +if step==1 %Transient Run #1 -%Mesh Generation #1 -if any(steps==1) + md = loadmodel('./Models/PIG.Control_drag'); - md.miscellaneous.name='PIG.Mesh_generation'; + md.inversion.iscontrol=0; + md.transient.ismasstransport=1; + md.transient.isstressbalance=1; + md.transient.isgroundingline=1; + md.transient.ismovingfront=0; + md.transient.isthermal=0; + + md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); + md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1); + + md.timestepping.time_step=0.1; + md.timestepping.final_time=10; + md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'} - %Mesh parameters - domain =['./DomainOutline.exp']; - hmax=40000; % maximum element size of the final mesh - hmin=5000; % minimum element size of the final mesh - hinit=10000; % element size for the initial mesh - gradation=1.7; % maximum size ratio between two neighboring elements - err=8; % maximum error between interpolated and control field + md=solve(md,TransientSolutionEnum); - % Generate an initial uniform mesh (resolution = hinit m) - md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1); - - %ploting - plotmodel(md,'data','mesh') - - % Load Velocities - % http://nsidc.org/data/nsidc-0484.html - nsidc_vel='../Data/Antarctica_ice_velocity.nc'; - - % Get necessary data to build up the velocity grid - xmin = ncreadatt(nsidc_vel,'/','xmin'); - xmin = strtrim(xmin); % this is a string, and we need to recover the double value - xmin = xmin(1:end-2); % get rid of the unit - xmin = str2num(xmin); % convert to double - - ymax = ncreadatt(nsidc_vel,'/','ymax'); - ymax = strtrim(ymax); - ymax = ymax(1:end-2); - ymax = str2num(ymax); - - nx = ncreadatt(nsidc_vel,'/','nx'); - ny = ncreadatt(nsidc_vel,'/','ny'); - - spacing = ncreadatt(nsidc_vel,'/','spacing'); - spacing = strtrim(spacing); - spacing = spacing(1:end-2); - spacing = str2num(spacing); - - % Get velocities (Note: You can use ncdisp('file') to see an ncdump) - vx = double(ncread(nsidc_vel,'vx')); - vy = double(ncread(nsidc_vel,'vy')); - - x=xmin+(0:1:nx)'*spacing; - x=double(x); - y=(ymax-ny*spacing)+(0:1:ny)'*spacing; - y=double(y); - - % Interpolate velocities onto coarse mesh - vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0); - vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0); - vel_obs=sqrt(vx_obs.^2+vy_obs.^2); - clear vx vy x y; - - % 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') - - % Convert x,y coordinates (Polar stereo) to lat/lon - [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1); - % Save model - save ./Models/PIG.Mesh_generation md; + save ./Models/PIG.Transient md; end -%Masks #2 -if any(steps==2) +if step==2 %High Melt #2 + md = loadmodel('./Models/PIG.Transient'); - md = loadmodel('./Models/PIG.Mesh_generation'); + 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; + md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'} - % 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)); - groundedice(groundedice<=0)=-1; - clear thkmask; + md=solve(md,TransientSolutionEnum); - %fill in the md.mask structure - md.mask.groundedice_levelset=groundedice; %ice is grounded for mask equal one - md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie - - %ploting - plotmodel(md,'data',md.mask.groundedice_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice') - - % Save model - save ./Models/PIG.SetMask md; + save ./Models/PIG.HighMelt md; end -%Parameterization #3 -if any(steps==3) +if step==3 %High surface mass balance #3 + md = loadmodel('./Models/PIG.Transient'); - md = loadmodel('./Models/PIG.SetMask'); - md = parameterize(md,'./Pig.par'); + md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); + md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1); - % Use a MacAyeal flow model - md = setflowequation(md,'SSA','all'); + md.smb.mass_balance=2*md.smb.mass_balance; - % Save model - save ./Models/PIG.Parameterization md; + md.timestepping.time_step=0.1; + md.timestepping.final_time=10; + + md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'} + + md=solve(md,TransientSolutionEnum); + + save ./Models/PIG.HighMelt md; end -%Control Method #4 -if any(steps==4) +if step==4 %Ice Front retreat + md = loadmodel('./Models/PIG.Transient'); - md = loadmodel('./Models/PIG.Parameterization'); + md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); + md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1); - % Control general - md.inversion.iscontrol=1; - md.inversion.maxsteps=20; - md.inversion.maxiter=40; - md.inversion.dxmin=0.1; - md.inversion.gttol=1.0e-4; - md.verbose=verbose('solution',true,'control',true); + md.timestepping.time_step=0.1; + md.timestepping.final_time=10; + md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'} - % Cost functions - md.inversion.cost_functions=[101 103 501]; - md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3); - md.inversion.cost_functions_coefficients(:,1)=1; - md.inversion.cost_functions_coefficients(:,2)=1; - md.inversion.cost_functions_coefficients(:,3)=8e-15; + md=solve(md,TransientSolutionEnum); - % Controls - md.inversion.control_parameters={'FrictionCoefficient'}; - md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1); - md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1); - - % Additional parameters - md.stressbalance.restol=0.01; - md.stressbalance.reltol=0.1; - md.stressbalance.abstol=NaN; - - % Solve - md.toolkits=toolkits; - md.cluster=generic('name',oshostname,'np',2); - md=solve(md,StressbalanceSolutionEnum); - - % Update model friction fields accordingly - md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient; - - plotmodel(md,'data',md.friction.coefficient) - - % Save model - save ./Models/PIG.Control_drag md; + save ./Models/PIG.HighMelt md; end - -%Plot #5 -if any(steps==5) - - md = loadmodel('./Models/PIG.Control_drag'); - - plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',... - 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,... - 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,... - 'FontSize#all',12,... - 'data',md.initialization.vel,'title','Observed velocity',... - 'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',... - 'data',md.geometry.base,'title','Bed elevation',... - 'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',... - 'colorbar#all','on','colorbartitle#1-2','[m/yr]',... - 'caxis#1-2',([1.5,4000]),... - 'colorbartitle#3','[m]', 'log#1-2',10); -end - -%HO #6 -if any(steps==6) - - % Load Model - - % Disable inversion - - % Extrude Mesh - - % Set Flowequation - - % Solve - - % Save Model - -end - -%Plot #7 -if any(steps==7) - - mdHO = loadmodel('./Models/PIG.ModelHO'); - mdSSA = loadmodel('./Models/PIG.Control_drag'); - - basal=find(mdHO.mesh.vertexonbase); - surf=find(mdHO.mesh.vertexonsurface); - - plotmodel(mdHO,'nlines',3,'ncols',2,'unit#all','km','axis#all','equal',... - 'xlim#all',[min(mdHO.mesh.x) max(mdHO.mesh.x)]/10^3,... - 'ylim#all',[min(mdHO.mesh.y) max(mdHO.mesh.y)]/10^3,... - 'FontSize#all',12,... - '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