norm_residual is the norm of arnoldi_basis_[xy], which is b[xy] - workspace_[xy].
- b[xy] is all zero
- workspace_[xy] is matvec(sol[xy]) and sol[xy] is the initial iterate which is zero because of icnone (no ice so velocity is zero everywhere...).
bx (and similarly for y) is bx(i,j) = bxfix(i,j) + taux + strintx (calc_bvec)
- strintx is dPr/dx which is zero if velocity is zero
- taux = vrel(i,j)*waterxU(i,j) is zero because vrel is zero
  - vrel is zero because ocean velocity is zero (calc_vrel_Cb)
- bxfix is umassdti(i,j)*uvel_init(i,j) + forcexU(i,j)
  - uvel_init is zero
  - forcexU is forcex(i,j) = strairx(i,j) + strtltx(i,j) (dyn_prep2)
    - strtltx is zero if uocn is zero (dyn_prep2)
    - strairx is also all zero.
strairx from dyn_prep2 is strairxU in implicit_solver
strairxU is call grid_average_X2Y('F',strairxT,'T',strairxU,'U') if calc_strair = .true. (default)
- strairxT is initialized to 0 in init_coupler_flux and init_flux_atm.
- strairxT is "merged" across categories from strairxn (merge_fluxes) IF aicen_init(n) > puny (icepack_step_therm1) (which is false since aice_init = 0
- strairxn is computed in icepack_atm_boundary if calc_strair AND aicen_init(n) > puny, which is false if starting with no ice

So all of these fields being zero seems correct and makes sense.

We have A x = b with b = 0 so we must have x = 0. It is thus the FMGRES algorithm which has to be tweaked to account for b = 0. This is e.g. what Scipy does:
https://github.com/scipy/scipy/blob/651a9b717deb68adde9416072c1e1d5aa14a58a1/scipy/sparse/linalg/_isolve/iterative.py#L620-L628

VP solver: division by zero in FGMRES if starting with no ice #42

Description