compute the residual and L2 norm
compute_residual(nxf, nyf, dx, dy, f, un, r)
rms = compute_l2norm(nxf, nyf, r)
if (rms/init_rms) <= 1e-12
break
end