use conduction timestep by default

[maxrudolph]

Per discussion on the aspect-devel mailing list, this pull request enables the conduction timestep by default. This results in improved accuracy for systems that are very close to the critical Rayleigh number.

[bangerth]

/run-tests

[gassmoeller]

Wow that are a lot of failed tests for such an innocent looking patch. Seems like some of the tests should have used the conduction timestep all the time (lots of smaller timesteps now). But I am more worried about the ones that now crash ("Test aborted"), or throw exceptions ("Product of rho * cp has to be non-negative"). Max if you have some time for this, could you take a look which tests crash exactly and maybe figure out why? After all a smaller timestep should only improve things, maybe there is something wrong in the computation of the new time step. If you do not have the time, or you would like some help, let us know.

[maxrudolph]

Yes, I was surprised at the sheer number of tests that were broken. I can find some time to go through the test output and try to figure out whether any of the tests were broken in a significant way.

[bangerth]

Debugging help with this would be much appreciated! It's a bit of a shocker that a five-character change breaks 73 tests :-)

[maxrudolph]

I've had a super busy week with no time to look into this, but it is on the 'to-do' list.

[maxrudolph]

@tjhei would you mind re-running the tests on this version? Looks like the tester timed out.

[tjhei]

several tests are timing out:

depth_average_03
depth_average_04
point_value_02
minimum_refinement_function_time_dependent
diffusion
point_value_02

I would assume they need to be set up to run without conduction timesteps.

[gassmoeller]

From what I can see there are two types of test failures for this PR.

1. A lot of the tests complain that rho*cp is negative in the compute_time_step function. In my view we should simply remove this assertion (helper_functions.cc:572). This check is likely carried over from the assembly of the advection system (where it would make no sense to have negative rho*cp), but in fact most of the failing tests only solve the incompressible stokes equation and there we do not care if the density is negative. Even if the density would be negative it would not affect the result of the compute_time_step function, because the conduction time step is only computed if the thermal diffusivity is positive, and that is checked right after the assertion. Therefore, it should be safe to remove the assertion.

2. A number of tests now have much smaller timesteps than before (these are the ones that should have used the conduction time step from the beginning). Likely the tests that timed out also fall into this category. I disagree that we should simply not use the conduction time step for them, because that kind of defeats the purpose of this PR. Instead we would need to go the tedious route of running each of these tests, figuring out the new timestep and the setting the new end time to a reasonable value (maybe so that they make approximately as many timesteps as before the PR). We also need to make sure they still test what they were supposed to test if they wanted to test something that is actually time dependent. Afterwards we can update the test results. I am willing to help with this, but I would like to wait until 1. is fixed so that we see the actual number of tests that fall into the second category.

[bangerth]

Would we also need to adjust some code to not use the conduction time step if we only solve the Stokes equations? We can, after all, switch off the temperature solve altogether in ASPECT, at which point there is no reason any more to even consider a conduction time step -- it's just an invalid value we should disregard.

[gassmoeller]

Hmm, that is a valid point, though the only nonlinear solver scheme that does not solve the temperature equation is Stokes only and for that it does not matter, because it quits after one Stokes solve anyway (i.e. the time step is computed, but not used). Should we still check for that solver scheme and skip the conduction time step computation?

This would become a problem if somebody implemented a nonlinear solver scheme that solves the stokes and composition equation, but not the temperature. Hopefully until then we have the rework of the nonlinear solver schemes done, and can somehow interrogate the solver scheme about this kind of thing (something to keep in mind when working on #972)

[bangerth]

Yes, let's not deal with this possibility at this time then. I had not remembered that "Stokes only" implies "one time step only".

I think I need to read up in the code how we actually do that... :-)

[jdannberg]

We had a discussion about this problem (@gassmoeller, @naliboff, @LudovicJnnt and I), and there is one problem we saw with using the conduction time step in all models by default.
The convection and the conduction time step scale differently with cell size: the conduction time step scales with cell_size^2/kappa, the convection time step with cell_size/velocity. That means for every given model solution, if you increase the resolution going to smaller and smaller cells, at some point you will reach a cell size where the conduction time step is smaller than the convection time step. We computed when this crossover would happen for typical model parameters (kappa = 1e-6 m²/s; velocity = 10 cm/yr). For these parameters, we got a cell size of 300 m. So for mantle convection models, this is not likely to be a problem. However, for lithosphere dynamics simulations, I can see this becoming an issue. They often need a high resolution, and often have even smaller velocities of only around 1 cm/yr (so the crossover would be at a 3 km cell size). So for these models, the time step would then always be dominated by the conduction time step, and it might not be clear to the users why their time steps are so small.
So I guess at some point we should discuss what would be the best solution here.

[bangerth]

I think that's a good point. Do you propose to just document this, or do you actually want to change the default again?

[jdannberg]

I don't think we ever changed the default (this pull request is still open). But because of this issue, it was not clear to us what would be the best solution: using or not using the conduction time step by default.

[maxrudolph]

I would think that the default behavior should be to guarantee a n accurate solution regardless of the nature of the problem. Thus, using the minimum of the advection and condution timesteps seems like a conservative default. If the user understands what they are doing, they can disable the conduction timestep.

[bangerth]

It's a good question. I tend to like erring on the safe side like @maxrudolph does. It's true that if you have models where the conduction time step dominates, you get small time steps. But you get them for good reasons -- because you won't get a very accurate solution otherwise.

At least that was true when we used an explicit time discretization. Does anyone have intuition what happens if your time step is too large with the implicit time discretization we use these days?

[jdannberg]

That was exactly the point we weren't sure about: As the conduction time step we use only comes from a stability criterion for explicit time stepping, it was not clear to us if it is a measure for accuracy and should always be considered, or if it would just make the time steps in lithospheric models really small without there being a good reason for it.

[maxrudolph]

The problem that I was working on when I opened this PR was growth of rayleigh-taylor instabilities near the conductive-convective transition. If you don't use the conductive timestep, and take advection-limited timesteps, you will suffer numerical diffusion and hence underpredict the growth rate of convective instabilities. This would suppress for instance the onset of sublithospheric small scale convection.

[bangerth]

I don't have the intuition to tell one way or the other.

@jdannberg -- you calculated that we are unlikely to run into problems for global models because there we would have to have a mesh size of 300 meters before the cross-over point, but that it may become an issue for lithosphere models. But we don't have to obey the time step limits with a factor of one, just like we are not bound to have a CFL of one or less. We use an implicit time stepping method, after all, and the worst that can happen is that we become inaccurate. So we could say that we choose the time step as the lesser of (i) the convection time step times the CFL number, and (ii) the conduction time step times another number. That latter number we could choose to be 5 or 10, for example. Would that be good enough for lithospheric models? How large would that number have to be so that the conduction time step never dominates?

[jdannberg]

That's exactly the problem, the conduction time step scales with the square of the cell size, while the convection time step scales only linearly will cell size. So there is no such value, if you assume that people can always go to smaller scales with their models and for example model crustal deformation, where the whole domain size would only be 10s of kilometres.

I agree with @maxrudolph that for some problems it is essential to use a conduction time step. But for other models the time step might just be 10 or 100 times smaller than what you would expect from the CFL criterion, and that is confusing for the user, because the CFL number is a very prominent input parameter used in many input files, while the conduction time step is not used very often.
One may argue that if we make using the conduction time step the default, it's the safe option, because the time step will never be too large for the modelled problem. But maybe we can find a way that works well for all models.

So whatever we decide, we have to make sure to carefully document it.

[bangerth]

On 05/23/2018 01:11 PM, Juliane Dannberg wrote: That's exactly the problem, the conduction time step scales with the square of the cell size, while the convection time step scales only linearly will cell size. So there is no such value, if you assume that people can always go to smaller scales with their models [...]

Well, true. But will people ever go to scales less than, say, 100 meters? I mean, surely the Stokes flow model ceases to be valid at some scale and using a code like ASPECT just doesn't make sense any more then.

[maxrudolph]

To offer some counterpoint, one can never anticipate fully what the code will be used for. There is a MS student at Portland State right now using ASPECT to model magma deformation on length scales of decimeters.

[bangerth]

On 05/23/2018 01:24 PM, Max Rudolph wrote: To offer some counterpoint, one can never anticipate fully what the code will be used for. There is a MS student at Portland State right now using ASPECT to model magma deformation on length scales of decimeters.

Bizarre. I love open source software!

[jdannberg]

I had to think about this again when I read @maxrudolph's comment about magma dynamics, and also discussed it with @gassmoeller. I think as long as the dominant physical process in the domain is convection and/or conduction, the conduction time step should be switched on.

If you make the domain smaller, but keep a similar Rayleigh number, this would imply that either the buoyancy forces increase, the viscosity decreases or the conductivity decreases. So as long as the system is still convecting, no matter how small your model domain, the convection time step will control the time step.

And even if the physics and the domain size stay the same and you just increase the resolution so that the conduction time step becomes dominant, this is probably also reasonable: On these small length scales, conduction will become the dominant process at some point, and if you decrease your cell size to that point, that implies that you're interested in the processes that are dominant on those scales.

The only problem comes in for the lithosphere dynamics models: They are not convection/diffusion dominated, they are usually controlled by the boundary conditions. So in some way, they are looking at a different type of physical behaviour. And so they may want to use a different way to choose their time step.

So I would say, let's move forward with this (maybe at the hack or whenever one of us has the time), and one of us can look at all of the failing tests (and for all the lithosphere models, we may want to switch off the conduction time step manually).
In addition, we should add a section to the manual about how the time step is chosen, and carefully document that we have a new default.

[naliboff]

Quoting from Juliane:

So I would say, let's move forward with this (maybe at the hack or whenever one of us has the time), and one of us can look at all of the failing tests (and for all the lithosphere models, we may want to switch off the conduction time step manually).
In addition, we should add a section to the manual about how the time step is chosen, and carefully document that we have a new default.

I agree with all of this. If we do have the conduction time step option enabled by default, it does just need to be carefully documented in the manual and in all of the lithospheric-scale cookbooks, benchmarks and tests we should also document in the parameter file why we don't use this option.

Quoting from Wolfgang:

Well, true. But will people ever go to scales less than, say, 100 meters? I mean, surely the Stokes flow model ceases to be valid at some scale and using a code like ASPECT just doesn't make sense any more then.

Like Max, I'm also working on projects with ASPECT that will examine deformation at quite small scales (1-10 meters). The context is ductile shear zones (diffusion + dislocation creep, grain-size evolution), where the viscous approximation is appropriated and ASPECT is ideally suited. Just wait until groups start to model short-term seismic deformation ;)

[jdannberg]

Well, true. But will people ever go to scales less than, say, 100 meters? I mean, surely the Stokes flow model ceases to be valid at some scale and using a code like ASPECT just doesn't make sense any more then.

Just to add to this: I've done models of shear bands in partially molten rocks on a laboratory length/time scale, so 1 x 4 mm 😄
It's useful to compare model outputs to experiments.

[bangerth]

Alright alright, I'm clearly not imaginative enough :-)

Either way, I agree as well that this patch should go in.