Solving a nonlinear System with mixed boundary conditions

Question

Hello,

I would like to solve (for a start) a small PDE-System consisting of two reaction-diffusion equations with mixed boundary conditions:
$ u_1 = g_D \text{ on } \Gamma_{D1} $
$ \frac{\partial u_1}{\partial n} = g_N \text{ on } \Gamma_{N1} $
and
$u_2 = g_D \text{ on } \Gamma_{D2} $
$\frac{\partial u_1}{\partial n} = g_N \text{ on } \Gamma_{N2}$
with $ \Omega = [0, x_{max}=12] $ being my domain.
$ \Gamma_{N1} $ and $\Gamma_{D2}$ should be applied to the left side, i.i. $ x = 0$, and $ \Gamma_{D1} $ and $\Gamma_{N2}$ to the right side $ x = x_{max}$.

I tried to define two FacetFunctions, one for each function $u_1$ and $u_2$:

exterior_facet_domains1 = FacetFunction("uint", mesh1D)
exterior_facet_domains2 = FacetFunction("uint", mesh1D)

but it's not working, I get as an error:

Found two domain data objects for same domain type 'exterior_facet', only one is allowed.

How do you define different boundary conditions for each function in such a system?

In the following code I didn't implement the Dirichlet conditions yet:

from dolfin import *
import numpy as np

# to shut off the output
set_log_level(WARNING)

nx = 120
xmax = 12
mesh1D = IntervalMesh(nx, 0, xmax)

# 0 = no plot
do_plot = 1;

# definition of solving method
V = FunctionSpace(mesh1D, 'CG', 1)
W = MixedFunctionSpace([V,V])

T = 20      
dt = 0.25

t = dt

# Define boundary segments
class left(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13   # tolerance for coordinate comparisons
        return on_boundary and x[0] < tol

class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13 
        return on_boundary and x[0] >= xmax - tol

Neumann_l1 = left()
Dirichlet_r1 = right();

Dirichlet_l2 = left();
Neumann_r2 = right();

# Create mesh function over cell facets
exterior_facet_domains1 = FacetFunction("uint", mesh1D)
exterior_facet_domains2 = FacetFunction("uint", mesh1D)

exterior_facet_domains1.set_all(1)
exterior_facet_domains2.set_all(1)

Neumann_l1.mark(exterior_facet_domains1, 2)
Neumann_r2.mark(exterior_facet_domains1, 5)

dss1 = ds[exterior_facet_domains1]
dss2 = ds[exterior_facet_domains2]


g = Expression('2.0')

# RHS
def rhs(u):
  return - u

# initial values
u1_start = Constant(0.0)
u2_start = Constant(0.0)

# Definition of the variational form
U = Function(W)
U_old = Function(W)
dU = TrialFunction(W)

(u1, u2) = split(U)
(v1, v2) = TestFunction(W)
(du1, du2) = split(dU)
(u1_old, u2_old) = split(U_old)


u1_old = interpolate(u1_start, V)
u2_old = interpolate(u2_start, V)

F1 = (u1 - u1_old)*v1*dx + dt*dot(grad(u1), grad(v1))*dx - dt*rhs(u1)*v1*dx - g*v1*dss1(2)
F2 = (u2 - u2_old)*v2*dx + dt*dot(grad(u2), grad(v2))*dx - dt*rhs(u2)*v2*dx - g*v2*dss2(5)
F = F1 + F2

dF = derivative(F,U, dU)

#>> to define: bcs!
#problem = NonlinearVariationalProblem(F, u_, bcs, J)
problem = NonlinearVariationalProblem(F, U, J = dF)
pdesys_newton = NonlinearVariationalSolver(problem)

#constants for ploting
x_points = []
for i in range(0,nx):
  x_points.append((xmax*i)/(1.0*nx)) 

x_array = np.array(x_points)

ymax = 0

solution1 = []
solution2 = []

#--------------- Solving the Equation
while t <= T:

    u1_array = u1_old.vector().array()
    u2_array = u2_old.vector().array()

    u1_nodal_values = []
    u2_nodal_values = []

    for i in range(0,nx):
      u1_nodal_values.append(u1_array[i])
      u2_nodal_values.append(u2_array[i])

    #--------------- Newton
    pdesys_newton.solve()
    u1, u2 = U.split()

    t += dt
    t_out += dt
    u1_old.assign(u1)
    u2_old.assign(u2)

Answer 1

In FEniCS you are not allowed to have two domains of same type (like you have dss1 and dss2 for type "ds")

But I think that the problem is easily solved by deleting dss2. And then in the F2 equation simply use dss1(5) instead of dss2(5). Since it is the only measure that you have set anything useful in.

Comments

Thanks, works like a charm :-)
to close this topic: by using another FacetFunctions

    exterior_facet_domains_D = FacetFunction("uint", mesh1D)

(and marking) I could define the Dirichlet boundaries

bcs = [DirichletBC(W.sub(0), g_D, exterior_facet_domains_D, 4),
       DirichletBC(W.sub(1), g_D, exterior_facet_domains_D, 5)]

so here is my complete source code for the described system. Sorry for the long post, the preview is displaying it properly as source code and thus making it much shorter.

from dolfin import *
import numpy as np

# to shut off the output
set_log_level(WARNING)

nx = 120
xmax = 12
mesh1D = IntervalMesh(nx, 0, xmax)

# 0 = no plot
do_plot = 1;

# definition of solving method
V = FunctionSpace(mesh1D, 'CG', 1)
W = MixedFunctionSpace([V,V])

T = 20      
dt = 0.25

t = dt

# ------------------- boundary conditions
# Define boundary segments
class left(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13   # tolerance for coordinate comparisons
        return on_boundary and abs(x[0]) < tol

class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13 
        return on_boundary and abs(x[0] - xmax) < tol

# instances of SubDomain
Neumann_l1 = left()
Dirichlet_r1 = right()

Dirichlet_l2 = left()
Neumann_r2 = right()

# FacetFunctions, one for Neumann and one for Dirichlet boundary conditions
exterior_facet_domains = FacetFunction("uint", mesh1D)
exterior_facet_domains_D = FacetFunction("uint", mesh1D)
exterior_facet_domains.set_all(1)
exterior_facet_domains_D.set_all(1)

# marking of the corresponding FaceFunctions
Neumann_l1.mark(exterior_facet_domains, 2)
Neumann_r2.mark(exterior_facet_domains, 3)
Dirichlet_r1.mark(exterior_facet_domains_D, 4)
Dirichlet_l2.mark(exterior_facet_domains_D, 5)

# definition of g_N and g_D, corresponding Functions
g_N = Expression('1.0')
g_D = Expression('2')

# definition of Dirichlet boundary conditions, Neumann boundaries in variational form
bcs = [DirichletBC(W.sub(0), g_D, exterior_facet_domains_D, 4),
       DirichletBC(W.sub(1), g_D, exterior_facet_domains_D, 5)]

# assignin markers to measure
dss = ds[exterior_facet_domains]
# ------------------- boundary conditions

# ------------------- problem
# RHS
def rhs(u):
  return - u

# initial values
u1_start = Constant(0.0)
u2_start = Constant(0.0)

# Definition of the variational form
U = Function(W)
U_old = Function(W)
dU = TrialFunction(W)

(u1, u2) = split(U)
(v1, v2) = TestFunction(W)
(du1, du2) = split(dU)
(u1_old, u2_old) = split(U_old)

u1_old = interpolate(u1_start, V)
u2_old = interpolate(u2_start, V)

F1 = (u1 - u1_old)*v1*dx + dt*dot(grad(u1), grad(v1))*dx - dt*rhs(u1)*v1*dx - g_N*v1*dss(2)
F2 = (u2 - u2_old)*v2*dx + dt*dot(grad(u2), grad(v2))*dx - dt*rhs(u2)*v2*dx - g_N*v2*dss(3)
F = F1 + F2

dF = derivative(F,U, dU)


problem = NonlinearVariationalProblem(F, U, bcs, J = dF)
pdesys_newton = NonlinearVariationalSolver(problem)
# ------------------- problem

#constants for ploting
x_points = []
for i in range(0,nx):
  x_points.append((xmax*i)/(1.0*nx)) 

x_array = np.array(x_points)

solution1 = []
solution2 = []

# ------------------- Solving the Equation
while t <= T:

    u1_array = u1_old.vector().array()
    u2_array = u2_old.vector().array()

    u1_nodal_values = []
    u2_nodal_values = []

    for i in range(0,nx):
      u1_nodal_values.append(u1_array[i])
      u2_nodal_values.append(u2_array[i])

    # ------------------- Newton
    pdesys_newton.solve()
    u1, u2 = U.split()

    t += dt
    t_out += dt
    u1_old.assign(u1)
    u2_old.assign(u2)