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

45 views

### David Jessop

Dec 6, 2016, 9:26:44 AM12/6/16
to sfepy-devel

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

Dec 7, 2016, 2:59:27 AM12/7/16
Hi David,

On 12/06/2016 03:17 PM, David Jessop wrote:
>
>
> 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

Dec 9, 2016, 10:21:31 AM12/9/16
to sfepy-devel
Hi Robert,

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

Regards,
David

### David Jessop

Dec 16, 2016, 4:41:40 AM12/16/16
to sfepy-devel
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

Dec 16, 2016, 4:46:32 AM12/16/16
to sfepy-devel
Here's the complete problem description file.

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

### Robert Cimrman

Dec 16, 2016, 4:54:39 AM12/16/16
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:
>>
>>
Message has been deleted

### David Jessop

Dec 21, 2016, 7:59:43 AM12/21/16
to sfepy-devel
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 ssfepy: creating regions...sfepy:     Rightsfepy:     Topsfepy:     Bottomsfepy:     Omegasfepy:     Leftsfepy: ...done in 0.00 ssfepy: 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: lssfepy: updating variables...sfepy: ...donesfepy: setting up dof connectivities...sfepy: ...done in 0.00 ssfepy: matrix shape: (399, 399)sfepy: assembling matrix graph...sfepy: ...done in 0.00 ssfepy: matrix structural nonzeros: 3355 (2.11e-02% fill)sfepy: ====== time 0.000000e+00 (step 1 of 2) =====sfepy: updating materials...sfepy:     Gsfepy:     condsfepy:     insulatedsfepy: ...done in 0.00 ssfepy: 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: ...donesfepy: updating materials...sfepy:     Gsfepy:     condsfepy:     insulatedsfepy: ...done in 0.00 ssfepy: 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:   Omegasfepy:   Leftsfepy:   Rightsfepy:   Bottomsfepy:   Topsfepy: ...donesfepy: updating variables...sfepy: ...donesfepy: setting up dof connectivities...sfepy: ...done in 0.00 ssfepy: matrix shape: (1599, 1599)sfepy: assembling matrix graph...sfepy: ...done in 0.00 ssfepy: matrix structural nonzeros: 24311 (9.51e-03% fill)sfepy: updating materials...sfepy:     condsfepy:     insulatedsfepy:     Gsfepy: ...done in 0.00 ssfepy: 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 lowersfepy: 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']`

myPoisson.py
myPoissonInteractive.py
myPoissonInteractive_solution-reduced.pdf
myPoisson_Soln-reduced.pdf

### Robert Cimrman

Dec 21, 2016, 8:09:03 AM12/21/16
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

Dec 21, 2016, 2:14:08 PM12/21/16
to sfepy-devel
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

Dec 21, 2016, 2:25:45 PM12/21/16
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