Hi Geoff,
On 08/06/2014 09:58 PM, Geoff Wright wrote:
Hi Robert,
I'm attaching the files. Also attaching a screenshot from paraview of the voltage solution.
Two questions:
- I know you were suggesting I use a higher order for the voltage field. But I haven't found that it makes much difference. For example, with approx_order==1 results are: Source flux: 5.70041929669 Sink flux: -3.1885130653 Block flux: -0.222003267731 Entire surf: 2.28990296366
My results with approx_order==1:
Source flux: 0.113939919911 Sink flux: -0.094333017518 Block flux: 1.73879971396e-05 Entire surf: 0.0196242903901
Much better, right? :) That's because I used the same resistivity in the laplace equation and in the flux evaluation. You were using the getResistivity() output in the equation and
m = Material('coef', sigma=nm.eye(3))
in flux evaluation. So use
out['val'] = mult
out['sigma'] = mult * np.eye(3)
in getResistivity() and
snkf = pb.evaluate('d_surface_flux.i2.Gamma_Left(p.sigma, v)', mode='eval')
instead of
snkf = pb.evaluate('d_surface_flux.i1.Gamma_Left(coef.sigma, v)', mode='eval', coef=m)
etc in post_process(),
and with approx_order==2, results are: Source flux: 6.47939779365 Sink flux: -3.90908956187 Block flux: 0.0399662212412 Entire surf: 2.61027445302
I will try that tomorrow on a bigger machine...
Its true that the amount of leakage in the block region has decreased, but the proportion of current that is observed at the Sink relative to the source is 56% in the first case and 60% in the second case... not much difference.
- I've been thinking more about your explanation for the flux error and I don't get it. If the solver is properly enforcing the Neuman boundary condition of zero flux on the surface then it should not reach a solution with residual approx equal to zero. I.e. the non zero flux on the surface should contribute to the residual and the system would not converge. However when I look at the residual I see something like e-10. Can you explain this? I'm wondering whether we just need another term in the equations that explicitly enforces the Neuman boundary condition in the "block" region. What do you think?
The Neumann boundary condition can be enforced by adding a term to the equation. When it is zero, the term is zero, so it disappears from the equation. Maybe you can try adding a dw_surface_flux(...) = 0 equation via a Lagrange multiplier, but I do not see how it could be better.
Everything is in the weak sense with respect to the test function space. And the gradient computation on surface is a sort of extrapolation from the inside of the boundary element - it is influenced by the inside of the domain (the normal component of the gradient) -> the smaller the boundary elements the better. That is all explanation I can come up with now.
r.