Hi all, I am playing around with the RT/BDM L2 error norm convergence for the Darcy equation: mu/K*v + grad(p) = rhob div(v) = 0 Where mu is viscosity, K is permeability tensor, and rhob is specific body force. The analytical solutions for a bi-unit square are: p(x,y) = sin(pi*x)*sin(pi*y) v_x(x,y) = sin(pi*x)*cos(pi*x) v_y(x,y) = -cos(pi*x)*sin(pi*y) And mu = K = 1. Theory states that RT0 has velocity/pressure slope of 1.0 and that BDM has velocity slope of 2.0. I get these or better. Now if I extend the Darcy formulation to include pressure-dependence for viscosity: mu(p)/K*v + grad(p) = rhob div(v) = f where mu(p) = mu*(1 + beta*p) So we now have a nonlinear formulation. By employing the same analytical solutions, the only thing that changes is rhob, but the velocity L2 errornorm slopes for both RT and BDM worsen. For example, BDM goes from 2.0 to 1.0. Pressure slopes remains unchanged. Why is this happening? Does this have something to do with the DG0 pressure functionspace that is now included in the mass matrix term? If I used a stabilized equal-order formulation for velocity/pressure, no such slope differences were observed. These trends were also observed with FEniCS. Attached is the code y'all can play around with: Run it as follows: python H_conv_barus.py Darcy RT 10 python H_conv_barus.py Darcy RT 20 python H_conv_barus.py Darcy RT 40 ... python H_conv_barus.py Darcy BDM 10 python H_conv_barus.py Darcy BDM 20 python H_conv_barus.py Darcy BDM 40 ... python H_conv_barus.py Barus RT 10 python H_conv_barus.py Barus RT 20 python H_conv_barus.py Barus RT 40 ... python H_conv_barus.py Barus BDM 10 python H_conv_barus.py Barus BDM 20 python H_conv_barus.py Barus BDM 40 ... I am trying to understand whether these H(div) elements are suitable for nonlinear mixed formulations or if there is a better way of implementing what I have. Thanks! Justin