From ccd919417fea265f64c8e32f9c08cbfd34f8c3e8 Mon Sep 17 00:00:00 2001 From: Antoine Jacquey Date: Thu, 12 Dec 2024 12:42:05 -0500 Subject: [PATCH 1/5] Implemented Blanco-Martin model for rock salt #24 --- include/materials/BVBlancoMartinModelUpdate.h | 62 ++++ src/materials/BVBlancoMartinModelUpdate.C | 293 ++++++++++++++++++ 2 files changed, 355 insertions(+) create mode 100644 include/materials/BVBlancoMartinModelUpdate.h create mode 100644 src/materials/BVBlancoMartinModelUpdate.C diff --git a/include/materials/BVBlancoMartinModelUpdate.h b/include/materials/BVBlancoMartinModelUpdate.h new file mode 100644 index 0000000..956f6ba --- /dev/null +++ b/include/materials/BVBlancoMartinModelUpdate.h @@ -0,0 +1,62 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#pragma once + +#include "BVTwoCreepUpdateBase.h" + +class BVBlancoMartinModelUpdate : public BVTwoCreepUpdateBase +{ +public: + static InputParameters validParams(); + BVBlancoMartinModelUpdate(const InputParameters & parameters); + +protected: + virtual void initQpStatefulProperties() override; + virtual ADReal creepRate(const std::vector & eqv_strain_incr, + const unsigned int i) override; + virtual ADReal creepRateLemaitre(const std::vector & eqv_strain_incr); + virtual ADReal creepRateMunsonDawson(const std::vector & eqv_strain_incr); + virtual ADReal creepRateDerivative(const std::vector & eqv_strain_incr, + const unsigned int i, + const unsigned int j) override; + virtual ADReal creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, + const unsigned int /*j*/); + virtual ADReal creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, + const unsigned int j); + virtual ADReal lemaitreCreepStrain(const std::vector & eqv_strain_incr); + virtual ADReal munsondawsonCreepStrain(const std::vector & eqv_strain_incr); + virtual void preReturnMap() override; + virtual void postReturnMap(const std::vector & eqv_strain_incr) override; + + // Lemaitre creep strain rate parameters + const Real _alpha; + const Real _kr1; + const Real _kr2; + const Real _beta1; + const Real _beta2; + + // Munson-Dawson creep strain rate parameters + const Real _A1; + const Real _n1; + const Real _A; + const Real _B; + const Real _m; + const Real _n; + + // Internal variable for Lemaitre and Munson-Dawson creep strain + ADMaterialProperty & _eqv_creep_strain_L; + const MaterialProperty & _eqv_creep_strain_L_old; + ADMaterialProperty & _eqv_creep_strain_R; + const MaterialProperty & _eqv_creep_strain_R_old; +}; diff --git a/src/materials/BVBlancoMartinModelUpdate.C b/src/materials/BVBlancoMartinModelUpdate.C new file mode 100644 index 0000000..e93bbb5 --- /dev/null +++ b/src/materials/BVBlancoMartinModelUpdate.C @@ -0,0 +1,293 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#include "BVBlancoMartinModelUpdate.h" + +registerMooseObject("BeaverApp", BVBlancoMartinModelUpdate); + +InputParameters +BVBlancoMartinModelUpdate::validParams() +{ + InputParameters params = BVTwoCreepUpdateBase::validParams(); + params.addClassDescription( + "Material for computing a RTL2020 creep update. See Azabou et al. (2021), Rock salt " + "behavior: From laboratory experiments to pertinent long-term predictions."); + params.addRequiredRangeCheckedParam("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); + params.addRequiredRangeCheckedParam("kr1", "kr1 > 0.0", "The kr1 parameter."); + params.addRequiredRangeCheckedParam("kr2", "kr2 > 0.0", "The kr2 parameter."); + params.addRequiredRangeCheckedParam("beta1", "beta1 > 0.0", "The beta1 parameter."); + params.addRequiredRangeCheckedParam("beta2", "beta2 > 0.0", "The beta2 parameter."); + params.addRequiredRangeCheckedParam("A1", "A1 > 0.0", "The A1 parameter."); + params.addRequiredRangeCheckedParam("n1", "n1 > 0.0", "The n1 parameter."); + params.addRequiredRangeCheckedParam("A", "A >= 0.0", "The A parameter."); + params.addRequiredRangeCheckedParam("B", "B >= 0.0", "The B parameter."); + params.addRequiredRangeCheckedParam("m", "m > 1.0", "The m parameter."); + params.addRequiredRangeCheckedParam("n", "n > 1.0", "The n parameter."); + return params; +} + +BVBlancoMartinModelUpdate::BVBlancoMartinModelUpdate(const InputParameters & parameters) + : BVTwoCreepUpdateBase(parameters), + // Modified Lemaitre creep strain rate parameters + _alpha(getParam("alpha")), + _kr1(getParam("kr1")), + _kr2(getParam("kr2")), + _beta1(getParam("beta1")), + _beta2(getParam("beta2")), + // Munson-Dawson creep strain rate parameters + _A1(getParam("A1")), + _n1(getParam("n1")), + _A(getParam("A")), + _B(getParam("B")), + _m(getParam("m")), + _n(getParam("n")), + // Internal variable for Lemaitre and Munson-Dawson creep strain + _eqv_creep_strain_L(declareADProperty(_base_name + "eqv_creep_strain_L")), + _eqv_creep_strain_L_old(getMaterialPropertyOld(_base_name + "eqv_creep_strain_L")), + _eqv_creep_strain_R(declareADProperty(_base_name + "eqv_creep_strain_R")), + _eqv_creep_strain_R_old(getMaterialPropertyOld(_base_name + "eqv_creep_strain_R")) +{ +} + +void +BVBlancoMartinModelUpdate::initQpStatefulProperties() +{ + _eqv_creep_strain_L[_qp] = 0.0; + _eqv_creep_strain_R[_qp] = 0.0; +} + +ADReal +BVBlancoMartinModelUpdate::creepRate(const std::vector & eqv_strain_incr, const unsigned int i) +{ + if (i == 0) // Lemaitre + return creepRateLemaitre(eqv_strain_incr); + else if (i == 1) // Munson-Dawson + return creepRateMunsonDawson(eqv_strain_incr); + else + throw MooseException("BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRate`!"); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateLemaitre(const std::vector & eqv_strain_incr) +{ + if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) + return _alpha * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, + _beta2 / _alpha)); + else + return _alpha * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, + _beta2 / _alpha)) * + std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateMunsonDawson(const std::vector & eqv_strain_incr) +{ + ADReal saturation_strain = + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + + ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + + if (gamma_ms < saturation_strain) + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, + _beta2 / _alpha)); + else + return -_B * std::pow(gamma_ms / saturation_strain, _m) * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, + _beta2 / _alpha)); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateDerivative(const std::vector & eqv_strain_incr, + const unsigned int i, + const unsigned int j) +{ + if (i == 0) // Lemaitre + return creepRateLemaitreDerivative(eqv_strain_incr, j); + else if (i == 1) // Munson-Dawson + return creepRateMunsonDawsonDerivative(eqv_strain_incr, j); + else + throw MooseException( + "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, + const unsigned int j) +{ + if (j == 0) // Lemaitre wrt Lemaitre + if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) + return -3.0 * _G * + (_beta1 / _kr1 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)); + else + return -3.0 * _G * + (_beta1 / _kr1 * + std::pow( + (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow( + (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)) + + (_alpha - 1.0) * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)); + else if (j == 1) // Lemaitre wrt Munson-Dawson + if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) + return -3.0 * _G * + (_beta1 / _kr1 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)); + else + return -3.0 * _G * + (_beta1 / _kr1 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)) * + std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + else + throw MooseException( + "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, + const unsigned int j) +{ + ADReal saturation_strain = + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + + ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + + if (j == 0) // Munson-Dawson wrt Lemaitre + if (gamma_ms < saturation_strain) + return -3.0 * _G / _alpha * _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + (_beta1 / _kr1 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)); + else + return 3.0 * _G / _alpha * _B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + (_beta1 / _kr1 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)); + else if (j == 1) // Munson-Dawson wrt Munson-Dawson + if (gamma_ms < saturation_strain) + return -3.0 * _G / _alpha * _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + (_beta1 / _kr1 * + std::pow( + (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow( + (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)) - + _A * _n / saturation_strain * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha)); + else + return 3.0 * _G / _alpha * _B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + (_beta1 / _kr1 * + std::pow( + (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha - 1.0) + + _beta2 / _kr2 * + std::pow( + (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha - 1.0)) - + _B * _m * std::pow(gamma_ms / saturation_strain - 1.0, _m - 1.0) * + (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr1, + _beta1 / _alpha) + + std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / + _kr2, + _beta2 / _alpha)); + else + throw MooseException( + "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); +} + +ADReal +BVBlancoMartinModelUpdate::lemaitreCreepStrain(const std::vector & eqv_strain_incr) +{ + return _eqv_creep_strain_L_old[_qp] + eqv_strain_incr[0]; +} + +ADReal +BVBlancoMartinModelUpdate::munsondawsonCreepStrain(const std::vector & eqv_strain_incr) +{ + return _eqv_creep_strain_R_old[_qp] + eqv_strain_incr[1]; +} + +void +BVBlancoMartinModelUpdate::preReturnMap() +{ + _eqv_creep_strain_L[_qp] = _eqv_creep_strain_L_old[_qp]; + _eqv_creep_strain_R[_qp] = _eqv_creep_strain_R_old[_qp]; +} + +void +BVBlancoMartinModelUpdate::postReturnMap(const std::vector & eqv_strain_incr) +{ + _eqv_creep_strain_L[_qp] = lemaitreCreepStrain(eqv_strain_incr); + _eqv_creep_strain_R[_qp] = munsondawsonCreepStrain(eqv_strain_incr); +} \ No newline at end of file From cb1fbcaa0c85e04528b64f710e1bec3b6dfc5b16 Mon Sep 17 00:00:00 2001 From: Antoine Jacquey Date: Thu, 19 Dec 2024 17:10:05 +0100 Subject: [PATCH 2/5] Implemented modified Lemaitre and modified RTL2020 #24 --- .../viscoelasticity/blanco-martin/RTL_S1-1.py | 31 ++ .../blanco-martin/blanco-martin-RTL-2024.i | 296 +++++++++++++++++ .../blanco-martin-lemaitre-2024.i | 300 ++++++++++++++++++ include/materials/BVBlancoMartinModelUpdate.h | 2 + .../materials/BVModifiedLemaitreModelUpdate.h | 44 +++ src/materials/BVBlancoMartinModelUpdate.C | 190 ++++------- src/materials/BVModifiedLemaitreModelUpdate.C | 117 +++++++ 7 files changed, 851 insertions(+), 129 deletions(-) create mode 100644 examples/viscoelasticity/blanco-martin/RTL_S1-1.py create mode 100644 examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i create mode 100644 examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i create mode 100644 include/materials/BVModifiedLemaitreModelUpdate.h create mode 100644 src/materials/BVModifiedLemaitreModelUpdate.C diff --git a/examples/viscoelasticity/blanco-martin/RTL_S1-1.py b/examples/viscoelasticity/blanco-martin/RTL_S1-1.py new file mode 100644 index 0000000..c374f24 --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/RTL_S1-1.py @@ -0,0 +1,31 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.style.use('publication.mplstyle') +plt.rcParams['text.usetex'] = False # Ensure LaTeX is not used for rendering text + +def read_csv(filename): + # Read data + time, EpsAx = np.loadtxt(filename, delimiter=',', skiprows=1, usecols=[0, 4], unpack=True) + return time, EpsAx + +def plot_figure(): + fig, ax = plt.subplots() + ax.set_title('Axial_strain (-) vs time (d)') # Use Unicode directly + +# Reading data + time, EpsAx = read_csv("RTL_S1-1.csv") +# Plot AE data + ax.plot(time, EpsAx, color="k", label="EpsAx") +# Additional plot settings + ax.set_xlabel('Time (days)') + ax.set_ylabel('EpsAx (-)') + # ax.set_yscale('log') # Set y-axis to log base 10 + ax.legend() +# Axes limit +# ax.set_xlim(4, 8) +# Save figure + fig.savefig("RTL_S1-1.pdf", format="pdf", dpi=300, bbox_inches="tight") + +if __name__ == '__main__': + plot_figure() diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i b/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i new file mode 100644 index 0000000..501cece --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i @@ -0,0 +1,296 @@ +# Modified Lemaitre creep model +# See Blanco-Martin et al. (2023) +# Parameters +# Units: stress in MPa, time in days, strain in m / m +E = 12000 +nu = 0.3 +alpha = 0.575 +kr1 = 1.302 +beta1 = 3.053 +kr2 = 0.091 +beta2 = 1.053 +A = 100 +n = 9 +A1 = 0.034 +n1 = 1.499 +P = 5.0 +Q = 5.0 + +[Mesh] + type = GeneratedMesh + dim = 3 + # nx = 5 + # ny = 10 + # nz = 5 + nx = 1 + ny = 1 + nz = 1 + xmin = 0 + xmax = 65e-03 + ymin = 0 + ymax = 130e-03 + zmin = 0 + zmax = 65e-03 +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] + [disp_z] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [stress_x] + type = BVStressDivergence + component = x + variable = disp_x + [] + [stress_y] + type = BVStressDivergence + component = y + variable = disp_y + [] + [stress_z] + type = BVStressDivergence + component = z + variable = disp_z + [] +[] + +[AuxVariables] + [eqv_stress] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain_rate] + order = CONSTANT + family = MONOMIAL + [] + [eqv_creep_strain_L] + order = CONSTANT + family = MONOMIAL + [] + [strain_yy] + order = CONSTANT + family = MONOMIAL + [] + [stress_yy] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [eqv_stress_aux] + type = BVMisesStressAux + variable = eqv_stress + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_aux] + type = BVEqvStrainAux + variable = eqv_strain + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_rate_aux] + type = BVEqvStrainRateAux + variable = eqv_strain_rate + execute_on = 'TIMESTEP_END' + [] + [eqv_creep_strain_L_aux] + type = ADMaterialRealAux + variable = eqv_creep_strain_L + property = eqv_creep_strain_L + execute_on = 'TIMESTEP_END' + [] + [strain_zz_aux] + type = BVStrainComponentAux + variable = strain_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] + [stress_zz_aux] + type = BVStressComponentAux + variable = stress_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] +[] + +[Functions] + [loading2] + type = ParsedFunction + expression = 'if(t<=20,5.5,if(t<=55,6,if(t<=76,7,if(t<=97,10,if(t<=118,15,if(t<=139,20,25))))))' + [] +[] + +[BCs] + [no_x] + type = DirichletBC + variable = disp_x + boundary = 'left' + value = 0.0 + [] + [no_y] + type = DirichletBC + variable = disp_y + boundary = 'bottom' + value = 0.0 + [] + [no_z] + type = DirichletBC + variable = disp_z + boundary = 'back' + value = 0.0 + [] + [BVPressure] + [pressure_right] + boundary = 'right' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_front] + boundary = 'front' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_top] + boundary = 'top' + displacement_vars = 'disp_x disp_y disp_z' + function = loading2 + [] + [] +[] + +[Materials] + [elasticity] + type = BVMechanicalMaterial + displacements = 'disp_x disp_y disp_z' + young_modulus = ${E} + poisson_ratio = ${nu} + initial_stress = '-${P} -${Q} -${P}' + inelastic_models = 'viscoelastic' + [] + [viscoelastic] + type = BVBlancoMartinModelUpdate + alpha = ${alpha} + kr1 = ${kr1} + beta1 = ${beta1} + kr2 = ${kr2} + beta2 = ${beta2} + A1 = ${A1} + n1 = ${n1} + A = ${A} + n = ${n} + B = 0.0 + m = ${n} + [] +[] + +[Postprocessors] + [eqv_strain_rate] + type = ElementAverageValue + variable = eqv_strain_rate + outputs = csv + [] + [eqv_strain] + type = ElementAverageValue + variable = eqv_strain + outputs = csv + [] + [eqv_strain_L] + type = ElementAverageValue + variable = eqv_creep_strain_L + outputs = csv + [] + [strain_zz] + type = ElementAverageValue + variable = strain_yy + outputs = csv + [] + [stress_zz] + type = ElementAverageValue + variable = stress_yy + outputs = csv + [] + [q] + type = ElementAverageValue + variable = eqv_stress + outputs = csv + [] +[] + +[Preconditioning] + active = 'hypre' + [hypre] + type = SMP + full = true + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type -pc_hypre_type + -snes_atol -snes_rtol -snes_stol -snes_max_it + -snes_linesearch_type' + petsc_options_value = 'hypre boomeramg + 1.0e-10 1.0e-12 0 20 + basic' + [] + [superlu] + type = SMP + full = true + petsc_options = '-snes_ksp_ew -snes_converged_reason -ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_diagonal_scale -ksp_diagonal_scale_fix' + petsc_options_iname = '-snes_type + -snes_atol -snes_rtol -snes_max_it + -pc_type -pc_factor_mat_solver_package + -snes_linesearch_type' + petsc_options_value = 'newtonls + 1e-10 1e-12 50 + lu superlu_dist + l2' + [] + [asm] + type = SMP + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-ksp_type + -pc_type + -sub_pc_type + -snes_type -snes_atol -snes_rtol -snes_max_it -snes_linesearch_type + -ksp_gmres_restart' + petsc_options_value = 'fgmres + asm + ilu + newtonls 1e-10 1e-10 120 basic + 201' + [] +[] + +[Executioner] + type = Transient + solve_type = 'NEWTON' + automatic_scaling = true + start_time = 0.0 + end_time = 200 # 200 days + dt = 0.02 + timestep_tolerance = 1.0e-10 +[] + +[Outputs] + perf_graph = true + exodus = true + [csv] + type = CSV + execute_on = 'timestep_end' + [] +[] \ No newline at end of file diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i new file mode 100644 index 0000000..ef9cefa --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i @@ -0,0 +1,300 @@ +# Modified Lemaitre creep model +# See Blanco-Martin et al. (2024) +# Parameters +# Units: stress in MPa, time in days, strain in m / m +E = 12000 +nu = 0.3 +alpha = 0.326 +kr1 = 0.7 +beta1 = 2.922 +kr2 = 0.009 +beta2 = 0.867 +P = 5.0 +Q = 5.0 + +[Mesh] + type = GeneratedMesh + dim = 3 + # nx = 5 + # ny = 10 + # nz = 5 + nx = 1 + ny = 1 + nz = 1 + xmin = 0 + xmax = 65e-03 + ymin = 0 + ymax = 130e-03 + zmin = 0 + zmax = 65e-03 +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] + [disp_z] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [stress_x] + type = BVStressDivergence + component = x + variable = disp_x + [] + [stress_y] + type = BVStressDivergence + component = y + variable = disp_y + [] + [stress_z] + type = BVStressDivergence + component = z + variable = disp_z + [] +[] + +[AuxVariables] + [eqv_stress] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain_rate] + order = CONSTANT + family = MONOMIAL + [] + [eqv_creep_strain_L] + order = CONSTANT + family = MONOMIAL + [] + [strain_yy] + order = CONSTANT + family = MONOMIAL + [] + [stress_yy] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [eqv_stress_aux] + type = BVMisesStressAux + variable = eqv_stress + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_aux] + type = BVEqvStrainAux + variable = eqv_strain + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_rate_aux] + type = BVEqvStrainRateAux + variable = eqv_strain_rate + execute_on = 'TIMESTEP_END' + [] + [eqv_creep_strain_L_aux] + type = ADMaterialRealAux + variable = eqv_creep_strain_L + property = eqv_creep_strain + execute_on = 'TIMESTEP_END' + [] + [strain_zz_aux] + type = BVStrainComponentAux + variable = strain_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] + [stress_zz_aux] + type = BVStressComponentAux + variable = stress_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] +[] + +[Functions] + [loading2] + type = ParsedFunction + expression = 'if(t<=20,5.5,if(t<=55,6,if(t<=76,7,if(t<=97,10,if(t<=118,15,if(t<=139,20,25))))))' + [] +[] + +[BCs] + [no_x] + type = DirichletBC + variable = disp_x + boundary = 'left' + value = 0.0 + [] + [no_y] + type = DirichletBC + variable = disp_y + boundary = 'bottom' + value = 0.0 + [] + [no_z] + type = DirichletBC + variable = disp_z + boundary = 'back' + value = 0.0 + [] + [BVPressure] + [pressure_right] + boundary = 'right' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_front] + boundary = 'front' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_top] + boundary = 'top' + displacement_vars = 'disp_x disp_y disp_z' + function = loading2 + [] + [] +[] + +[Materials] + [elasticity] + type = BVMechanicalMaterial + displacements = 'disp_x disp_y disp_z' + young_modulus = ${E} + poisson_ratio = ${nu} + initial_stress = '-${P} -${Q} -${P}' + inelastic_models = 'viscoelastic' + [] + # [viscoelastic] + # type = BVBlancoMartinModelUpdate + # alpha = ${alpha} + # kr1 = ${kr1} + # beta1 = ${beta1} + # kr2 = ${kr2} + # beta2 = ${beta2} + # A1 = ${A1} + # n1 = ${n1} + # A = ${A} + # n = ${n} + # B = 0.0 + # m = ${n} + # [] + [viscoelastic] + type = BVModifiedLemaitreModelUpdate + alpha = ${alpha} + kr1 = ${kr1} + beta1 = ${beta1} + kr2 = ${kr2} + beta2 = ${beta2} + [] +[] + +[Postprocessors] + [eqv_strain_rate] + type = ElementAverageValue + variable = eqv_strain_rate + outputs = csv + [] + [eqv_strain] + type = ElementAverageValue + variable = eqv_strain + outputs = csv + [] + [eqv_strain_L] + type = ElementAverageValue + variable = eqv_creep_strain_L + outputs = csv + [] + [strain_zz] + type = ElementAverageValue + variable = strain_yy + outputs = csv + [] + [stress_zz] + type = ElementAverageValue + variable = stress_yy + outputs = csv + [] + [q] + type = ElementAverageValue + variable = eqv_stress + outputs = csv + [] +[] + +[Preconditioning] + active = 'hypre' + [hypre] + type = SMP + full = true + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type -pc_hypre_type + -snes_atol -snes_rtol -snes_stol -snes_max_it + -snes_linesearch_type' + petsc_options_value = 'hypre boomeramg + 1.0e-10 1.0e-12 0 20 + basic' + [] + [superlu] + type = SMP + full = true + petsc_options = '-snes_ksp_ew -snes_converged_reason -ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_diagonal_scale -ksp_diagonal_scale_fix' + petsc_options_iname = '-snes_type + -snes_atol -snes_rtol -snes_max_it + -pc_type -pc_factor_mat_solver_package + -snes_linesearch_type' + petsc_options_value = 'newtonls + 1e-10 1e-12 50 + lu superlu_dist + l2' + [] + [asm] + type = SMP + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-ksp_type + -pc_type + -sub_pc_type + -snes_type -snes_atol -snes_rtol -snes_max_it -snes_linesearch_type + -ksp_gmres_restart' + petsc_options_value = 'fgmres + asm + ilu + newtonls 1e-10 1e-10 120 basic + 201' + [] +[] + +[Executioner] + type = Transient + solve_type = 'NEWTON' + automatic_scaling = true + start_time = 0.0 + end_time = 200 # 200 days + dt = 0.02 + timestep_tolerance = 1.0e-10 +[] + +[Outputs] + perf_graph = true + exodus = true + [csv] + type = CSV + execute_on = 'timestep_end' + [] +[] \ No newline at end of file diff --git a/include/materials/BVBlancoMartinModelUpdate.h b/include/materials/BVBlancoMartinModelUpdate.h index 956f6ba..6d5b8f4 100644 --- a/include/materials/BVBlancoMartinModelUpdate.h +++ b/include/materials/BVBlancoMartinModelUpdate.h @@ -25,11 +25,13 @@ class BVBlancoMartinModelUpdate : public BVTwoCreepUpdateBase virtual void initQpStatefulProperties() override; virtual ADReal creepRate(const std::vector & eqv_strain_incr, const unsigned int i) override; + virtual ADReal creepRateR(const std::vector & eqv_strain_incr); virtual ADReal creepRateLemaitre(const std::vector & eqv_strain_incr); virtual ADReal creepRateMunsonDawson(const std::vector & eqv_strain_incr); virtual ADReal creepRateDerivative(const std::vector & eqv_strain_incr, const unsigned int i, const unsigned int j) override; + virtual ADReal creepRateRDerivative(const std::vector & eqv_strain_incr); virtual ADReal creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, const unsigned int /*j*/); virtual ADReal creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, diff --git a/include/materials/BVModifiedLemaitreModelUpdate.h b/include/materials/BVModifiedLemaitreModelUpdate.h new file mode 100644 index 0000000..148d340 --- /dev/null +++ b/include/materials/BVModifiedLemaitreModelUpdate.h @@ -0,0 +1,44 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#pragma once + +#include "BVCreepUpdateBase.h" + +class BVModifiedLemaitreModelUpdate : public BVCreepUpdateBase +{ +public: + static InputParameters validParams(); + BVModifiedLemaitreModelUpdate(const InputParameters & parameters); + +protected: + virtual void initQpStatefulProperties() override; + virtual ADReal creepRateR(const ADReal & eqv_strain_incr); + virtual ADReal creepRate(const ADReal & eqv_strain_incr) override; + virtual ADReal creepRateRDerivative(const ADReal & eqv_strain_incr); + virtual ADReal creepRateDerivative(const ADReal & eqv_strain_incr) override; + virtual ADReal lemaitreCreepStrain(const ADReal & eqv_strain_incr); + virtual void preReturnMap() override; + virtual void postReturnMap(const ADReal & eqv_strain_incr) override; + + // Lemaitre creep strain rate parameters + const Real _alpha; + const Real _kr1; + const Real _kr2; + const Real _beta1; + const Real _beta2; + + // Internal variable for creep strain + ADMaterialProperty & _eqv_creep_strain; + const MaterialProperty & _eqv_creep_strain_old; +}; diff --git a/src/materials/BVBlancoMartinModelUpdate.C b/src/materials/BVBlancoMartinModelUpdate.C index e93bbb5..7c41bd8 100644 --- a/src/materials/BVBlancoMartinModelUpdate.C +++ b/src/materials/BVBlancoMartinModelUpdate.C @@ -77,44 +77,42 @@ BVBlancoMartinModelUpdate::creepRate(const std::vector & eqv_strain_incr throw MooseException("BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRate`!"); } +ADReal +BVBlancoMartinModelUpdate::creepRateR(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 0.0; + else + return 1.0e-06 * std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha); +} + ADReal BVBlancoMartinModelUpdate::creepRateLemaitre(const std::vector & eqv_strain_incr) { - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return _alpha * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, - _beta2 / _alpha)); + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateR(eqv_strain_incr); else - return _alpha * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, - _beta2 / _alpha)) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + return _alpha * creepRateR(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); } ADReal BVBlancoMartinModelUpdate::creepRateMunsonDawson(const std::vector & eqv_strain_incr) { - ADReal saturation_strain = - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; - ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); if (gamma_ms < saturation_strain) return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, - _beta2 / _alpha)); + creepRateR(eqv_strain_incr); else - return -_B * std::pow(gamma_ms / saturation_strain, _m) * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _kr2, - _beta2 / _alpha)); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateR(eqv_strain_incr); } ADReal @@ -131,62 +129,39 @@ BVBlancoMartinModelUpdate::creepRateDerivative(const std::vector & eqv_s "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); } +ADReal +BVBlancoMartinModelUpdate::creepRateRDerivative(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 1.0; + else + return -1.0e-06 * 3.0 * _G / _alpha * + std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha - 1.0) * + (_beta1 / _kr1 * std::pow(q / _kr1, _beta1 - 1.0) + + _beta2 / _kr2 * std::pow(q / _kr2, _beta2 - 1.0)); +} + ADReal BVBlancoMartinModelUpdate::creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, const unsigned int j) { + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + if (j == 0) // Lemaitre wrt Lemaitre - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return -3.0 * _G * - (_beta1 / _kr1 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)); + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); else - return -3.0 * _G * - (_beta1 / _kr1 * - std::pow( - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow( - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)) + - (_alpha - 1.0) * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)); + return std::pow(gamma_l, -1.0 / _alpha) * + (_alpha * gamma_l * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * (_alpha - 1.0) * creepRateR(eqv_strain_incr)); + else if (j == 1) // Lemaitre wrt Munson-Dawson - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return -3.0 * _G * - (_beta1 / _kr1 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)); + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); else - return -3.0 * _G * - (_beta1 / _kr1 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + return _alpha * creepRateRDerivative(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); else throw MooseException( "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); @@ -196,71 +171,28 @@ ADReal BVBlancoMartinModelUpdate::creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, const unsigned int j) { - ADReal saturation_strain = - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; - ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); if (j == 0) // Munson-Dawson wrt Lemaitre if (gamma_ms < saturation_strain) - return -3.0 * _G / _alpha * _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - (_beta1 / _kr1 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)); + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + creepRateRDerivative(eqv_strain_incr); else - return 3.0 * _G / _alpha * _B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * - (_beta1 / _kr1 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateRDerivative(eqv_strain_incr); + else if (j == 1) // Munson-Dawson wrt Munson-Dawson if (gamma_ms < saturation_strain) - return -3.0 * _G / _alpha * _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - (_beta1 / _kr1 * - std::pow( - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow( - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)) - - _A * _n / saturation_strain * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha)); + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * + ((1.0 - gamma_ms / saturation_strain) * creepRateRDerivative(eqv_strain_incr) - + 1.0e+06 * _n / saturation_strain * creepRateR(eqv_strain_incr)); else - return 3.0 * _G / _alpha * _B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * - (_beta1 / _kr1 * - std::pow( - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha - 1.0) + - _beta2 / _kr2 * - std::pow( - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha - 1.0)) - - _B * _m * std::pow(gamma_ms / saturation_strain - 1.0, _m - 1.0) * - (std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr1, - _beta1 / _alpha) + - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / - _kr2, - _beta2 / _alpha)); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m - 1.0) * + ((gamma_ms / saturation_strain - 1.0) * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * _m / saturation_strain * creepRateR(eqv_strain_incr)); else throw MooseException( "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); diff --git a/src/materials/BVModifiedLemaitreModelUpdate.C b/src/materials/BVModifiedLemaitreModelUpdate.C new file mode 100644 index 0000000..56562f3 --- /dev/null +++ b/src/materials/BVModifiedLemaitreModelUpdate.C @@ -0,0 +1,117 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#include "BVModifiedLemaitreModelUpdate.h" + +registerMooseObject("BeaverApp", BVModifiedLemaitreModelUpdate); + +InputParameters +BVModifiedLemaitreModelUpdate::validParams() +{ + InputParameters params = BVCreepUpdateBase::validParams(); + params.addClassDescription( + "Material for computing a modified Lemaitre creep update (see Blanco-Martin 2024)."); + params.addRequiredRangeCheckedParam("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); + params.addRequiredRangeCheckedParam("kr1", "kr1 > 0.0", "The kr1 parameter."); + params.addRequiredRangeCheckedParam("kr2", "kr2 > 0.0", "The kr2 parameter."); + params.addRequiredRangeCheckedParam("beta1", "beta1 > 0.0", "The beta1 parameter."); + params.addRequiredRangeCheckedParam("beta2", "beta2 > 0.0", "The beta2 parameter."); + return params; +} + +BVModifiedLemaitreModelUpdate::BVModifiedLemaitreModelUpdate(const InputParameters & parameters) + : BVCreepUpdateBase(parameters), + // Modified Lemaitre creep strain rate parameters + _alpha(getParam("alpha")), + _kr1(getParam("kr1")), + _kr2(getParam("kr2")), + _beta1(getParam("beta1")), + _beta2(getParam("beta2")), + // Internal variable for creep strain + _eqv_creep_strain(declareADProperty(_base_name + "eqv_creep_strain")), + _eqv_creep_strain_old(getMaterialPropertyOld(_base_name + "eqv_creep_strain")) +{ +} + +void +BVModifiedLemaitreModelUpdate::initQpStatefulProperties() +{ + _eqv_creep_strain[_qp] = 0.0; +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRateR(const ADReal & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * eqv_strain_incr; + + if (q == 0.0) + return 0.0; + else + return 1.0e-06 * std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha); +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRate(const ADReal & eqv_strain_incr) +{ + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateR(eqv_strain_incr); + else + return _alpha * creepRateR(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRateRDerivative(const ADReal & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * eqv_strain_incr; + + if (q == 0.0) + return 1.0; + else + return -1.0e-06 * 3.0 * _G / _alpha * + std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha - 1.0) * + (_beta1 / _kr1 * std::pow(q / _kr1, _beta1 - 1.0) + + _beta2 / _kr2 * std::pow(q / _kr2, _beta2 - 1.0)); +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRateDerivative(const ADReal & eqv_strain_incr) +{ + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); + else + return std::pow(gamma_l, -1.0 / _alpha) * + (_alpha * gamma_l * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * (_alpha - 1.0) * creepRateR(eqv_strain_incr)); +} + +ADReal +BVModifiedLemaitreModelUpdate::lemaitreCreepStrain(const ADReal & eqv_strain_incr) +{ + return _eqv_creep_strain_old[_qp] + eqv_strain_incr; +} + +void +BVModifiedLemaitreModelUpdate::preReturnMap() +{ + _eqv_creep_strain[_qp] = _eqv_creep_strain_old[_qp]; +} + +void +BVModifiedLemaitreModelUpdate::postReturnMap(const ADReal & eqv_strain_incr) +{ + _eqv_creep_strain[_qp] = lemaitreCreepStrain(eqv_strain_incr); +} \ No newline at end of file From 2bf43c7febc8e6fa7bd88050f4df98e39f8fde93 Mon Sep 17 00:00:00 2001 From: Antoine Jacquey Date: Thu, 19 Dec 2024 17:13:19 +0100 Subject: [PATCH 3/5] Update CI --- .github/workflows/CI.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e8f1e0b..459f0fb 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,7 @@ jobs: - name: Install MOOSE packages shell: bash -l {0} run: | - mamba install moose-dev=2024.10.17=mpich + mamba install moose-dev=2024.12.02=mpich - name: Compile BEAVER shell: bash -l {0} run: | @@ -63,7 +63,7 @@ jobs: - name: Install MOOSE packages shell: bash -l {0} run: | - mamba install moose-dev=2024.10.17=mpich lcov + mamba install moose-dev=2024.12.02=mpich lcov - name: Compile and test BEAVER (with coverage) shell: bash -l {0} run: | From 7e9cc9cd46c3586e4f073654b9a2cc685c800cad Mon Sep 17 00:00:00 2001 From: Antoine Jacquey Date: Mon, 13 Jan 2025 21:59:25 -0500 Subject: [PATCH 4/5] Blanco-Martin model finalized --- .../blanco-martin/blanco-martin-RTL-2024.py | 41 +++++++ .../blanco-martin-lemaitre-2024.py | 41 +++++++ include/materials/BVRTL2020ModelUpdate.h | 2 + src/materials/BVBlancoMartinModelUpdate.C | 4 + src/materials/BVRTL2020ModelUpdate.C | 115 ++++++++++-------- 5 files changed, 149 insertions(+), 54 deletions(-) create mode 100644 examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py create mode 100644 examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.py diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py b/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py new file mode 100644 index 0000000..bf4569d --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py @@ -0,0 +1,41 @@ +import os, csv +import numpy as np +import matplotlib.pyplot as plt +plt.style.use('../../publication.mplstyle') + +def plot_experimental_data(ax): + # Read data from file + t, stress, strain = np.loadtxt("S1-1.txt", skiprows=7, usecols=[0, 1, 3], unpack=True) + + # ax.plot(t, stress, color="k") + ax[0].plot(t, strain / 1.0e+04, color="k") + ax[1].plot(t, stress, color="k", label=r"Experimental data") + +def plot_simulation_data(ax): + # Read data from file + t, stress, strain = np.loadtxt("blanco-martin-RTL-2024_csv.csv", delimiter=',', skiprows=1, usecols=[0, 6, 5], unpack=True) + + # ax.plot(t, -stress, color="xkcd:red") + ax[0].plot(t, -strain * 1.0e+02, color="xkcd:red") + ax[1].plot(t, -stress, color="xkcd:red", label=r"Numerical fit") + +if __name__ == "__main__": + + # Figure stress + fig, axes = plt.subplots(2, 1, figsize=(3.25, 5), sharex=True) + # plt.subplots_adjust(hspace=0.12) + + plot_experimental_data(axes) + plot_simulation_data(axes) + + axes[1].set_xlim(0, 200) + axes[0].set_ylim(0, 8) + axes[1].set_ylim(0, 30) + axes[1].set_xlabel(r"Time, (days)") + axes[0].set_ylabel(r"Axial strain, $\varepsilon_{\text{ax}}$ (\%)") + axes[1].set_ylabel(r"Axial stress, $\sigma_{\text{ax}}$ (MPa)") + axes[0].set_title(r"Modified RTL2020 model (Blanco-Martìn et al., 2024)") + axes[1].legend(loc="lower right") + + # plt.show() + fig.savefig("blanco-martin-RTL-2024.png", format="PNG", dpi=300, bbox_inches="tight") diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.py b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.py new file mode 100644 index 0000000..8d447f8 --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.py @@ -0,0 +1,41 @@ +import os, csv +import numpy as np +import matplotlib.pyplot as plt +plt.style.use('../../publication.mplstyle') + +def plot_experimental_data(ax): + # Read data from file + t, stress, strain = np.loadtxt("S1-1.txt", skiprows=7, usecols=[0, 1, 3], unpack=True) + + # ax.plot(t, stress, color="k") + ax[0].plot(t, strain / 1.0e+04, color="k") + ax[1].plot(t, stress, color="k", label=r"Experimental data") + +def plot_simulation_data(ax): + # Read data from file + t, stress, strain = np.loadtxt("blanco-martin-lemaitre-2024_csv.csv", delimiter=',', skiprows=1, usecols=[0, 6, 5], unpack=True) + + # ax.plot(t, -stress, color="xkcd:red") + ax[0].plot(t, -strain * 1.0e+02, color="xkcd:red") + ax[1].plot(t, -stress, color="xkcd:red", label=r"Numerical fit") + +if __name__ == "__main__": + + # Figure stress + fig, axes = plt.subplots(2, 1, figsize=(3.25, 5), sharex=True) + # plt.subplots_adjust(hspace=0.12) + + plot_experimental_data(axes) + plot_simulation_data(axes) + + axes[1].set_xlim(0, 200) + axes[0].set_ylim(0, 8) + axes[1].set_ylim(0, 30) + axes[1].set_xlabel(r"Time, (days)") + axes[0].set_ylabel(r"Axial strain, $\varepsilon_{\text{ax}}$ (\%)") + axes[1].set_ylabel(r"Axial stress, $\sigma_{\text{ax}}$ (MPa)") + axes[0].set_title(r"Modified Lemaitre model (Blanco-Martìn et al., 2024)") + axes[1].legend(loc="lower right") + + # plt.show() + fig.savefig("blanco-martin-lemaitre-2024.png", format="PNG", dpi=300, bbox_inches="tight") diff --git a/include/materials/BVRTL2020ModelUpdate.h b/include/materials/BVRTL2020ModelUpdate.h index fba7bbb..93f134d 100644 --- a/include/materials/BVRTL2020ModelUpdate.h +++ b/include/materials/BVRTL2020ModelUpdate.h @@ -25,11 +25,13 @@ class BVRTL2020ModelUpdate : public BVTwoCreepUpdateBase virtual void initQpStatefulProperties() override; virtual ADReal creepRate(const std::vector & eqv_strain_incr, const unsigned int i) override; + virtual ADReal creepRateR(const std::vector & eqv_strain_incr); virtual ADReal creepRateLemaitre(const std::vector & eqv_strain_incr); virtual ADReal creepRateMunsonDawson(const std::vector & eqv_strain_incr); virtual ADReal creepRateDerivative(const std::vector & eqv_strain_incr, const unsigned int i, const unsigned int j) override; + virtual ADReal creepRateRDerivative(const std::vector & eqv_strain_incr); virtual ADReal creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, const unsigned int /*j*/); virtual ADReal creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, diff --git a/src/materials/BVBlancoMartinModelUpdate.C b/src/materials/BVBlancoMartinModelUpdate.C index 7c41bd8..6fe41e7 100644 --- a/src/materials/BVBlancoMartinModelUpdate.C +++ b/src/materials/BVBlancoMartinModelUpdate.C @@ -22,11 +22,13 @@ BVBlancoMartinModelUpdate::validParams() params.addClassDescription( "Material for computing a RTL2020 creep update. See Azabou et al. (2021), Rock salt " "behavior: From laboratory experiments to pertinent long-term predictions."); + // Modified Lemaitre creep strain rate parameters params.addRequiredRangeCheckedParam("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); params.addRequiredRangeCheckedParam("kr1", "kr1 > 0.0", "The kr1 parameter."); params.addRequiredRangeCheckedParam("kr2", "kr2 > 0.0", "The kr2 parameter."); params.addRequiredRangeCheckedParam("beta1", "beta1 > 0.0", "The beta1 parameter."); params.addRequiredRangeCheckedParam("beta2", "beta2 > 0.0", "The beta2 parameter."); + // Munson-Dawson creep strain rate parameters params.addRequiredRangeCheckedParam("A1", "A1 > 0.0", "The A1 parameter."); params.addRequiredRangeCheckedParam("n1", "n1 > 0.0", "The n1 parameter."); params.addRequiredRangeCheckedParam("A", "A >= 0.0", "The A parameter."); @@ -162,6 +164,7 @@ BVBlancoMartinModelUpdate::creepRateLemaitreDerivative(const std::vector return _alpha * creepRateRDerivative(eqv_strain_incr); else return _alpha * creepRateRDerivative(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); + else throw MooseException( "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); @@ -193,6 +196,7 @@ BVBlancoMartinModelUpdate::creepRateMunsonDawsonDerivative(const std::vector("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); params.addRequiredRangeCheckedParam("A2", "A2 > 0.0", "The A2 parameter."); params.addRequiredRangeCheckedParam("n2", "n2 > 0.0", "The n2 parameter."); + // Munson-Dawson creep strain rate parameters params.addRequiredRangeCheckedParam("A1", "A1 > 0.0", "The A1 parameter."); params.addRequiredRangeCheckedParam("n1", "n1 > 0.0", "The n1 parameter."); params.addRequiredRangeCheckedParam("A", "A >= 0.0", "The A parameter."); @@ -73,36 +75,42 @@ BVRTL2020ModelUpdate::creepRate(const std::vector & eqv_strain_incr, con throw MooseException("BVRTL2020ModelUpdate: error, unknow creep model called in `creepRate`!"); } +ADReal +BVRTL2020ModelUpdate::creepRateR(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 0.0; + else + return 1.0e-06 * std::pow(q / _A2, _n2); +} + ADReal BVRTL2020ModelUpdate::creepRateLemaitre(const std::vector & eqv_strain_incr) { - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2); + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateR(eqv_strain_incr); else - return _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + return _alpha * creepRateR(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); } ADReal BVRTL2020ModelUpdate::creepRateMunsonDawson(const std::vector & eqv_strain_incr) { - ADReal saturation_strain = - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; - ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); if (gamma_ms < saturation_strain) return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2); + creepRateR(eqv_strain_incr); else - return -_B * std::pow(gamma_ms / saturation_strain, _m) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateR(eqv_strain_incr); } ADReal @@ -119,33 +127,37 @@ BVRTL2020ModelUpdate::creepRateDerivative(const std::vector & eqv_strain "BVRTL2020ModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); } +ADReal +BVRTL2020ModelUpdate::creepRateRDerivative(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 1.0; + else + return -1.0e-06 * 3.0 * _G * _n2 / _A2 * std::pow(q / _A2, _n2 - 1.0); +} + ADReal BVRTL2020ModelUpdate::creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, const unsigned int j) { + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + if (j == 0) // Lemaitre wrt Lemaitre - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return -3.0 * _G / _A2 * _n2 * _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); else - return _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), -1.0 / _alpha) * - (-3.0 * _G / _A2 * _n2 * lemaitreCreepStrain(eqv_strain_incr) + - (1.0 - 1.0 / _alpha) * - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2); + return std::pow(gamma_l, -1.0 / _alpha) * + (_alpha * gamma_l * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * (_alpha - 1.0) * creepRateR(eqv_strain_incr)); + else if (j == 1) // Lemaitre wrt Munson-Dawson - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return -3.0 * _G / _A2 * _n2 * _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); else - return -3.0 * _G / _A2 * _n2 * _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + return _alpha * creepRateRDerivative(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); + else throw MooseException( "BVRTL2020ModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); @@ -155,34 +167,29 @@ ADReal BVRTL2020ModelUpdate::creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, const unsigned int j) { - ADReal saturation_strain = - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; - ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); if (j == 0) // Munson-Dawson wrt Lemaitre if (gamma_ms < saturation_strain) - return -3.0 * _G / _A2 * _n2 * _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + creepRateRDerivative(eqv_strain_incr); else - return 3.0 * _G / _A2 * _n2 * _B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateRDerivative(eqv_strain_incr); + else if (j == 1) // Munson-Dawson wrt Munson-Dawson if (gamma_ms < saturation_strain) - return -_A * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - (_n / saturation_strain * - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2 + - 3.0 * _G / _A2 * _n2 * (1.0 - gamma_ms / saturation_strain)); + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * + ((1.0 - gamma_ms / saturation_strain) * creepRateRDerivative(eqv_strain_incr) - + 1.0e+06 * _n / saturation_strain * creepRateR(eqv_strain_incr)); else return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m - 1.0) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - (_m * (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2 - - 3.0 * _G / _A2 * _n2 * (gamma_ms / saturation_strain - 1.0)); + ((gamma_ms / saturation_strain - 1.0) * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * _m / saturation_strain * creepRateR(eqv_strain_incr)); + else throw MooseException( "BVRTL2020ModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); From dfb490c4534a2319bba9e74cbcee2e4694921506 Mon Sep 17 00:00:00 2001 From: Antoine Jacquey Date: Tue, 14 Jan 2025 12:33:18 -0500 Subject: [PATCH 5/5] Tests and CI --- .github/workflows/CI.yml | 4 +- .../viscoelasticity/blanco-martin-lemaitre.i | 257 ++++++++++++++++++ .../tests/viscoelasticity/blanco-martin-rtl.i | 253 +++++++++++++++++ .../gold/blanco-martin-lemaitre_out.e | Bin 0 -> 106788 bytes .../gold/blanco-martin-rtl_out.e | Bin 0 -> 107268 bytes test/tests/viscoelasticity/tests | 10 + 6 files changed, 522 insertions(+), 2 deletions(-) create mode 100644 test/tests/viscoelasticity/blanco-martin-lemaitre.i create mode 100644 test/tests/viscoelasticity/blanco-martin-rtl.i create mode 100644 test/tests/viscoelasticity/gold/blanco-martin-lemaitre_out.e create mode 100644 test/tests/viscoelasticity/gold/blanco-martin-rtl_out.e diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 459f0fb..5fcc64a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,7 @@ jobs: - name: Install MOOSE packages shell: bash -l {0} run: | - mamba install moose-dev=2024.12.02=mpich + mamba install moose-dev=2024.12.23=mpich - name: Compile BEAVER shell: bash -l {0} run: | @@ -63,7 +63,7 @@ jobs: - name: Install MOOSE packages shell: bash -l {0} run: | - mamba install moose-dev=2024.12.02=mpich lcov + mamba install moose-dev=2024.12.23=mpich lcov - name: Compile and test BEAVER (with coverage) shell: bash -l {0} run: | diff --git a/test/tests/viscoelasticity/blanco-martin-lemaitre.i b/test/tests/viscoelasticity/blanco-martin-lemaitre.i new file mode 100644 index 0000000..52ce1cd --- /dev/null +++ b/test/tests/viscoelasticity/blanco-martin-lemaitre.i @@ -0,0 +1,257 @@ +# Modified Lemaitre creep model +# See Blanco-Martin et al. (2024) +# Parameters +# Units: stress in MPa, time in days, strain in m / m +E = 12000 +nu = 0.3 +alpha = 0.326 +kr1 = 0.7 +beta1 = 2.922 +kr2 = 0.009 +beta2 = 0.867 +P = 5.0 +Q = 5.0 +Q1 = 5.5 + +[Mesh] + type = GeneratedMesh + dim = 3 + # nx = 5 + # ny = 10 + # nz = 5 + nx = 1 + ny = 1 + nz = 1 + xmin = 0 + xmax = 65e-03 + ymin = 0 + ymax = 130e-03 + zmin = 0 + zmax = 65e-03 +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] + [disp_z] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [stress_x] + type = BVStressDivergence + component = x + variable = disp_x + [] + [stress_y] + type = BVStressDivergence + component = y + variable = disp_y + [] + [stress_z] + type = BVStressDivergence + component = z + variable = disp_z + [] +[] + +[AuxVariables] + [eqv_stress] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain_rate] + order = CONSTANT + family = MONOMIAL + [] + [eqv_creep_strain_L] + order = CONSTANT + family = MONOMIAL + [] + [strain_yy] + order = CONSTANT + family = MONOMIAL + [] + [stress_yy] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [eqv_stress_aux] + type = BVMisesStressAux + variable = eqv_stress + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_aux] + type = BVEqvStrainAux + variable = eqv_strain + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_rate_aux] + type = BVEqvStrainRateAux + variable = eqv_strain_rate + execute_on = 'TIMESTEP_END' + [] + [eqv_creep_strain_L_aux] + type = ADMaterialRealAux + variable = eqv_creep_strain_L + property = eqv_creep_strain + execute_on = 'TIMESTEP_END' + [] + [strain_zz_aux] + type = BVStrainComponentAux + variable = strain_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] + [stress_zz_aux] + type = BVStressComponentAux + variable = stress_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] +[] + +[BCs] + [no_x] + type = DirichletBC + variable = disp_x + boundary = 'left' + value = 0.0 + [] + [no_y] + type = DirichletBC + variable = disp_y + boundary = 'bottom' + value = 0.0 + [] + [no_z] + type = DirichletBC + variable = disp_z + boundary = 'back' + value = 0.0 + [] + [BVPressure] + [pressure_right] + boundary = 'right' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_front] + boundary = 'front' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_top] + boundary = 'top' + displacement_vars = 'disp_x disp_y disp_z' + function = ${Q1} + [] + [] +[] + +[Materials] + [elasticity] + type = BVMechanicalMaterial + displacements = 'disp_x disp_y disp_z' + young_modulus = ${E} + poisson_ratio = ${nu} + initial_stress = '-${P} -${Q} -${P}' + inelastic_models = 'viscoelastic' + [] + # [viscoelastic] + # type = BVBlancoMartinModelUpdate + # alpha = ${alpha} + # kr1 = ${kr1} + # beta1 = ${beta1} + # kr2 = ${kr2} + # beta2 = ${beta2} + # A1 = ${A1} + # n1 = ${n1} + # A = ${A} + # n = ${n} + # B = 0.0 + # m = ${n} + # [] + [viscoelastic] + type = BVModifiedLemaitreModelUpdate + alpha = ${alpha} + kr1 = ${kr1} + beta1 = ${beta1} + kr2 = ${kr2} + beta2 = ${beta2} + [] +[] + +[Preconditioning] + active = 'hypre' + [hypre] + type = SMP + full = true + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type -pc_hypre_type + -snes_atol -snes_rtol -snes_stol -snes_max_it + -snes_linesearch_type' + petsc_options_value = 'hypre boomeramg + 1.0e-10 1.0e-12 0 20 + basic' + [] + [superlu] + type = SMP + full = true + petsc_options = '-snes_ksp_ew -snes_converged_reason -ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_diagonal_scale -ksp_diagonal_scale_fix' + petsc_options_iname = '-snes_type + -snes_atol -snes_rtol -snes_max_it + -pc_type -pc_factor_mat_solver_package + -snes_linesearch_type' + petsc_options_value = 'newtonls + 1e-10 1e-12 50 + lu superlu_dist + l2' + [] + [asm] + type = SMP + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-ksp_type + -pc_type + -sub_pc_type + -snes_type -snes_atol -snes_rtol -snes_max_it -snes_linesearch_type + -ksp_gmres_restart' + petsc_options_value = 'fgmres + asm + ilu + newtonls 1e-10 1e-10 120 basic + 201' + [] +[] + +[Executioner] + type = Transient + solve_type = 'NEWTON' + automatic_scaling = true + start_time = 0.0 + end_time = 1 + dt = 0.02 + timestep_tolerance = 1.0e-10 +[] + +[Outputs] + perf_graph = true + exodus = true +[] \ No newline at end of file diff --git a/test/tests/viscoelasticity/blanco-martin-rtl.i b/test/tests/viscoelasticity/blanco-martin-rtl.i new file mode 100644 index 0000000..0fc483f --- /dev/null +++ b/test/tests/viscoelasticity/blanco-martin-rtl.i @@ -0,0 +1,253 @@ +# Modified Lemaitre creep model +# See Blanco-Martin et al. (2023) +# Parameters +# Units: stress in MPa, time in days, strain in m / m +E = 12000 +nu = 0.3 +alpha = 0.575 +kr1 = 1.302 +beta1 = 3.053 +kr2 = 0.091 +beta2 = 1.053 +A = 100 +n = 9 +A1 = 0.034 +n1 = 1.499 +P = 5.0 +Q = 5.0 +Q1 = 5.5 + +[Mesh] + type = GeneratedMesh + dim = 3 + # nx = 5 + # ny = 10 + # nz = 5 + nx = 1 + ny = 1 + nz = 1 + xmin = 0 + xmax = 65e-03 + ymin = 0 + ymax = 130e-03 + zmin = 0 + zmax = 65e-03 +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] + [disp_z] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [stress_x] + type = BVStressDivergence + component = x + variable = disp_x + [] + [stress_y] + type = BVStressDivergence + component = y + variable = disp_y + [] + [stress_z] + type = BVStressDivergence + component = z + variable = disp_z + [] +[] + +[AuxVariables] + [eqv_stress] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain] + order = CONSTANT + family = MONOMIAL + [] + [eqv_strain_rate] + order = CONSTANT + family = MONOMIAL + [] + [eqv_creep_strain_L] + order = CONSTANT + family = MONOMIAL + [] + [strain_yy] + order = CONSTANT + family = MONOMIAL + [] + [stress_yy] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [eqv_stress_aux] + type = BVMisesStressAux + variable = eqv_stress + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_aux] + type = BVEqvStrainAux + variable = eqv_strain + execute_on = 'TIMESTEP_END' + [] + [eqv_strain_rate_aux] + type = BVEqvStrainRateAux + variable = eqv_strain_rate + execute_on = 'TIMESTEP_END' + [] + [eqv_creep_strain_L_aux] + type = ADMaterialRealAux + variable = eqv_creep_strain_L + property = eqv_creep_strain_L + execute_on = 'TIMESTEP_END' + [] + [strain_zz_aux] + type = BVStrainComponentAux + variable = strain_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] + [stress_zz_aux] + type = BVStressComponentAux + variable = stress_yy + index_i = y + index_j = y + execute_on = 'TIMESTEP_END' + [] +[] + +[BCs] + [no_x] + type = DirichletBC + variable = disp_x + boundary = 'left' + value = 0.0 + [] + [no_y] + type = DirichletBC + variable = disp_y + boundary = 'bottom' + value = 0.0 + [] + [no_z] + type = DirichletBC + variable = disp_z + boundary = 'back' + value = 0.0 + [] + [BVPressure] + [pressure_right] + boundary = 'right' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_front] + boundary = 'front' + displacement_vars = 'disp_x disp_y disp_z' + value = ${P} + [] + [pressure_top] + boundary = 'top' + displacement_vars = 'disp_x disp_y disp_z' + function = ${Q1} + [] + [] +[] + +[Materials] + [elasticity] + type = BVMechanicalMaterial + displacements = 'disp_x disp_y disp_z' + young_modulus = ${E} + poisson_ratio = ${nu} + initial_stress = '-${P} -${Q} -${P}' + inelastic_models = 'viscoelastic' + [] + [viscoelastic] + type = BVBlancoMartinModelUpdate + alpha = ${alpha} + kr1 = ${kr1} + beta1 = ${beta1} + kr2 = ${kr2} + beta2 = ${beta2} + A1 = ${A1} + n1 = ${n1} + A = ${A} + n = ${n} + B = 0.0 + m = ${n} + [] +[] + +[Preconditioning] + active = 'hypre' + [hypre] + type = SMP + full = true + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type -pc_hypre_type + -snes_atol -snes_rtol -snes_stol -snes_max_it + -snes_linesearch_type' + petsc_options_value = 'hypre boomeramg + 1.0e-10 1.0e-12 0 20 + basic' + [] + [superlu] + type = SMP + full = true + petsc_options = '-snes_ksp_ew -snes_converged_reason -ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_diagonal_scale -ksp_diagonal_scale_fix' + petsc_options_iname = '-snes_type + -snes_atol -snes_rtol -snes_max_it + -pc_type -pc_factor_mat_solver_package + -snes_linesearch_type' + petsc_options_value = 'newtonls + 1e-10 1e-12 50 + lu superlu_dist + l2' + [] + [asm] + type = SMP + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-ksp_type + -pc_type + -sub_pc_type + -snes_type -snes_atol -snes_rtol -snes_max_it -snes_linesearch_type + -ksp_gmres_restart' + petsc_options_value = 'fgmres + asm + ilu + newtonls 1e-10 1e-10 120 basic + 201' + [] +[] + +[Executioner] + type = Transient + solve_type = 'NEWTON' + automatic_scaling = true + start_time = 0.0 + end_time = 1 # + dt = 0.02 + timestep_tolerance = 1.0e-10 +[] + +[Outputs] + perf_graph = true + exodus = true +[] \ No newline at end of file diff --git a/test/tests/viscoelasticity/gold/blanco-martin-lemaitre_out.e b/test/tests/viscoelasticity/gold/blanco-martin-lemaitre_out.e new file mode 100644 index 0000000000000000000000000000000000000000..922e286d76106a759c1f42b2933351b4bd765939 GIT binary patch literal 106788 zcmeG_2YeJo_o4R^LI>ffDw+rh5Q^C)Khof{LItQYa$ZY)Dc&L(xWHYh82W(w>ZpJ z7>C&k<1kxc9A+zw!)!g{GE^pv!?Gd#VmT=JNR%%r8?8Rs0H2k>FqXMatA{cV>xtjM zUlI5#p)zWf8Jx+iHRvH2RtWy`#pm1l0l*?Y;wvZX3@VFOpQ1Dxw0etjxWS}>GRQ{U zC| z0v+2z4F4_t9_cHlH+$uWG|>m&t;`2XVFfA*@4M%R?*__`e%}S&EoS!459JJIXgm|TB5&FkC7x&9XWK2MsP zTFte@RRg8Sij+U0m}x)kI4q5iBrR9!wg~k z2Ggs+AARzN=$7R1lilLLiY_n)C)ZHW~EKvH(K;f2jX|mg|40{}1*5&vN~L zc1J!T&O;I)5Urr)iW9 zYM*hZiP~$BM*1B|*J#Z~WhTTvqYvSqUuSzJB|!M+*AAaa{<&cpO0|jOjIi`%*2TZlq_S`{`k8?a{y8(Q6g-2Eqt$0Ut11Bahq6~bv;0Gq$MKZc zUmcMrPtBu_qK%1z9!8@`%2IUg-o1B%qE6lr2A#q>K%*f&U8UD3x`4}^BHol@P3QC$ zvjVb`L961GR6{!VifX7zJ7zYRuG=ZIZ=cQm000xD2hDQ>cwrpjpwE_wZA6gcVUS zCfb?6B>?>R(#RX8ba3x8=v&BOm43n(tRveN3fGgOwJ6jeTmZ4LvB~P_#^tq7u^zM?- z!#6^QUw{;xZLB5})H+2XxaYd~1yLqiIYlRxUeTzrqEVw3O=4O^$Ds5A(wvaVsjaYn zuHy(!aI8!>62zVe8#isyFn0*0aP#lHipsN8Y6mg1pjiwi;BaPggsx@|1=~XG5$ye?QXs~*V$)E$@ixjgv zCdcnMMaOP^NkgHppKtE@W9i>?oS{X&FVi&mvm zYAm8CrBkRg3NJ-VuIc|_{c|o^niG4W@pZ%=XB;^cVzO3^7vAr z7!dV}w|PUUvS!Mo7|Ydf7fU)=UMTUT3sz|~n4L0>2+AcfAPa|cLdy%j>%oH508;3; z&z9bn1N7O_dmBJ#?w*9>U|vQ@v)d!RKU66EHa;knwq0S~VFESLi&N>K6#R-6X(4^9 z8&MuV<;ale6O}2;0}&0#L8wmU)M>n$0<_@l2=iF28fRDV$dMtZC*wm(lP8pgkPJFu zCyw3|7YW5F5<7NH=$)9*Lz&R6t>5w>iTF`63>{$dBmT&t(2$%Ao6yAZfeaf4Cd=&& zM1z}SlqsAADx1=5O-|Nkx*?a-QV=0$viPbt!p_Z&9Fd<$i<#qVh#Ot5pK@3>gRaMF zryOQ$JRBWr6gZ+nQqfVs;tgulhm0|-ZheAXM8_GD7(EuY~z=`|HsIHcS5S)>< zKHEy1LRcvN4eF?*qF#}wk3`D&zYyhAqYg_U=%~Dri{%+)HR>WQC;~tVLuwSI;CkW| zt@|bn%jiuOi2~_RCD`&x0J09C07A0>m*iZ(LMk9& z>qBAvuc!cMCIR_{7oHLo;IOo0XXjl50Lwu$1eq2Y5DfrX-k;8HThsLf3gPEiC zy0yyef7d3GL<HxkJQ99w?OR#t)+}Z3-vm9O;mube4lYoJ-Ciw21$cQz&pOJM763 z6QIm$V>aku6P(hVC5AIjq1Q=~BJxQFOIDRmlxOKIKV{y1^7S$8El_U3(967G0& zkQxoCQOSGF0`3b4sZrW!eL!q_AvKCpa6NJqeMpVESZAvGFOqaiiwrZGi-VoM3BQIvw~k=|p6)Tm2({I8IeShTsh zN0QZBr>j`rwP+YhKCf!jNe@zr2t#T#q((z()Vj|mRK(!_eh4qZ7(Wk(D zx^iLo8AVUhPIlp(x4qsCg+GMCi9zaPKXx17am3pe2mpQ}{W2G)^Jy zPFEhGpK=BY1J@Ndnezzgfsf?WXmBKUbLNawJWsA@%2W`L#e?kU-EesPEJYf(I-w3@I-$vmY?i}IAqhkrI3rXpOBWqU+T>lgx{7T z?e=dTVW_3>a^U}mT8c-z6$_xC+*=A2$0OY1gnvhwv4WvEJmDYzBR zFeQU9Wd@vJZZOeHZsd#){q%+?LmemO-f*N)qgACS&2UK=99@oKk4H|SxR0>iM%O%E zIV=Z#uNjvKM@uW=1~QX2la_^(ZCMl&gr0?)CxqH`qfV=a3(UkGIr{>9Fw3Kdi4W{* zvvBDO+^XqwzT>Rt_!Uk?{oH6ET?P26f`mjl_3+6*EZWb+Opb>n&h-4=7r@VMoy8d7 z2MXRtl_q5Zy$#d>tp;{52EgPX$TpqH~n zisZlK6dIZ@d)9L~Dj>joIZ35X6D?5DNNjo`y&TfZq4~0C^KgsAmJ;;C0RI)eEZ*O{ zu}Xg;%bxue#<)HI`h58*aF5~dkB6h&I8#t`eSqmVnpxVIIJjWQEo0tD?iN8y7$KN+ ztrmk4hS@L;=TDvUJdaJfboO=b=q)O>MG0pr0#;6CQl~0aI-OxSl(tTlo}^JJ)p&?p zoFdw1<#S8mZF53|I*Wl;M(Ks3wwj?-l(3GZbfIwx5m@FTo8(@OKp~INC~}D|qLodd ziPQO}uwgSdDf$12BW&ngElQ)2Uo>+u%m0$l>PSwA(2Td}dbDC{x>o89rgSB#skv>> zE`*L&oCjQ#Lh+Rz$~v$d)i75zke4k`Xn}5ve(nvi+PoAarP=^Ok@Z0*W>OCu<(|X< zdxPcZ6i%;E;`)`7CeI(zI|t5Wg0))EDdN0KN9e9aL>e@ap`^C#fl!xXFljBR=`Ob8 z6nq8s-{4%SaPhgDo=Y~E)EsFZ&Q)CLP0t-6DpIZIiJWEBET=CP)s*D zG4e9Zy{9Q4bm?O!7PlX93JsU6vg&d-I>bS{8$zWKY(A~0j5375wMAvJSd9QVqq9?>G)nIr^nA_MXfs6X4_2F( zB=v_5gh-wkM%fhls@>9?9#H6uU6dnXxsjBn2M9?=jdKaFB0+??M~0OON}2B+(5+nr zA5ti>b8sLOC6BNyp@=Dk*ii?yltp@eED!RKQkJ9wa7+p3v}ePd0=>V$oB^etoXn3k zld@nyfpQGoh#hZ0p;*(EDJHE(tVDoDDaa5fk5IG~QC zC^9s3Ae7qMI1u_8k4UfB&MBGrm_&M8hMqx{N$fpfRYLDgzJ<{L(kUbe!M@-eAW~F3 zlb=LyH`QoJy00a|^KSYB&qn{-F9erI<&Ad_Un% zZXOXr!GQhW-S8~EP>>v@kuzxx8YP&mgrgDZZP^_m;n)*po)CJX;>$%ILONm>$PZBx z8A_V?fD6CgtfcM5SzYK4>7`5%!pJ6_tn`>Ia3!o7_P`q)TZE;zWp{*L{WY>IwR7N2 z4ZPVwWg*5+BGM!-1#Pp^S;l@f9n2_N6&TyS2N3cQf)Z8s?%Ex?ZQePYmIx_jsfr$M z+)08PsN8UP{FF-JQ?LTwnKS7*-4lFB1lcv}xF<0v8`!W!13Ss{a^ctJaRXZ#k&{DO zuTXl3hDL;b7uBVAUw*gc%W2F1|C`t7HMMQE(5T=#Fty*i+)#$W-@K+XsM8R$*dwP< zoEN2uBNX4;C>X^#gNP&K-S-4~_a*DBTJXLAg@Py22#Q*&3c%0l%tnBU4@#fW<7$BD zUP6kH-(BadXaRoUSHZuJ-p3A4LZbAc2uyuAPQ&vj@?TD&;9!LA9n~n|%E(kV#BvdW z!wzRM8gwwX_^)4ZAN1{Yuy%t3`vBF2CO;+H>l~UoC|m5Xf&^*gz;C}mgrSiGO2PHu zNl?AXwb;+l$RXFeaY*Y!VTDExp^<}xH9f31SZzWh2b6;Ak*nyR{*gnjyBri&fQ<(K z%Oi)NYyw0X25IELZxBI*p^*bh!Sy_oBZpk?7YdCWeCR$oh9Le(A0NG<9E2g81Gw(G znzD@53NR6bY!3XgrbCW$#x``}$gL<`Kfw=HYc?u#-R0nyKH3P66fUUxcWn+&!2(TJ zuqAQ;m zS5BR((u0qkj*8DIGOr)_T~e^uQv#fpAl2m(D75O*l<5YIRcH0~Baow)9a{jY)I6aa z(>>Hf@R)#ooo+rSPORV&)V*7Rm+XB;;D;Jcty94oxY?TQsX{SPj-{PWZSlx*FCoe~ zv^?ksDh7Q$1$Hd({TJ39SPe(fIiEO1hSsb$INtDC7S06xP^;I1mjgPD<0_I975@Wd z7}+GMp~Cquz!8Y}74M!5i&6&gBLfEYg-Ya?~*{qBa z?hgz_qGFs93j71V4gQ)`u$Pd`skqU_M|UHqw}I-S0$+Ag3Ejm1U6aZ(SoJBsS#A_9Z70`LadNVI^Uv!k5m=J)$2RViI2sWm42^A23a*EYZ49sl zM$2isaOo;7^uk@SWA8&Y&v|t}_@!ghf{3pmBuaM|1cW+cs>(%^=P%M9MmP1M#8(7< zjULHL;jRygR!(n|#9360C;}k_XYr`!(NcC;1BHMvT8M~-k>zGhV|{Vx6_4LN`2)@; ze>nH#kG|4OK5HTIGL6c8rQwGnX(pewkVn|CVRNtg<|_j0A-`Y9aSn^e&%pcFy z7VMllm73G^Fd5(|mvo1b3CM9rhrNy8u@@qg>UHmT`i%&sT1EUG>nVpqzn*i5Q0hJ@ ztzNCOYB(jBmtf8e6QGn-16UL3*#Xk+*pAt0 zQl_RfA4E#8hfajzL`p|UhDs(Y9JK@+HK4GyX4p_^g2SftFf=DYU_194HWUg{yx?c) zpj2QK0c28}o;O0MjLG24$!`#oUdnzijFHn)%;XqArNN}}*e@!-Vv|jjFf-!F9!e7j z+iopjx&k|nDj8M3iAx6KSgs2O7U6JLrW1>^4d;PE-1A6w-clSMPb1C!i^|G_JP1gs)QFKk$&_sM6r$q8(h=eY7Ho8Ig$nc! z!T~H$C{?GpBgFnefKnzxYmx>|1jJL8gf!+Ap%;xZ5aKy~P{wdJMiMN@QRPI7NV?V| z1_GtY02@)=H=w!@N-Y$ee*4BJqm*IWyjkhKDo<&n^-!v!CkV4VvABh#%a9<te_f|~vr42RL^F_%Qzd)R#?uL{B=^HtqUiN-B8FKjP@XQP zrOQz7ySGNPQra=L9yqwkcW*5dq3_;W9zr}QMyJEmU_}qJgWH-C+Qs+jlIW{95X8f} zqpM#qQREQ7-Az!XZ8VRmpSdLx=>0;+7z58-Qh^M5ct31b@(tB*RK_@%_5GkXNbe)28JvdN` zA`n7Z2w^q~gD;imc1#Cm>DD>NdWujVQ%a>5irSn(KNf@15o!#{=_()B89)xw5ptOv z-1@+YdV~^$&v)$BF|lKO7e!*nt_i&p6M85Uy0!J`GzhcNOz?p9B|$l}fMhLPXW)Gz z9+ju`ik*@k{S8}S;1GhtmAr6Uab^wzPir%fi65#}s#or>GoU3xtPXWC@2U8s!ihkv-UoO`<$ zDjAqOIuXK(R)$v%lHS1z5Q60}9WbX;lLjD^=zJnXNJ`sD8GwUjFoyG2+JG}G z6bU)*NduO06{(;ySawH9k9Bi(XL1pu%@t0&ROwS-%@Ef0478l3S8OAp&~?_m7uV8g zL|=dPR8JbfK!|B%|PQDwnC_S8M`T3n8UWl|_?T0~^{Aidg9gjziqGralEBSP;Pp9~o~CXcXoD`Xexl~17& zT<}2YJI*9_k_wyGB+o|q=`A=H!KIm;$?U{M1S%gP9#IcNPLCqa6U#h80evV#xh(Ss zF_OV1T7#&Yhx|hEMo5H-)rMFIeDu)5XfzScI?%?Ti1kcxmnx8JGqj`6hLe|EUVf$u2!I9bvV%)7HnXPmcwKRr7UGL zTo`|A7*uE{LP7s1y>k#!C^(&7Izn6lHyMVzCS#rj11w9LJ;Kfk@b)r;)q!`i^fEn2 zSEG6koV37E%lq7XLMl+mMM}Cj(tA?oQP9#$Sp=bbK}#)gCl9r8pwyFR~v_AiP%BF zw@-=pkvaQo><^_5MVtunGCXe@WmD)ZYynQ$2!&lg#Vfl|%xV+IiC8A3BQ)ZuUFc{@ z>BB;4?7yA}QXSSg*t9L~JuRb@eHvJJDy0_+MF?XC7%F+#UlFAYgr37VaOiZ7!{i|} ztMQ(6*jvlnMtDYD|D{uih9K9q5MH|tWFYjdPk5FaMa4FcZQ2~KB+kV}(csf=&}HE6 z6mL5k;ahs=z<19+=TdSJ;)K~%5^ld_AVk|2+~p~>qB)StC={nL&0!yLBM*Mcp%CPk z?qSGvM+r(Q8SF3IaDypLLELTeVt7~Y$q#^GMDEaH3I%U&fD07UMLI8}P^yw52w`>E zhXqW!3<*M!g=J@&a&VqVnj6&hEeTnJ-sE3leFW==aGlf6zi!3n*XKyA} zIlo3z=_bjnu`3 z0Cr;O-4}EZq8_z`KBO>#)B&GD6C~hgeyy-_Wa#E`;~>ioNNUJHEq>~T(A$M44-^`p zn=q@=jh^!rITS)~h*eRmw;H^CJmesB8B3uXBxb}s%FwUjSx!ClYj~DZDDLxZH)i>t zvJg7VX)am1`DO$MA)Pye?M0Rh9TbXMvh$44)1-?ON;OEjOp?4M+f9C0EP0WksA)Ve zEDLC()kdn7*u}D7IKt!43(L+qo_oRM1)->|Jufm8wYBGo!Yt=}z-4RC3x%SVY^iMz z>jr;Vw4-aly{%xJRPy#Ac@^82*^5*}IWlyvNb^N~e!b`9jBSLLY!6#Il(Ixi zdjMOz9F~O?x|?8Sbzgo><8l?pL?L1w5kctEta4;1g3wKj5-}o$4l69| z_vG|mQ(Wmgxts{Snlxqxr7U&BixZ)^OQW2Y;xs&ikD!+&TcQZVkSixbF-vwHScU~r za?;=_LjMp+cXd}6psYZ0~g$wBDX)-HWRG}zh= zX5w)}#pVI!)I-nG@U+>;YI^);J<(&1<+WwUAW38z1iZ`cv$c_f&$k4CT zrJUY|H{Of3WO-v7a_E4sIj$@j`dbZvwk*k*hO+l-yq8rd{;s%}h0vezUX~1<#(P@% z+@P##<3uR38X$*dVtB>{6|sdgIX$*tOEw^sXvYp<$(BTia1o|FWeW*&!Ql z_+h0`NJwGG(vDJaJ#mU2AhhyeFie3jNk@(zSAxAF<-74ao z-Ue&F?8luEhs2=)r|(8WQn00%o#R6c2_8B1(1kEG;6!oC=^TPQoD1^k>nS(}qz8D~ z81(!`Brwh8eLNSSB}q;u1J2Ye@CXVLZCKR;TzAr&vsjJplHp9ePlkIONb8l^zPcjpU?zp_@&wa2Qa|K>J$_F@EC_cqMRN(=mXAVtIp?0 z1og6jDgQxaISuZ@)Wb0^_}lA=gJXfAt~%R0rwNQOJvrI;F_&?Q;YR-`JX4V(J)An4 zB`gd|svHx6Avqar4XPB_=kM(&nJ0v>XM+;vgmMrX)#RX6Ck~Gva#nndoW%^Tutu-( zHDZoaM4I)SS&3DF8}3;EVj=v-10k5Jz)2OO`3KxqBY_3td4)enXY2@)oK)al>BtoO zLz62Ao8D;(ZqtDiB6!+xVbTnn5a3KrL>pk&29n`PAvmF1Jt-taDUu9^bT}tFJq3oe z4Z$4N&;g?r4dF%;EYAE9tIeSfKw=y3qJ&*A@eIdB_8}B!k~{Fl2z$i zUrzVX+jNO<-z&ab`&?cTiWCOl%hYkTrBkzEF7%Ppmi=CmAV)9zx*dTK8K-F1u~+Zh z>=@+>8$ZHzJwYEs8AcT#)k!{}V{KQ1!OYRR*BY*U^Qd>i0FpxQ8zuK?5N{q_zi@0Y z!D^ET=Y|PC^c5CInBm&SCWmE}$)xhJnMB@f)aUu)NWL5j@v2{gUhL9af>8SSkkh7! z?Sn)5MP=7#3`px$jAY1AQK(M!a>J6DIPZ{cXD}=)%{#EjCDNH zSq}2Z!B>-%9q&@3LAaP5sX$5PQH^GKaDxa#YSeR$K~75vsZp#+t|v|rQllOX%{zd? z3aL?7mBC+f2As(uHR>V6|7A6r{lv^hvppI#7C=aihSX?iHtNwVL;~1SLbFkng6omq zW5c#4JsJ>QI6QuZ)Tl_!^Nd1yJ=wF-Cx1!ZQ=(69h@F;wvi-91BPfN;it0afRKM93 zVTUrUqG&vYBwqjbO!8j&3GdN(A+#NdX;Puaye0&`IPbn_dXrIeNqSVjO~Q0xnB?2~ z{m%%zFnte!M-$&5bW~lii#f7=!R24tzRVVGKi$0e+zwWia%t;;Tjkl~j~_oiySJ6S zSff`@)f`u4O9co&eytO`K*0CP)vtDlS{%i!Tpjz(?h5aF<-c-u%jg3a?T{Yx)0t6Q zh3PbGJGyhRQ9r#iz$FzM47o+)DTMehUp+P~7GcJVOIp!*Ul8(Z(QxF%B+U2HH>VW< z(7?hmlBP4VDA%&D;*ays>@%oh1AvC@&xZ?knGyBAEKzneV~7IwT#ba#94_4gMI zAHI4^#jvm{n}>8`m#yo0aZHVMG(YCs9e0yw?&$#V_O@2H2h(^8F}?P+j&T5Qo`3Pr zKWV%W+K&C1GML%_9l{^)-z+bZo;s=O+}|NRbyfIEn=oA%hWzfN9TV1Is zUTZP)CcwWxS~Z{=b8gY&zxA`Gu|*G!8a?<`EnDRmL-OGV4Uq2%@$I|WOKg3y(85Z4 zknb{gKhv%QPDKF4lhLa@I!b`Fffkv-(8EP6ZQl z09Gizc(n-MCn~0^_o zv~DU)7lx6pVkJsd#(WE%m>UJKQk_u;Phxn%73Egu>Z*J5vXqb6;>W^voQUSwlJzxb zk3P7He77k*{GPqkw#BV>_UMOv7x_-nmi-P}ch5^SA7+iFW67Rj5dhm?pZir!8c!jn zcZ=)Y1Yo=PD&Nnc@j_@jp1p2l>JQ%|Tz#d}1d;Uh^;bNq2kGlu*e|9E(}iKkufvhu znf)-|Qe`b`kd8n8bMFfbe;rkP8*}Z@+SbQ*j$}&=sF8nvT1B?#zSZx1TkcKdyG5Bl z`q@jj-Ikr4aR&J=ys}MW_6I@xUURgHecRGL%$lQ>ZAFz{^g?{q!b5j{2(VJ*$qqi( zEp5lIPw198bs@sH*Oc^ftH$^Xm!~8me3)AKrjSoz81k#~@vH|!0am=Ra_oJiW5Sfv zZIF&>gACd1EL&10LrPmfrYP&uI`7U(f z;cj+E@i^oB-;2|Hn6>RW-Q1790H}Fm%_c9O3QS+UYe+0Wb>Tm{enr0%LfeU`3*&C+ zj0orGvQ~%Z#k_dVCbA=LHW9o5=q1!mX_z419s+$E+)I?orG4Mj@=X z@amg1oYKxsyfS0l+@FWA`G(fo z8b6=3@B06gNU)bV)at>IA$_pC3tHd4%I+rp`C6I}vu^P66mC%^fXhBw_n%JLuqTom<&u<1*Sp3cr8c!j_4{or9y#sK-*~D7) zXuJ^GPMk_U+`s!JgqQYi*d>ylQe{T*?;t%zvGQh7VY)Dk^lZCbKNj6moq z`TYO~zFJ3Lp1IX#!6%za7iUXNp0=Zt@iOyhS=fq6TNTK6t>Y1M?PcR?Jh$lE_Q-dE zKfYYX?jd^YBF%?c_i%2p-Wvx4oYS__qct?1LQL=d`HSxZ{H&JZ4X=3%eb08Xc#-R~ zyJY|@T5#CO1R75v#y9MG^a{e*7oJa|@xCC^b9MPIp8=e`OR-UhbgW7)`vB>%7An(# zxpgRZ>fSbo*iuLSY*S&(bIikns|NJWxrTh#{Db+xUN))q(Ug?_gztF|nz8$0$V2mE z*0&G4+pO>ugv}Zc^IV~Xmmof0>xXSVVb-^Mtcc&ApMEEVwv$oobc@(!2;aJx>E*DA z@k<`QaTDP;a|cZkelH9oJ;zsdT?_d>`f1tm#z@E9?;6!aI!xJh>NB^w#8J(@oynH2 zrCax0x8IoiA1qv3{O3Bzca0juJKM{Sd9T7w!z9AD{-t{Ceu4g8Kk>%+1Fyf&te<#& zXP=Kxd*?s##@>A8M?w0v%^m086{gd$?PSjtv9+`<5O(=zBjKR6Um_-?gkI(q+4i0RpJ z|5OIJ#}u*AO5=slcG6tyKj+5IMmUG9G(;@@r^4}vF#XrA(O&ve_&(-a=~THQm~W-h z-`B_Rz9RW`82;{^I~wL6o$1?!ZOX8v7u?_7w?jAPZjJY5oOvM(`L5O?a+$sCj-GSB z9XXEhT{-(Vmh@`Bm!SDD8_L}1^k{KgfY+z|Gqw~ zsg0_Cx@-u-$jN%AU4y=d@qKs4gd-eyWpY0n?+aqS4G%-ThJ2fskL-@&^B=DJ1jDDh z_0D4cIk&at7n5IPOXtl0u<+=O%$*A3OCH?&8S-7Vz}N@&vS;g!y0x`A@}2+hkLIz5 z$n?86&5zm8b!vn0w~_!(u`C|uHEzcAVcV2?F)-Dm#cYS!Ty)a!E#(c{>n0yBFRam<3#qjpYwH{&k7ww61%$>pQQc`MmXUjZy zX5ot0GMInr{Z?*cWzw%zzVz~E_V5DRC&teBl<>X3e;f8N8Lv6V4a|n|S4S2$Oaypk zXV^Q=`cLV=^lfuWEd+S^n09VU8ZU&llfN$CP<7^Ygg>iev|{Nu*3OT@^gACF+9*sH zhB04Fp}og2-!}?1{R-f<%F9zm0=)9V*T!t-&TsmtiKTCi_^Jw~?>I0!QJ5|aW4>XF&RxcQYfe6q ziQ#WQd3_&-e>t_f2y>U2ko?1_ zBz%9mF^N4!#;GsSe3%V)XDM604fW{m%w^3A)TQwhV)}yV_9%ceb{_BPHLjuW*>Z{& z`npJ`D+o7jtyGf6Q;6}WU(4Eo@ciQuooKu-i21f@e*R6&H+4+H2n_!{+tv}_%!(rt zdNX&IREx46oX(c{zDMB~S}U2qcPS&!M%_ogE5@cawuko}v8napLxk^lC+4t!;A||; zO!H$lzVd2&t-3P6mLJdU_zR7v5Ys2VlrR@yiw!^YQPX%KwB^)3SiA3uGYA{)OT8+V z-oL_>Fiao3aL-s_x-g9SHaO9HBIet-yg38IM~+q~0kCDkMO*qZ_sYI{BWlrlw#@Al z=jv}*#{6CP<5S6x`yk&HGP3vD!;NiUUS;1$_|^}OW&b4Ob7%WvHmc(iBY!Ok@YRJ` zz3zGEua5gZ@FvpAHK&eIji!XMRy3e6`}O(XE(!7UuqexB9VV%NuWf)Zqkkb38YtLxCd5cljR{ zRMbA&HldS0Zw`&R5VYv=(?@0V{`;CTt;TlK?wD>2_5NmFYf9Y21H-=-#4^z(Fp+eZtl{rEsduT)4%RKBMac3D>+4w;-BA&~n~bl-lJ5LuM;HFvGiUAm zdMC5-_Kd3ww3}%>g%CgUPvz*BF}>a!wRh5ZA++U;9n)jv=T{LbBd>VbQ!svvQokAD zd-Jvp6n-xZW4>=z`SnwPpKOYGbQa*O7sl0Ch;-;4Gce%UG=i9@7d+zs3S$e{Ey_PZTNirY1O!Hwj)qB2Lz1QjitZ!KIpY1fB zLQEgjBKu2#^_GnJdmoJ#LR-%4os%Ez+Kq6|$j7V1(zj(*t%B)4e>K3XzZbrb^xVI2 z;A6~}W!6_gIu;F`Q5VD4#*G}x+&}+9Sn}Jo*|NjFYNDDxk-0wN@Ha7ZzZXEh!|z1J+QT38ywh-rk?>vQXES?> zOhcQ{oS02-9j$jPu_(YJBS!CDN8>5P^pyOw8Gwg3{b=}`#tWe>=jVkV|2h@g_sj52i0! z*l1ZJw(O@nKPwVBi@CPIlyB$F8-(w7Ca$xW3!l^X)PRM^cfNbWwzH?m#KhU3GMm1v zbGPQL{0O_p54b_|r4ZulHb1a20$}Z_f6psRqP;b&!Q6_Cr~LRU@?CaW z#v*&U7Ta@nRUb(B&M9wV&%kQQ#WA&Le$1waM^|XSnTN3Gz!f<(o;65vun`xaJKB>BVP=d(}lTO<=YeW|H>wOciC6lUas4+A+@(1 zCVU^fG@Cs`><-&#KFsF&>RVgpK98_#$;;2vcnTq2(}aod2QXYj=_FGS^I>+RG(1_-EA3mk8g#DF0#4CgY=}`7oQi zRVg{^Kbrwo)-+3)O5-WS^v>(nnf-uGT!-pT-NJ?Nq{sFKgZ}gmCSQ$nIk4$6Bm1VEW0Dg?!HA>35i~B6Ug< z=KI!=FZN=1+aP8qz=AQgK3&2*GL%-2&`o2*BkKOrAv~VB99L)6r;ExEz7u{fY%e$N z&btps&m(-VI~&fPgW31V_A6ZYXGPcQ6Vs5{JgZUGh^5YPyL-BF)4As$z2UE`?|P-v z_iU&7?Wtnk7m2X`z7mURJcStFxp&%hggvTF+(6@fL8Ry5(v}A>-@#qBRtFgUdDVq~ zVSI$X!!hR3y3D~}z4$vDUT|J2`q@M}j%THBg z{~~s#J~TgO^OnU~9qSzj_{E+thkMzPFn#Uok8POVAo*;4`kfHkPNhvuI6kp8!rt}6 zo)b$SK6YOdOn;+Yyt7{wzDL8DZ=a~z9|Byo?)5VJF?{#Axb{d#`r)^~VjhRR({J>G zxSxxQp3*#`xx& zZ8`%q9E_OZMK7jrd{t8m(;GEU^67_UJN4$Z4%{!{2>VQ%v4ql3A;xF@{`zr*Z%rD% zm&W^om~X!uA58?Pw@o^)#PEaEZMTpPi+0mI=5c#AV#$t`YUVD};|KoSaGVQwZ_-YPY;P7h(O29mdgkUl8e8eD=y6%(v{3 zHwqygoi9{cj^Q1pI~QjjPu$S__viy``0F2h*C??rbFpsKFK5Q9NIkFFKi*#MS9R9D z(v?Yhf84hYdjS^5elH$R^J8q05o3UsG54MH(~svw1eLR+*%-F^F#WbFNpaL$zoFghX1^DOn!jn3zXjX64Eg$ zws17_cu#)o%}1Bn@U(fqg`a7_Tv!^*9}e2@PwlI52Lo#RBtX6n`H+f_*bd(OXn zxg5=xLQEgmwtsa@FLHdYnZ^sD?bHvdhNH`MMHn~W?sc(rRZPqlOy_pBuO>_vhB4o| zquSR8*kj5~Lo9~3cd~3lIy%)oxttAqsbb@t6UUD1GGwji(UfFYGQl2jNxK*MHG?Ul8+c`rzI%fK|RQ-fV>7L&*#O zgLDkOagbxfI!~FaTGEpZpWR^efY)AT&Z}a!S5a&ue82czwml;JQHx8fiGNGJ5vo7g zOTXZ=hvvuFRu8;g@$7pDyT8B7%bo)9gQ9D`{{g@OYfVdEqu&Xk?bMmGV~>UZiE!%V z5q7ckwatc(!}N{YE_;n*h3_Lhzdkn|>i6ma@uOlEAsxn-ep-RyyFKTHv0-TwccjhS z$cBIY{B-;B+st40^<946v<3Mtd1KtC_J~^dyLDBYNO@bU++r^iJF;_qkFg!B-lWu( zT>z_1nl(a%?}OEi2kx(h^s09>nV$;jqhZ^r+wJNv`g$b7n2xDl^%&#(9u4n-a6r;J zFTFwWP>6JuTGg}{<~!=P>Te7$8@j0&z-lcb#+_lq-Y@n^V*Xle_@*C#{59j9 zg46296TUYOZ*Pxiw0_NJ>-hfM5T40i!QtAnDlYuD2VFkiYBOWIJ*fTYrOt7?Tl(M^ zQx9>FKJeQbOZL%tA+((?oKj#w&*BKnT+aWiSbB^$E*aBfFLpX6Oc#ccp6$<9_z?4b zHFiZCfP<4KSSn!n?>T!Kuwh@8TJvye0XF=X)KZdk4s-6rYx&fjh#o81zxzUaM2~mAn$k3a@U31xiM@*EwGUJ@ zAI9D!d*g7<*#|^4VoW?nD5k1H#%azC3h}(2yjG9>9xHuyzPhWQW<(d%y2aO)p*vsDr z)79_(v=`vOn{9@8&DSV4+vz4LwWc*_hA{GWr}i|SLX7X(^oOzt`z%TKveVP|C`5W5 z=k$9I^BqyR+Zm)|&VYml7~eQy#>Z^Fg2t50KaQ~BcV}>o4}HO$YaCIc#FAOacZmms UD%m6S>Cdff^(NuFzv8F=10JZY`~Uy| literal 0 HcmV?d00001 diff --git a/test/tests/viscoelasticity/gold/blanco-martin-rtl_out.e b/test/tests/viscoelasticity/gold/blanco-martin-rtl_out.e new file mode 100644 index 0000000000000000000000000000000000000000..629c636c5f16de6f14f3158a12c0b9d92a9d770f GIT binary patch literal 107268 zcmeG_2YeJo_n`@dmO@hyj*bW+^rG3LB@jpmC3Nw;TsFy(Yg`ElMFa#zqzX#$vm&At z0YUm5NH2nPq?gb;hzO$m-uh|6ip_132~+9cM210wtxEq#hcWr5U(4)w4H`phh!Q#=#;V3 zQ5aYw#7P2>hmwfr=$H_Rak}Cf3*b-|{P&s);F<~Gn%i+DC_g$9ZAf9jyWlWdJ`SVh z<1ku24x{DcFj{w9mePng%o~zk%m)RBMER1uQR~w5@L4!MW1ef&I>_^o?)VM-io&mi z(x6sku||_xuY<>sBJeBFw!qJm0T%Uwub8IME6r+My27AW>&%KFdSfc&!BC`)azNQ) z-V{$$YaoBQ0GBWfkJKu&)mki5P!1KqT)H^}(2nnAa9xAln=d-AN;{pswVsugZB+ay@8fPd>7@6-xcQS8{xH3d?~#MzCeQe=Owh|UzBd3_4BmH*58use zd}es}1(hC3N;V;nGQqM-FXX)q_fYCXfZV$|1mo-7-HvyM`g^xD@``zcS%KvV%MHy> znvPIKoLZylt<+dZ(IM|({CL6p^tN;YJfr+wY0q?8}gGgHCkX9(l z*+GpjB0~+5&(y z?&74aGUpgL5p(&_DYkD*cn2?6dK!2xN#1vlAK#%D-rJw6d2yLNu5 zZYu>idb!e_4@gSa-@6O(T{~T>A4^lB0A1-WaiQz)-Q_NHsg5nhwpOb z*HIqAcy#U0RW`W1OY;qNC%rIj^xa?Kp3xAlccrqF8fdjRdRW=vLYn;c z!P&MK7(5rk-M&LX!{zV0TzJRrI}|!xe!er-gm+4jbdy4B?=^AXRUj0c?UJ9g)AqW$ z4898|X{R(Bu_?UtW1B}W{E&MtU`&O3_ir2L`CDuY-M?*;=Wop(-!{$jx5GWYZI=!Y*N9%*4EC=5U3L+D+^;U)Aw36L-O`bcPuVJ?LP z>GDVM6uJbt-O4xE>DZxLAUyaB_jsSj7mc4n{`=3{bBa>MjclK{gxk}bqBooM#3b;1 z3c#pN&%lQCdA-bf!}FSSl5P!~_B>~QNSaZvd!7XlZ2yDp|MOh_gYAE?{ePaz|8w2> z2_`-$0b!d4+P84N0&QGEHV=L{zQebjpj`#=Xj?&u_7#NqEkakj3w=f!7*4}zJk&np z$`iHMAdmDr5}v9y85G%oeNHcuer}%axs(7&KR0*yT&ABLoTX41S=Ink&qgJ z0$0AOM-5qfBMU7|sw^c(mY9^(Ggek3e*^>0VGZC>uhlAbsj>ucm6NqKrdzbE&TNuF zR1$cVHO$a!*@ntNO4VQso6``Srv`yvjHzJ?s|4S=1}3G(Oz-suv#fzy=HbQ`rIs~l zz#5JEsHmtAqy664Ko8#=!F>`C8g+_UZ&YP8NM*CsDpsa4Sb$wwZ_cSt77vr5UiptR zf=dAS?`4wLj9PH()9YGGNGkn=%UFBdt!2(DU2T@BK)3)}v}lp0YMj=rag$~;y-wCR zx^dG6jT$wGZb6B-!UT0nB5TT^$#If&eM!uHx7g&KDo=7US(B*7jiS)!LgdvmA-0=$ zf_A?EF*Ms-j7BJRvSe_}b*2|YnQUQYos~LSqo%S(jaoNr(z`t(0^VZQ#@LEt{P|d~<^SJdj>I}(V1Of+$rKiy z8X6jfQiau+cmO&ppiCAcYnz(dgVm5qZpu&_5|uh-I&1XChYTvL&TQ0c!1p5EI)qewLi zk`&LOMn42R%v1ChPyzy%cn*~sAgR!?=^*wjRxx-8i?hfMamL8RwGfag78A_EgL}MQ zmFXbkWK}$ec7ot711iHwh~9~pA-Y9D1EwGpDH$h(a@Z?5$y=U!;$<=~geXbz1c?`= zcdFtEqBkpLM*zmiI=1cAv!`d;ykG*}ICzOiI(k4`OLiC5sADyzKsRL4nR5G5CL0jd z(9^u3P+GF3ag6!uvxy}QEGLw<4o)LY`eN@$I2Z=CkrargTVzafcDqL=RS}fn4w@R^)}Wtu9Qt z-Aa)m(I-k{jvFERAcdeRgH>g6Y6{4Lvm=aS)hZla!6ik8td8^#DNpWD5<$}GWEIfE z_q<|c$?=J?J(FX*DPp^J@R=VZ5GkZV*FkUeq(v}NIMk=5!QwNo1f%`4JYm)NP9AN|dY^ea8D3m^ovNsO@ z9Th-X%ZuQBy@S;?r~*(5&dW^&^nfzpM-AYYK^5SXAKnF3fUjCd!ca{ALn>eAZtGgAUFzeO3rx~Q~`c#9}4S# zMFl`L@taN!=MU4DL5}w0h|>G ze5DL}9ZdB4R0H^Ba1`K_AKpC`DnPAEWwRA(ch2Zjj7*$ED1OHg97IXDLn#DXT08aw}hwD9S?gB0?v5SJqfT{_@6?8kqMOx=WXVKFgf|b z!DGPJKO%<#*UM>2ty&}O4}`;t#H6lCiScdoxP}us^y!f;#b=3nWJA(Z!8hC8{Y#2R zP{7>f)dCWA8`FJ)dIg$aUN6b9fWCCN&oQR%H*13~dx zFW@!@SiQuVI|O{>fkUxw{LuSSr?W!Nk#-r1PjbM!b7@(G7vVoCnF4#V!H3(q7AUts=w;3@0aTG-P>lxF zsOYt3e%A#A)hMmBULaPzpc+LfI4>!RKBz{W^259Tc{NIWe)6SDQA`0-qd2_P<;4%l zOFJNlMh=5|IjEO|dRee~xFE5n1c$dM1?MGR1^gFgD}s92nIU*5K2x5h;Vn!dAtXY7LMK!SdH4nRflQ(+_$ z45)Lma5VYAB8TABZm?&Hf*g;UpR7TGKC=ba-qJE4lrubDLI^#XkIhG8Wa7?rt46Wlr!|j$Mm}b$D4S2aHw|^6eoF^EK8|@-TVcw7@6N4{XLTWbWg-#GF@v5IY|2{sVUr{-dslbtSRD7|K4Yh= zU@g%yl`DA|Lb9le~BzmBj|94IJ>dRcf=EF>7z%P0ls6(f84^>T{d3@52M z3x@AIOj}SdyB`D?)XToh41{8_3F_tipX2kd=w;#dULemgm!DI;JoG8h%R>dS0~xA!hdrMH-5 z_x1`gF4w<4UVaALW4QC?JwMrql%b@sN=1vY)T3pYjb!!jvEIxaZ zabG-~OtJF^`R$f;2^{KmW~IukfISs~DyuZAG89UUMn42{Tcgybq$(9E+(a%$*4S&| za|z&Vb9{mtvz``4@tLBsm>^daFps2gl5q+VlFT+}mUlS@JtU)IwUzNI!`5#lPz#)hGvU??g_EjyxK%eu>pbt^Mek`q#QQNJBxnS220T? ztS(i7^H&a@+;7A;4y@4#bG4vTgwrmbpsN-UsL(`)qUy36K~1{es5WP4oym@oaRt#0&!RIb|E=_M#v7~x9W^u(=Jy(LLNY(B;a^_ElP$C0MP@oMXnIfwBQl~?cf)oDS zT%1V(f`lWC%zCO-#Wk@qnZ|4DPqJEkvC+WXw?|gE7Mn)Va5>9`pF8AGNH;ksax%=j zrzsG0>SG5Lmzx+_Dx0RXX!6!O#7?^lL4^TqKFV~rb^>Y&Yc{F;B-9x!lc5B&(^##- zVn}r!^6uaMX(XgT>}mQSG6f zAhIV!Q8I_#O1Jo`2ON4!F3OQ4xsjNv2Lwq&jbjOqEJ1>KCx(Rza+zx#(5+p79&#wM zaj+85yHv7ruVE(^5$m>=XKxhzTrV4o89X&(w>3iSR0a|YykS{m2aOv!-` z1W= zc7kGC8#_U7{Som6+c6{)?vseG%g{0?vx&V2tV-y;$<+}0FP=j}5bO(%4k9JRLGi9` z#1jwO6yZy}3LxWIMV3i1Na?TzL9fZEwQDcqATMTR(A;cC?dzMPj74m}By4w=0-hUyTkN_pxWk`|lVywZ(& z=z)0@9`TKX35*`@EtWIAa2F$^PN)r7i-tAYb&*#9i;XkY1_hgqc2m@NoK3|^Q@jjK zI@Y93Qy5tkrBE*CU4^=eDOfkq%tWP_f zJ;hl~g}a>Mn=c^f%W@+jQ%J@0YG{!{kX5GU;Mf8zyQ!=hY=2OE98%0DL9U%}WjCJ$ zA!ES)?`n7!pDDCKi=(*E8}YeJ2tv#{3BkM{QbKsU=Ai-eI0j1!)a3`pq zr%g;iC|K*-T zKo$X_3i6@omE-F zYz3GIf;I~whCZ1ZIyoU4j}Cl-oh7cy(SPhLC9@VJ;LxJURA}|77LCQzk3foEwl4vsQgeq=O!rU@ z!D9l}b-MVRIH-b0P*T@e57~P~;D@QKN~451aFZp?U4>$x6iYiD+Txbw9zv9|X?egm zR1Elh3anV*+Al0SFdL4db6hd9EVW6cx4+>vFB}2*p<1T~F9&oQ$5|u?DeeYj*q~Wd zw9N76XA4C9igQnfNhv+}@tD=D8%-&>?9<4;3QW- z=cID<7G1jcBsXd-ZY9?VJr9xU1$VS8rT%ptn`{MF|sC6 z(arL%rjQXN*qCy}haf}o1fTjo_S4^kefm3`nMSXrkW{^>CL;D_PuP}x~~ujOH)KQuuh+Scjm~G6b^Mh@&wz% zsHUy*o^0ge8R)dbctQcoqQJ!DIe?ThtgtU^4V+GS?)<(2+$DY%ES%tbm*0N1Yt!SW zhmEZbw+UMZux4wNDmJy7Q4bsTXziv%V8@jlc13*lu#uqHwG}>%EJ#qS?aSxptx`Dj z=^u#%^R7=&nrIVIX)LL%0y;Qg8G&JLdWIfs#&m}W@jY03^Xg6p+g4vO zrs6Y&Bj~LID)3fpfHo-36kES(E=N4!EzkPY_$oN%^G1#F&zxT$Ri{> z@5^#zGLpUEO>5_rXZiwSVym7fK`4xA;Ihjt0~DXjK92>GQd3N1vqy#AnCiAFSbD)G zOJ!jU%#xL|Miv%Vo54T{);eVxs(vGz1}4xvXADfW;@VpW6-OP;2Zy+_m8`F(G~A&0 zT=u@T4IvyO>&U{UpXjHZ=+WN0xBlb_wwKkZakUWI->V6!$Ei zATCqFN(Tq!K>NV&VgrX_b&4xNY#(?eB@(oxq{6H_P7nFP~K=c6szb7!5nug zE>GfRNC*yfA>n>2fuJuV3)W@PnPO3EQgJI!yVippJz=|5$D$JPby)?IZ$>>XL!wPJ zG}GKpl~|0#bLiB|4n@rAWrvyuqIu70>=J{1Qh1#bmck(t^xAI6$olr{FVHf3L!pbj zA0tNIO?_O84GL82ok0ni;?qzirRR{Y)M|aS4&pN4`6YN^$N>MzA8djHX&cHbp!5wMD6GTImhE=9{@W$~8wIuJ|S132I`4#MY zVp8*zr!&#wWvKJsS|eF8^_VNzI`6HuM1tO1YiR^=J0Oh)_xu&yvJlSXij8aAD+n?9eE(*l+~xn=s(au$fJF$tRejXJKVEFJayvWisPX z1=naRvxQ1yjI5vcuLVL}2of#eQNYDlJqL#lP3J~59|U2qTG+J0jT&T)c!Ff$!wr`d zLo!&M!o3$SzUsl|N{+-j3XAhbJcrrPemTGQxP#UvD!RFtX}D~kJ28U{MV0FinJmXY zhmAz2W+doquV@Jz4)xEWI2jHNEQimbQ8$=mfhPFhMh!@^nBjzk!0jHB7t_67tJUkc zh&hA3oe4oFNq+GMxpbibYFC9rz+k$XP-~u>to8L_)JlmvgpS`P&~m@eVSJ3< zVSi_`had3`g`Iz$6vYz^`0&;~ucuBUg%QRqaA=N8-(7q;1efVFuw5t(u3NmgCOPkR zFH|xxcythi8Lcdj5+uHXsYwF}igb_R6C}Cqq!=U1`v}0EWH5$v2N!|M zGh_+b)GHM%DgQZ&*FCPuX&k>&l8QKSrFn7yOF7hj# zLj$$0J@xRS)8vVa|C((p$w%a zncHrYbT-i%MBO~3XNo65B1|kc#6sYuh32EtM8IW=eG>-v%+61-5hRFAijyLq0B9!5 z9Z4ba{ELyrcI_SCBdKd*Y}dTZ?cvF5DP#yB)E|~qZ0rYU(7rnCTMiR8un^B~vV&X} zv*a;^I~xuPw1Xh8eH7m~@Hyn2K_H$W&VU>BL!5&#_l$ugiyA$`+7IydGJ(~BbF%a> zJ&0GMIu@L?z){QdbOgdGaL6_gbv&r&pv*0!#pkjBLD!5HpD6+aT{BudK`^su)ZkS_ zcH{^5OP;SO)a2->ma@iPEs(tChAej>1Zh*v>o9=~d7t=xbgZAjp0s-j+k4DkIAZSSH02G+?it z?`TQ+!%WHhGFdE8cIzB0%@_8bmXOO{4J;g&;xmO3ggygwmE7#FNKyhp_ih|GbUONB z(g>PVYUo<4)8KS=9H4Otb^eRz5Dh`j!#_M$8%QANU7v7CZWPs`Rg31W@Zjk@90v|Q z?RrfXu1@i^qmgurZydPh+3Q$JDnX(Akjf$2zTm`Gz7@@$(GoJnVNA2zN9>dw1jRP6 zE?r@FQrr_LsWh;^utW65Oc`;vRa&xfQ!%bjWbXh&kKC@sWHQd&0LNr%1sX3YU1CBI zW|zHKz{Jat5EPhLcE}Lctf_Urb*D?HJ;`|1oZ3XkCa$s+Uk(X}Ubun0u8Ea)^_SUOgSe*A!g!h&Ex*IYNvK zYKqv&1wVFT@x3o-9t1sV34h38EGYwCg~kiO&wN^8rO435Ej`Vpwn0i%^)!&=2M0~4bM`_p-;oJluU6Q zZ@VOxyOl)HZccN`(#1C;Fa+t?8Eh|-WN7D5(2|{Rg6<|=+ z3R<$o);%m6+?L?>rUBQsf^ky8*@vVRY;R^SQV^xc&@m&;6=nB-DFl6b&q?Xq@GaSH zwsy#6k(PEpwst8==5y$3f|b;K`8177QS^?%uXBdMMT{bah-E~8pi{F-k)Z%V7cq*Y zh#cCju&~~f)p-nY#joUY5cH_h7#ZZU*aa^Rg2FD1Qfi9B@C-hJ9+qs8A`HP)IS2|_ zvh$H-m=Gn~m%DFWn-7AH(t}((u6z)5SId;56rCd6{kKz;#376)h8Wdm5Bn?9ltR#_ ztsVA1ycev=5;?@ZfxL}zCCSj&Yyi|{QN}csyHDf2q)b7FFTu$a zsk!XScrQtY4&yxvfOknMHV%RUvjI|MNDR;3JYl7@*gh@UKv1L|+m9t%N~WM`92~T8 zdpf8vybg{O83xoej)lahC0k0>BN|VUk($ybJ@pp`ve{C?S+s+MVk)@^5WN zG5yPyc1#XgaKm*g`BFg+gO+xbg7b=zbpxSw#8pa3WX^&wG~P3Q8E{Mcb8OW0V(yEJ z@7gwjCwlLXm?_|Sk_EdI!g>HBY(dXkawfj{av=!&0;0E)n_r$Jg1%NgN>TJ0*onmJ z`U)v}Ip9_i$M7~#<7FT2jHsl79Zv6+goI#CF**8&_{6)UltU-MV22aMDW!1;aCa`i zy{~6r8<1|`X`|P13z5Jym;e5ppPD2(m~=QZa=;@fK)hj5^K+<8Pu6TPxJrgKay}WZ zIFMGujHsxn5admCvV47j90og_CJ*+_?2aAW%w+ z9qLs@mycQ(~OL2X)^_iZj?WJ3(T zNw|k1M>^PbG>4xVlvODv0)1K<*cz1Su+HDpPcmNwVa*0*%mJkkG^of%s}34&H&SMN z46NA%uCNA={xwpLku@;sSd#*a0z1S#1B9ON8#jV5Uky&Gh~_SETWtVIAf8uvLmFfI zBZ)}~-j()9u{|_8A7Rlu4Z&qPutNmL8%|7`U=afBsflF$tlB^}+&Kg%bc;KOgeh5y zUay6HvbE{ZrHux2ShO8BmPNyvEi%|HMy-O#kRazk?Zur;N-mp?78pipbD$hTv%^~> zAii(&{B=coE8*l|FdYVNgD#=so(~-9ir`roitNDjMX((D>570;D3uG>1))!-h8!nu z=|O@4?>VF?wQ6rp_t4rTwC&iVZP$)@JO~y!47``A{Qy&kX2V$MC8aL=JS0JiUiNl6 z0%0;n78l>6XI@r}a+!@E;k+KNjiCuf6(Dxyf?vnlM7`d`(z4eMPP%g|cl-xrg`O%) z&eI^?JUD-0Utog8CK>h(TFZL(g_xTy_{` zk6CtUMX~W6X8~MOapt-~6i*@GxBju9QaZrz&VN_NLh*cPJ^oK!dCzy(0DPrv{{f!} zg(qA)9R=a_W~j7N`QiL$h;My`LH!B={5Gc4&Q1U~m21_|0B}vs1U8oWJf?Qh%4N&S zi-nvUdS+QSdCG<*x1LN$lAnI^DAzo?y)E?W_*;gF@f*rLd9rg!k&E)VymYTf_^D+p zRApGk>U28`2!|Aj2l&ZqYxJowxzS#b#@e##S2>8lH z#?)$zaNh9wemnIL-@KNan|+6TeV#b62*4%3maf$T;M}VtW)@+-Zx?fGV2>m6(5`Lz zw!S`IUVmG;u9uD+lAj4Fa3b8+%T}tz<$7Pe^}>eoAt6kUNtNV_S9QNOw#q8X5A)N{ zHLo?k`4hkwJAQwVq8$+MvRzyHwFX%6pP!1p;t1f~SWh3F*!=pqRsbi}4=MU1FE|A7 zlQ;dndCc{wt5*05rn_3X5;Nt?dHKC^`)$>qb6#Us zZok}n(^u8_+yMUOvXiNY0bWY{8_v2=_$o z8P5;rKSO+1A1Olz0lX;twgSd~yJ*TwIf!riQSWb<6%EGB-hTHZdC6tP|JYOHJ$c1N zsZnqB8-VHlETO2zR%Xu9_p=rc#&nk}a{aM<^(?$*(|DLwQ?Cx$k+~9JZq-XawV-$k z0e>a8$?cN>&+YxRu*w*Sfx=RZSy z7gd+uo&oUew3PUX0IzOrXzm2?!m<5v#hF#_joo!++DUmSQ|IPuz6+H{oQ`8EJvfNz z9^Ca&vMnqs`_$u(Z7|)X&+j=Zw~}(cn#RMdZd&T-)J44j7Bei5*h}#gV)*f~M>zG6{+3Yq@3*Gygm9**@~Dg-&VPpZLWd2wRvBQCbtBp@ z23UH?{jzNkUo&Ro7G~|NGh3BqUz3+E(`eQ3b9d#XO)dJ|N}7b}?o(~J&K5RY^XY;3 znwahq-Rr$0--xxl`*j)*vwE$qRk!^U0j^28KG0)^4#Phxz4kP~Rrk8B>O;TdL+e>} z`NnUD-a=@ccjf zgxBXYR!i}^57TDKOHW;py0Xm~dGTSJzE0~Kjp>f-Jg<{2Z2fnG3-tO9(_Or1+%5TL zUj1GZt=ju_=^Mzp4>7_ zPVp21ULAG4{BeYu2SYAVyf+B(DXT3n^&`MHCe~Xs31I5l;oBk*Uv7lz3ueRAM?YLL zeF%eRHb7`{FdO5e-so6@W$z#leD`*|+K^P%&505`C)6fnM@5 zsm-2rURq3Es70e06Yl+p>2BO6yPqw*_RV#rf0&5rexcXK#_}DcU7JPuV%FR_V;eY2 z4)DZ^U8Ozr58%&^zx8_-z~gaf9-F>t+~*OdyU6fe3+21(-v4>>8+9G@*LHa2y#@vEF>5>2 z`n*%|t`tup;OnKtzBK`0&G&4Nt09HRyB>HFMV@Ac|!x5j|G)Kp5ZuW8XlA+b*VMRu^xjIkF{B!AX zQkS|rFx|E9-F<8eAKK$_rFG*m-GvWrEg|1S+Vhn(9%k+2#`QPUtO_vd&BTr!DV{>W zH;G&ja{^$4@cD)%6wim&b3J~IJT-F#!cIN*{UsDWp>V4n7``pL$_#!u{~6+Ic)Q%d zjR0T$qOxix^1%$Pehgr={80Ix%$6Q=7j7Ha;tct(P`ROo(|DM*2b-?yK3EPg`tPyB4^cdYfN%EgiSBBE4OeHh zDoycxXgz1?o>Hge0E9EfT>vS&Ux+`PaPu{QJt9W5J>eee`2_JbelJD;CBUc-3+pdK zK6cl8I1TZaYClM3wzBVk)4Jvad01!Lmj1JDFb}VN^4^_Wq<&Z3_IHde-1gDPJ(V_M zx(i+SXQ6x_(PP(WOw78n)yI5jstmBo(dLys^Lt&{8aYFI`~l%Ff7<%mtMofQw4R&M zsq4JXUn3kPD^oxyd_%XD-(mRUxcYth;rwTaPd31KI~ri+n(KQ%0$5`vd-^-X@5rqG z53^1GW?Y|Lt>t0rVSn^|l+65D=!JF%lI~)Dzg+1`dt3Oaq$Q(QeNWOoar7SffhOdl z@i6O>Rvo=FO#^U6tNYU#il-3ps}IZ{69I5}k?%KtPVszbJ@;ee$b$duL%3wkK#%?y z;%|P^{w~0>_4*a{=oj+oL3}^0{Af`&z-1FQ4}o+it-e`IRSn?E#;-Ry%WNMUcDLY@ z3-YjW)pj*`w3PWnIVmY5z6hrKrQ=WP+rl5+ANz`R1W9*$y+wYIr@z;|H)+?z_iHoj z-kZ>T+At4(06%%hxuk&rKPdk0n{V*>qtDiJd;S?wV)e%e7gx4kp?C@rfA#3YMF5}F z))u`@@!lZB_u&r%mo5hQ{=^ac${`32?{u22+YqJcWSY)7`ZH zI>7DUZa8p+;`z{ezTnGSzMn80;r&Iueh~<-68%wl7=%~JJUgTYKb-#z@$G8$(uPX_ zx796i<89#P8p@v?Q}D;hla{KHMu0+h1+pj%snt-S1x9_~w?vnBV1&56`qk z#ND6J?#)gl-Jh>|EI&$y-#uwO%=*fA(ngK%3-HFph6oH=li) z;`z{e{>2{-9{cn~fJ+MNuN4vqU&4H`yBmh*Y^%SKAI^V<_-^&wJF+&wYh5BLh9e)@ zYuesN{2#4vbY*ro?Y=MW_g(U^%e!`0J#&}2vwYV_Th^b#beH|)#%s2SfuD~ZUSl0e z_mH>e$bTdKwK+5%W_|qbD@{!w0opfpf!E{IF=dalXY8JO;fp>H0vMZCB9uO|E3lXN#3X_22G{Zxk@X4d~_ z@}0F!y8--g#I$G+dmrGZJd8a*0N@1Eydf98@o7EZ?@6`zPg)~PI#l~H|2@RN-~P%e zgqvO**p460e}?!bo}Uul5#acni)vN}I5niN;x~Yk<~%G^mD%&lz9Um6E|iDY{_@;E zY)9tiroZ0);{6Sn?(iA8(`*rQET2X`FpzYYf6z{T5=Ub(CK?a3e)C~@zZz2r%Q$8L9aNqXZ>H!qbht~7jui{TGZ-X$q+7b_^NyIOy7_tQ6)g$3$`QP)O zA-;o+ON6cl_-niKo3A1t2|r{{NBnKunKzidUu`{CWKs)xc&jl>ixeKf+<0f-)xHz^ zV!Fe&SD0&ySlnk&u?pib-39*c*I#~$^v~DQc$oE9r@y7#cLU&*JIil(pm+)aKOT zy{~Q!rMc!=rf~t-@T>IC;6PX*CtI4bDRH+q%=`Ouv=I6GEb9>ok>D5TOt@Q@WNppTS zp2o*mdmh@|Hfkoo{i>u!P#fGYz#l$TqF+OR`*hzla*W&E-@^P{QN7%js&x^zuNi$u zC_HCQjaM*yd3YuN^9I(Q2R2MP`2)ne_d>zgEy%~g(VhMTcyRi+pPyv*WlXKKDrbT` z{O$H5Rvt`et{9s4J}~+>lJ3`5{a}l@SE2c+5o8=%VBwFu;~nFv6^QY;&M&` z{J7YrUN(xS5b(2eN}XE|a7Kl-+tyP&A6n0E8FTqaxk&)uEnoDK=ePv$35~wfV)(e3 zpBLtT&wqyaK1p4$>l=U{P3)7cLq3Kid{PPE%$RyiBabYTzqIqaVQ?9f+?b$Xu%KYm^mVDeOh#vu>N z2jKgKWVRasu;*;7vI{wlRW%C9Ya!URbwt4T~K4!i#)4(jS0gVLA&vztag&o(ayVs*PLE)Y)1&M z`g-ZoqbQyat><&M|Cw-e9l#y6UM+WAAbdyvosa8b_>A(uZ{vsapCP_#BM%&?2k_-3 zoyIQ*SaZtDD=Gu5{?_d`7cu+C7nw1z=*RN#h2G6=cRQ;8S556M(^OE~2b+$#y zE?=&8a|KEFJF&0I&y#`4Z!{jpy07H0A)S5)SiFM#$O4L|5b&kO+}W@JV6o`B{qIsd zA6n1fTR*sR?=*z}Iiz&-59oV{KhPp`89-~PAuT=T73g~uLVTh6Wp$SUeBtX2H|roD zZ{Ilu^50h?`NxYxnf(_PqrM!$$ivrUW%NJ#3UhJD<@4`U9)40*O$JTz+5=GykE*EZ;*5s53OK}jES|(U)qqQJN?`q z`30CQxi+>Mjfb&?>QD5)@D@N#&hp}p`6HKbL+(1-0u0W_SPU3eZq-8##Z>CA!u zw2SLBznMJzaH)?ZOIdA>O}F|7E&@d@TOv+5mtBr~Y&l(*4J!*=tIjlZT&A zTJ?Uv5GHqJ|0zA&eofNdf5&xOWZxQtBR?%b^1H{HF7k_Mc&TYTjIGV3J0+$R1(-Ye z<1LFRoVa&m#%y&C)`dc1;J0|O$mi3t2CfNrP zmJ&O1i4q4@wn#R-Vc`euNV?lly)JEu??t^pc{6Z;wc1to#p2? zRRvhx)a~?36wim&+?G2h_Is%}!qwTEHVcK{IF{QP!gtlV+$4z~&VPpZYJM@ag#n-} zw~3`Q^08oKQh$K8K2|L%&K#Wg+i#JBbLHWWkAKnf$ZY2P%4*R$RW4$Fhi?3~g)MUM z4}DjTd`!~artw?yOECIA+i@w4hq0M=M1K7jlymczA1kc=gW@R!{LYuotIGh~e8Mup zVV|Pm6k2m*=WO^V<4=S)3XfA!JcWS&wAyPo<{=!=X8-pT?+rqH+g?x?iUYVQ;Y_;+ z$Vb?^k#!OOByMLl=HTBkU(Ge^hDKW1Qv^-v93SvpXTY*XqY*Yf(HOT65n_jVe+{gRtL;s2M`x z6F2{O2gA3msXUV(&VPpZ`kcH_u>rv3SH>ld1K9uKy+(Hd_C3)#bQk|?-+o)&J2dhz zz@C?4`tG233IX3Yu}|eO7+&FnrVhPN!zr}p4oW%x*8n|0Rnc`{9;SE-5&zmdUsp$H z(OveKm!a=b2=Vo5IlBk))qUK8au|Q#wh{NX1MJ;u-sxk^p+^%x>sRM@c|`pQpH6JL ziaA##Eu%%rZ%Dc~{5{MT`K2QHLFkVp-PPg_$*+>>@3v0*hYrg(PTI%V4(&gaI(-<$ zQwaDYaRuXE!0?FQmX)D+KD6d$RquODu>j$_;~umY3SU*X@l6cBTW0A-emMUb;yc*r z)vHYb{yJ(#MnmMIRr$*|03K>;E6|5IoLDz>Z>mNf(c*}^w{E_*er;_joz9 z=L;9qdiT(BlJ1vg?3G_71JMgK9>#Y4#QhBI7J$Ef95d(@il-3pC;v)65B2N%(Vg2L zSty_0j?_A zHup{D@U35ZN8GO^k4W5{to-69=FA_VHzzJ6{qy40D{Qqz{yej!b;~G{?$WIm%C8f9 zY6BV%v!Ud!dlMFQ1-P?AugK%rTV$=BF^PiEQYs(ka#&|14-`s$F=&rtA1K_T~Mc!M<9FbrAt(*Cd zJfiQQA11!hojLRMO>?hmTQR?j-HrXk7WoU4esN_jlI~)oa^%-{;usNj-BKVDFu!_73Lrf#H9h zt5+PuGot_SZg0t5G?|XK>JVINec*2}b%$d;fue?=z63Opb?RVQE z&;GITvwOsz8uIs8XpMe!73_)7}!pt$d9PK835AmNF`@e+nXp0X#Yy?wP%W=B*2kIkc^PiEQEe~_H0<3<&|0inXqf_q+-vF$YG;~`9 z=4iIHbiFqH(*1k0cjdRxyf#@$ z<6$