From 213fff3c8db75f7b48c0b8f1766026e6d5037e58 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Tue, 10 Sep 2019 22:07:27 +0200
Subject: [PATCH] [MatL] CreepBGRa; Faster residual/jac computation.
Also simplify the logic by removing the side effects.
---
MaterialLib/SolidModels/CreepBGRa.cpp | 33 ++++++++++++---------------
1 file changed, 14 insertions(+), 19 deletions(-)
diff --git a/MaterialLib/SolidModels/CreepBGRa.cpp b/MaterialLib/SolidModels/CreepBGRa.cpp
index 46ca8ca412f..2443ddcfe58 100644
--- a/MaterialLib/SolidModels/CreepBGRa.cpp
+++ b/MaterialLib/SolidModels/CreepBGRa.cpp
@@ -74,31 +74,26 @@ CreepBGRa<DisplacementDim>::integrateStress(
dt * constant_coefficient *
std::exp(-Q / (MaterialLib::PhysicalConstant::IdealGasConstant * T));
- // In newton_solver.solve(), the Jacobian is calculated first, and then
- // then comes the assembly of the residue vector. In order to save
- // computation time, the following two variables are set visible in this
- // function. The two variables keep the computation results of the
- // norm of the deviatoric stress and 2bG||s_{n+1}||^{n-1} in
- // update_jacobian, and then they are directly used in update_residual
- // without repeating the computation. Of course, update_jacobian is not a
- // pure function anymore, but has side effects.
- double pow_norm_s_n1_n_minus_one_2b_G = 0.;
-
- KelvinVector s_n1;
- auto const update_jacobian = [&](JacobianMatrix& jacobian) {
- // side effect
- s_n1 = deviatoric_matrix * solution;
+ double const G2b = 2.0 * b * this->_mp.mu(t, x);
+
+ auto const update_jacobian = [&solution, &G2b, &n](JacobianMatrix& jacobian) {
+ auto const& D = Invariants::deviatoric_projection;
+ KelvinVector const s_n1 = D * solution;
double const norm_s_n1 = Invariants::FrobeniusNorm(s_n1);
- double const G2b = 2.0 * b * this->_mp.mu(t, x);
- // side effect
- pow_norm_s_n1_n_minus_one_2b_G = G2b * std::pow(norm_s_n1, n - 1);
+ double const pow_norm_s_n1_n_minus_one_2b_G =
+ G2b * std::pow(norm_s_n1, n - 1);
jacobian = KelvinMatrix::Identity() +
- (pow_norm_s_n1_n_minus_one_2b_G * deviatoric_matrix +
+ (pow_norm_s_n1_n_minus_one_2b_G * D +
(n - 1) * G2b * std::pow(norm_s_n1, n - 3) * s_n1 *
s_n1.transpose());
};
- auto const update_residual = [&](ResidualVectorType& r) {
+ auto const update_residual = [&solution, &G2b, &n,
+ &sigma_try](ResidualVectorType& r) {
+ KelvinVector const s_n1 = Invariants::deviatoric_projection * solution;
+ double const norm_s_n1 = Invariants::FrobeniusNorm(s_n1);
+ double const pow_norm_s_n1_n_minus_one_2b_G =
+ G2b * std::pow(norm_s_n1, n - 1);
r = solution - sigma_try + pow_norm_s_n1_n_minus_one_2b_G * s_n1;
};
--
GitLab