Issue with SetIceSheetBC

DominoAJones

Hello!

I've been using SetIceSheetBC for a while now and everything was fine when I checked it (which was regularly, to make sure instability wasn't related to issues with the boundary conditions). Until last week when I checked the spcvx and it looked like this (the same gaps show up when I plot spcvy):

I relaunched everything and ran through my setup again to make sure nothing had gone wrong. There's nothing fancy going on, I just set md=SetIceSheetBC(md);

Any reason this might be happening all of a sudden?! I thought my files might have gotten corrupted so re-installed ISSM just in case but I get the same result. The vx and vy components are correctly interpolated on the mesh under md.initialization and md.inversion. Should I just set the BC manually?

Thanks!

mathieumorlighem

Hello @DominoAJones something must be happening with how you set up the boundary condition because this looks like they are set "per element" and not "per vertex". Can you check the size of spcvx ? SetIceSheetBC definitely sets Dirichlet conditions per vertex so I think there must be something else in your code after you call SetIceSheetBC.
Cheers
Mathieu

DominoAJones

Hi @mathieumorlighem

I was seeing the same thing but even on a "per element" basis, it's not showing all the elements along the boundary. I really simplified my code to make sure there wasn't anything after I set the BC that could be causing this issue, let me know what you think:

md=model;
% Name and Coordinate system
md.miscellaneous.name='KNS';
md.miscellaneous.notes='Setting up initialization';
md.mesh.epsg=3413;

% Creating a mesh
md=triangle(model,'.\KNS-setup\domain\domain.exp',1000000);

disp(' Refining Mesh');
ncdatavx = '.\trunk\datasets\Velocity\Clipped\vx_vel-CL.nc';
ncdatavy = '.\trunk\datasets\Velocity\Clipped\vy_vel-CL.nc' ;

velx= ncread(ncdatavx,'Band1');
xx= ncread(ncdatavx,'x');
yx= ncread(ncdatavx,'y');

vely= ncread(ncdatavy,'Band1');
xy= ncread(ncdatavx,'x');
yy= ncread(ncdatavx,'y');

vx = InterpFromGridToMesh(xx,yx,velx',md.mesh.x,md.mesh.y,0);
vy = InterpFromGridToMesh(xy,yy,vely',md.mesh.x,md.mesh.y,0);
vel = sqrt(vx.²+vy.²);
in=ContourToNodes(md.mesh.x,md.mesh.y,'.\KNS-setup\domain\no-ice-mask.exp',1);
vx(find(in))=0;
vy(find(in))=0;
vel(find(in))=0;

% Refine mesh using surface velocities as metric
md=bamg(md,'hmin',500,'hmax',1000000,'field',vel,'err',5);
[md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,39,71);

disp('Refine mesh at TW margin');
md.private.bamg=struct();

h=NaN*ones(md.mesh.numberofvertices,1);
in=ContourToNodes(md.mesh.x,md.mesh.y,'.\KNS-setup\domain\refinement.exp',1);
h(find(in))=500; % may be too fine

md=bamg(md,'hmin',500,'hmax',50000,'field',vel,'err',5,'hVertices',h);
% % vel=shock(md.mesh.x,md.mesh.y);

disp('Interpolate geometry');
nctopg='.\trunk\datasets\Bedmachine\Clipped\bed-bedmachine-CL.nc';
x2= ncread(nctopg,'x');
y2= ncread(nctopg,'y');
topg = ncread(nctopg,'bed')';
disp(' Bedrock topography');
md.geometry.base = InterpFromGridToMesh(x2,y2,topg,md.mesh.x,md.mesh.y,0);
md.geometry.bed=md.geometry.base;

ncsurf='.\trunk\datasets\Bedmachine\Clipped\surface-bedmachine-CL.nc';
x3= ncread(ncsurf,'x');
y3= ncread(ncsurf,'y');
surf = ncread(ncsurf, 'surface')';
disp(' Surface elevation');
md.geometry.surface=InterpFromGridToMesh(x3,y3,surf,md.mesh.x,md.mesh.y,0);
in=ContourToNodes(md.mesh.x,md.mesh.y,'.\KNS-setup\domain\no-ice-mask.exp',1);
md.geometry.surface(find(in))=md.geometry.base(find(in));

md.geometry.thickness=md.geometry.surface-md.geometry.base;
md.inversion.thickness_obs=md.geometry.thickness;

md.inversion.surface_obs=md.geometry.surface;

disp('Interpolate other parameters');
disp(' Velocity');
ncdatavv = '.\trunk/datasets/Velocity/Clipped/cal_vel-CL.nc';
velv= ncread(ncdatavv,'Band1');
xv= ncread(ncdatavv,'x');
yv= ncread(ncdatavv,'y');
md.inversion.vel_obs = InterpFromGridToMesh(xv,yv,velv',md.mesh.x,md.mesh.y,0);
md.inversion.vel_obs(find(in))=0;
md.initialization.vel= md.inversion.vel_obs;

vx = InterpFromGridToMesh(xx,yx,velx',md.mesh.x,md.mesh.y,0);
vy = InterpFromGridToMesh(xy,yy,vely',md.mesh.x,md.mesh.y,0);
md.inversion.vx_obs = vx;
md.inversion.vx_obs(find(in))=0;
md.inversion.vy_obs = vy;
md.inversion.vy_obs(find(in))=0;
md.initialization.vx = md.inversion.vx_obs;
md.initialization.vy = md.inversion.vy_obs;
md.initialization.vz=zeros(md.mesh.numberofvertices,1);

disp(' Temperature');
ncdatatemp='C:/Users/domjon/Desktop/KNS-ISSM/trunk/datasets/RACMO23p2/Clipped/t2m_clipped.nc';
xtemp= ncread(ncdatatemp,'x');
ytemp= ncread(ncdatatemp,'y');
temp_month = [];
for i = 25:385 % 30 year average from 1960 to 1990, note data is given in bands each month since 1958-01-15
varName = sprintf('Band%d', i);
tempBand = ncread(ncdatatemp, varName);
md.initialization.temperature = InterpFromGridToMesh(xtemp, ytemp, tempBand', md.mesh.x, md.mesh.y, 0);
temp_month = cat(3, temp_month, md.initialization.temperature);
end
md.initialization.temperature = mean(temp_month, 3);

disp(' Construct ice rheological properties'); %#md.materials
md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
md.materials.rheology_B=cuffey(md.initialization.temperature);
md.damage.D=zeros(md.mesh.numberofvertices,1);
md.materials.rheology_law='None';
% law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeH2O', or 'NyeCO2'

md=setflowequation(md,'SSA','all');

disp('Set mask');
in=ContourToNodes(md.mesh.x,md.mesh.y,'.\KNS-setup\domain\no-ice-mask.exp',1);
md=setmask(md,'','');
md.mask.ice_levelset(find(in))=1;

disp('Set boundary conditions');
md=SetIceSheetBC(md);

mathieumorlighem

Your code looks good, can you please check the size of md.stressbalance.spcvx and also plot plot model(md,'data','BC') to check the boundary conditions? (maybe it is just a plotting artifact of MATLAB)
Mathieu

DominoAJones

Plotting the BCs, everything looks fine:

I'll try updating MATLAB. What an odd issue to have crop up out of nowhere, hopefully it's ok to proceed regardless. 🙂

mathieumorlighem

Out of curiosity, how did you generate the first plot where the triangles show up? I get something similar if I use the option 'shading','flat' but not if I use the default plotmodel(md,'data',md.stressbalance.spcvx). Anyway I think all is well with your model.

DominoAJones

Yeah, I did run it with option 'shading', 'flat'. When I just use the default, nothing gets plotted. Hopefully this isn't a problem that pops up with other plots.

I initially thought there was something wrong with the BC because there's instability and sometimes total collapse if I run transients for longer than 50 years. If I run something and do a transientrestart, I can buy myself some more time but sooner or later it crashes. Thought I'd mention it in case you have any advice or suggestions of parameters I should be paying attention to.

Cheers!

mathieumorlighem

oh ok so then all is good!! and nothing has changed, just the way you plot things. NaN will turn "white" so that's why MATLAB plots nothing if you plot md.stressbalance.spcvx. You can try the plot option 'nan',0, which will replace the NaNs by 0 when you plot.
If you have a crash, the first thing I would look into is time stepping. If you are using a fixed time step make sure it is small enough. I recommend using timesteppingadaptive so that it adjusts automatically.

DominoAJones

Right, but there aren't any NaN values at this stage of the boundary conditions. I add NaN values later on, beyond where the tidewater margin advances into the fjord (on the left side of the domain) but they shouldn't show up all over (e.g. the big gaps that show up upstream, on the right side). I did run the plot option 'nan', 0 and it shows NaN values where I need spc for the upstream ice.

Haha, I shouldn't be surprised! Time steps are the bane of my existence - I am using adaptive timestepping but the minimum never seems to be small enough.

mathieumorlighem

SetIceSheetBC is setting boundary conditions along all sides as 0 or vx_obs. Take a look at the function for yourself

%node on Dirichlet
pos=find(md.mesh.vertexonboundary);
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);
md.stressbalance.spcvx(pos)=0;
md.stressbalance.spcvy(pos)=0;
md.stressbalance.spcvz(pos)=0;

if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices)
   disp('      boundary conditions for stressbalance model: spc set as observed velocities');
   md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
   md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
else
   disp('      boundary conditions for stressbalance model: spc set as zero');
end

If your time steps become super small that's because the velocity is blowing up. You have to identify why it becomes so large (plot the velocity at different steps, look at the geometry, etc). If you can't figure it out send me your model and I will take a look when I get a chance.
Cheers
Mathieu

DominoAJones

Right, and that's exactly the problem. I know it is supposed to set the boundary conditions as either 0 or vx_obs/vy_obs on all sides but it doesn't. For some reason, it sets some BCs as NaN and I do not understand why or how to fix it.

I'll send you the model, I would really appreciate some help sorting out the instability and BC issue. Cheers!

mathieumorlighem

Hi @DominoAJones

Maybe I am missing something, but I think the model is perfectly ok and consistent with what you said. spcvx/vy is of size numberofvertices and we use NaN for any node that should have "no constraint". If you look at your plotmodel(md,'data','BC') you will see that all the nodes along the boundary have a Dirichlet condition and no constraint inside.
Mathieu

DominoAJones

Any chance you'd have time for a chat? I want to make sure we're on the same page. Potentially everything's fine but I have a few follow up questions.