Hello,
I would like to use fenics to obtain the potential in a quantum dot structure. Below is the code to generate the layered quantum dot structure by pygmsh.
import numpy as np
import pygmsh
import meshio
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_min = 0
geom.characteristic_length_max = 5
g=2.0
sio2 = geom.add_box([0.0, 0.0, 0.0], [70, 70, 20])
hbn = geom.add_box([0.0, 0.0, 20.0], [70, 70, 5])
gate1 = geom.add_box([0.0, 30.0, 25.0], [30-g, 10, 5])
gate2 = geom.add_box([30.0, 40.0+g, 25.0], [10, 30-g, 5])
gate3 = geom.add_box([40.0+g, 30.0, 25.0], [30-g, 10, 5])
gate4 = geom.add_box([30.0, 0.0, 25.0], [10, 30-g, 5])
group1=[sio2,gate1,gate2,gate3,gate4]
x=np.linspace(0,70,281)
points=[[xi,35,20] for xi in x]
group2=[hbn]
for i in range(281):
group2.append(geom.add_point(points[i], mesh_size=0.1))
geom.boolean_fragments(group1, group2)
mesh = geom.generate_mesh(3)
for cell in mesh.cells:
if cell.type == "triangle":
triangle_cells = cell.data
elif cell.type == "tetra":
tetra_cells = cell.data
meshio.write("mesh3d_for_potential.xdmf", tetra_mesh)
In the above code, I also embed a line below into the mesh which is located at the interface of SiO2 and hBN (z=20). Then I would like to obtain the potential on the line.
x=np.linspace(0,70,281)
points=[[xi,35,20] for xi in x]
The code for solving the Poisson equation in the quantum dot structure is shown below.
from fenics import *
from mshr import *
import numpy as np
mesh = Mesh(MPI.comm_self)
with XDMFFile(MPI.comm_self, "mesh3d_for_potential.xdmf") as infile:
infile.read(mesh)
vtkfile=File('mesh_3D.pvd')
vtkfile<<mesh
V=FunctionSpace(mesh,'P',1)
g=2
tol=1e-14
class Gate1(SubDomain):
def inside(self,x,on_boundary):
return between(x[0],(0,30-g)) and between(x[1],(30,40)) and between(x[2],(25-tol,30+tol))
class Gate2(SubDomain):
def inside(self,x,on_boundary):
return between(x[0],(30,40)) and between(x[1],(40+g,70)) and between(x[2],(25-tol,30+tol))
class Gate3(SubDomain):
def inside(self,x,on_boundary):
return between(x[0],(40+g,70)) and between(x[1],(30,40)) and between(x[2],(25-tol,30+tol))
class Gate4(SubDomain):
def inside(self,x,on_boundary):
return between(x[0],(30,40)) and between(x[1],(0,30-g)) and between(x[2],(25-tol,30+tol))
class SiO2(SubDomain):
def inside(self, x, on_boundary):
return between(x[2],(0,20+tol))
class hBN(SubDomain):
def inside(self, x, on_boundary):
return between(x[2],(20-tol,25+tol))
def boundary(x,on_boundary):
return on_boundary and near(x[2],0,tol)
gate1=Gate1()
gate2=Gate2()
gate3=Gate3()
gate4=Gate4()
hbn=hBN()
sio2=SiO2()
bc1=DirichletBC(V,-1.5,gate1)
bc2=DirichletBC(V,-1.5,gate2)
bc3=DirichletBC(V,-1.5,gate3)
bc4=DirichletBC(V,-1.5,gate4)
bc5=DirichletBC(V,0.0,boundary)
bcs=[bc1,bc2,bc3,bc4,bc5]
materials = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure('dx', domain=mesh, subdomain_data=materials)
gate1.mark(materials,2)
gate2.mark(materials,2)
gate3.mark(materials,2)
gate4.mark(materials,2)
hbn.mark(materials,1)
sio2.mark(materials,0)
vtkfile=File('materials.pvd')
vtkfile<<materials
class Permeability(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(**kwargs)
self.materials = materials
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = 3.9 # sio2
elif self.materials[cell.index] == 1:
values[0] = 5.1 # hBN
def value_shape(self):
return ()
er = Permeability(materials, degree=2)
u=TrialFunction(V)
v=TestFunction(V)
f=Constant(0.0)
a=er*inner(grad(u),grad(v))*dx
L=f*v*dx
u=Function(V)
solve(a==L,u,bcs)
vtkfile=File('TMDC_4gate_2.pvd')
vtkfile<<u
# output potential profiles
x=np.linspace(0+tol,70-tol,281)
points=[(xi,35,20) for xi in x]
potential_line=np.array([[point[0],u(point)] for point in points])
np.savetxt('potential_line.txt',potential_line)
After solving the Poisson equation, the potential in the line is obtained and saved as “potential_line.txt”. Below is the result.
The potential in the line is not symmetric. However, because the quantum dot structure is symmetric, the potential in the line is supposed to be symmetric.
I have embedded the line into the mesh and used a denser mesh along the line. I have no idea on how to obtain higher accuracy in the balance of computational time and mesh size? Appreciate any idea on this issue.
Thank you very much!