Fix visualization of flux densities

[jfeist]

Hi Homer, here's my fix for the plots of flux densities. For one example calculation, this changes the plot of power flux density from looking like this: to looking like this:

These plots are done with GMSH by turning off lighting, so you don't explicitly see the triangles of the grid, but I found that this best shows how smooth the flux plots come out in the new version, and how there are clearly errors in the previous one.

I included my arguments of why I think it should be the formulas I am using in the comments, but I'm not at all sure about whether even my assumptions true. What I do know is that like this, it works.

BTW, I'm not convinced that the absolute scale of the results is correct in scuff units (i.e., that this is power/area/frequency in scuff units for the spectral power flux density we are plotting). But there should only be a global scaling factor missing. Víctor is checking this for a simple case where we can compare with analytical results for the power flux density. I'll update this pull request (or send a new one if this one is merged in the meantime) once we figure out the correct absolute scaling.

EDIT: Forgot to say: this includes two commits, where the first one uses a slightly different approach that gives almost identical results. I think the second one is more correct, but I don't have proof. So I'm including the first one as well, just in case you want to compare.

[HomerReid]

Thanks for looking into this and for contributing a fix that works. A considerable piece of detective work! And thanks for taking the time to explain your reasoning in the comments.

I spent some time going over this today and I still can't see what was wrong about the previous version, nor why your proposed fix is correct (though I acknowledge it is). I will need to write up some notes on this in an attempt to pin it down precisely. What is being plotted in this case is actually not a tangential-vector-valued quantity approximated by an expansion in RWG basis functions, but a scalar quantity that can be decomposed into a sum of contributions associated with each basis function. The total quantity is the sum of the contributions of each edge, and each edge is connected to two vertices, so it seems natural to assign half the contribution of each edge to each of the two vertices to which it is connected. Then the density at a given vertex should be the total value assigned to that vertex, divided by the total area assigned to that vertex. The total area assigned to the vertex should just be the area of the region lying closer to that vertex than to any other vertex (the Voronoi cell?) That was the reasoning behind the original calculation.

Your calculation seems instead to be (a) divide the contribution of each edge by the total area of the edge to obtain a "density-per-edge," (b) for each vertex, average the density-per-edge of all edges connected to the vertex and take that to be the density associated to the vertex. I have to think a little about why this is the right thing to do. I'll try to get this straight in conjunction with the revamping of the units of output quantities reported by scuff-neq.

[jfeist]

Thanks for the additional details! I think the "trick" here is that for plotting the flux density (or any function), we need to find the best approximation for the "raw" function value at the vertex positions. Your argument sounds to me like it would possibly conserve the sum (i.e., whereas before the total flux was the sum of all edge contributions, you are reexpressing it as a sum over vertex contributions). That means that the values assigned to each vertex implicitly contain the function value times some weight related to the size of the triangles.

However, for plotting, we need the actual function value at the vertex, not the function times some weight (and we are not interested in conserving the total sum, in the sense that the sum over vertex values by itself is not an approximation to the integral of the function over the area).

My argument is then that each edge contribution corresponds to the integral of the "raw" function over the area associated to the edge. The average value of the function in that area is therefore the edge contribution divided by the area. For each vertex, there are a number of edges touching it, each of which have different values of the function. A reasonable approximation for the actual function value at the vertex is thus the average of all these values.

Now, some caveats:

1. I think that in a really correct argument, the shape of the RWG function must also come in. I guess the flux is calculated as the integral of the dot product (or normal part of the cross product?) of some function with the vectorial RWG basis function, and the value associated with that edge thus corresponds to the prefactor of a "new" scalar basis function that is maximum along the edge and 0 at the opposing vertices. Therefore, the correct normalization to go from edge contribution to raw function value might not use just the triangle area, but the integral over that new basis function - fortunately, that is just the area times some prefactor (2/3 if its linear as I think). Note that part of this argument is already taken into account, in that each edge only contributes to the two vertices it touches, and not to the "opposing" vertices at the tips of the two associated triangles.
2. Additionally, I think my approach could be improved by doing the final averaging with nonequal weights (e.g., according to the opening angle of the triangle associated with each edge around that vertex, which would correspond to what you might get from a small "ring" integral around the vertex).

For sure, I agree with you that this should be worked out correctly, which is hard to do here since you can not even write equations... Also, since I do not even know the starting equation that gives the edge contributions, my arguments include quite a bit of guesswork. I'm happy I got lucky this time!

[HomerReid]

The fix you proposed (the "better" one) has been implemented in 5fd499bdae3bd232e16ef54f27bf198ff87a8b11. Thanks again for this contribution.

In thinking through this I realized that the current scheme for spatial visualization might not be optimal, because it computes density values at mesh vertices (which GMSH then interpolates over triangles to yield color plots), but the densities in question are discontinuous at mesh vertices. An alternative might be instead to compute density values at panel centroids and ask GMSH to interpolate these over the elements of the "dual" mesh, i.e. the Voronoi cell corresponding to each vertex. The resulting plot might be a better way to visualize the density in question (I'm not sure). I will investigate.

[jfeist]

Great, thanks! I'll close this pull request.