specifying an initial guess for non-linear problems with no ebcs

[David Jessop]

Hi,

I have a non-linear form of the Poisson equation where the "diffusion" coefficient depends on the derivatives of the state variable (see image).  There are no Dirichlet-type boundary conditions (i.e. no ebcs) so I really need to specify a good initial guess at the solution. 

My questions are:
1) How are initial guesses passed to the solver?  From the documentation, it doesn't seem like initial guesses can be directly passed to the nls.newton solver as an argument (though they are an argument of the subroutine __call__).
2) In what form should the initial guess vector/array be written?  Can a function be passed?  Do numerical values have be be specified at the mesh nodes?

Thanks for any help on this matter.

D

[Robert Cimrman]

Hi David,

On 12/06/2016 03:17 PM, David Jessop wrote:
>
>
> <https://lh3.googleusercontent.com/-nqPKBodss14/WEbGCb-rgXI/AAAAAAAAA70/sQ8dedCI87g1P1S984VDxaC8wnJ26NsLwCLcB/s1600/non-linear_Poisson.png>

> Hi,
>
> I have a non-linear form of the Poisson equation where the "diffusion"
> coefficient depends on the derivatives of the state variable (see image).
> There are no Dirichlet-type boundary conditions (i.e. no ebcs) so I really
> need to specify a good initial guess at the solution.
>
> My questions are:
> 1) How are initial guesses passed to the solver? From the documentation,
> it doesn't seem like initial guesses can be directly passed to the
> nls.newton solver as an argument (though they are an argument of the
> subroutine __call__).

You are right that it is not possible to influence the initial guess of the
(non)linear solver from the problem description file. However, you can use the
interactive approach, where everything can be controlled in fine detail. See
[1], for example, where the line `vec = pb.solve()` [2] calls the solver. Its
first (optional) argument is the initial state guess (instance of State [3] -
State.vec is the numpy array holding the DOFs - you can set it in-place to
whatever you want, see below). You can also call the solvers directly (see the
parallel examples).

> 2) In what form should the initial guess vector/array be written? Can a
> function be passed? Do numerical values have be be specified at the mesh
> nodes?

For the Lagrange basis (the default), you can use a function of nodal
coordinates - just evaluate what you need in the coordinates returned by
Field.get_coor().

Or, more general, you could misuse the initial conditions facilities for time
dependent problems - create the initial conditions as in [4], and call
State.apply_ic().

Let me know if you need a more detailed help.

r.

[1]
http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_interactive.html
[2]
http://sfepy.org/doc-devel/src/sfepy/discrete/problem.html?highlight=problem.solve#sfepy.discrete.problem.Problem.solve
[3]
http://sfepy.org/doc-devel/src/sfepy/discrete/state.html?#sfepy.discrete.state.State
[4] http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html

[David Jessop]

Hi Robert,

Thanks for those hints.  I'll try them out and let you know how it goes.

Regards,
David

[David Jessop]

Hi,

I'm afraid that I'll need some specific help.  I've tried recasting my code to run in interactive mode but am getting the following error:

ValueError: wrong arguments shapes for "dw_laplace.1.Omega(coef.val, s, T)" term! (see above)

I can't see how to resolve this issue from the linear_elasticity_interactive.py case. 

I've written the respective variables as:

coef = Material('coef', val=[2.0])
T = FieldVariable('T', 'unknown', field)
s = FieldVariable('s', 'test', field, primary_var_name='T')

I've also tried writing:

coef = Material('coef', values={'.val' : 0.0})
as per its2D_interactive.py, but with the same effect.

Thanks,

[David Jessop]

Here's the complete problem description file.

On Tuesday, 6 December 2016 15:26:44 UTC+1, David Jessop wrote:

[Robert Cimrman]

Hi David,

The main problem is having 'vector' in the field creation (line 119) - this
means, that your variables are vectors, and not scalars, as required by the
Laplace term. Try replacing that with 'scalar' or 1.

As for materials, you can use simply:

coef = Material('coef', val=2.0)

etc.

r.

On 12/16/2016 10:46 AM, David Jessop wrote:
> Here's the complete problem description file.
>
> On Tuesday, 6 December 2016 15:26:44 UTC+1, David Jessop wrote:
>>
>>

>> <https://lh3.googleusercontent.com/-nqPKBodss14/WEbGCb-rgXI/AAAAAAAAA70/sQ8dedCI87g1P1S984VDxaC8wnJ26NsLwCLcB/s1600/non-linear_Poisson.png>

[David Jessop]

I'm getting some very strange behaviour in my "interactive" solution.  For the time being I'm just solving a linear Poisson equation on a square region with a constant source term everywhere and Dirichelet BCs on the left and right boundaries.  I'm set up and solved the same problem using a problem description file with the simple.py routine (see myPoisson_Soln.pdf) yet the interactive form has massive oscillations (see myPoissonInteractive_solution.

pdf) in the solution.  The error seems to be in that the matrix for the interactive case is too large (1599x1599 elements for a domain of 21x21 cells, myPoisson.py gives 399x399 elements for the same number of cells in the domain) and so the inversion is singular, yet I can't find where the error in my code could be.  Would someone please mind pointing it out to me?

Thanks.
D

Output of ./simple.py myPoisson.py:

sfepy: left over: ['verbose', '__builtins__', 'n_step', 'dims', 'shape', '__file__', '__name__', 't1', 'center', 'UserMeshIO', 'gen_block_mesh', 't0', '__package__', 'output_dir', '_filename', 'np', 'output', '__doc__', 'mesh_hook']
sfepy: reading mesh [user] (function:mesh_hook)...
sfepy: ...done in 0.00 s
sfepy: creating regions...
sfepy:     Right
sfepy:     Top
sfepy:     Bottom
sfepy:     Omega
sfepy:     Left
sfepy: ...done in 0.00 s
sfepy: equation "Temperature":
sfepy:
           dw_laplace.i.Omega(cond.val, s, T)
         - dw_surface_integrate.2.Top(insulated.val, s)
         - dw_volume_lvf.2.Omega(G.val, s)

sfepy: using solvers:
                ts: ts
               nls: newton
                ls: ls
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: matrix shape: (399, 399)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 3355 (2.11e-02% fill)
sfepy: ====== time 0.000000e+00 (step 1 of 2) =====
sfepy: updating materials...
sfepy:     G
sfepy:     cond
sfepy:     insulated
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 2.530862e+01 (rel: 1.000000e+00)
sfepy:   rezidual:    0.00 [s]
sfepy:      solve:    0.00 [s]
sfepy:     matrix:    0.00 [s]
sfepy: nls: iter: 1, residual: 6.515061e-14 (rel: 2.574246e-15)
sfepy: ====== time 1.000000e-01 (step 2 of 2) =====
sfepy: updating variables...
sfepy: ...done
sfepy: updating materials...
sfepy:     G
sfepy:     cond
sfepy:     insulated
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 6.515061e-14 (rel: 1.000000e+00)

Output of python myPoissonInteractive.py:

Enter code here...sfepy: saving regions as groups...
sfepy:   Omega
sfepy:   Left
sfepy:   Right
sfepy:   Bottom
sfepy:   Top
sfepy: ...done
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: matrix shape: (1599, 1599)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 24311 (9.51e-03% fill)
sfepy: updating materials...
sfepy:     cond
sfepy:     insulated
sfepy:     G
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 6.169742e+01 (rel: 1.000000e+00)
sfepy:   rezidual:    0.00 [s]
sfepy:      solve:    0.01 [s]
sfepy:     matrix:    0.00 [s]
sfepy: warning: linear system solution precision is lower
sfepy: then the value set in solver options! (err = 2.856021e+01 < 1.000000e-10)
sfepy: nls: iter: 1, residual: 2.946994e+01 (rel: 4.776527e-01)
IndexedStruct
  condition:
    1
  err:
    29.4699421316
  err0:
    61.6974210411
  n_iter:
    1
  time_stats:
    dict with keys: ['rezidual', 'solve', 'matrix']

[Robert Cimrman]

Hi David,

note "approx_order=2" in the field definition in your "interactive" script
(line 77). Use 1 to have the smaller matrix as in the problem description file.
Or, if you do want to use the bi-quadratic order, increase the order of
numerical integration (line 93) to two. Currently you are under-integrating and
that is why you keep getting a singular matrix.

r.

[David Jessop]

Hi Robert.

I've now got the interactive form working nicely, with an initial guess passed to the problem solver as you describe in your first reply.  I'm having some issues with passing function values to the Material properties, but that might be best as a separate thread.

Thanks for all your help so far!
David

[Robert Cimrman]

On 12/21/2016 08:14 PM, David Jessop wrote:
> Hi Robert.
>
> I've now got the interactive form working nicely, with an initial guess
> passed to the problem solver as you describe in your first reply. I'm
> having some issues with passing function values to the Material properties,
> but that might be best as a separate thread.

OK, you may also want to check [1].

Cheers,
r.
[1] http://sfepy.org/doc-devel/users_guide.html#functions