"""
Flow between rotating concentric cylinders
"""
from dolfin import *

r0 = 0.2
r1 = 1
mesh = Interval(100, r0, r1)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
x = Expression("x[0]")
w0 = Constant(1.)
w1 = Constant(0.1)
T0 = Constant(8)
T1 = Constant(10)
mu = Constant(0.25)
k  = Constant(0.001)
u_ = Function(V)
T_ = Function(V)
F  = inner(grad(v), grad(u))*x*dx + u*v/x*dx
FT = inner(grad(v), grad(u))*x*dx - mu/k*v*x*(grad(u_)-u_/x)**2*dx

def inner_bnd(x, on_bnd):
    return on_bnd and x[0] < r0+DOLFIN_EPS

def outer_bnd(x, on_bnd):
    return on_bnd and x[0] > 1-DOLFIN_EPS    

bci = DirichletBC(V, r0*w0, inner_bnd)
bco = DirichletBC(V, r1*w1, outer_bnd)
bcTi = DirichletBC(V, T0, inner_bnd)
bcTo = DirichletBC(V, T1, outer_bnd)

solve(lhs(F) == rhs(F), u_, bcs=[bci, bco])
solve(lhs(FT) == rhs(FT), T_, bcs=[bcTi, bcTo])

# Compute exact solution given in (3.22)
ue = project(r0*w0*(r1/x-x/r1)/(r1/r0-r0/r1) + 
             r1*w1*(x/r0-r0/x)/(r1/r0-r0/r1), V)

# Compute exact solution given in (3.23)
Tn = project((T_-T0)/(T1-T0), V)
PrEc = mu*r0**2*w0**2/k/(T1-T0)
Te = PrEc*r1**4/(r1**2-r0**2)**2*((1-r0**2/x**2) - (1-r0**2/r1**2)*ln(x/r0)/ln(r1/r0)) + ln(x/r0)/ln(r1/r0)
Te = project(Te, V)

print "\nError velocity    = ", errornorm(u_, ue, degree=0)
print "Error temperature = ", errornorm(Tn, Te, degree=0)
print "??? Should be close to zero.\n"

plot(u_, title="Computed axial velocity")
plot(ue, title="Exact axial velocity (3.22)")
plot(Tn, title="Computed temperature")
plot(Te, title="Exact temperature (3.23)")

# Since I cannot reproduce the temperature I have recomputed
# (3.23) analytically. I plot that result here
C1 = w0-r1**2*(w1-w0)/(r0**2-r1**2)
C2 = r1**2*(w1-w0)/(1-r1**2/r0**2)
C3 = (T1-T0+C2**2*mu/k*(1/r1**2-1/r0**2))/ln(r1/r0)
C4 = T0 + C2**2*mu/k/r0**2-C3*ln(r0)
uee = C1*x+C2/x
Tee = -C2**2*mu/k/x**2+C3*ln(x)+C4
Tee = project((Tee-T0)/(T1-T0), V)
plot(uee, mesh=mesh, title="Recomputed velocity")
plot(Tee, mesh=mesh, title="Recomputed temperature")

print "\nError recomputed temperature = ", errornorm(Tn, Tee, degree=0)
print "Which is a strong indication that (3.23) is wrong!!"