Add viscoelastic material model and elastic outputs functionality #1593

naliboff · 2017-05-12T20:58:55Z

The next step is to incorporate the right-hand side terms.

tjhei · 2017-05-12T21:19:53Z

include/aspect/elastic.h

+    {
+      public:
+        ElasticOutputs (const unsigned int n_points,
+                     const unsigned int /*n_comp*/)


remove n_comp

tjhei · 2017-05-12T21:20:02Z

include/aspect/elastic.h

+                     const unsigned int /*n_comp*/)
+        {
+          elastic_viscosities.resize(n_points);
+          elastic_evolutions.resize(n_points);


reinit elastic_evolutions

tjhei · 2017-05-12T21:26:39Z

do you want to derive from AdditionalMaterialOutputsStokesRHS instead?

naliboff · 2017-05-16T17:52:24Z

The viscoelastic code currently matches a time-dependent analytical solution for a few time steps , but eventually the solution blows up. No need for a review until I have a better handle on what exactly is going awry.

naliboff · 2017-08-21T21:53:17Z

With commit 2643312 the code is now passing a time-dependent analytical benchmark. Still some more work to do before a formal review (additional tests, particles), but I think the code structure is now largely correct if anyone would like to try it with a specific test case.

naliboff · 2017-09-08T20:57:00Z

Update and question:

The current implementation only matches analytical solutions for problems with a fixed strain-rate. I think I know what the issue is, but the fix will require accessing basis function values (grad_phi_u) within the material model.
Is there any major pitfall or impediment to accessing grad_phi_u through finite_element_values as in stokes.cc, or is there already a way to access these values through the material model interface?

bangerth · 2017-09-08T21:04:05Z

The MaterialModelInputs structure, as far as I know, only gives you access to the strain rate (i.e., the symmetric gradients). Do you really need the real gradient?

naliboff · 2017-09-08T21:20:59Z

I'm pretty sure what other people have used is the real gradient, but I will start with the symmetric gradients to get the new code structure in place and see how this affects the results.

bangerth · 2017-09-08T21:58:25Z

Hm, interesting. The reason why the strain is used in most material laws (and in the Stokes and the elasticity equations) is that it is a tensor that rotates appropriately with the coordinate system. In other words, quantities such as the scalar strain rate are independent of the coordinate system. This is not the case for the gradient of the velocity -- it depends on the coordinate system, and consequently is not usually used in material descriptions.

naliboff · 2017-11-18T00:37:32Z

@bangerth Everything you outlined above makes sense. For now I've reverted use the symmetric gradient velocity grades (grads_phi_u) as this is what was used in other deal.II elasticity examples and I don't have intuitive sense yet for why one would use the full gradient. I'll revisit this after a discussion with @cedrict.

I've added a number of changes and can now reproduce additional benchmarks (ex: viscoelastic bending beam, ..). However, still some work to do condensing commits, adding additional features and improving the benchmark cases. I actually removed a few benchmark cases that will be added in with separate pull requests.

For now, is it a reasonable strategy to keep this pull request as simple as possible and then continue adding complexity with additional pull requests?

The large section of code that remains to be added is an option to track the elastic stresses with tracers.

naliboff · 2018-02-08T00:07:15Z

This should be ready for a full review now. There have been a significant number of changes since the last pull request and I condensed everything down to 2 commits to help clean things up. I'll add in two benchmarks with a separate commit.

jdannberg

Thanks for contributing this! It'll be really cool to have elasticity in ASPECT, and it looks like you've put a lot of effort into this.
I left you a few comments, and I guess you will add a benchmark/cookbook in a later pull request?

jdannberg · 2018-02-08T17:52:24Z

include/aspect/material_model/viscoelastic.h

+     *
+     * The first 3 (2D) or 6 (3D) compositional fields are used to track
+     * individual components of the viscoelastic stress tensor, while additional
+     * fields represent distinct rock types.


I feel it would make a lot of sense to give those compositional fields specific names (like strain_xx, strain_xy, etc.), and then asking for those names (like it's done for example in the melt material models with the porosity field), because then you don't rely on the order of the fields. It's also safer, because then the user would have to define those fields explicitly in the input file, and you could make sure that they are there and otherwise throw an error (otherwise it could just happen that someone uses 3 different rock types, doesn't think about it, and then those fields will be used to compute elastic stresses).

I agree the naming convention is much safer and I actually meant to do this before updating the pull request. The only advantage to the current method is that it only requires two lines of codes, but it certainly could lead to all sorts of issues. However, the user will still have to be quite careful when assigning properties to a large number of fields in 3D (>10 in some cases?). Would it be worth requiring both a specific naming scheme and ordering as an extra 'safety' feature?

It would certainly be much easier to code that way. You could basically leave most of the code as it is, and just add an AssertThrow that makes sure that the first 3 (or 6) components are named stress_xx, stress_xy, etc. I think that would make a lot of sense.

For now I am fine either way. I can imagine another plugin requiring that its special properties have to be fields 0..2, then we have a problem. But that could be resolved when it finally crops up. Make sure you require the existence of the names though, otherwise allowing for variable ordering but requiring the names in the future will break compatibility.

update documentation

jdannberg · 2018-02-08T19:38:04Z

source/material_model/viscoelastic.cc

+                                   "listed compositional fields in the parameter file."
+                                   "\n\n "
+                                   "Expanding the model to include non-linear viscous flow (e.g., "
+                                   "diffusion/dislocation creep) and plasticity would produce a "


There's a word missing here.

jdannberg · 2018-02-08T19:51:09Z

include/aspect/material_model/viscoelastic.h

+     * are given that value. Other lengths of lists are not allowed. For a given
+     * compositional field the material parameters are treated as constant,
+     * except density, which varies linearly with temperature according to the
+     * thermal expansivity.


I think that piece of information should rather go into the description of the class at the end of the cc file, so that it shows up in the manual. I suspect that's where users would look first when they want to find out how to write their input file.

My view on this is to make the description of the class at the end of the .cc file an identical copy of this documentation. It is meant for the same audience.

jdannberg · 2018-02-08T19:52:42Z

include/aspect/material_model/viscoelastic.h

+     * When more than one field is present at a point, they are averaged
+     * arithmetically. An exception is viscosity, which may be averaged
+     * arithmetically, harmonically, geometrically, or by selecting the
+     * viscosity of the composition with the greatest volume fraction.


Same as above: I think it would be really useful to have that part show up in the manual.

jdannberg · 2018-02-08T19:55:42Z

include/aspect/material_model/viscoelastic.h

+    class Viscoelastic : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
+    {
+      public:
+


Remove empty line.

jdannberg · 2018-02-08T23:17:51Z

source/simulator/assemblers/stokes.cc

+      Assert(!this->get_parameters().enable_elasticity
+             ||
+             outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()->rhs_e.size()
+             == n_points, ExcInternalError());


You can combine all of your new lines with the ones below for the right-hand side, as you do exactly the same things for elasticity as for additional_stokes_rhs. So no need to copy all of that, just combine the if-statements and the Asserts.

jdannberg · 2018-02-08T23:19:15Z

tests/viscoelastic_fixed_elastic_time_step.prm

@@ -0,0 +1,132 @@
+#  Global parameters


Can you add a short description of what the test is supposed to do?

jdannberg · 2018-02-08T23:20:13Z

tests/viscoelastic_fixed_elastic_time_step.prm

+set End time                               = 1e3
+set Use years in output instead of seconds = true
+set Linear solver tolerance                = 1e-7
+set Nonlinear solver scheme                = IMPES


Why do you use IMPES for your test cases, isn't the problem nonlinear?

In this case I believe the problem is linear as there is no plasticity or non-linear viscous flow terms. Would one need to use the iterative solver simply due to having the stress (strain-rate dependent) in the constitutive law?

I guess I never thought about that, as models with elasticity, but without plasticity don't come up very often. But I think as soon as your viscosity depends on the solution variables of the Stokes system, the problem is nonlinear.
You can easily test this by just using the iterated IMPES solver, and if the nonlinear solver converges after the second iteration with a very small residual, the problem is linear; otherwise you probably need the nonlinear solver.

jdannberg · 2018-02-08T23:31:24Z

tests/viscoelastic_numerical_elastic_time_step.prm

@@ -0,0 +1,132 @@
+#  Global parameters


Can you add a description here too?
Something like, this test is exactly like the test XXX, with the only difference being...

I also wondered: The end time of the model is exactly the same as your elastic time step, so both tests do exactly the same thing. Or is this what you're testing? Wouldn't it make sense to have a model where the elastic time step is different from the real time step, and then see if the stress averaging option works and accounts for the differences?

I also think it would make sense to add the statistics files to the test outputs, because they include much more information (and digits), so it's much easier to see when things have changed and something is different than before.

jdannberg · 2018-02-08T23:32:02Z

tests/viscoelastic_stress_averaging.prm

@@ -0,0 +1,132 @@
+#  Global parameters


naliboff · 2018-02-09T00:06:29Z

@jdannberg - Thanks for the review! I'll likely go through most of the changes tomorrow, but a few initial replies are forthcoming.

naliboff · 2018-02-15T19:36:37Z

@jdannberg, @bangerth - Thanks for the additional discussion on various points above. @cedrict and I are having a follow-up conversation and comparing results, which I will summarize here when we are satisfied. I'll start pushing various updates in the next few days for the benchmarks, tests and stylistics/re-organizational requests.

naliboff · 2018-03-27T03:12:39Z

@jdannberg - I think I've addressed most of the comments and I've also added two benchmarks. I need to revisit the FEFieldFunction method for obtaining previous solutions, but I'm not sure now whether that is necessary. Something to discuss.

The last major item on the to-do list is the particle functionality, but I think that should go in a separate pull request.

jdannberg · 2018-03-27T03:17:40Z

Okay, great, I'll have a look later this week!
(Before, I still have to look at @gassmoeller's open pull requests I promised to review weeks ago...)

gassmoeller

Only half-through, but have to leave. Rest of the comments coming later. Overall very clean code, just a few places where we need to make a decision about which way to go is best.

gassmoeller · 2018-03-29T00:02:04Z

source/material_model/viscoelastic.cc

+      if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && in.strain_rate.size() > 0)
+        {
+          // Get old (previous time step) velocity gradients
+          const QGauss<dim> quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.velocities).degree+1);


This is dangerous because it assumes the quadrature never changes. But it does, e.g. for the postprocessors you will have positions at the vertices instead of the quadrature points (as for the assembly). The solution would be to use the mapping and cell information to convert the global positions into reference cell coordinates and create a quadrature at those positions (we do not care about the quadrature weights here). You could replace this and the following line by this (ignore indentation):

std::vector<Point<dim> > quadrature_positions(in.position.size()); for (unsigned int i=0; i < in.position.size(); ++i) quadrature_positions[i] = this->get_mapping().transform_real_to_unit_cell(in.current_cell, in.position[i]); FEValues<dim> fe_values (this->get_mapping(), this->get_fe(), Quadrature<dim>(quadrature_positions), update_gradients);

This will be somewhat slower, but hopefully not dramatically, but it will be much safer. And actually I suspect that your current output it slightly wrong, and will be more accurate with this change.

Great, thanks for the straightforward fix! In terms of the current output being wrong, this should only apply to the post-processor data, yes?

The only 'output' derived from this section of code is this the compositional field reaction terms and the stokes RHS terms. My interpretation was that these values would only be used/updated on the quadrature points. If this is not the case, then there certainly could be some errors in the solution.

gassmoeller · 2018-03-29T00:06:06Z

benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm

+# compositional fields. Significantly, the time step is limited to 1 Kyr through
+# the "Maximum time step" parameter and a relatively high (0.5) CFL value. Using
+# a constant time step ensures the effective (viscoelastic) viscosity of the
+# viscoelastic beam and viscous fluid remains constants throughout the model run.


Does that mean the solution is time step dependent otherwise? Is there a separate elastic timestep from the convection timestep?

There is an option in viscoelastic.cc to use a separate elastic time step from the convection time step (I call it 'numerical' time step). However, in the documentation it is strongly recommended that one uses the averaging scheme option to account for the difference between the elastic and numerical time step.

The ideal case is to set the elastic time step equal to the numerical time step, which I have done in this benchmark. For simplicity and testing purposes, I simply restricted the numerical time step so the effective (viscoelastic) viscosities are constant through time.

One only really runs into trouble if the convection time step becomes much smaller than realistic elastic time steps (~ 10^3-10^5 years). Unfortunately, this can happen when using a really fine mesh.

This situation is effectively the opposite of what the operating splitting was designed for (reaction time scale <<<<< numerical time step). Maybe others will have some insight on how to bridge these time scales.

Most time steps in lithospheric deformation simulations are on the order of 10^3-10^4 years, so this issue does not occur too often.

gassmoeller · 2018-03-29T00:07:48Z

benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm

+set Output directory                       = output
+set Timing output frequency                = 1
+set Pressure normalization                 = surface
+set Surface pressure                       = 0.


Can you remove those parameters from the list that are set to default values and not very interesting for this problem? I am thinking of use direct solver, linear solver tolerance, nonlinear solver scheme, nonlinear solver tolerance, timing output frequency.

gassmoeller · 2018-03-29T00:09:41Z

benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm

+  set Free surface boundary indicators        = 
+  set Tangential velocity boundary indicators = bottom, top, right
+  set Zero velocity boundary indicators       = left
+  set Enable elasticity                       = true


I would like to remove this subsection in #2105. Would you be ok with moving the parameter into the Formulation subsection that handles all modifications to the equations?

Yes, Formulation seems like the logical place to put this parameter.

gassmoeller · 2018-03-29T00:10:34Z

benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm

+  set Temperature polynomial degree           = 1
+  set Use locally conservative discretization = false
+  set Use discontinuous temperature discretization = false
+  set Use discontinuous composition discretization = true


Oh, so you are using discontinuous elements for the stresses as well. Did that cause any problems? Just curious.

I recall it did not have much an effect on the stresses, but I would need to check again. As expected, the DG elements + limiter improved tracking of the beam interface.

Conceptually, is there any reason not to use the DG elements with the stresses as long as the limiter values are far greater than the stress magnitudes?

gassmoeller · 2018-03-29T00:23:22Z

include/aspect/material_model/viscoelastic.h

+     *
+     * The first 3 (2D) or 6 (3D) compositional fields are used to track
+     * individual components of the viscoelastic stress tensor, while additional
+     * fields represent distinct rock types.


For now I am fine either way. I can imagine another plugin requiring that its special properties have to be fields 0..2, then we have a problem. But that could be resolved when it finally crops up. Make sure you require the existence of the names though, otherwise allowing for variable ordering but requiring the names in the future will break compatibility.

gassmoeller · 2018-03-29T00:24:38Z

include/aspect/material_model/viscoelastic.h

+     * are given that value. Other lengths of lists are not allowed. For a given
+     * compositional field the material parameters are treated as constant,
+     * except density, which varies linearly with temperature according to the
+     * thermal expansivity.


My view on this is to make the description of the class at the end of the .cc file an identical copy of this documentation. It is meant for the same audience.

gassmoeller · 2018-03-29T00:49:29Z

source/material_model/viscoelastic.cc

+#include <numeric>
+#include <deal.II/base/signaling_nan.h>
+
+//using namespace dealii;


Remove this line

gassmoeller · 2018-03-29T00:51:19Z

source/material_model/viscoelastic.cc

+                                 "compositional field is called stress_xz."));
+          AssertThrow(this->introspection().compositional_index_for_name("stress_yz") == 5,
+                      ExcMessage("Material model Viscoelastic only works if the sixth "
+                                 "compositional field is called stress_yz."));


Please move all of these asserts into parse_parameters or initialize. We do not need to check this every time the model is called.

gassmoeller · 2018-03-29T00:52:42Z

source/material_model/viscoelastic.cc

+                           ?
+                           this->get_timestep()
+                           :
+                           fixed_elastic_time_step * year_in_seconds );


If you use the fixed elastic time step is the solution still dependent on the model time step as we discussed a few minutes ago?

naliboff · 2018-03-29T18:30:22Z

@gassmoeller - Thanks for the review!

jdannberg · 2018-04-05T23:25:08Z

include/aspect/material_model/interface.h

+         * momentum equation (first part of the Stokes equation) in each
+         * quadrature point.
+         */
+        std::vector<Tensor<2,dim> > rhs_e;


Can you make a separate class for this (something like ElasticOutputs) that derives from AdditionalMaterialOutputs?

jdannberg · 2018-04-05T23:25:57Z

include/aspect/material_model/viscoelastic.h

@@ -0,0 +1,273 @@
+/*
+  Copyright (C) 2014 - 2017 by the authors of the ASPECT code.


update copyright

jdannberg · 2018-04-05T23:30:06Z

include/aspect/material_model/viscoelastic.h

+     *
+     * The first 3 (2D) or 6 (3D) compositional fields are used to track
+     * individual components of the viscoelastic stress tensor, while additional
+     * fields represent distinct rock types.


update documentation

jdannberg · 2018-04-05T23:31:00Z

include/aspect/material_model/viscoelastic.h

+     * number of compositional fields. The additional field corresponds to
+     * the value for background material. They should be ordered
+     * ``background, composition1, composition2...''. However, the first 3 (2D)
+     * or 6 (3D) composition fields correspond to components of the elastic


jdannberg · 2018-04-05T23:32:03Z

include/aspect/material_model/viscoelastic.h

+        /**
+         * The first 3 (2D) or 6 (3D) compositional fields are assumed
+         * to be components of the viscoelastic stress tensor and
+         * assigned volume fractions of zero.


Or just write "All compositional fields that are stresses are assigned ..."

jdannberg · 2018-04-06T00:29:25Z

source/material_model/viscoelastic.cc

+              const double average_elastic_shear_modulus =
+                calculate_average_elastic_shear_modulus(composition,
+                                                        elastic_shear_moduli,
+                                                        viscosity_averaging);


You can use the one from elastic_out->elastic_shear_moduli[i]

@jdannberg - For some reason this broke the tests. I am going to leave the current structure for now, but happy to revisit this later.

jdannberg · 2018-04-06T00:33:25Z

source/material_model/viscoelastic.cc

+              const double dt = this->get_timestep();
+              if (use_fixed_elastic_time_step == true && use_stress_averaging == true)
+                {
+                  stress_new = ( ( 1. - ( dt / dte ) ) * stress_old ) + ( ( dt / dte ) * stress_new ) ;


Can you add an Assert in parse_parameters that dte is not zero?

jdannberg · 2018-04-06T00:40:19Z

source/material_model/viscoelastic.cc

+              // Fill reaction terms
+              for (unsigned int j = 0; j < SymmetricTensor<2,dim>::n_independent_components ; ++j)
+                out.reaction_terms[i][j] = -stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)]
+                                           + stress_new[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)];


This will not work with the iterated IMPES scheme. Can you add an Assert in Parse_parameters that makes sure we do not use the material model with iterated IMPES.

jdannberg · 2018-04-06T00:42:57Z

source/material_model/viscoelastic.cc

+                             "constant value, which can be enforced through the CFL condition and"
+                             "maximum time step parameter. When using a fixed elastic time step that"
+                             "differs the numerical time step it is strongly recommended that the stress"
+                             "averaging scheme is also applied");


Dot at the end of sentence

jdannberg · 2018-04-06T00:51:25Z

source/material_model/viscoelastic.cc

+                             "maximum time step parameter. When using a fixed elastic time step that"
+                             "differs the numerical time step it is strongly recommended that the stress"
+                             "averaging scheme is also applied");
+          prm.declare_entry ("Fixed elastic time step", "1.e3",


Lets think about what exactly this does and then document this parameter a bit more.

gassmoeller

A few more comments

gassmoeller · 2018-04-06T00:43:23Z

source/material_model/viscoelastic.cc

+      }
+      prm.leave_subsection();
+    }
+


usually 3 empty lines between functions

gassmoeller · 2018-04-06T00:44:14Z

source/material_model/viscoelastic.cc

+    Viscoelastic<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      //not pretty, but we need to get the number of compositional fields before
+      //simulatoraccess has been initialized here...


This comment is outdated, and you can remove the subsection below. this->n_compositional_fields() should work just fine anywhere you need it.

gassmoeller · 2018-04-06T00:45:11Z

source/material_model/viscoelastic.cc

+            (new MaterialModel::ElasticAdditionalOutputs<dim> (n_points)));
+        }
+    }
+


remove empty line

gassmoeller · 2018-04-06T00:49:49Z

source/simulator/assemblers/stokes.cc

                data.local_rhs(i) += (force->rhs_u[q] * scratch.phi_u[i]
                                      + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i])
                                     * JxW;

+              if (force != NULL && this->get_parameters().enable_elasticity)


If you move rhs_e out of AdditionalMaterialOutputsStokesRHS and into ElasticOutputs then you can create a pointer elastic_outputs and replace this condition by if (elastic_outputs != NULL)

jdannberg · 2018-04-06T00:55:42Z

source/simulator/assemblers/stokes.cc

@@ -229,11 +229,15 @@ namespace aspect
              data.local_rhs(i) += (density * gravity * scratch.phi_u[i])
                                   * JxW;

-              if (force != NULL)
+              if (force != NULL && this->get_parameters().enable_additional_stokes_rhs)


If you do what @gassmoeller describes below, you don't need this change anymore.

replace force with elastic_outputs and remove the second part of the if

naliboff · 2018-05-02T21:30:13Z

@jdannberg @gassmoeller - Ok, I think all of the comments have been addresses except for making a new ElasticOutputs class. The branch has also been updated to rebased to the current master. Would it be possible to do a review before implementing the ElasticOutputs restructuring? I'm hoping to simplify the review process a bit.

…ial model

jdannberg

Alright, this is more or less good to go. There are two points I would like to discuss before merging, and all of my other comments are only cleanup and can be addressed in a follow-up pull request.
(1) My comments in stokes.cc, line 218. @gassmoeller, would you agree that this should be changed?
(2) At the moment, Use fixed elastic time step = false by default. Is this what we want? This is the option that makes the viscosity time-step dependent. I just wanted to make sure that this is really the right default before we merge.

In addition, you should add an issue that lists all of the things you want to improve in the future. The things I remember are

adding a cookbook
use tracers to track the stress
determine if we can use operator splitting to enable this material model for all solver schemes
add a comment that the elasticity can be switched off by setting the shear modulus to a very large value
cleanup (remaining comments I had as part of this review)
You might have other things on a list somewhere, but it would make sense to record them on github.

jdannberg · 2018-05-24T02:51:58Z

benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm

+  end
+end
+
+# Model geometry (7.5x5 km, 0.025 km spacing)


Can you update the comment, this should be 0.1 km, correct? (7.5 km, and 75 repetitions?)

jdannberg · 2018-05-24T03:37:10Z

source/simulator/assemblers/stokes.cc

                data.local_rhs(i) += (force->rhs_u[q] * scratch.phi_u[i]
                                      + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i])
                                     * JxW;

+              if (elastic_outputs != NULL )


This implementation means that in the assembly, if the material model fills the elastic outputs, they will be assembled (no matter what the enable_elasticity parameter is set to). I just want to make sure that this is really what we want (as for the force terms, they are only assembled if this option is selected in the input file, see above). Shouldn't this rather be
if (elastic_outputs != NULL && this->get_parameters().enable_elasticity)

Yes, from our last conversation I recall we did decide on using if (elastic_outputs != NULL && this->get_parameters().enable_elasticity) . While enable_elasticity could be removed entirely, having this parameter is a safeguard to make sure users realize they are adding additional terms to the RHS side of the momentum equation.

However, this could lead to some confusion if people don't read the documentation for enable_elasticity. If enable_elasticity is set to false, perhaps the material model should provide the background viscosities to the material model outputs instead of the viscoelastic viscosity? This would allow one to run simulations with the viscoelastic material model without any elastic effects.

At the moment, the model will just crash if you set enable_elasticity = false, and the error message says that you have to switch it on if you want to use the elasticity material model. So even if people don't read the documentation, they will notice.
This would only be to make sure if people add other material models with elasticity that the elastic terms are only added if enable_elasticity = true.

And I think we agreed that if you want to use the material model without elasticity, you just have to set the elastic moduli to a large value. So no need to do any big changes, just this one line. I think we want to merge this as soon as possible.

That's right, there is already an Assert for Enable elasticity. Probably a sign that is indeed to time to merge as soon as possible :) I've added a list of future features in a new github issue and will address the remaining points in a little bit.

jdannberg · 2018-05-24T03:42:33Z

tests/viscoelastic_fixed_elastic_time_step.prm

@@ -0,0 +1,159 @@
+# This test checks viscoelastic stress buildup with a


If this is mostly a copy of an existing benchmark, it would make sense to just include the benchmark in this prm file instead of copying everything (see for example the mid_ocean_ridge test).

This test does derive from the viscoelastic stress build-up benchmark, but there are a number of differences and the other viscoelastic tests are variants from this test.

My preference is to have one test parameter file with all of the parameters specified. However, I could derive the other two viscoelastic tests from viscoelastic_fixed_elastic_time_step.prm to clean things up.

I think it's fine to have separate parameter files that specify all of the parameters. My only point was, if there is one parameter file that basically represents the benchmark case (maybe just with a lower resolution and shorter end time), it makes sense to include the benchmark prm file rather than copying it, as we will immediately notice when a pull request breaks that test (while the benchmarks are not run for every pull request), so the benchmark prm file will always be working and up to date.

I actually used the benchmark case as a template for the test for exactly that reason, but having it derive directly from the benchmark would make that perfectly clear for anyone else who runs into issues. Definitely worth the change 👍

Yes, and the other big point is that in case we ever do changes to the structure of the parameter file again, we will notice if we have to update the benchmark parameter file (in the past, it happened that we just updated the test prm file, and forgot about the file in the benchmarks folder, as it wasn't used directly in any test).

naliboff · 2018-05-24T16:14:09Z

@jdannberg - My preference is to have Use fixed elastic time step = false by default as this is the option most commonly used in the community.

gassmoeller · 2018-05-24T16:51:04Z

Hmm, I feel like the Use fixed elastic time step parameter is a controversial point. I just discussed this with @jdannberg and I can understand her point that this makes models time-step dependent and nearly impossible to reproduce. On the other hand if it is the commonly used option in the community it would be nice to allow for it. So my suggestion is: Let us force the user to make this decision. Instead of a Patterns::Bool() use a Patterns::Selection("true|false|unspecified") with unspecified as the default setting. Then parse this in the parse_parameters function and add AssertThrow(parameter != "unspecified",ExcMessage("... (this parameter has to be specified, because ... )). Is that a reasonable workaround?

naliboff · 2018-05-24T17:13:39Z

@gassmoeller @jdannberg - Yes, forcing the user to at least make a decision sounds like a good solution. Hopefully at some point soon we can find a new formulation that handles the time-dependence in a more reasonable matter or simply removes it.

jdannberg · 2018-05-24T21:18:01Z

Okay, looks good to me. You can merge yourself once the tester is happy :)
And good job, it's great that we can finally merge this!

naliboff · 2018-05-24T21:21:58Z

Yahoo! Thanks to you and Rene for all the reviews, feedback and patience!

tjhei reviewed May 12, 2017

View reviewed changes

naliboff force-pushed the viscoelastic branch 2 times, most recently from 358b628 to afd2255 Compare May 16, 2017 17:48

naliboff force-pushed the viscoelastic branch from afd2255 to c41bf2c Compare August 21, 2017 21:47

naliboff force-pushed the viscoelastic branch from 6b1740c to 3bd3098 Compare November 18, 2017 00:12

naliboff force-pushed the viscoelastic branch 2 times, most recently from a088779 to 58f4350 Compare February 8, 2018 00:04

jdannberg requested changes Feb 8, 2018

View reviewed changes

naliboff force-pushed the viscoelastic branch 2 times, most recently from 0881424 to bc57436 Compare March 20, 2018 01:17

gassmoeller reviewed Mar 29, 2018

View reviewed changes

jdannberg reviewed Apr 6, 2018

View reviewed changes

gassmoeller reviewed Apr 6, 2018

View reviewed changes

jdannberg reviewed Apr 6, 2018

View reviewed changes

naliboff force-pushed the viscoelastic branch from 8562575 to 989fed2 Compare May 2, 2018 21:22

naliboff added 15 commits May 22, 2018 17:21

Update viscoelastic documentation and formatting requirements

393a8c2

Restructure and cleanup viscoelastic implementation

f5cc347

Add viscoelastic additional outputs

ca59260

Check for access to strain rate

39bb765

Fix assert error messages

37d6e2c

Add viscoelastic stress build-up benchmark

a56428f

Add viscoelastic bending beam benchmark

3523e6a

Move enable_elasticity to Formulation subsection

5401f9b

Fix method for obtaining old velocity gradients in viscoelastic mater…

5f00c16

…ial model

Address comments

ebae58b

Add ElasticOutputs for Stokes RHS

b5371cc

Address comments

dc3669e

Minor cleanup of viscoelastic material model

4cc9f61

Update viscoelastic tests

1a804fd

Update viscoelastic bending beam benchmark

06aa82b

naliboff force-pushed the viscoelastic branch from cea8e83 to 06aa82b Compare May 23, 2018 00:23

naliboff changed the title ~~[WIP] Add viscoelastic material model and elastic outputs functionality~~ Add viscoelastic material model and elastic outputs functionality May 23, 2018

jdannberg reviewed May 24, 2018

View reviewed changes

naliboff force-pushed the viscoelastic branch 2 times, most recently from d14744c to f583bf0 Compare May 25, 2018 17:04

Minor updates to viscoelastic material model

1c4781e

naliboff force-pushed the viscoelastic branch from f583bf0 to 1c4781e Compare May 25, 2018 17:25

naliboff merged commit c757cc1 into geodynamics:master May 25, 2018

naliboff deleted the viscoelastic branch June 8, 2018 19:26

jdannberg mentioned this pull request Jun 13, 2019

move elasticity to rheology model #3107

Merged

		@@ -0,0 +1,273 @@
		/*
		Copyright (C) 2014 - 2017 by the authors of the ASPECT code.

		@@ -0,0 +1,159 @@
		# This test checks viscoelastic stress buildup with a

Add viscoelastic material model and elastic outputs functionality #1593

Add viscoelastic material model and elastic outputs functionality #1593

Conversation

naliboff commented May 12, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tjhei commented May 12, 2017

naliboff commented May 16, 2017

naliboff commented Aug 21, 2017

naliboff commented Sep 8, 2017

bangerth commented Sep 8, 2017

naliboff commented Sep 8, 2017

bangerth commented Sep 8, 2017

naliboff commented Nov 18, 2017

naliboff commented Feb 8, 2018

jdannberg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Feb 9, 2018

naliboff commented Feb 15, 2018

naliboff commented Mar 27, 2018

jdannberg commented Mar 27, 2018

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Mar 29, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented May 2, 2018

jdannberg left a comment