Hi Robert,

You're right! I think I must have been trying to process the wrong output file before, so I thought it wasn't producing a solution. The output file does indeed look like the correct solution.

The message is a bit scarey, so I'm glad you changed it.

Many thanks,

• Phil

On Tuesday, March 25, 2014 6:04:12 PM UTC, Robert Cimrman wrote:

Hi Phil,

On 03/25/2014 05:51 PM, seismo_phil wrote:

...

Hi,

Thanks for inviting me to contribute usage issues, Robert. I'm using sfepy to solve a relatively simple, 2D, thermal diffusion problem. i can solve it useing the following ebcs and equations (i.e., implied zero heat flux at left and right boundaries):

ebcs = {

't1' : ('Gamma_Top', {'t.0' : 0.0}),

't2' : ('Gamma_Bottom', {'t.0' : 1000.0}),

}

equations = {

'Temperature' : """dw_laplace.2.Omega( coef.val, s, t )

= 0 """

}

This works perfectly. I actually need to solve the same problem with an ebc at the top but a Neumann boundary condition at the bottom. So I use:

ebcs = {

't1' : ('Gamma_Top', {'t.0' : 0.0}),

}

equations = {

'Temperature' : """dw_laplace.2.Omega( coef.val, s, t )

= dw_surface_integrate.2.Gamma_Bottom(flux.val, s)"""

}

where flux is defined in materials as:

materials = {

 'flux' : ({'val' : 1.0},), 

 'coef' : ({'val' : 1.0},), 

}

This seems to be exactly like what is done in the example poisson_neumann.py But now, I find that the solver fails:

sfepy: updating materials...

sfepy: coef

sfepy: flux

sfepy: ...done in 0.01 s

sfepy: nls: iter: 0, residual: 9.332878e+04 (rel: 1.000000e+00)

sfepy: rezidual: 0.00 [s]

sfepy: solve: 0.00 [s]

sfepy: matrix: 0.00 [s]

sfepy: linear system not solved! (err = 5.659775e-10 < 1.000000e-10)

sfepy: nls: iter: 1, residual: 5.748516e-10 (rel: 6.159425e-15)

I'm confused by this, because I thought that this is still a well-posed problem. Can anyone please suggest to me what I may be doing wrong?

The solution is perfectly fine. Look at the sizes of residual norms - the absolute one is about 5e-10, the relative one about 6e-15, which is as close to machine precision in double-precision floats as you can get. The message is just a warning, and it can be prevented by setting the options 'lin_red' or 'eps_a' of the Newton solver to higher values. There is a check that the linear system solution error is smaller than (eps_a * lin_red). Anyway, it almost passes with your settings.

Cheers, r.

On 03/26/2014 04:23 PM, seismo_phil wrote:

BTW, my next job is to calculate the heat flux out the top of the model. Any hints for that?

For this you can use the d_surface_flux term, for example in a postprocess hook function [1] - something like (your region/variable/other names may differ...):

flux = problem.evaluate('d_surface_flux.2.Gamma_Top(coef.val, t)', mode='eval', copy_materials=False, verbose=False) print flux

r.

[1] http://sfepy.org/doc-devel/users_guide.html#solution-postprocessing