Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/python/tamaas/nonlinear_solvers.py
9f1e5f8c78d4features/declarative

Last Change
Raw File

nonlinear_solvers.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/>.
import tamaas as tm
from scipy.optimize import root
import numpy as np
import scipy.interpolate as ip
from scipy.sparse.linalg import LinearOperator, lgmres
from scipy.linalg import norm
from scipy.optimize import newton_krylov
def makePeriodic(data):
data = np.append(data, np.expand_dims(data[0, ...], axis=0), axis=0)
data = np.append(data, np.expand_dims(data[:, 0, ...], axis=1), axis=1)
return data
class NonSolver(tm.EPSolver):
def __init__(self, residual, model):
tm.EPSolver.__init__(self, residual)
self.model = model
self.elastic_strain_increment = self.getStrainIncrement()
self.last_increment = self.elastic_strain_increment.copy()
def solve(self):
self.elastic_strain_increment[...] = 0
self.getResidual().computeResidual(self.elastic_strain_increment)
self.elastic_strain_increment[...] = self.getResidual().getVector()
self.last_increment *= -1
self.last_increment += self.elastic_strain_increment
self.getResidual().computeResidualDisplacement(
self.elastic_strain_increment)
class ScipySolver(tm.EPSolver):
def __init__(self, residual, model):
tm.EPSolver.__init__(self, 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):
# 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):
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-10, 'fatol': 1e-10}
def scipy_solve(self, compute_residual):
solution = root(compute_residual, self._x, method='df-sane',
options=self.options)
print("DF-SANE: {} ({} iterations)".format(solution.message,
solution.nit))
return solution.x.copy()
class NewtonRaphsonSolver(tm.EPSolver):
def __init__(self, residual, model):
tm.EPSolver.__init__(self, residual)
self.model = model
self._residual = self.getResidual()
self._x = self.getStrainIncrement()
self.last_solution = np.zeros_like(self._x)
def applyTangent(self, _input, strain):
_input = np.reshape(_input, strain.shape)
out = np.zeros_like(_input)
self._residual.applyTangent(out, _input, strain)
return out
def solve(self):
self._x[...] = 0
# For initial guess, compute the strain due to boundary tractions
self._residual.computeResidual(self._x)
self._x[...] = self._residual.getVector()
res_vec = np.ravel(self._residual.getVector())
initial_norm = norm(res_vec)
self._residual.computeResidual(self._x)
n_max = 100
i = 0
while norm(res_vec) / initial_norm > 1e-10 and i < n_max:
# Making linear tangent
tangent = LinearOperator(
(self._x.size, self._x.size),
matvec=lambda x: self.applyTangent(x, self._x))
res_vec *= -1
correction, ok = lgmres(tangent, res_vec, tol=1e-12, maxiter=100)
if ok != 0:
raise Exception("LGMRES not converged")
self._x += np.reshape(correction, self._x.shape)
self._residual.computeResidual(self._x)
i += 1
print("Solved Newton-Raphson in {} iterations".format(i))
# Computing the strain correction
self.last_solution *= -1
self.last_solution += self._x
# Computing displacements
self._residual.computeResidualDisplacement(self.last_solution)
self.last_solution[...] = self._x[...]
class FromTraction:
"""Functor for Neumann boundary conditions"""
def __init__(self, interpolation):
self.interpolation = interpolation
def operator(self, quad_point, dual, coord, normals):
val = self.interpolation(*(coord[:-1]))
dual[2] = val[0]
try:
import akantu as ak
class PoincareFEM_plast(tm.Residual):
def __init__(self, tm_model, ak_model):
tm.Residual.__init__(self, tm_model)
self.tm_model = tm_model
self.discretization = tm_model.getDiscretization().copy()
self.previous_pressure = np.zeros(self.discretization[1:])
self.residual_displacement = np.zeros(self.discretization[1:])
self.discretization[1] += 1
self.discretization[2] += 1
self.coordinates = [np.linspace(0, L, N)
for L, N in zip(tm_model.getSystemSize(),
self.discretization)]
self.model = ak_model
self.bot_face = np.where(self.model.getMesh().getNodes()[:, 2] == 0)
self.ids = \
{'inelastic_strain':
'previous_solid_mechanics_model:0:'
+ 'plastic_linear_isotropic_hardening:inelastic_strain',
'iso_hardening':
'previous_solid_mechanics_model:0:'
+ 'plastic_linear_isotropic_hardening:iso_hardening'}
self.savePreviousState()
def apply(self, field, out=None):
interpolation = ip.interp2d(self.coordinates[1],
self.coordinates[2],
makePeriodic(field))
functor = FromTraction(interpolation)
self.model.getExternalForce()[...] = 0
self.model.applyNeumannBC(functor, "bot")
self.model.solveStep()
displacements = self.model.getDisplacement()[self.bot_face, 2]
if out is not None:
# All these shenanigans because of node ordering
bot_points = self.model.getMesh().getNodes()[self.bot_face, :-1]
bot_points = bot_points.reshape(
(self.discretization[1] * self.discretization[2], 2))
x, y = np.meshgrid(*self.coordinates[1:], indexing='ij')
disp = ip.griddata(bot_points, displacements.T, (x, y))
out[...] = disp.reshape(self.discretization[1:])[:-1, :-1]
out -= out.mean()
def savePreviousState(self):
self.inelastic_strain = \
self.model.getMaterial(0).getArrayReal(
self.ids['inelastic_strain'], ak._tetrahedron_4).copy()
self.hardening = \
self.model.getMaterial(0).getArrayReal(
self.ids['iso_hardening'], ak._tetrahedron_4).copy()
self.previous_pressure[...] = self.tm_model.getTraction()[:, :, 2]
def restorePreviousState(self):
inelastic_strain = \
self.model.getMaterial(0).getArrayReal(
self.ids['inelastic_strain'], ak._tetrahedron_4)
hardening = \
self.model.getMaterial(0).getArrayReal(
self.ids['iso_hardening'], ak._tetrahedron_4)
inelastic_strain[...] = self.inelastic_strain
hardening[...] = self.hardening
def updateState(self, increment):
"""
Update the plastic deformation state
- Here we disregard the increment
"""
self.restorePreviousState()
self.apply(self.tm_model.getTraction()[:, :, 2])
self.apply(np.zeros_like(self.residual_displacement),
self.residual_displacement)
self.apply(self.tm_model.getTraction()[:, :, 2])
self.savePreviousState()
def getStrainIncrement(self):
return np.array()
class AkantuSolver(tm.EPSolver):
"""Non-linear finite elements solver"""
def __init__(self, tm_model, material_file, mesh_name):
"""Initialize Akantu SolidMechanicsModel"""
ak.parseInput(material_file)
self.tm_model = tm_model
self.name = mesh_name
# Loading mesh
self.mesh = ak.Mesh(3)
self.mesh.read(mesh_name)
# Creating Akantu model
self.model = ak.SolidMechanicsModel(self.mesh)
self.model.initFull(_analysis_method=ak._static)
# Poincare operator/dummy_residual
self.poincare = PoincareFEM_plast(tm_model, self.model)
tm.EPSolver.__init__(self, self.poincare)
# Setting boundary conditions
position = self.mesh.getNodes()
blocked_dofs = self.model.getBlockedDOFs()
epsilon = 1e-10
for i in range(self.mesh.getNbNodes()):
if np.abs(position[i, 2] - 1) < epsilon: # Blocking top surface
blocked_dofs[i, 2] = True
if np.abs(position[i, 0]) < epsilon:
blocked_dofs[i, 1] = True
if np.abs(position[i, 1]) < epsilon:
blocked_dofs[i, 0] = True
if np.abs(position[i, 0]) < epsilon \
or np.abs(position[i, 0] - 1.) < epsilon:
blocked_dofs[i, 0] = True
if np.abs(position[i, 1]) < epsilon \
or np.abs(position[i, 1] - 1.) < epsilon:
blocked_dofs[i, 1] = True
solver = self.model.getNonLinearSolver("static")
solver.set('max_iterations', 200)
solver.set('threshold', 1e-6)
solver.set('convergence_type', ak._scc_residual)
self.surface_correction = tm_model.getDisplacement()[0, :, :, 2]
self.previous_solution = np.zeros_like(self.surface_correction)
def solve(self):
"""Apply and release pressure to get the residual displacements"""
# this is the pressure increment
pressure = self.tm_model.getTraction()[:, :, 2]
self.poincare.restorePreviousState()
self.poincare.apply(pressure + self.poincare.previous_pressure)
residual_disp = np.zeros_like(self.surface_correction)
self.poincare.apply(np.zeros_like(pressure), residual_disp)
self.surface_correction[...] = \
residual_disp - self.poincare.residual_displacement
def reset(self):
self.previous_solution[...] = 0
except ImportError:
print("Akantu module not found: AkantuSolver unavailable")
c4science · Help