cd2010 viscosity dalphadt update

[RoguePotato]

In Headers/Sph.h, at the end of the ComputeCullenAndDehnenViscosity() function, parti.dalphadt is always being set, but should only be when alpha_loc > parti.alpha (see section 3.2 of Cullen & Dehnen, 2010).

[rbooth200]

Have you noticed whether this makes much of a difference? You could try replacing line 456 with

parti.dalphadt = (FLOAT)0.1*parti.sound*(max(alpha_visc_min, alpha_loc) - parti.alpha)*invh;

which I think should give the intended behaviour.

[rbooth200]

The fix I made was by moving the if-statement on line 453 to line 455 and putting the update within the if-statement. This keeps part.alpha = 0 unless the particle is in a pre-shock region (the check on line 447).

I'm not quite sure what you mean, are you setting alpha to 0 when not in a shock region? I don't think this is the correct behaviour. I think the fix I proposed is the correct one, it ensures that alpha decays exponentially to the desired new value.

Prior to this change, a low-mass disc was warmer for cd2010 viscosity vs mm97 viscosity. Making the change results in a cooler disc for cd2010, which makes sense, as there is less viscosity and less viscous heating.

Ok, can you open a pull request since it seems imporant.

[RoguePotato]

Sorry, that previous comment was posted by mistake, the fix I proposed was wrong. I'm going to run some more tests with your change and see what happens.

[RoguePotato]

Okay I've noticed a few things. I'm using my own initial conditions from a file, and before the simulation runs, all the particle alphas are set to alpha_visc, because alpha_loc is incorrectly high (line 449). After the simulation starts running, alpha_loc seems to have reasonable values, but part.alpha never changes thereafter in the cd2010 function because the if-statement on line 453 is never met. One would expect however, that the particle alpha would evolve due to dalphadt being set. After noticing that the viscosity remained at alpha_visc for a long amount of time, I did some investigating. Normally, the flow is:

1. Zero values before calculation parti.dalphadt = 0 (performed here))
2. Calculate part.dalphadt based on viscosity method.
3. Update such that parti.alpha += parti.dalphadt * dt (performed here)

However, the placement of the cd2010 function means that the dalphadt is reset to zero before alpha is updated. This means that this line has no effect on the code, because it gets zeroed before it gets used.

Also, where does the factor of 10 come from on line 448?

[rbooth200]

Ok, so I think you're close to getting to the bottom of this. I've just checked, and dalphadt should not be set to zero in Sph::ZeroAccelerations. The ordering meant this did not affect anything for the M&M switch, but as you pointed out its wrong for the C&D switch, so you can remove it (your step 1). Does removing this line fix the problem?

The C&D implementation we use is similar to phantom one. The factor of 10 is just a scaling factor and I haven't really played around with it so much, but both phantom and Rosswog (2015) use a similar factor. Reducing it might cause problems at weak shocks though.

[RoguePotato]

Yeah this gives expected results now, the cd2010 viscosity is generally lower than the mm97 method.

[rbooth200]

Ok great. Can you open a pull-request with these (and only these) changes?

@giovanni-rosotti: once this bug is fixed we should re-run the viscous spreading test for the paper. Can you do this and update the overleaf?

[rbooth200]

Fixed in #162