Changeset 20737


Ignore:
Timestamp:
06/14/16 14:24:14 (9 years ago)
Author:
seroussi
Message:

NEW: starting to design tutorial 07 on melting sensitivities

Location:
issm/trunk-jpl/examples/MeltingSensitivity
Files:
2 deleted
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/examples/MeltingSensitivity/runme.m

    r20725 r20737  
    1 %Which steps to be performed
    2 steps=[1] ;
     1step=[1];
    32
    4 %Run Steps
     3if step==1 %Transient Run #1
    54
    6 %Mesh Generation #1
    7 if any(steps==1)
     5        md = loadmodel('./Models/PIG.Control_drag');   
    86
    9         md.miscellaneous.name='PIG.Mesh_generation';
     7        md.inversion.iscontrol=0;
     8        md.transient.ismasstransport=1;
     9        md.transient.isstressbalance=1;
     10        md.transient.isgroundingline=1;
     11        md.transient.ismovingfront=0;
     12        md.transient.isthermal=0;
     13       
     14        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     15        md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1);
     16       
     17        md.timestepping.time_step=0.1;
     18        md.timestepping.final_time=10;
     19        md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
    1020
    11         %Mesh parameters
    12         domain =['./DomainOutline.exp'];
    13         hmax=40000;    % maximum element size of the final mesh
    14         hmin=5000;     % minimum element size of the final mesh
    15         hinit=10000;   % element size for the initial mesh
    16         gradation=1.7; % maximum size ratio between two neighboring elements
    17         err=8;         % maximum error between interpolated and control field
     21        md=solve(md,TransientSolutionEnum);
    1822
    19         % Generate an initial uniform mesh (resolution = hinit m)
    20         md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1);
    21 
    22         %ploting
    23         plotmodel(md,'data','mesh')
    24 
    25         % Load Velocities
    26         % http://nsidc.org/data/nsidc-0484.html
    27         nsidc_vel='../Data/Antarctica_ice_velocity.nc';         
    28 
    29         % Get necessary data to build up the velocity grid
    30         xmin = ncreadatt(nsidc_vel,'/','xmin');
    31         xmin = strtrim(xmin);  % this is a string, and we need to recover the double value
    32         xmin = xmin(1:end-2);  % get rid of the unit
    33         xmin = str2num(xmin);  % convert to double
    34        
    35         ymax = ncreadatt(nsidc_vel,'/','ymax');
    36         ymax = strtrim(ymax); 
    37         ymax = ymax(1:end-2); 
    38         ymax = str2num(ymax);
    39        
    40         nx = ncreadatt(nsidc_vel,'/','nx');
    41         ny = ncreadatt(nsidc_vel,'/','ny');
    42        
    43         spacing = ncreadatt(nsidc_vel,'/','spacing');
    44         spacing = strtrim(spacing);
    45         spacing = spacing(1:end-2); 
    46         spacing = str2num(spacing);
    47        
    48         % Get velocities (Note: You can use ncdisp('file') to see an ncdump)
    49         vx = double(ncread(nsidc_vel,'vx'));
    50         vy = double(ncread(nsidc_vel,'vy'));
    51        
    52         x=xmin+(0:1:nx)'*spacing;
    53         x=double(x);
    54         y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
    55         y=double(y);
    56 
    57         % Interpolate velocities onto coarse mesh
    58         vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
    59         vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
    60         vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
    61         clear vx vy x y;
    62 
    63         % Adapt the mesh to minimize error in velocity interpolation
    64         md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err);
    65        
    66         %ploting
    67         plotmodel(md,'data','mesh')
    68 
    69         % Convert x,y coordinates (Polar stereo) to lat/lon
    70         [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1);
    71        
    7223        % Save model
    73         save ./Models/PIG.Mesh_generation md;
     24        save ./Models/PIG.Transient md;
    7425end
    7526
    76 %Masks #2
    77 if any(steps==2)
     27if step==2 %High Melt #2
     28        md = loadmodel('./Models/PIG.Transient');       
    7829
    79         md = loadmodel('./Models/PIG.Mesh_generation');
     30        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     31        md.basalforcings.floatingice_melting_rate=60*ones(md.mesh.numberofvertices,1);
     32       
     33        md.timestepping.time_step=0.1;
     34        md.timestepping.final_time=10;
     35        md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
    8036
    81         % Load SeaRISe dataset for Antarctica 
    82         % http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
    83         searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
    84        
    85         %read thickness mask from SeaRISE
    86         x1=double(ncread(searise,'x1'));
    87         y1=double(ncread(searise,'y1'));
    88         thkmask=double(ncread(searise,'thkmask'));
    89        
    90         %interpolate onto our mesh vertices
    91         groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
    92         groundedice(groundedice<=0)=-1;
    93         clear thkmask;
     37        md=solve(md,TransientSolutionEnum);
    9438
    95         %fill in the md.mask structure
    96         md.mask.groundedice_levelset=groundedice; %ice is grounded for mask equal one
    97         md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie
    98 
    99         %ploting
    100         plotmodel(md,'data',md.mask.groundedice_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
    101        
    102         % Save model
    103         save ./Models/PIG.SetMask md;
     39        save ./Models/PIG.HighMelt md;
    10440end
    10541
    106 %Parameterization #3
    107 if any(steps==3)
     42if step==3 %High surface mass balance #3
     43        md = loadmodel('./Models/PIG.Transient');       
    10844
    109         md = loadmodel('./Models/PIG.SetMask');
    110         md = parameterize(md,'./Pig.par');
     45        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     46        md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1);
    11147
    112         % Use a MacAyeal flow model
    113         md = setflowequation(md,'SSA','all');
     48        md.smb.mass_balance=2*md.smb.mass_balance;
    11449       
    115         % Save model
    116         save ./Models/PIG.Parameterization md;
     50        md.timestepping.time_step=0.1;
     51        md.timestepping.final_time=10;
     52
     53        md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
     54
     55        md=solve(md,TransientSolutionEnum);
     56
     57        save ./Models/PIG.HighMelt md;
    11758end
    11859
    119 %Control Method #4
    120 if any(steps==4)
     60if step==4 %Ice Front retreat
     61        md = loadmodel('./Models/PIG.Transient');       
    12162
    122         md = loadmodel('./Models/PIG.Parameterization');
     63        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     64        md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1);
    12365
    124         % Control general
    125         md.inversion.iscontrol=1;
    126         md.inversion.maxsteps=20;
    127         md.inversion.maxiter=40;
    128         md.inversion.dxmin=0.1;
    129         md.inversion.gttol=1.0e-4;
    130         md.verbose=verbose('solution',true,'control',true);
     66        md.timestepping.time_step=0.1;
     67        md.timestepping.final_time=10;
     68        md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
    13169
    132         % Cost functions
    133         md.inversion.cost_functions=[101 103 501];
    134         md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
    135         md.inversion.cost_functions_coefficients(:,1)=1;
    136         md.inversion.cost_functions_coefficients(:,2)=1;
    137         md.inversion.cost_functions_coefficients(:,3)=8e-15;
     70        md=solve(md,TransientSolutionEnum);
    13871
    139         % Controls
    140         md.inversion.control_parameters={'FrictionCoefficient'};
    141         md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
    142         md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
    143 
    144         % Additional parameters
    145         md.stressbalance.restol=0.01;
    146         md.stressbalance.reltol=0.1;
    147         md.stressbalance.abstol=NaN;
    148 
    149         % Solve
    150         md.toolkits=toolkits;
    151         md.cluster=generic('name',oshostname,'np',2);
    152         md=solve(md,StressbalanceSolutionEnum);
    153 
    154         % Update model friction fields accordingly
    155         md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
    156 
    157         plotmodel(md,'data',md.friction.coefficient)
    158 
    159         % Save model
    160         save ./Models/PIG.Control_drag md;
     72        save ./Models/PIG.HighMelt md;
    16173end
    162 
    163 %Plot #5
    164 if any(steps==5)
    165 
    166         md = loadmodel('./Models/PIG.Control_drag');
    167 
    168         plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',...
    169                 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
    170                 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
    171                 'FontSize#all',12,...
    172                 'data',md.initialization.vel,'title','Observed velocity',...
    173                 'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',...
    174                 'data',md.geometry.base,'title','Bed elevation',...
    175                 'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',...
    176                 'colorbar#all','on','colorbartitle#1-2','[m/yr]',...
    177                 'caxis#1-2',([1.5,4000]),...
    178                 'colorbartitle#3','[m]', 'log#1-2',10);
    179 end
    180 
    181 %HO #6
    182 if any(steps==6)
    183 
    184         % Load Model
    185 
    186         % Disable inversion
    187 
    188         % Extrude Mesh
    189 
    190         % Set Flowequation
    191 
    192         % Solve
    193 
    194         % Save Model
    195 
    196 end
    197 
    198 %Plot #7
    199 if any(steps==7)
    200 
    201         mdHO = loadmodel('./Models/PIG.ModelHO');
    202         mdSSA = loadmodel('./Models/PIG.Control_drag');
    203 
    204         basal=find(mdHO.mesh.vertexonbase);
    205         surf=find(mdHO.mesh.vertexonsurface);
    206 
    207         plotmodel(mdHO,'nlines',3,'ncols',2,'unit#all','km','axis#all','equal',...
    208                                                 'xlim#all',[min(mdHO.mesh.x) max(mdHO.mesh.x)]/10^3,...
    209                                                 'ylim#all',[min(mdHO.mesh.y) max(mdHO.mesh.y)]/10^3,...
    210                                                 'FontSize#all',12,...
    211                                                 'data',mdHO.initialization.vel,'title','Observed velocity',...
    212                                                 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
    213                                                 'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
    214                                                 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
    215                                                 'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
    216                                                 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
    217                                                 'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
    218                                                 'colorbar#all','on','view#all',2,...
    219                                                 'colorbartitle#all','[m/yr]',...
    220                                                 'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
    221 end
Note: See TracChangeset for help on using the changeset viewer.