Neumann BCs for a system of PDEs

Question

So I'm currently solving a system of coupled PDEs using fenics. I'm a little confused as to what I'm supposed to do with the Neumann BCs for the code though. I'm trying to solve

$$\frac{\partial P}{\partial t} - D_{p} \Delta P = - \lambda M $$
$$\frac{\partial M}{\partial t} - D_{M} \Delta M = kM $$
$$\frac{\partial P}{\partial n} = 0 $$
$$\frac{\partial M}{\partial n} = 0 $$

I have the weak formulation of this as follows, after using a Crank-Nicholson like scheme to discretize in time:

$$\int_\Omega \nu_P \frac{P^{(n+1)} - P^{(n)}}{\Delta t} - D_P \nabla P^{(n)} \cdot \nabla \nu_P + \lambda \mathbb{M} \nu_P d\Omega = 0$$
$$\int_\Omega \nu_M \frac{M^{n+1} - M^{n}}{\Delta t} - D_M \nabla M^{(n)} \cdot \nabla \nu_M - k \mathbb{M} \nu_M d\Omega = 0$$
$$\mathbb{M} = \frac{1}{2} (M^{(n+1)} + M^{(n)})$$

The tutorial for systems of PDEs doesn't really seem to discuss what I need to do for this. I'm solving this on the unit disk. I know that I next need to create the class for the system. After that I'm lost. Below is the code that I have so far; any help/advice that you can give would be greatly appreciated.

from cbc.pdesys import *
import numpy

# Set up the mesh to solve on
mesh = UnitCircle(20)
Q = FunctionSpace(mesh,'Q',1)

#establish time parameters
problem_parameters['time_integration'] = "Transient"
problem = Problem(mesh, problem_parameters)

# Set up first PDESystem
solver_parameters['space']['u'] = VectorFunctionSpace
solver_parameters['degree']['u'] = 2 
Plaque = PDESystem([['u', 'p']], problem, solver_parameters)

boundary = AutoSubDomain[lambda x, on_boundary: on_boundary]

Answer 1

Hi.

Sorry I missed that you was trying to use Crank-Nicolson method.
But here is a solution where I use Continuous Galerkin, hope you can use it anyway.

from dolfin import *
import numpy
import time

# set parameters
t = 0.0;            #Start time
dt = 0.1;          #Step size
T = 5             #Total simulation time

# define piecewise const diffusion coefficients
Dp = 1
Dm = 1

# Define reaction constants
Lambda = 1
k = 1


# Set up the mesh to solve on
circ = Circle(0.0,0.0,1)
mesh = Mesh(circ,80)

# Create function spaces
W1 = FunctionSpace(mesh, "CG", 1)
V = MixedFunctionSpace([W1,W1])

# Create intial conditions and functions
(P, M) = TrialFunctions(V)
(vp, vm) = TestFunctions(V)


## Class representing the intial conditions
class InitialConditions(Expression):
    def eval_cell(self, value, x, ufc_cell):
    # P0
        value[0] = 1.0
    # M0
        value[1] = 1.0
    def value_shape(self):
        return (2,)

c0 = Function(V)
c_init = InitialConditions()
c0.interpolate(c_init)
(P0, M0) = as_vector((c0[0], c0[1]))

c1 = Function(V)
(P1, M1) = as_vector((c1[0], c1[1]))


FP = (1/dt)*(P - P0)*vp *dx \
    + Dp * inner(grad(P), grad(vp)) * dx \
    + Lambda*M*vp* dx

FM = (1/dt)*(P - P0)*vm *dx \
    + Dm * inner(grad(M), grad(vm)) * dx \
    + k*M*vm* dx


F = FP + FM

a = lhs(F); L = rhs(F)

# preassembly
A = assemble(a)
b = None

t += dt
while t < T:
    b = assemble(L, tensor=b)
    solve(A, c1.vector(), b)

    c0.assign(c1)
    (P1, M1) = c1.split()

    print " time = %g, |c_trial| = %g" %(t, c1.vector().norm("l2"))

    plot(P1)
    time.sleep(5)
    plot(M1)
    time.sleep(5)

    t += dt

Comments

Thanks for this! I think I finally understand how to do this. The one problem I have when I run this code is that the plots are blank. Do you happen to know why that is? I get the scale but everything else is just blank.

---

Nope. I can see all the plots.
But perhaps you could try to only show one of the plots, and set interactive to true.
I hope it helps :)

---

The only thing I changed from your code to mine was using

mesh = UnitCircle(40)

Could this be causing the problem? When I ran your code, I got the error

Traceback (most recent call last):
  File "UnitCircle.py", line 20, in <module>
    circ = Circle(0.0,0.0,1)
NameError: name 'Circle' is not defined.