Python solvers.Solver类代码示例

    def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None, mtx=None, **kwargs):
        ls_conf = get_default(ls_conf, self.ls_conf, "you must set linear solver!")

        nls_conf = get_default(nls_conf, self.nls_conf, "you must set nonlinear solver!")

        ev = self.get_evaluator(**kwargs)
        ls = Solver.any_from_conf(ls_conf, mtx=mtx)
        nls = Solver.any_from_conf(nls_conf, evaluator=ev, lin_solver=ls, status=nls_status)
        self.solvers = Struct(name="solvers", ls=ls, nls=nls)

    def test_eigenvalue_solvers(self):
        eig_confs = self._list_eigenvalue_solvers(self.conf.solvers)

        n_eigs = 5

        ok = True
        tt = []
        for eig_conf in eig_confs:
            self.report(eig_conf.name)

            eig_solver = Solver.any_from_conf(eig_conf)

            t0 = time.clock()
            eigs, vecs = eig_solver(self.mtx, n_eigs=n_eigs, eigenvectors=True)
            tt.append((' '.join((eig_conf.name, eig_conf.kind)),
                       time.clock() - t0))

            self.report(eigs)

            _ok = nm.allclose(eigs.real, eigs_expected, rtol=0.0, atol=1e-8)
            ok = ok and (_ok or (eig_conf.kind in self.can_fail))

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

        return ok

    def presolve(self, mtx):
        """Prepare A^{-1} B^T for the Schur complement."""

        mtx_a = mtx['A']
        mtx_bt = mtx['BT']
        output('full A size: %.3f MB' % (8.0 * nm.prod(mtx_a.shape) / 1e6))
        output('full B size: %.3f MB' % (8.0 * nm.prod(mtx_bt.shape) / 1e6))

        ls = Solver.any_from_conf(self.problem.ls_conf,
                                   presolve=True, mtx=mtx_a)
        if self.mode == 'explicit':
            tt = time.clock()
            mtx_aibt = nm.zeros(mtx_bt.shape, dtype=mtx_bt.dtype)
            for ic in xrange(mtx_bt.shape[1]):
                mtx_aibt[:,ic] = ls(mtx_bt[:,ic].toarray().squeeze())
            output('mtx_aibt: %.2f s' % (time.clock() - tt))
            action_aibt = MatrixAction.from_array(mtx_aibt)
        else:
            ##
            # c: 30.08.2007, r: 13.02.2008
            def fun_aibt(vec):
                # Fix me for sparse mtx_bt...
                rhs = sc.dot(mtx_bt, vec)
                out = ls(rhs)
                return out
            action_aibt = MatrixAction.from_function(fun_aibt,
                                                    (mtx_a.shape[0],
                                                     mtx_bt.shape[1]),
                                                    nm.float64)
        mtx['action_aibt'] = action_aibt

def solve_optimize(conf, options):
    opts = conf.options
    trunk = io.get_trunk(conf.filename_mesh)
    data = {}

    dpb = ProblemDefinition.from_conf(conf, init_equations=False)
    equations = getattr(conf, "_".join(("equations_direct", opts.problem)))

    dpb.set_equations(equations)

    dpb.name = "direct"
    dpb.time_update(None)

    apb = dpb.copy("adjoint")
    equations = getattr(conf, "_".join(("equations_adjoint", opts.problem, opts.objective_function)))

    apb.set_equations(equations)
    apb.time_update(None)
    apb.ebcs.zero_dofs()
    apb.update_equations(None, ebcs=apb.ebcs)

    ls_conf = dpb.get_solver_conf(opts.ls)
    dnls_conf = dpb.get_solver_conf(opts.nls_direct)
    anls_conf = dpb.get_solver_conf(opts.nls_adjoint)
    opt_conf = dpb.get_solver_conf(opts.optimizer)

    dpb.init_solvers(ls_conf=ls_conf, nls_conf=dnls_conf)

    apb.init_solvers(ls_conf=ls_conf, nls_conf=anls_conf)

    shape_opt = so.ShapeOptimFlowCase.from_conf(conf, dpb, apb)
    design0 = shape_opt.dsg_vars.val
    shape_opt.cache = Struct(design=design0 + 100, state=None, i_mesh=-1)

    opt_status = IndexedStruct()
    optimizer = Solver.any_from_conf(
        opt_conf, obj_fun=so.obj_fun, obj_fun_grad=so.obj_fun_grad, status=opt_status, obj_args=(shape_opt, opts)
    )

    ##
    # State problem solution for the initial design.
    vec_dp0 = so.solve_problem_for_design(dpb, design0, shape_opt, opts)

    dpb.save_state(trunk + "_direct_initial.vtk", vec_dp0)

    ##
    # Optimize.
    des = optimizer(design0)
    print opt_status

    ##
    # Save final state (for "optimal" design).
    dpb.domain.mesh.write(trunk + "_opt.mesh", io="auto")
    dpb.save_state(trunk + "_direct_current.vtk", shape_opt.cache.state)

    print des

    def __init__(self, problem, options):
        self.mtx_mass = problem.evaluate(options.mass, mode='weak',
                                         auto_init=True, dw_mode='matrix')

        if options.lumped:
            raise NotImplementedError

        else:
            # Initialize solvers (and possibly presolve the matrix).
            self.ls = Solver.any_from_conf(problem.ls_conf, mtx=self.mtx_mass,
                                           presolve=True)

    def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None,
                     force=False):
        """
        Create and initialize solver instances.

        Parameters
        ----------
        nls_status : dict-like, IndexedStruct, optional
            The user-supplied object to hold nonlinear solver convergence
            statistics.
        ls_conf : Struct, optional
            The linear solver options.
        nls_conf : Struct, optional
            The nonlinear solver options.
        force : bool
            If True, re-create the solver instances even if they already exist
            in `self.nls` attribute.
        """
        if (self.nls is None) or force:
            ls_conf = get_default(ls_conf, self.ls_conf,
                                  'you must set linear solver!')

            nls_conf = get_default(nls_conf, self.nls_conf,
                                   'you must set nonlinear solver!')

            ls = Solver.any_from_conf(ls_conf, problem=self)

            ev = self.get_evaluator()

            if self.conf.options.get('ulf', False):
                self.nls_iter_hook = ev.new_ulf_iteration

            nls = Solver.any_from_conf(nls_conf, fun=ev.eval_residual,
                                       fun_grad=ev.eval_tangent_matrix,
                                       lin_solver=ls,
                                       iter_hook=self.nls_iter_hook,
                                       status=nls_status, problem=self)

            self.set_solver(nls)

def eig( mtx_a, mtx_b = None, num = None, eigenvectors = True,
         return_time = None, method = 'eig.scipy', **ckwargs ):

    kwargs = {'name' : 'aux', 'kind' : method}
    kwargs.update( ckwargs )
    conf = Struct( **kwargs )
    solver = Solver.any_from_conf( conf )

    status = {}
    out = solver( mtx_a, mtx_b, num, eigenvectors, status )
    if return_time is not None:
        return_time[0] = status['time']
        
    return out

    def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None,
                     mtx=None, presolve=False):
        """Create and initialize solvers."""
        ls_conf = get_default( ls_conf, self.ls_conf,
                               'you must set linear solver!' )

        nls_conf = get_default( nls_conf, self.nls_conf,
                              'you must set nonlinear solver!' )

        if presolve:
            tt = time.clock()
        if get_default_attr(ls_conf, 'needs_problem_instance', False):
            extra_args = {'problem' : self}
        else:
            extra_args = {}

        ls = Solver.any_from_conf(ls_conf, mtx=mtx, presolve=presolve,
                                  **extra_args)
        if presolve:
            tt = time.clock() - tt
            output('presolve: %.2f [s]' % tt)

        if get_default_attr(nls_conf, 'needs_problem_instance', False):
            extra_args = {'problem' : self}
        else:
            extra_args = {}
        ev = self.get_evaluator()

        if get_default_attr(self.conf.options, 'ulf', False):
            self.nls_iter_hook = ev.new_ulf_iteration

        nls = Solver.any_from_conf(nls_conf, fun=ev.eval_residual,
                                   fun_grad=ev.eval_tangent_matrix,
                                   lin_solver=ls, iter_hook=self.nls_iter_hook,
                                   status=nls_status, **extra_args)

        self.solvers = Struct( name = 'solvers', ls = ls, nls = nls )

    def prepare_matrices(self, problem):
        """
        A = K + B^T D^{-1} B
        """
        equations = problem.equations
        mtx = equations.eval_tangent_matrices(None, problem.mtx_a,
                                              by_blocks=True)

        ls = Solver.any_from_conf(problem.ls_conf,
                                  presolve=True, mtx=mtx['D'])

        mtx_b, mtx_m = mtx['B'], mtx['M']
        mtx_dib = nm.empty(mtx_b.shape, dtype=mtx_b.dtype)
        for ic in xrange(mtx_b.shape[1]):
            mtx_dib[:,ic] = ls(mtx_b[:,ic].toarray().squeeze())
        mtx_a = mtx['K'] + mtx_b.T * mtx_dib

        return mtx_a, mtx_m, mtx_dib

    def solve_eigen_problem(self):
        opts = self.app_options
        pb = self.problem

        pb.set_equations(pb.conf.equations)
        pb.time_update()

        output('assembling lhs...')
        tt = time.clock()
        mtx_a = pb.evaluate(pb.conf.equations['lhs'], mode='weak',
                            auto_init=True, dw_mode='matrix')
        output('...done in %.2f s' % (time.clock() - tt))

        if 'rhs' in pb.conf.equations:
            output('assembling rhs...')
            tt = time.clock()
            mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
                                dw_mode='matrix')
            output('...done in %.2f s' % (time.clock() - tt))

        else:
            mtx_b = None

        _n_eigs = get_default(opts.n_eigs, mtx_a.shape[0])

        output('solving eigenvalue problem for {} values...'.format(_n_eigs))
        eig = Solver.any_from_conf(pb.get_solver_conf(opts.evps))
        if opts.eigs_only:
            eigs = eig(mtx_a, mtx_b, opts.n_eigs, eigenvectors=False)
            svecs = None

        else:
            eigs, svecs = eig(mtx_a, mtx_b, opts.n_eigs, eigenvectors=True)

        output('...done')

        vecs = self.make_full(svecs)
        self.save_results(eigs, vecs)

        return Struct(pb=pb, eigs=eigs, vecs=vecs)

    def test_eigenvalue_solvers(self):
        eig_confs = self._list_eigenvalue_solvers(self.conf.solvers)

        n_eigs = 5

        ok = True
        tt = []
        for eig_conf in eig_confs:
            self.report(eig_conf.name)

            try:
                eig_solver = Solver.any_from_conf(eig_conf)

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

                else:
                    raise

            t0 = time.clock()
            eigs, vecs = eig_solver(self.mtx, n_eigs=n_eigs, eigenvectors=True)
            tt.append([' '.join((eig_conf.name, eig_conf.kind)),
                       time.clock() - t0])

            self.report(eigs)

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

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

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

        return ok

    def solve_eigen_problem_n( self ):
        opts = self.app_options
        pb = self.problem

        dim = pb.domain.mesh.dim

        pb.set_equations( pb.conf.equations )
        pb.select_bcs( ebc_names = ['ZeroSurface'] )

        output( 'assembling rhs...' )
        tt = time.clock()
        mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
                            auto_init=True, dw_mode='matrix')
        output( '...done in %.2f s' % (time.clock() - tt) )
        assert_( nm.alltrue( nm.isfinite( mtx_b.data ) ) )

        ## mtx_b.save( 'b.txt', format='%d %d %.12f\n' )

        aux = pb.create_evaluable(pb.conf.equations['lhs'], mode='weak',
                                  dw_mode='matrix')
        mtx_a_equations, mtx_a_variables = aux

        if self.options.plot:
            log_conf = {
                'is_plot' : True,
                'aggregate' : 1,
                'yscales' : ['linear', 'log'],
            }
        else:
            log_conf = {
                'is_plot' : False,
            }
        log =  Log.from_conf( log_conf, ([r'$|F(x)|$'], [r'$|F(x)-x|$']) )

        file_output = Output('', opts.log_filename, combined = True)

        eig_conf = pb.get_solver_conf( opts.eigen_solver )
        eig_solver = Solver.any_from_conf( eig_conf )

        # Just to get the shape. Assumes one element group only!!!
        v_hxc_qp = pb.evaluate('dq_state_in_volume_qp.i1.Omega(Psi)')
        v_hxc_qp.fill(0.0)
        self.qp_shape = v_hxc_qp.shape
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)

        self.norm_v_hxc0 = nla.norm(vec_v_hxc)
        self.itercount = 0
        aux = wrap_function(self.iterate,
                            (eig_solver,
                             mtx_a_equations, mtx_a_variables,
                             mtx_b, log, file_output))
        ncalls, times, nonlin_v, results = aux

        # Create and call the DFT solver.
        dft_conf = pb.get_solver_conf(opts.dft_solver)
        dft_status = {}
        dft_solver = Solver.any_from_conf(dft_conf,
                                          fun = nonlin_v,
                                          status = dft_status)
        v_hxc_qp = dft_solver(v_hxc_qp.ravel())

        v_hxc_qp = nm.array(v_hxc_qp, dtype=nm.float64)
        v_hxc_qp.shape = self.qp_shape
        eigs, mtx_s_phi, vec_n, vec_v_h, v_ion_qp, v_xc_qp, v_hxc_qp = results
        output( 'DFT iteration time [s]:', dft_status['time_stats'] )

        fun = pb.materials['mat_v'].function
        variable = self.problem.create_variables(['scalar'])['scalar']
        vec_v_ion = fun(None, variable.field.get_coor(),
                        mode='qp')['V_ion'].squeeze()

        vec_v_xc = self._interp_to_nodes(v_xc_qp)
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)
        vec_v_sum = self._interp_to_nodes(v_hxc_qp + v_ion_qp)

        coor = pb.domain.get_mesh_coors()
        r2 = norm_l2_along_axis(coor, squared=True)
        vec_nr2 = vec_n * r2

        pb.select_bcs( ebc_names = ['ZeroSurface'] )
        mtx_phi = self.make_full( mtx_s_phi )

        out = {}
        update_state_to_output(out, pb, vec_n, 'n')
        update_state_to_output(out, pb, vec_nr2, 'nr2')
        update_state_to_output(out, pb, vec_v_h, 'V_h')
        update_state_to_output(out, pb, vec_v_xc, 'V_xc')
        update_state_to_output(out, pb, vec_v_ion, 'V_ion')
        update_state_to_output(out, pb, vec_v_hxc, 'V_hxc')
        update_state_to_output(out, pb, vec_v_sum, 'V_sum')
        self.save_results(eigs, mtx_phi, out=out)

        if self.options.plot:
            log( save_figure = opts.iter_fig_name )
            pause()
            log(finished=True)

        return Struct( pb = pb, eigs = eigs, mtx_phi = mtx_phi,
                       vec_n = vec_n, vec_nr2 = vec_nr2,
                       vec_v_h = vec_v_h, vec_v_xc = vec_v_xc )

    def test_semismooth_newton(self):
        import numpy as nm
        from sfepy.solvers import Solver

        ns = [0, 6, 2, 2]

        offsets = nm.cumsum(ns)
        nx = offsets[-1]

        iw = slice(offsets[0], offsets[1])
        ig = slice(offsets[1], offsets[2])
        il = slice(offsets[2], offsets[3])

        def fun_smooth(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]
            xl = vec_x[il]

            m = self.ms
            rw = m['A'] * xw - m['Bb'].T * xg - self.m['fs'] 
            rg = m['Bb'] * xw + xl * xg

            rwg = nm.r_[rw, rg]

            return rwg

        def fun_smooth_grad(vec_x):
            xw = vec_x[iw]
            xl = vec_x[il]
            xg = vec_x[ig]

            m = self.m

            mzl = nm.zeros((6, 2), dtype=nm.float64)

            mw = nm.c_[m['A'], -m['Bb'].T, mzl]
            mg = nm.c_[m['Bb'], nm.diag(xl), nm.diag(xg)]

            mx = nm.r_[mw, mg]

            mx = sps.csr_matrix(mx)
            return mx

        def fun_a(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]

            subsd = {}
            for ii, key in enumerate(self.w_names):
                subsd[key] = xw[ii]

            sn = eval_matrix(self.m['sn'], **subsd).squeeze()

            ra = nm.abs(xg) - fc * nm.abs(sn)

            return -ra

        def fun_a_grad(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]
            xl = vec_x[il]

            subsd = {}
            for ii, key in enumerate(self.w_names):
                subsd[key] = xw[ii]

            md = eval_matrix(self.m['D'], **subsd)
            sn = eval_matrix(self.m['sn'], **subsd).squeeze()

            ma = nm.zeros((xl.shape[0], nx), dtype=nm.float64)
            ma[:,iw] = - fc * nm.sign(sn)[:,None] * md
            ma[:,ig] = nm.sign(xg)[:,None] * self.m['C']

            return -sps.csr_matrix(ma)

        def fun_b(vec_x):
            xl = vec_x[il]

            return xl

        def fun_b_grad(vec_x):
            xl = vec_x[il]

            mb = nm.zeros((xl.shape[0], nx), dtype=nm.float64)
            mb[:,il] = self.m['C']

            return sps.csr_matrix(mb)

        vec_x0 = 0.1 * nm.ones((nx,), dtype=nm.float64)

        lin_solver = Solver.any_from_conf(dict_to_struct(ls_conf))
        status = {}
        solver = Solver.any_from_conf(dict_to_struct(conf),
                                      fun_smooth=fun_smooth,
                                      fun_smooth_grad=fun_smooth_grad,
                                      fun_a=fun_a,
                                      fun_a_grad=fun_a_grad,
                                      fun_b=fun_b,
                                      fun_b_grad=fun_b_grad,
                                      lin_solver=lin_solver,
#.........这里部分代码省略.........

def setup_precond(mtx, problem):
    """
    Setup a preconditioner for `mtx`.
    """
    import scipy.sparse.linalg as spla
    from sfepy.solvers import Solver

    # Get active DOF indices for u, p.
    adi = problem.get_variables().adi
    iu = adi.indx['u']
    ip = adi.indx['p']

    # Get the diagonal blocks of the linear system matrix.
    K = mtx[iu, iu]
    M = mtx[ip, ip]

    # Create solvers for K, M blocks to be used in matvec_bj(). A different
    # solver for each block could be used.
    conf = problem.solver_confs['direct']
    # conf = problem.solver_confs['cg-s']
    # conf = problem.solver_confs['cg-p']
    # conf = problem.solver_confs['pyamg']
    ls1 = Solver.any_from_conf(conf, mtx=K, context=problem)
    ls2 = Solver.any_from_conf(conf, mtx=M, context=problem)

    def matvec_bj(vec):
        """
        The application of the Block Jacobi preconditioner.

        The exact version (as with the `'direct'` solver) can be obtained also
        by using the following PETSs command-line options, together with the
        `'iterative-p'` solver::

          -ksp_monitor -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_u_ksp_type preonly -fieldsplit_u_pc_type lu -fieldsplit_p_ksp_type preonly -fieldsplit_p_pc_type lu

        The inexact version (20 iterations of a CG solver for each block, as
        with the `'cg-s'` or `'cg-p'` solvers) can be obtained also by using
        the following PETSs command-line options, together with the
        `'iterative-p'` solver::

          -ksp_monitor -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_u_ksp_type cg -fieldsplit_u_pc_type none -fieldsplit_p_ksp_type cg -fieldsplit_p_pc_type none -fieldsplit_u_ksp_max_it 20 -fieldsplit_p_ksp_max_it 20
        """
        vu = ls1(vec[iu])
        vp = ls2(vec[ip])

        return nm.r_[vu, vp]

    def matvec_j(vec):
        """
        The application of the Jacobi (diagonal) preconditioner.

        The same effect can be obtained also by using the following PETSs
        command-line options, together with the `'iterative-p'` solver::

          -ksp_monitor -pc_type jacobi
        """
        D = mtx.diagonal()

        return vec / D

    # Create the preconditioner, using one of matvec_bj() or matvec_j().
    precond = Struct(name='precond', shape=mtx.shape, matvec=matvec_bj)
    precond = spla.aslinearoperator(precond)

    return precond

def solve_navier_stokes(conf, options):
    opts = conf.options

    dpb = ProblemDefinition.from_conf(conf, init_equations=False)
    equations = getattr(conf, "_".join(("equations_direct", opts.problem)))
    dpb.set_equations(equations)

    ls_conf = dpb.get_solver_conf(opts.ls)
    nls_conf = dpb.get_solver_conf(opts.nls_direct)

    method = opts.direct_method
    if method == "stationary":
        data = {}
        dpb.time_update(None)
        state_dp = dpb.solve(nls_conf=nls_conf)

    elif method == "transient":
        ls = Solver.any_from_conf(ls_conf)
        ts_conf = dpb.get_solver_conf(opts.ts_direct)

        data = {"ts": Struct(dt=ts_conf.dt)}

        # Plug in mass term.
        mequations = {}
        for key, eq in equations.iteritems():
            if "dw_div_grad" in eq:
                eq = "+".join((ts_conf.mass_term, eq)).replace("++", "+")
            mequations[key] = eq

        if ts_conf.stokes_init:
            state_dp0 = solve_stokes(dpb, conf.equations_direct_stokes, nls_conf)
            dpb.set_equations(mequations)
        else:
            dpb.set_equations(mequations)
            state_dp0 = dpb.create_state()
            dpb.time_update(None)
            state_dp0.apply_ebc()

        from sfepy.base.log import Log

        log = Log.from_conf(Struct(is_plot=True), ([r"$||u||$"], [r"$||p||$"]))

        output("Navier-Stokes...")
        ev = BasicEvaluator(dpb, ts=Struct(dt=ts_conf.dt))
        nls = Solver.any_from_conf(nls_conf, evaluator=ev, lin_solver=ls)

        n_step = ts_conf.n_step
        step = 0
        while 1:
            for ii in xrange(n_step):
                output(step)

                vec_u = state_dp0("w")
                vec_p = state_dp0("r")
                log(nm.linalg.norm(vec_u), nm.linalg.norm(vec_p))

                dpb.variables.set_data_from_state("w_0", state_dp0(), "w")
                vec_dp = nls(state_dp0())

                step += 1
                state_dp = state_dp0.copy()
                state_dp.set_reduced(vec_dp)

                state_dp0 = state_dp

            if ts_conf.interactive:
                try:
                    n_step = int(raw_input("continue: "))
                    if n_step <= 0:
                        break
                except:
                    break

        vec_u = state_dp("w")
        vec_p = state_dp("r")
        log(nm.linalg.norm(vec_u), nm.linalg.norm(vec_p), finished=True)

    else:
        raise "unknown Navier-Stokes solution method (%s)!" % method

    return dpb, state_dp, data