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