https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: 04c710c35371ded12eb6f1cd7e7fa6ddfe0d7f5b authored by Lars Bilke on 05 March 2021, 08:21:32 UTC
Merge branch 'snakemake-vtkdiff-followup' into 'master'
Merge branch 'snakemake-vtkdiff-followup' into 'master'
Tip revision: 04c710c
Lubby2.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#include "Lubby2.h"
#include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
namespace MPL = MaterialPropertyLib;
namespace MaterialLib
{
namespace Solids
{
namespace Lubby2
{
/// Calculates the 18x6 derivative of the residuals with respect to total
/// strain.
template <int DisplacementDim>
Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
Lubby2<DisplacementDim>::KelvinVectorSize>
calculatedGdEBurgers()
{
Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
Lubby2<DisplacementDim>::KelvinVectorSize>
dGdE =
Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
Lubby2<DisplacementDim>::KelvinVectorSize>::Zero();
dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
Lubby2<DisplacementDim>::KelvinVectorSize>()
.diagonal()
.setConstant(-2.);
return dGdE;
}
template <int DisplacementDim, typename LinearSolver>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> tangentStiffnessA(
double const GM0, double const KM0, LinearSolver const& linear_solver)
{
// Calculate dGdE for time step
auto const dGdE = calculatedGdEBurgers<DisplacementDim>();
// Consistent tangent from local Newton iteration of material
// functionals.
// Only the upper left block is relevant for the global tangent.
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
KelvinMatrix const dzdE =
linear_solver.solve(-dGdE)
.template topLeftCorner<KelvinVectorSize, KelvinVectorSize>();
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
auto const& P_sph = Invariants::spherical_projection;
auto const& P_dev = Invariants::deviatoric_projection;
KelvinMatrix C = GM0 * dzdE * P_dev + 3. * KM0 * P_sph;
return C;
};
template <int DisplacementDim>
std::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
std::unique_ptr<typename MechanicsBase<
DisplacementDim>::MaterialStateVariables>,
typename Lubby2<DisplacementDim>::KelvinMatrix>>
Lubby2<DisplacementDim>::integrateStress(
MaterialPropertyLib::VariableArray const& variable_array_prev,
MaterialPropertyLib::VariableArray const& variable_array, double const t,
ParameterLib::SpatialPosition const& x, double const dt,
typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
material_state_variables) const
{
auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]);
auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
variable_array_prev[static_cast<int>(
MPL::Variable::mechanical_strain)]);
auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
variable_array_prev[static_cast<int>(MPL::Variable::stress)]);
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
assert(dynamic_cast<MaterialStateVariables const*>(
&material_state_variables) != nullptr);
MaterialStateVariables state(
static_cast<MaterialStateVariables const&>(material_state_variables));
state.setInitialConditions();
auto local_lubby2_properties =
detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
// calculation of deviatoric parts
auto const& P_dev = Invariants::deviatoric_projection;
KelvinVector const eps_m_d_i = P_dev * eps_m;
KelvinVector const eps_m_d_t = P_dev * eps_m_prev;
// initial guess as elastic predictor.
KelvinVector sigd_j = 2.0 * (eps_m_d_i - state.eps_M_t - state.eps_K_t);
// Note: sigd_t contains dimensionless stresses!
KelvinVector sigd_t = P_dev * sigma_prev / local_lubby2_properties.GM0;
// Calculate effective stress and update material properties
double sig_eff = Invariants::equivalentStress(sigd_j);
local_lubby2_properties.update(sig_eff);
using LocalJacobianMatrix =
Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
Eigen::RowMajor>;
// Linear solver for the newton loop is required after the loop with the
// same matrix. This saves one decomposition.
Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
// Different solvers are available for the solution of the local system.
// TODO Make the following choice of linear solvers available from the
// input file configuration:
// K_loc.partialPivLu().solve(-res_loc);
// K_loc.fullPivLu().solve(-res_loc);
// K_loc.householderQr().solve(-res_loc);
// K_loc.colPivHouseholderQr().solve(res_loc);
// K_loc.fullPivHouseholderQr().solve(-res_loc);
// K_loc.llt().solve(-res_loc);
// K_loc.ldlt().solve(-res_loc);
LocalJacobianMatrix K_loc;
{ // Local Newton solver
using LocalResidualVector =
Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
auto const update_residual = [&](LocalResidualVector& residual) {
calculateResidualBurgers(dt, eps_m_d_i, eps_m_d_t, sigd_j, sigd_t,
state.eps_K_j, state.eps_K_t,
state.eps_M_j, state.eps_M_t, residual,
local_lubby2_properties);
};
auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
calculateJacobianBurgers(
t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
local_lubby2_properties); // for solution dependent Jacobians
};
auto const update_solution = [&](LocalResidualVector const& increment) {
// increment solution vectors
sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
KelvinVectorSize * 0);
state.eps_K_j.noalias() +=
increment.template segment<KelvinVectorSize>(KelvinVectorSize *
1);
state.eps_M_j.noalias() +=
increment.template segment<KelvinVectorSize>(KelvinVectorSize *
2);
// Calculate effective stress and update material properties
sig_eff = MathLib::KelvinVector::Invariants<
KelvinVectorSize>::equivalentStress(sigd_j);
local_lubby2_properties.update(sig_eff);
};
auto newton_solver = NumLib::NewtonRaphson<
decltype(linear_solver), LocalJacobianMatrix,
decltype(update_jacobian), LocalResidualVector,
decltype(update_residual), decltype(update_solution)>(
linear_solver, update_jacobian, update_residual, update_solution,
_nonlinear_solver_parameters);
auto const success_iterations = newton_solver.solve(K_loc);
if (!success_iterations)
{
return {};
}
// If the Newton loop didn't run, the linear solver will not be
// initialized.
// This happens usually for the first iteration of the first timestep.
if (*success_iterations == 0)
{
linear_solver.compute(K_loc);
}
}
KelvinMatrix C =
tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
local_lubby2_properties.KM0,
linear_solver);
// Hydrostatic part for the stress and the tangent.
double const delta_eps_m_trace = Invariants::trace(eps_m - eps_m_prev);
double const sigma_trace_prev = Invariants::trace(sigma_prev);
KelvinVector const sigma =
local_lubby2_properties.GM0 * sigd_j +
(local_lubby2_properties.KM0 * delta_eps_m_trace +
sigma_trace_prev / 3.) *
Invariants::identity2;
return {std::make_tuple(
sigma,
std::unique_ptr<
typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
new MaterialStateVariables{state}},
C)};
}
template <int DisplacementDim>
void Lubby2<DisplacementDim>::calculateResidualBurgers(
const double dt,
const KelvinVector& strain_curr,
const KelvinVector& strain_t,
const KelvinVector& stress_curr,
const KelvinVector& stress_t,
const KelvinVector& strain_Kel_curr,
const KelvinVector& strain_Kel_t,
const KelvinVector& strain_Max_curr,
const KelvinVector& strain_Max_t,
ResidualVector& res,
detail::LocalLubby2Properties<DisplacementDim> const& properties) const
{
// calculate stress residual
res.template segment<KelvinVectorSize>(0).noalias() =
(stress_curr - stress_t) -
2. * ((strain_curr - strain_t) - (strain_Kel_curr - strain_Kel_t) -
(strain_Max_curr - strain_Max_t));
// calculate Kelvin strain residual
res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
(strain_Kel_curr - strain_Kel_t) -
dt / (2. * properties.etaK) *
(properties.GM0 * stress_curr -
2. * properties.GK * strain_Kel_curr);
// calculate Maxwell strain residual
res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
(strain_Max_curr - strain_Max_t) -
dt * 0.5 * properties.GM0 / properties.etaM * stress_curr;
}
template <int DisplacementDim>
void Lubby2<DisplacementDim>::calculateJacobianBurgers(
double const t,
ParameterLib::SpatialPosition const& x,
const double dt,
JacobianMatrix& Jac,
double s_eff,
const KelvinVector& sig_i,
const KelvinVector& eps_K_i,
detail::LocalLubby2Properties<DisplacementDim> const& properties) const
{
Jac.setZero();
// build G_11
Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
.diagonal()
.setConstant(1.);
// build G_12
Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
.diagonal()
.setConstant(2.);
// build G_13
Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
2 * KelvinVectorSize)
.diagonal()
.setConstant(2.);
// build G_21
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
.noalias() =
-0.5 * dt * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
if (s_eff > 0.)
{
KelvinVector const eps_K_aid =
1. / (properties.etaK * properties.etaK) *
(properties.GM0 * sig_i - 2. * properties.GK * eps_K_i);
KelvinVector const dG_K = 1.5 * _mp.mK(t, x)[0] * properties.GK *
properties.GM0 / s_eff * sig_i;
KelvinVector const dmu_vK = 1.5 * _mp.mvK(t, x)[0] * properties.GM0 *
properties.etaK / s_eff * sig_i;
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
0)
.noalias() += 0.5 * dt * eps_K_aid * dmu_vK.transpose() +
dt / properties.etaK * eps_K_i * dG_K.transpose();
}
// build G_22
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
KelvinVectorSize)
.diagonal()
.setConstant(1. + dt * properties.GK / properties.etaK);
// nothing to do for G_23
// build G_31
Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
0)
.noalias() =
-0.5 * dt * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
if (s_eff > 0.)
{
KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
properties.etaM / s_eff * sig_i;
Jac.template block<KelvinVectorSize, KelvinVectorSize>(
2 * KelvinVectorSize, 0)
.noalias() += 0.5 * dt * properties.GM0 /
(properties.etaM * properties.etaM) * sig_i *
dmu_vM.transpose();
}
// nothing to do for G_32
// build G_33
Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
2 * KelvinVectorSize)
.diagonal()
.setConstant(1.);
}
template class Lubby2<2>;
template class Lubby2<3>;
} // namespace Lubby2
} // namespace Solids
} // namespace MaterialLib