Hi, everybody. I am a newbie in numerical simulation. I performed the simulation in Transient, but I found that the glacier front position in the first year in the result was different from the mask I set. What is the reason?
` md.inversion.iscontrol=0;
md.transient.ismasstransport=1;
md.transient.isstressbalance=1;
md.transient.isgroundingline=1;
md.transient.ismovingfront=1;
md.transient.isthermal=0;
md.transient.issmb = 1;
md.groundingline.migration='SubelementMigration';
md.groundingline.friction_interpolation='SubelementFriction1';
md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
%Dont use damage model
md.damage.D=zeros(md.mesh.numberofvertices,1);
md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
%Front and GL options
md.groundingline.migration='SubelementMigration';
md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
md.groundingline.friction_interpolation='SubelementFriction1';
md.transient.amr_frequency = timestep;
md.verbose.solution=1;
md.timestepping.time_step=timestep;
md.timestepping.start_time = 0;
md.timestepping.final_time= md.timestepping.start_time + Transient_year;
md.settings.output_frequency=1/md.timestepping.time_step;
md.cluster=generic('name',oshostname,'np',ncore);
md.toolkits=toolkits;
%% thicknessevolution(md), which will give you -div(Hv).
% linearbasalforcings
% floating ice shelf
geothermalflux=md.basalforcings.geothermalflux;
md.basalforcings=linearbasalforcings();
md.basalforcings.upperwater_melting_rate=0;
md.basalforcings.deepwater_melting_rate=30; % melt rate 30
md.basalforcings.deepwater_elevation = -450;
md.basalforcings.upperwater_elevation = -100;
% grounded ice
md.basalforcings.groundedice_melting_rate= 0 * ones(md.mesh.numberofvertices,1);
md.basalforcings.geothermalflux = geothermalflux; % nan * zeros(md.mesh.numberofvertices,1);
disp(' Set geothermal heat flux');
ncdata='/home/ISSM/trunk/examples/Data/Greenland_5km_dev1.2.nc';
x1 = ncread(ncdata,'x1');
y1 = ncread(ncdata,'y1');
gflux = ncread(ncdata,'bheatflx')';
md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux,md.mesh.x,md.mesh.y,0);
%% sin frontalforcing
fmelt_rate =2* ones(md.mesh.numberofvertices,1);
md.frontalforcings.meltingrate = fmelt_rate;
clear fmelt_rate melt_active
%%
md.calving = calvingvonmises();
md.calving.stress_threshold_floatingice = 1.5 * 1e5;
md.calving.stress_threshold_groundedice = 1 * 1e6;
md.calving.min_thickness=100; %m, default NaN
%% for permanent CF prediction
md.levelset.spclevelset = nan * ones(md.mesh.numberofvertices,1);
md.levelset.stabilization = 1;
md.levelset.kill_icebergs=1;
md.levelset.migration_max=10000; % -- maximum allowed migration rate (m/a)
%Constrain areas above sea level
pos = find(md.geometry.bed>0);
md.levelset.spclevelset(pos) = md.mask.ice_levelset(pos);
% Reset levelset mask on domain boundary
pos = find(md.mesh.vertexonboundary==1);
md.levelset.spclevelset(pos) = md.mask.ice_levelset(pos);
%Set BCs along domain boundaries
idxB=find(md.mesh.vertexonboundary==1); %vertices on boundary
md.masstransport.spcthickness=NaN(md.mesh.numberofvertices,1);
%%
% md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'};
md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation',...
'MaskOceanLevelset','MaskIceLevelset', 'IceMass','TotalSmb','SmbMassBalance',...
'BasalforcingsFloatingiceMeltingRate',...
'TotalCalvingFluxLevelset',... %Gt/r
'GroundinglineMassFlux',... %Gt/yr
'TotalCalvingMeltingFluxLevelset','IcefrontMassFluxLevelset',...
'TotalCalvingFluxLevelset','TotalGroundedBmb'};`