...
 
Commits (2)
......@@ -215,7 +215,10 @@
"mprealsupport": "cpp",
"LynxADDarcyVelocityAux.C": "cpp",
"LynxADPorosityAux.C": "cpp",
"LynxADHeatFluxAux.C": "cpp"
"LynxADHeatFluxAux.C": "cpp",
"LynxADDamageBase.C": "cpp",
"LynxADDamageModelBase.C": "cpp",
"LynxADLyakhovskyDamage.C": "cpp"
},
"C_Cpp.errorSquiggles": "Disabled",
"python.pythonPath": "/opt/moose/miniconda/bin/python"
......
......@@ -26,8 +26,8 @@ protected:
const unsigned int _component;
const VariableGradient & _grad_pf;
const ADMaterialProperty<Real> & _fluid_mobility;
const MaterialProperty<Real> & _fluid_mobility;
const bool _coupled_grav;
const ADMaterialProperty<RealVectorValue> * _gravity;
const MaterialProperty<RealVectorValue> * _gravity;
const ADMaterialProperty<Real> * _rho_f;
};
\ No newline at end of file
......@@ -34,7 +34,7 @@ protected:
// velocities
unsigned int _nvel;
std::vector<const ADVariableValue *> _vel;
std::vector<const VariableValue *> _vel;
std::vector<const VariableValue *> _vel_old;
std::vector<const VariableValue *> _vel_older;
......
......@@ -26,13 +26,13 @@ protected:
const Real _coeff_Hs;
unsigned int _nvel;
std::vector<const ADVariableValue *> _vel;
std::vector<const VariableValue *> _vel;
const ADVariableGradient & _grad_pressure;
const ADMaterialProperty<Real> & _rhoC_b;
const MaterialProperty<Real> & _radiogenic_heat;
const bool _has_inelastic_heat_mat;
const ADMaterialProperty<Real> * _radiogenic_heat;
const ADMaterialProperty<Real> * _inelastic_heat_mat;
const bool _coupled_inelastic_heat;
const ADVariableValue & _inelastic_heat;
const ADMaterialProperty<Real> & _thermal_exp;
const VariableValue & _inelastic_heat;
const MaterialProperty<Real> & _thermal_exp;
};
\ No newline at end of file
......@@ -24,8 +24,8 @@ public:
protected:
virtual ADReal computeQpResidual() override;
const ADMaterialProperty<Real> & _fluid_mobility;
const MaterialProperty<Real> & _fluid_mobility;
const bool _coupled_grav;
const ADMaterialProperty<RealVectorValue> * _gravity;
const MaterialProperty<RealVectorValue> * _gravity;
const ADMaterialProperty<Real> * _rho_f;
};
\ No newline at end of file
......@@ -25,5 +25,5 @@ protected:
virtual ADReal computeQpResidual() override;
const ADMaterialProperty<Real> & _bulk_density;
const ADMaterialProperty<RealVectorValue> & _gravity;
const MaterialProperty<RealVectorValue> & _gravity;
};
\ No newline at end of file
......@@ -30,9 +30,9 @@ protected:
const bool _vol_locking_correction;
const ADMaterialProperty<RankTwoTensor> & _stress;
const bool _coupled_pf;
const ADMaterialProperty<Real> * _biot;
const MaterialProperty<Real> * _biot;
const bool _coupled_grav;
const ADMaterialProperty<RealVectorValue> * _gravity;
const MaterialProperty<RealVectorValue> * _gravity;
const ADMaterialProperty<Real> * _rho_b;
std::vector<ADReal> _avg_grad_test;
};
\ No newline at end of file
......@@ -23,13 +23,13 @@ public:
void setQp(unsigned int qp);
virtual void creepUpdate(ADRankTwoTensor & stress_dev,
const ADReal & pressure,
const ADReal & G,
const Real & G,
ADRankTwoTensor & elastic_strain_incr);
void resetQpProperties() final {}
void resetProperties() final {}
protected:
virtual ADReal viscousIncrement(const ADReal & pressure, const ADReal & eqv_stress, const ADReal & G);
virtual ADReal viscousIncrement(const ADReal & pressure, const ADReal & eqv_stress, const Real & G);
virtual ADReal computeQpEffectiveViscosity(const ADReal & pressure);
virtual ADReal computeQpOneOnDiffViscosity(const ADReal A);
virtual ADReal computeQpOneOnDislViscosity(const ADReal A, const ADReal n, const ADReal eII);
......@@ -70,15 +70,15 @@ protected:
// Creep effective parameters
ADReal _A_diff;
ADReal _E_diff;
ADReal _V_diff;
Real _E_diff;
Real _V_diff;
ADReal _A_disl;
ADReal _n_disl;
ADReal _E_disl;
ADReal _V_disl;
Real _n_disl;
Real _E_disl;
Real _V_disl;
// For iterative update
ADReal _tau_II_tr;
ADReal _G;
Real _G;
};
\ No newline at end of file
/******************************************************************************/
/* This file is part of */
/* LYNX, a MOOSE-based application */
/* Lithosphere dYnamic Numerical toolboX */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* Licensed under GNU General Public License 3, */
/* please see LICENSE for details */
/* or http://www.gnu.org/licenses/gpl.html */
/******************************************************************************/
#pragma once
#include "LynxADElasticDeformation.h"
#include "LynxADDamageModelBase.h"
class LynxADDamageDeformation : public LynxADElasticDeformation
{
public:
static InputParameters validParams();
LynxADDamageDeformation(const InputParameters & parameters);
protected:
virtual void initializeQpDeformation() override;
virtual void computeQpStress() override;
// Damage rheology
LynxADDamageModelBase * _damage_model;
// Strain properties
ADMaterialProperty<RankTwoTensor> & _elastic_strain;
const MaterialProperty<RankTwoTensor> & _elastic_strain_old;
// Elasticity tensor
ADRankFourTensor _Cijkl;
};
\ No newline at end of file
......@@ -58,7 +58,7 @@ protected:
// Strain properties
ADMaterialProperty<RankTwoTensor> & _strain_increment;
ADMaterialProperty<RankTwoTensor> & _spin_increment;
const ADMaterialProperty<Real> * _thermal_exp;
const MaterialProperty<Real> * _thermal_exp;
// Stress properties
ADMaterialProperty<RankTwoTensor> & _stress;
......
......@@ -24,14 +24,14 @@ public:
protected:
virtual void computeQpGravity();
const ADVariableValue & _porosity;
const VariableValue & _porosity;
bool _has_gravity;
Real _g;
const std::vector<Real> _fluid_density;
const std::vector<Real> _solid_density;
ADMaterialProperty<RealVectorValue> & _gravity;
MaterialProperty<RealVectorValue> & _gravity;
ADMaterialProperty<Real> & _rho_f;
ADMaterialProperty<Real> & _rho_s;
ADMaterialProperty<Real> & _rho_b;
......
......@@ -44,8 +44,8 @@ protected:
// Elastic properties
const VariableValue & _plith_old;
ADMaterialProperty<RankTwoTensor> & _elastic_strain_incr;
ADMaterialProperty<Real> & _K;
ADMaterialProperty<Real> & _G;
MaterialProperty<Real> & _K;
MaterialProperty<Real> & _G;
const MaterialProperty<RankTwoTensor> & _stress_old;
const ADMaterialProperty<RankTwoTensor> * _viscous_strain_incr;
const ADMaterialProperty<RankTwoTensor> * _plastic_strain_incr;
......
......@@ -31,22 +31,22 @@ protected:
virtual void computeQpPermeability() = 0;
virtual void computeQpFluidViscosity() = 0;
const ADVariableValue & _porosity;
const VariableValue & _porosity;
const bool _coupled_mech;
const ADMaterialProperty<Real> * _K;
const MaterialProperty<Real> * _K;
const ADMaterialProperty<RankTwoTensor> * _strain_increment;
const bool _has_viscous;
const ADMaterialProperty<RankTwoTensor> * _viscous_strain_incr;
const bool _has_plastic;
const ADMaterialProperty<RankTwoTensor> * _plastic_strain_incr;
ADMaterialProperty<Real> & _biot;
ADMaterialProperty<Real> & _C_d;
ADMaterialProperty<Real> & _C_biot;
ADMaterialProperty<Real> & _fluid_mobility;
MaterialProperty<Real> & _biot;
MaterialProperty<Real> & _C_d;
MaterialProperty<Real> & _C_biot;
MaterialProperty<Real> & _fluid_mobility;
ADMaterialProperty<Real> & _poro_mech;
std::vector<ADReal> _C_f;
std::vector<ADReal> _C_s;
std::vector<ADReal> _k;
std::vector<ADReal> _eta_f;
std::vector<Real> _C_f;
std::vector<Real> _C_s;
std::vector<Real> _k;
std::vector<Real> _eta_f;
};
\ No newline at end of file
......@@ -24,13 +24,13 @@ public:
const std::vector<T> & getLynxParam(const std::string & name) const;
protected:
virtual ADReal averageProperty(const std::vector<Real> & properties);
virtual ADReal arithmetic_average(const std::vector<Real> & properties);
virtual ADReal harmonic_average(const std::vector<Real> & properties);
virtual ADReal max_average(const std::vector<Real> & properties);
virtual Real averageProperty(const std::vector<Real> & properties);
virtual Real arithmetic_average(const std::vector<Real> & properties);
virtual Real harmonic_average(const std::vector<Real> & properties);
virtual Real max_average(const std::vector<Real> & properties);
const bool _has_compositional_phases;
const unsigned int _n_composition;
const unsigned int _average_type;
std::vector<const ADVariableValue *> _compositional_phases;
std::vector<const VariableValue *> _compositional_phases;
};
\ No newline at end of file
......@@ -23,16 +23,16 @@ public:
void setQp(unsigned int qp);
virtual void plasticUpdate(ADRankTwoTensor & stress_dev,
ADReal & pressure,
const ADReal & G,
const ADReal & K,
const Real & G,
const Real & K,
ADRankTwoTensor & elastic_strain_incr);
void resetQpProperties() final {}
void resetProperties() final {}
protected:
virtual void initQpStatefulProperties() override;
virtual ADReal plasticIncrement(const ADReal & eqv_stress, const ADReal & pressure, const ADReal G, const ADReal K);
virtual void initPlasticParameters(const ADReal & pressure, const ADReal & K);
virtual ADReal plasticIncrement(const ADReal & eqv_stress, const ADReal & pressure, const Real G, const Real K);
virtual void initPlasticParameters(const ADReal & pressure, const Real & K);
virtual ADReal plasticYieldFunction(const ADReal & eqv_stress, const ADReal & pressure);
// Plastic parameters
......@@ -53,14 +53,14 @@ protected:
const MaterialProperty<Real> * _intnl_old;
// Plastic effective parameters
ADReal _alpha_0;
ADReal _alpha_res;
ADReal _k_0;
ADReal _k_res;
ADReal _beta;
ADReal _intnl_0;
ADReal _intnl_lim;
ADReal _one_on_eta;
Real _alpha_0;
Real _alpha_res;
Real _k_0;
Real _k_res;
Real _beta;
Real _intnl_0;
Real _intnl_lim;
Real _one_on_eta;
ADReal _alpha;
ADReal _k;
ADReal _H;
......
......@@ -55,11 +55,11 @@ protected:
ADMaterialProperty<Real> & _eta_eff;
ADReal _A_diff;
ADReal _E_diff;
ADReal _V_diff;
Real _E_diff;
Real _V_diff;
ADReal _A_disl;
ADReal _n_disl;
ADReal _E_disl;
ADReal _V_disl;
Real _n_disl;
Real _E_disl;
Real _V_disl;
};
\ No newline at end of file
......@@ -30,10 +30,11 @@ protected:
virtual void computeQpHeatCap() = 0;
virtual void computeQpThermalCond() = 0;
virtual void computeQpThermalExp() = 0;
virtual ADReal computeMixtureProperty(const ADReal fluid_prop, const ADReal solid_prop);
virtual Real computeMixtureProperty(const Real fluid_prop, const Real solid_prop);
virtual ADReal computeADMixtureProperty(const ADReal fluid_prop, const ADReal solid_prop);
// const bool _coupled_porosity;
const ADVariableValue & _porosity;
const VariableValue & _porosity;
const std::vector<Real> _heat_source;
const bool _coupled_dens;
......@@ -42,13 +43,13 @@ protected:
ADMaterialProperty<Real> & _thermal_diff;
ADMaterialProperty<Real> & _rhoC_b;
ADMaterialProperty<Real> & _rhoC_f;
ADMaterialProperty<Real> & _thermal_exp;
ADMaterialProperty<Real> & _radiogenic_heat;
MaterialProperty<Real> & _thermal_exp;
MaterialProperty<Real> & _radiogenic_heat;
std::vector<ADReal> _c_f;
std::vector<ADReal> _c_s;
std::vector<ADReal> _lambda_f;
std::vector<ADReal> _lambda_s;
std::vector<ADReal> _beta_f;
std::vector<ADReal> _beta_s;
std::vector<Real> _c_f;
std::vector<Real> _c_s;
std::vector<Real> _lambda_f;
std::vector<Real> _lambda_s;
std::vector<Real> _beta_f;
std::vector<Real> _beta_s;
};
\ No newline at end of file
......@@ -33,9 +33,9 @@ LynxADDarcyVelocityAux::LynxADDarcyVelocityAux(const InputParameters & parameter
: AuxKernel(parameters),
_component(getParam<unsigned int>("component")),
_grad_pf(coupledGradient("fluid_pressure")),
_fluid_mobility(getADMaterialProperty<Real>("fluid_mobility")),
_coupled_grav(hasADMaterialProperty<Real>("fluid_density")),
_gravity(_coupled_grav ? &getADMaterialProperty<RealVectorValue>("gravity_vector") : nullptr),
_fluid_mobility(getMaterialProperty<Real>("fluid_mobility")),
_coupled_grav(hasMaterialProperty<RealVectorValue>("gravity_vector")),
_gravity(_coupled_grav ? &getMaterialProperty<RealVectorValue>("gravity_vector") : nullptr),
_rho_f(_coupled_grav ? &getADMaterialProperty<Real>("fluid_density") : nullptr)
{
}
......@@ -44,9 +44,9 @@ Real
LynxADDarcyVelocityAux::computeValue()
{
RealVectorValue grav_term = _coupled_grav ? MetaPhysicL::raw_value(-(*_rho_f)[_qp]) *
MetaPhysicL::raw_value((*_gravity)[_qp])
(*_gravity)[_qp]
: RealVectorValue();
return -MetaPhysicL::raw_value(_fluid_mobility[_qp]) *
return -_fluid_mobility[_qp] *
(_grad_pf[_qp](_component) + grav_term(_component));
}
......@@ -31,10 +31,10 @@ LynxADPorosityAux::LynxADPorosityAux(const InputParameters & parameters)
_biot(getADMaterialProperty<Real>("biot_coefficient")),
_C_d(getADMaterialProperty<Real>("bulk_compressibility")),
_strain_increment(getADMaterialProperty<RankTwoTensor>("strain_increment")),
_has_viscous(hasADMaterialProperty<RankTwoTensor>("viscous_strain_increment")),
_has_viscous(hasMaterialProperty<RankTwoTensor>("viscous_strain_increment")),
_viscous_strain_incr(
_has_viscous ? &getADMaterialProperty<RankTwoTensor>("viscous_strain_increment") : nullptr),
_has_plastic(hasADMaterialProperty<RankTwoTensor>("plastic_strain_increment")),
_has_plastic(hasMaterialProperty<RankTwoTensor>("plastic_strain_increment")),
_plastic_strain_incr(
_has_viscous ? &getADMaterialProperty<RankTwoTensor>("plastic_strain_increment") : nullptr)
{
......
......@@ -80,13 +80,13 @@ LynxADAdvectionBase::LynxADAdvectionBase(const InputParameters & parameters)
for (unsigned i = 0; i < _nvel; ++i)
{
_vel[i] = &adCoupledValue("velocities", i);
_vel[i] = &coupledValue("velocities", i);
_vel_old[i] = &coupledValueOld("velocities", i);
_vel_older[i] = &coupledValueOlder("velocities", i);
}
for (unsigned i = _nvel; i < 3; ++i)
{
_vel[i] = &adZeroValue();
_vel[i] = &_zero;
_vel_old[i] = &_zero;
_vel_older[i] = &_zero;
}
......@@ -128,7 +128,7 @@ LynxADAdvectionBase::precalculateResidual()
ADReal
LynxADAdvectionBase::computeQpResidual()
{
ADRealVectorValue u((*_vel[0])[_qp], (*_vel[1])[_qp], (*_vel[2])[_qp]);
RealVectorValue u((*_vel[0])[_qp], (*_vel[1])[_qp], (*_vel[2])[_qp]);
ADReal res = u * _grad_u[_qp] * _test[_i][_qp];
Real a2 = this->_t_step > 1 ? _dt / this->_dt_old : 0.0;
Real a1 = 1.0 + a2;
......
......@@ -36,7 +36,7 @@ LynxADAdvectionComposition::computeArtificialViscosity()
computeEntropyResidual();
ADReal max_residual = 0.0;
ADReal max_velocity = 0.0;
Real max_velocity = 0.0;
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
{
......@@ -49,8 +49,8 @@ LynxADAdvectionComposition::computeArtificialViscosity()
max_velocity = std::max(u.norm(), max_velocity);
}
ADReal max_viscosity = _beta_stabilization * max_velocity * diameter;
ADReal entropy_variation =
Real max_viscosity = _beta_stabilization * max_velocity * diameter;
Real entropy_variation =
std::max(_pp_max_entropy - _pp_avg_entropy, _pp_avg_entropy - _pp_min_entropy);
if (this->_t_step <= 1 || std::abs(entropy_variation) < 1e-50)
......
......@@ -35,7 +35,7 @@ LynxADAdvectionTemperature::LynxADAdvectionTemperature(const InputParameters & p
_grad_pressure(isCoupled("pressure") ? adCoupledGradient("pressure") : adZeroGradient()),
_thermal_diff(getADMaterialProperty<Real>("thermal_diffusivity")),
_rhoC(getADMaterialProperty<Real>("bulk_specific_heat")),
_has_inelastic_heat_mat(hasADMaterialProperty<Real>("inelastic_heat")),
_has_inelastic_heat_mat(hasMaterialProperty<Real>("inelastic_heat")),
_radiogenic_heat(_has_inelastic_heat_mat
? &getADMaterialProperty<Real>("radiogenic_heat_production")
: nullptr),
......
......@@ -38,26 +38,25 @@ LynxADHeatSources::LynxADHeatSources(const InputParameters & parameters)
_vel(3),
_grad_pressure(isCoupled("pressure") ? adCoupledGradient("pressure") : adZeroGradient()),
_rhoC_b(getADMaterialProperty<Real>("bulk_specific_heat")),
_has_inelastic_heat_mat(hasADMaterialProperty<Real>("inelastic_heat")),
_radiogenic_heat(_has_inelastic_heat_mat
? &getADMaterialProperty<Real>("radiogenic_heat_production")
: nullptr),
_radiogenic_heat(getMaterialProperty<Real>("radiogenic_heat_production")),
// Here she need hasADMaterialProperty...
_has_inelastic_heat_mat(hasMaterialProperty<Real>("inelastic_heat")),
_inelastic_heat_mat(_has_inelastic_heat_mat ? &getADMaterialProperty<Real>("inelastic_heat")
: nullptr),
_coupled_inelastic_heat(isCoupled("inelastic_heat")),
_inelastic_heat(_coupled_inelastic_heat ? adCoupledValue("inelastic_heat") : adZeroValue()),
_thermal_exp(getADMaterialProperty<Real>("thermal_expansion_coefficient"))
_inelastic_heat(_coupled_inelastic_heat ? coupledValue("inelastic_heat") : _zero),
_thermal_exp(getMaterialProperty<Real>("thermal_expansion_coefficient"))
{
for (unsigned i = 0; i < _nvel; ++i)
_vel[i] = &adCoupledValue("velocities", i);
_vel[i] = &coupledValue("velocities", i);
for (unsigned i = _nvel; i < 3; ++i)
_vel[i] = &adZeroValue();
_vel[i] = &_zero;
}
ADReal
LynxADHeatSources::computeQpResidual()
{
ADReal Hr = _has_inelastic_heat_mat ? (*_radiogenic_heat)[_qp] : 0.0;
Real Hr = _radiogenic_heat[_qp];
ADReal Hs = _coeff_Hs;
if (_coupled_inelastic_heat)
Hs *= _inelastic_heat[_qp];
......@@ -65,7 +64,7 @@ LynxADHeatSources::computeQpResidual()
Hs *= (*_inelastic_heat_mat)[_qp];
else
Hs *= 0.0;
ADRealVectorValue vel((*_vel[0])[_qp], (*_vel[1])[_qp], (*_vel[2])[_qp]);
RealVectorValue vel((*_vel[0])[_qp], (*_vel[1])[_qp], (*_vel[2])[_qp]);
ADReal Ha = _thermal_exp[_qp] * _u[_qp] * vel * _grad_pressure[_qp];
ADReal heat_sources = Hr + Hs + Ha;
......
......@@ -25,9 +25,9 @@ LynxADHydroDarcy::validParams()
LynxADHydroDarcy::LynxADHydroDarcy(const InputParameters & parameters)
: ADKernel(parameters),
_fluid_mobility(getADMaterialProperty<Real>("fluid_mobility")),
_coupled_grav(hasADMaterialProperty<Real>("fluid_density")),
_gravity(_coupled_grav ? &getADMaterialProperty<RealVectorValue>("gravity_vector") : nullptr),
_fluid_mobility(getMaterialProperty<Real>("fluid_mobility")),
_coupled_grav(hasMaterialProperty<RealVectorValue>("gravity_vector")),
_gravity(_coupled_grav ? &getMaterialProperty<RealVectorValue>("gravity_vector") : nullptr),
_rho_f(_coupled_grav ? &getADMaterialProperty<Real>("fluid_density") : nullptr)
{
}
......
......@@ -27,7 +27,7 @@ LynxADPressureLoad::validParams()
LynxADPressureLoad::LynxADPressureLoad(const InputParameters & parameters)
: ADKernel(parameters),
_bulk_density(getADMaterialProperty<Real>("reference_bulk_density")),
_gravity(getADMaterialProperty<RealVectorValue>("gravity_vector"))
_gravity(getMaterialProperty<RealVectorValue>("gravity_vector"))
{
}
......
......@@ -38,10 +38,10 @@ LynxADSolidMomentum::LynxADSolidMomentum(const InputParameters & parameters)
_component(getParam<unsigned int>("component")),
_vol_locking_correction(getParam<bool>("volumetric_locking_correction")),
_stress(getADMaterialProperty<RankTwoTensor>("stress")),
_coupled_pf(hasADMaterialProperty<Real>("biot_coefficient")),
_biot(_coupled_pf ? &getADMaterialProperty<Real>("biot_coefficient") : nullptr),
_coupled_grav(hasADMaterialProperty<Real>("bulk_density")),
_gravity(_coupled_grav ? &getADMaterialProperty<RealVectorValue>("gravity_vector") : nullptr),
_coupled_pf(hasMaterialProperty<Real>("biot_coefficient")),
_biot(_coupled_pf ? &getMaterialProperty<Real>("biot_coefficient") : nullptr),
_coupled_grav(hasMaterialProperty<RealVectorValue>("gravity_vector")),
_gravity(_coupled_grav ? &getMaterialProperty<RealVectorValue>("gravity_vector") : nullptr),
_rho_b(_coupled_grav ? &getADMaterialProperty<Real>("bulk_density") : nullptr),
_avg_grad_test()
{
......
......@@ -117,7 +117,7 @@ LynxADCreepModel::setQp(unsigned int qp)
void
LynxADCreepModel::creepUpdate(ADRankTwoTensor & stress_dev,
const ADReal & pressure,
const ADReal & G,
const Real & G,
ADRankTwoTensor & elastic_strain_incr)
{
if (G != 0.0) // Visco elastic update
......@@ -145,7 +145,7 @@ LynxADCreepModel::creepUpdate(ADRankTwoTensor & stress_dev,
ADReal
LynxADCreepModel::viscousIncrement(const ADReal & pressure,
const ADReal & tau_II_tr,
const ADReal & G)
const Real & G)
{
if (tau_II_tr == 0.0)
{
......
/******************************************************************************/
/* This file is part of */
/* LYNX, a MOOSE-based application */
/* Lithosphere dYnamic Numerical toolboX */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* Licensed under GNU General Public License 3, */
/* please see LICENSE for details */
/* or http://www.gnu.org/licenses/gpl.html */
/******************************************************************************/
#include "LynxADDamageDeformation.h"
#include "ElasticityTensorTools.h"
#include "LynxADCreepModel.h"
#include "LynxADPlasticModel.h"
registerMooseObject("LynxApp", LynxADDamageDeformation);
InputParameters
LynxADDamageDeformation::validParams()
{
InputParameters params = LynxADElasticDeformation::validParams();
params.addClassDescription("Class calculating the deformation of a material following a "
"damage viscoelasto-(visco)plastic rheology.");
params.addRequiredParam<MaterialName>(
"damage_model",
"The material object to use for the damage rheology.");
return params;
}
LynxADDamageDeformation::LynxADDamageDeformation(const InputParameters & parameters)
: LynxADElasticDeformation(parameters),
// Damage model
// _damage_model(&getMaterialByName("damage_model")),
// Strain properties
_elastic_strain(declareADProperty<RankTwoTensor>("elastic_strain")),
_elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("elastic_strain"))
{
LynxADDamageModelBase * damage_r =
dynamic_cast<LynxADDamageModelBase *>(&this->getMaterialByName("damage_model"));
_damage_model = damage_r;
}
void
LynxADDamageDeformation::initializeQpDeformation()
{
// Initialize elastic properties
_elastic_strain_incr[_qp] = _strain_increment[_qp];
_K[_qp] = averageProperty(_bulk_modulus);
_G[_qp] = averageProperty(_shear_modulus);
_Cijkl = ElasticityTensorTools::elasticityTensorKandG(_K[_qp], _G[_qp]);
_elastic_strain[_qp] = spinRotation(_elastic_strain_old[_qp]) + _elastic_strain_incr[_qp];
}
void
LynxADDamageDeformation::computeQpStress()
{
ADRankTwoTensor stress_old = spinRotation(_stress_old[_qp]);
ADRankTwoTensor elastic_strain_old = spinRotation(_elastic_strain_old[_qp]);
// Damage elastic guess
_damage_model->setQp(_qp);
_damage_model->elasticGuess(_stress[_qp], stress_old, _Cijkl, elastic_strain_old, _elastic_strain_incr[_qp]);
_stress[_qp].addIa(_plith_old[_qp] - _plith[_qp]);
// Damage - plastic correction
_damage_model->damageUpdate(_stress[_qp], elastic_strain_old);
}
// void
// LynxADDamageDeformation::computeQpThermalSources()
// {
// LynxADDeformationBase::computeQpThermalSources();
// if (_has_plasticity)
// {
// Real I2 = Utility::pow<2>(_elastic_strain[_qp].L2norm());
// Real xi = strainRatio(_elastic_strain[_qp]);
// Real dPsi_dalpha = _damage_plasticity->_gamma * I2 * (xi - _damage_plasticity->_xi0);
// _damage_heat[_qp] = dPsi_dalpha * _damage_rate[_qp];
// }
// else
// _damage_heat[_qp] = 0.0;
// }
......@@ -57,7 +57,7 @@ LynxADDeformationBase::LynxADDeformationBase(const InputParameters & parameters)
_strain_increment(declareADProperty<RankTwoTensor>("strain_increment")),
_spin_increment(declareADProperty<RankTwoTensor>("spin_increment")),
_thermal_exp(_coupled_temp || _coupled_temp_aux
? &getADMaterialProperty<Real>("thermal_expansion_coefficient")
? &getMaterialProperty<Real>("thermal_expansion_coefficient")
: nullptr),
// Stress properties
_stress(declareADProperty<RankTwoTensor>("stress")),
......
......@@ -19,7 +19,7 @@ LynxADDensityBase::validParams()
{
InputParameters params = LynxADMaterialBase::validParams();
params.addClassDescription("Base class for calculating densities and gravity.");
params.addCoupledVar("porosity", "The porosity auxiliary variable.");
params.addCoupledVar("porosity", 0.0, "The porosity auxiliary variable.");
params.addParam<bool>("has_gravity", false, "Model with gravity on?");
params.addParam<Real>("gravity_acceleration", 9.81, "The magnitude of the gravity acceleration.");
params.addParam<std::vector<Real>>("fluid_density", "The fluid density.");
......@@ -29,14 +29,14 @@ LynxADDensityBase::validParams()
LynxADDensityBase::LynxADDensityBase(const InputParameters & parameters)
: LynxADMaterialBase(parameters),
_porosity(isCoupled("porosity") ? adCoupledValue("porosity") : adZeroValue()),
_porosity(coupledValue("porosity")),
_has_gravity(getParam<bool>("has_gravity")),
_g(_has_gravity ? getParam<Real>("gravity_acceleration") : 0.0),
_fluid_density(isParamValid("fluid_density") ? getLynxParam<Real>("fluid_density")
: std::vector<Real>(_n_composition, 0.0)),
_solid_density(isParamValid("solid_density") ? getLynxParam<Real>("solid_density")
: std::vector<Real>(_n_composition, 0.0)),
_gravity(declareADProperty<RealVectorValue>("gravity_vector")),
_gravity(declareProperty<RealVectorValue>("gravity_vector")),
_rho_f(declareADProperty<Real>("fluid_density")),
_rho_s(declareADProperty<Real>("solid_density")),
_rho_b(declareADProperty<Real>("bulk_density")),
......@@ -48,9 +48,9 @@ void
LynxADDensityBase::computeQpGravity()
{
if (_mesh.dimension() == 3)
_gravity[_qp] = ADRealVectorValue(0., 0., -_g);
_gravity[_qp] = RealVectorValue(0., 0., -_g);
else if (_mesh.dimension() == 2)
_gravity[_qp] = ADRealVectorValue(0., -_g, 0.);
_gravity[_qp] = RealVectorValue(0., -_g, 0.);
else if (_mesh.dimension() == 1)
_gravity[_qp] = ADRealVectorValue(-_g, 0., 0.);
_gravity[_qp] = RealVectorValue(-_g, 0., 0.);
}
\ No newline at end of file
......@@ -49,8 +49,8 @@ LynxADElasticDeformation::LynxADElasticDeformation(const InputParameters & param
// Elastic properties
_plith_old(isCoupled("lithostatic_pressure") ? coupledValueOld("lithostatic_pressure") : _zero),
_elastic_strain_incr(declareADProperty<RankTwoTensor>("elastic_strain_increment")),
_K(declareADProperty<Real>("bulk_modulus")),
_G(declareADProperty<Real>("shear_modulus")),
_K(declareProperty<Real>("bulk_modulus")),
_G(declareProperty<Real>("shear_modulus")),
_stress_old(getMaterialPropertyOld<RankTwoTensor>("stress")),
_viscous_strain_incr(
_has_creep ? &getADMaterialProperty<RankTwoTensor>("viscous_strain_increment") : nullptr),
......
......@@ -24,21 +24,21 @@ LynxADHydroBase::validParams()
LynxADHydroBase::LynxADHydroBase(const InputParameters & parameters)
: LynxADMaterialBase(parameters),
_porosity(adCoupledValue("porosity")),
_coupled_mech(hasADMaterialProperty<Real>("bulk_modulus")),
_K(_coupled_mech ? &getADMaterialProperty<Real>("bulk_modulus") : nullptr),
_porosity(coupledValue("porosity")),
_coupled_mech(hasMaterialProperty<Real>("bulk_modulus")),
_K(_coupled_mech ? &getMaterialProperty<Real>("bulk_modulus") : nullptr),
_strain_increment(_coupled_mech ? &getADMaterialProperty<RankTwoTensor>("strain_increment")
: nullptr),
_has_viscous(hasADMaterialProperty<RankTwoTensor>("viscous_strain_increment")),
_has_viscous(hasMaterialProperty<RankTwoTensor>("viscous_strain_increment")),
_viscous_strain_incr(
_has_viscous ? &getADMaterialProperty<RankTwoTensor>("viscous_strain_increment") : nullptr),
_has_plastic(hasADMaterialProperty<RankTwoTensor>("plastic_strain_increment")),
_has_plastic(hasMaterialProperty<RankTwoTensor>("plastic_strain_increment")),
_plastic_strain_incr(
_has_viscous ? &getADMaterialProperty<RankTwoTensor>("plastic_strain_increment") : nullptr),
_biot(declareADProperty<Real>("biot_coefficient")),
_C_d(declareADProperty<Real>("bulk_compressibility")),
_C_biot(declareADProperty<Real>("biot_compressibility")),
_fluid_mobility(declareADProperty<Real>("fluid_mobility")),
_biot(declareProperty<Real>("biot_coefficient")),
_C_d(declareProperty<Real>("bulk_compressibility")),
_C_biot(declareProperty<Real>("biot_compressibility")),
_fluid_mobility(declareProperty<Real>("fluid_mobility")),
_poro_mech(declareADProperty<Real>("poro_mechanical")),
_C_f(_fe_problem.getMaxQps()),
_C_s(_fe_problem.getMaxQps()),
......@@ -70,7 +70,7 @@ LynxADHydroBase::computeQpCompressibilities()
_biot[_qp] -= _C_s[_qp] / _C_d[_qp];
// Pore compressibility
ADReal C_phi = (_biot[_qp] - _porosity[_qp]) * _C_d[_qp];
Real C_phi = (_biot[_qp] - _porosity[_qp]) * _C_d[_qp];
// Biot compressibility
_C_biot[_qp] = _porosity[_qp] * _C_f[_qp] + (1.0 - _biot[_qp]) * C_phi;
......
......@@ -36,7 +36,7 @@ LynxADMaterialBase::LynxADMaterialBase(const InputParameters & parameters)
{
if (_has_compositional_phases)
for (unsigned i = 0; i < _n_composition; ++i)
_compositional_phases[i] = &adCoupledValue("compositional_phases", i);
_compositional_phases[i] = &coupledValue("compositional_phases", i);
}
template <typename T>
......@@ -52,12 +52,12 @@ LynxADMaterialBase::getLynxParam(const std::string & name) const
return prop;
}
ADReal
Real
LynxADMaterialBase::averageProperty(const std::vector<Real> & properties)
{
if (_n_composition == 1)
return properties[0];
ADReal average_property = 0.0;
Real average_property = 0.0;
switch (_average_type)
{
case 0: // ARITHMETIC
......@@ -74,12 +74,12 @@ LynxADMaterialBase::averageProperty(const std::vector<Real> & properties)
return average_property;
}
ADReal
Real
LynxADMaterialBase::arithmetic_average(const std::vector<Real> & properties)
{
ADReal phase = 0;
ADReal sum = 0;
ADReal value = 0;
Real phase = 0;
Real sum = 0;
Real value = 0;
for (unsigned i = 0; i < _n_composition; ++i)
{
phase = std::min(std::max(0.0, (*_compositional_phases[i])[_qp]), 1.0);
......@@ -91,12 +91,12 @@ LynxADMaterialBase::arithmetic_average(const std::vector<Real> & properties)
return value;
}
ADReal
Real
LynxADMaterialBase::harmonic_average(const std::vector<Real> & properties)
{
ADReal phase = 0;
ADReal sum = 0;
ADReal value = 0;
Real phase = 0;
Real sum = 0;
Real value = 0;
for (unsigned i = 0; i < _n_composition; ++i)
{
phase = std::min(std::max(0.0, (*_compositional_phases[i])[_qp]), 1.0);
......@@ -108,11 +108,11 @@ LynxADMaterialBase::harmonic_average(const std::vector<Real> & properties)
return (value != 0) ? 1.0 / value : 0;
}
ADReal
Real
LynxADMaterialBase::max_average(const std::vector<Real> & properties)
{
ADReal max_phase = (*_compositional_phases[0])[_qp];
ADReal value = properties[0];
Real max_phase = (*_compositional_phases[0])[_qp];
Real value = properties[0];
for (unsigned i = 1; i < _n_composition; ++i)
if ((*_compositional_phases[i])[_qp] > max_phase)
{
......
......@@ -100,8 +100,8 @@ LynxADPlasticModel::setQp(unsigned int qp)
void
LynxADPlasticModel::plasticUpdate(ADRankTwoTensor & stress_dev,
ADReal & pressure,
const ADReal & G,
const ADReal & K,
const Real & G,
const Real & K,
ADRankTwoTensor & elastic_strain_incr)
{
const ADReal eqv_stress = std::sqrt(1.5) * stress_dev.L2norm();
......@@ -124,8 +124,8 @@ LynxADPlasticModel::plasticUpdate(ADRankTwoTensor & stress_dev,
ADReal
LynxADPlasticModel::plasticIncrement(const ADReal & eqv_stress,
const ADReal & pressure,
const ADReal G,
const ADReal K)
const Real G,
const Real K)
{
// Initialize hardening
if (_has_hardening)
......@@ -168,7 +168,7 @@ LynxADPlasticModel::plasticIncrement(const ADReal & eqv_stress,
}
void
LynxADPlasticModel::initPlasticParameters(const ADReal & pressure, const ADReal & K)
LynxADPlasticModel::initPlasticParameters(const ADReal & pressure, const Real & K)
{
_alpha_0 = std::sqrt(3.0) * std::sin(averageProperty(_friction_angle_0) * libMesh::pi / 180.0);
_alpha_res =
......
......@@ -18,7 +18,7 @@ LynxADThermalBase::validParams()
{
InputParameters params = LynxADMaterialBase::validParams();
params.addClassDescription("Base class for calculating the thermal properties.");
params.addCoupledVar("porosity", "The porosity auxiliary variable.");
params.addCoupledVar("porosity", 0.0, "The porosity auxiliary variable.");
params.addParam<std::vector<Real>>("radiogenic_heat_production",
"The radiogenic heat production value.");
return params;
......@@ -27,18 +27,18 @@ LynxADThermalBase::validParams()
LynxADThermalBase::LynxADThermalBase(const InputParameters & parameters)
: LynxADMaterialBase(parameters),
// _coupled_porosity(isCoupled("porosity")),
_porosity(isCoupled("porosity") ? adCoupledValue("porosity") : adZeroValue()),
_porosity(coupledValue("porosity")),
_heat_source(isParamValid("radiogenic_heat_production")
? getLynxParam<Real>("radiogenic_heat_production")
: std::vector<Real>(_n_composition, 0.0)),
_coupled_dens(hasADMaterialProperty<Real>("solid_density")),
_coupled_dens(hasMaterialProperty<RealVectorValue>("gravity_vector")),
_rho_f(_coupled_dens ? &getADMaterialProperty<Real>("fluid_density") : nullptr),
_rho_s(_coupled_dens ? &getADMaterialProperty<Real>("solid_density") : nullptr),
_thermal_diff(declareADProperty<Real>("thermal_diffusivity")),
_rhoC_b(declareADProperty<Real>("bulk_specific_heat")),
_rhoC_f(declareADProperty<Real>("fluid_specific_heat")),
_thermal_exp(declareADProperty<Real>("thermal_expansion_coefficient")),
_radiogenic_heat(declareADProperty<Real>("radiogenic_heat_production")),
_thermal_exp(declareProperty<Real>("thermal_expansion_coefficient")),
_radiogenic_heat(declareProperty<Real>("radiogenic_heat_production")),
_c_f(_fe_problem.getMaxQps()),
_c_s(_fe_problem.getMaxQps()),
_lambda_f(_fe_problem.getMaxQps()),
......@@ -65,7 +65,7 @@ LynxADThermalBase::computeQpSpecificHeat()
_rhoC_f[_qp] = _coupled_dens ? (*_rho_f)[_qp] * _c_f[_qp] : 0.0;
_rhoC_f[_qp] = 0.0;
ADReal rhoC_s = _coupled_dens ? (*_rho_s)[_qp] * _c_s[_qp] : 0.0;
_rhoC_b[_qp] = computeMixtureProperty(_rhoC_f[_qp], rhoC_s);
_rhoC_b[_qp] = computeADMixtureProperty(_rhoC_f[_qp], rhoC_s);
}
void
......@@ -94,12 +94,14 @@ LynxADThermalBase::computeQpThermalSource()
_radiogenic_heat[_qp] = averageProperty(_heat_source);
}
Real
LynxADThermalBase::computeMixtureProperty(const Real fluid_prop, const Real solid_prop)
{
return _porosity[_qp] * fluid_prop + (1.0 - _porosity[_qp]) * solid_prop;
}
ADReal
LynxADThermalBase::computeMixtureProperty(const ADReal fluid_prop, const ADReal solid_prop)
LynxADThermalBase::computeADMixtureProperty(const ADReal fluid_prop, const ADReal solid_prop)
{
// if (_coupled_porosity)
// return (*_porosity)[_qp] * fluid_prop + (1.0 - (*_porosity)[_qp]) * solid_prop;
// else
// return solid_prop;
return _porosity[_qp] * fluid_prop + (1.0 - _porosity[_qp]) * solid_prop;
}
\ No newline at end of file
......@@ -46,31 +46,6 @@
[../]
[]
[AuxVariables]
[./vol_strain]
order = CONSTANT
family = MONOMIAL
[../]
[./pressure]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./vol_strain_aux]
type = LynxADVolStrainRateAux
variable = vol_strain
strain_type = elastic
execute_on = 'TIMESTEP_END'
[../]
[./pressure_aux]
type = LynxADEffectivePressureAux
variable = pressure
execute_on = 'TIMESTEP_END'
[../]
[]
[BCs]
[./T1_everywhere]
type = PresetBC
......