Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/src/solvers/maxwell_viscoelastic.cpp
b0c44a420f99archive/legacy_mpi

Last Change
Raw File

maxwell_viscoelastic.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides),
* Laboratory (IJLRDA - Institut Jean Le Rond d'Alembert)
* Copyright (©) 2020-2024 Lucas Frérot
*
* 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "maxwell_viscoelastic.hh"
#include "logger.hh"
#include "loop.hh"
#include "model.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
MaxwellViscoelastic::MaxwellViscoelastic(Model& model,
const GridBase<Real>& surface,
Real tolerance, Real time_step,
std::vector<Real> shear_moduli,
std::vector<Real> characteristic_times)
: PolonskyKeerRey(model, surface, tolerance, PolonskyKeerRey::pressure,
PolonskyKeerRey::pressure),
time_step_(time_step), shear_moduli_(shear_moduli),
characteristic_times_(characteristic_times),
M(characteristic_times.size()),
U(model.getDisplacement().getNbPoints(),
model.getDisplacement().getNbComponents()) {
// Check that the number of shear moduli and characteristic times is correct
if (shear_moduli_.size() != characteristic_times_.size()) {
throw std::invalid_argument(
"The number of shear moduli and characteristic times must be equal.");
}
for (UInt k = 0; k < M.size(); ++k) {
M[k] = GridBase<Real>(model.getDisplacement().getNbPoints(),
model.getDisplacement().getNbComponents());
}
}
/* ------------------------------------------------------------------------- */
Real MaxwellViscoelastic::computeGtilde(
const Real& time_step, const std::vector<Real>& shear_moduli,
const std::vector<Real>& characteristic_times) const {
Real gtilde = 0;
for (size_t i = 0; i < shear_moduli.size(); ++i) {
gtilde += shear_moduli[i] * characteristic_times[i] /
(characteristic_times[i] + time_step);
}
return gtilde;
}
/* ------------------------------------------------------------------------- */
std::vector<Real> MaxwellViscoelastic::computeGamma(
const Real& time_step, const std::vector<Real>& shear_moduli,
const std::vector<Real>& characteristic_times) const {
std::vector<Real> gamma(shear_moduli.size());
for (size_t i = 0; i < shear_moduli.size(); ++i) {
gamma[i] = characteristic_times[i] / (characteristic_times[i] + time_step);
}
return gamma;
}
/* ------------------------------------------------------------------------- */
void MaxwellViscoelastic::reset() {
for (auto& m : M) {
m = 0;
}
U = 0;
}
/* ------------------------------------------------------------------------- */
Real MaxwellViscoelastic::solve(std::vector<Real> target_viscoelastic) {
auto& M_new = model.getDisplacement();
auto gtilde = computeGtilde(time_step_, shear_moduli_, characteristic_times_);
auto gamma = computeGamma(time_step_, shear_moduli_, characteristic_times_);
// Elastic branch of the generalized maxwell model is determined from model
Real alpha = gtilde + model.getShearModulus();
Real beta = gtilde;
GridBase<Real> M_maxwell(M_new.getNbPoints(), M_new.getNbComponents());
GridBase<Real> H_old(this->surface);
for (size_t k = 0; k < M.size(); ++k) {
auto gamma_k = gamma[k];
Loop::loop(
[gamma_k](Real& M_maxwell, const Real& M) { M_maxwell += gamma_k * M; },
M_maxwell, M[k]);
}
Loop::loop(
[alpha, beta, gamma](Real& h_new, const Real& surface, const Real& u,
const Real& maxwell) {
h_new = alpha * surface - beta * u + maxwell;
},
this->surface, H_old, U, M_maxwell);
// Normalize integral operator for the solve operation
this->integral_op->field<Complex>("influence") *= model.getShearModulus();
Real error = PolonskyKeerRey::solve(target_viscoelastic);
model.solveNeumann();
this->integral_op->field<Complex>("influence") *= 1 / model.getShearModulus();
GridBase<Real> U_new(U.getNbPoints(), U.getNbComponents());
Loop::loop(
[alpha, beta](Real& U_new, const Real& M_new, const Real& M_maxwell,
const Real& U) {
U_new = (1 / alpha) * (M_new - M_maxwell + beta * U);
},
U_new, M_new, M_maxwell, U);
for (size_t k = 0; k < shear_moduli_.size(); ++k) {
auto gamma_k = gamma[k];
auto shear_moduli_k = shear_moduli_[k];
Loop::loop(
[gamma_k, shear_moduli_k](Real& M, const Real& U_new, const Real& U) {
M = gamma_k * (M + shear_moduli_k * (U_new - U));
},
M[k], U_new, U);
}
// Reset surface
this->surface = H_old;
// Save displacements
U = U_new;
// Reconstruct true displacements
*dual = U;
*dual *= model.getShearModulus();
*dual -= this->surface;
enforceAdmissibleState();
return error;
}
} // namespace tamaas
/* ------------------------------------------------------------------------- */
c4science · Help