[aspect-devel] Problem with iterative advection solver

[Wolfgang Bangerth]

Hi Luuk,
welcome to ASPECT!

> I am Luuk Schuurmans, a master student from Utrecht University, and at the 
> moment I am writing my master thesis on subduction modelling, for which I use 
> the ASPECT code. During one of the initial tests, in which I try to simulate 
> the subduction of a oceanic plate below a continental plate, I adjusted 
> (lowered a bit) the value of the activation energy for dislocation creep of 
> the overriding plate crust. When I ran the model with a composite viscosity, 
> after ~1000 timesteps I ended up with the error message that the iterative 
> advection solver did not converge, and that the residual in the last step was 
> -nan. The time steps around this time have all been equal, because they are at 
> the maximum time step set in the input file. For a higher value of the 
> activation energy, the model did not crash around this timestep.
> 
> The statistics file revealed that there is a strange thing going on with the 
> minimum temperature of the model at certain timesteps.

Good detective work already!


> I have also tried to run the model with shear heating turned off, since the 
> viscosity terms plays and important role in that factor. The model just runs, 
> and a crash of the advection solver does not occur. However, as you can see in 
> noshearheating.png, the problem with the low temperatures also occurs many 
> times in this case, but without a crash as result.

The fact that the model doesn't crash seems like coincidence. The temperature 
still looks wrong, and whether or not the solver succeeds, figuring out the 
reason for the temperature issue sounds like a worth-while goal.


> Besides the error in temperature, it also seems that the composition is 
> assigned a value which is not totally correct (see compositionerror.png, 
> composition should be between 0 and 1). The location where this occurs is the 
> cell directly beneath the one where the temperature anomaly is present.

That too is a great observation.

Here's a wild shot into the dark: Are you using a compressible model? At the 
location of the problem, the velocity vectors vary wildly. I would not be 
surprised if in a compressible model, you have a large negative divergence 
here because you simply don't resolve the dynamics of the flow with your mesh. 
If there is a large negative divergence, then in compressible models you will 
have adiabatic expansion -- cooling! -- happening at these locations. If the 
neighboring cell has a large positive divergence, then you will have 
compression that might lead to a larger value of the compositional field.

I'm just speculating, but at least the first of these two could be tested by 
removing adiabatic heating from the list of heating models (somewhere in the 
input file).

Best
  W.


-- 
------------------------------------------------------------------------
Wolfgang Bangerth          email:                 bangerth at colostate.edu
                            www: http://www.math.colostate.edu/~bangerth/