Hi all,
I have a problem with numerical instabilities causing spurious strainrates. The model setup is as follows:
- isoviscous, closed, free-slip box, 660x2640 km
- cold lithosphere at surface
- 250-K thermal instability at the bottom boundary at x=1000 km (to develop mantle plume)
I’ve stripped the model setup down to a small model with very few ingredients, but the problem still remains. The following images illustrate the problem. The strainrate field at t=0 looks like:
with velocity field:
and a corresponding temperature field that looks smooth:
The corresponding input file can be found here:
The mesh is 128x32 elements, without mesh refinement, and solver tolerance is set low (1e-9). The visco-plastic material model is used (it is isoviscous now, nbut the original problem was rheologically more complex). Local mesh refinement seems to make it worse.
The problem occurs (or can be noticed) only when overall flow velocities are low. Velocities at t=0 are small, but not beyond the point of resolving, I would have thought (up to 0.4 mm/yr, although this is probably at least in part due to the numerical instability). The problem remains for for later model times, until velocities are large enough to overshadow this problem.
Has anyone encountered this problem, and if so, is there a solution?
Many thanks,
Jeroen van Hunen
Hi Jereon,
Thanks for posting this to the forum! The short answer is I have seen this type of instability (banding) in brittle layers before deformation localizes and in some cases throughout the lithosphere (brittle and viscous) for the same scenario (no localization).
However, I don’t recall ever seeing this in purely viscous models and the instabilities typically disappear after 1 or 2 time steps.
I’ll try experimenting with the parameter file later this afternoon to see if I can pinpoint any commonality with previous models exhibiting this behavior.
Cheers,
John
Hi Jeroen,
my guess would be that at the start of the model, there is just no buoyancy to drive the flow, so it’s really hard for the solver to figure out an accurate solution. I am not quite sure what the numerical reason is, but @rgassmoeller suggested it might have to do with the solver tolerance: we use the right-hand side of the equation to compute the solver tolerance, and if there’s nothing there on the right-hand side, that may be a problem.
You prescribe your plume only as part of the temperature boundary condition, not as part of the temperature initial condition. So in the first time step, the only lateral temperature variations are along the bottom of the model. But the material there is constrained to only move horizontally by the velocity boundary conditions.
If you add the term that describes the plume influence ( + dTplume * exp(-((x-xplume)/wplume)^2) * exp(-(y/wplume)^2)) to the temperature initial condition, the problem goes away.
Hope that helps!
Juliane
my guess would be that at the start of the model, there is just no buoyancy to drive the flow, so it’s really hard for the solver to figure out an accurate solution.
I agree, this is effectively the same underlying issue (lack of heterogeneity to drive localization) that lead to similar behavior in the first few time steps of crustal or lithospheric-scale models.
John
Hi John and Juliane,
Many thanks for your replies. I agree that the velocities in that model were low, but even when I add the thermal plume anomaly in the initial condition (as Juliane suggests), and flow velocities are significantly larger immediately from the first time step (up to 1.1 cm/yr), the problem doesn’t disappear:
And the problem doesn’t quickly fade away either: the solution after 3 Myrs (10 timesteps) looks like this (showing the instability remaining in the top right corner):
and after 4 Myrs (80 timesteps):
when strain rates (and velocities reach large values).
Thanks again for your help.
Best wishes,
Jeroen
You are right, I noticed this just now that the oscillations in the strain rate are still there, I just didn’t see them because the starting plume causes a larger signal in the velocity.
I discussed this with @rgassmoeller, and one other idea we have is that this could be caused by the ‘Interpolate output’ functionality. If I set Interpolate output = false, I still see a signal in the strain rate, but it doesn’t show these oscillations any more, and it also becomes resolved more clearly with increasing resolution. This could mean that this is a problem with the interpolate output function we need to look into.
Hi @jdannberg, @gassmoeller, and @jbnaliboff,
Thanks again. I agree that switching off the ‘interpolate output’ reduces the problem. But given the fact that the problem is visible in both the velocity and the strainrate suggests to me that the cause is not in the postprocessing.
I tried a few other things that might give a hint about the causes of this issue:
-
the velocity banding only occurs in the y-velocity, whereas the x-velocity looks perfectly smooth. Since the direction of gravity is an obvious difference between the x- and y-direction, I was wondering if the lithostatic pressure could play a role: most other codes use pressure deviation from the lithostatic pressure in the Stokes solver, but the manual states that ASPECT uses a clever method that allows using the full pressure. Could this be an issue? And if perhaps so, is there a way to instruct the Stokes solver to use pressure deviation from the lithostatic one instead of the full pressure? I tried scaling down the dimension of the model by a factor 1000, which doesn’t make a difference. I also tried reducing rho0 (while increasing thermal expansion to keep the same buoyancy force), which also didn’t make a difference, so perhaps lithostatic pressure isn’t an issue?
-
Also, if y-velocity shows this banding, while x-velocity doesn’t, this should imply that the velocity field is not divergence-free. How is the constraint of a divergence-free velocity implemented in the models? In old codes with direct solvers, the ‘penalty function method’ meant that the flow is never entirely divergence-free, but the user could specify how divergence-free they wanted to have it, but this is probably irrelevant here. I tried explicitly setting the Boussinesq approximation, which didn’t seem to have any effect. I also tried a direct solver instead of iterative one, again without any effect.
-
the horizontal banding in y-velocity (and strainrate) only occurs at depths where the temperature varies, i.e. in and around the thermal boundary layer near the top. If I add a thermal boundary layer at the bottom, then banding occurs there as well. Interestingly, if I choose a linear geotherm all the way from top to bottom of the box, the problem entirely disappears. I expected it to reduce the problem, since the temperature profile is smoother, but not that it disappeared entirely.
Is any of this information helpful?
Once again, thanks for your support.
Jeroen
Hi Jeroen,
Thanks for the tests and observations, especially the observation that the banding disappears for linear temperature helped a lot! I found a solution and preliminary explanation though I am not sure I understand it completely.
First the solution: Switching from a Q2Q1 element (continuous pressure) to a Q2P-1 element (discontinuous pressure) eliminates the oscillations. You can do so using the parameter set Use locally conservative discretization = true in the Discretization subsection.
Preliminary explanation: ASPECT’s default finite element space can only represent continuous linear pressures, which does not ensure a divergence-free velocity when integrated over a single cell (although it does so when integrated globally, see Section 3.2.3. of the first ASPECT paper). This is in particular a problem when the pressure gradient (rho g) changes significantly (e.g. in your thermal boundary layers). I can not give you the exact mathematical explanation, but I would imagine it like this: Since p can not fulfill the momentum equation exactly, the solver ‘makes up’ a velocity that reduces the residual of the equation further, even if that means that the continuity equation is no longer exactly fulfilled. And since (rho g) is purely in y-direction the created velocity is also in the y-component. Because ASPECT does not enforce mass conservation separately, but instead tries to minimize the residual of the combined Stokes block for pressure and velocity this seems like the best solution to the solver.
Now the question I can not answer: I have seen this behavior occasionally in test models, but it was never an issue in my production models, and I can not explain when and why it becomes so prominent. I assume it happens when strong density gradients appear in a model with very little flow in the region of the density gradients, but that is just a guess. Increasing resolution and Stokes degree seem to decrease their wavelength and amplitude, but only Q2P-1 elements completely removes them.
If anyone else has insights I am looking for further explanations. I know that some codes (e.g. PTatin) use Q2P-1 by default, maybe we should consider this as well? We would need to fix issue https://github.com/geodynamics/aspect/issues/714 first though.
Hope that helps
Hi Rene,
Thanks so much for taking the time to look into this further, and for offering a solution for this problem.
I compared the model results for the two cases (for the Q2Q1 and Q2P1 elements), and the model evolution is quite different, with a spurious temperature anomaly entering the model in the Q2P1 case. So perhaps this element might lead to other problems in this model setup (but then again, this model setup was very simple just to illustrate the velocity/strainrate ‘banding’ problem). Do you know of other disadvantages of this element besides a problem with calculating pressure gradients?
However, regardless of this, it is good to know what caused the ‘banding’ problem, and how it can be resolved. And your explanation makes a lot of sense.
Many thanks!
Jeroen
Great, glad that it helped. I am not aware of any problems with this element except for the non-intuitive jumps in pressure at element boundaries. What do you mean by spurious temperature anomaly entering the model at the model boundary? There might be some oscillations at low resolution, but I am not sure why anything in temperature should change when changing the Stokes element (this model has no shear heating, correct?).
Best,
Rene
