Page MenuHomec4science
Files F65263899

scipy_solver.py
No OneTemporary
Actions

File Metadata

Created
Sun, Jun 2, 09:13

scipy_solver.py
View Options

"""
@author Lucas Frérot <lucas.frerot@epfl.ch>
@section LICENSE
Copyright (©) 2018 EPFL (Ecole Polytechnique Fédérale de
Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
Solides)
Tamaas is free software: you can redistribute it and/or modify it under the
terms of the GNU Lesser General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option) any
later version.
Tamaas 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 Lesser General Public License for more
details.
You should have received a copy of the GNU Lesser General Public License
along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
"""
from .. import EPSolver
import numpy as np
from scipy.optimize import newton_krylov, root
__all__ = ["DFSANESolver", "NewtonKrylovSolver"]
class ScipySolver(EPSolver):
"""
Base class for solvers wrapping SciPy routines
"""
def __init__(self, residual, model):
super().__init__(residual)
self.model = model
self._x = self.getStrainIncrement()
self._residual = self.getResidual()
self.surface_correction = model.getDisplacement()[0, :, :, 2]
self.last_residual_disp = np.zeros_like(self.surface_correction)
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.last_residual_disp[...] = 0
self._x[...] = 0
class NewtonKrylovSolver(ScipySolver):
"""
Solve using a finite-difference Newton-Krylov method
"""
def __init__(self, residual, model):
ScipySolver.__init__(self, residual, model)
self.options = {'ftol': 1e-8, 'fatol': 1e-15}
def scipy_solve(self, compute_residual):
return newton_krylov(compute_residual, self._x,
f_tol=1e-8, verbose=True)
class DFSANESolver(ScipySolver):
def __init__(self, residual, model):
ScipySolver.__init__(self, residual, model)
self.options = {'ftol': 1e-4, 'fatol': 1e-15}
def scipy_solve(self, compute_residual):
solution = root(compute_residual, self._x, method='df-sane',
options=self.options)
self.options['ftol'] /= 2
if self.options['ftol'] < 1e-8:
self.options['ftol'] = 1e-8
print("DF-SANE: {} ({} iterations, {})".format(solution.message,
solution.nit,
self.options['ftol']))
return solution.x.copy()

Event Timeline

c4science · Help