import pytest
from sfepy import data_dir
import sfepy.base.testing as tst
filename_mesh = data_dir + '/meshes/3d/special/cube_cylinder.mesh'
if 0:
from sfepy.discrete.fem.utils import refine_mesh
refinement_level = 1
filename_mesh = refine_mesh(filename_mesh, refinement_level)
material_2 = {
'name' : 'coef',
'values' : {'val' : 1.0},
}
field_1 = {
'name' : 'temperature',
'dtype' : 'real',
'shape' : (1,),
'region' : 'Omega',
'approx_order' : 1,
}
variables = {
't' : ('unknown field', 'temperature', 0),
's' : ('test field', 'temperature', 't'),
}
regions = {
'Omega' : 'all',
'Gamma_Left' : ('vertices in (x < 0.0001)', 'facet'),
'Gamma_Right' : ('vertices in (x > 0.999)', 'facet'),
}
ebcs = {
't1' : ('Gamma_Left', {'t.0' : 2.0}),
't2' : ('Gamma_Right', {'t.0' : -2.0}),
}
integral_1 = {
'name' : 'i',
'order' : 1,
}
equations = {
'Temperature' : """dw_laplace.i.Omega(coef.val, s, t) = 0"""
}
[docs]
class DiagPC(object):
"""
Diagonal (Jacobi) preconditioner.
Equivalent to setting `'precond' : 'jacobi'`.
"""
[docs]
def setUp(self, pc):
A = pc.getOperators()[0]
self.idiag = 1.0 / A.getDiagonal()
[docs]
def apply(self, pc, x, y):
y.pointwiseMult(x, self.idiag)
[docs]
def setup_petsc_precond(mtx, problem):
return DiagPC()
solvers = {
'd00' : ('ls.scipy_direct',
{}
),
'd01' : ('ls.scipy_direct',
{'method' : 'umfpack',
'warn' : True,}
),
'd02' : ('ls.scipy_direct',
{'method' : 'superlu',
'warn' : True,}
),
'd10' : ('ls.mumps', {}),
'i00' : ('ls.pyamg',
{'method' : 'ruge_stuben_solver',
'accel' : 'cg',
'eps_r' : 1e-12,
'method:max_levels' : 5,
'solve:cycle' : 'V',}
),
'i01' : ('ls.pyamg',
{'method' : 'smoothed_aggregation_solver',
'accel' : 'cg',
'eps_r' : 1e-12,}
),
'i02' : ('ls.pyamg_krylov',
{'method' : 'cg',
'eps_r' : 1e-12,
'i_max' : 1000,}
),
'i10' : ('ls.petsc',
{'method' : 'cg', # ksp_type
'precond' : 'none', # pc_type
'eps_a' : 1e-12, # abstol
'eps_r' : 1e-12, # rtol
'i_max' : 1000,} # maxits
),
'i11' : ('ls.petsc',
{'method' : 'cg', # ksp_type
'precond' : 'python', # just for output (unused)
'setup_precond' : setup_petsc_precond, # user-defined pc
'eps_a' : 1e-12, # abstol
'eps_r' : 1e-12, # rtol
'i_max' : 1000,} # maxits
),
'i12' : ('ls.petsc',
{'method' : 'cg', # ksp_type
'precond' : 'jacobi', # pc_type
'eps_a' : 1e-12, # abstol
'eps_r' : 1e-12, # rtol
'i_max' : 1000,} # maxits
),
'i13' : ('ls.petsc',
{'method' : 'cg', # ksp_type
'precond' : 'icc', # pc_type
'eps_a' : 1e-12, # abstol
'eps_r' : 1e-12, # rtol
'i_max' : 1000,} # maxits
),
'i20' : ('ls.scipy_iterative',
{'method' : 'cg',
'i_max' : 1000,
'eps_a' : 1e-12,
'eps_r' : 1e-12,}
),
'i21' : ('ls.scipy_iterative',
{'method' : 'bicgstab',
'i_max' : 1000,
'eps_a' : 1e-12,
'eps_r' : 1e-12,}
),
'i22' : ('ls.scipy_iterative',
{'method' : 'qmr',
'i_max' : 1000,
'eps_a' : 1e-12,
'eps_r' : 1e-12,}
),
'newton' : ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-10,
}),
}
options = {
'nls' : 'newton',
}
can_fail = ['ls.pyamg', 'ls.pyamg_krylov', 'ls.petsc', 'ls.mumps',
'ls.scipy_direct']
[docs]
@pytest.fixture(scope='module')
def problem():
import sys
from sfepy.discrete import Problem
from sfepy.base.conf import ProblemConf
conf = ProblemConf.from_dict(globals(), sys.modules[__name__])
problem = Problem.from_conf(conf)
problem.time_update()
return problem
def _list_linear_solvers(confs):
d = []
for key, val in confs.items():
if val.kind.find('ls.') == 0:
d.append(val)
d.sort(key=lambda a: a.name)
return d
[docs]
def test_solvers(problem, output_dir):
from sfepy.base.base import IndexedStruct
import os.path as op
solver_confs = _list_linear_solvers(problem.solver_confs)
ok = True
tt = []
for solver_conf in solver_confs:
method = solver_conf.get('method', '')
precond = solver_conf.get('precond', '')
name = ' '.join((solver_conf.name, solver_conf.kind,
method, precond)).rstrip()
tst.report(name)
tst.report('matrix size:', problem.mtx_a.shape)
tst.report(' nnz:', problem.mtx_a.nnz)
status = IndexedStruct()
try:
problem.init_solvers(status=status,
ls_conf=solver_conf,
force=True)
state = problem.solve(save_results=False)
failed = status.nls_status.condition != 0
except Exception as aux:
failed = True
status = None
exc = aux
ok = ok and ((not failed) or (solver_conf.kind in can_fail))
if status is not None:
status = status.nls_status
for kv in status.time_stats.items():
tst.report('%10s: %7.2f [s]' % kv)
tst.report('condition: %d, err0: %.3e, err: %.3e'
% (status.condition, status.err0, status.err))
tt.append([name,
status.time_stats['solve'],
status.ls_n_iter,
status.err])
aux = name.replace(' ', '_')
fname = op.join(output_dir, 'test_linear_solvers_%s.vtk') % aux
problem.save_state(fname, state)
else:
tst.report('solver failed:')
tst.report(exc)
tt.append([name, -1, 1e10, 1e10])
tt.sort(key=lambda a: a[1])
tst.report('solution times / numbers of iterations (residual norms):')
for row in tt:
tst.report('%.2f [s] / % 4d' % (row[1], row[2]),
'(%.3e)' % row[3], ':', row[0])
assert ok
[docs]
def test_ls_reuse(problem):
import numpy as nm
from sfepy.solvers import Solver
problem.init_solvers(ls_conf=problem.solver_confs['d00'])
nls = problem.get_nls()
state0 = problem.get_initial_state()
state0.apply_ebc()
vec0 = state0.get_state(problem.active_only)
problem.update_materials()
rhs = nls.fun(vec0)
mtx = nls.fun_grad(vec0)
ok = True
for name in ['i12', 'i01']:
solver_conf = problem.solver_confs[name]
method = solver_conf.get('method', '')
precond = solver_conf.get('precond', '')
name = ' '.join((solver_conf.name, solver_conf.kind,
method, precond)).rstrip()
tst.report(name)
try:
ls = Solver.any_from_conf(solver_conf)
except:
tst.report('skipped!')
continue
conf = ls.conf.copy()
conf.force_reuse = True
sol00 = ls(rhs, mtx=mtx, conf=conf)
digest00 = ls.mtx_digest
sol0 = ls(rhs, mtx=mtx)
digest0 = ls.mtx_digest
sol1 = ls(rhs, mtx=2*mtx, conf=conf)
digest1 = ls.mtx_digest
sol2 = ls(rhs, mtx=2*mtx)
digest2 = ls.mtx_digest
ls(rhs, mtx=2*mtx)
digest3 = ls.mtx_digest
_ok = digest00 != digest0
tst.report(digest00, '!=', digest0, ':', _ok); ok = ok and _ok
_ok = digest0 == digest1
tst.report(digest0, '==', digest1, ':', _ok); ok = ok and _ok
_ok = digest1 != digest2
tst.report(digest1, '!=', digest2, ':', _ok); ok = ok and _ok
_ok = digest2[1] == digest3[1]
tst.report(digest2[1], '==', digest3[1], ':', _ok); ok = ok and _ok
_ok = nm.allclose(sol00, sol0, atol=1e-12, rtol=0.0)
tst.report('sol00 == sol0:', _ok); ok = ok and _ok
_ok = nm.allclose(sol0, sol1, atol=1e-12, rtol=0.0)
tst.report('sol0 == sol1:', _ok); ok = ok and _ok
_ok = nm.allclose(sol0, 2 * sol2, atol=1e-12, rtol=0.0)
tst.report('sol0 == 2 * sol2:', _ok); ok = ok and _ok
assert ok