diff --git a/python/tamaas/dumpers/numpy_dumper.py b/python/tamaas/dumpers/numpy_dumper.py index 218927a..0d16223 100644 --- a/python/tamaas/dumpers/numpy_dumper.py +++ b/python/tamaas/dumpers/numpy_dumper.py @@ -1,54 +1,56 @@ -# @author Lucas Frérot -# -# @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 . +""" +@author Lucas Frérot + +@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 . +""" from ._helper import FieldDumper import os import numpy as np __all__ = ['NumpyDumper'] class NumpyDumper(FieldDumper): """Dumper to compressed numpy files""" def __init__(self, basename, *fields): default_fields = ['traction', 'displacement', 'stress', 'plastic_strain'] super().__init__(basename, default_fields, *fields) if not os.path.exists('numpys'): os.mkdir('numpys') self.count = 0 def dump(self, model): dump_fields = { field: self.make_periodic[field](model.getField(field)) for field in self.fields } np.savez_compressed('numpys/{}_{:04d}.npz'.format(self.basename, self.count), **dump_fields) self.count += 1 diff --git a/python/tamaas/dumpers/uvw_dumper.py b/python/tamaas/dumpers/uvw_dumper.py index d6a1844..646c51f 100644 --- a/python/tamaas/dumpers/uvw_dumper.py +++ b/python/tamaas/dumpers/uvw_dumper.py @@ -1,71 +1,73 @@ # @author Lucas Frérot # # @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 . from ._helper import FieldDumper import os import uvw import numpy as np __all__ = ["UVWDumper"] class UVWDumper(FieldDumper): """Dumper for elasto-plastic calculations""" def __init__(self, basename, *fields): default_fields = ['displacement', 'stress', 'plastic_strain'] super().__init__(basename, default_fields, *fields) if not os.path.exists('paraview'): os.mkdir('paraview') self.count = 0 def dump(self, model): """Dump displacements, plastic deformations and stresses""" discretization = model.getDiscretization().copy() discretization[1] += 1 discretization[2] += 1 coordinates = [np.linspace(0, L, N) for L, N in zip(model.getSystemSize(), discretization)] # Correct order of coordinate dimensions dimension_indices = [1, 2, 0] # Creating rectilinear grid with correct order for components grid = uvw.RectilinearGrid( "paraview/{}_{:04d}.vtr".format(self.basename, self.count), (coordinates[i] for i in dimension_indices)) dump_fields = { field: uvw.DataArray(self.make_periodic[field](model.getField(field)), dimension_indices, field) for field in self.fields } for array in dump_fields.values(): grid.addPointData(array) grid.write() + + self.count += 1 diff --git a/python/tamaas/nonlinear_solvers/scipy_solver.py b/python/tamaas/nonlinear_solvers/scipy_solver.py index 9cb38ba..c9882f6 100644 --- a/python/tamaas/nonlinear_solvers/scipy_solver.py +++ b/python/tamaas/nonlinear_solvers/scipy_solver.py @@ -1,96 +1,132 @@ """ @author Lucas Frérot @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 . """ from .. import EPSolver import numpy as np from scipy.optimize import newton_krylov, root -__all__ = ["DFSANESolver", "NewtonKrylovSolver"] +__all__ = ["DFSANESolver", "NewtonKrylovSolver", "ToleranceManager"] 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): + """ + Solve using a spectral residual jacobianless method + """ def __init__(self, residual, model): ScipySolver.__init__(self, residual, model) - self.options = {'ftol': 1e-4, 'fatol': 1e-15} + self.options = {'ftol': 1e-9, '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'])) + print("DF-SANE: {} ({} iterations, {:.3e})".format( + solution.message, + solution.nit, + self.options['ftol'])) return solution.x.copy() + + +class ToleranceManager: + """ + Manage tolerance of non-linear solver + """ + def __init__(self, start, end, rate): + self.start = start / rate # just anticipating first multiplication + self.end = end + self.rate = rate + + 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['ftol'] = self.start + + def scipy_solve(obj, *args, **kwargs): + ftol = obj.options['ftol'] + ftol *= self.rate + + obj.options['ftol'] = max(ftol, self.end) + return orig_solve(obj, *args, **kwargs) + + def updateState(obj, *args, **kwargs): + obj.options['ftol'] = self.start + return orig_updateState(obj, *args, **kwargs) + + cls.__init__ = __init__ + cls.scipy_solve = scipy_solve + cls.updateState = updateState + return cls diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp index ab6f160..b5813fb 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,129 +1,133 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 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 . * */ /* -------------------------------------------------------------------------- */ #include "wrap.hh" #include "contact_solver.hh" #include "polonsky_keer_rey.hh" #include "kato.hh" #include "beck_teboulle.hh" #include "condat.hh" #include "polonsky_keer_tan.hh" #include "epic.hh" #include "ep_solver.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace wrap { /* -------------------------------------------------------------------------- */ using namespace py::literals; class PyEPSolver : public EPSolver { public: using EPSolver::EPSolver; void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); } + + void updateState() override { + PYBIND11_OVERLOAD(void, EPSolver, updateState); + } }; /* -------------------------------------------------------------------------- */ void wrapSolvers(py::module& mod) { py::class_(mod, "ContactSolver") // .def(py::init&, Real>(), "model"_a, // "surface"_a, "tolerance"_a) .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a) .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a) .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm) .def("solve", py::overload_cast>(&ContactSolver::solve), "target_force"_a) .def("solve", py::overload_cast(&ContactSolver::solve), "target_normal_pressure"_a); py::class_ pkr(mod, "PolonskyKeerRey"); pkr.def(py::init&, Real, PolonskyKeerRey::type, PolonskyKeerRey::type>(), "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a, "constraint_type"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("computeError", &PolonskyKeerRey::computeError); py::enum_(pkr, "type") .value("gap", PolonskyKeerRey::gap) .value("pressure", PolonskyKeerRey::pressure) .export_values(); py::class_ kato(mod, "Kato"); kato.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a=50) .def("solveRelaxed", &Kato::solveRelaxed, "g0"_a) .def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a=0.01) .def("computeCost", &Kato::computeCost, "use_tresca"_a=false); py::class_ bt(mod, "BeckTeboulle"); bt.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &BeckTeboulle::solve, "p0"_a) .def("computeCost", &BeckTeboulle::computeCost); py::class_ cd(mod, "Condat"); cd.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &Condat::solve, "p0"_a, "grad_step"_a=0.9) .def("computeCost", &Condat::computeCost); py::class_ pkt(mod, "PolonskyKeerTan"); pkt.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &PolonskyKeerTan::solve, "p0"_a) .def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a) .def("computeCost", &PolonskyKeerTan::computeCost, "use_tresca"_a=false); py::class_(mod, "EPSolver") .def(py::init()) .def("solve", &EPSolver::solve) .def("getStrainIncrement", &EPSolver::getStrainIncrement) .def("getResidual", &EPSolver::getResidual) .def("updateState", &EPSolver::updateState); py::class_(mod, "EPICSolver") .def(py::init(), "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10) .def("solve", &EPICSolver::solve, "force"_a) .def("solve", [](EPICSolver& solver, Real pressure) { return solver.solve(std::vector{pressure}); }, "normal_pressure"_a); } } // namespace wrap __END_TAMAAS__ diff --git a/src/solvers/ep_solver.hh b/src/solvers/ep_solver.hh index 2eea028..03c2def 100644 --- a/src/solvers/ep_solver.hh +++ b/src/solvers/ep_solver.hh @@ -1,62 +1,62 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2016-2017 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef EP_SOLVER_HH #define EP_SOLVER_HH /* -------------------------------------------------------------------------- */ #include "grid_base.hh" #include "residual.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ class EPSolver { public: /// Constructor EPSolver(Residual& residual); /// Destructor virtual ~EPSolver() = default; // Pure virtual methods public: virtual void solve() = 0; public: - void updateState(); + virtual void updateState(); // Accessors public: GridBase& getStrainIncrement() { return *_x; } Residual& getResidual() { return _residual; } protected: std::unique_ptr> _x; Residual& _residual; }; /* -------------------------------------------------------------------------- */ } // namespace tamaas #endif