Implement the Poynting vector calculation

[lucasgrjn]

Regarding the issue #106:

• Implement the function calculate_poynting() for the class Mode() to calculate the Poynting vector of a mode
• Improve the docs by adding the graph Fig.1.h et Fig.1.f from the article used in the effective_area example.
  

[HelgeGehring]

Looks amazing!

Two small ideas:

• we might want to use the discontinuous variant of the element we use for Ez, that would probably allow higher accuracy in case we use a higher order and also quad meshes? (Instead of hardcoding ElementTriP1)
• what about having it as a cached property instead? Once calculated there's no need to calculate it again?

What do you think?

[lucasgrjn]

Thanks for the comments/ideas :)

• I was able to gather the element used from Ez but my knowledges on discontinuous element / mesh / quadrature remain really limited actually, so I really appreciate your insight!
• For the cached_property, it definitely make sense. I started like this but I switched to a normal function because of the poynting basis. But returns the basis + the property seems a good way! (I also add Px, Py and Pz as properties of Mode).

[HelgeGehring]

Great progress!

[HelgeGehring]

I meant as this line we should reuse the element for Ez of the basis, i.e.

Suggested change

	poynting_basis = self.basis.with_element(ElementDG(ElementTriP1()))
	poynting_basis = self.basis.with_element(ElementDG(self.basis.elems[1]))

would brobably do it.

Even better would be to have here a

Suggested change

	poynting_basis = self.basis.with_element(ElementDG(ElementTriP1()))
	poynting_basis = self.basis.with_element(ElementDG(ElementVector(self.basis.elems[1],3)))

as this way we could directly project using just one basis. This way we'd just have one field to return.

Some background:
ElementDG converts Elements in their discontinous counterpart. For example, ElementTriP1 has its degrees of freedom (DOFs) on the nodes around (Tri->Triangles->3 nodes/dofs).
As those nodes are shared with the neighboring triangles, DOFs are shared with the neighbors (just imagine the triangles's node-heights being the values and the triangles the surfaces in-between the nodes. A continuous element would lead to a continuous surface as the values at the nodes are reused.)

Here, we don't expect it to be continuous as epsilon makes discrete jumps as the boundaries, thus we have to consider that in how we choose the element.
For maxwell's equations we know that the tangential component of the field is continuous while the normal component is not -> best choice for Ex/Ey is a H(curl)-conforming element like the Nedelec-element. For Ez we know it's continuous so we're best of with a Lagrange element.

For the pointing vector, i'd guess it's save to use a vector of Lagrangian elements. As we see in the plots that it's discontinuous at boundaries, we need the discontinuous version of it.
Here, the question would be if there's any continuity across boundaries of the poynting vector which we could have in our choice of element as that would reduce the needed number of DOFs and would probably also be a bit more precise. (@DorisReiter do you know if there's any component of the pointing vector continuous across boundaries?)

For now I think we can just get it in with the Vector of Lagrangian elements, as we don't do a solve using those elements.

[lucasgrjn]

I see for the use of basis.elem.elems.
(In the beginning I wanted/was thinking to use split() and get the elements from En_basis, but this way is much cleaner!)
(Corrected on 98cb853)

[lucasgrjn]

Thanks for the explanation! :)

I don't get the second part of your message to use a function like self.basis.with_element(ElementDG(ElementVector(self.basis.elems[1],3))). Why will we be allowed to project using one basis? Don't we already make that with the poynting_basis.project() ?

For the Poynting vector, from what I remember, you have the Poynting theorem which is an expression for the conservation of energy:
$$-\frac{\partial u}{\partial t} = \nabla\cdot\mathbf{S} + \mathbf{J}\cdot\mathbf{E}$$
with $u$ the energy density, $S$ the Poynting vector, $J$ the current density and $E$ the electric field. If we consider a mode (i.e. no source), we should obtain $\nabla\cdot\mathbf{S}=0$ but I think @DorisReiter will have better insights!

[HelgeGehring]

at the moment, we have three sets of DOFs as the element is a scalar element.
if we use a vector-element we can have all three components in one set of DOFs (which we can split if we need)

something like

poynting_basis = self.basis.with_element(ElementDG(ElementVector(self.basis.elems[1],3)))
... = poynting_basis.project(np.stack([Px,Py,Pz]))

should do it (I always forget if stack/hstack/vstack, it should be stacked as a new last axis)

alternatively, there should also be a way to combine the DOFs into one set and use it with the vector basis, but I think it's better readable if we treat it directly as a vector.

About the continuity: that equation doesn't tell us about the interface-conditions of a single component (as a discrete change in one component can be compensated in another one)
So I guess it's better to assume it's discontinuous (the pictures you posted/in the paper also show that it's not continuous across interfaces)

[HelgeGehring]

I'm also wondering if treating it as a vector makes sense at all 🤔
we're calculating eigenmodes, which should mean that the mode propagates only in z direction...
That lead to $S_x$=$S_y$=0?
Or is it only the integral over the whole domain which needs to be zero for those two components?

Edit: seems like it can be non-zero https://doi.org/10.1364/OE.472148

[lucasgrjn]

I'm also wondering if treating it as a vector makes sense at all 🤔 we're calculating eigenmodes, which should mean that the mode propagates only in z direction... That lead to Sx=$S_y$=0? Or is it only the integral over the whole domain which needs to be zero for those two components?

Honestly that's a good question... I also found an article which states about that. So I think it is safer to keep the 3 components :)