Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78058610
solvers.cc
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
Sun, Aug 18, 04:00
Size
17 KB
Mime Type
text/x-c
Expires
Tue, Aug 20, 04:00 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19977421
Attached To
rMUSPECTRE µSpectre
solvers.cc
View Options
/**
* file solvers.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 24 Apr 2018
*
* @brief implementation of dynamic newton-cg solver
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with µSpectre; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* * Boston, MA 02111-1307, USA.
*
* Additional permission under GNU GPL version 3 section 7
*
* If you modify this Program, or any covered work, by linking or combining it
* with proprietary FFT implementations or numerical libraries, containing parts
* covered by the terms of those libraries' licenses, the licensors of this
* Program grant you additional permission to convey the resulting work.
*/
#include "solver/solvers.hh"
#include "common/iterators.hh"
#include <Eigen/Dense>
#include <iomanip>
namespace
muSpectre
{
Eigen
::
IOFormat
format
(
Eigen
::
FullPrecision
,
0
,
", "
,
",
\n
"
,
"["
,
"]"
,
"["
,
"]"
);
//--------------------------------------------------------------------------//
std
::
vector
<
OptimizeResult
>
newton_cg
(
Cell
&
cell
,
const
LoadSteps_t
&
load_steps
,
SolverBase
&
solver
,
Real
newton_tol
,
Real
equil_tol
,
Dim_t
verbose
)
{
const
Communicator
&
comm
=
cell
.
get_communicator
();
using
Vector_t
=
Eigen
::
Matrix
<
Real
,
Eigen
::
Dynamic
,
1
>
;
using
Matrix_t
=
Eigen
::
Matrix
<
Real
,
Eigen
::
Dynamic
,
Eigen
::
Dynamic
>
;
// Corresponds to symbol δF or δε
Vector_t
incrF
(
cell
.
get_nb_dof
());
// field to store the rhs for cg calculations
Vector_t
rhs
(
cell
.
get_nb_dof
());
solver
.
initialise
();
size_t
count_width
{};
const
auto
form
{
cell
.
get_formulation
()};
std
::
string
strain_symb
{};
if
(
verbose
>
0
&&
comm
.
rank
()
==
0
)
{
// setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
std
::
cout
<<
"Newton-"
<<
solver
.
get_name
()
<<
" for "
;
switch
(
form
)
{
case
Formulation
::
small_strain:
{
strain_symb
=
"ε"
;
std
::
cout
<<
"small"
;
break
;
}
case
Formulation
::
finite_strain:
{
strain_symb
=
"F"
;
std
::
cout
<<
"finite"
;
break
;
}
default
:
throw
SolverError
(
"unknown formulation"
);
break
;
}
std
::
cout
<<
" strain with"
<<
std
::
endl
<<
"newton_tol = "
<<
newton_tol
<<
", cg_tol = "
<<
solver
.
get_tol
()
<<
" maxiter = "
<<
solver
.
get_maxiter
()
<<
" and "
<<
strain_symb
<<
" from "
<<
strain_symb
<<
"₁ ="
<<
std
::
endl
<<
load_steps
.
front
()
<<
std
::
endl
<<
" to "
<<
strain_symb
<<
"ₙ ="
<<
std
::
endl
<<
load_steps
.
back
()
<<
std
::
endl
<<
"in increments of Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
(
load_steps
.
back
()
-
load_steps
.
front
())
/
load_steps
.
size
()
<<
std
::
endl
;
count_width
=
size_t
(
std
::
log10
(
solver
.
get_maxiter
()))
+
1
;
}
auto
shape
{
cell
.
get_strain_shape
()};
switch
(
form
)
{
case
Formulation
::
finite_strain:
{
cell
.
set_uniform_strain
(
Matrix_t
::
Identity
(
shape
[
0
],
shape
[
1
]));
for
(
const
auto
&
delF
:
load_steps
)
{
if
(
not
((
delF
.
rows
()
==
shape
[
0
])
and
(
delF
.
cols
()
==
shape
[
1
])))
{
std
::
stringstream
err
{};
err
<<
"Load increments need to be given in "
<<
shape
[
0
]
<<
"×"
<<
shape
[
1
]
<<
" matrices, but I got a "
<<
delF
.
rows
()
<<
"×"
<<
delF
.
cols
()
<<
" matrix:"
<<
std
::
endl
<<
delF
;
throw
SolverError
(
err
.
str
());
}
}
break
;
}
case
Formulation
::
small_strain:
{
cell
.
set_uniform_strain
(
Matrix_t
::
Zero
(
shape
[
0
],
shape
[
1
]));
for
(
const
auto
&
delF
:
load_steps
)
{
if
(
not
((
delF
.
rows
()
==
shape
[
0
])
and
(
delF
.
cols
()
==
shape
[
1
])))
{
std
::
stringstream
err
{};
err
<<
"Load increments need to be given in "
<<
shape
[
0
]
<<
"×"
<<
shape
[
1
]
<<
" matrices, but I got a "
<<
delF
.
rows
()
<<
"×"
<<
delF
.
cols
()
<<
" matrix:"
<<
std
::
endl
<<
delF
;
throw
SolverError
(
err
.
str
());
}
if
(
not
check_symmetry
(
delF
))
{
throw
SolverError
(
"all Δε must be symmetric!"
);
}
}
break
;
}
default
:
throw
SolverError
(
"Unknown strain measure"
);
break
;
}
// initialise return value
std
::
vector
<
OptimizeResult
>
ret_val
{};
// storage for the previous mean strain (to compute ΔF or Δε)
Matrix_t
previous_macro_strain
{
load_steps
.
back
().
Zero
(
shape
[
0
],
shape
[
1
])};
auto
F
{
cell
.
get_strain_vector
()};
//! incremental loop
for
(
const
auto
&
tup
:
akantu
::
enumerate
(
load_steps
))
{
const
auto
&
strain_step
{
std
::
get
<
0
>
(
tup
)};
const
auto
&
macro_strain
{
std
::
get
<
1
>
(
tup
)};
if
((
verbose
>
0
)
and
(
comm
.
rank
()
==
0
))
{
std
::
cout
<<
"at Load step "
<<
std
::
setw
(
count_width
)
<<
strain_step
+
1
<<
std
::
endl
;
}
using
StrainMap_t
=
RawFieldMap
<
Eigen
::
Map
<
Eigen
::
MatrixXd
>>
;
for
(
auto
&&
strain
:
StrainMap_t
(
F
,
shape
[
0
],
shape
[
1
]))
{
strain
+=
macro_strain
-
previous_macro_strain
;
}
std
::
string
message
{
"Has not converged"
};
Real
incr_norm
{
2
*
newton_tol
},
grad_norm
{
1
};
Real
stress_norm
{
2
*
equil_tol
};
bool
has_converged
{
false
};
auto
convergence_test
=
[
&
incr_norm
,
&
grad_norm
,
&
newton_tol
,
&
stress_norm
,
&
equil_tol
,
&
message
,
&
has_converged
]()
{
bool
incr_test
=
incr_norm
/
grad_norm
<=
newton_tol
;
bool
stress_test
=
stress_norm
<
equil_tol
;
if
(
incr_test
)
{
message
=
"Residual tolerance reached"
;
}
else
if
(
stress_test
)
{
message
=
"Reached stress divergence tolerance"
;
}
has_converged
=
incr_test
||
stress_test
;
return
has_converged
;
};
Uint
newt_iter
{
0
};
for
(;
newt_iter
<
solver
.
get_maxiter
()
&&
!
has_converged
;
++
newt_iter
)
{
// obtain material response
auto
res_tup
{
cell
.
evaluate_stress_tangent
()};
auto
&
P
{
std
::
get
<
0
>
(
res_tup
)};
rhs
=
-
P
;
cell
.
apply_projection
(
rhs
);
stress_norm
=
std
::
sqrt
(
comm
.
sum
(
rhs
.
squaredNorm
()));
if
(
convergence_test
())
{
break
;
}
try
{
incrF
=
solver
.
solve
(
rhs
);
}
catch
(
ConvergenceError
&
error
)
{
std
::
stringstream
err
{};
err
<<
"Failure at load step "
<<
strain_step
+
1
<<
" of "
<<
load_steps
.
size
()
<<
". In Newton-Raphson step "
<<
newt_iter
<<
":"
<<
std
::
endl
<<
error
.
what
()
<<
std
::
endl
<<
"The applied boundary condition is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
<<
std
::
endl
<<
"and the load increment is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
-
previous_macro_strain
<<
std
::
endl
;
throw
ConvergenceError
(
err
.
str
());
}
F
+=
incrF
;
incr_norm
=
std
::
sqrt
(
comm
.
sum
(
incrF
.
squaredNorm
()));
grad_norm
=
std
::
sqrt
(
comm
.
sum
(
F
.
squaredNorm
()));
if
((
verbose
>
1
)
and
(
comm
.
rank
()
==
0
))
{
std
::
cout
<<
"at Newton step "
<<
std
::
setw
(
count_width
)
<<
newt_iter
<<
", |δ"
<<
strain_symb
<<
"|/|Δ"
<<
strain_symb
<<
"| = "
<<
std
::
setw
(
17
)
<<
incr_norm
/
grad_norm
<<
", tol = "
<<
newton_tol
<<
std
::
endl
;
if
(
verbose
-
1
>
1
)
{
std
::
cout
<<
"<"
<<
strain_symb
<<
"> ="
<<
std
::
endl
<<
StrainMap_t
(
F
,
shape
[
0
],
shape
[
1
]).
mean
()
<<
std
::
endl
;
}
}
convergence_test
();
}
if
(
newt_iter
==
solver
.
get_maxiter
())
{
std
::
stringstream
err
{};
err
<<
"Failure at load step "
<<
strain_step
+
1
<<
" of "
<<
load_steps
.
size
()
<<
". Newton-Raphson failed to converge. "
<<
"The applied boundary condition is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
<<
std
::
endl
<<
"and the load increment is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
-
previous_macro_strain
<<
std
::
endl
;
throw
ConvergenceError
(
err
.
str
());
}
// update previous macroscopic strain
previous_macro_strain
=
macro_strain
;
// store results
ret_val
.
emplace_back
(
OptimizeResult
{
F
,
cell
.
get_stress_vector
(),
convergence_test
(),
Int
(
convergence_test
()),
message
,
newt_iter
,
solver
.
get_counter
(),
form
});
// store history variables for next load increment
cell
.
save_history_variables
();
}
return
ret_val
;
}
//----------------------------------------------------------------------------//
std
::
vector
<
OptimizeResult
>
de_geus
(
Cell
&
cell
,
const
LoadSteps_t
&
load_steps
,
SolverBase
&
solver
,
Real
newton_tol
,
Real
equil_tol
,
Dim_t
verbose
)
{
const
Communicator
&
comm
=
cell
.
get_communicator
();
using
Vector_t
=
Eigen
::
Matrix
<
Real
,
Eigen
::
Dynamic
,
1
>
;
using
Matrix_t
=
Eigen
::
Matrix
<
Real
,
Eigen
::
Dynamic
,
Eigen
::
Dynamic
>
;
// Corresponds to symbol δF or δε
Vector_t
incrF
(
cell
.
get_nb_dof
());
// Corresponds to symbol ΔF or Δε
Vector_t
DeltaF
(
cell
.
get_nb_dof
());
// field to store the rhs for cg calculations
Vector_t
rhs
(
cell
.
get_nb_dof
());
solver
.
initialise
();
size_t
count_width
{};
const
auto
form
{
cell
.
get_formulation
()};
std
::
string
strain_symb
{};
if
(
verbose
>
0
&&
comm
.
rank
()
==
0
)
{
// setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
std
::
cout
<<
"de Geus-"
<<
solver
.
get_name
()
<<
" for "
;
switch
(
form
)
{
case
Formulation
::
small_strain:
{
strain_symb
=
"ε"
;
std
::
cout
<<
"small"
;
break
;
}
case
Formulation
::
finite_strain:
{
strain_symb
=
"F"
;
std
::
cout
<<
"finite"
;
break
;
}
default
:
throw
SolverError
(
"unknown formulation"
);
break
;
}
std
::
cout
<<
" strain with"
<<
std
::
endl
<<
"newton_tol = "
<<
newton_tol
<<
", cg_tol = "
<<
solver
.
get_tol
()
<<
" maxiter = "
<<
solver
.
get_maxiter
()
<<
" and "
<<
strain_symb
<<
" from "
<<
strain_symb
<<
"₁ ="
<<
std
::
endl
<<
load_steps
.
front
()
<<
std
::
endl
<<
" to "
<<
strain_symb
<<
"ₙ ="
<<
std
::
endl
<<
load_steps
.
back
()
<<
std
::
endl
<<
"in increments of Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
(
load_steps
.
back
()
-
load_steps
.
front
())
/
load_steps
.
size
()
<<
std
::
endl
;
count_width
=
size_t
(
std
::
log10
(
solver
.
get_maxiter
()))
+
1
;
}
auto
shape
{
cell
.
get_strain_shape
()};
Matrix_t
default_strain_val
{};
switch
(
form
)
{
case
Formulation
::
finite_strain:
{
cell
.
set_uniform_strain
(
Matrix_t
::
Identity
(
shape
[
0
],
shape
[
1
]));
for
(
const
auto
&
delF
:
load_steps
)
{
auto
rows
=
delF
.
rows
();
auto
cols
=
delF
.
cols
();
if
(
not
((
rows
==
shape
[
0
])
and
(
cols
==
shape
[
1
])))
{
std
::
stringstream
err
{};
err
<<
"Load increments need to be given in "
<<
shape
[
0
]
<<
"×"
<<
shape
[
1
]
<<
" matrices, but I got a "
<<
delF
.
rows
()
<<
"×"
<<
delF
.
cols
()
<<
" matrix:"
<<
std
::
endl
<<
delF
;
throw
SolverError
(
err
.
str
());
}
}
break
;
}
case
Formulation
::
small_strain:
{
cell
.
set_uniform_strain
(
Matrix_t
::
Zero
(
shape
[
0
],
shape
[
1
]));
for
(
const
auto
&
delF
:
load_steps
)
{
if
(
not
((
delF
.
rows
()
==
shape
[
0
])
and
(
delF
.
cols
()
==
shape
[
1
])))
{
std
::
stringstream
err
{};
err
<<
"Load increments need to be given in "
<<
shape
[
0
]
<<
"×"
<<
shape
[
1
]
<<
" matrices, but I got a "
<<
delF
.
rows
()
<<
"×"
<<
delF
.
cols
()
<<
" matrix:"
<<
std
::
endl
<<
delF
;
throw
SolverError
(
err
.
str
());
}
if
(
not
check_symmetry
(
delF
))
{
throw
SolverError
(
"all Δε must be symmetric!"
);
}
}
break
;
}
default
:
throw
SolverError
(
"Unknown strain measure"
);
break
;
}
// initialise return value
std
::
vector
<
OptimizeResult
>
ret_val
{};
// storage for the previous mean strain (to compute ΔF or Δε)
Matrix_t
previous_macro_strain
{
load_steps
.
back
().
Zero
(
shape
[
0
],
shape
[
1
])};
auto
F
{
cell
.
get_strain_vector
()};
//! incremental loop
for
(
const
auto
&
tup
:
akantu
::
enumerate
(
load_steps
))
{
const
auto
&
strain_step
{
std
::
get
<
0
>
(
tup
)};
const
auto
&
macro_strain
{
std
::
get
<
1
>
(
tup
)};
using
StrainMap_t
=
RawFieldMap
<
Eigen
::
Map
<
Eigen
::
MatrixXd
>>
;
if
((
verbose
>
0
)
and
(
comm
.
rank
()
==
0
))
{
std
::
cout
<<
"at Load step "
<<
std
::
setw
(
count_width
)
<<
strain_step
+
1
<<
", "
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
(
macro_strain
+
default_strain_val
).
format
(
format
)
<<
std
::
endl
;
}
std
::
string
message
{
"Has not converged"
};
Real
incr_norm
{
2
*
newton_tol
},
grad_norm
{
1
};
Real
stress_norm
{
2
*
equil_tol
};
bool
has_converged
{
false
};
auto
convergence_test
=
[
&
incr_norm
,
&
grad_norm
,
&
newton_tol
,
&
stress_norm
,
&
equil_tol
,
&
message
,
&
has_converged
]()
{
bool
incr_test
=
incr_norm
/
grad_norm
<=
newton_tol
;
bool
stress_test
=
stress_norm
<
equil_tol
;
if
(
incr_test
)
{
message
=
"Residual tolerance reached"
;
}
else
if
(
stress_test
)
{
message
=
"Reached stress divergence tolerance"
;
}
has_converged
=
incr_test
||
stress_test
;
return
has_converged
;
};
Uint
newt_iter
{
0
};
for
(;
((
newt_iter
<
solver
.
get_maxiter
())
and
(
!
has_converged
))
or
(
newt_iter
<
2
);
++
newt_iter
)
{
// obtain material response
auto
res_tup
{
cell
.
evaluate_stress_tangent
()};
auto
&
P
{
std
::
get
<
0
>
(
res_tup
)};
try
{
if
(
newt_iter
==
0
)
{
for
(
auto
&&
strain
:
StrainMap_t
(
DeltaF
,
shape
[
0
],
shape
[
1
]))
{
strain
=
macro_strain
-
previous_macro_strain
;
}
rhs
=
-
cell
.
evaluate_projected_directional_stiffness
(
DeltaF
);
F
+=
DeltaF
;
stress_norm
=
std
::
sqrt
(
comm
.
sum
(
rhs
.
matrix
().
squaredNorm
()));
if
(
stress_norm
<
equil_tol
)
{
incrF
.
setZero
();
}
else
{
incrF
=
solver
.
solve
(
rhs
);
}
}
else
{
rhs
=
-
P
;
cell
.
apply_projection
(
rhs
);
stress_norm
=
std
::
sqrt
(
comm
.
sum
(
rhs
.
matrix
().
squaredNorm
()));
if
(
stress_norm
<
equil_tol
)
{
incrF
.
setZero
();
}
else
{
incrF
=
solver
.
solve
(
rhs
);
}
}
}
catch
(
ConvergenceError
&
error
)
{
std
::
stringstream
err
{};
err
<<
"Failure at load step "
<<
strain_step
+
1
<<
" of "
<<
load_steps
.
size
()
<<
". In Newton-Raphson step "
<<
newt_iter
<<
":"
<<
std
::
endl
<<
error
.
what
()
<<
std
::
endl
<<
"The applied boundary condition is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
<<
std
::
endl
<<
"and the load increment is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
-
previous_macro_strain
<<
std
::
endl
;
throw
ConvergenceError
(
err
.
str
());
}
F
+=
incrF
;
incr_norm
=
std
::
sqrt
(
comm
.
sum
(
incrF
.
squaredNorm
()));
grad_norm
=
std
::
sqrt
(
comm
.
sum
(
F
.
squaredNorm
()));
if
((
verbose
>
0
)
and
(
comm
.
rank
()
==
0
))
{
std
::
cout
<<
"at Newton step "
<<
std
::
setw
(
count_width
)
<<
newt_iter
<<
", |δ"
<<
strain_symb
<<
"|/|Δ"
<<
strain_symb
<<
"| = "
<<
std
::
setw
(
17
)
<<
incr_norm
/
grad_norm
<<
", tol = "
<<
newton_tol
<<
std
::
endl
;
if
(
verbose
-
1
>
1
)
{
std
::
cout
<<
"<"
<<
strain_symb
<<
"> ="
<<
std
::
endl
<<
StrainMap_t
(
F
,
shape
[
0
],
shape
[
1
]).
mean
()
<<
std
::
endl
;
}
}
convergence_test
();
}
if
(
newt_iter
==
solver
.
get_maxiter
())
{
std
::
stringstream
err
{};
err
<<
"Failure at load step "
<<
strain_step
+
1
<<
" of "
<<
load_steps
.
size
()
<<
". Newton-Raphson failed to converge. "
<<
"The applied boundary condition is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
<<
std
::
endl
<<
"and the load increment is Δ"
<<
strain_symb
<<
" ="
<<
std
::
endl
<<
macro_strain
-
previous_macro_strain
<<
std
::
endl
;
throw
ConvergenceError
(
err
.
str
());
}
// update previous macroscopic strain
previous_macro_strain
=
macro_strain
;
// store results
ret_val
.
emplace_back
(
OptimizeResult
{
F
,
cell
.
get_stress_vector
(),
convergence_test
(),
Int
(
convergence_test
()),
message
,
newt_iter
,
solver
.
get_counter
(),
form
});
// store history variables for next load increment
cell
.
save_history_variables
();
}
return
ret_val
;
}
}
// namespace muSpectre
Event Timeline
Log In to Comment