Hi there,
I'm not sure if I could post a really detailed question here, but I'm really stuck here now.
This time-dependent non-linear problem doesn't converge. I don't know what kind of techniques I can use to solve this problem. It will be greatly appreciated if you happen to have time to look at it, and give me some advice. Thanks.
The equations are like:
d(s*p,t) - div(p*s*grad(p)) = f1
d(s,t) - div(p*grad(s)) = f2
where s and p are unknowns, and essential boundary conditions are applied on the boundary.
I created a manufactured solution to this system and it is expressed as:
p_exact = 1 + x*x + y*y + t
s_exact = 1 + x*x + y*y + t
and the RHS source terms f1 and f2 are calculated according to this manufactured solution.
I tried the following code to calculate it, but it failed to converge at time around 254(s). Even before this break point, there is always a small difference between the analytical solution and calculation solution.
from dolfin import *
import numpy
mesh = UnitSquareMesh(12,12)
def boundaries(x, on_boundary):
return on_boundary
Q = FunctionSpace(mesh, "CG", 1)
S = FunctionSpace(mesh, "CG", 1)
W = Q*S
# exact solution
g = Expression(("1 + x[0]*x[0] + x[1]*x[1] + t","1 + x[0]*x[0] + x[1]*x[1] + t"), t = 0)
# previous and current solution
w0 = interpolate(g, W)
w = interpolate(g, W)
(p,s) = split(w)
(p0,s0)= split(w0)
# time variables
dt = Constant(0.1); t = float(dt); T = 2000.0
# test functions:
q, r = TestFunctions(W)
# right hand side source term
f1 = Expression("2*t - 8*x[0]*x[0]*(x[0]*x[0] + x[1]*x[1] + t +1)-8*x[1]*(x[0]*x[0] + x[1]*x[1] + t +1)+2*x[0]*x[0]+2*x[1]*x[1]-4*(x[0]*x[0] + x[1]*x[1] + t +1)*(x[0]*x[0] + x[1]*x[1] + t +1)+2", t = 0)
f2 = Expression("-8*x[0]*x[0] - 8*x[1]*x[1]-4*t-3", t = 0)
# essential boundary condition
bc = DirichletBC(W, g, boundaries)
# weak form
F1 = s*p*q*dx+ p*s*q*dx - s*p0*q*dx - p*s0*q*dx + dt*p*s*inner(nabla_grad(p), nabla_grad(q))*dx -dt*f1*q*dx
F2 = p*s*r*dx- p*s0*r*dx + dt*p*inner(nabla_grad(s), nabla_grad(r))*dx -dt*f2*r*dx
F = F1 + F2
while t <= T:
print t
g.t = t
f1.t = t
f2.t = t
solve(F == 0, w, bc, solver_parameters={"newton_solver":{"maximum_iterations": 20, "relative_tolerance": 1e-9,"absolute_tolerance":1e-9}})
# verify accuracy
w_e = interpolate(g, W)
maxdiff = numpy.abs(w_e.vector().array() - w.vector().array()).max()
print 'Max error, t=%.2f: %-10.3f' % (t, maxdiff)
# time increase
w0.vector()[:] = w.vector()
t += float(dt)
