45 views

Skip to first unread message

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

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

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

to sfepy...@googlegroups.com

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>

(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

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
>

> 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__).

(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?

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

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

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

Regards,

David

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:

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

I've written the respective variables as:

I've also tried writing:

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

to sfepy...@googlegroups.com

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>

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

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:

Output of python myPoissonInteractive.py:

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']

Dec 21, 2016, 8:09:03 AM12/21/16

to sfepy...@googlegroups.com

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.

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.

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

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

Dec 21, 2016, 2:25:45 PM12/21/16

to sfepy...@googlegroups.com

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].
> 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.

Cheers,

r.

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

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu