diff --git a/python/tamaas/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py index 9d2708b..db969c2 100644 --- a/python/tamaas/nonlinear_solvers/__init__.py +++ b/python/tamaas/nonlinear_solvers/__init__.py @@ -1,176 +1,176 @@ # -*- mode:python; 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 . """Nonlinear solvers for plasticity problems. Solvers in this module use :py:mod:`scipy.optimize` to solve the implicit non-linear equation for plastic deformations with fixed contact pressures. """ from functools import wraps from scipy.optimize import newton_krylov, root -from scipy.optimize.nonlin import NoConvergence +from scipy.optimize._nonlin import NoConvergence from .._tamaas import ( EPSolver, Logger, LogLevel, mpi, _tolerance_manager, _DFSANESolver as DFSANECXXSolver, ) __all__ = ['NLNoConvergence', 'DFSANESolver', 'DFSANECXXSolver', 'NewtonKrylovSolver', 'ToleranceManager'] class NLNoConvergence(Exception): """Convergence not reached exception.""" class ScipySolver(EPSolver): """Base class for solvers wrapping SciPy routines.""" def __init__(self, residual, model=None, callback=None): """Construct nonlinear solver with residual. :param residual: plasticity residual object :param model: Deprecated :param callback: passed on to the Scipy solver """ super(ScipySolver, self).__init__(residual) if mpi.size() > 1: raise RuntimeError("Scipy solvers cannot be used with MPI; " "DFSANECXXSolver can be used instead") self.callback = callback self._x = self.getStrainIncrement() self._residual = self.getResidual() self.options = {'ftol': 0, 'fatol': 1e-9} def solve(self): """Solve R(δε) = 0 using a Scipy function.""" # For initial guess, compute the strain due to boundary tractions # self._residual.computeResidual(self._x) # self._x[...] = self._residual.getVector() EPSolver.beforeSolve(self) # Scipy root callback def compute_residual(vec): self._residual.computeResidual(vec) return self._residual.vector.copy() # Solve Logger().get(LogLevel.debug) << \ "Entering non-linear solve\n" self._x[...] = self._scipy_solve(compute_residual) Logger().get(LogLevel.debug) << \ "Non-linear solve returned" # Computing displacements self._residual.computeResidualDisplacement(self._x) def reset(self): """Set solution vector to zero.""" self._x[...] = 0 def _scipy_solve(self, compute_residual): """Actually call the scipy solver. :param compute_residual: function returning residual for given variable """ raise NotImplementedError() class NewtonKrylovSolver(ScipySolver): """Solve using a finite-difference Newton-Krylov method.""" def _scipy_solve(self, compute_residual): try: return newton_krylov(compute_residual, self._x, f_tol=self.tolerance, verbose=True, callback=self.callback) except NoConvergence: raise NLNoConvergence("Newton-Krylov did not converge") class DFSANESolver(ScipySolver): """Solve using a spectral residual jacobianless method. See :doi:`10.1090/S0025-5718-06-01840-0` for details on method and the relevant Scipy `documentation `_ for details on parameters. """ def _scipy_solve(self, compute_residual): solution = root(compute_residual, self._x, method='df-sane', options={'ftol': 0, 'fatol': self.tolerance}, callback=self.callback) Logger().get(LogLevel.info) << \ "DF-SANE/Scipy: {} ({} iterations, {})".format( solution.message, solution.nit, self.tolerance) if not solution.success: raise NLNoConvergence("DF-SANE/Scipy did not converge") return solution.x.copy() def ToleranceManager(start, end, rate): """Decorate solver to manage tolerance of non-linear solve step.""" def actual_decorator(cls): orig_init = cls.__init__ orig_solve = cls.solve orig_update_state = cls.updateState @wraps(cls.__init__) def __init__(obj, *args, **kwargs): orig_init(obj, *args, **kwargs) obj.setToleranceManager(_tolerance_manager(start, end, rate)) @wraps(cls.solve) def new_solve(obj, *args, **kwargs): ftol = obj.tolerance ftol *= rate obj.tolerance = max(ftol, end) return orig_solve(obj, *args, **kwargs) @wraps(cls.updateState) def updateState(obj, *args, **kwargs): obj.tolerance = start return orig_update_state(obj, *args, **kwargs) cls.__init__ = __init__ # cls.solve = new_solve # cls.updateState = updateState return cls return actual_decorator diff --git a/tests/pytest.ini b/tests/pytest.ini index b0e5a94..1a5c089 100644 --- a/tests/pytest.ini +++ b/tests/pytest.ini @@ -1,3 +1 @@ -[pytest] -filterwarnings = - ignore::DeprecationWarning \ No newline at end of file +[pytest] \ No newline at end of file diff --git a/tests/test_flood_fill.py b/tests/test_flood_fill.py index f0fc859..c3ebbf4 100644 --- a/tests/test_flood_fill.py +++ b/tests/test_flood_fill.py @@ -1,104 +1,104 @@ # -*- 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 print_function, division import tamaas as tm import numpy as np def test_flood_fill2d(): - field = np.zeros((18, 18), np.bool) + field = np.zeros((18, 18), bool) # Cluster of permeter 18 area 13 field[2:4, 3:6] = True field[4, 4:6] = True # noqa field[5, 5] = True # noqa field[3:5, 6:8] = True # noqa # Cluster of perimeter 8 area 4 field[14:16, 3:5] = True # Cluster of perimeter 20 area 11 field[9:11, 8:11] = True # noqa field[7:9, 9] = True # noqa field[10, 11:14] = True # noqa # Cluster of perimeter 18 area 9 field[3:5, 14:16] = True field[5:10, 15] = True # Cluster spanning periodic boundary field[0:3, 12:14] = True field[-2:, 12:14] = True # Percolating cluster field[12, :] = True if False: # Show contact map import matplotlib.pyplot as plt plt.imshow(field, interpolation='none', cmap='gray_r', origin='lower') plt.xlim(0, 17) plt.ylim(0, 17) plt.show() clusters = tm.FloodFill.getClusters(field, False) # Dummy class for clusters class dummy: def __init__(self, area, perimeter, extent=None): self.area = area self.perimeter = perimeter self.extent = extent # Solution ref = [ dummy(10, 14, [5, 2]), # cluster spanning pbc dummy(13, 18), dummy(9, 18), dummy(11, 20), dummy(18, 36, [1, 18]), # percolating cluster dummy(4, 8) ] assert len(ref) == len(clusters) for x, y in zip(clusters, ref): assert x.area == y.area assert x.perimeter == y.perimeter if y.extent is not None: assert x.extent == y.extent field[:] = True field[2:4, 2:4] = False clusters = tm.FloodFill.getClusters(field, False) assert clusters[0].extent == [18, 18] def test_flood_fill3d(): - field = np.zeros((5, 5, 5), np.bool) + field = np.zeros((5, 5, 5), bool) field[2:4, 3, 1:3] = True clusters = tm.FloodFill.getVolumes(field, False) assert len(clusters) == 1 assert clusters[0].area == 4 diff --git a/tests/test_hertz.py b/tests/test_hertz.py index d6c152f..404861e 100644 --- a/tests/test_hertz.py +++ b/tests/test_hertz.py @@ -1,72 +1,72 @@ # -*- 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 numpy as np from numpy.linalg import norm import tamaas as tm def compute_perimeter(tractions): - contact = np.zeros_like(tractions, np.bool) + contact = np.zeros_like(tractions, bool) contact[tractions > 0] = True # Computing clusters clusters = tm.FloodFill.getClusters(contact, False) return np.sum([c.perimeter for c in clusters]) def test_contact_area(pkr): "Testing contact area against Hertz solution for solids of revolution" model, hertz = pkr def error(erzt, true): return np.abs(erzt - true) / true tractions = model['traction'] true_area = hertz.a**2 * np.pi contact_area = tm.Statistics2D.contact(tractions) hertz_area = np.mean(hertz.pressure > 0) assert error(contact_area, hertz_area) < 1e-2 assert error(contact_area, true_area) < 5e-3 contact_area = tm.Statistics2D.contact(tractions, compute_perimeter(tractions)) hertz_perimeter = compute_perimeter(hertz.pressure) hertz_area -= (np.pi - 1 + np.log(2)) / 24 \ * hertz_perimeter / tractions.size assert error(contact_area, hertz_area) < 1e-2 assert error(contact_area, true_area) < 3e-3 def test_pressure(pkr): model, hertz = pkr error = norm(model['traction'] - hertz.pressure) / norm(hertz.pressure) assert error < 5e-3 def test_displacement(pkr): model, hertz = pkr error = norm(model['displacement'] - hertz.displacement)\ / norm(hertz.displacement) assert error < 8e-3 diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py index 9d1ba58..10cbf54 100644 --- a/tests/test_integral_operators.py +++ b/tests/test_integral_operators.py @@ -1,230 +1,230 @@ # -*- 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 . import tamaas as tm import numpy as np import pytest from numpy.linalg import norm from register_integral_operators import register_kelvin_force, \ set_integration_method @pytest.fixture(params=[tm.integration_method.linear, tm.integration_method.cutoff], ids=['linear', 'cutoff']) def integral_fixture(request): return request.param def test_kelvin_volume_force(integral_fixture): N = 65 E = 1. nu = 0.3 mu = E / (2*(1+nu)) domain = np.array([1.] * 3) omega = 2 * np.pi * np.array([1, 1]) / domain[:2] omega_ = norm(omega) discretization = [N] * 3 model = tm.ModelFactory.createModel(tm.model_type.volume_2d, domain, discretization) model.E = E model.nu = nu register_kelvin_force(model) kelvin = model.operators['kelvin_force'] set_integration_method(kelvin, integral_fixture, 1e-12) coords = [np.linspace(0, domain[i], discretization[i], endpoint=False, dtype=tm.dtype) for i in range(2)]\ + [np.linspace(0, domain[2], discretization[2])] x, y = np.meshgrid(*coords[:2], indexing='ij') displacement = model['displacement'] source = np.zeros_like(displacement) # The integral of forces should stay constant source[N//2, :, :, 2] = np.sin(omega[0]*x) * np.sin(omega[1]*y) * (N-1) kelvin(source, displacement) z = coords[2] - 0.5 z, x, y = np.meshgrid(z, *coords[:2], indexing='ij') solution = np.zeros_like(source) solution[:, :, :, 0] = -np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \ * omega[0]*z*np.cos(omega[0]*x)*np.sin(omega[1]*y) solution[:, :, :, 1] = -np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \ * omega[1]*z*np.sin(omega[0]*x)*np.cos(omega[1]*y) solution[:, :, :, 2] = np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \ * (3-4*nu + omega_*np.abs(z))*np.sin(omega[0]*x)*np.sin(omega[1]*y) error = norm(displacement - solution) / norm(solution) assert error < 5e-2 def test_mindlin(integral_fixture): # Definition of modeled domain # tm.set_log_level(tm.LogLevel.debug) model_type = tm.model_type.volume_2d discretization = [126, 128, 128] system_size = [1., 3., 3.] integration_method = integral_fixture print(integration_method) # Material contants E = 1. # Young's modulus nu = 0.3 # Poisson's ratio h = 0.01 * E # Hardening modulus sigma_y = 0.3 * E # Initial yield stress # Creation of model model = tm.ModelFactory.createModel(model_type, system_size, discretization) model.E = E model.nu = nu mu = E / (2*(1+nu)) lamda = E * nu / ((1+nu) * (1-2*nu)) # Setup for integral operators residual = tm.ModelFactory.createResidual(model, sigma_y, h) - residual.setIntegrationMethod(integration_method, 1e-12) + model.setIntegrationMethod(integration_method, 1e-12) # Coordinates x = np.linspace(0, system_size[1], discretization[1], endpoint=False, dtype=tm.dtype) y = np.linspace(0, system_size[2], discretization[2], endpoint=False, dtype=tm.dtype) z = np.linspace(0, system_size[0], discretization[0], endpoint=True, dtype=tm.dtype) z, x, y = np.meshgrid(z, x, y, indexing='ij') # Inclusion definition a, c = 0.1, 0.2 center = [system_size[1] / 2, system_size[2] / 2, c] r = np.sqrt((x-center[0])**2 + (y-center[1])**2 + (z-center[2])**2) ball = r < a # Eigenstrain definition alpha = 1 # constant isotropic strain beta = (3 * lamda + 2 * mu) * alpha * np.eye(3) eigenstress = np.zeros(discretization + [3, 3], dtype=tm.dtype) eigenstress[ball, ...] = beta eigenstress = tm.compute.to_voigt(eigenstress.reshape(discretization + [9])) # Array references - gradient = residual.getVector() - stress = residual.getStress() + gradient = residual.vector + stress = model["stress"] # Applying operator # mindlin = model.operators["mindlin"] mindlin_gradient = model.operators["mindlin_gradient"] # Not testing displacements yet # mindlin(eigenstress, model['displacement']) # Applying gradient mindlin_gradient(eigenstress, gradient) model.operators['hooke'](gradient, stress) stress -= eigenstress # Normalizing stess as in Mura (1976) T, alpha = 1., 1. beta = alpha * T * (1+nu) / (1-nu) stress *= 1. / (2 * mu * beta) n = discretization[1] // 2 vertical_stress = stress[:, n, n, :] z_all = np.linspace(0, 1, discretization[0], dtype=tm.dtype) sigma_z = np.zeros_like(z_all) sigma_t = np.zeros_like(z_all) inclusion = np.abs(z_all - c) < a # Computing stresses for exterior points z = z_all[~inclusion] R_1 = np.abs(z - c) R_2 = np.abs(z + c) sigma_z[~inclusion] = 2*mu*beta*a**3/3 * ( 1 / R_1**3 - 1 / R_2**3 - 18*z*(z+c) / R_2**5 + 3*(z+c)**2 / R_2**5 - 3*(z-c)**2 / R_1**5 + 30*z*(z+c)**3 / R_2**7 ) sigma_t[~inclusion] = 2*mu*beta*a**3/3 * ( 1 / R_1**3 + (3-8*nu) / R_2**3 - 6*z*(z+c) / R_2**5 + 12*nu*(z+c)**2 / R_2**5 ) # Computing stresses for interior points z = z_all[inclusion] R_1 = np.abs(z - c) R_2 = np.abs(z + c) sigma_z[inclusion] = 2*mu*beta*a**3/3 * ( - 2 / a**3 - 1 / R_2**3 - 18*z*(z+c) / R_2**5 + 3*(z+c)**2 / R_2**5 + 30*z*(z+c)**3 / R_2**7 ) sigma_t[inclusion] = 2*mu*beta*a**3/3 * ( - 2 / a**3 + (3-8*nu) / R_2**3 - 6*z*(z+c) / R_2**5 + 12*nu*(z+c)**2 / R_2**5 ) # This test can be used to debug if False: import matplotlib.pyplot as plt plt.plot(z_all, vertical_stress[:, 2]) plt.plot(z_all, sigma_z / (2 * mu * beta)) plt.figure() plt.plot(z_all, vertical_stress[:, 0]) plt.plot(z_all, sigma_t / (2 * mu * beta)) plt.show() z_error = norm(vertical_stress[:, 2] - sigma_z / (2 * mu * beta)) \ / discretization[0] t_error = norm(vertical_stress[:, 0] - sigma_t / (2 * mu * beta)) \ / discretization[0] assert z_error < 1e-3 and t_error < 1e-3