Scalar Mixed Function Space Problem

Question

Hello
I am trying to solve two coupled equations (Nernst-Planck and Poisson) for a simple domain. First of all I would like to explain my equations and boundary conditions. The equations are:

For the above equations c and phi are the unknowns corresponding to the concentration and electrical potential respectively. In addition, epsilon is constant = 0.0082
The domain is a square as shown below:

The Dirichlet boundary condition should be applied on the top and the bottom boundaries:

For both unknowns, zero Neumann boundary condition should be applied on all of the boundaries.
In order to derive the weak form of the Eq.1 we can multiply the Eq.1 by a test function (v_1) :

After integration by parts and applying zero Neumann boundary conditions:

After discretizing the first term in the Eq.4:

The c_0 in the above equation corresponds to the solution in the previous time step.Finally the weak form of the Eq.1 can be derived as:

We can derive the weak form of the Eq.2 easily by multiplying it by another test function (v_2):

After integration by parts (and considering zero Neumann B.C) the weak form of the Eq.2 can be derived as:

The weak form of the whole system is obtained by summing up the Eq.6 and Eq.8:

As far as I know, I should use a mixed function space for the two unknowns in this problem (c & phi). I have implemented it in a code but I am getting an error. This is my complete code:

from dolfin import *
import mshr

# Define boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Domain
domain = mshr.Rectangle(Point(0,0), Point(1,1))

# Mesh
mesh = mshr.generate_mesh(domain, 20)

# Marking boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)

#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, 'CG', 1)

#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, 'CG', 1)

#Defining the mixed function space
W = MixedFunctionSpace((V_c, V_phi))

#Defining the Trial  functions
(c, phi) = TrialFunctions(W)

#Defining the test  functions
(v_1, v_2) = TestFunctions(W)

# Time variables
dt = 0.05; t = 0; T = 1.1

# Previous solution
W_const = Constant((1200))
C_previous= interpolate(W_const, V_c)

# Define Dirichlet boundary conditions at top boundary (V=4.0 volt)
Voltage = Constant((4.0))
bc_top = DirichletBC(W.sub(1), Voltage, 1)

# Define Dirichlet boundary conditions at bottom boundary (V=0.0 volt)
bc_bottom = DirichletBC(W.sub(1), 0.0, 3)

#Collecting the boundary conditions
bcs = [bc_top, bc_bottom]

#Define eps in the weak form
eps=0.0082

# Define variational form
a = dot(grad(phi),grad(v_2))*dx - (1./(2*(pow(eps,2))))*c*v_2*dx\
    + (1./(2*(pow(eps,2))))*v_2*dx+c*v_1*dx\
     +dt *eps*dot(grad(c),grad(v_1))*dx\
     +dt*eps*c*dot(grad(phi),grad(v_1))*dx

L = C_previous*v_1*dx

u = Function(W)

while t <= T:
    A, b = assemble_system(a, L, bcs)

    solve(A, u.vector(), b)
    (c, phi) = u.split()

    C_previous.assign(c)

    t += dt
    plot(c, interactive=False)

This is the error I get:

#line 87, in <module>
    A, b = assemble_system(a, L, bcs)

I think there is something wrong in the defining the variational form that should be exactly the Eq.9. In addition to that I am not sure about this line:

    solve(A, u.vector(), b)

As you can see in the code, finally I want to split the mixed space and for example plot the concentration (c) over the time.
Could you please help get this code to work? Thanks in advance for your time.

Comments

Your system of equations is nonlinear if solve simultaneously. C.f. dot(c*grad(phi), grad(v1)).

---

I did what you said but still the error exists!

Answer 1

Hi,

you have several issues there. As nate pointed, the problem in nonlinear. Also correct BC application is needed.
Moreover, your splitting into bilinear and (uni)linear functional is wrong. If it was a linear problem then the term + (1./(2 * (pow(eps,2)))) * v_2 * dx has to be in functional L.

Anyway, I have edited your code for the nonlinear solver. The problem with assign is solved
setting deepcopy=True in u.split(True).

from dolfin import *
import mshr


# Define boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)


class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)


# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Domain
domain = mshr.Rectangle(Point(0,0), Point(1,1))

# Mesh
mesh = mshr.generate_mesh(domain, 20)

# Marking boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)

CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)

#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, CG1_elem)

#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, CG1_elem)

#Defining the mixed function space
W_elem = MixedElement([CG1_elem, CG1_elem])
W = FunctionSpace(mesh, W_elem)

#Defining the "Trial" functions
u = Function(W)
c, phi = split(u)

#Defining the test  functions
(v_1, v_2) = TestFunctions(W)

# Time variables
dt = 0.05; t = 0; T = 10

# Previous solution
W_const = Constant(1.0)
C_previous= interpolate(W_const, V_c)

# Define Dirichlet boundary conditions at top boundary (V=4.0 volt)
Voltage = Constant(4.0)
bc_top = DirichletBC(W.sub(1), Voltage, boundaries, 1)

# Define Dirichlet boundary conditions at bottom boundary (V=0.0 volt)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)

#Collecting the boundary conditions
bcs = [bc_top, bc_bottom]

#Define eps in the weak form
eps=0.0082

# Define variational form
a = dot(grad(phi), grad(v_2)) * dx \
    + c*v_1 * dx \
    - (1./(2 * (pow(eps,2)))) * c * v_2 * dx \
    + dt * eps * dot(grad(c),grad(v_1)) * dx \
    + dt * eps * c * dot(grad(phi),grad(v_1)) * dx

L = C_previous * v_1 * dx - (1./(2 * (pow(eps,2)))) * v_2 * dx

# For nonlinear solver
F = a - L

while t <= T:
    solve(F == 0, u, bcs)

    (c, phi) = u.split(True)
    C_previous.assign(c)

    t += dt
    plot(c, key="c", interactive=False)

Comments

Thank you for the time you spent on this code. When I run your code I get this error:

 V_c = FunctionSpace(mesh, CG1_elem)
TypeError: __init__() takes at least 4 arguments (3 given)

Anyway I changed the V_c and V_phi this way:

V_c = FunctionSpace(mesh, CG1_elem,1)
V_phi = FunctionSpace(mesh, CG1_elem,1)

I run the code again and Unfortunately this time I got another error:

*** Error:   Unable to create function space.
*** Reason:  Illegal argument for finite element family, not a    string: <CG1 on a <Domain built from <triangle cell in 2D> with label None>>.
*** Where:   This error was encountered inside functionspace.py.

Do you have any idea to fix these errors?
Thanks!

---

What version of FEniCS do you use? Just run command dolfin-version in terminal to see.

You should probably use the old way of spaces definition - as you had in the original code.

from dolfin import *
import mshr


# Define boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)


class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)


# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Domain
domain = mshr.Rectangle(Point(0,0), Point(1,1))

# Mesh
mesh = mshr.generate_mesh(domain, 20)

# Marking boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)

#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, 'CG', 1)

#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, 'CG', 1)

#Defining the mixed function space
W = MixedFunctionSpace((V_c, V_phi))

#Defining the "Trial" functions
u = Function(W)
c, phi = split(u)

#Defining the test  functions
(v_1, v_2) = TestFunctions(W)

# Time variables
dt = 0.05; t = 0; T = 10

# Previous solution
W_const = Constant(1.0)
C_previous= interpolate(W_const, V_c)

# Define Dirichlet boundary conditions at top boundary (V=4.0 volt)
Voltage = Constant(4.0)
bc_top = DirichletBC(W.sub(1), Voltage, boundaries, 1)

# Define Dirichlet boundary conditions at bottom boundary (V=0.0 volt)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)

#Collecting the boundary conditions
bcs = [bc_top, bc_bottom]

#Define eps in the weak form
eps=0.0082

# Define variational form
a = dot(grad(phi), grad(v_2)) * dx \
    + c*v_1 * dx \
    - (1./(2 * (pow(eps,2)))) * c * v_2 * dx \
    + dt * eps * dot(grad(c),grad(v_1)) * dx \
    + dt * eps * c * dot(grad(phi),grad(v_1)) * dx

L = C_previous * v_1 * dx - (1./(2 * (pow(eps,2)))) * v_2 * dx

# For nonlinear solver
F = a - L

while t <= T:
    solve(F == 0, u, bcs)

    (c, phi) = u.split(True)
    C_previous.assign(c)

    t += dt
    plot(c, key="c", interactive=False)

---

Thank you so much. It fixed the errors.
I am using the version 1.6.0. You are absolutely right. I had to bring the term (1./(2 * (pow(eps,2)))) * v_2 * dx in functional L (not a).
As the last question: How can I extract the results in an arbitrary point of the domain? Lets say I want to know the value of the "c" (or "phi") at the time = 5 in the center of the square domain [point:(0.5,0.5)]. How can I do that?

---

Please read the section 5.4.3 in https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf

In short, you can use

c(0.5, 0.5)

but this is very expensive and I would not use it for any other finite element space than CG1.

---

That was really helpful. Thanks again.

Answer 2

It is non-linear (as pointed out by nate), but you also need to apply the BC correctly:

 bc_top = DirichletBC(W.sub(1), Voltage, boundaries, 1)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)

Comments

I modified the boundary condition based on what you suggested. The problem here is related to the writing of the variational form. There are two problems here:

1- Variational form :
I know this problem is non-linear. I derived the weak form equations (Eq.6 and Eq.8) that correspond to the Eq.1 and Eq.2 respectively and added them to find the final weak form (Eq.9):

a = dot(grad(phi),grad(v_2))*dx \
    - (1./(2*(pow(eps,2))))*c*v_2*dx\
    + (1./(2*(pow(eps,2))))*v_2*dx\
    +c*v_1*dx\
     +dt *eps*dot(grad(c),grad(v_1))*dx\
     +dt*eps*c*dot(grad(phi),grad(v_1))*dx 

I removed the terms in the above form one by one and I realized if I only keep the first and the 4th term, I do not get that error again, like:

a = dot(grad(phi),grad(v_2))*dx \
    +c*v_1*dx

2- Defining the initial condition
I define the c_0 (C_previous in the code) as:

# Previous solution
W_const = Constant((1200))
C_previous= interpolate(W_const, V_c)

And when I assign it to the unknown c in the time loop:

C_previous.assign(c)

I get this error:

Error:   Unable to initialize vector of degrees of freedom for function.
Reason:  Cannot re-initialize a non-empty vector. Consider creating a new function.

I would like to mention again that I want to solve these two equarions over the time simultaneously and plot the concentration (unknown c) in time.
Could you please help me fix the two problems I mentioned above to get it to work? Once thanks again for your time.

---

You can probably solve the assignment problem by doing a projection.

Furthermore, a has to be bilinear (a product of test and trial functions) and to me it looks like you forgot to multiply with dt on some of the terms.

---

I have multiplied the Eq.5 by "dt" and it was reduced to Eq.6 which is a weak form of the Eq.1.
In addition there is no time derivative in the Eq.7. So "dt" does not apprear in the weak form of the Eq.2. After all, I think the weak form is correct.
In addition the projection did not solve the assignment problem. I still get the same error!