Drift-diffusion problem involving submeshes

variational formulation
1

Hi,

I have a problem similar to the one discussed below except that I need to adapt methods from the legacy version to my conda installation of 0.4.1:

To sum up: I have Poisson equation in the whole domain and diffusion equations active only in certain subdomains. The diffusion equation determines charge concentration which feeds as source into Poisson equation.
Likewise, Poisson equation calculates potential that feeds into diffusion equation’s advection term.

I could not figure how to express the source term for poisson equation which depends upon diffused charges inside the submesh but is zero outside. Further, I get error in trying to define the nonlinear problem which suggests I have mistakes even where I do not get explicit errors.

Below is an MWE of how I am trying to do it. I will appreciate if someone could direct me to relevant tutorials or point out mistakes in my attempt specially towards the end where I just cannot proceed due to errors:

import dolfinx, ufl, petsc4py
import mpi4py
import numpy as np

#create the geometry
A = 1 # side length

#define a simple field of view
mesh = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, A, A, dolfinx.mesh.CellType.quadrilateral)

#left half is a submesh
submesh_entities = dolfinx.mesh.locate_entities(mesh, dim=2, marker=lambda x: (x[0] < 0.5))
submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, submesh_entities)

#function spaces, taken from How to correctly assign values to mixed elements in FEniCS-X?
Omega1_FE = ufl.FiniteElement(“Lagrange”, mesh.ufl_cell(), 2) # complete mesh
Omega2_p_FE = ufl.FiniteElement(“Lagrange”, submesh.ufl_cell(), 2) # submesh for positive ions
Omega1 = dolfinx.fem.FunctionSpace(mesh, Omega1_FE) # for poisson equation/whole mesh
Omega2_p = dolfinx.fem.FunctionSpace(mesh, Omega2_p_FE) # for diffusion equation/submesh
ME = ufl.MixedElement([Omega1_FE, Omega2_p_FE])
Omega = dolfinx.fem.FunctionSpace(mesh, ME)

#trial and test functions
u = dolfinx.fem.Function(Omega)
V, cp = ufl.split(u)
p, q = ufl.TestFunctions(Omega)

#define teh dielectric constant for poisson equation
epsr = [1, 4] # values of dielectric constant in left and right half of the mesh
kappa = dolfinx.fem.Function(Omega1) # dielectric constat, defined on the whole mesh
cells_inside_submesh = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, lambda x: x[0] <= 0.5)
cells_outside_submesh = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, lambda x: x[0] >= 0.5)
kappa.x.array[cells_inside_submesh] = epsr[0]
kappa.x.array[cells_outside_submesh] = epsr[1]

#define boundary condition
#left boundary, at 10V
dofs = dolfinx.fem.locate_dofs_geometrical(Omega1, lambda x: np.isclose(x[0], 0))
bcs = [dolfinx.fem.dirichletbc(petsc4py.PETSc.ScalarType(10), dofs, Omega1)]
#right boundary, at 0V
dofs = dolfinx.fem.locate_dofs_geometrical(Omega1, lambda x: np.isclose(x[0], 1))
bcs = [dolfinx.fem.dirichletbc(petsc4py.PETSc.ScalarType(0), dofs, Omega1)]

#finally define the weak form
#first define the current terms of positive charges
Jp = -ufl.grad(cp) - cp * ufl.grad(V) # current of positive ions

aV = kappa * ufl.dot(ufl.grad(V), ufl.grad(p)) * ufl.dx
acp = - ufl.dot(Jp, ufl.grad(q)) * ufl.dx

LV = cp * ufl.dx
Lcp = dolfinx.fem.Constant(submesh, petsc4py.PETSc.ScalarType(0.0)) * ufl.dx

F = aV + acp - LV - Lcp

problem = dolfinx.fem.petsc.NonlinearProblem(F, u=u, bcs = bcs) # this step runs into error