Diffusion solving in two material domains

General
1

Hi everyone, I am solving the diffusion equation and the model contains two materials with different diffusion coefficients. However, it seems that the diffusion coefficients don’t work after defining them as below. I guess it should be a function space problem (diffusion coefficients and solution function use two different function spaces), do you have any better suggestions, thanks!

D_g_0 = 8.22e-8  
E_g   = 0.5  
D_gb_0 = 8.36e-4
E_gb  = 0.06  
k = 8.6e-5  
T = 300.0  
mesh = Mesh()
with XDMFFile("mesh/grain_model_triangles.xdmf") as xdmf:
    xdmf.read(mesh)

mvc_surface = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh/grain_model_triangles.xdmf") as xdmf:
    xdmf.read(mvc_surface, "gmsh:physical")
facet_markers = MeshFunction("size_t", mesh, mvc_surface)

mvc_line = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("mesh/grain_model_lines.xdmf") as xdmf:
    xdmf.read(mvc_line, "gmsh:physical")
line_markers = MeshFunction("size_t", mesh, mvc_line)

P = FunctionSpace(mesh, "DG", 0)  # Discontinuous Galerkin
D_diffusion = Function(P)

# Define conductivity values
kappa_grain = D_g_0 * exp(-E_g / (k * T))
kappa_gb = D_gb_0 * exp(-E_gb / (k * T))

D_diffusion.vector().set_local(np.zeros(D_diffusion.vector().local_size()))

for cell in cells(mesh):
    if facet_markers[cell] == 1:  # Assuming physical group 1 corresponds to value 1
        D_diffusion.vector()[cell.index()] = kappa_grain
    elif facet_markers[cell] == 2:  # Assuming physical group 2 corresponds to value 2
        D_diffusion.vector()[cell.index()] = kappa_gb

with XDMFFile("D_diffusion.xdmf") as file:
    file.write(D_diffusion)

D_left = 1e2  # Example diffusion coefficient for left boundary
D_right = 1e0  # Example diffusion coefficient for right boundary

V = FunctionSpace(mesh, "DG", 1)

bc_left = DirichletBC(V, Constant(D_left), line_markers, 4)  # Left boundary (physical group 4)
bc_right = DirichletBC(V, Constant(D_right), line_markers, 6)  # Right boundary (physical group 6)

u = TrialFunction(V)
v = TestFunction(V)
a = D_diffusion*dot(grad(u), grad(v)) * dx  # Diffusion operator
L = Constant(0) * v * dx  # Assuming no source term

u = Function(V)
solve(a == L, u, [bc_left, bc_right])
with XDMFFile("solution.xdmf") as file:
    file.write(u)
2

Your code is not MWE, as it is missing the includes. This also means it is difficult to ascertain which version you are running. And you’re not saying what “doesn’t work”.

The ‘poisson with piecewise diffusivity’ is a demo in nearly all Fenics versions (old and new), so you’d just have to look for the right one.

For FenicsX:
https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html

1 Like