[aspect-devel] Problem with iterative advection solver

Tue Dec 12 02:35:28 PST 2017

Hello everyone,

I have been working on my project to try to find the cause of the
aforementioned problem with the iterative advection solver (for those who
did not see the initial problem statement of the 24th of November, I added
it at the end of my email). It became clear that the crash is probably
caused by the implementation of a free surface. During the last tests, I
added a sticky air layer on top of the modelling domain to check whether
the free surface caused the problems. Running models with the sticky air
did not cause the iterative advection solver to crash, and also the strange
temperature fluctuations disappeared. Therefore, I suspect that the models
which make use of the free surface try to use values from outside the
modelling domain at certain timesteps. As a consequence, the model crashes
due to this obtained undefined value. This does not happen in the case of
sticky air, since this material flows out through an open top boundary.
Does anyone know how to get rid of the free surface problem, without having
to use the sticky air layer on top of the domain?

If you need any further information, do not hesitate to send me an email.
Thanks in advance!

Kind regards,

Luuk Schuurmans

Initial problem:
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. The temperature
at the top boundary is prescribed to be 273 K, but at certain timesteps the
minimum temperature in the modelling domain is quite far below this value
of 273 K (see Temperature timesteps.png). From the paraview files it
becomes clear that the cell which is exactly at the interface between the
subducting plate and the overriding plate corresponds to this lower than
expected temperature (see temperatureerror.png and velocities.png). Note:
the paraview screenshots correspond to timestep 196, at which the heaviest
anomaly occurs in Temperature timesteps.png. What is remarkable, is that
the minimum temperature recovers after a few timesteps. However, for the
moment of crashing, we presume the error in the temperature to be too large
to recover, resulting in a -nan and a crash of the model.

I have looked at this problem with Menno Fraters, but we were not able to
figure out the exact cause. We presume that the error is caused by
something in the free surface plugin, which possibly can not handle the
interpolation in the case of a big temperature gradient in combination with
opposing velocities at the boundary between the two compositions (possibly
it tries to use material from outside the modelling domain and fails to
assign good values to it?). We use a vertical surface anomaly projection
for the model, and the prescribed convergence velocities are in the order a
few of centimeters/year.

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.

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.

Our question now is, do you maybe know the exact cause of the crash (is it
in ASPECT, or is it caused by dealii, and could it be caused by the thing
we mentioned) and how can it possibly be solved?

Thanks in advance for your collaboration!

2017-11-29 16:28 GMT+01:00 Luuk Schuurmans <l.f.j.schuurmans at students.uu.nl>
:

> Dear Wolfgang,
>
> First of all, thanks a lot for your reply! Since I was not subscribed to
> the mailing list, initially I did not receive your email, but now I
> subscribed to the list. I discussed your proposed solution with Menno, and
> tried to run the model with all heating models turned off.
>
> For the screenshots I sent to you in my previous email, I used an
> incompressible model. However, since I am using the extended Boussinesq
> approximation, both shear heating and adiabatic heating were turned on. In
> my previous mail I already showed that turning off the shear heating did
> not solve the problem, and unfortunately switching off the adiabatic
> heating also did not have any effect on the problem. The strange low
> temperature peaks still are present, but like the model without shear
> heating, the model run did not crash. What is remarkable, is that the
> errors in the temperature seem to occur at or around the same (random)
> timesteps, e.g. 18 (for model 1, 2 and 3) and 194 (model 3), 195 (model 2)
> and 196 (model 1). To clarify, model 1 corresponds to the model with both
> adiabatic and shear heating on, model 2 to the one with only adiabatic
> heating, and model 3 to the one with both heating models switched off.
>
> When the value of the composition field is analyzed, it becomes clear that
> again values slightly larger than 1 are obtained for model 3 around the
> location of the erroneous temperature, which still is the boundary between
> the subducting and overriding plate. The value is definitely highest (e.g.
> 1.09) at the problematic gridpoint, but also at the ones where there seems
> to be no temperature problem the composition value can be slightly larger
> than 1 (e.g. 1.003).
>
> As a conclusion, we still do not know the exact cause of the problem.
> Hopefully you can make something of this additional information. If you
> need any further information, please let me know!
>
> Kind regards,
>
> Luuk Schuurmans
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20171212/cf1562c4/attachment-0001.html>