Index: /issm/trunk-jpl/examples/Greenland/CheatSheet.m =================================================================== --- /issm/trunk-jpl/examples/Greenland/CheatSheet.m (revision 20740) +++ /issm/trunk-jpl/examples/Greenland/CheatSheet.m (revision 20741) @@ -1,35 +1,70 @@ +if any(steps==7) + disp(' Step 7: Historical Relaxation run'); + md = loadmodel('./Models/Greenland.Control_drag'); + + load smbbox + + %convert mesh x,y into the Box projection + [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,39,71); + [xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70); + + %Interpolate and set surface mass balance + index = BamgTriangulate(x1(:),y1(:)); + smb_mo = InterpFromMeshToMesh2d(index,x1(:),y1(:),smbmean(:),xi,yi); + smb = smb_mo*12/1000*md.materials.rho_freshwater/md.materials.rho_ice; + md.smb.mass_balance = [smb;1 ]; + + %Set transient options, run for 20 years, saving every 5 timesteps + md.timestepping.time_step=0.2; + md.timestepping.final_time=200; + md.settings.output_frequency=5; + + %Additional options + md.inversion.iscontrol=0; + md.transient.requested_outputs={'IceVolume','TotalSmb', ... + 'SmbMassBalance'}; + md.verbose=verbose('solution',true,'module',true); + + %Go solve + md.cluster=generic('name',oshostname,'np',2); + md=solve(md,TransientSolutionEnum); + + save ./Models/Greenland.HistoricTransient_200yr md; +end + if any(steps==8) disp(' Step 8: Plotting exercise'); %Load historic transient model - md = loadmodel('./Models/Greenland.HistoricTransient'); + md = loadmodel('./Models/Greenland.HistoricTransient_200yr'); %Create Line Plots of relaxation run. Create a figure. figure - %Save surface mass balance, by looping through 200 years, or 1000 steps - surfmb=[]; for i=1:1000; surfmb=[surfmb ... + %Save surface mass balance, by looping through 200 years (1000 steps) + % Note, the first output will always contain output from time step 1 + surfmb=[]; for i=2:201; surfmb=[surfmb ... md.results.TransientSolution(i).SmbMassBalance]; end %Plot surface mass balance time series in first subplot - subplot(3,1,1); plot([0.2:0.2:200],mean(surfmb)); + subplot(3,1,1); plot([1:200],mean(surfmb)); %Title this plot Mean surface mass balance title('Mean Surface mass balance'); - %Save velocity by looping through 1000 steps - vel=[]; for i=1:1000; vel=[vel md.results.TransientSolution(i).Vel]; end + %Save velocity by looping through 200 years + vel=[]; for i=2:201; vel=[vel md.results.TransientSolution(i).Vel]; end %Plot velocity time series in second subplot - subplot(3,1,2); plot([0.2:0.2:200],mean(vel)); + subplot(3,1,2); plot([1:200],mean(vel)); %Title this plot Mean Velocity title('Mean Velocity'); - %Save Ice Volume by looping through 1000 steps - volume=[]; for i=1:1000; volume=[volume md.results.TransientSolution(i).IceVolume]; end + %Save Ice Volume by looping through 200 years + volume=[]; for i=2:201; volume=[volume md.results.TransientSolution(i).IceVolume]; end %Plot volume time series in third subplot - subplot(3,1,3); plot([0.2:0.2:200],volume); + subplot(3,1,3); plot([1:200],volume); %Title this plot Mean Velocity and add an x label of years Index: /issm/trunk-jpl/examples/Greenland/Greenland.par =================================================================== --- /issm/trunk-jpl/examples/Greenland/Greenland.par (revision 20740) +++ /issm/trunk-jpl/examples/Greenland/Greenland.par (revision 20741) @@ -69,14 +69,9 @@ md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness; -disp(' Single point constraint for continental model'); -%Set at least one vertex to velocity 0 so as to not get a singular problem (point on the wet peninsula) +disp(' Single point constraints'); +%Initialize single point constraint arrarys md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6); 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); -location = 1.0e+06 *[.32011 -2.2039]; -[dist pos]=min(sqrt((md.mesh.x - location(1)).^2 + (md.mesh.y - location(2)).^2)); -md.stressbalance.spcvx(pos) = 0; -md.stressbalance.spcvy(pos) = 0; -md.stressbalance.spcvz(pos) = 0; Index: /issm/trunk-jpl/examples/Greenland/runme.m =================================================================== --- /issm/trunk-jpl/examples/Greenland/runme.m (revision 20740) +++ /issm/trunk-jpl/examples/Greenland/runme.m (revision 20741) @@ -117,5 +117,13 @@ %Planview plots - plotmodel(md,'data',md.results.TransientSolution(end).Vel,'caxis',[1e-1 6000],'log', 10, 'title', 'Velocity (m/y)','gridded',11); + plotmodel(md,'data',md.results.TransientSolution(end).Vel,'log#1', 10, ... + 'caxis#1', [1e-1 1e4], ... + 'title#1', 'Velocity (m/y)',... + 'data',md.results.TransientSolution(1).SmbMassBalance,... + 'title#2', 'Surface Mass Balance (m/y)',... + 'data',md.results.TransientSolution(end).Thickness,... + 'title', 'Thickness (m)',... + 'data',md.results.TransientSolution(end).Surface,... + 'title', 'Surface (m)'); %Line Plots @@ -161,5 +169,5 @@ %plot a time series of total SMB - plot(t,totalsmb/1e9); title('Total Surface mass balance, Gt'); xlabel('year'); ylabel('Gt/yr'); + figure; plot(t,totalsmb/1e9); title('Total Surface mass balance, Gt'); xlabel('year'); ylabel('Gt/yr'); clear smbbox @@ -210,5 +218,6 @@ - %Save surface mass balance, by looping through 200 years, or 1000 steps + %Save surface mass balance, by looping through 200 years (1000 steps) + % Note, the first output will always contain output from time step 1 @@ -219,5 +228,5 @@ - %Save velocity by looping through 1000 steps + %Save velocity by looping through 200 years @@ -228,5 +237,5 @@ - %Save Ice Volume by looping through 1000 steps + %Save Ice Volume by looping through 200 years @@ -240,5 +249,5 @@ if any(steps==9) disp(' Step 9: Box Transient run'); - md = loadmodel('./Models/Greenland.HistoricTransient'); + md = loadmodel('./Models/Greenland.HistoricTransient_200yr'); %load past transient results