WIP: analytic jacobian

This resolves #15.

The numeric Jacobian implementation will be replaced by an explicit analytic implementation.

• Find expressions for the derivatives of the conductivity K(\Theta) and the saturation \Theta(h):

  \Theta(h) = (1 + ( \alpha h )^n)^{-m}

  \rightarrow \quad \frac{\partial}{\partial h} \Theta (h) = - m n \alpha^n h^{n-1} ( 1 + ( \alpha h )^n )^{-m-1}

  K(\Theta) = K_0 \Theta^\tau ( 1 - ( 1 - \Theta^{1/m} )^m )^2

  \rightarrow \quad \frac{\partial}{\partial \Theta} K (\Theta) = K_0 \left[ \tau \Theta^{\tau-1} ( 1 - ( 1 - \Theta^{1/m} )^m )^2 + 2 \Theta^\tau \Theta^{- 1 + 1/m} ( 1 - \Theta^{1/m} )^{m-1} ( 1 - ( 1 - \Theta^{1/m} )^m ) \right]

• add jacobian_ implementations

  • jacobian_volume:

    \alpha_i = K ( \Theta ( h ) ) \left[ (\nabla h - \mathbf{\hat{g}}) \nabla \varphi_i \right] \quad \rvert \, h = \sum_k u_k \varphi_k , \, \nabla h = \sum_k u_k \nabla \varphi_k

    \rightarrow \quad \frac{\partial}{\partial u_j} \alpha_i = \frac{\partial K}{\partial \Theta} \frac{\partial \Theta}{\partial h} \varphi_j \left[ (\nabla h - \mathbf{\hat{g}}) \nabla \varphi_i \right] + K ( \Theta ( h ) ) \left[ \nabla \varphi_j \nabla \varphi_i \right]

  • jacobian_skeleton:

    \alpha^s_i = K_1 ( \Theta ( h^u ) ) j_N \varphi^s_i + \theta \cdot \omega^s K_0^s K_1 ( \Theta ( h^u ) ) ( \nabla \varphi^s_i \cdot \mathbf{\hat{n}} ) [ h ]

    \alpha^n_i = - K_1 ( \Theta ( h^u ) ) j_N \varphi^n_i + \theta \cdot \omega^n K_0^n K_1 ( \Theta ( h^u ) ) ( \nabla \varphi^n_i \cdot \mathbf{\hat{n}} ) [ h ]

    [ h ] = h^s - h^n \quad \rightarrow \quad \frac{ \partial }{ \partial u^s_j } [ h ] = \varphi^s_j \, , \, \frac{ \partial }{ \partial u^n_j } [ h ] = - \varphi^n_j

    \text{if } j_N > 0 \, : \, h^u = h^s \quad \rightarrow \quad \frac{ \partial }{ \partial u^s_j } h^u = \varphi^s_j

    \text{else } \, : \, h^u = h^n \quad \rightarrow \quad \frac{ \partial }{ \partial u^n_j } h^u = \varphi^n_j

    K_1 ( \Theta ) = K ( \Theta ) / K_0

    \rightarrow \quad \frac{ \partial }{ \partial u^s_j } K_1 ( h^u ) = \frac{ 1 }{ K_0 } \frac{ \partial K }{ \partial \Theta^u } \frac{ \partial \Theta^u }{\partial h^u} \frac{ \partial h^u }{ \partial u^s_j }

    j_N = \omega^s j^s + \omega^n j^n \, , \quad j^s = K_0^s ( - ( ( \nabla h^s - \mathbf{\hat{g}} ) \cdot \mathbf{\hat{n}} ) + p [ h ] )

    \rightarrow \quad \frac{ \partial }{ \partial u^s_j } j_N = \omega^s K_0^s ( - ( \nabla \varphi^s_j \cdot \mathbf{\hat{n}} ) + p \varphi^s_j ) + \omega^n K_0^n p \varphi^s_j

    \rightarrow \quad \frac{ \partial }{ \partial u^n_j } j_N = \omega^s K_0^s ( - p \varphi^n_j ) + \omega^n K_0^n ( - ( \nabla \varphi^n_j \cdot \mathbf{\hat{n}} ) - p \varphi^n_j )

    \frac{\partial}{\partial u^s_j} \alpha^s_i = \frac{\partial}{\partial u^s_j} K_1 ( h^u ) \cdot j_N \varphi^s_i + K_1 ( h^u ) \cdot \frac{\partial}{\partial u^s_j} j_N \cdot \varphi^s_i \, + \, \theta \omega^s K_0^s ( \nabla \varphi^s_i \cdot \mathbf{\hat{n}} ) ( \frac{\partial}{\partial u^s_j} K_1 ( h^u ) + \frac{\partial}{\partial u^s_j} [ h ] )

  • jacobian_volume (temporal OP):

    \alpha_i = \theta_w (h) \varphi_i

    \rightarrow \quad \frac{\partial}{\partial u_j} \alpha_i = \frac{\partial \theta_w}{\partial h} \varphi_j \varphi_i = (\theta_s - \theta_r) \frac{\partial \Theta}{\partial h} \varphi_j \varphi_i

• check residuals and convergence rates