sfepy/sfepy — Issue #43: Schrodinger equation solver doesn't work in the latest sfepy
Author: rc
Created: 2012-01-30 10:24:39+00:00
State: open

Schrodinger equation solver doesn't work in the latest sfepy:

$ ./schroedinger.py -s
sfe: warning: other missing: ['functions', 'modules', 'epbc_[0-9]+|epbcs',
'lcbc_[0-9]+|lcbcs', 'nbc_[0-9]+|nbcs']
sfe: warning: left over: ['funV', 'solver_2']
sfe: reading mesh...
nodes: 100% |############################################| Time:
00:00:00
elements: 100% |#########################################| Time:
00:00:00
sfe: ...done in 1.20 s
sfe: setting up domain edges...
sfe: ...done in 0.89 s
sfe: setting up domain faces...
sfe: ...done in 0.66 s
sfe: creating regions...
sfe: leaf Omega region_1000
sfe: leaf Surface region_2
sfe: ...done in 1.04 s
sfe: equation "rhs":
sfe: dw_mass_scalar.i1.Omega( v, Psi )
sfe: equation "lhs":
sfe: dw_laplace.i1.Omega( m.val, v, Psi )
+ dw_mass_scalar_variable.i1.Omega( matV.V, v, Psi )
sfe: describing geometries...
sfe: ...done in 0.16 s
sfe: setting up dof connectivities...
sfe: ...done in 0.00 s
sfe: matrix shape: (7964, 7964)
sfe: assembling matrix graph...
sfe: ...done in 0.10 s
sfe: matrix structural nonzeros: 117024 (1.85e-03% fill)
sfe: updating materials...
sfe: m
sfe: matV
sfe: ...done in 0.01 s
sfe: assembling lhs...
sfe: setting up dof connectivities...
sfe: ...done in 0.00 s
Traceback (most recent call last):
File "./schroedinger.py", line 317, in <module>
main()
File "./schroedinger.py", line 312, in main
evp = solveEigenProblem1( conf, options )
File "./schroedinger.py", line 211, in solveEigenProblem1
dwMode = 'matrix', tangentMatrix = pb.mtxA )
File "/home/ondra/sfepy/sfe/fem/evaluate.py", line 186, in evalTermOP
chunkSize = problem.domain.shape.nEl, **kwargs )
File "/home/ondra/sfepy/sfe/fem/evaluate.py", line 237, in evalTerm
chunkSize, groupCanFail = False, **kwargs )
File "/home/ondra/sfepy/sfe/fem/evaluate.py", line 173, in assembleMatrix
**args ):
File "/home/ondra/sfepy/sfe/terms/termsMass.py", line 129, in __call__
fargs = self.buildCFunArgs( state, ap, vg, **kwargs )
File "/home/ondra/sfepy/sfe/terms/termsMass.py", line 172, in buildCFunArgs
modeIn = 'vertex' )
File "/home/ondra/sfepy/sfe/terms/cache.py", line 178, in __call__
self.update( key, groupIndx, ih, **kwargs )
File "/home/ondra/sfepy/sfe/terms/cachesBasic.py", line 350, in update
matQP = self.data[key][ckey][ih].reshape( vshape )
ValueError: total size of new array must be unchanged

Migrated from http://code.google.com/p/sfepy/issues/detail?id=39

earlier comments

ondrej.certik said, at 2008-06-11T15:25:29.000Z:

Whatever the problem is, let's make sfepy more robust to raise a meaningful exception
at the right place. + a test for it.

ondrej.certik said, at 2008-06-11T23:00:07.000Z:

I used "hg bisect" and here is what I got:

The first bad revision is:
changeset: 180:c73e7db416dd
user: Robert Cimrman
date: Tue May 06 18:03:32 2008 +0200
summary: fixed MatInQPDataCache to work with diffusion tensor in vertices, in

well, I think it doesn't need any more comments. :) Do you know how to fix it?

ondrej.certik said, at 2008-06-11T23:12:15.000Z:

This patch makes it run with the latest sfepy:

diff --git a/sfe/terms/cachesBasic.py b/sfe/terms/cachesBasic.py
--- a/sfe/terms/cachesBasic.py
+++ b/sfe/terms/cachesBasic.py
@@ -341,7 +341,7 @@ class MatInQPDataCache( DataCache ):

# dq_state_in_qp() works for vectors -> make a view of
# shape (nEl, nQP, nRow * nCol, 1).
- vshape = shape[0:2] + (nm.prod( mat.shape[1:] ), 1)
+ vshape = shape[0:2] + (mat.shape[1], 1)
## print self
## print self.shape, ckey
## print vshape

Could you please explain me (or even better in the comments in the sources) what is
this supposed to do? And then let's write a test for this, as clearly it can be
broken sometimes.

robert.cimrman said, at 2008-06-12T11:20:30.000Z:

I will check it. The problem is that the tests do not cover special solvers as
schroedinger.py - I will make a test veryfing if it runs.

robert.cimrman said, at 2008-06-12T11:58:29.000Z:

Fixed (the input files).

Explanation: MatInQPDataCache is a class that can, for example, interpolate material
data from vertices into Gauss points of the elements. The changeset 180 broadened the
possible input data shapes, but I forgot to fix also the schreodinger input files.

Test it please.

ondrej.certik said, at 2008-06-12T13:29:10.000Z:

Works well now. Thanks!