Numerical stability.

[Search-Enemy]

Hello professor boris, I have another question to ask you. Because I just read the source code of Example2D.jl and related functions, I haven’t had time to look at the MTK_GMG examples yet.

My question is: CFL criteria mentioned in the literature to limit dt:

If you only keep kappa on the left, it is equivalent to limiting a maximum value for thermal conductivity. In Example2D.jl, you provide the setting of dt:

Then you can calculate the maximum upper limit of kappa, and the calculation is 2.040816326530612e-6

Then, the actual kappa of the model only needs to meet this condition kappa<2.040816326530612e-6 to ensure numerical stability.
However, when I take the values Kc = 2.899, Cp = 980, Rho = 2800,

Now, kappa = 1e-6 (approximately) ,

this is obviously within the convergence conditions, but Example2D running in this way has a large negative temperature anomaly in just 5kyr.

Of course, I not only tested this set of values. I also used custom thermal conductivity temperature-dependent parameters:
K[i,j] = 0.70 + 770/(350+T[i,j] - 273.15),K in range(1.3,2.9)(approximately) when T is 273.15-1173.15℃ [Clauser,1995]
When Cp is 980 and below (such as 700, 800, 900, 950, etc.), numerical instability will occur, and the time for the abnormal value of T to appear increases as Cp increases. When Cp is 1000 and above, stability can be maintained for a longer period of time.(within 90kyr). At first I thought there was a problem with the temperature-related functions and thermal diffusion I wrote myself.
Later, I tested the MTK source program Example2D.jl (just like the set of values ​​exampled above), and encountered Same situation.
Therefore, what is the reason for the situation that the parameters comply with the CFL criteria, but the model runs unstable?Maybe you can give me some important tips.
Thank you very much.

[boriskaus]

Not really sure but I guess it may be related to the effect of latent heating, so you can try switching that off. The effect you are seeing are indeed the usual blowups.

Please note that you don't have to program temperature-dependent conductivity yourself - we have a range of those already available in GeoParams so you can just select one.