Using the solution from one problem as a boundary condition in another (decoupled)

[krei]

Hi,

I have solved a Laplace problem in a domain (according to the early tutorial steps 3 and 4) and have received the solution for electric potential. Now, in a neighboring domain I want to solve a similar PDE (dependent variable is also electric potential, but the equation is slightly different, it calculates the current distribution due to emission currents, but this is not important) while using the electric field from the first problem as a boundary condition on one boundary.

In the second problem, to impose the non-homogeneous Neumann BC-s, I have to make a class extending Function and implementing a method value(), right? And somehow in this method value() I have to obtain the electric field from the first solution.

How would I evaluate the gradient (electric field) of the first problem's solution (electric potential) in a point on the boundary?

Best regards,
krei

[Simon Sticko]

Hi,
it sounds like you could use

VectorTools::point_gradient

/Simon

[Wolfgang Bangerth]

On 01/07/2016 07:19 PM, krei wrote:
>
> I have solved a Laplace problem in a domain (according to the early
> tutorial steps 3 and 4) and have received the solution for electric
> potential. Now, in a neighboring domain I want to solve a similar PDE
> (dependent variable is also electric potential, but the equation is
> slightly different, it calculates the current distribution due to
> emission currents, but this is not important) while using the electric
> field from the first problem as a boundary condition on one boundary.

Can you put this into formulas?

> In the second problem, to impose the non-homogeneous Neumann BC-s, I
> have to make a class extending Function and implementing a method
> value(), right? And somehow in this method value() I have to obtain the
> electric field from the first solution.
>
> How would I evaluate the gradient (electric field) of the first
> problem's solution (electric potential) in a point on the boundary?

As already suggested, you could use VectorTools::point_gradient. But
this may be more work than necessary if the meshes for the two parts of
your domain meet in some sort of regular way. Do their vertices and node
locations match?

Best
W.

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

[krei]

Hi,

The domains of the two problems share a boundary, but I made two different meshes for 1) Vacuum over a metal surface (http://i.imgur.com/R2KOiXB.png); and 2) the metal (http://i.imgur.com/XIBAsg2.png).

I vacuum, I want to solve the Laplace equation for electric potential such that the boundary condition at the surface of the metal is just potential = 0. This I have already managed to do.

In the metal part I want to solve a stationary equation for currents (hopefully this latex is displayed correctly): [; \nabla \cdot \sigma(x) \nabla \varphi = 0 ;]

with Neumann boundary condition on the copper-vacuum boundary (represent emission currents): [; \vec{n} \cdot \vec{J}(E_\text{vac}) = \vec{n} \cdot \sigma(x) \nabla \varphi  ;]

From my initial testing with VectorTools::point_gradient as Simon suggested, It seems to work to implement this BC, but I would love to hear suggestions how to solve the problem in a better way.

Eventually, I also need to couple the currents equation with the heat equation [; \nabla \cdot \kappa(x,T) \nabla T = - \frac{J^2}{\sigma(x,T)} ;] and also the emission currents and the electric conductivity become temperature dependent. I am currently studying step-21 on how to couple multiple PDE-s, but any suggestions on this are also welcome.

Best regards,
krei

[Wolfgang Bangerth]

On 01/08/2016 06:31 PM, krei wrote:
>
> The domains of the two problems share a boundary, but I made two
> different meshes for 1) Vacuum over a metal surface
> (http://i.imgur.com/R2KOiXB.png); and 2) the metal
> (http://i.imgur.com/XIBAsg2.png).
>
> I vacuum, I want to solve the Laplace equation for electric potential
> such that the boundary condition at the surface of the metal is just
> potential = 0. This I have already managed to do.
>
> In the metal part I want to solve a stationary equation for currents
> (hopefully this latex is displayed correctly): [; \nabla \cdot \sigma(x)
> \nabla \varphi = 0 ;]
>
> with Neumann boundary condition on the copper-vacuum boundary (represent
> emission currents): [; \vec{n} \cdot \vec{J}(E_\text{vac}) = \vec{n}
> \cdot \sigma(x) \nabla \varphi ;]
>
> From my initial testing with VectorTools::point_gradient as Simon
> suggested, It seems to work to implement this BC, but I would love to
> hear suggestions how to solve the problem in a better way.

You could try to go into the implementation of
VectorTools::interpolate_b_v to see what this function does. You will
find that it requires the evaluation of the boundary values at the
"support points", i.e., the points where the finite element shape
functions are defined. For example, for a Q1 element, these are just the
vertices of faces at the boundary, whereas for Q2 elements there are
also the face and edge midpoints. If your two meshes match, then the
points where you need to evaluate the gradient of one solution happen to
be also the vertices of the other solution, and this can make many
things simpler. Specifically, you would only need to evaluate the first
solution (as boundary conditions for the second solution) at points that
are quite well defined and easy to find, rather than having to do what
VectorTools::point_gradient does (namely, loop over all cells, find the
appropriate cell, evaluate the solution there). The result should of
course be the same, but it would be more efficient. If, on the other
hand, your assembly is currently not the bottleneck, then there is no
reason for you to go to a more complicated scheme than you currently
already have.

> Eventually, I also need to couple the currents equation with the heat
> equation [; \nabla \cdot \kappa(x,T) \nabla T = -
> \frac{J^2}{\sigma(x,T)} ;] and also the emission currents and the
> electric conductivity become temperature dependent. I am currently
> studying step-21 on how to couple multiple PDE-s, but any suggestions on
> this are also welcome.

This is a classical problem in solving coupled problems. step-21 is an
example of solving something like this. step-15 also provides a number
of tools that you will need (specifically on how to evaluate solutions
at quadrature points).

[krei]

Thanks for the reply. Unfortunately, the assemble system is a bottleneck currently (takes about 4 times more time than the solver) with using VectorTools::point_gradient. The mesh on the corresponding boundaries of the two domains (http://i.imgur.com/R2KOiXB.png and http://i.imgur.com/XIBAsg2.png) doesn't match up and I think I would eventually want to solve the system with Newton's method by starting initial iterations on a coarse grid and refining it after each step. In this case, a solution which works even without matching boundary meshes would be very convenient. I also checked the code of interpolate_boundary_values, but I don't really see how it could be adapted to the present case.

Best regards,
krei

[Wolfgang Bangerth]

On 01/14/2016 04:46 PM, krei wrote:
>
> Thanks for the reply. Unfortunately, the assemble system is a bottleneck
> currently (takes about 4 times more time than the solver) with using
> VectorTools::point_gradient. The mesh on the corresponding boundaries of the
> two domains (http://i.imgur.com/R2KOiXB.png

> <http://www.google.com/url?q=http%3A%2F%2Fi.imgur.com%2FR2KOiXB.png&sa=D&sntz=1&usg=AFQjCNE-xZo2OfDOfM1NEiYswbiBGuQwxA>
> and http://i.imgur.com/XIBAsg2.png
> <http://www.google.com/url?q=http%3A%2F%2Fi.imgur.com%2FXIBAsg2.png&sa=D&sntz=1&usg=AFQjCNF6IfyUsoDoWT9Y5m5WCzW7Ccnz4w>)

> doesn't match up and I think I would eventually want to solve the system with
> Newton's method by starting initial iterations on a coarse grid and refining
> it after each step. In this case, a solution which works even without matching
> boundary meshes would be very convenient. I also checked the code of
> interpolate_boundary_values, but I don't really see how it could be adapted to
> the present case.

It's more complicated indeed. The best you can do in that case is to compute
once and for all which faces from one subdomain are adjacent to the faces of
the other subdomain. This way at least you won't have to search for the
appropriate cell every time you want to evaluate a point in point_gradient().

It is true that the algorithm would have to be more sophisticated -- working
with multiple, non-matching meshes is non-trivial.

[krei]

Hi, thanks for the idea. To implement it, I have to define a custom point_gradient method and change the line

const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
cell_point
= GridTools::find_active_cell_around_point (mapping, dof, point);

to something that would only loop through the elements on the defined boundary. One way would be to somehow modify the DofHandler dof object to only contain the boundary cells, but I am not sure if this is possible. Another way is to go deeper, and define a custom  find_active_cell_around_point and find_closest_vertex methods, where in latter instead of

const std::vector< Point<spacedim> > &vertices = tria.get_vertices();

a precalculated subset of vertices is given. But in this case, there is a lot of code duplication and I am not sure about the details.

Could you perhaps give me some advice on how you would implement this?

[Wolfgang Bangerth]

On 01/24/2016 07:23 AM, krei wrote:
>
> |
> conststd::pair<typenameDoFHandler<dim,spacedim>::active_cell_iterator,Point<spacedim>>
> cell_point
> =GridTools::find_active_cell_around_point (mapping,dof,point);

> |
>
> to something that would only loop through the elements on the defined
> boundary. One way would be to somehow modify the DofHandler dof object to only
> contain the boundary cells, but I am not sure if this is possible. Another way

> is to go deeper, and define a cu|stom |||find_active_cell_around_point|and

> ||find_closest_vertex |||methods, where in latter instead of
>
> |
> conststd::vector<Point<spacedim>>&vertices =tria.get_vertices();
> |
>
> a precalculated subset of vertices is given. But in this case, there is a lot
> of code duplication and I am not sure about the details.
>
> Could you perhaps give me some advice on how you would implement this?

Since you're going to call the function that finds the appropriate boundary
cell/face many times, you want to pull as much of the work out to a place
where you can do it once rather than every time the function is called.

For example, you could just create a std::vector<active_cell_iterator> into
which you put all of the cells on the boundary of interest. In the function
that computes the point values or gradients, you would then only have to loop
over those when you try to find the appropriate cell.