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.

Sign up for GitHub

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

Jump to bottom

Add option to use full finite strain tensor for strain weakening #1644

Merged
merged 5 commits into from May 17, 2017
Merged

Add option to use full finite strain tensor for strain weakening #1644

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f135f60
Add option to use full finite strain tensor for strain weakening
naliboff May 13, 2017
e29c9f5
Address comments, fix astyle
naliboff May 14, 2017
2df6e70
Add fe_values access
naliboff May 15, 2017
12d71cd
Add test output
naliboff May 16, 2017
23cbf2a
Add Dannberg as co-author for finite strain cookbook and plugin
naliboff May 16, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions include/aspect/material_model/visco_plastic.h
View file
Open in desktop
Expand Up @@ -186,6 +186,7 @@ namespace aspect
const YieldScheme &yield_type) const;

bool use_strain_weakening;
bool use_finite_strain_tensor;
std::vector<double> start_strain_weakening_intervals;
std::vector<double> end_strain_weakening_intervals;
std::vector<double> viscous_strain_weakening_factors;
Expand Down
128 changes: 113 additions & 15 deletions source/material_model/visco_plastic.cc
View file
Open in desktop
Expand Up @@ -20,6 +20,7 @@

#include <aspect/material_model/visco_plastic.h>
#include <aspect/utilities.h>
#include <deal.II/fe/fe_values.h>
using namespace dealii;

namespace aspect
Expand All @@ -41,7 +42,17 @@ namespace aspect

// assign compositional fields associated with strain a value of 0
if (use_strain_weakening == true)
x_comp[0] = 0.0;
{
if (use_finite_strain_tensor == false)
{
x_comp[0] = 0.0;
}
else
{
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
x_comp[i] = 0.;
}
}

//sum the compositional fields for normalization purposes
double sum_composition = 0.0;
Expand Down Expand Up @@ -200,11 +211,25 @@ namespace aspect
double coh = cohesions[j];

// Strain weakening
double strain_ii = 0.;
if (use_strain_weakening == true)
{

// Constrain the second strain invariant of the previous timestep
const double strain_ii = std::max(std::min(composition[0],end_strain_weakening_intervals[j]),start_strain_weakening_intervals[j]);
// Calculate and/or constrain the strain invariant of the previous timestep
if ( use_finite_strain_tensor == true )
{
// Calculate second invariant of left stretching tensor "L"
Tensor<2,dim> strain;
for (unsigned int q = 0; q < Tensor<2,dim>::n_independent_components ; ++q)
strain[Tensor<2,dim>::unrolled_to_component_indices(q)] = composition[q];
const SymmetricTensor<2,dim> L = symmetrize( strain * transpose(strain) );
strain_ii = std::max(std::min(std::fabs(second_invariant(L)),end_strain_weakening_intervals[j]),start_strain_weakening_intervals[j]);
}
else
{
// Here the compositional field already contains the finite strain invariant magnitude
strain_ii = std::max(std::min(composition[0],end_strain_weakening_intervals[j]),start_strain_weakening_intervals[j]);
}

// Linear strain weakening of cohesion and internal friction angle between specified strain values
const double strain_fraction = ( strain_ii - start_strain_weakening_intervals[j] ) /
Expand Down Expand Up @@ -274,6 +299,10 @@ namespace aspect
evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
MaterialModel::MaterialModelOutputs<dim> &out) const
{



// Loop through points
for (unsigned int i=0; i < in.temperature.size(); ++i)
{
const double temperature = in.temperature[i];
Expand Down Expand Up @@ -347,13 +376,55 @@ namespace aspect
out.reaction_terms[i][c] = 0.0;
// If strain weakening is used, overwrite the first reaction term,
// which represents the second invariant of the strain tensor
if (use_strain_weakening && this->get_timestep_number() > 0)
double edot_ii = 0.;
double e_ii = 0.;
if (use_strain_weakening == true && use_finite_strain_tensor == false && this->get_timestep_number() > 0)
Copy link
Contributor

Choose a reason for hiding this comment

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

Why should the timestep be bigger than 0?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This is to test if we are not on the first time step.

{
const double edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate);

edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate);
e_ii = edot_ii*this->get_timestep();
// Update reaction term
out.reaction_terms[i][0] = edot_ii*this->get_timestep();
out.reaction_terms[i][0] = e_ii;
}
}

// We need the velocity gradient for the finite strain (they are not included in material model inputs),
// so we get them from the finite element.
if (in.cell && use_strain_weakening == true && use_finite_strain_tensor == true && this->get_timestep_number() > 0)
{
const QGauss<dim> quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.velocities).degree+1);
FEValues<dim> fe_values (this->get_mapping(),
this->get_fe(),
quadrature_formula,
update_gradients);

std::vector<Tensor<2,dim> > velocity_gradients (quadrature_formula.size(), Tensor<2,dim>());

fe_values.reinit (*in.cell);
fe_values[this->introspection().extractors.velocities].get_function_gradients (this->get_solution(),
velocity_gradients);
Copy link
Contributor Author

Choose a reason for hiding this comment

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

The gcc-hd compiler build error seems to be related to retrieving the fe_values here. Any suggestions for a fix I can apply?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sorry, should have quoted lines 394-397, but the build error there propagates down to these lines as well.

Copy link
Contributor

Choose a reason for hiding this comment

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

I think you need #include <deal.II/fe/fe_values.h> at the top of the file.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Doh, of course. Thank you.


// Assign the strain components to the compositional fields reaction terms.
// If there are too many fields, we simply fill only the first fields with the
// existing strain rate tensor components.
for (unsigned int q=0; q < in.position.size(); ++q)
{
// Convert the compositional fields into the tensor quantity they represent.
Tensor<2,dim> strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
{
strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];
}

// Compute the strain accumulated in this timestep.
const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);

// Output the strain increment component-wise to its respective compositional field's reaction terms.
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
{
out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
}
}

}
}

Expand Down Expand Up @@ -429,6 +500,10 @@ namespace aspect
Patterns::Bool (),
"Apply strain weakening to viscosity, cohesion and internal angle "
"of friction based on accumulated finite strain. Units: None");
prm.declare_entry ("Use finite strain tensor", "false",
Patterns::Bool (),
"Track and use the full finite strain tensor for strain weakening. "
"Units: None");
prm.declare_entry ("Start strain weakening intervals", "0.",
Patterns::List(Patterns::Double(0)),
"List of strain weakening interval initial strains "
Expand Down Expand Up @@ -568,6 +643,9 @@ namespace aspect
//increment by one for background:
const unsigned int n_fields = this->n_compositional_fields() + 1;

// number of required compositional fields for full finite strain tensor
const unsigned int s = Tensor<2,dim>::n_independent_components;

prm.enter_subsection("Material model");
{
prm.enter_subsection ("Visco Plastic");
Expand Down Expand Up @@ -603,7 +681,10 @@ namespace aspect
if (use_strain_weakening)
AssertThrow(this->n_compositional_fields() >= 1,
ExcMessage("There must be at least one compositional field. "));

use_finite_strain_tensor = prm.get_bool ("Use finite strain tensor");
if (use_finite_strain_tensor)
AssertThrow(this->n_compositional_fields() >= s,
ExcMessage("There must be enough compositional fields to track all components of the finite strain tensor (4 in 2D, 9 in 3D). "));
start_strain_weakening_intervals = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Start strain weakening intervals"))),
n_fields,
"Start strain weakening intervals");
Expand Down Expand Up @@ -758,12 +839,29 @@ namespace aspect
"See, for example, Thieulot, C. (2011), PEPI 188, pp. 47-68. "
"\n\n"
"The user has the option to linearly reduce the cohesion and "
"internal friction angle as a function of the finite strain invariant. "
"The finite strain invariant may be calculated through particles or the compositional "
"field finite strain invariant plugin (see \\cite{buiter2008dissipation} benchmark). "
"In either case, the user specifies a compositional field for the finite "
"strain invariant, which must be the first listed compositional field "
"in the parameter file."
"internal friction angle as a function of the finite strain magnitude. "
"The finite strain invariant or full strain tensor is calculated through "
"compositional fields within the material model. This implementation is "
"identical to the compositional field finite strain plugin and cookbook "
"described in the manual (author: Gassmoeller). If the user selects to track "
Copy link
Member

Choose a reason for hiding this comment

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

Please add Juliane, we contributed the section together.

"the finite strain invariant ($e_ii$), a single compositional field tracks "
"the value derived from $e_ii^t = (e_ii)^(t-1) + edot_ii*dt, where $t$ and $t-1$ "
"are the current and prior time steps, $edot_ii$ is the second invariant of the "
"strain rate tensor and $dt$ is the time step size. In the case of the "
"full strain tensor $F$, the finite strain magnitude is derived from the "
"second invariant of the symmetric stretching tensor $L$, where "
"$L = F * [F]^T$. The user must specify a single compositional "
"field for the finite strain invariant or multiple fields (4 in 2D, 9 in 3D) "
"for the finite strain tensor. These field(s) must be the first lised "
"compositional fields in the parameter file. Note that one or more of the finite strain "
"tensor components must be assigned a non-zero value intially. This value can be "
"be quite small (ex: 1.e-8), but still non-zero. While the option to track and use "
"the full finite strain tensor exists, tracking the associated compositional "
"is computationally expensive in 3D. Similarly, the finite strain magnitudes "
"may in fact decrease if the orientation of the deformation field switches "
"through time. Consequently, the ideal solution is track the finite strain "
"invariant (single compositional) field within the material and track "
"the full finite strain tensor through tracers."
""
"\n\n"
"Viscous stress may also be limited by a non-linear stress limiter "
Expand All @@ -775,7 +873,7 @@ namespace aspect
"reference strain rate, $\\varepsilon_{ii}$ is the strain rate "
"and $n_y$ is the stress limiter exponent. The yield stress, "
"$\\tau_y$, is defined through the Drucker Prager yield criterion "
"formulation. This method of limiting viscous stress has been used "
"formulation. This method of limiting viscous stress has been used "
"in various forms within the geodynamic literature, including "
"Christensen (1992), JGR, 97(B2), pp. 2015-2036; "
"Cizkova and Bina (2013), EPSL, 379, pp. 95-103; "
Expand Down
116 changes: 116 additions & 0 deletions tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm
View file
Open in desktop
@@ -0,0 +1,116 @@
# Global parameters
Copy link
Member

Choose a reason for hiding this comment

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

I think you forgot the output for this test.

set Dimension = 2
set Start time = 0
set End time = 1e2
set Use years in output instead of seconds = true
set Nonlinear solver scheme = iterated Stokes
set Max nonlinear iterations = 1
set Number of cheap Stokes solver steps = 0
set Output directory = visco_plastic_yield_strain_weakening_full_strain_tensor
set Timing output frequency = 1

# Model geometry (100x100 km, 10 km spacing)
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 10
set Y repetitions = 10
set X extent = 100e3
set Y extent = 100e3
end
end

# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end


# Boundary classifications (fixed T boundaries, prescribed velocity)
subsection Model settings
set Fixed temperature boundary indicators = bottom, top, left, right
set Prescribed velocity boundary indicators = bottom y: function, top y: function, left x: function, right x: function
end

# Velocity on boundaries characterized by functions
subsection Boundary velocity model
subsection Function
set Variable names = x,y
set Function constants = m=0.0005, year=1
set Function expression = if (x<50e3 , -1*m/year, 1*m/year); if (y<50e3 , 1*m/year, -1*m/year);
end
end

# Temperature boundary and initial conditions
subsection Boundary temperature model
set Model name = box
subsection Box
set Bottom temperature = 273
set Left temperature = 273
set Right temperature = 273
set Top temperature = 273
end
end
subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 273
end
end

# Compositional fields used to track finite strain
subsection Compositional fields
set Number of fields = 4
set Names of fields = s11, s12, s21, s22
end

# The compositional fields have an initial value 0 or 0.5
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0.5; 0; 0; 0.5;
end
end

# Boundary composition specification
subsection Boundary composition model
set Model name = initial composition
end

# Material model (values for background material and strain fields)
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Reference strain rate = 1.e-16
set Viscous flow law = dislocation
set Prefactors for dislocation creep = 5.e-23
set Stress exponents for dislocation creep = 1.0
set Activation energies for dislocation creep = 0.
set Activation volumes for dislocation creep = 0.
set Yield mechanism = drucker
set Angles of internal friction = 0.
set Cohesions = 1.e6
set Use strain weakening = true
set Use finite strain tensor = true
set Start strain weakening intervals = 0.
set End strain weakening intervals = 1.0
set Cohesion strain weakening factors = 0.5
set Friction strain weakening factors = 0.5
end
end

# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.0
end
end

# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, mass flux statistics
end