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