Poisson equation, implementation of boundary values, mixed finite element, semiconductor devices

[Konrad Wiśniewski]

 Hi dear deal.ii community!

I am trying to solve 2D transient problem in semiconductor devices and I'm stuck with application of Dirichlet boundary values via ConstrainMatrix or AffineMatrix (in the last version of this library)

The program (which I am rewriting) solves consecutively two equation: (i) for a given charge density it solves Poisson equation, and than (ii) the continuity equation to find the densities in next step.

The poisson equation looks like this:

This is a formulation of mixed finite element problem with phi as electric potential and E = (Ex, Ey) as a vector of electric field. 

Program uses a Raviart-Thomas finite element to describe E, and Discontinous Galerkin (FE_DGQ) for potential. 

Weak fe formulation:

The geometry (as well as boundary conditions) is quite simple:

I have choosen this kind of boundary conditions because:

- the potential is applied on the left side of device

- Interesting things (non-zero carrier densities) shound appeard only near the interface of two domains  so I assumed that Ex=0, Ey=0 on Dirichlet boundaries

- there should be no current through top and bottom of domains, so  I decided in most simple case to have zero electric field perpendicular to the boundary.

I hope those are reasonable bc. 

My implementation of boundary conditions goes like this:

For a given FESystem fe and DoFHanfler dof_handler:

dof_handler.distribute_dofs(fe);

DoFRenumbering::component_wise(dof_handler);

// make hanging node constraints

constraints.clear();

DoFTools::make_hanging_node_constraints(dof_handler, constraints);

// add constaints to Poisson equation

const FEValuesExtractors::Vector ElectricField(0);

ComponentMask electric_field_mask = fe.component_mask(ElectricField);

DoFTools::make_zero_boundary_constraints(dof_handler,

Dirichlet, // DIRICHLET BOUNDARY INDICATOR

constraints,

electric_field_mask);

std::set<unsigned int> no_flux_boundary;

no_flux_boundary.insert(Neumann);

VectorTools::compute_no_normal_flux_constraints(dof_handler,

                                0,

                                no_flux_boundary,

                                constraints);

                constraints.close();

during the procedure of assembling the global matrix I use in each loop:

constraints.distribute_local_to_global( data.local_matrix,

data.local_dof_indices,

Poisson_object.system_matrix);

and for right hand side:

constraints.distribute_local_to_global(data.local_rhs,

       data.local_dof_indices,

       Poisson_object.system_rhs);

Unfortunately when I look at the output the E value on Dirichlet boundary only x-component Ex is 0, but not Ey. On Neumann boundary non of the components are 0, but of course Ey should be. Additionally the value of a potential (which comes from face integral on right hand side) on Dirichlet boundaries is not longer as it should be (eg. I ascribed 1V and in the output I read 1e-5V)

What did I do wrong? Is it because the Raviart-Thomas should be treated differently? Maybe the value that I look at in Paraview is extrapolated from the values from the middle of the cells face and that is why the values are incorrect?

Another problem is that for the simplest case all Ey on all of the domain should be zero. Maybe this is the problem, but I kind of don't want to write a code especially for the simplest case when Ey=0 -> is there a better way to do this?

I will appreciate any help! 

[Wolfgang Bangerth]

Konrad,

> I am trying to solve 2D transient problem in semiconductor devices and
> I'm stuck with application of Dirichlet boundary values via
> ConstrainMatrix or AffineMatrix (in the last version of this library)
>
>
> The program (which I am rewriting) solves consecutively two equation:
> (i) for a given charge density it solves Poisson equation, and than (ii)
> the continuity equation to find the densities in next step.
>
> The poisson equation looks like this:
>
>

> eq.png
>
> This is a formulation of mixed finite element problem with phi as
> electric potential and E = (Ex, Ey) as a vector of electric field.

I suspect you have already found that out (even though you don't state
it in your email) that this is *exactly* the problem step-20 solves.

> The geometry (as well as boundary conditions) is quite simple:

This doesn't work. You can't impose boundary conditions on *both* Phi
and E at the left and right sides. It would be equivalent to imposing
both Dirichlet and Neumann values for the non-mixed Laplace equation.

> I have choosen this kind of boundary conditions because:
>
> - the potential is applied on the left side of device
>
> - Interesting things (non-zero carrier densities) shound appeard only
> near the interface of two domains  so I assumed that Ex=0, Ey=0 on
> Dirichlet boundaries

This assumption is unphysical. From your first equation, it is clear
that (up to some constant), E=-grad Phi. So if E=0 at the boundary,
you'd have grad Phi=0 AND ALSO Phi=something on these boundaries.

> Unfortunately when I look at the output the E value on Dirichlet
> boundary only x-component Ex is 0, but not Ey.

That's because for the R-T element, the degrees of freedom at a face are
only the normal component of the vector. That's the only thing you can
constrain, i.e., you can only prescribe E.n=something, but not the
tangential component.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Konrad Wiśniewski]

Than you for your response!

and if I may ask for some advise:

I tried to give some additional constrains for electirc field E (I must admit that it wasn't good idea) to avoid another unphysical results. General physical poin of view is that I should have this kind of a situation:

The Ex in the middle should go from 0 to -1e7 or -1e6 in the middle of device and then go back to zero 

but my results (for Ex) look rather like this:

Where all Ex values, over all domain is almost -9e7 and it drops a little in the middle (as it should, but not from such a big value). 

I couldn't find any mistake in my implementation so I've dicided to change boundary condition to enforce electric filed E to be 0. 

Do you have an intuition what I'm making wrong? 

Probably I can subtract the values on the boundary from all of the results but I don't think it is proper way to deal with this problem.

[Wolfgang Bangerth]

Konrad,

> I tried to give some additional constrains for electirc field E (I must admit
> that it wasn't good idea) to avoid another unphysical results. General
> physical poin of view is that I should have this kind of a situation:
>

> str.png
>
>
> The Ex in the middle should go from 0 to -1e7 or -1e6 in the middle of device
> and then go back to zero
> but my results (for Ex) look rather like this:
>

> results.png
>
>
> Where all Ex values, over all domain is almost -9e7 and it drops a little in
> the middle (as it should, but not from such a big value).

What is it that you show in your figure? The x-component of the electrical
field? Can you also show the potential? Recall that the x-component of the
electrical field is something like the x-derivative of the potential Phi, so
it seems to me like you should have a linearly decreasing potential with a
large gradient. Does the potential satisfy the boundary conditions you impose?

> I couldn't find any mistake in my implementation so I've dicided to change
> boundary condition to enforce electric filed E to be 0.

This is some more general advice: If you get to a situation that you don't
understand (electric field looks wrong) and you "fix" it by adding conditions
you have no idea whether they are correct/physical/mathematically allowed
(adding boundary conditions for E), then in all likelihood you've just added a
second bug to the first one. You really should get in the habit of
investigating and *understanding* what is going wrong when when something
looks odd, rather than *papering over* the issue. In the long run, this will
save you time; it will also make you a better programmer and computational
scientists if you build a mental toolbox for how to debug problems.

[Konrad Wiśniewski]

Wolfgang,

thank you very much once again! 

Yes, on the picture is Ex, and potential was exactly as you described. Today I was able to correct my calculations. I changed the right hand side boundary condition to be only E=0 (what is, I think, physically correct thing to do), and after I get the proper shape of a solution I was able to deduce where my mistake was. 

It looks like even the writing about a problem on forum is beneficial - it somehow structuring the problem. 

Thank you very much for your time! Your advises and answers really helped me!

Now I will start to test continuity equation so it is not excluded that I will ask for help in the near future on this forum (although I hope it won't be necessary :). 

best regards,

Konrad

[Konrad]

Hi Konrad,

I may add as a comment (without going into details): Your problem is - as you correctly describe above - a mixed Poisson problem. You essentially have two options to solve it and the way you choose influences the way you need to treat the boundary conditions: in primal form and in mixed form. In your formulation you have natural boundary conditions (Dirichlet) on Phi which means that you need to prescribe them in the variational form. This will lead to a boundary integral on the relevant part of the boundary and result in an additional term on the right-hand side. The boundary condition for E is a flux condition and in your mixed form this (Neumann) condition is essential. This means that you need to enforce them on the Raviart-Thomas space as a hard constraints. There are essentially two options for this:

1. If you have given a vector field that describes E on the relevant part of the boundary you can use the function

VectorTools::project_boundary_values_div_conforming. This will project the given vector field onto its normal component (since this is the one you must describe with RT elements in H(div)) and use the projected value as boundary condition for the normal flux of E

2. If you only are given the scalar value of the normal flux of E then you may have to set the relevant dofs in the RT-space by hand. This you can do through constraints that way:

    AffineConstraints<double>           constraints;

    constraints.clear ();

    DoFTools::make_hanging_node_constraints (dof_handler, constraints);

    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
    const unsigned int   dofs_per_face   = fe.dofs_per_face;
    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
    std::vector<types::global_dof_index> local_dof_face_indices (dofs_per_face);

   

    typename DoFHandler<dim>::active_cell_iterator
                                    cell = dof_handler.begin_active(),
                                    endc = dof_handler.end();
    for (; cell!=endc; ++cell)
    {
        cell->get_dof_indices(local_dof_indices);

        for (unsigned int face_n=0;
                face_n<GeometryInfo<dim>::faces_per_cell;
                ++face_n)
        {
            cell->face(face_n)->get_dof_indices(local_dof_face_indices);
            if (cell->at_boundary(face_n))
            {
                if (cell->face(face_n)->boundary_id()==id_bottom_boundary)
                {
                    for (unsigned int i = 0; i<dofs_per_face; ++i)
                    {

                        // see documentation for these functions
                        constraints.add_line (local_dof_face_indices.at(i));
                        constraints.set_inhomogeneity (local_dof_face_indices.at(i),
                                value_of_dof);
                    }
                }
                else // add some constraints on other boundaries. Here: zero-flux
                {
                    for (unsigned int i = 0; i<dofs_per_face; ++i)
                    {
                        constraints.add_line (local_dof_face_indices.at(i));
                    }
                }
            }
        }
        ++cell_n;
    }

    constraints.close();

Once you use the DoFTools::make_sparsity_pattern function you need to make sure to keep the inhomogenities:

        DoFTools::make_sparsity_pattern(dof_handler,
                                              dsp,
                                              constraints,
                                              /*keep_constrained_dofs = */ true);

The same in the assembly routine:

            cell->get_dof_indices(local_dof_indices);
            constraints.distribute_local_to_global(local_matrix,
                                                                local_rhs_local,
                                                                local_dof_indices,
                                                                system_matrix,
                                                                system_rhs,
                                                                /* use inhomogeneities for rhs */ true);

Keep in mind that the dofs for RT elements are edge moments in 2D and face moments in 3D (and not nodal values).

Best,
Konrad

[Konrad Wiśniewski]

Hi Konrad,

Thank you very much for your comment! I will try to go into the detail of it in next days. I must admit that to this point I've just used ::make_zero_boundary_constraints with regard to the conditions on E field. I didn't have a better idea, and your solutions seem much more better. I will probably use the second option, but I need to read more about the first function that you mentioned. Anyway - it will help me for sure, thanks! We will see if I will be able to proceed without any further advice :). 

Best

Konrad