Seperated Domains within problem

[Cyph3r]

Hi all,

A couple of days ago I posted implementing internal boundaries in fenics and got some really helpful feedback. However it has made me realise I now have another issue. For now, I am considering a 2D convection-diffusion problem. However what I would like is to solve the problem in two decoupled domain simaltaneously, with a method that would allow me to add coupling to the BC’s at a later date. Right now I use no flux boundary condtions at the extenral boundaries and this gives me the expected results for the steady state solution of the problem (expontetial function which builds up at the end we convect towards). This is image 1.

I then added an internal boundary which also has the no flux condtion, and what I expected to see was the domain split in half, with the original steady state solution essentially copied onto each half. However, what I see is some weird solution which ssems somewhat unphysical (image 2).

Having tried to look around on these forums, I understand that this is because I’m using continuous elements. I get the impression to get the behaviour I want I need to either use DG elements, or use the mixed dimesnional brach of fenics. I wanted to confirm that this is right and get the opinions of others on what to do next.

To add some context, the real problem I want to solve involves an elastic medium, which is deformed by a chemical undergoing reaction-diffusion on top, and this is in turn advected by the medium. I have a version of this code working already, but want to add internal boundaries which the chemical cannot move across, (but the elastic medium can interact accros the boundaries, which is why I want to do this all as one simulation). So I’m using this problem to try and understand how to do advection-diffusion with internal boundaries. Any advice on how you would proceed would be appreciated.

Here’s my code for the simplified problem:

from fenics import *
from ufl import Index

import numpy as np
from os.path import expanduser,join

# import gmsh
# import meshio

############### mesh conversion funcs ################
def load_gmsh(path_to_mesh):
    space_dim=2
    mesh = Mesh()
    with XDMFFile(path_to_mesh + "_mesh.xdmf") as infile:
        infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 1)
    with XDMFFile(path_to_mesh + "_mf.xdmf") as infile:
        infile.read(mvc, "name_to_read")
    mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
    return mesh, mf

def create_meshio_mesh(mesh, cell_type, prune_z=True):
    import meshio
    cells = mesh.get_cells_type(cell_type)

    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

############### Make mesh ################
p=0.01
gmsh.initialize()
gmsh.model.add("square")
gmsh.model.occ.addPoint(0,0,0,p,1)
gmsh.model.occ.addPoint(1,0,0,p,2)
gmsh.model.occ.addPoint(1,1,0,p,3)
gmsh.model.occ.addPoint(0,1,0,p,4)

gmsh.model.occ.addLine(1,2,1)
gmsh.model.occ.addLine(2,3,2)
gmsh.model.occ.addLine(3,4,3)
gmsh.model.occ.addLine(4,1,4)
gmsh.model.occ.addCurveLoop([1,2,3,4],1)
gmsh.model.occ.addPlaneSurface([1],1)


gmsh.model.occ.addPoint(2,0,0,p,5)
gmsh.model.occ.addPoint(2,1,0,p,6)

# old line 2 here
gmsh.model.occ.addLine(3,6,5)
gmsh.model.occ.addLine(6,5,6)
gmsh.model.occ.addLine(5,2,7)
gmsh.model.occ.addCurveLoop([2,5,6,7],2)
gmsh.model.occ.addPlaneSurface([2],2)

gmsh.model.occ.synchronize()
gmsh.model.add_physical_group(1,[1,3,4,5,6,7],1)
gmsh.model.add_physical_group(1,[2],2)
gmsh.model.add_physical_group(2,[1],3)
gmsh.model.add_physical_group(2,[2],4)


gmsh.model.mesh.generate(2)
gmsh.write(expanduser("./mesh_files/square.msh"))
gmsh.fltk.run()
gmsh.finalize()

msh = meshio.read(expanduser("./mesh_files/square.msh"))
line_mesh = create_meshio_mesh(msh, "line", prune_z=True)
triangle_mesh = create_meshio_mesh(msh, "triangle", prune_z=True)

meshio.write(join(expanduser("./mesh_files/square_mf.xdmf")), line_mesh)
meshio.write(join(expanduser("./mesh_files/square_mesh.xdmf")), triangle_mesh)

############### Run simulation ################

set_log_level(30)
parameters["ghost_mode"] = "shared_facet"
T=10
dt=10**-2
num_steps=int(T/dt)

mesh,mf=load_gmsh(expanduser("./mesh_files/square"))

P1=FiniteElement('P',mesh.ufl_cell(),1)
func_space = FunctionSpace(mesh,P1)

rho=Function(func_space)
rho_prev=Function(func_space)
test_rho=TestFunction(func_space)

rho_0=interpolate(Constant(1),func_space)
assign(rho_prev,rho_0)
assign(rho,rho_0)

out_file1 = File(expanduser("./con_diff/output_rho.pvd"))

n=FacetNormal(mesh)
dx=Measure('dx',domain=mesh,subdomain_data=mf)
ds_ext=Measure('ds',domain=mesh,subdomain_data=mf,subdomain_id=1)
ds_int=Measure('dS',domain=mesh,subdomain_data=mf,subdomain_id=2)
j=Index()
v=Expression(('1','0'),degree=2,domain=mesh)
D=Constant(1)
dt=Constant(dt)

F =(
    (rho-rho_prev)*test_rho*dx + dt*Dx(rho*v[j],j)*test_rho*dx
    + dt*D*Dx(rho,j)*Dx(test_rho,j)*(dx)
     - dt*rho*v[j]*n[j]*test_rho*ds_ext - dt('+')*rho('+')*v[j]('+')*n[j]('+')*test_rho('+')*ds_int
    )

for step in range(num_steps):
    t=step*dt
    solve(F==0,rho)
    rho_prev.assign(rho)
    out_file1 << (rho_prev,t)
    print(t)

[IngeborgGjerde]

Hi,

In order to get a function that is discontinuous across the interface you can combine two functions that are P1 on each submesh

From your code I see that you are assembling and solving the problem in the form F==0. Is this because you want to solve a non-linear problem eventually? The reason I ask is that there is a bug causing non-linear problems to fail with fenics-mixed-dim.

We can fix that by switching to fenics_ii, with the caveat that the syntax is a bit more involved.

sketch of fenics_ii code

from fenics import *
from ufl import Index
set_log_level(30)
parameters["ghost_mode"] = "shared_facet"
from xii import *

# Make 2d mesh
mesh = UnitSquareMesh(30,30)


marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]<0.5

bound = Boundary()
bound.mark(marker, 1)
meshleft = SubMesh(mesh, marker, 1)
meshright = SubMesh(mesh, marker, 0)


# Run simulation
T=1
dt=0.1
num_steps=int(T/dt)

W = [FunctionSpace(mesh, 'CG', 1) for mesh in (meshleft, meshright)]

rho = ii_Function(W)
rhol, rhor = rho

rho_prev = ii_Function(W)
rhol_prev, rhor_prev = rho_prev

test_rho = list(map(TestFunction, W))

assign(rhol_prev,interpolate(Constant(1),W[0]))
assign(rhor_prev,interpolate(Constant(1),W[1]))

j=Index()
v=Expression(('1','0'),degree=2)
D=Constant(1)
dt=Constant(dt)

F = [0,0]
for i in range(len(F)):

    rhoi = rho[i]
    rhoi_prev = rho_prev[i]
    test_rhoi = test_rho[i]
    vi = interpolate(v, VectorFunctionSpace(W[i].mesh(), 'CG', 1))

    dxi = Measure('dx', domain=W[i].mesh())
    dsi = Measure('ds', domain=W[i].mesh())
    n = FacetNormal(W[i].mesh())

    # note: I changed some of the syntax here to fit fenics_ii, please double check
    F[i] +=(
        (rhoi-rhoi_prev)*test_rhoi*dxi + dt*dot(grad(rhoi),vi)*test_rhoi*dxi
        + dt*dot(grad(rhoi),vi)*test_rhoi*dxi
        + dt*D*rhoi.dx(j)*test_rhoi.dx(j)*dxi
        - dt*rhoi*dot(vi,n)*test_rhoi*dsi
        )

    
from xii.nonlin.jacobian import block_jacobian
# Solution
dF = block_jacobian(F, rho)

# Newton
eps = 1.0
tol = 1.0E-10
niter = 0
maxiter = 125

drho = ii_Function(W)

for step in range(num_steps):
    t=step*dt
    print('timestep: ',step)
    
    # newton solver at this time step
    while eps > tol and niter < maxiter:
        niter += 1

        A, b = list(map(ii_assemble, (dF, F)))
        A, b = list(map(ii_convert, (A, b)))

        solve(A, drho.vector(), b)

        eps = sqrt(sum(x.norm('l2')**2 for x in drho.vectors()))

        for i in range(len(W)):
            rho[i].vector().axpy(-1, drho[i].vector())

    rhol.assign(rho[0])
    rhor.assign(rho[1])
    
File('shared/mixed-dim-comparisons/rhol.pvd') << rhol
File('shared/mixed-dim-comparisons/rhor.pvd') << rhor

Hope that helps, feel free to ask more questions and I will try my best to help

[Cyph3r]

Hey, thanks so much for your reply, this looks really helpful!

To answer your question in the end yes the problem will be highly non-linear so I’ll absolutely try fenics_ii.

I’m working on something else temporarily but will come back to this post when I go back to this problem. Thanks again for your help!

[Cyph3r]

Hey,

So I am finally back to this problem and the code you provided is really helpful. Thank you!

However my overall goal is to eventually have the velocity field v be defined by another coupled field which has its own evolution in time. However I’m not sure how to solve this problem in fenics_ii as I want this second field to be continuous across the whole mesh, unlike the density field I’m solving for here. I wondered if you would be able to point me in the right direction here? Or link some examples of similar problems? (IE a field which has an internal boundary coupled to a field which doesn’t have this constraint.)