Index: /issm/trunk-jpl/examples/Pig2/DomainOutline.exp =================================================================== --- /issm/trunk-jpl/examples/Pig2/DomainOutline.exp (revision 23388) +++ /issm/trunk-jpl/examples/Pig2/DomainOutline.exp (revision 23388) @@ -0,0 +1,101 @@ +## 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: /issm/trunk-jpl/examples/Pig2/FrontRetreat.exp =================================================================== --- /issm/trunk-jpl/examples/Pig2/FrontRetreat.exp (revision 23388) +++ /issm/trunk-jpl/examples/Pig2/FrontRetreat.exp (revision 23388) @@ -0,0 +1,17 @@ +## Name:FrontRetreat +## Icon:0 +# Points Count Value +11 1.000000 +# X pos Y pos +-1626858.6376043619 -327176.4356578548 +-1605190.4936934854 -325911.1425827671 +-1604716.0087903275 -316263.2828852235 +-1607088.4333061168 -312625.5652943464 +-1612465.9288752396 -311676.5954880306 +-1618476.0709819060 -312467.4036599604 +-1624011.7281854146 -313890.8583694341 +-1626858.6376043619 -315788.7979820656 +-1627807.6074106777 -319268.3539385567 +-1628756.5772169936 -322589.7482606619 +-1626858.6376043619 -327176.4356578548 + Index: /issm/trunk-jpl/examples/Pig2/Pig.par =================================================================== --- /issm/trunk-jpl/examples/Pig2/Pig.par (revision 23388) +++ /issm/trunk-jpl/examples/Pig2/Pig.par (revision 23388) @@ -0,0 +1,118 @@ +% 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')'; + +%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); %For Dakota + +%Set min thickness to 1 meter +pos0=find(md.geometry.thickness<=1); +md.geometry.thickness(pos0)=1; +md.geometry.surface=md.geometry.thickness+md.geometry.base; +md.geometry.bed=md.geometry.base; +md.geometry.bed(pos)=md.geometry.base(pos)-1000; + +%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(' Loading velocities data from NetCDF'); +nsidc_vel='../Data/Antarctica_ice_velocity.nc'; +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')); +velx = double(ncread(nsidc_vel,'vx')); +vely = 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)); +spacing = strtrim(spacing); +spacing = str2num(spacing(1:end-2)); +% Build the coordinates +x2=xmin+(0:1:nx)'*spacing; +y2=(ymax-ny*spacing)+(0:1:ny)'*spacing; + +disp(' Set observed velocities') +md.initialization.vx=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0); +md.initialization.vy=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0); +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.vel=sqrt(md.initialization.vx.^2+md.initialization.vy.^2); +clear velx vely; + +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'); +mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0); +md.smb.mass_balance=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.groundingline.migration='SubelementMigration'; + +md.inversion=m1qn3inversion(); +md.inversion.vx_obs=md.initialization.vx; +md.inversion.vy_obs=md.initialization.vy; +md.inversion.vel_obs=md.initialization.vel; + +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; +md.masstransport.spcthickness = NaN*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/examples/Pig2/runme.m =================================================================== --- /issm/trunk-jpl/examples/Pig2/runme.m (revision 23388) +++ /issm/trunk-jpl/examples/Pig2/runme.m (revision 23388) @@ -0,0 +1,231 @@ +steps=[1:8]; + +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 + + % Generate an initial uniform mesh (resolution = hinit m) + md=bamg(model,'domain',domain,'hmax',hinit); + + % 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')); + + % 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); + 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); + + % Plot and save model + plotmodel(md,'data','mesh') + save ./Models/PIG_Mesh_generation md; +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)); + groundedice(groundedice<=0)=-1; + clear thkmask; + + %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 + + plotmodel(md,'data',md.mask.groundedice_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 + + 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 + + md = loadmodel('./Models/PIG_Parameterization'); + + % Control general + md.inversion.iscontrol=1; + md.inversion.maxsteps=40; + md.inversion.maxiter=40; + md.inversion.dxmin=0.1; + md.inversion.gttol=1.0e-6; + md.verbose=verbose('control',true); + + % Cost functions + md.inversion.cost_functions=[101 103 502]; + md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3); + md.inversion.cost_functions_coefficients(:,1)=1000; + md.inversion.cost_functions_coefficients(:,2)=1; + md.inversion.cost_functions_coefficients(:,3)=1.e-16; + + % Controls + md.inversion.control_parameters={'MaterialsRheologyBbar'}; + md.inversion.min_parameters=md.materials.rheology_B; + md.inversion.max_parameters=md.materials.rheology_B; + pos = find(md.mask.groundedice_levelset<0); + md.inversion.min_parameters(pos) = cuffey(273); + md.inversion.max_parameters(pos) = cuffey(200); + + % Additional parameters + md.stressbalance.restol=0.01; + md.stressbalance.reltol=0.1; + md.stressbalance.abstol=NaN; + + % Solve + md.cluster=generic('name',oshostname,'np',2); + mds=extract(md,md.mask.groundedice_levelset<0); + mds=solve(mds,'Stressbalance'); + + % Update model rheology_B accordingly + md.materials.rheology_B(mds.mesh.extractedvertices)=mds.results.StressbalanceSolution.MaterialsRheologyBbar; + plotmodel(md,'data',md.materials.rheology_B) + + % Save model + save ./Models/PIG_Control_B md; +end + +if any(steps==5) %drag inversion + + md = loadmodel('./Models/PIG_Control_B'); + + % 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)=2000; + md.inversion.cost_functions_coefficients(:,2)=1; + md.inversion.cost_functions_coefficients(:,3)=8e-7; + + % 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); + + % Solve + md=solve(md,'Stressbalance'); + + % Update model friction fields accordingly + md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient; + + %Plot and save + plotmodel(md,... + '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); + save ./Models/PIG_Control_drag md; +end + +if any(steps==6) %Transient Run #1 + + md = loadmodel('../Pig/Models/PIG_Control_drag'); + + 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.verbose.solution=1; + md.timestepping.time_step=0.1; + md.timestepping.final_time=10; + md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}; + + %Set melt to 25 m/yr under floating ice + md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); + md.basalforcings.floatingice_melting_rate=25*ones(md.mesh.numberofvertices,1); + + + md=solve(md,'Transient'); + time = cell2mat({md.results.TransientSolution(:).time}); + V_control = cell2mat({md.results.TransientSolution(:).IceVolumeAboveFloatation}); + plot(time,V_control); legend({'Control run'}); + + % Save model + save ./Models/PIG_Transient md; +end + +if any(steps==7) %High Melt #2 + md = loadmodel('./Models/PIG_Transient'); + + %Set melt to 60 m/yr under floating ice + md.basalforcings.floatingice_melting_rate=60*ones(md.mesh.numberofvertices,1); + + md=solve(md,'Transient'); + V_melt = cell2mat({md.results.TransientSolution(:).IceVolumeAboveFloatation}); + plot(time,[V_control',V_melt']); legend({'Control run','High melt'}); + + save ./Models/PIG_HighMelt md; +end + +if any(steps==8) %Ice Front retreat + md = loadmodel('./Models/PIG_Transient'); + + pos = find(ContourToNodes(md.mesh.x,md.mesh.y,'FrontRetreat.exp',2)); + md.mask.ice_levelset(pos) = +1; %Deactivate nodes in retreated area + + md=solve(md,'Transient'); + V_retreat = cell2mat({md.results.TransientSolution(:).IceVolumeAboveFloatation}); + plot(time,[V_control',V_melt',V_retreat']); legend({'Control run','High melt','Retreat'}); + + save ./Models/PIG_FrontRetreat md; +end + +if any(steps==9) %High surface mass balance #3 + %Load model from PIG_Transient + %... + + %Change surface mass balance (x2) + %... + + %Solve + %... + + %Get volume time series + %... + + %plot + plot(time,[V_control',V_melt',V_retreat',V_smb']); legend({'Control run','High melt','Retreat','SMB'}); + + %Save model + save ./Models/PIG_SMB md; +end