Index: /issm/trunk-jpl/examples/StISSM/DomainOutline.exp =================================================================== --- /issm/trunk-jpl/examples/StISSM/DomainOutline.exp (revision 27531) +++ /issm/trunk-jpl/examples/StISSM/DomainOutline.exp (revision 27531) @@ -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/StISSM/PigStISSM.par =================================================================== --- /issm/trunk-jpl/examples/StISSM/PigStISSM.par (revision 27531) +++ /issm/trunk-jpl/examples/StISSM/PigStISSM.par (revision 27531) @@ -0,0 +1,117 @@ + +% Name and Coordinate system +md.miscellaneous.name = 'PigStISSM'; +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.ocean_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'); +Temp_change = 0; +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'); +friction_coefficient = 10.0; +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.ocean_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.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1); +md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1); +md = SetMarineIceSheetBC(md); +md.thermal.spctemperature = md.initialization.temperature; +md.masstransport.spcthickness = NaN*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/examples/StISSM/runme.m =================================================================== --- /issm/trunk-jpl/examples/StISSM/runme.m (revision 27531) +++ /issm/trunk-jpl/examples/StISSM/runme.m (revision 27531) @@ -0,0 +1,288 @@ + +%%% +% Tutorial for StISSM +%%% + +steps = [5]; + +if any(steps==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_StISSM_Mesh_generation md; + +end + +if any(steps==2) %Masks #2 + md = loadmodel('./Models/PIG_StISSM_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.ocean_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.ocean_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice') + + save ./Models/PIG_StISSM_SetMask md; +end + +if any(steps==3) %Parameterization #3 + + md = loadmodel('./Models/PIG_StISSM_SetMask'); + md = setflowequation(md,'SSA','all'); + md = parameterize(md,'./PigStISSM.par'); + + save ./Models/PIG_StISSM_Parameterization md; +end + +if any(steps==4) %Stochastic SMB + md = loadmodel('./Models/PIG_StISSM_Parameterization'); + + % Create the different subdomains for SMB % + ymax = max(md.mesh.y); + ymin = min(md.mesh.y); + xmax = max(md.mesh.x); + xmin = min(md.mesh.x); + idsmb = zeros(md.mesh.numberofelements,1); %subdomain ID is defined for each element + iid1 = find(md.mesh.x>=(xmax-2/3*(xmax-xmin))); %vertices in subdomain 1 + iid2 = find(md.mesh.x<(xmax-2/3*(xmax-xmin))); %vertices in subdomain 2 + for ii=1:md.mesh.numberofelements + for vertex=1:3 + if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in subdomain 1 + idsmb(ii) = 1; + end + end + if idsmb(ii)==0 %no vertex was found in subdomain 1 + idsmb(ii) = 2; + end + end + + % SMBarma implementation + md.smb = SMBarma(); + md.smb.num_basins = 2; %we use two different subdomains + md.smb.basin_id = idsmb; %element subdomain IDs + md.smb.num_breaks = 1; %1 breakpoint in the piecewise polynomial + md.smb.datebreaks = [5;5]; %breakpoint occurs at year 5 in both subdomains + md.smb.num_params = 2; %use a constant and a linear trend for the piecwise polynomial + constsmb = [0.5,0.2;0.3,0.5]; %constant SMB term for subdomains pre- and post-breakpoint [m yr^-1] + trendsmb = [0,0;0.01,0.001]; %trend is SMB for subdomains pre- and post-breakpoint [m yr^-2] + md.smb.polynomialparams = cat(3,constsmb,trendsmb); %concatenating const and trend along a 3rd dimension + md.smb.ar_order = 1; %first-order AR + md.smb.ma_order = 1; %first-order MA + md.smb.arlag_coefs = [0.3;0]; %AR coefficients in each subdomain + md.smb.malag_coefs = zeros(md.smb.num_basins,md.smb.ma_order); %all zeros is equivalent to MA order 0 + md.smb.arma_timestep = 1.0; %yearly ARMA + md.smb.elevationbins = [300,1000;300,1000]; %elevations separating different lapse rate values [m] + md.smb.lapserates = 1e-2*[0.03,0.01,-1e-4;0.02,0.02,-1e-5]; %lapse rate values [m ice eq. m^-1 yr^-1] + md.smb.refelevation = [500,500]; %elevations at which the SMBarma calculated values apply (i.e., before using lapse rates) + + % Set-up the covariance matrix + sdevSMB1 = 0.01; %low standard deviation in subdomain 1 [m ice eq. yr^-1] + sdevSMB2 = 0.2; %higher variability in subdomain 2 [m ice eq. yr^-1] + correlationSMB = [1.0,0.5;0.5,1.0]; %moderate correlation between the subdomains + covarianceSMB = diag([sdevSMB1,sdevSMB2])*correlationSMB*diag([sdevSMB1,sdevSMB2]); %covariance matrix [(m ice eq. yr^-1)^2] + + % Stochasticforcing implementation + md.stochasticforcing.isstochasticforcing = 1; %activate stochasticity + md.stochasticforcing.fields = [{'SMBarma'}]; + md.stochasticforcing.covariance = covarianceSMB; %prescribe the SMB covariance + md.stochasticforcing.stochastictimestep = 1.0; %yearly stochasticity + + save ./Models/PIG_StISSM_StochSMB md; + +end + +if any(steps==5) % Transient Run #1 + + md = loadmodel('./Models/PIG_StISSM_StochSMB'); + + 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','SmbMassBalance','MaskIceLevelset'}; + + md = solve(md,'Transient'); + timing = cell2mat({md.results.TransientSolution(:).time}); + fullsmb = cell2mat({md.results.TransientSolution(:).SmbMassBalance}); + timemeansmb = mean(fullsmb,2); + timesdevsmb = std(fullsmb,0,2); + maxsmb = max(fullsmb,[],'all'); + minsmb = min(fullsmb,[],'all'); + + + plotmodel(md,'figure',1,... + 'data',md.results.TransientSolution(1).SmbMassBalance,'title','SMB at timestep 1','caxis#1',[minsmb,maxsmb],... + 'data',md.results.TransientSolution(end).SmbMassBalance,'title','SMB at last timestep','caxis#2',[minsmb,maxsmb]) + plotmodel(md,'figure',2,... + 'data',timemeansmb,'title','SMB mean over simulation',... + 'data',timesdevsmb,'title','SMB standard deviation over simulation') + + seriesSMB1 = fullsmb(1,:); + seriesSMB71 = fullsmb(71,:); + figure(3); + plot(timing,seriesSMB1,'b'); + hold on + plot(timing,seriesSMB71,'r'); + title('SMB at two different vertices'); + hold off + + save ./Models/PIG_StISSM_Transient1 md; + +end + +if any(steps==6) % Set up stochastic calving also + + md = loadmodel('./Models/PIG_StISSM_Transient1'); + + md.calving.calvingrate = 20*ones(md.mesh.numberofvertices,1); + md.frontalforcings.meltingrate = zeros(md.mesh.numberofvertices,1); + md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); + md.levelset.migration_max = 100.0; %avoid too fast advance/retreat of the front + md.transient.ismovingfront = 1; + + % Parameterize stochastic calving + idbasinsCalving = md.smb.basin_id; %same subdomains as for SMB + sdevClv1 = 0.5; %low standard deviation in subdomain 1 [m yr^-1] + sdevClv2 = 5; %higher variability in subdomain 2 [m yr^-1] + correlationClv = [1.0,0.0;0.0,1.0]; %no correlation between the subdomains + covarianceClv = diag([sdevClv1,sdevClv2])*correlationClv*diag([sdevClv1,sdevClv2]); %covariance matrix [(m yr^-1)^2] + + % Adjust stochastic forcing class + md.stochasticforcing.fields = [{'SMBarma'},{'DefaultCalving'}]; + oldcovarianceSMB = md.stochasticforcing.covariance; + covarianceGlobal = blkdiag(oldcovarianceSMB,covarianceClv); %independence between SMB and calving + md.stochasticforcing.covariance = covarianceGlobal; + md.stochasticforcing.defaultdimension = md.smb.num_basins; %2 basins for stochastic default (used for calving, same as SMB) + md.stochasticforcing.default_id = idbasinsCalving; + + save ./Models/PIG_StISSM_StochCalving md; + +end + +if any(steps==7) % Transient Run #2 + + md = loadmodel('./Models/PIG_StISSM_StochCalving'); + + md.verbose.solution = 1; + md.timestepping.start_time = md.timestepping.final_time; + md.timestepping.time_step = 0.1; + md.timestepping.final_time = md.timestepping.start_time+5; + md.transient.requested_outputs = {'default','SmbMassBalance','CalvingCalvingrate'}; + + % Set up initial conditions from previous transient results + md.geometry.thickness = md.results.TransientSolution(end).Thickness; + md.initialization.vx = md.results.TransientSolution(end).Vx; + md.initialization.vy = md.results.TransientSolution(end).Vy; + md.initialization.vel = md.results.TransientSolution(end).Vel; + md.mask.ocean_levelset = md.results.TransientSolution(end).MaskOceanLevelset; + md.initialization.pressure = md.results.TransientSolution(end).Pressure; + md.geometry.base = md.results.TransientSolution(end).Base; + md.geometry.surface = md.results.TransientSolution(end).Surface; + md.mask.ice_levelset = md.results.TransientSolution(end).MaskIceLevelset; + + md = solve(md,'Transient'); + md = loadresultsfromcluster(md); + timing = cell2mat({md.results.TransientSolution(:).time}); + fullsmb = cell2mat({md.results.TransientSolution(:).SmbMassBalance}); + fullclv = cell2mat({md.results.TransientSolution(:).CalvingCalvingrate}); + + timemeansmb = mean(fullsmb,2); + timesdevsmb = std(fullsmb,0,2); + maxsmb = max(fullsmb,[],'all'); + minsmb = min(fullsmb,[],'all'); + timemeanclv = mean(fullclv,2); + timesdevclv = std(fullclv,0,2); + maxclv = max(fullclv,[],'all'); + minclv = min(fullclv,[],'all'); + + plotmodel(md,'figure',1,... + 'data',md.results.TransientSolution(1).SmbMassBalance,'title','SMB at timestep 1','caxis#1',[minsmb,maxsmb],... + 'data',md.results.TransientSolution(end).SmbMassBalance,'title','SMB at last timestep','caxis#2',[minsmb,maxsmb]) + plotmodel(md,'figure',2,... + 'data',timemeansmb,'title','SMB mean over simulation',... + 'data',timesdevsmb,'title','SMB standard deviation over simulation') + plotmodel(md,'figure',3,... + 'data',timemeanclv,'title','Calving mean over simulation',... + 'data',timesdevclv,'title','Calving standard deviation over simulation') + + seriesSMB1 = fullsmb(1,:); + seriesSMB71 = fullsmb(71,:); + figure(4); + plot(timing,seriesSMB1,'b'); + hold on + plot(timing,seriesSMB71,'r'); + title('SMB at two different vertices'); + hold off + + seriesclv1 = fullclv(1,:); + seriesclv71 = fullclv(71,:); + figure(5); + plot(timing,seriesclv1,'b'); + hold on + plot(timing,seriesclv71,'r'); + title('Calving at two different vertices'); + hold off + +end + + + + + + + + Index: /issm/trunk-jpl/src/m/classes/SMBarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27530) +++ /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27531) @@ -20,4 +20,5 @@ self.polynomialparams = np.nan self.ar_order = 0 + self.ma_order = 0 self.arlag_coefs = np.nan self.ma_order = 0 @@ -29,4 +30,5 @@ self.elevationbins = np.nan self.refelevation = np.nan + self.datebreaks = np.nan self.steps_per_step = 1 self.averaging = 0 Index: /issm/trunk-jpl/src/m/classes/clusters/pace.m =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/pace.m (revision 27530) +++ /issm/trunk-jpl/src/m/classes/clusters/pace.m (revision 27531) @@ -9,5 +9,5 @@ properties (SetAccess=public) % {{{ - name = 'login-phoenix-4.pace.gatech.edu' %Phoenix cluster name + name = 'login-phoenix-slurm.pace.gatech.edu' %Phoenix cluster name login = ''; %personal login numnodes = 1; %number of nodes requested @@ -17,5 +17,5 @@ queue = 'inferno'; %queue time = 60; %time requested per run [minutes] - accountname = 'GT-arobel3-atlas'; %group account name + accountname = 'gts-arobel3-atlas'; %group account name codepath = ''; %path to issm binaries executionpath = ''; %path for execution folder @@ -62,17 +62,20 @@ fid=fopen([modelname '.queue'],'w'); fprintf(fid,'#!/bin/sh\n'); - fprintf(fid,'#PBS -l walltime=%i\n',cluster.time*60); %conversion of walltime from minutes to seconds. - fprintf(fid,'#PBS -N %s\n',modelname); - fprintf(fid,'#PBS -l nodes=1:ppn=%i\n',cluster.np); - fprintf(fid,'#PBS -l pmem=%igb\n',cluster.mem); - fprintf(fid,'#PBS -q %s\n',cluster.queue); - fprintf(fid,'#PBS -A %s\n',cluster.accountname); - - fprintf(fid,'#PBS -o %s/%s/%s.outlog \n',cluster.executionpath,dirname,modelname); - fprintf(fid,'#PBS -e %s/%s/%s.errlog \n\n',cluster.executionpath,dirname,modelname); - fprintf(fid,'export PBS_O_WORKDIR=%s\n',[cluster.executionpath '/' dirname]); - fprintf(fid,'cd $PBS_O_WORKDIR\n'); - fprintf(fid,'mpiexec -np %i %s/%s %s %s %s \n',cluster.np,cluster.codepath,executable,solution,[cluster.executionpath '/' dirname],modelname); + fprintf(fid,'#SBATCH -t%i\n',cluster.time); + fprintf(fid,'#SBATCH -J%s\n',modelname); + fprintf(fid,'#SBATCH -N 1 --ntasks-per-node=%i\n',cluster.np); + %fprintf(fid,'#SBATCH -N %i\n',cluster.numnodes); + %fprintf(fid,'#SBATCH --ntasks=1\n'); + %fprintf(fid,'#SBATCH --cpus-per-task=%i\n',cluster.np); + fprintf(fid,'#SBATCH --mem-per-cpu=%iG\n',cluster.mem); + fprintf(fid,'#SBATCH -p%s\n',cluster.queue); + fprintf(fid,'#SBATCH -A %s\n',cluster.accountname); + fprintf(fid,'#SBATCH -o%s/%s/%s.outlog \n',cluster.executionpath,dirname,modelname); + fprintf(fid,'#SBATCH -e%s/%s/%s.errlog \n\n',cluster.executionpath,dirname,modelname); + fprintf(fid,'export SLURM_SUBMIT_DIR=%s\n',[cluster.executionpath '/' dirname]); + fprintf(fid,'cd $SLURM_SUBMIT_DIR\n'); + fprintf(fid,'export LD_LIBRARY_PATH=/opt/slurm/current/lib:/opt/pmix/current/lib:$LD_LIBRARY_PATH \n'); + fprintf(fid,'srun --mpi=pmi2 -n %i %s/%s %s %s %s \n',cluster.np,cluster.codepath,executable,solution,[cluster.executionpath '/' dirname],modelname); fclose(fid); @@ -97,8 +100,8 @@ disp('launching solution sequence on remote cluster'); if ~isempty(restart) - launchcommand=['cd ' cluster.executionpath ' && cd ' dirname ' && qsub ' modelname '.queue ']; + launchcommand=['cd ' cluster.executionpath ' && cd ' dirname ' && sbatch ' modelname '.queue ']; else launchcommand=['cd ' cluster.executionpath ' && rm -rf ./' dirname ' && mkdir ' dirname ... - ' && cd ' dirname ' && mv ../' dirname '.tar.gz ./ && tar -zxf ' dirname '.tar.gz && qsub ' modelname '.queue ']; + ' && cd ' dirname ' && mv ../' dirname '.tar.gz ./ && tar -zxf ' dirname '.tar.gz && sbatch ' modelname '.queue ']; end issmssh(cluster.name,cluster.login,cluster.port,launchcommand); Index: /issm/trunk-jpl/src/m/classes/hydrologyarmapw.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyarmapw.m (revision 27530) +++ /issm/trunk-jpl/src/m/classes/hydrologyarmapw.m (revision 27531) @@ -18,4 +18,5 @@ malag_coefs = NaN; datebreaks = NaN; + requested_outputs = {}; end methods @@ -184,5 +185,13 @@ WriteData(fid,prefix,'data',dbreaks,'name','md.hydrology.datebreaks','format','DoubleMat','scale',yts); WriteData(fid,prefix,'data',tempmonthlyfactors,'name','md.hydrology.monthlyfactors','format','DoubleMat'); - WriteData(fid,prefix,'data',{'FrictionWaterPressure'},'name','md.hydrology.requested_outputs','format','StringArray'); + + %process requested outputs + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); end % }}} end Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 27530) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 27531) @@ -79,5 +79,5 @@ structstoch = self.structstochforcing() # Check if hydrologyarmapw is used - if type(md.hydrology).__name__ == 'hydrologyarmapw' and md.transient.ishydrology: + if((type(md.hydrology).__name__ == 'hydrologyarmapw') and md.transient.ishydrology==1): ispwHydroarma = 1 else: @@ -241,5 +241,5 @@ md.stochasticforcing.stochastictimestep = md.timestepping.time_step #by default: stochastictimestep set to ISSM time step # Check if hydroarmapw is used - if type(md.hydrology).__name__ == 'hydrologyarmapw' and md.transient.ishydrology: + if((type(md.hydrology).__name__ == 'hydrologyarmapw') and md.transient.ishydrology==1): ispwHydroarma = 1 else: Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m (revision 27530) +++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m (revision 27531) @@ -45,4 +45,5 @@ ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],... ['./' basename '-' ncdate '.nc'],... + '/media/vincent/TOSH4TB/GeorgiaTech/DataSearch/BedMachine/BedMachineGreenland-2021-04-20.nc',... }; @@ -89,7 +90,9 @@ if strcmp(string,'mask') | strcmp(string,'source'), %Need nearest neighbor to avoid interpolation between 0 and 2 - output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); + %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); + output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV else - output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); + %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); + output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV end Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpSeaRISE.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpSeaRISE.m (revision 27530) +++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpSeaRISE.m (revision 27531) @@ -58,4 +58,8 @@ searisenc='/totten_1/ModelData/SeaRISE/Antarctica_5km_dev1.0.nc'; end + case {'vincentxps159500'} + if hemisphere==1 + searisenc = '/media/vincent/TOSH4TB/GeorgiaTech/ISSM/fromSeroussi/GrIS04Feb2022/seariseData/Greenland_5km_dev1.2.nc'; + end otherwise error('hostname not supported yet'); @@ -80,7 +84,9 @@ if verbose, disp([' -- SeaRISE: interpolating ' string]); end if strcmpi(string,'LandMask'); - dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj,'nearest'); + %dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj,'nearest'); + dataout = InterpFromGridToMesh(xdata,ydata,data,xproj,yproj,NaN); %VV else - dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj); + %dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj); + dataout = InterpFromGridToMesh(xdata,ydata,data,xproj,yproj,NaN); %VV end dataout = reshape(dataout,size(X,1),size(X,2));