Solving on subdomains with robin condition between subdomains

Question

Hello,

I'm trying to solve the heat equation on part of my domain ($ \Omega_1 $) with a RobinBC at the boundary of the subdomain ($ x=l $). On the other part ($ \Omega_2 $), i want to set the solution to a fixed value ($ u_0 $) this is needed as initial condition for further time steps). I plan to do this via DirichletBC. How do i implement the interior condition (find the facets and use dS measure)? Do i need DG methods for this?

Similar question

Answer 1

Generally speaking, Robin BCs are for conforming discretizations of second-order elliptic PDEs treated as natural condition, i.e. adding an appropriate surface integral into weak formulation. There would be no usage of DirichletBC.

But it is not clear from your formulation what do you want. Should the solution be continuous at the interface? If yes, then the paragraph above applies.

find the facets

It depends how your program looks like - what should be input for searching for facets. Please, start with some minimal example (at least mesh) and state clearly what do you want to achieve.

Comments

state clearly what do you want to achieve

Solve the PDE on a subdomain independently of the other subdomain (where the solution may be constant) with a robin boundary condition for the left solution at the interface.

adding an appropriate surface integral into weak formulation.

I have a working solution for one of the subdomains. I think die RobinBC is working the way it's supposed to work but if i solve on the other subdomain too, the solution is different.

My idea is to set the values in the right subdomain manually (I only need it as initial condition). I found the command ident_zeros(), which basically does what i want without touching the solution on the left side.

The problem is, i want to use dolfin_adjoint and i don't think it's able to annotate manually set functions or ident_zeros(), so i try to use DirichletBC. But how can i do this without changing my solution on the left subdomain? And there can be a discontinuity at the interface (although it would be possible to work around that with a rapidly decreasing function)

Here is a running example.

from dolfin import * 

## Mesh and function space
mesh = UnitIntervalMesh(100) 
V = FunctionSpace(mesh, 'CG', 2) 

u0 = Constant(10.0)

## Sourceterm and constant for RobinBC
f = Constant(10.0)
alpha = Constant(1.0) 

## Define domains
class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < 0.5 

class Omega2(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 0.5 

## Define boundary between subdomains for selection boundary facet
class Grenze(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.5) 

## Set DirichletBC
class Boundary(SubDomain):
  def inside(self, x, on_boundary):
                return x[0]<DOLFIN_EPS 
bc = DirichletBC(V, u0, Boundary())#, 'pointwise')

u = TrialFunction(V)
v = TestFunction(V)


## Cell and Facetfunctions
subdom = CellFunction("size_t", mesh)
omega1 = Omega1()
omega2 = Omega2()

omega1.mark(subdom, 0)
omega2.mark(subdom, 1)

grenze = Grenze()
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
grenze.mark(boundaries, 1)

## Define Measures
dx = Measure("dx")[subdom]
dSs = Measure("dS")[boundaries]

## Define form
a = dot(grad(v), grad(u))*dx(0) + alpha*avg(u)*avg(v)*dSs(1)
L = f*v*dx(0) + alpha*u0*avg(v)*dSs(1)

## Assemble
A, b = assemble_system(a, L , bc) 
#A.ident_zeros() # This gives a solution with zeros on the unsolved subdomain

## Solving
u = Function(V)   
solve(A, u.vector(), b)

## Plots
plot(u)

plot(boundaries)
interactive()