Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84212753
__init__.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 21, 09:54
Size
6 KB
Mime Type
text/x-python
Expires
Mon, Sep 23, 09:54 (2 d)
Engine
blob
Format
Raw Data
Handle
20959663
Attached To
rTAMAAS tamaas
__init__.py
View Options
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# 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 <https://www.gnu.org/licenses/>.
"""
Pulling solvers to nonlinear_solvers module
"""
from
..
import
EPSolver
import
numpy
as
np
from
scipy.sparse.linalg
import
LinearOperator
,
lgmres
from
scipy.linalg
import
norm
from
scipy.optimize
import
newton_krylov
,
root
try
:
from
.akantu_solver
import
*
# noqa
except
ImportError
:
pass
__all__
=
[
'NLNoConvergence'
,
'DFSANESolver'
,
'NewtonKrylovSolver'
,
'ToleranceManager'
]
class
NLNoConvergence
(
Exception
):
"""Convergence not reached exception"""
pass
class
NewtonRaphsonSolver
(
EPSolver
):
"""
Newton-Raphson based on CTO
"""
def
__init__
(
self
,
residual
,
model
):
super
(
NewtonRaphsonSolver
,
self
)
.
__init__
(
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
ScipySolver
(
EPSolver
):
"""
Base class for solvers wrapping SciPy routines
"""
def
__init__
(
self
,
residual
,
model
,
callback
=
None
):
super
(
ScipySolver
,
self
)
.
__init__
(
residual
)
self
.
model
=
model
self
.
callback
=
callback
self
.
_x
=
self
.
getStrainIncrement
()
self
.
_residual
=
self
.
getResidual
()
self
.
options
=
{
'ftol'
:
0
,
'fatol'
:
1e-9
}
@property
def
tolerance
(
self
):
return
self
.
options
[
'fatol'
]
@tolerance.setter
def
set_tolerance
(
self
,
v
):
self
.
options
[
'fatol'
]
=
v
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
.
_x
[
...
]
=
0
class
NewtonKrylovSolver
(
ScipySolver
):
"""
Solve using a finite-difference Newton-Krylov method
"""
def
__init__
(
self
,
residual
,
model
,
callback
=
None
):
ScipySolver
.
__init__
(
self
,
residual
,
model
,
callback
)
def
scipy_solve
(
self
,
compute_residual
):
try
:
return
newton_krylov
(
compute_residual
,
self
.
_x
,
f_tol
=
self
.
options
[
'fatol'
],
verbose
=
True
,
callback
=
self
.
callback
)
except
Exception
as
e
:
raise
NLNoConvergence
(
e
.
what
())
class
DFSANESolver
(
ScipySolver
):
"""
Solve using a spectral residual jacobianless method
"""
def
__init__
(
self
,
residual
,
model
,
callback
=
None
):
ScipySolver
.
__init__
(
self
,
residual
,
model
,
callback
)
def
scipy_solve
(
self
,
compute_residual
):
solution
=
root
(
compute_residual
,
self
.
_x
,
method
=
'df-sane'
,
options
=
self
.
options
,
callback
=
self
.
callback
)
print
(
"DF-SANE: {} ({} iterations, {})"
.
format
(
solution
.
message
,
solution
.
nit
,
self
.
options
))
if
not
solution
.
success
:
raise
NLNoConvergence
(
"DF-SANE did not converge"
)
return
solution
.
x
.
copy
()
class
ToleranceManager
:
"""
Decorator to manage tolerance of non-linear solver
"""
def
__init__
(
self
,
start
,
end
,
rate
,
key
=
'fatol'
):
self
.
start
=
start
/
rate
# just anticipating first multiplication
self
.
end
=
end
self
.
rate
=
rate
self
.
key
=
key
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
[
self
.
key
]
=
self
.
start
def
scipy_solve
(
obj
,
*
args
,
**
kwargs
):
ftol
=
obj
.
options
[
self
.
key
]
ftol
*=
self
.
rate
obj
.
options
[
self
.
key
]
=
max
(
ftol
,
self
.
end
)
return
orig_solve
(
obj
,
*
args
,
**
kwargs
)
def
updateState
(
obj
,
*
args
,
**
kwargs
):
obj
.
options
[
self
.
key
]
=
self
.
start
return
orig_updateState
(
obj
,
*
args
,
**
kwargs
)
cls
.
__init__
=
__init__
cls
.
scipy_solve
=
scipy_solve
cls
.
updateState
=
updateState
return
cls
Event Timeline
Log In to Comment