2 unknowns solved in 2 different subdomains

Question

Hi everyone,
I have a problem with 2 unknowns (C_0 and C_1) and a domain from x = 0 to 200. I want to split the domain into two subdomains and solve only one of the variables on each:
C_0 on omega0 -> x=0 to 100
C_1 on omega1 -> x = 100 to 200
Ideally, I don't even want C_0 defined on omega1, nor C_1 on omega0.

I read through http://fenicsproject.org/documentation/tutorial/materials.html
and
http://fenicsproject.org/qa/544/solving-pde-on-submesh

and came up with the following code:

from dolfin import *

code = '''
class MyFunc : public Expression
{
public:

  MyFunc() : Expression()
  {}

double C_0_0, C_1_0;
void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = C_0_0;
    values[1] = C_1_0;
  }
std::size_t value_rank() const
    { return 1;    }
std::size_t value_dimension(std::size_t i) const
    { return 2; }  

};'''


set_log_level(PROGRESS)

mesh = IntervalMesh(40, 0.0, 200)

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] <= 100 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] >= 100 else False

# Mark subdomains with numbers 0 and 1
subdomain0 = Omega0()
subdomain1 = Omega1()

subdomains = CellFunction('size_t', mesh)

subdomain0.mark(subdomains, 1)
subdomain1.mark(subdomains, 2) 

V1 = FunctionSpace(mesh, "Lagrange",1)
V = MixedFunctionSpace([V1, V1])

dU = TrialFunction(V)

test_C_0, test_C_1 = TestFunctions(V)
U = Function(V)
fcns = split(U)
C_0 = fcns[0]
C_1 = fcns[1]



##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Domain marking ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~##
left, right, interface= compile_subdomains([
    "(std::abs( x[0])                    < DOLFIN_EPS) && on_boundary",
    "(std::abs( x[0] - width )           < DOLFIN_EPS) && on_boundary",
    "(std::abs( x[0] - location ) < DOLFIN_EPS)"])
right.width = 200
interface.location = 100

dx = Measure('dx')[subdomains]

bcs = [

    DirichletBC(V.sub(0), .5, left),
    DirichletBC(V.sub(1), .3, left),

    DirichletBC(V.sub(0), .4, interface),
    DirichletBC(V.sub(1), .4, interface),

    DirichletBC(V.sub(0), .2, right),
    DirichletBC(V.sub(1), .4, right),

    ]

F_C_0 = inner(grad(C_0), grad(test_C_0) )*( dx(1))
F_C_1 = inner(grad(C_1), grad(test_C_1) )*( dx(2) )

F = F_C_0 + F_C_1


# ~~~~~~~~~~~~~~~~ Solver ~~~~~~~~~~~~~~~~~~~ #


f = Expression(code,
    C_0_0 = .5, C_1_0 = .3)
U.interpolate(f)

Jac = derivative(F, U, dU)
problem = NonlinearVariationalProblem(F,U, bcs, Jac)
solver = NonlinearVariationalSolver(problem)

solver.solve()

plot(U.split()[0])
plot(U.split()[1])
interactive()

It seems to work (in that the plots show separate ranges for C_0 and C_1), but I think, judging from the PETSc output, that both variables are being defined on the whole mesh and they are just set to zero or something on the other subdomain.

Is there a better way to do this? I tried a subMesh model based on the second hyperlink above but got an error about some data structure being defined twice for the cell level or something like that.
Thanks!

Comments

Here is my code with the submesh concept in case it is useful:

 from dolfin import *


code = '''
class MyFunc : public Expression
{
public:

  MyFunc() : Expression()
  {}

double C_0_0, C_1_0;
void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = C_0_0;
    values[1] = C_1_0;
  }
std::size_t value_rank() const
    { return 1;    }
std::size_t value_dimension(std::size_t i) const
    { return 2; }  

};'''


mesh = IntervalMesh(40, 0.0, 200)

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] <= 100 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] >= 100 else False

subdomain0 = Omega0()
subdomain1 = Omega1()

subdomains = CellFunction('size_t', mesh)

subdomain0.mark(subdomains, 0)
subdomain1.mark(subdomains, 1) 

submesh_0 = SubMesh(mesh, subdomains, 0)
submesh_0_domain = CellFunction('size_t', submesh_0)
submesh_0_domain.set_all(0)
dx_submesh_0 = Measure('dx')[submesh_0_domain]

V1_submesh = FunctionSpace(submesh_0, "Lagrange", 1)
V1 = FunctionSpace(mesh, "Lagrange",1)
V = MixedFunctionSpace([V1_submesh, V1])

dU = TrialFunction(V)

test_C_0, test_C_1 = TestFunctions(V)
U = Function(V)
fcns = split(U)
C_0 = fcns[0]
C_1 = fcns[1]



##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Domain marking ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~##
left, right, interface= compile_subdomains([
    "(std::abs( x[0])                    < DOLFIN_EPS) && on_boundary",
    "(std::abs( x[0] - width )           < DOLFIN_EPS) && on_boundary",
    "(std::abs( x[0] - location ) < DOLFIN_EPS)"])
right.width = 200
interface.location = 100

dx = Measure('dx')[subdomains]


bcs = [

    DirichletBC(V.sub(0), .5, left),
    DirichletBC(V.sub(1), .3, left),

    DirichletBC(V.sub(0), .4, interface),

    DirichletBC(V.sub(0), .2, right),
    #DirichletBC(V.sub(1), .4, right),

    ]

F_C_0 = inner(grad(C_0), grad(test_C_0) )*( dx_submesh_0(0) )
F_C_1 = inner(grad(C_1), grad(test_C_1) )*( dx(2) )

F = F_C_0 + F_C_1


# ~~~~~~~~~~~~~~~~ Solver ~~~~~~~~~~~~~~~~~~~ #


f = Expression(code,
    C_0_0 = .5, C_1_0 = .3)
U.interpolate(f)

Jac = derivative(F, U, dU)
problem = NonlinearVariationalProblem(F,U, bcs, Jac)
solver = NonlinearVariationalSolver(problem)

solver.solve()

plot(U.split()[0])
plot(U.split()[1])
interactive()

Which dies with the message: 'ufl.log.UFLException: Found two domain data objects for same domain type 'cell', only one is allowed.'

---

Update:

I found this link too:
https://answers.launchpad.net/fenics/+question/188087

Somebody please correct me if I'm wrong, but my current impression is that there isn't a neat way to use completely separated domains (not yet anyway), and so the problem comes down to what to do with the variables on the domains for which they are not defined. (C_0 on omega1 for example). From playing around, it looks like some issues associated with, viz empty lines in the system matrix as described in the link above, have already been addressed.

My actual problem is non-linear and contains a logarithm of C_0 and C_1 in the flux term. I found these terms were causing problem which were alleviated by using the snes solver instead of the NonlinearVariationalSolver, which I guess includes some techniques by default to avoid error of this sort.

---

Hey, did you continue working on this problem? I want to do something similar and, as you pointed out, defining and solving the two equations over the whole domain would be quite some additional work...

---

Hi Bobby,

Haven't looked at it in a while, but no, I don't think I found a better way. As I remember I started running into huge ill-conditioning problems (makes sense), so I had to use the PETSc solver with these options:

--petsc.ksp_type preonly --petsc.pc_type lu --petsc.pc_factor_shift_type NONZERO

the first two of which were to get it to use a direct solver, and the last one to make the entries that you are not defining anything for non-zero.
Good luck and let me know if you come up with something better!

Answer 1

The usual approach is to restrict fields by applying DirichletBC to spare DOFs. This is, of course, inefficient (but working) workaround.

You may also like to check experimental Restriction functionality. See restriction demo in development version. But I'm affraid that this does not alllow you to mix different domains in one form and will raise

ufl.log.UFLException: Found two domain data objects for same domain type 'cell', only one is allowed.

or something related. This will be fixed in a future but requires vast amount of work.

Comments

Thanks Jan, I tried fixing the other DOFs but found that in doing so, continuity is imposed in the variable between the desired and undesired domains, implying the dirichlet condition also applies at the boundary (which is definitely not desired!).

---

This is a tricky business sometimes. You may utilize

 DirichetBC(..., points_of_restricted_dofs, method="pointwise")

and you need to search for points to be restricted and create SubDomain points_of_restricted_dofs which can be very tricky thing. But there are tools in DOLFIN to do this.