check that stress computations follow Lemieux and Dupont 2020

[phil-blain]

Note:

stressp  = σ_1  = σ_11 + σ_22
stressm  = σ_2  = σ_11 - σ_22
stress12 = σ_12

[phil-blain]

• The strain rates are computed in stress_vp using uvel,vvel, which contains the converged velocities (u^k in Lemieux & Dupont 2020)
• zetaD is computed at the start of anderson_solver using {u,v}prev_k, i.e. u^{k-1} in the paper, and this value is passed in stress_vp to compute the stresses.
• the normalized stresses (σ^n_1 and σ^n_2 in paper, sig1, sig2 in the code) are computed in principal_stress using stressp,stressm,stress12
• the ice strength (strength in the code, P_p in the paper) is computed in icepack_ice_strength, which is called at the beginning of imp_solver, but does not depend on the velocity... maybe this is an error in the paper.

[phil-blain]

So the full formula is

σ_ij  = 2η(u^{k−1})\dot{ε}_ij(u^k) + 
        [ζ(u^{k−1}) − η(u^{k−1})]\dot{ε}_kk(u^k)δ_ij  
        − P(u^{k−1})δ_ij/2,

and here P is the replacement pressure, not the ice_strength (P_p).

The replacement pressure does depend on the velocity:

P = P_p Δ/Δ^*

[phil-blain]

The computations are already correct as P = P_p Δ/Δ^* = ζ * 2 Δ^* Δ / Δ^* = 2ζΔ