Hi all,
Hope you’re all doing good. I have been trying to reproduce a paper “A 3-D numerical study of the thermal evolution of the Moon after cumulate mantle overturn: The importance of rheology and core solidification” (Nan Zhang et al. 2013) but always get different results (see figures at the end, the initial ~40K offset is a printing issue). The authors used CitCom. Their model includes mantle convection and core solidification, and the coupling between the core and mantle is the heat flux at CMB: for a given time step dt, and CMB heat flux Q, the energy across CMB is Q*dt*S which is provided by secular cooling and latent heat of the core (S is area of the core at CMB). They monitored the CMB temperature and heat flux.
Before this post, I used entropy viscosity for energy equation stabilization. Yesterday I discussed with Nan and he reminded me that the difference probably came from the stabilization term in energy equation in ASPECT. So to reproduce the paper, we need to use the same calculation of CMB heat flux Q and temperature as that in CitCom. I searched the forum and found a post answered by Max (SUPG Stabilization - #10 by maxrudolph) and I happened to find the following line in source code ./source/postprocess/heat_flux_map.cc:
const double diffusion_constant =(simulator_access.get_parameters().advection_stabilization_method ==
Parameters<dim>::AdvectionStabilizationMethod::supg) ?
out.thermal_conductivities[q]
:
std::max(out.thermal_conductivities[q],
artificial_viscosity_cell);
which apparently shows the effective thermal conductivity is not literally the one read from material model, a constant. The entropy viscosity will probably not only smooth the temperature (and thus temperature profile) but also change the heat flux. And a recent paper A comparison of 3-D spherical shell thermal convection results at low to moderate Rayleigh number using ASPECT (version 2.2.0) and CitcomS (version 3.3.1) by Euen (2023) shows that entropy viscosity stabilization is pretty diffusive relative to SUPG, esp. for coarse meshes (though they also note that entropy viscosity could be as good or even better than SUP for some selected cases). So my tentative conclusion is that the entropy viscosity suppresses temperature → the temperature gradient near CMB is smaller → less heat flux is pouring into the core → lower temperature of the core. I’m not sure if this makes sense to you, additionally I’m not sure if CitCom uses simply q = -k \frac{\partial T}{\partial z} to calculate heat flux or not (I’ll check this later). Any ideas or comments?
Best,
Mingming
Hi Mingming,
I’m glad to see that someone else is interested in the advection stabilization issue. I tried to reproduce thermal evolution results from CitcomS several years ago and was unsuccessful. In order to solve an advection-diffusion equation in the advection-dominated regime, you will always need to introduce some numerical diffusion. CitcomS and other mantle convection codes use SUPG, which introduces numerical diffusion only along streamlines. This is advantageous for a problem like mantle convection where heat transfer by diffusion occurs across boundary layers where the temperature gradient is nearly orthogonal to the streamlines. ASPECT with SUPG should produce more accurate thermal evolution models than with entropy viscosity. However, my experience has been that using SUPG leads to random convergence failures, often after many thousands of timesteps. Perhaps the stabilization parameter needs to be chosen differently, or a limiter needs to be applied. In CitcomS, there is an option to limit the temperature field (applying a cutoff to keep the dimensionless temperature [0,1] for instance). Another path forward with ASPECT would be to use DG for temperature with the bound preserving limiter but (1) the bound preserving limiter is disabled in spherical coordinates and (2) I’ve seen seemingly random convergence failures with DG, again after thousands of timesteps. It would be great to try to make some progress on these long-standing issues.
For your lunar thermal evolution problem, the Rayleigh number may be low enough that you could disable stabilization altogether. If I recall correctly, Scott King did some benchmarking with ASPECT and found that the results in Zhong et al. (2008) could be reproduced when the entropy viscosity stabilization parameters were set to 0. Doing this for high Ra convection will result in numerical instabilities that eventually cause negative densities and convergence failures.
Best,
Max
Hi Max,
Thanks for sharing your experience. I should have replied much earlier but I have to wait for some tests to finish.
[1] You mentioned you tried to reproduce a CitComS model, may I ask what it is? I also want to take a look and do you want to try it again? What issues did you find in this paper?
[2] Yes, Euen2023 paper said in general SUPG in ASPECT should provide better accuracy. From the tests I have just finished, I found that at least for these cases, SUPG did not deliver better accuracy (all tests show exactly the same temperature profile and CMB temperature evolution). I guess this is probably due to the fact that the temperature gradient is not large enough, and the SUPG does not play a role.
[3] For the tests using SUPG (temperature is P2), the code runs smoothly up to 5000 steps (after which I killed it), I’m not sure if it will fail in another few thousand steps. But this is practically important for freshmen and even for experienced users. I’d like to try DG+SUPG first to see the random failure issue.
[4] For the paper you referred to, I guess it is Grant T. Euan 2023 paper “A comparison of 3-D spherical shell thermal convection results at low to moderate Rayleigh number using ASPECT (version 2.2.0) and CitcomS (version 3.3.1)”, right? I’m working on it now, thanks!
Finally, for the random failure of convergence from (DG+) SUPG, I’d like to make contribution honestly. But I cannot do it now because my supervisor is pushing me recently, besides I’m not a smart guy, so I’m afraid I cannot make any meaningful progress.
Best,
Mingming
