diff --git a/python/tamaas/nonlinear_solvers/scipy_solver.py b/python/tamaas/nonlinear_solvers/scipy_solver.py index e706961..9fcc523 100644 --- a/python/tamaas/nonlinear_solvers/scipy_solver.py +++ b/python/tamaas/nonlinear_solvers/scipy_solver.py @@ -1,152 +1,152 @@ # @file # @section LICENSE # # Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . """ Solvers using scipy.optimize routines """ from .. import EPSolver from scipy.optimize import newton_krylov, root __all__ = ["DFSANESolver", "NewtonKrylovSolver", "ToleranceManager"] class NLNoConvergence(Exception): """Convergence not reached exception""" pass class ScipySolver(EPSolver): """ Base class for solvers wrapping SciPy routines """ def __init__(self, residual, model, callback=None): super().__init__(residual) self.model = model self.callback = callback self._x = self.getStrainIncrement() self._residual = self.getResidual() self.options = {'ftol': 0, 'fatol': 1e-9} @property def tolerance(self): return self.options['fatol'] @tolerance.setter def set_tolerance(self, v): self.options['fatol'] = v def solve(self): """ Solve the nonlinear plasticity equation using the scipy_solve routine """ # For initial guess, compute the strain due to boundary tractions # self._residual.computeResidual(self._x) # self._x[...] = self._residual.getVector() # Scipy root callback def compute_residual(x): self._residual.computeResidual(x) return self._residual.getVector().copy() # Solve self._x[...] = self.scipy_solve(compute_residual) # Computing displacements self._residual.computeResidualDisplacement(self._x) def reset(self): self._x[...] = 0 class NewtonKrylovSolver(ScipySolver): """ Solve using a finite-difference Newton-Krylov method """ def __init__(self, residual, model, callback=None): ScipySolver.__init__(self, residual, model, callback) def scipy_solve(self, compute_residual): try: return newton_krylov(compute_residual, self._x, f_tol=self.options['fatol'], f_rtol=self.options['ftol'], verbose=True, callback=self.callback) except Exception as e: raise NLNoConvergence(e.what()) class DFSANESolver(ScipySolver): """ Solve using a spectral residual jacobianless method """ def __init__(self, residual, model, callback=None): ScipySolver.__init__(self, residual, model, callback) def scipy_solve(self, compute_residual): solution = root(compute_residual, self._x, method='df-sane', options=self.options, callback=self.callback) print("DF-SANE: {} ({} iterations, {})".format( solution.message, solution.nit, self.options)) - if not solution.sucess: + if not solution.success: raise NLNoConvergence("DF-SANE did not converge") return solution.x.copy() class ToleranceManager: """ Decorator to manage tolerance of non-linear solver """ def __init__(self, start, end, rate, key='fatol'): self.start = start / rate # just anticipating first multiplication self.end = end self.rate = rate self.key = key def __call__(self, cls): orig_init = cls.__init__ orig_solve = cls.scipy_solve orig_updateState = cls.updateState def __init__(obj, *args, **kwargs): orig_init(obj, *args, **kwargs) obj.options[self.key] = self.start def scipy_solve(obj, *args, **kwargs): ftol = obj.options[self.key] ftol *= self.rate obj.options[self.key] = max(ftol, self.end) return orig_solve(obj, *args, **kwargs) def updateState(obj, *args, **kwargs): obj.options[self.key] = self.start return orig_updateState(obj, *args, **kwargs) cls.__init__ = __init__ cls.scipy_solve = scipy_solve cls.updateState = updateState return cls