Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SmallDeformation nonlocal #2294

Merged
merged 12 commits into from
Dec 19, 2018
Merged

SmallDeformation nonlocal #2294

merged 12 commits into from
Dec 19, 2018

Conversation

endJunction
Copy link
Member

@endJunction endJunction commented Dec 11, 2018

Don't worry, most changes are the test cases, the code is around 2kLOC.

Based on Francesco's paper (published, not yet available online) "Experimental characterization and numerical modelling of fracture processes in granite".

Main difference to the SmallDeformation process is the additional, so called non-local integration loop.

@endJunction
Copy link
Member Author

@wenqing @TomFischer @renchao-lu Let's meet on Friday at 11 for discussion.

Copy link
Member

@wenqing wenqing left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good except some small stuffs.

@@ -53,6 +53,19 @@ template <int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> sOdotS(
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& v);

template <int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
elasticTangentStiffness(double const K, double const G)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a similar function in MaterialLib/SolidModels/LinearElasticIsotropic.cpp, and it is called getElasticTensor. It might be good to use the same name.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. I chose the version with the first Lame parameter instead of the bulk modulus.

config.getConfigParameter<double>("internal_length");

SmallDeformationNonlocalProcessData<DisplacementDim> process_data{
materialIDs(mesh), std::move(solid_constitutive_relations),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this line formatted by using clang-format?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. Table alignment.

{
double const damage = (1 - beta_d) * (1 - std::exp(-kappa_d / alpha_d));

if (damage < 0. || damage > 1.)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it necessary that let damage = damage < 0. ? 0 : damage and damage = damage > 1. ? 1 : damage?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Put an assert there and check for half open interval [0, 1).

// Compute sigma_eff from damage total stress sigma
using KelvinVectorType = typename BMatricesType::KelvinVectorType;
KelvinVectorType const sigma_eff_prev =
sigma_prev / (1. - damage_prev);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a possibility that damage_prev = 1.0?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Damage can not be 1, see calculateDamage(). Comment added.

@@ -18,7 +18,7 @@ double Invariants<KelvinVectorSize>::equivalentStress(
Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
{
assert(std::abs(trace(deviatoric_v)) <=
5e-14 * diagonal(deviatoric_v).norm());
2e-13 * diagonal(deviatoric_v).norm());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How did you find these magic number?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

reduce it until it works ;)

{//! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__process_variables__process_variable}
"process_variable"});

DBUG("Associate displacement with process variable \'%s\'.",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is the symbol ' escaped?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just to be on the safe side ;)
Removing.

if (damage < 0. || damage >= 1.)
{
ERR("Damage value %g outside of [0,1) interval.", damage);
assert(false && "damage not in [0,1) interval.");
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not OGS_FATAL?

typename BMatricesType::KelvinVectorType eps, eps_prev;
double free_energy_density = 0;
double damage = 0; ///< isotropic damage
double damage_prev = 0; ///< \copydoc damage
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the '\copydoc' intended?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes.

NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;

// TODO move to NumLib::ExtrapolatableElement
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this TODO easy to fulfil?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No. Is related to output and the TODO is in all processes.

{
namespace SmallDeformationNonlocal
{
/// Used by for extrapolation of the integration point values. It is ordered
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Either 'by' or 'for'.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixing.

if (ehlers_solid_material == nullptr)
{
OGS_FATAL(
"The SmallDeformationNonlocalLocal supports only Ehlers "
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMHO the word 'process' is missing.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixing.

_element.getID());
}

if (name == "sigma_ip")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This string comparison is done for every integration point. Is this computational expensive? Would an enum an option?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Enum is not an option because each process might call this differently. The function call (and so the comparison) happens once for each element when setting initial conditions.

BMatricesType, ShapeMatricesType, DisplacementDim>>>
_ip_data;

IntegrationMethod _integration_method;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not const?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed.

private:
int const _integration_order;
std::function<std::vector<std::vector<double>>()> _callback;
};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please move the writer to separate files.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Moved into single separate file.

@endJunction
Copy link
Member Author

@fparisio Once the paper is published we need to add the corresponding test descriptions to the web benchmark documentation page.

@endJunction endJunction merged commit 0fba06a into ufz:master Dec 19, 2018
@endJunction endJunction deleted the SDNonlocal branch December 19, 2018 18:07
@ogsbot
Copy link
Member

ogsbot commented Jun 19, 2020

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
5 participants