Hi, I noticed something strange when I ran the enthalpy model. When I request Vx from the thermal solution via md.thermal.requested_outputs, the Vx does not match the Vx field when solving the Stressbalance. Everything else is kept the same. When solving Thermal, I am curious how exactly it uses the stress balance information. Given the discrepancy, it doesn't look like it solves stress balance.
The code is here:
% where ice velocity < 10m/a: make it no-slip at the base
epsilon = 1e-9;
v_idx_base = md2.initialization.vel < 10 & md2.mesh.vertexonbase;
md2.stressbalance.spcvx(v_idx_base & ~logical(md2.mesh.vertexonboundary)) = epsilon;
md2.stressbalance.spcvy(v_idx_base & ~logical(md2.mesh.vertexonboundary)) = epsilon;
% solve stress balance for comparison
md2 = solve(md2,'sb');
% solve thermal with the velocity field in output
md2.thermal.requested_outputs = {'default','Vz','Vx'};
md2 = solve(md2,'thermal');
% make a 2D plot looking at Vx at the base
plotmodel(md2d, 'data', project2d(md2, md2.results.StressbalanceSolution.Vx,1) -project2d(md2, md2.results.ThermalSolution.Vx, 1),...
'caxis', [-10,10]);
The figure for the difference in Vx at the ice base is attached
