Source code for sfepy.tests.test_eigenvalue_solvers

import numpy as nm
import pytest

from sfepy.discrete.fem.meshio import UserMeshIO
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.solvers import Solver
import sfepy.base.testing as tst

# Mesh dimensions.
dims = nm.array([1, 1])

# Mesh resolution.
shape = [21, 21]



def mesh_hook(mesh, mode):
    """
    Generate the block mesh.
    """
    if mode == 'read':
        mesh = gen_block_mesh(dims, shape, [0, 0], name='user_block',
                              verbose=False)
        return mesh

    elif mode == 'write':
        pass


filename_mesh = UserMeshIO(mesh_hook)

regions = {
    'Omega' : 'all',
    'Gamma' : ('vertices of surface', 'facet'),
}

fields = {
    'temperature' : ('real', 1, 'Omega', 1),
}

variables = {
    't' : ('unknown field', 'temperature', 0),
    's' : ('test field',    'temperature', 't'),
}

ebcs = {
    't' : ('Gamma', {'t.0' : 0.0}),
}

integrals = {
    'i' : 2,
}

equations = {
    'eq' : 'dw_laplace.i.Omega(s, t) = 0'
}

solvers = {
    'evp0' : ('eig.scipy', {
        'method' : 'eig',
    }),
    'evp0h' : ('eig.scipy', {
        'method' : 'eigh',
    }),
    'evp0s' : ('eig.scipy', {
        'method' : 'eigs',
        'maxiter' : 100,
    }),
    'evp0sh' : ('eig.scipy', {
        'method' : 'eigsh',
        'maxiter' : 100,
    }),
    'evp1' : ('eig.sgscipy', {}),
    'evp2' : ('eig.scipy_lobpcg', {
        'i_max' : 100,
        'largest' : False,
    }),
    'evp4' : ('eig.slepc', {
        'method' : 'arnoldi',
        'problem' : 'hep',
        'i_max' : 100,
        'eps' : 1e-10,
        'which' : 'smallest_real',
    }),
    'evp5' : ('eig.matlab', {
        'maxit' : 100,
        'tol' : 1e-10,
        'which' : 'sr',
    }),
    'evp6' : ('eig.primme', {
        'maxiter' : 200,
        'tol' : 1e-10,
        'which' : 'sa',
    }),
    'evp7' : ('eig.octave', {
        'maxit' : 100,
        'tol' : 1e-10,
        'which' : 'sm',
    }),
}

eigs_expected = [nm.array([0.04904454, 0.12170685, 0.12170685,
                           0.19257998, 0.24082108]),
                 []]

can_fail = ['eig.slepc', 'eig.matlab', 'eig.octave', 'eig.primme']
can_miss = ['evp0'] # Depending on scipy version, evp0 can miss an
                    # eigenvalue.



@pytest.fixture(scope='module')
def data():
    import sys
    from sfepy.base.base import Struct
    from sfepy.discrete import Problem
    from sfepy.base.conf import ProblemConf

    conf = ProblemConf.from_dict(globals(), sys.modules[__name__])
    pb = Problem.from_conf(conf, init_solvers=False)
    pb.time_update()
    mtx = pb.equations.evaluate(mode='weak', dw_mode='matrix',
                                asm_obj=pb.mtx_a)

    return Struct(conf=conf, mtx=mtx)


def _list_eigenvalue_solvers(confs):
    d = []
    for key, val in confs.items():
        if val.kind.find('eig.') == 0:
            d.append(val)
    d.sort(key=lambda a: a.name)

    return d



def test_eigenvalue_solvers(data):
    from sfepy.base.base import IndexedStruct

    eig_confs = _list_eigenvalue_solvers(data.conf.solvers)

    all_n_eigs = [5, 0]

    ok = True
    tt = []
    for ii, n_eigs in enumerate(all_n_eigs):
        for eig_conf in eig_confs:
            tst.report(eig_conf.name)

            try:
                eig_solver = Solver.any_from_conf(eig_conf)

            except (ValueError, ImportError):
                if eig_conf.kind in can_fail:
                    continue

                else:
                    raise

            status = IndexedStruct()
            try:
                eigs, vecs = eig_solver(data.mtx, n_eigs=n_eigs,
                                        eigenvectors=True, status=status)

            except KeyboardInterrupt:
                raise

            except:
                if eig_conf.kind in can_fail:
                    continue

                else:
                    raise
            tt.append([' '.join((eig_conf.name, eig_conf.kind)),
                       status.time, n_eigs])

            tst.report(eigs)

            _ok = nm.allclose(eigs.real, eigs_expected[ii],
                              rtol=0.0, atol=1e-8)
            tt[-1].append(_ok)

            ok = ok and (_ok or (eig_conf.kind in can_fail)
                         or (eig_conf.name in can_miss))

    tt.sort(key=lambda x: x[1])
    tst.report('solution times:')
    for row in tt:
        tst.report('%.2f [s] : %s (%d) (ok: %s)'
                   % (row[1], row[0], row[2], row[3]))

    assert ok