I have a solved a Poisson problem with internal discontinuity in the field where I want the flux at the discontinuity to be proportional to the jump in the field.
-\kappa \nabla u \cdot \hat{n} = R(u^+ - u^- - \alpha) where \alpha and R are positive real numbers
I have Dirichlet boundary conditions set at the left and right external boundaries, and have Neumann boundary conditions set at the top and bottom external boundaries. The solution is implemented using DG with interior penalty method.
Everything works fine when I don’t have the flux proportional to the discontinuity. The expectation is that for a larger R, and constant flux, the jump in the solution field should decrease and vice-versa.
My main question is on the dolfinx snippet that tries to implement the flux condition on the internal boundary. Is the syntax correct, and how do I include the \alpha term.
a += -R * inner(avg(grad(v)), jump(u, n)) * dS(markers['pe_se']) - inner(jump(v, n), avg(grad(u))) * dS(markers['pe_se'])
I am including the rest of the MWE that’s split into two parts for mesh creation and for the solution.
import os
import gmsh
import meshio
import numpy as np
import pyvista
import ufl
import warnings
from dolfinx import cpp, fem, io, mesh, nls, plot
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
div, dx, ds, dS, grad, inner, grad, avg, jump)
import commons, geometry, utils
warnings.simplefilter('ignore')
adaptive_refine = True
markers = {
'ne_pcc': 0, # boundary not included yet
'ne_se': 1, # left boundary
'pe_se': 2, # internal boundary
'pe_pcc': 3, # left boundary
'ne': 4, # domain not included yet
'se': 5, # left domain
'pe': 6, # right domain
"insulated": 7, # top and bottom boundary
}
points = [
(5, 0, 0),
(10, 0, 0),
(20, 0, 0),
(20, 5, 0),
(10, 5, 0),
(5, 5, 0),
]
gpoints = []
lines = []
workdir = "output/full-cell"
utils.make_dir_if_missing(workdir)
output_meshfile = os.path.join(workdir, 'mesh.msh')
tria_meshfile = os.path.join(workdir, "tria.xdmf")
line_meshfile = os.path.join(workdir, "line.xdmf")
potential_resultsfile = os.path.join(workdir, "potential.xdmf")
current_resultsfile = os.path.join(workdir, "current.xdmf")
gmsh.initialize()
gmsh.model.add('full-cell')
for idx, p in enumerate(points):
gpoints.append(
gmsh.model.occ.addPoint(*p)
)
for idx in range(0, len(points)-1):
lines.append(
gmsh.model.occ.addLine(gpoints[idx], gpoints[idx+1])
)
lines.append(
gmsh.model.occ.addLine(gpoints[-1], gpoints[0])
)
lines.append(
gmsh.model.occ.addLine(gpoints[4], gpoints[1])
)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(1, [lines[5]], markers['ne_se'])
gmsh.model.addPhysicalGroup(1, [lines[6]], markers['pe_se'])
gmsh.model.addPhysicalGroup(1, [lines[2]], markers['pe_pcc'])
gmsh.model.addPhysicalGroup(1, [lines[idx] for idx in [0, 1, 3, 4]], markers['insulated'])
gmsh.model.occ.synchronize()
se_loop = gmsh.model.occ.addCurveLoop([lines[idx] for idx in [0, 6, 4, 5]])
pe_loop = gmsh.model.occ.addCurveLoop([lines[idx] for idx in [1, 2, 3, 6]])
gmsh.model.occ.synchronize()
se_phase = gmsh.model.occ.addPlaneSurface([se_loop])
pe_phase = gmsh.model.occ.addPlaneSurface([pe_loop])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [se_phase], markers['se'])
gmsh.model.addPhysicalGroup(2, [pe_phase], markers['pe'])
gmsh.model.occ.synchronize()
# adaptive refinement
if adaptive_refine:
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [lines[idx] for idx in [6]])
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", 0.01)
gmsh.model.mesh.field.setNumber(2, "LcMax", 1)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0)
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.5)
gmsh.model.mesh.field.add("Max", 5)
gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
gmsh.model.mesh.field.setAsBackgroundMesh(5)
gmsh.model.mesh.generate(2)
gmsh.write(output_meshfile)
gmsh.finalize()
mesh_2d = meshio.read(output_meshfile)
tria_mesh = geometry.create_mesh(mesh_2d, "triangle")
meshio.write(tria_meshfile, tria_mesh)
line_mesh = geometry.create_mesh(mesh_2d, "line")
meshio.write(line_meshfile, line_mesh)
the dolfinx implementation is shown below:
import os
import gmsh
import meshio
import numpy as np
import pyvista
import ufl
import warnings
from dolfinx import cpp, fem, io, mesh, nls, plot
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
div, dx, ds, dS, grad, inner, grad, avg, jump)
import commons, geometry, utils
warnings.simplefilter('ignore')
comm = MPI.COMM_WORLD
with io.XDMFFile(comm, tria_meshfile, "r") as infile3:
domain = infile3.read_mesh(cpp.mesh.GhostMode.none, 'Grid')
ct = infile3.read_meshtags(domain, name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
with io.XDMFFile(comm, line_meshfile, "r") as infile2:
ft = infile2.read_meshtags(domain, name="Grid")
meshtags = mesh.meshtags(domain, domain.topology.dim - 1, ft.indices, ft.values)
domaintags = mesh.meshtags(domain, domain.topology.dim, ct.indices, ct.values)
ds = ufl.Measure("ds", domain=domain, subdomain_data=meshtags)
dS = ufl.Measure("dS", domain=domain, subdomain_data=meshtags)
V = fem.FunctionSpace(domain, ("DG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
n = ufl.FacetNormal(domain)
x = ufl.SpatialCoordinate(domain)
alpha = 10
gamma = 10
h = 2 * ufl.Circumradius(domain)
h_avg = avg(h)
x = SpatialCoordinate(domain)
f = fem.Constant(domain, PETSc.ScalarType(0))
g = fem.Constant(domain, PETSc.ScalarType(0))
voltage = 100e-3 # 100 mV
u_left = fem.Function(V)
with u_left.vector.localForm() as u0_loc:
u0_loc.set(0)
u_right = fem.Function(V)
with u_right.vector.localForm() as u1_loc:
u1_loc.set(voltage)
R = fem.Constant(domain, PETSc.ScalarType(8.314))
a = inner(grad(u), grad(v)) * dx - inner(n, grad(u)) * v * ds
# Add DG/IP terms
a += - inner(avg(grad(v)), jump(u, n)) * dS - inner(jump(v, n), avg(grad(u))) * dS
a += + (gamma / h_avg) * inner(jump(v, n), jump(u, n)) * dS
# internal boundary
a += -R * inner(avg(grad(v)), jump(u, n)) * dS(markers['pe_se']) - inner(jump(v, n), avg(grad(u))) * dS(markers['pe_se'])
# Add Nitsche terms
a += - inner(n, grad(v)) * u * ds(markers['ne_se']) + alpha / h * inner(u, v) * ds(markers['ne_se'])
a += - inner(n, grad(v)) * u * ds(markers['pe_pcc']) + alpha / h * inner(u, v) * ds(markers['pe_pcc'])
L = inner(f, v) * dx
L += - inner(n, grad(v)) * u_left * ds(markers['ne_se']) + alpha / h * inner(u_left, v) * ds(markers['ne_se'])
L += - inner(n, grad(v)) * u_right * ds(markers['pe_pcc']) + alpha / h * inner(u_right, v) * ds(markers['pe_pcc'])
# Neumann BC on insulated boundary
L += - (h / alpha) * inner(n, grad(v)) * g * ds(markers['insulated'])
L += + inner(g, v) * ds(markers['insulated'])
problem = fem.petsc.LinearProblem(a, L)
uh = problem.solve()
uh.name = 'potential'
pyvista.start_xvfb()
# We create a mesh consisting of the degrees of freedom for visualization
topology, cell_types, x = plot.create_vtk_mesh(V)
num_dofs_local = V.dofmap.index_map.size_local
grid = pyvista.UnstructuredGrid(topology, cell_types, V.tabulate_dof_coordinates()[:num_dofs_local])
# To make this function work in parallel, we only consider dofs owned by the current process
grid.point_data["u"] = uh.vector.array.real
grid.set_active_scalars("u")
warped = grid.warp_by_scalar("u", factor=1)
# Create plotter
plotter = pyvista.Plotter()
pyvista.start_xvfb(wait=0.05)
plotter = pyvista.Plotter()
plotter.add_mesh(warped, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("DG.png")