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
Hi Max,
I have successfully benchmarked against Euen2023 and Zhong2008 paper. However, these papers are using classical Rayleigh-Benard convection. In modeling realistic problems, there could be a lot of extensions where our benchmark can diverge. Benchmark is no easy task for application papers, esp. those having errors, typos and unspecified conventions or assumptions.
Regards,
Mingming
Hi Mingming,
It’s good to hear that you can reproduce the benchmark values from Euen et al. 2023 and Zhong 2008. The Euen et al. 2023 models are pretty low Ra, several orders of magnitude smaller than Earth. We have still been unable to reproduce the thermal evolution of high-Ra (2x10^8 using Earth’s radius as the length scale) in older CitcomS models (Rudolph and Zhong, 2014). At higher Ra, one cannot disable the stabilization parameters and using SUPG or DG for temperature produces convergence failures after many thousands of timesteps.
Max
Just want to add, for myself and new comers who want to reproduce this paper: in Material model, you have to use project to Q1 only viscosity for Material averaging entry, otherwise the viscosity range would be smaller than the given value and thus Vrms (mean), T (mean), |Nu| will be larger.
Hi Max,
That is an interesting paper. As you said before, the bound preserving limiter is not available in spherical shells, maybe we could write a simple plugin to cut off the temperature? For SUPG or DG failure, what is the failure exactly? Have you tracked the evolution of the range of all fields and residuals? I’m curious and I’d like to try later.
Regards,
Mingming
