Solving diffusion equation on two domains with different D applying interface permeability condition

General
1

Hello all,

I hope someone can help me out with the following problem I am having. I am recently quite new to finite elements and Fenics is my first open-source solver. I have been looking through tutorials and done a lot of research online.

My aim is: Solve the diffusion equation D div(grad(u)) = du/dt on a domain that has two different D values (D1,D2). There is a large 2D squared domain (we will call it DOMAIN 1) and then a lot of polygons inside this 2D squared domain which we will refer to DOMAIN 2. DOMAIN 1 has diffusivitiy value D1 and DOMAIN 2 has diffusviity value D2.

My code works fine solving the temporal diffusion equation without any interface boundary condition between the two domains. Here is a brief of the code.

I want to add now a interface boundary condition that should set D dU/dX(‘+) = kappa*(u(’+‘) - u(’-')) . The interface boundary should set the flux based on the jump variation in the normal direction.

I am imposing the interface condition in the variational form but seems to not be working. Apart from that, it seems that the Dirichlet condition is also not applied. I have checked my boundary markers plotting and they seem fine. Does anyone have any idea?

V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

D_intra = Constant(1.)
D_extra = Constant(2.5)
diffusivity = Function(V)
diffusivity_expression = DiffusivityExpression(triangle_domains, 1., 2.5, degree=1)
diffusivity.interpolate(diffusivity_expression)
bc = DirichletBC(V, Constant(0), boundary_markers, 4)

u0 = Expression('x[0] >= -2. && x[0] <= 2. && x[1] >= -2. && x[1] <=2. ? 1.0 : 0.0', degree=1)
u_n = Function(V)
u_n.interpolate(u0)
# Normalize so that the integral over the domain is 1
u_n.vector()[:] /= assemble(u_n*dx)

# Define time-stepping parameters
T = 500
dt = Constant(.25)
f = Constant(0.)
kappa = Constant(0.01)
n = FacetNormal(mesh)
n_2d = as_vector([n[0], n[1]])  # Project normal onto z=0 plane

# Define variational form
dS_internal = Measure('dS', domain=mesh, subdomain_data=facet_domains, subdomain_id=6)
D = Constant(2.)

un_plus_normal = dot(grad(u('+')), n_2d('+'))
un_minus_normal = dot(grad(u('-')), n_2d('-'))

# Incorporate spatially varying diffusivity and boundary conditions
F = dt * (diffusivity * dot(grad(u), grad(v)) * dx - kappa * (un_plus_normal - un_minus_normal) * jump(v) * dS_internal) + u * v * dx - u_n * v * dx
a,L = lhs(F), rhs(F)

output_file = XDMFFile(comm, "output.xdmf")
output_file.parameters["flush_output"] = True
output_file.parameters["functions_share_mesh"] = True
output_file.write(u_n, 0.0)

# Define linear solver
solver = PETScKrylovSolver()
#solver.parameters.update(solver_parameters)
# Time loop
t = 0
step = 0
u = Function(V)
dtss=0.25
while t <= T:
    t += dtss
    step += 1
    A, b = assemble_system(a, L, bc)
    # Solve linear system
    solver.solve(A, u.vector(), b)

    # Update solution at previous time step
    u_n.assign(u)

    # Write solution to file
    output_file.write(u_n, t)

With class

class BoundaryConditionExpressionBetweenDomains(UserExpression):
    def __init__(self, cell_domains, kappa, **kwargs):
        super().__init__(**kwargs)
        self.cell_domains = cell_domains
        self.kappa = kappa
    def eval_cell(self, values, x, cell):
        if self.cell_domains[cell.index] == 1:
            values[0] = self.kappa * (1/2.5)  # your boundary condition for subdomain 1
        elif self.cell_domains[cell.index] == 2:
            values[0] = self.kappa * (1) # your boundary condition for subdomain 2
        else:
            raise RuntimeError("This cell does not have a subdomain attributed to it")
2

By using a single function space to solve the problem in both domains, you cannot apply such a condition with Lagrange elements, as u_one_side=u_other_side.
You would either have to use:

  1. Discontinuous galerkin method and apply different terms on different interior boundaries (most places enforce continuity)
  2. Use software such as fenics_ii, multiphenics or the mixed function space in fenics (C. Daversin has a paper on this: https://dl.acm.org/doi/abs/10.1145/3471138) to couple two sub-domains.
1 Like
3

Hello, thank you for the fast reply. I see, I tried using ‘DG’ but then my code crashes. What do you mean by applying different terms on different interior boundaries? In my case I do have flux continuity so D(‘+’) dU/dX(‘+) = D(‘-’) dU/dX(‘-).

How could I use two dS one for the interior and one for the exterior if they both share the same edge?

4

You cannot just change the element to DG, and expect the code to work. To use DG you would need to change your variational form to something similar to Nitsches method. See for instance: Bitbucket
If you want to go in that direction, I strongly recommend you to look at some introductory textbooks on the subject.

5

The DOFs in DG are local to each element, there is no continuity across elements,except if you enforce it with a coupling term (Interior penalty/Nitsche based method).

6

I think I will use the mixedFunctionSpace with two subdomains so that I can enforce the boundary condition. Thank you for the resource!