Thank you! It works beautifully now. (I made the boundary conditions zero to simplify the diagnosis.)

On Wednesday, December 11, 2013 3:58:54 PM UTC, Robert Cimrman wrote:
Hello Ben,

On 12/11/2013 01:41 PM, ben wrote:
> Hello,
>
> Thank you for the time you've put into this project.

Thanks!

> I've been solving a particular elliptic PDE by running SfePy on a short
> problem description file with: python simple.py problem_desc.py. Now, since
> I'd like to do this many times, with different parameters in the problem
> description file each time, I'd like to construct the problem directly, by
> creating a ProblemDefinition instance, as is done in the
> linear_elasticity.py example.
>
> (I'm working with the latest version of SfePy.) I've attempted the
> conversion, but clearly something is wrong, as I get a very large residual:
>
>> python direct.py
> sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8]
> (../rectangle_fine_quad_r.vtk)...
> sfepy: ...done in 0.72 s
> sfepy: setting up dof connectivities...
> sfepy: ...done in 0.00 s
> sfepy: updating variables...
> sfepy: ...done
> sfepy: matrix shape: (40602, 40602)
> sfepy: assembling matrix graph...
> sfepy: ...done in 0.09 s
> sfepy: matrix structural nonzeros: 723604 (4.39e-04% fill)
> sfepy: updating materials...
> sfepy:     mat_diffusion
> sfepy: ...done in 0.06 s
> sfepy: nls: iter: 0, residual: 5.209283e+08 (rel: 1.000000e+00)
> fish: Job 1, “python direct.py ” terminated by signal SIGTERM (Polite quit
> request)

The main problem was that you declared your variables to have 2 components
instead of one:

bad:

p = FieldVariable('p', 'unknown', field, mesh.dim)
s = FieldVariable('s', 'test', field, mesh.dim, primary_var_name='p')

good:

p = FieldVariable('p', 'unknown', field, 1)
s = FieldVariable('s', 'test', field, 1, primary_var_name='p')

This is related to a long-time-known quirk that the dimension has to be given
in both the field and variables.

Each term should really test this and fail gracefully - it's on my todo list
for some time...

Then you should actually use the boundary conditions you declared (and make
some nonzero(?)):

bad:

pb.time_update()

good:

pb.time_update(ebcs=Conditions([expiry_condition]))

Then the direct file uses two terms and the problem description file just one,
but I guess it's work in progress. :)

> , whereas, using the problem description file:
>
>> python ../../sfepy/simple.py problem_desc.py
> sfepy: left over: ['verbose', '__builtins__', 'pdb', '__doc__', '__name__',
> 'sys', 'data_dir', '__package__', 'refine_mesh', '_filename', 'np',
> '__file__', 'math']
> sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8]
> (../rectangle_fine_quad_r.vtk)...
> sfepy: ...done in 0.69 s
> sfepy: creating regions...
> sfepy:     Omega
> sfepy:     Expiry
> sfepy: ...done in 0.06 s
> sfepy: equation "Price":
> sfepy:  dw_diffusion.2.Omega(mat_diffusion.f, s, p)
>                    = 0
> sfepy: setting up dof connectivities...
> sfepy: ...done in 0.00 s
> sfepy: using solvers:
>                  ts: no ts
>                 nls: newton
>                  ls: ls
> sfepy: updating variables...
> sfepy: ...done
> sfepy: matrix shape: (20200, 20200)
> sfepy: assembling matrix graph...
> sfepy: ...done in 0.02 s
> sfepy: matrix structural nonzeros: 179998 (4.41e-04% fill)
> sfepy: updating materials...
> sfepy:     mat_diffusion
> sfepy: ...done in 0.06 s
> sfepy: nls: iter: 0, residual: 0.000000e+00 (rel: 0.000000e+00)

After the changes described above I get the same output as here.

>
> I've pared my problem description file and direct approach file down to the
> bare minimum. Would you mind having a quick look to see if you can spot the
> problem? Otherwise, is there some documentation that would help with the
> transition?

The linear elasticity example is described in the tutorial [1]. That's all we
have for now, as this interface is likely to change. Certainly it drifts much
more than the problem description file format.

Let us know if there are any other problems.

Cheers,
r.

[1] http://sfepy.org/doc-devel/tutorial.html#interactive-example-linear-elasticity