Extending ex13 / divergence theorem / interfacial resistance

[gatling-nrl]

There's been several previous discussions about calculating integrated flux, and in particular that it would be more stable/accurate to use the divergence theorem to replace the integral over facets with an integral over cells. The previous conversations were in the context of flux in/out of a subdomain. This discussion is a little different (I think) because I'm working with a subset of the exterior boundary facets instead of a subset of the cells.

Is it possible to change ex13 to use this divergence theorem idea when it is calculating port flux? I can't see how to do that. Is ex13 the best way to do ex13?

I'm also interested in how I might introduce an interfacial resistivity between the "port" and the interior. This would cause the value on the facet to be different than the value just inside the domain (a discontinuity at the facet). This is perhaps part of #867, but I'm asking exactly in the context of ex13.

My application has more than 2 ports, but I assume that if I can get this working for the 2 ports in ex13, I can do it for N ports.

[gdmcbain]

G'day. I'm quite unwell at the moment (yep, covid) but I'll see if I can set you in the right direction. First, can we split this? It sounds like two separate problems. I'll treat the divergence theorem here.

Is it possible to change ex13 to use this divergence theorem idea when it is calculating port flux?

ex13 does use this divergence theorem idea:

scikit-fem/docs/examples/ex13.py

Line 51 in 4084c41

conductance = {'skfem': u @ A @ u,

for an explanation, see ‘Three ways to compute multiport inertance’; this is the third way.

Is ex13 the best way to do ex13?

Well, depends on what you mean by best. The second of the three ways is slightly more accurate but the third is so much simpler that I pretty much always use the third.

My application has more than 2 ports, but I assume that if I can get this working for the 2 ports in ex13, I can do it for N ports.

Yep, this is covered in the paper. It does involve solving N − 1 boundary value problems (which reduces to 2 − 1 = 1 in the two-port case of ex13). Basically charge one port and earth the rest in each run.

This discussion is a little different (I think) because I'm working with a subset of the exterior boundary facets instead of a subset of the cells.

I don't think this should matter, but I'm a little befogged today. Could you provide a sketch? Is it like you've got a plane region with holes in it except that the holes are filled with highly conductive stuff? What's on the outer boundary? Insulated or earthed? Say earthed. The idea then is to solve a b.v.p. for each hole with the interior set to one and the interior of all the other holes set to zero. (Here I'm approximating the highly conductive stuff as perfect; is that a reasonable simplification to get you going?) You'll find that the matrix for all N problems is the same, only the r.h.s. varies. Combine the solutions using that matrix as a bilinear form, like U.T @ A @ U, where the columns of U are the solutions of the individual BVPs. The result is the overall conductance matrix between the N ports; the diagonals are the driving-point conductance, the offs minus the transfer-conductances.

[gdmcbain]

(∇v, ρ∇u)_Ω = (v n, ρ∇u)_Γ − (v, ∇⋅{ρ∇u})_Ω

Edit: Sorry, that was written without looking at the paper. That's correct in itself but not what's required for inertance since it's the specific volume ρ⁻¹ that's analogous to conductivity.

(∇v, ρ⁻¹∇u)_Ω = (v n, ρ⁻¹∇u)_Γ − (v, ∇⋅{ρ⁻¹∇u})_Ω

[gatling-nrl]

First of all, how did you get math in markdown?!

[gdmcbain]

In (11), what's written there as U^(j) is a flux, so it's the product of a coefficient (here the specific volume ρ⁻¹) and a gradient ∇Π^(j). Just keep that all together, and (12) follows from (11); for the test flux U⁽ⁱ⁾ though on the left-hand side of (11), just use the gradient, or divide the flux by the coeffcient (multiply by density here). So:

(∇v, ρ⁻¹∇u)_Ω = (v n, ρ⁻¹∇u)_Γ − (v, ∇⋅{ρ⁻¹∇u})_Ω

Then to get to (13) the inner product over the boundary Γ needs to be split up. The idea is that the integrand will vanish except over the port of interest and therefore this can be related to the ‘reciprocal inertance coefficient’ s (more generally a conductance coefficient). This vanishing is brought about by the choice of test potential v and its boundary conditions; i.e. (6) and (7).

Thank you, I am feeling better; maybe still a bit ‘foggy’, I'm not sure.

[gdmcbain]

First of all, how did you get math in markdown?!

Markdown can include HTML, including entities like ∇ and tags like ⁻¹.

[gdmcbain]

i.e. in (11) and (12), replace U⁽ⁱ⁾ with ∇v (and so Π⁽ⁱ⁾/ρ with v) and U^(j) with ρ⁻¹∇u.

[gatling-nrl]

So I'm still quite struggling with getting accurate flux measurements out of skfem. Context is still Poisson equation with a variable coefficient. I'm using some output from COMSOL solving a toy problem as a reference. I get extremely poor accuracy using what I now call the "naïve" method of first solving for u and then integrating sigma * grad(u) over a facet subset. Here is how the hunt for a "better" method as gone so far:

I tried to extend your paper @gdmcbain to handle the variable coef, but without any success so far. I also tried to tinker with ex13.py, but got nowhere when the coefficient was non-constant.

Googling harder, I found this four year old post on scicomp which I could have written word-for-word
https://scicomp.stackexchange.com/questions/26774/computing-accurate-fluxes-with-fem
and that put me on the path of looking into "mixed Poisson problems"

That led me right into ex37.py which seemed promising, but also doesn't include variable coefficients. I did successfully extend the example though, but then discovered that Dirichlet boundary conditions are somehow hardcoded in... I can't figure by what mechanism though. So i kept digging....

From a 2012 book by Dhatt, I found

Which made me think I could ditch the RT0 element altogether and treat the (x,y) components of the flux (2d problem) separately. Basically example 37, but with three basis sets: u, flux_x, and flux_y. But this didn't work either.

I'm really beating my head against the wall at this point, because this seems like a well-solved, much discussed problem (calculating accurate fluxes) and I just can't get it to work. I also must admit to being very much a "user" of FEM and not well versed on its guts, or some of the seemingly important math concepts (like notation that reads H(div, Omega) or the concept of H-div conforming). I have no idea how I would extend ElementTriRT0 to ElementTriRT1 or why ElementTriRT0 and ElementTriP2 can't work together. At the moment, all these questions make a deep, deep rabbit hole.

What i really need to do first is specify a P0 sigma in the domain, specify a voltage on a subset of the external boundary facets, specify ground on some other disjoint subset of the boundary, and require zero normal flux everywhere else on the external boundary. Then I need an accurate calculation of the total flux between the two "electrodes" by any practical means. The potential solution inside the domain is irrelevant to me, except to sanity check things.

[gatling-nrl]

[gatling-nrl]

Edit:
The small mesh size gets quite sluggish (mostly due to gmsh i think). Here are the extremes used in the above convergence plots:

Tomorrow I will try to get some similar data out of COMSOL, because I think it is converging much faster. Obviously I would prefer to work with as coarse a mesh as possible...

[gdmcbain]

The small mesh size gets quite sluggish (mostly due to gmsh i think)

Could be, but also it will get sluggish due to the use of the direct solver.

Obviously I would prefer to work with as coarse a mesh as possible...

Yes, of course. The mesh-convergence study here is mainly to get an idea of whether everything is working correctly. We can already see that the method not marked ‘naive’ is doing much bettter on the coarser meshes than the one that is. Is that one the second or third way? Anyway, yes, when I first started using the second it was very much because itt gave reasonable engineering results (±1–10%) on reasonable engineering (i.e. coarse) meshes.

If you do find that you want to get quicker results on finer meshes, one idea is switching from a direct to an iterative solver; e.g. algebraic multigrid #264. Fortunately this is really easy in scikit-fem #273.

But once you've got decent results which seem to be converging, as they do, you can try extrapolating from the last few to get the zero-mesh size limit. A nice tool for this is NASA's VERIFY, wrapped for Python as convergence. Then instead of the flux, you can plot the discrepancy and plot that and the mesh parameter logarithmically, as in figure 3 of the ‘Three ways’. The solid orange, green, and red curves above are already looking pretty good; I can't tell from the graph whether they're ±1% or ±0.1%.

[gdmcbain]

Well, in the long run, I really want for (2)

For this, i.e. changing from Dirichlet to Robin conditions on the electrodes, what I'd think about doing is first solve the complete set of n Dirichlet problems to construct the conductance matrix for the n electrodes and then think of that as a purely resistive circuit and then add another resistor branch out from the Dirichlet electrode to the real terminal. That is, I'd work at the lumped-element circuit level rather than the finite-element level.

[gatling-nrl]

Is that one the second or third way?

naïve in that plot is

    @skfem.Functional
    def measure_flux_naive(w):
        return dot(w.sigma * grad(w.u), w.n)

    measure_flux_naive.assemble(fb2, u=fb2.interpolate(u1), sigma=fb0.interpolate(sigma))

non-naïve is

    @skfem.Functional
    def measure_flux(w):
        return w.sigma * dot(grad(w.u1), grad(w.u2))

   measure_flux.assemble(cb2, u1=u1, u2=u2, sigma=cb0.interpolate(sigma))