That's ok - this can happen, when the initial residual is high. Try increasing 'lin_red' option of the Newton solver to 1e1.

More details: Usually, one wants the linear solver precision match the nonlinear solver precision, so sfepy warns you when the linear system solution rezidual norm is > than eps_a * lin_red of the nonlinear solver, and default values are eps_a = 1e-10, lin_red = 1. In this case the relative rezidual is 2.823174e-15 which is about as good as it can get in double precision.

Conclusion: ignore the message :)

r.

On 06/06/2012 05:35 PM, steve wrote:

Hmm.. also I just noticed that simple.py complains that all is not well. The output really

*looks*reasonable though!sfepy: reading mesh (/Users/steve/Development/sfepy/es-lens.mesh)... sfepy: ...done in 0.10 s sfepy: warning: bad element orientation, trying to correct... sfepy: ...corrected sfepy: creating regions... sfepy: Omega_Source sfepy: Omega_Lens2 sfepy: Omega_Lens1 sfepy: Omega_Target sfepy: Omega sfepy: Omega_Needle sfepy: ...done in 0.04 s sfepy: equation "Laplace equation": sfepy: dw_laplace.i1.Omega( m.val, v, u ) = 0 sfepy: setting up dof connectivities... sfepy: ...done in 0.00 s sfepy: using solvers: nls: newton ls: ls sfepy: updating variables... sfepy: ...done sfepy: matrix shape: (4465, 4465) sfepy: assembling matrix graph... sfepy: ...done in 0.00 s sfepy: matrix structural nonzeros: 30483 (1.53e-03% fill) sfepy: updating materials... sfepy: m sfepy: ...done in 0.00 s sfepy: nls: iter: 0, residual: 4.936809e+04 (rel: 1.000000e+00) sfepy: rezidual: 0.00 [s] sfepy: solve: 0.03 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.335327e-10< 1.000000e-10) sfepy: nls: iter: 1, residual: 1.393747e-10 (rel: 2.823174e-15)

-steve

On Wednesday, June 6, 2012 9:04:39 AM UTC-6, Robert Cimrman wrote:

On 06/06/2012 04:33 PM, steve wrote:

Hi Robert,

OK. So it seems that for the above, you could use directly dw_laplace

and

dw_volume_dot terms, with material arguments equal to 2\pi r and 2\pi

k^2

r, respectively. Or just r, and put the known constants directly into

the

equations string.

Hmm.. OK. I'll try to wrap my head around that. ;-)

I will try to support you, don't worry.

Well... I've had some success here, I think!

Glad to hear that!

here's the setup file 'es-test.py' which is trying to model a

cylindrically

symmetric electrostatic lens. The results compare favorably with a relaxation code with the same setup. Now of course I have a couple of issues that I hope you can help me with.

Nice! Could you send me also the mesh? (off-list, if it's big...) So this just solves the direct problem, not the eigenvalue problem, right?

- In the relaxation version I have to impose the du/dr=0 boundary
condition at r=0 at each time step. I made no mention of that BC in the setup file, but the solution seems to respect it more or less automatically. I remember in the examples seeing mention of "zero flux" over the boundary. Is this just a consequence of leaving out the

boundary

integral in the weak formulation?

Yes, by leaving out the boundary integral you impose the zero flux.

- I've developed a simple wxpython app that allows a user to specify
voltages on, and position/geometry of the lenses and then compute

particle

trajectories through the system. My relaxation algorithm is currently limited to a uniform rectangular mesh. Part of my motivation for

pursuing

sfepy is the possibility of using variable/arbitrary mesh density.

However,

to run this test I had to generate a mesh in gmsh and then hand-edit the output to get a .mesh file that was acceptable to sfepy. I'd like to increase the density of nodes around the corners of the conductors in

the

problem (especially around the needle where the particles are

introduced).

Any thoughts on dynamically building meshes in code to accomplish this?

What had to be changed in the gmsh-generated mesh? INRIA medit format (.mesh) should work - if it does not, it should be fixed.

As sfepy cannot do the meshing/adaptive refinement, so you will need to call an external mesher (e.g. gmsh) to do that. Skimming the gmsh docs seems to suggest that it's possible to drive it from Python. If you choose to go that path, I would really like to see the result - having a MeshIO class taking a gmsh mesh description file and returning the mesh would be great.

Anyway.. thanks for the help! I think it's working pretty well after I finally groked how variable materials work in the setup file.

Hth, r.