ice_dyn_vp: fgmres: exit early if right-hand-side vector is zero

If starting a run with with "ice_ic='none'" (no ice), the linearized
problem for the ice velocity A x = b will have b = 0, since all terms in
the right hand side vector will be zero:

- strint[xy] is zero because the velocity is zero
- tau[xy] is zero because the ocean velocity is also zero
- [uv]vel_init is zero
- strair[xy] is zero because the concentration is zero
- strtlt[xy] is zero because the ocean velocity is zero

We thus have a linear system A x = b with b=0, so we
must have x=0.

In the FGMRES linear solver, this special case is not taken into
account, and so we end up with an all-zero initial residual since
workspace_[xy] is also zero because of the all-zero initial guess
'sol[xy]', which corresponds to the initial ice velocity. This then
leads to a division by zero when normalizing the first Arnoldi vector.

Fix this special case by computing the norm of the right-hand-side
vector before starting the iterations, and exiting early if it is zero.
This is in line with the GMRES implementation in SciPy [1].

[1] https://github.com/scipy/scipy/blob/651a9b717deb68adde9416072c1e1d5aa14a58a1/scipy/sparse/linalg/_isolve/iterative.py#L620-L628

Close: #42

cicecore/cicedynB/dynamics/ice_dyn_vp.F90

@@ -2780,6 +2780,7 @@ subroutine fgmres (zetax2   , etax2   , &
       real (kind=dbl_kind) :: &
          norm_residual   , & ! current L^2 norm of residual vector
          inverse_norm    , & ! inverse of the norm of a vector
+         norm_rhs        , & ! L^2 norm of right-hand-side vector
          nu, t               ! local temporary values
 
       integer (kind=int_kind) :: &
@@ -2819,6 +2820,25 @@ subroutine fgmres (zetax2   , etax2   , &
       arnoldi_basis_x = c0
       arnoldi_basis_y = c0
 
+      ! solution is zero if RHS is zero
+      !$OMP PARALLEL DO PRIVATE(iblk)
+      do iblk = 1, nblocks
+         call calc_L2norm_squared(nx_block       , ny_block      , &
+                                  icellu   (iblk),                 &
+                                  indxui (:,iblk), indxuj(:,iblk), &
+                                  bx(:,:,iblk)   , &
+                                  by(:,:,iblk)   , &
+                                  norm_squared(iblk))
+
+      enddo
+      !$OMP END PARALLEL DO
+      norm_rhs = sqrt(global_sum(sum(norm_squared), distrb_info))
+      if (norm_rhs == c0) then
+         solx = bx
+         soly = by
+         return
+      endif
+
       ! Residual of the initial iterate
 
       !$OMP PARALLEL DO PRIVATE(iblk)