Enhancing Numerical Robustness For Convective Flux In PROTEUS Simulations

Hey guys! Today, we're diving deep into a fascinating issue we've encountered while comparing coupled time evolution simulations in PROTEUS between Spider and Aragog. We've noticed some significant differences in the results, particularly a 'bump' in the melt fraction profile that indicates solidification from above. This is super intriguing, and we need to get to the bottom of it. Let's break it down and see what's going on.

The Curious Case of the Melt Fraction 'Bump'

In our comparative analysis, as highlighted here by @planetmariana, we've observed this peculiar behavior. The melt fraction profile isn't behaving as we'd expect, showing this unexpected solidification from the top down. This is a critical issue because understanding melt fraction is crucial for modeling planetary evolution and thermal behavior. We need to ensure our simulations accurately represent these processes, and this 'bump' is throwing a wrench in the works. The team's hard work in running these complex simulations deserves accurate interpretation, so let's dig into the potential causes and solutions.

This deviation suggests that there might be some numerical instability or sensitivity within our implementation. Specifically, we're zeroing in on the convective flux implementation in Aragog, which we suspect might be more susceptible to numerical precision issues than its counterpart in Spider. You see, the convective flux in Aragog is proportional to the difference between the temperature gradient and the adiabatic temperature gradient, while in Spider, it’s directly proportional to the entropy gradient. This difference in formulation might be the root cause of our problems. We're essentially trying to model how heat is transferred within these planetary bodies, and if our numerical methods aren't up to snuff, we'll get these kinds of artifacts. Therefore, numerical robustness in convective flux is paramount for reliable simulation results.

To put it simply, imagine trying to bake a cake with a faulty oven. If the oven's temperature fluctuates wildly, your cake isn't going to turn out right. Similarly, if our simulations have numerical instabilities, the results won't accurately reflect the physical processes we're trying to model. We're aiming for a smooth, consistent baking process—err, I mean, simulation—so we can trust the outcome. We're not just chasing pretty pictures; we're after scientifically sound conclusions. The numerical discrepancies can lead to misinterpretations about the planetary thermal evolution, potentially affecting our understanding of planetary habitability and geological processes. This is a serious business, guys!

Potential Culprits and Solutions

Okay, so we think the convective flux implementation might be the issue. What can we do about it? Let's explore some potential solutions we've brainstormed.

1. Adiabatic Temperature Gradient Evaluation

One area we're looking at is how the adiabatic temperature gradient is evaluated. Currently, it's derived from thermophysical properties. This can be a bit clunky and might introduce some numerical noise. The adiabatic temperature gradient, in essence, dictates how temperature changes with pressure in a system where no heat is exchanged with the surroundings. If we calculate this gradient using thermophysical properties on the fly, we might be introducing small errors that accumulate and cause issues down the line. So, what's the alternative?

Instead of calculating it dynamically, we could use a lookup table for this quantity. Think of it as a pre-calculated cheat sheet. We'd have a table of adiabatic temperature gradient values for different conditions, and we'd simply look up the appropriate value during the simulation. The hope here is that this will result in a smoother behavior, as we're bypassing the direct calculation, which can be prone to numerical fluctuations. It's like using a pre-mixed cake batter instead of measuring out all the ingredients yourself – less room for error!

Furthermore, we're also considering a correction suggested in issue #66. This correction might address some specific inaccuracies in our current method of calculating the adiabatic temperature gradient. It's all about refining our approach and making sure we're using the most accurate values possible. By implementing these changes, we aim to reduce the numerical noise and achieve a more stable and reliable simulation. We're essentially trying to iron out any wrinkles in our process to ensure a smoother, more accurate outcome. The goal is to represent the adiabatic temperature gradient accurately, as it's a fundamental component in modeling thermal convection within planetary interiors.

2. Smoothing Thermophysical Quantities

Another aspect we're focusing on is the smoothing applied to thermophysical quantities during phase evaluation. When we're evaluating a property close to the liquidus or solidus (the temperatures at which a material melts or freezes), we use a hyperbolic tangent weighing method. This combines the single-phase property and the mixed-phase property. This is a common technique to handle phase transitions smoothly in numerical simulations, but it might not be enough on its own.

On top of this, we're contemplating applying a spatial smoothing—a spatial filter—to the convective flux. This would help ensure a continuous profile and prevent oscillations. Spatial smoothing is like applying a blur filter to an image; it helps to smooth out sharp transitions and reduce noise. In our case, it could help to dampen any oscillations in the convective flux that might be contributing to the 'bump' in the melt fraction profile. Think of it as adding a stabilizer to our cake batter to prevent it from curdling. The goal is to make the simulation more robust and less prone to numerical artifacts.

Spatial smoothing the convective flux is a critical consideration because sudden changes in the flux can introduce instability. We want a nice, gradual transition, and smoothing can help us achieve that. By ensuring that our thermophysical quantities are well-behaved, we can improve the overall stability and accuracy of the simulation. It's all about creating a more continuous and predictable flow of heat within the model. The application of smoothing techniques in numerical modeling is essential for managing the sharp transitions and discontinuities associated with phase changes, which are inherently complex physical phenomena. We're essentially trying to simplify the problem for the computer, making it easier to solve without sacrificing accuracy.

3. Eddy Diffusivity Refinement

The eddy diffusivity is also a key player here. This is a sensitive quantity because it has different scalings depending on the flow regime. In essence, eddy diffusivity represents how effectively heat is mixed by turbulent eddies within the fluid. The scaling laws that govern eddy diffusivity can change dramatically depending on whether the flow is laminar (smooth) or turbulent (chaotic). This means that small changes in the flow conditions can lead to large changes in the eddy diffusivity, which can, in turn, affect the simulation results.

We're thinking about applying spatial smoothing to this quantity, either instead of or in addition to smoothing the convective flux. If the eddy diffusivity is fluctuating wildly, it could be causing some of the issues we're seeing. Smoothing it out could help stabilize the simulation and give us more reliable results. It's like adjusting the sensitivity of a thermostat to prevent it from overreacting to small temperature changes. We want the eddy diffusivity to reflect the overall flow behavior, not just random fluctuations.

The choice of whether to smooth the eddy diffusivity, the convective flux, or both is something we'll need to investigate further. Each approach has its own potential benefits and drawbacks. We might find that smoothing one quantity is more effective than the other, or that a combination of both techniques yields the best results. It's all about experimenting and finding the right balance to ensure the stability and accuracy of our simulation. Therefore, spatial smoothing on eddy diffusivity is a key area of focus for improvement.

Beyond Convective Flux: Addressing the Dilatation Source Term

It's worth noting that other terms in our simulation might also benefit from spatial smoothing, particularly the dilatation source term. This term, which describes the rate of volume change, can exhibit stiff and sharp behavior. This is like having a sudden burst of expansion or contraction in our model, which can cause numerical instability. These behaviors are detailed here.

The dilatation source term is particularly tricky because it's directly related to changes in density, which can be quite abrupt during phase transitions or under certain pressure conditions. If this term is fluctuating rapidly, it can introduce significant noise into the simulation and make it harder to converge to a stable solution. Think of it as having a hiccup in the system – it disrupts the smooth flow of calculations. Smoothing this term could help to dampen these fluctuations and lead to a more stable and accurate simulation.

By applying spatial smoothing to the dilatation source term, we can effectively filter out some of the high-frequency noise and prevent it from propagating through the simulation. This is particularly important in regions where the density is changing rapidly, such as near phase boundaries. The goal is to capture the overall trend without being overly sensitive to small, localized fluctuations. Just as a conductor manages the orchestra, ensuring each instrument contributes harmoniously, spatial smoothing helps manage the numerical terms, allowing the simulation to proceed smoothly and accurately. Therefore, addressing the dilatation source term is crucial for improving the numerical stability of our models.

Next Steps and Conclusion

So, where do we go from here? Well, we've identified a few key areas to investigate. We need to test these potential solutions—the lookup table for the adiabatic temperature gradient, the spatial smoothing of convective flux and eddy diffusivity, and the handling of the dilatation source term. This will likely involve running a series of simulations with different configurations and comparing the results.

Our goal is to eliminate this 'bump' in the melt fraction profile and ensure our simulations are robust and reliable. This is crucial for accurately modeling planetary evolution and understanding the thermal processes that shape these fascinating worlds. We're not just tweaking numbers; we're striving for a deeper understanding of how planets work. This is the core of our mission, guys! Let's keep digging, keep testing, and keep pushing the boundaries of what we can simulate. We're on the verge of some exciting breakthroughs, and by addressing these numerical challenges, we're paving the way for more accurate and insightful models of planetary interiors. Keep an eye out for updates as we progress! This adventure into enhancing numerical robustness is a testament to our commitment to high-quality science and the relentless pursuit of accurate planetary models.