[aspect-devel] Problem with iterative advection solver

[Luuk Schuurmans]

Hi John,

Thanks for your response! From the start, I was using a low CFL number, to
that probably has not caused my problem, of which I think it indeed is
different from the problem of your student. At the moment I am trying to
work out a different initial temperature profile, hopefully this will give
a new insight in why the temperature solver crashes.

Concerning the DG options for composition and temperature, I do not use
them, so that will not give any problem.

Again, thanks for your response, and hopefully I can fix the problem soon!

Cheers,

Luuk

2017-12-18 22:20 GMT+01:00 John Naliboff <jbnaliboff at ucdavis.edu>:

> Hi Luuk,
>
> Sorry for the delayed response and glad you found a solution with the
> sticky air!
>
> A student working with me has been running ocean-continent subduction
> models with a free surface. As far as I know this problem has not occurred
> so far. We are also using the vertical surface advection option.
>
> At some stage we had convergence issues with the compositional field
> advection, but this problem appeared relatively early on in the simulation
> (< 2 Myr). We resolved the issue by going to a smaller CFL number (0.5 —>
> 0.1 or 0.2). However, it sounds like this might be a different issue than
> what we saw? I’ll check with the student to see if minimum temperature
> errors ever appeared.
>
> Are you by chance using the DG option for the composition or temperature?
> Anne Glerum and I have both found errors when using the DG Limiter with a
> vertical free surface advection scheme. This is on the list of things to
> fix.
>
> Cheers,
> John
>
> *************************************************
> John Naliboff
> Assistant Project Scientist, CIG
> Earth & Planetary Sciences Dept., UC Davis
>
>
>
>
> On Dec 12, 2017, at 5:35 AM, Luuk Schuurmans <
> l.f.j.schuurmans at students.uu.nl> wrote:
>
> 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
>>
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20171221/fe96e76e/attachment.html>