function R = calc_velocity_residual( grid, ice)
% Calculate the L2-norm of the two consecutive velocity solutions

res1 = sum( (ice.u_bk - ice.u_prev_bk).^2 + (ice.v_bk - ice.v_prev_bk).^2);
res2 = sum( (ice.u_bk + ice.u_prev_bk).^2 + (ice.v_bk + ice.v_prev_bk).^2);

R = 2 * res1 / max( res2, 1E-8);

end