3D Diffusion with Zero Neumann Boundary Conditions

Asked by Michael Trogdon on 2013-02-01

I am trying to run a time dependent diffusion model on a 3D mesh (and E. Coli bacteria) with Neumann boundary conditions equal to zero and no forcing term, f=0. I have drawn from the d2_d2D.py example problem and tried to adapt it to my problem. I have the mesh imported just fine and I think I have made the necessary adjustments for Neumann boundary conditions but when I run the model I am getting no evolution in time (i.e. my solution at t=T is the same as t=0). So I think one problem is that I can not come up with a simple example problem (that will satisfy the boundary conditions) to check if my code has an error or my problem set up. Any advice on how to go about this would be greatly appreciated. I have also put my code below if there are any glaring errors I would appreciate the feedback on that too. I am new to this forum so I am not sure how it works to share files but I can provide the .xml file with the mesh to anyone who wants to actually run the model. Thanks so much in advance for the help, I really appreciate it!!

"""This demo program solves the diffusion equation

     div grad u1 = du1/dt

with only BULK diffusion of one species on the ecoli mesh with homogeneous Neumann boundary conditions and no forcing term (i.e. no reactions).

from dolfin import *
import scipy.io
import numpy
#parameters["linear_algebra_backend"] = "uBLAS"

#Create mesh from Coli problem
mesh = Mesh('coli.xml')
V = FunctionSpace(mesh, "Lagrange", 1)
print mesh

#Diffusion constants for membrane and cytoplasm
Dmem = 1E-14
Dcyt = 2.5E-12

# Define initial conditions
alpha = 2; beta = 1; gamma=8;
u0 = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*x[2]*x[2] + gamma*t',
                alpha=alpha, beta=beta, gamma=gamma, t=0)
u0i = interpolate(u0, V)

T = 20 # total simulation time
dt = 0.3 # time step

# Define variational problem

# Laplace term
u = TrialFunction(V)
v = TestFunction(V)
a_K = Dcyt*inner(nabla_grad(u), nabla_grad(v))*dx

# "Mass matrix" term
a_M = u*v*dx

M = assemble(a_M)
K = assemble(a_K)
A = M + dt*K

# f term
#f = Expression('0')

# Compute solution
u = Function(V)
t = dt
while t <= T:
    #f_k = interpolate(f, V)
    #F_k = f_k.vector()
    b = M*u0i.vector() #+ dt*M*F_k
    u0.t = t
    solve(A, u.vector(), b)

    t += dt

#Try to get solution evaluated at vertex points instead of cell centers
VertexSolution = numpy.zeros(mesh.num_vertices(), dtype='d')
u.compute_vertex_values(VertexSolution, mesh)

Mbulk = M.array()
Kbulk = K.array()

# export to MATLAB as a dense matrix
scipy.io.savemat("rdBULKdense.mat", {"Mbulkdense": Mbulk, "Kbulkdense": Kbulk, "VertexSolution" : VertexSolution})

file = File("BulkDiffusion.pvd")
file << u


Question information

English Edit question
DOLFIN Edit question
No assignee Edit question
Solved by:
Michael Trogdon
Last query:
Last reply:
M Waste (michael-waste) said : #1


I'm not sure, but maybe change u0i.assign(u) to u0i.assign(u0).

Then I get a different result.



Michael Trogdon (mindgame252) said : #2

Thanks for the reply but I am not sure this solves the problem directly. While it does create in evolution in time it is not iterating to solve the problem it is simply updating an expression (the initial conditions) and plotting those. I think we need to update u0 to u, while in the loop, which is the solution from the first solve in the iteration. I think this is correct but let me know if I am missing something.

Thanks again for the reply.

Michael Trogdon (mindgame252) said : #3

Well I figured it out. The diffusion coefficients were too small because they came from a model that was in actual units while the gmsh file is using relative numbers. So change the Dcyt from 2.5E-12 to 2.5 and it diffuses just fine. So I think the script is working just fine.

Thanks for the help!