Draft: Streamlined MAGMA2 implementation

This is a reimplementation of the MAGMA2 branch with:

• Simpler symmetric matrix handling
• Only 1st order reconstruction of v and u.
• Only basic viscosity and diffusion terms as in the Rosswog 2020 implementation

Questions:

• Do we need more advanced terms? SPHENIX-like? Or the Rosswog entropy-based viscosity switch?
• What implementation do we want? Traditional or Gasoline?

Todo:

• Document the matrix stuff
• 1D and 2D cases
• Scale-factors
• Strange corner cases (e.g. no ngb)
• Tests in galaxy settings
• Properly implement the symmetric functions.
• Fix documentation. It still mentions Minimal and SPHENIX in places.

added SPH label

@jborrow @dc-sand2 would be good to have your thoughts on the points above.

I am not fully convinced by the Gasoline-like implementation. I would say to go without it... In Zhen's tests did this not cause problems? I remember when we swapped to the 'regular' version everything suddenly worked. To me, I worry that kind of EoM leads to 'blobbing' as particles with vastly different phase-space properties can easily inter-mix.

As for advanced terms, let's try some galaxy formation simulations with the basics first. No point in fixing something that's not broke (or rather, it's probably easier to debug weird things that happen in those settings if you don't have fancier switches). For the entropy switch, do we not need to consider how sub-grid cooling impacts it?

I was thinking of the tests from your paper on hydro+feedback combination.

Yes, those would certainly be helpful. It is interesting to know how much energy is leaked (iirc for GIZMO it is 50%!). This is obviously less important for kinetic feedback, though.

The viscosity reconstruction is (unhelpfully) too smooth in cases where we actually want viscosity e.g. shock while still being overdissipative in shear flows e.g. KH. In particular the overdissipation in the KH tests I've been running has been pretty dominant for the evolution (when it should be fully switched off) to the point where I'm not sure whether I trust the results for instability growth from Rosswog 2020 or Frontiere et al. 2017 (making things slightly viscous will make it look nice and laminar, but not be accurately simulating an inviscid fluid). I've found with the hydro scheme that I've been working on, that increasing the alpha and beta to 2 and 4 respectively helps reduce the ringing in shocks and that multiplying Q by a Balsara switch prevents overdissipation in my KH simulations, even with these higher constants. Note also that with the Q viscosity formalism and velocity reconstruction the alpha, beta aren't really directly compatible to the alpha, beta we're used to (see Frontiere et al. 2017 for some weirder choices for these). A side note that might be relevant: In my hydro scheme I find that using a normal Balsara switch is actually underdissipative as you evolve to get small-scale compression between the shearing layers in the KH so I've been testing using a pairwise divergence estimate to make a pairwise Balsara switch (similar to the pairwise divergence estimate used in the viscosity in the first place) and this seems to work really well. This might be irrelevant if there are other dominant effects at play e.g. using normal, smoothed SPH density estimate, but if you end up with something looking like a fingering instability between the shearing layers on the particle lengthscale, then it's probably the same problem.

In my tests I've consistently been finding way better KH evolution with the Gasoline-style EoM. Note that this is with different gradients and density estimates so not a one-to-one comparison to the magama2 scheme. It might be the case that smoothed initial densities, equal particle spacings (rather than equal mass) and too much viscosity (as discussed above) in the Rosswog hydro tests hide some of the differences, but I'm not sure. It could also be that with enough artificial conduction it wouldn't matter. I reckon it would at least be worth testing again to see if e.g. including a Balsara switch gives a bigger difference between the normal and gasoline-like schemes. Something to note with the gasoline-like scheme: Wadsley et al. 2017 have those "f" terms in the gradients which are not in the Rosswog version. My understanding of these is that they're a correction introduced because the new internal energy evolution equation and the normal SPH density estimate are not thermodynamically consistent (see section 3.2 in Read et al. 2012). I don't know how big of an effect not including something like that has (maybe comparatively small given that artificial dissipation is breaking thermodynamic consistency anyway and there aren't any grad-h terms?). At one point Sergio was testing a normal gasoline-style scheme for our planetary simulations and I seem to remember that the "f" terms were causing some problems at vacuum boundaries, so maybe something to avoid if it can easily cause weird things to happen even if you don't have any vacuum boundaries in your simulations?

This is all from the perspective of improving purely the hydro, so I have no idea how any of this is interplaying with e.g. feedback

changed the description

@dc-sand2 what signal velocity for diffusion were you using in your test when saying that the GDF-like EoMs perform better on KH tests?

The gravity version from Rosswog 2020 (which is significantly less than the non-gravity version for the same alpha)

I've been doing the KH in 3D (similar to Rosswog 2020), but with equal mass particles and without a smoothed initial density. I still keep the smoothed initial x-velocity from Rosswog 2020, so that there's less conduction. In these tests I've also been using a Wendland C6 kernel with eta=1.866 for validation.

But the main thing I've been focussing on is getting good behaviour with little/no conduction across the boundary (i.e. for a material boundary where we can't have conduction between different materials), so not sure how things will behave if you're happy to have more conduction. Though it could be that the different density estimates and gradient estimates might make things behave differently enough from my tests that the amount of conduction isn't the key difference.

@yvesrevaz your NFW hydrostatic test could be interesting with this branch too.

@matthieu I'm preparing the test:

could you put : https://drive.switch.ch/index.php/s/mXrYj7n9XRR5alg here: http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/NFW_Hydrostatic/nfw.hdf5

The test is done and ready to be merged with the master (!1807 (merged)). It will be simple to use it with the MAGMA2 solver.

If anyone tests this, keep an eye on particles near the vacuum boundary as well. We found the need to switch to a normal, spherically symmetric kernel gradient function near vacuum boundaries as magma2 gradients can send particles flying outwards. Something to keep in mind if the intended use of this hydro scheme has sharp vacuum boundaries.

added 28 commits

• 9475c3a5...ce511ffa - 26 commits from branch master
• 1c1edf7e - Merge branch 'master' into MAGMA2_matthieu
• 22fd117e - Applied formatting script

Compare with previous version

marked the checklist item Document the matrix stuff as completed