Problem with solving temperature after inverse method

Extension of Elmer in computational glaciology
ygong
Posts: 15
Joined: 08 May 2013, 16:51
Antispam: Yes

Problem with solving temperature after inverse method

Post by ygong »

Hi,

I am trying to do the inverse method run in this way : (1) first relax the topography ; (2)then restart from the relaxation and do the inverse method run with an initial guess for temperature; (3) then restart from the previous results and add temperature solver to calculate temperature; (4) then do the inverse method again.

I have got these errors when I do (3) if I restart from the relaxation which is shorter than 50yrs:
Max/min values Temperature: 574.677894395332 / -3.64321112064340
EffectiveViscosity: Positive Temperature detected in Glen - limiting to zero!
...
TemperateIceSolver (temperature): Assembly done
ERROR:: ComputeChange: Norm of solution appears to be NaN

Well, the temperature solver is happy if I restart from the longer period of relaxation (>50yrs), but I think the topography might be relaxed too much.
I have tried to use a small timesteps (1/52, weekly) in my relaxation,but it didn't work.
I have tried to unfreeze the vertical update of the mesh of side wall by comment "Mesh Update 3 = Real 0.0" out in the boundary condition it still doesn't work

Then I found out that the only huge differences between the results of (2) is the mesh velocity (attached a pic ). It seems something wrong with the mesh update.
Dose anyone have any idea how to solve this problem? I attach my sif for relaxation here.
Attachments
ASF_relaxation5yr.sif
relaxation files
(15.08 KiB) Downloaded 925 times
inverse method results.png
The mesh velocity results after doing (2) restarting from different relaxations
(719.55 KiB) Not downloaded yet
fgillet
Posts: 46
Joined: 30 Sep 2010, 16:58

Re: Problem with solving temperature after inverse method

Post by fgillet »

Hi,

I think there is no problem with your relaxation.
When starting from a given observed topography, the incompressibility equation usually leads to very high and unphysical ds/dt (i.e. vertical mesh velocity).
This is why relaxation is needed, and at the beginning there is very high mesh velocities and the maximum value is decreasing as you relax the surface and tends to a steady state.

I think the problem with the temperature might come from the fact that the mesh velocity is used to correct the advection term in the temperature solver (ALE method); so that the high and unphysical mesh velocities when you restart your solution can affect your temperature solution.

Maybe the best is to reset the mesh velocity to 0 after the restart and before the temperature computation.
You can write a solver to do the job but maybe someone has a better solution to reset the mesh velocity to 0?
An alternative would be to put a flag in the Temperate ice solver to switch off the ALE method?

Fabien
rgladstone
Posts: 64
Joined: 15 Apr 2013, 16:23
Antispam: Yes

Re: Problem with solving temperature after inverse method

Post by rgladstone »

Fab, thanks for your reply. Would it work to set the following in the simulation section?

Restart Before Initial Conditions = Logical True

Then in the initial conditions section set mesh velocity to zero? Can the restart variables be overwritten like this?

Cheers,
Rupert
fgillet
Posts: 46
Joined: 30 Sep 2010, 16:58

Re: Problem with solving temperature after inverse method

Post by fgillet »

Hi,

I didn't know about this keyword; Thanks Rupert
If it works as expected that would be an elegant solution.

cheers

Fab
ygong
Posts: 15
Joined: 08 May 2013, 16:51
Antispam: Yes

Re: Problem with solving temperature after inverse method

Post by ygong »

Hi, Thank you very much, Rupert and Fabien! I have set my simulation run again using the way mentioned by Rupert.

As what was said by Fabien, I wonder whether it means that the relaxation is not long enough if ds/dt is still rather high after a couple of years.

Cheers,
Yongmei
fgillet wrote:
When starting from a given observed topography, the incompressibility equation usually leads to very high and unphysical ds/dt (i.e. vertical mesh velocity).
This is why relaxation is needed, and at the beginning there is very high mesh velocities and the maximum value is decreasing as you relax the surface and tends to a steady state.

Fabien
Martina
Posts: 39
Joined: 17 Apr 2013, 14:06
Antispam: Yes

Re: Problem with solving temperature after inverse method

Post by Martina »

Hi,

I didn't had Yongmei's problem when doing exactly the same simulations on another ice-cap, this puzzles me a bit, but didn't had yet the time to look deeper in it.

The "Restart Before..." will do the job, I've used that for this purpose in another occasions, you just have to be careful in setting your initial conditions (i.e. not putting sth to zero that should be read from the restart).

I had in similar simulations the observed problems with the temperature solver, and I could solve it by carefully setting my initial conditions for the temperature, but since your observations depend on the length of your relaxation run, it's maybe not the solution here? What are your inital conditions for the temperature?

Other temperature problems occured to me in the case of problems with the MeshUpdate Solver, so I'd recommend to really make sure that everything is well in your relaxation run (for example Depth>=0 and good convergence behaviour of the FreeSurfaceSolver).

Also an option would be to try directly the relaxation with sliding, maybe you end up with something strange by doing it first without sliding.

Martina
Martina
Posts: 39
Joined: 17 Apr 2013, 14:06
Antispam: Yes

Mesh velocities, ALE method in the temperature Solver

Post by Martina »

Hi all,

I'm using a very similar set up as described in the initial post, so far with no real convergence problems. By curiosity I tried now Fabien's suggestion to set the Mesh Velocities to 0 after relaxation before calculating temperature fields.
I observe very different temperature fields if I do so (it is much colder when putting the Mesh Velocities to 0). The temperature field that I get when setting Mesh Velocities to zero corresponds closely to the one calculated without relaxation, so setting the Mesh Velocities to 0 seems to take away the "benefit" from the relaxation (if it is a benefit and not an error).

Now I'm wondering what is right: should I set Mesh Velocity=0 everywhere since the Mesh Velocity rather introduces errors? Is it wrong to set the Mesh Velocities to zero since they are important? Is there a possibility to switch off the ALE method as suggested by Fabien (I couldn't find such a flag so far in the documentation) to make a comparison? Does it imply that my relaxation run is not well defined since setting the Mesh Velocities to 0 or not should lead roughly to the same solution?

I would be happy if somebody could give me some more insights in the role of these Mesh Velocities in the context of relaxation and advection to find the right way calculating temperature fields.

Thanks in advance,
Martina
Martina
Posts: 39
Joined: 17 Apr 2013, 14:06
Antispam: Yes

UPDATE on Mesh Velocity

Post by Martina »

Hi all,

after discussion with Fabien and Thomas here a quick update:

- There is no flag to turn off the ALE method, and probably there is no need for it since it can be done with the simply initial condition "Mesh Velocity =0".

- Mesh Velocities result from the transient relaxation.
When computing the steady state temperature fields (implicitly assuming a steady-state free-surface elevation), MeshVelocities should always be set to 0 after restart.
For further transient simulations, with arbitrary moving meshs, ALE is the correct method and has to be used including the MeshVelocities.

Martina

P.S. In the meanwhile I have some simulations (restart from restart...) where setting MeshVelocity=0 even in combination of "Restart Before..." doesn't work (no idea why, ElmerPost says MeshVelocity is -inf, ParaView says MeshVelocity doesn't exist). For the moment I use a modified Temperature Solver that sets the MeshVelocities to zero and not to their actual values. This leads to small modifications in the temperature in some simulations and in others the simulation doesn't converge at all if I don't use this modification.
denis.cohen
Posts: 56
Joined: 15 Dec 2010, 13:50

Re: Problem with solving temperature after inverse method

Post by denis.cohen »

Hi Everyone,

I am writing to this post instead of opening a new one since I think this is a related problem.

I am solving a paleo ice-flow problem, trying to reconstruct part of an ice field at the LGM so I have no temperature or ice velocity data, only a bed and ice surface DEM.
First, I start with an initial guess for temperature and solve Stokes equations with the free surface keeping the temperature constant. This works fine if no sliding but with sliding I have, at some nodes, sliding speeds > surface speeds. This happens where there are large gradients in bed and/or ice surface topography.
Then I do a restart and solve for the temperature only, keeping the ice flow and free surface (from the previous run) constant. If I solve for the steady state temperature, this only works if I have no sliding. If I have sliding the code does not converge at the first steady state temperature iteration probably because the velocities at the bed are not very good. If I solve for the temperature in transient mode no problem, the code converges at each iterations but it takes too much time to reach a steady state for temperature.

I am wondering if anyone had similar experience they could share and ways of solving this. I've tried several mesh refinements to better resolve the bed topography without much success. Smaller time steps does not help either, at least to the extent I have tried. I can give more details if needed.

Also I have tried without success to implement and use a temperature-dependent sliding law. The code always crashes after some time iterations. Has anyone experienced this?

Thanks

Denis
Martina
Posts: 39
Joined: 17 Apr 2013, 14:06
Antispam: Yes

Re: Problem with solving temperature after inverse method

Post by Martina »

Hi Denis,

I didn't had exactly your problems, but I'll still share some quite some temperture-convergence-issus, maybe that inspires you.

My steady-state temperature fields converges well in most cases, however it depends a lot on the initial conditions. Did you try different ones (for me the best is a linear depth dependency like a steady-state profile at the ELA)? When I run in transient I managed to do around 50years fully theromechanical coupled, then it crashes and so far I have no clue why.
I can also say that sometimes enforcing the conditions T<T_pm strictly (including Loop While Unconstrained Nodes = Logical True) does not work out. Somehow Elmer can't solve the problem with the strict condition. Did you try to be less strict on that, maybe even putting a constant limit of 273 or 274 to see if it helps? That does help me in some cases where I can't find a converged steady-state otherwise.
Did you isolate your different heat sources (advection, strain, friction, heat flux...) to see to which one the problem relates and to check if all of them are correctly calculated? In my case it's always friction heat which is the problematic one. You could also imagine to take the temperature field (if you can get that) without friction heat as initial condition for the run with it or that kind of setups. I had at some point plenty of runs with different restarts, including doing some timesteps in transient and taking that as initial conditions and so on. For me it was possible to calculate a temperature-only-steady-state by doing a transient run and I took that as check to see if my directly calculated steady-state was finally right, but if I understood you right, you can't do that, even not with constant velocities?

I have also tried a temperature sliding dependent sliding law, which did not work out, but it's a long time ago, and never went back to it with my newer setup. Are you aware of the Greenland paper where they did it in Elmer? Is your sliding law similar to that? I do not know enough details about that work, but I was wondering if the fact that it worked for them was linked to their well defined initial condition for the temperature field (but that's really only guessing from myside).

One last thing, I imagine you checked at each timestep that you calculate your temperature steady-state (sorry, I haven't quite understood at which moment of your simulation you are doing that) that your mesh is well defined, i.e. now surface points below the bed? I had that earlier creating a lot of troubles with the temperature solver and finally switched to the Internal Mesh Extrusion which avoided that problem.

Good luck,
Martina
Post Reply