High temperatures in funnel zone of blondin simulations

[Higginbottom]

This behaviour can be reproduced by running either of the example simulations in the PY_PLUTO/blondin subdirectory of the Test_Problems directory.

Essentially what happens is that the cells close to the rotational axis heat up to very high temperatures. This causes the simulation to be rather slow, since the sound speed is high and the cells are small meaning that the dt to cross a cell is very short.

Here is a little figure showing the high temp...

[Higginbottom]

Tracing the energy changes in cell 2,2 - the innermost 'real' cell in both directions, I see that the main energy input when the cell is rapidly heating appears to be the term

rhs[j][ENG] -= phi_c*rhs[j][RHO];

in rhs_source.c - around line 340.

This is where the energy term in the right hand side of the equations is made, and this term in particular is just the density x the gravitational potential.

phi_c = BodyForcePotential(x1[i], x2[j], x3[k])

where I,j,k is the cell in question. This function is defined in init.c

Equation 6.3 in the Pluto user guide is the one to look at here.

[Higginbottom]

Well, this is odd.
When the cells are heating up - what appears to be the density in the equation

rhs[j][ENG] -= phi_c*rhs[j][RHO];

is negative.

so this makes a positive energy. Need to find out what rhs[j][RHO] actually is. I bet it isn't density - more like density times velocity or something.....

[jhmatthews]

Hey Nick - I had forgotten this, but I think rho*phi is right for the energy equation the way it is written in pluto. It's because the energy equation is really written as total energy to be fully conservative, i.e we want to calculate the time derivative of internal energy + kinetic energy + GPE in each cell. See pg 59 here: http://plutocode.ph.unito.it/userguide.pdf

but that's interesting...is there already a density floor in the simulation?

[jhmatthews]

Actually, I'm not entirely sure that's right, at least according to the docs. The rho phi should be in the advection term but nowhere else - but, it doesn't appear in the "flux" bit of the code, so maybe they rewrite it to appear on the RHS.

[Higginbottom]

Interesting - there is a density floor - but I'm just going to check that something isn't going wrong here. In particular I want to work out what the term rhs[j][RHO] actually is - if it is supposed to be density, then there is a problem...

[jhmatthews]

Ah ok, it's set in rhs.c. It looks right.

 rhs[j][RHO] = -dtdV*(fA[j][RHO] - fA[j-1][RHO]);

So that's essentially the 1D form of a div(rho * v) i.e. the RHS of the mass conservation equation, multiplying it by phi then gives me my div (rho * phi * v) which is what I want.

[jhmatthews]

I've had numerical problems before where the code was kinetically dominated, meaning that when you subtract the kinetic energy from the total energy you could get negative/incorrect thermal energies. I could imagine similar things happening here (is rho * phi greater than E_t or something because the gas is cooling quickly?) but I don't have a specific feel for which way it would go.

[Higginbottom]

Note:
I've found that the term rhs[j][RHO] - mass conservation - is only negative in the J-direction sweep. This is going in theta.

[Higginbottom]

For the cell I'm looking at in detail - that ids cell 2,2 - at the time when the temperature is very high, the split between internal energy and kinetic energy is

KE=1.31e-14
U=8.94e-14

for a total of 1.03e-13

[jhmatthews]

How does that compare to rho*phi?