[aspect-devel] Problem with iterative advection solver

Fri Dec 22 05:51:09 PST 2017

Hi Luuk,

For reference, we are juxtaposing a 70 Myr oceanic plate (plate cooling model geotherm) against a continent with a surface heat flow of 55 mW/m2 (leads to roughly the same plate thickness). 

Looking back through your original message, another difference is that we are running incompressible simulations (Boussinesq, no shear or adiabatic heating). 

Sorry, none of this is likely very helpful. If we can try something out as an extra test case please do let us know!

Cheers,
John

*************************************************
John Naliboff
Assistant Project Scientist, CIG
Earth & Planetary Sciences Dept., UC Davis

> On Dec 21, 2017, at 9:24 AM, Luuk Schuurmans <l.f.j.schuurmans at students.uu.nl> wrote:
> 
> 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 <mailto: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 <mailto: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 <mailto: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 <mailto:Aspect-devel at geodynamics.org>
>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel>
> 
> _______________________________________________
> 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/20171222/f4065991/attachment-0001.html>