Solving 3 time dependent PDE's simultaneously!!

Question

Hey, I have the following code for solving three different PDE's which are time dependent.
They are subjected to same mesh and all have same initial and boundary condition. During each iteration, I generate new vector element by applying conditions to solution of three PDE's and use this vector as variable values for t-1 seconds in all the PDE's.

My code is as follows:

from dolfin import*
import numpy
nx = ny = 50
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "CG", 1)

#plot(mesh)

Tc = Constant(310)
Tp = Constant(77.0)

def left(x, on_boundary):
    return on_boundary and x[0] <= 0.05 and x[1] <= 0.05

bcs = [DirichletBC(V, Tp, left)]

# Initial condition
Tinitial = 'x[0] <= 0.05 && x[1] <= 0.05 ? Tp : Tc'
Temp = Expression(Tinitial, Tp = Tp, Tc = Tc)
T_1 = interpolate(Temp, V)

#print T_1.compute_vertex_values()

dt = 1    # time step

# Parameters
omegabcb = 40 * 1E3
tb = 310
qmet = 33.8 * 1E3

k1 = Constant(0.5 * 1E4)
k2 = Constant(15.98 * 1E4)
k3 = 1005 * 1E4 / T_1**1.5

c1 = Constant(3.6 * 1E6)
c2 = Constant(880 * 1E6) 
c3 = Constant(0.00415 * 1E6)

# Define variational problem
T1 = TrialFunction(V)
v1 = TestFunction(V)
T2 = TrialFunction(V)
v2 = TestFunction(V)
T3 = TrialFunction(V)
v3 = TestFunction(V)

a1 = (c1 + omegabcb*dt)*T1*v1*dx + dt*k1*inner(nabla_grad(T1), nabla_grad(v1))*dx
L1 = (c1*T_1 + dt*omegabcb*tb + qmet*dt)*v1*dx
a2 = (c2)*T2*v2*dx + dt*k2*inner(nabla_grad(T2), nabla_grad(v2))*dx - 3.21*1E6*T2*T_1*v2*dx
L2 = (c2*T_1)*v2*dx - 3.21*1E6*T_1*T_1*v2*dx
a3 = dt*k3*inner(nabla_grad(T3), nabla_grad(v3))*dx + c3*T3*T_1*v3*dx
L3 = c3*T_1*T_1*v3*dx

A1 = None          
A2 = None          
A3 = None          
b1 = None         
b2 = None          
b3 = None          

# Compute solution
T1 = Function(V)  
T2 = Function(V)   
T3 = Function(V)   
duration = 70  # total simulation time
t = dt

while t <= duration:
    print 'time =', t
    A1 = assemble(a1, tensor = A1)
    A2 = assemble(a2, tensor = A2)
    A3 = assemble(a3, tensor = A3)
    b1 = assemble(L1, tensor=b1)
    b2 = assemble(L2, tensor=b2)
    b3 = assemble(L3, tensor=b3)
    for bc in bcs : bc.apply(A1, b1)
    for bc in bcs : bc.apply(A2, b2)
    for bc in bcs : bc.apply(A3, b3)
    solve(A1, T1.vector(), b1)
    solve(A2, T2.vector(), b2)
    solve(A3, T3.vector(), b3)
    T1mat = T1.vector().array()
    T2mat = T2.vector().array()
    T3mat = T3.vector().array()
    T_1mat = T_1.vector().array()
    T = []
    for i in range(0, len(T1mat)):
        if(T_1mat[i] >= 273):
            T.append(T1mat[i])
        if(T_1mat[i] >= 251 and T_1mat[i] < 273):
            T.append(T2mat[i])
        if(T_1mat[i] < 251):
            T.append(T3mat[i])
    T_1.vector()[:] = numpy.array(T)
    t += dt   

plot (T_1, title = 'T plot')
interactive()

When I run this code for duration greater than or equal to 60 secs. I get an error saying :

File "TestCode14.py", line 94, in
T_1.vector()[:] = numpy.array(T) File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1183, in
setslice
self.setitem(slice(i, j, 1), values) File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1214, in
setitem
_set_vector_items_array_of_float(self, indices, values)
RuntimeError: non matching dimensions on input

Please help in solving this problem.
Thanks

Answer 1

Hi, putting

T = []
#  loop ...
print len(T)

you'll see that before the code crashes len(T) is not equal V.dim() as it should. Inspecting
T1_mat reveals that it contains quite a few NaNs. You should check stability of your numerical scheme.

Comments

Hey,
I think this is occurring because I am generating T_1 vector (function at time t-1) using solution of three different PDE's each time which gives discontinuity and hence returns nan value after some iterations depending on the number of elements.
I am not sure if this is the correct way to address such problems. Can you suggest any other method to solve such problems where each time function at time t-1 is made using solution of three different PDE's over the domain.

Thanks