gandalfcode / gandalf

GANDALF (Graphical Astrophysics code for N-body Dynamics And Lagrangian Fluids)
GNU General Public License v2.0
44 stars 12 forks source link

RadWS Fix #159

Closed RoguePotato closed 7 years ago

RoguePotato commented 7 years ago

Fixed Radws EoS integration scheme. Requires EoS table with gamma. Added gpot_hydro to ignore sink potential.

The EoS pressure function is not inherited by the polytropic nor the radws EoS regimes. For now I have change the default pressure function to use part.gamma.

dhubber commented 7 years ago

A couple of quick comments. First of all, congratulations for making your first pull request :-) . Second, without going through the changes in a huge amount of detail, we spent a while optimising the code for both serial and parallel performance mainly by reducing or optimising the memory usage. We did this largely by removing as many redundant fields from the Particle and TreeCell data structure. However, here there are two new fields. I would suggest these should only be added if they are only 100% necessary and cannot be computed some other way (e.g. by a function using other particle properties). So, if you can think of a way to get rid of them then great, but if there is no way to do this then it's fine to leave them.

RoguePotato commented 7 years ago

The gpot_hydro field will be able to be removed if we implement the Lombardi variant of the RadWS scheme (Lombardi et. al., 2015), as this uses pressure instead of gravitational potential. Maybe it is possible to remove it now, and within the RadWS scheme, deduct the gravitational potential of stars within the simulation (if it is that simple anyway).

The particle gamma may need to stay as it is set in Hydrodynamics/EnergyRadws.cpp but used in the EoS implementation i.e. Thermal/RadwsEOS.cpp. The pressure calculation cannot be done within EnergyRadws.cpp as the only opportunity to do so is at the start or at the end of the timestep. The eos->Pressure call is in between the two and would overwrite the calculation.

The main concern right now is to why the EOS.h Pressure function is not inheritable, and whether this can be changed. Having a single pressure function using a constant gamma is not going to work for the RadWS method.

giovanni-rosotti commented 7 years ago

The main concern right now is to why the EOS.h Pressure function is not inheritable, and whether this can be changed.

I did this change some months ago. I don't think there was any fundamental reason why Pressure should be defined in the base class; it's just that I noticed that all subclasses had exactly the same implementation, and therefore I just got rid of useless code. You are welcome to make it virtual (so that you don't have to reimplement it in every base class), and override it in RadWS. I suppose that the effect on the computational time should be quite small, but it would be great if you could time the code before and after this change, and let us how much it has changed. If it's a few % we can live with that.

RoguePotato commented 7 years ago

In changing the EOS::Pressure function to be virtual and to take a Particle address, I noticed an issue whereby two calls were passing in HydroForcesParticle (see the GradhSph.cpp changes in latest commit), but HydroForcesParticle doesn't inherit from Particle. I couldn't get around this issue so I kept the original template function but renamed it to PressureHydroForces. If there is a solution such that only one EOS::Pressure function is needed, that would be great.

RoguePotato commented 7 years ago

Is there a place in the SPH code specifically where gpot from sinks is added to the particles? I feel I have not set gpot_hydro to gpot in the correct place (see GradhSph.cpp diff in the first commit). I may think of another way of doing this as to remove the gpot_hydro particle data field.

dhubber commented 7 years ago

Is there a place in the SPH code specifically where gpot from sinks is added to the particles? I feel I have not set gpot_hydro to gpot in the correct place (see GradhSph.cpp diff in the first commit). I may think of another way of doing this as to remove the gpot_hydro particle data field.

I think it's added in GradhSph::ComputeStarGravForces (and similar routines for the meshless I presume). This is called from the GradhSphTree::ComputeAllSphForces routine, when all other forces are also computed.

RoguePotato commented 7 years ago

Okay so I found the issue: the gpot_hydro was being added to, but not being reset, so it kept increasing. Should be fixed now and is used by RadWS.

I've also been working on a separate radiative feedback branch (watch this space!) and have noticed some possible issues with this RadWS implementation. I'll have to do some more tests to see.

rbooth200 commented 7 years ago

Now included in radiative feed back