specifying an initial guess for non-linear problems with no ebcs
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:
- 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
Hi David,
On 12/06/2016 03:17 PM, David Jessop wrote:
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:
- 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).
- 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_interac... [2] http://sfepy.org/doc-devel/src/sfepy/discrete/problem.html?highlight=problem... [3] http://sfepy.org/doc-devel/src/sfepy/discrete/state.html?#sfepy.discrete.sta... [4] http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html
Hi Robert,
Thanks for those hints. I'll try them out and let you know how it goes.
Regards, David
On Tuesday, 6 December 2016 15:26:44 UTC+1, David Jessop wrote:
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:
- 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
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
On Tuesday, 6 December 2016 15:26:44 UTC+1, David Jessop wrote:
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:
- 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
Here's the complete problem description file.
On Tuesday, 6 December 2016 15:26:44 UTC+1, David Jessop wrote:
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:
- 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
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:
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:
- 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
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']
On Friday, 16 December 2016 10:54:39 UTC+1, Robert Cimrman wrote:
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/sQ8de...
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:
- 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
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.
On 12/21/2016 01:59 PM, David Jessop wrote:
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']
On Friday, 16 December 2016 10:54:39 UTC+1, Robert Cimrman wrote:
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/sQ8de...
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:
- 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
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
On Wednesday, 21 December 2016 14:09:03 UTC+1, Robert Cimrman wrote:
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.
On 12/21/2016 01:59 PM, David Jessop wrote:
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']
On Friday, 16 December 2016 10:54:39 UTC+1, Robert Cimrman wrote:
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/sQ8de...
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:
- 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
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
Thanks for all your help so far! David
participants (2)
-
David Jessop -
Robert Cimrman