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