Closed phil-blain closed 1 year ago
norm_residual
is the norm of arnoldi_basis_[xy]
, which is b[xy] - workspace_[xy]
.
b[xy]
is all zeroworkspace_[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 zerotaux = vrel(i,j)*waterxU(i,j)
is zero because vrel
is zerovrel
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 zeroforcexU
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 iceSo 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
52fd683 is a missing piece for a complete fix.
Reproducer;
It is the first normalization by
norm_residual
that causes the abort.