diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp index 6cbd075..39eb938 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,246 +1,252 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 Lucas Frérot * * 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 . * */ /* -------------------------------------------------------------------------- */ +#include "anderson.hh" #include "beck_teboulle.hh" #include "condat.hh" #include "contact_solver.hh" #include "dfsane_solver.hh" #include "ep_solver.hh" #include "epic.hh" #include "kato.hh" #include "kato_saturated.hh" #include "polonsky_keer_rey.hh" #include "polonsky_keer_tan.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { /* -------------------------------------------------------------------------- */ using namespace py::literals; class PyEPSolver : public EPSolver { public: using EPSolver::EPSolver; // NOLINTNEXTLINE(readability-else-after-return) void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); } void updateState() override { // NOLINTNEXTLINE(readability-else-after-return) PYBIND11_OVERLOAD(void, EPSolver, updateState); } }; class PyContactSolver : public ContactSolver { public: using ContactSolver::ContactSolver; Real solve(std::vector load) override { PYBIND11_OVERLOAD(Real, ContactSolver, solve, load); } Real solve(Real load) override { PYBIND11_OVERLOAD(Real, ContactSolver, solve, load); } }; /* -------------------------------------------------------------------------- */ void wrapSolvers(py::module& mod) { py::class_(mod, "ContactSolver") .def(py::init&, Real>(), py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def( "setMaxIterations", [](ContactSolver& m, UInt iter) { TAMAAS_DEPRECATE("setMaxIterations()", "the max_iter property"); m.setMaxIterations(iter); }, "max_iter"_a) .def( "setDumpFrequency", [](ContactSolver& m, UInt freq) { TAMAAS_DEPRECATE("setDumpFrequency()", "the dump_freq property"); m.setDumpFrequency(freq); }, "dump_freq"_a) .def_property("tolerance", &ContactSolver::getTolerance, &ContactSolver::setTolerance, "Solver tolerance") .def_property("max_iter", &ContactSolver::getMaxIterations, &ContactSolver::setMaxIterations, "Maximum number of iterations") .def_property("dump_freq", &ContactSolver::getDumpFrequency, &ContactSolver::setDumpFrequency, "Frequency of displaying solver info") .def_property_readonly("model", &ContactSolver::getModel) .def_property_readonly("surface", &ContactSolver::getSurface) .def_property("functional", &ContactSolver::getFunctional, &ContactSolver::setFunctional) .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm, "Add a term to the contact functional to minimize") .def("solve", py::overload_cast>(&ContactSolver::solve), "target_force"_a, py::call_guard(), "Solve the contact for a mean traction/gap vector") .def("solve", py::overload_cast(&ContactSolver::solve), "target_normal_pressure"_a, py::call_guard(), "Solve the contact for a mean normal pressure/gap"); py::class_ pkr( mod, "PolonskyKeerRey", "Main solver class for normal elastic contact problems. Its functional " "can be customized to add an adhesion term, and its primal variable can " "be set to either the gap or the pressure."); // Need to export enum values before defining PKR constructor py::enum_(pkr, "type") .value("gap", PolonskyKeerRey::gap) .value("pressure", PolonskyKeerRey::pressure) .export_values(); pkr.def(py::init&, Real, PolonskyKeerRey::type, PolonskyKeerRey::type>(), "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a = PolonskyKeerRey::type::pressure, "constraint_type"_a = PolonskyKeerRey::type::pressure, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("computeError", &PolonskyKeerRey::computeError) .def("setIntegralOperator", &PolonskyKeerRey::setIntegralOperator); py::class_( mod, "KatoSaturated", "Solver for pseudo-plasticity problems where the normal pressure is " "constrained above by a saturation pressure \"pmax\"") .def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "pmax"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def_property("pmax", &KatoSaturated::getPMax, &KatoSaturated::setPMax, "Saturation normal pressure"); 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", "Main solver for frictional contact problems. It has no restraint on the " "material properties or friction coefficient values, but solves an " "associated version of the Coulomb friction law, which differs from the " "traditional Coulomb friction in that the normal and tangential slip " "components are coupled."); 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, "_tolerance_manager", "Manager object for the tolereance of nonlinear plasticity solvers. " "Decreases the solver tolerance by geometric progression.") .def(py::init(), "start_tol"_a, "end_tol"_a, "rate"_a); py::class_( mod, "EPSolver", "Mother class for nonlinear plasticity solvers") .def(py::init(), "residual"_a, py::keep_alive<1, 2>()) .def("solve", &EPSolver::solve) .def("getStrainIncrement", &EPSolver::getStrainIncrement, py::return_value_policy::reference_internal) .def("getResidual", &EPSolver::getResidual, py::return_value_policy::reference_internal) .def("updateState", &EPSolver::updateState) .def_property("tolerance", &EPSolver::getTolerance, &EPSolver::setTolerance) .def("setToleranceManager", &EPSolver::setToleranceManager) .def("beforeSolve", &EPSolver::beforeSolve); py::class_(mod, "_DFSANESolver") .def(py::init(), "residual"_a, py::keep_alive<1, 2>()) .def(py::init([](Residual& res, Model&) { TAMAAS_DEPRECATE("Solver(residual, model)", "Solver(residual)"); return std::make_unique(res); }), "residual"_a, "model"_a, py::keep_alive<1, 2>()); py::class_( mod, "EPICSolver", "Main solver class for elastic-plastic contact problems") .def(py::init(), "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10, "relaxation"_a = 0.3, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def( "solve", [](EPICSolver& solver, Real pressure) { return solver.solve(std::vector{pressure}); }, "normal_pressure"_a, py::call_guard(), "Solves the EP contact with a relaxed fixed-point scheme. Adjust " "the relaxation parameter to help convergence.") .def( "acceleratedSolve", [](EPICSolver& solver, Real pressure) { return solver.acceleratedSolve(std::vector{pressure}); }, "normal_pressure"_a, py::call_guard(), "Solves the EP contact with an accelerated fixed-point scheme. May " "not converge!") .def_property("tolerance", &EPICSolver::getTolerance, &EPICSolver::setTolerance) .def_property("relaxation", &EPICSolver::getRelaxation, &EPICSolver::setRelaxation) .def_property_readonly("model", &EPICSolver::getModel); + + py::class_(mod, "AndersonMixing") + .def(py::init(), + "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10, + "memory"_a = 5, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()); } } // namespace wrap } // namespace tamaas diff --git a/src/SConscript b/src/SConscript index 2ed290c..fdd85ae 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,147 +1,148 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . from SCons.Script import Import, Export def prepend(env, path, list): return [env.File(x, path) for x in list] Import('main_env') env = main_env.Clone() env.AddMethod(prepend, 'PrependDir') # Core core_list = """ fft_engine.cpp grid.cpp grid_hermitian.cpp statistics.cpp tamaas.cpp loop.cpp computes.cpp logger.cpp mpi_interface.cpp """.split() core_list = env.PrependDir('core', core_list) if env['use_fftw']: core_list += ['core/fftw/fftw_engine.cpp'] if env['use_mpi']: core_list += ['core/fftw/mpi/fftw_mpi_engine.cpp'] if env['use_cuda']: core_list += ['core/cuda/cufft_engine.cpp'] info_file = env.Substfile('tamaas_info.cpp', 'tamaas_info.cpp.in') env.Depends(info_file, Value(env.subst(env['SUBST_DICT'].values()))) core_list.append(info_file) # Lib roughcontact generator_list = """ surface_generator.cpp surface_generator_filter.cpp surface_generator_random_phase.cpp isopowerlaw.cpp regularized_powerlaw.cpp """.split() generator_list = env.PrependDir('surface', generator_list) # Lib PERCOLATION percolation_list = """ flood_fill.cpp """.split() percolation_list = env.PrependDir('percolation', percolation_list) # Model model_list = """ model.cpp model_factory.cpp model_type.cpp model_template.cpp integral_operator.cpp be_engine.cpp westergaard.cpp dcfft.cpp elastic_functional.cpp meta_functional.cpp adhesion_functional.cpp volume_potential.cpp kelvin.cpp mindlin.cpp boussinesq.cpp hooke.cpp residual.cpp materials/isotropic_hardening.cpp materials/mfront_material.cpp integration/element.cpp """.split() model_list = env.PrependDir('model', model_list) # Solvers solvers_list = """ contact_solver.cpp polonsky_keer_rey.cpp kato_saturated.cpp kato.cpp beck_teboulle.cpp condat.cpp polonsky_keer_tan.cpp ep_solver.cpp dfsane_solver.cpp epic.cpp +anderson.cpp """.split() solvers_list = env.PrependDir('solvers', solvers_list) # Assembling total list rough_contact_list = \ core_list + generator_list + percolation_list + model_list + solvers_list # Adding extra warnings for Tamaas base lib env.AppendUnique(CXXFLAGS=['-Wextra']) # Allowing libTamaas.so to find libs in the same directory env.AppendUnique(RPATH=["'$$$$ORIGIN'"]) # Build static library for packaging if env['build_static_lib']: env.AppendUnique(CXXFLAGS='-fPIC') libTamaas = env.StaticLibrary('Tamaas', rough_contact_list) # Build shared library (default) else: libTamaas = env.SharedLibrary('Tamaas', rough_contact_list, SHLIBVERSION=env['version']) # Specify install target to install lib lib_prefix = env.Dir('lib', env['prefix']) lib_install = env.InstallVersionedLib(target=lib_prefix, source=libTamaas) # Defining alias targets main_env.Alias('build-cpp', libTamaas) main_env.Alias('install-lib', lib_install) # Export target for use in python builds Export('libTamaas') diff --git a/src/solvers/anderson.cpp b/src/solvers/anderson.cpp new file mode 100644 index 0000000..529d1eb --- /dev/null +++ b/src/solvers/anderson.cpp @@ -0,0 +1,192 @@ +/* + * SPDX-License-Indentifier: AGPL-3.0-or-later + * + * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * Copyright (©) 2020-2023 Lucas Frérot + * + * 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 . + * + */ +/* -------------------------------------------------------------------------- */ +#include "anderson.hh" +#include +#include +/* -------------------------------------------------------------------------- */ +namespace tamaas { + +/* -------------------------------------------------------------------------- */ +AndersonMixing::AndersonMixing(ContactSolver& csolver, EPSolver& epsolver, + Real tolerance, UInt memory) + : EPICSolver(csolver, epsolver, tolerance, 0.5), M(memory) {} + +/* -------------------------------------------------------------------------- */ +Real AndersonMixing::solve(const std::vector& load) { + UInt n = 0, nmax = 1000; + Real error = 0.; + const GridBase initial_surface(surface); + Real normalizing_factor = std::sqrt(initial_surface.var()); + GridBase x(*residual_disp), x_prev(*residual_disp); + memory_t memory, residual_memory; + + do { + { + // Compute a residual vector for this iteration + GridBase F(*residual_disp); + fixedPoint(F, x, initial_surface, load); + F -= x; // F = g(x) - x + memory.push_front(x); + residual_memory.push_front(std::move(F)); + } + + // Keep M previous iterations: memory has M + 1 vectors + if (memory.size() > M + 1) + memory.pop_back(); + if (residual_memory.size() > M + 1) + residual_memory.pop_back(); + + x = mixingUpdate(std::move(x), computeGamma(residual_memory), memory, + residual_memory, relaxation); + error = computeError(x, x_prev, normalizing_factor); + x_prev = x; + + Logger().get(LogLevel::info) << "[Anderson] error = " << std::scientific + << error << std::fixed << std::endl; + } while (error > tolerance && n++ < nmax); + + std::cout << "Converged in " << n << " iterations\n"; + + // computing full pressure + pressure += *pressure_inc; + // setting the model pressure to full converged pressure + *pressure_inc = pressure; + // updating plastic state + epsolver.updateState(); + return error; +} + +/* -------------------------------------------------------------------------- */ +GridBase AndersonMixing::mixingUpdate(GridBase x, + std::vector gamma, + const memory_t& memory, + const memory_t& residual_memory, + Real relaxation) { + Loop::loop([beta = relaxation](Real& x, Real Fl) { x += beta * Fl; }, x, + residual_memory.front()); + + for (UInt m = 1; m < memory.size(); ++m) { + Loop::loop( + [gm = gamma[m - 1], beta = relaxation](Real& x, Real xmp1, Real xm, + Real Fmp1, Real Fm) { + // Eyert (10.1006/jcph.1996.0059) p282 eq (7.7) + x -= gm * ((xmp1 - xm) + beta * (Fmp1 - Fm)); + }, + x, memory[m - 1], memory[m], residual_memory[m - 1], + residual_memory[m]); + } + + return x; +} + +/* -------------------------------------------------------------------------- */ +/// Crout(e) (Antoine... et Daniel...) factorization from +/// https://en.wikipedia.org/wiki/Crout_matrix_decomposition +auto factorize(Grid A) { + TAMAAS_ASSERT(A.sizes()[0] == A.sizes()[1], "Matrix is not square"); + const auto N = A.sizes()[0]; + auto LU = + std::make_pair(Grid{A.sizes(), 1}, Grid{A.sizes(), 1}); + auto& L = LU.first; + auto& U = LU.second; + + for (UInt i = 0; i < N; ++i) { + U(i, i) = 1; + } + + for (UInt j = 0; j < N; ++j) { + for (UInt i = j; i < N; ++i) { + double sum = 0; + for (UInt k = 0; k < j; ++k) { + sum += L(i, k) * U(k, j); + } + + L(i, j) = A(i, j) - sum; + } + + for (UInt i = j; i < N; ++i) { + double sum = 0; + for (UInt k = 0; k < j; ++k) { + sum += L(j, k) * U(k, i); + } + + U(j, i) = (A(j, i) - sum) / L(j, j); + } + } + return LU; +} + +/* -------------------------------------------------------------------------- */ +auto LUsubstitute(std::pair, Grid> LU, Grid b) { + const auto& L = LU.first; + const auto& U = LU.second; + const auto N = b.sizes()[0]; + + // Forward substitution + for (UInt m = 0; m < N; ++m) { + for (UInt i = 0; i < m; ++i) + b(m) -= L(m, i) * b(i); + b(m) /= L(m, m); + } + + // Backward substitution + for (UInt l = 0; l < N; ++l) { + auto m = N - 1 - l; + for (UInt i = m + 1; i < N; ++i) + b(m) -= U(m, i) * b(i); + b(m) /= U(m, m); + } + + return b; +} + +/* -------------------------------------------------------------------------- */ +std::vector AndersonMixing::computeGamma(const memory_t& residual) { + auto N = residual.size() - 1; + Grid A({N, N}, 1); + Grid b({N}, 1); + const auto& Fl = residual.front(); + + // Fill RHS vector + for (UInt n = 1; n < residual.size(); ++n) { + b(n - 1) = residual[n - 1].dot(Fl) - residual[n].dot(Fl); + } + + /// Regularization from Eyert + auto w0 = 0; + // Fill matrix + for (UInt n = 1; n < residual.size(); ++n) { + for (UInt m = 1; m < residual.size(); ++m) { + A(n - 1, m - 1) = residual[n - 1].dot(residual[m - 1]) - + residual[n - 1].dot(residual[m]) - + residual[n].dot(residual[m - 1]) + + residual[n].dot(residual[m]); + A(n - 1, m - 1) *= (1 + w0 * w0 * (n == m)); + } + } + + auto x = LUsubstitute(factorize(A), b); + return std::vector(x.begin(), x.end()); +} + +} // namespace tamaas diff --git a/src/solvers/anderson.hh b/src/solvers/anderson.hh new file mode 100644 index 0000000..87f415d --- /dev/null +++ b/src/solvers/anderson.hh @@ -0,0 +1,53 @@ +/* + * SPDX-License-Indentifier: AGPL-3.0-or-later + * + * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * Copyright (©) 2020-2023 Lucas Frérot + * + * 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 . + * + */ +/* -------------------------------------------------------------------------- */ +#ifndef ANDERSON_HH +#define ANDERSON_HH +/* -------------------------------------------------------------------------- */ +#include "epic.hh" +#include +/* -------------------------------------------------------------------------- */ +namespace tamaas { + +class AndersonMixing : public EPICSolver { + using memory_t = std::deque>; + +public: + AndersonMixing(ContactSolver& csolver, EPSolver& epsolver, + Real tolerance = 1e-10, UInt memory = 5); + + Real solve(const std::vector& load) override; + +private: + static std::vector computeGamma(const memory_t& residual); + static GridBase mixingUpdate(GridBase x, std::vector gamma, + const memory_t& memory, + const memory_t& residual_memory, + Real relaxation); + +private: + UInt M; +}; + +} // namespace tamaas +/* -------------------------------------------------------------------------- */ +#endif diff --git a/src/solvers/epic.hh b/src/solvers/epic.hh index 4b6f9c4..1f13d7e 100644 --- a/src/solvers/epic.hh +++ b/src/solvers/epic.hh @@ -1,88 +1,88 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 Lucas Frérot * * 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef EPIC_HH #define EPIC_HH /* -------------------------------------------------------------------------- */ #include "contact_solver.hh" #include "ep_solver.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ class EPICSolver { public: /// Constructor EPICSolver(ContactSolver& csolver, EPSolver& epsolver, Real tolerance = 1e-10, Real relaxation = 0.3); - Real solve(const std::vector& load); + virtual Real solve(const std::vector& load); Real acceleratedSolve(const std::vector& load); Real computeError(const GridBase& current, const GridBase& prev, Real factor) const; void fixedPoint(GridBase& result, const GridBase& x, const GridBase& initial_surface, std::vector load); template void setViews(); TAMAAS_ACCESSOR(tolerance, Real, Tolerance); TAMAAS_ACCESSOR(relaxation, Real, Relaxation); const Model& getModel() const { return csolver.getModel(); } protected: GridBase surface; ///< corrected surface GridBase pressure; ///< current pressure std::unique_ptr> residual_disp; ///< plastic residual disp std::unique_ptr> pressure_inc; ///< pressure increment ContactSolver& csolver; EPSolver& epsolver; Real tolerance, relaxation; }; /* -------------------------------------------------------------------------- */ /* Template impl. */ /* -------------------------------------------------------------------------- */ template void EPICSolver::setViews() { constexpr UInt dim = model_type_traits::dimension; constexpr UInt bdim = model_type_traits::boundary_dimension; constexpr UInt comp = model_type_traits::components; pressure_inc = std::unique_ptr>{new GridView( csolver.getModel().getTraction(), {}, comp - 1)}; residual_disp = std::unique_ptr>{new GridView( csolver.getModel().getDisplacement(), model_type_traits::indices, comp - 1)}; } /* -------------------------------------------------------------------------- */ } // namespace tamaas #endif // EPIC_HH diff --git a/tests/test_epic.py b/tests/test_epic.py index cb3f87f..ed24b09 100644 --- a/tests/test_epic.py +++ b/tests/test_epic.py @@ -1,78 +1,77 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . from __future__ import division, print_function import pytest from conftest import pkridfn import tamaas as tm from numpy.linalg import norm from tamaas.nonlinear_solvers import DFSANESolver @pytest.fixture(scope="module", - params=[tm.PolonskyKeerRey.pressure], - ids=pkridfn) + params=[tm.EPICSolver, tm.AndersonMixing], + ids=["EPIC", "Anderson"]) def model_setup(hertz_coarse, request): models, solvers = [], [] for mtype in [tm.model_type.basic_2d, tm.model_type.volume_2d]: dim = tm.type_traits[mtype].dimension n = [hertz_coarse.n] * dim d = [hertz_coarse.domain_size] * dim if dim == 3: n[0] = 2 d[0] = 0.001 model = tm.ModelFactory.createModel(mtype, d, n) model.E, model.nu = hertz_coarse.e_star, 0 - solver = tm.PolonskyKeerRey(model, hertz_coarse.surface, 1e-12, - request.param, request.param) + solver = tm.PolonskyKeerRey(model, hertz_coarse.surface, 1e-12) models.append(model) solvers.append(solver) mat = tm.materials.IsotropicHardening(model, sigma_y=1e2 * model.E, hardening=0) epsolver = DFSANESolver(tm.Residual(model, mat)) - solvers[1] = tm.EPICSolver(solvers[1], epsolver, 1e-12) + solvers[1] = request.param(solvers[1], epsolver, 1e-12) for solver in solvers: solver.solve(hertz_coarse.load) return models, hertz_coarse @pytest.mark.xfail(tm.TamaasInfo.backend == 'tbb', reason='TBB reductions are unstable?') def test_energy(model_setup): """Test the full elastic-plastic solve step in elasticity""" # We use computed values from basic PKR to test (model_elastic, model_plastic), hertz = model_setup error = norm( (model_plastic.traction[..., 2] - model_elastic.traction) * (model_plastic.displacement[0, ..., 2] - model_elastic.displacement)) error /= norm(hertz.pressure * hertz.displacement) assert error < 1e-16