Incorrect solution for Poisson equation

variational formulation
1

I’m solving a 3D static electric field in a box of 1000*1000*50 um^3. The boundary between two areas with different charge density somehow shifted.

How could I restate the formal problem to get a correct solution?

The problem is

\nabla^2U(\vec{r})=-\rho/\varepsilon

in which

\rho =\left\{ \begin{aligned} e*4*10^{17}\mathrm{cm}^{-3},0<z<1,\\ e*1*10^{13}\mathrm{cm}^{-3},1<z<50 \end{aligned} \right. , \varepsilon = 11.7 \varepsilon_0,U(z=0)=-1000V,U(z=50)=0

Python codes are below.

import fenics
import mshr

mesh = mshr.generate_mesh(mshr.Box(fenics.Point(0, 0, 0), fenics.Point(1000, 1000, 50)),32)
V = fenics.FunctionSpace(mesh, 'P', 1)
u_D = fenics.Expression('x[2]<tol ? p_1:p_2', degree=2, tol=1E-14, p_1=-2000, p_2=0)
def boundary(x, on_boundary):
    return abs(x[2])<1E-14 or abs(x[2]-50)<1E-14
bc_l = fenics.DirichletBC(V, u_D, boundary)

e0 = 1.60217733e-19
perm_mat = 9.76
perm0 = 8.854187817e-12
f_1 = e0*4e5*1e6/perm0/perm_mat
f_2 = e0*10*1e6/perm0/perm_mat

u = fenics.TrialFunction(V)
v = fenics.TestFunction(V)
f = fenics.Expression('x[2] < 1 + tol ? f_1 : f_2',
    degree=1,width=1,f_1=f_1,f_2=f_2, tol = 1E-14)
a = fenics.dot(fenics.grad(u), fenics.grad(v))*fenics.dx
L = f*v*fenics.dx
U = fenics.Function(V)
fenics.solve(a == L, U, bc_l, solver_parameters=dict(linear_solver='gmres', preconditioner='ilu'))

W = fenics.VectorFunctionSpace(mesh, 'P', 1)
E_field = fenics.project(fenics.as_vector((U.dx(0), U.dx(1), U.dx(2))),W)
2

lgad3D
lgad3D1440×1404 108 KB

The returned solution in z direction are in the upright corner of the given picture.

3

Your code is not complete (it cannot be executed by others).
Also it is lacking markdown formatting, i.e.

```python
def test():
    assert True
```

It would also help to have the mathematical description of the problem you are solving.

1 Like
4

Thank you for you remind and I’ve updated the main post.

5

First of all, you are using an unstructured grid (since you are using mshr to generate the mesh).

This means that the internal boundary you are specifying

is not aligned with the mesh facets.

Also, You should not use a continuous finite element (which is the behavior when using degree=1)
https://bitbucket.org/fenics-project/dolfin/src/74a6ac2095059ea6582fe4ea11686c1c723e14f5/python/dolfin/function/expression.py#lines-22:34
You should instead send in degree=0 to get a DG 0 interpolation of the expression.

2 Likes
6

Thank you, and how could I define a structured grid?

7

You can use BoxMesh, see for instance: https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/built-in-meshes/demo_built-in-meshes.py.html?highlight=boxmesh

1 Like