B. Just iterate over neighbours - order NM
Entries up to 04 Oct 2004 are reformatted versions of
~/CLOUDY/Results-2004/README.sonic
Attempts to get the transsonic models working.
ZS016.in starts off on the R-branch but has an explicit antishock set at a certain depth. I added the code for this to dynamics.c and parse_set.c as mirror images of the "set dynamic shock XXXX" code. The antishock is unphysical but the idea is that it should be a small jump. How it works is that the pressure solver is initially looking for a supersonic solution but after the antishock it switches over to looking for a subsonic solution.
I tried to choose the antishock depth to be the minimum in the velocity from a weak-R test. This works quite well for the first few iterations: the antishock occurs at the H i-front and is a jump of < 10% in the velocity. However, eventually the solution ends up hitting the sonic point before getting to the antishock depth - and there it stops with a pressure failure.
The ZX series have the problem that the antishock is too strong in the initial iterations (because I chose its position according to where it should be in the later iterations). As a result, the subsequent iterations develop an echo of the antishock that is one advection length downstream. Late, we get an echo of the echo, etc. As a result, it never converges. It might be better to have the antishock at a given Mach number, rather than a given radius. That way, it would happily do a weak-R solution for the first few iterations.
In theory, this should work with the "strongd" pressure model. In that case, it shouldn't be necessary to specify a depth at which the solver switches from the supersonic to subsonic branch - it should do it automatically when it fails to find a supersonic solution. At the moment this is not working beacause we get a pressure failure at the sonic point.
ConvPresTempEdenIoniz() repeatedly calls PressureChange(), which calls lgConvPres() (in the same file). It exits the loop once conv.lgConvPres is TRUE (or it aborts or exceeds max iterations).
For the wind case, lgConvPres() calls DynaPresChngFactor() to calculate pressure_change_factor. This has branches for the subsonic and supersonic cases and uses the history of previous pressure evaluations from previous iterations.
Then lgConvPres() sets its return value and conv.lgConvPres to TRUE if pressure_change_factor is in range 1.0 +/- conv.PressureErrorAllowed. Otherwise, FALSE.
PressureChange() also sets conv.lgConvPres to the return value of lgConvPres(), which seems redundant...
Then PressureChange() returns FALSE if conv.lgConvPres is TRUE. Otherwise, it carries on to calculate some extra stuff (mainly to do with abundance fluctuations). The value of pressure_change_factor gets capped to some value, depending on conditions. Then it returns the value of lgChange, which looks like it will be TRUE in general.
With the strongd or the antishock-by-mach pressure modes, then it is possible to hit a region where no pressure solution is possible for a while. There is some logic in DynaPresChngFactor() to detect this. However, it sets conv.lgConvPres to FALSE so that ConvPresTempEdenIoniz() sees it as a simple pressure failure. If this happens repeatedly, then it gives up and calls ConvFail("pres"), which see that the gas and ram pressures are nearly equal and so sets pressure.lgSonicPoint. When this flag is seen by lgEndFun(), it will stop the program. We need some way to indicate to ConvPresTempEdenIoniz() that we are in such a region and that it should just get over it and move on to the next zone.
Question is, do we want it to stop looking for a pressure solution immediately when this happens? Assuming yes, then I have added a new variable, pressure.lgStrongDLimbo, that gets set when this happens.
Actually, I've found a better way - we now have a new flag, pressure.lgSonicPointAbortOK, which is set to FALSE by the pressure modes that might want to go through the sonic point, thereby stopping lgEndFun() from ending it all. There may still be a use for lgStrongDLimbo, so I'll leave it in for now.
Things sort of work for a few iterations if I choose the right value of the various convergence tolerances, but it always falls over adter a bit. One problem is that the model I am doing is starting off very close to sonic.
In most of the models, there are 3 maxima in c-squared:
The third one is slightly higher than the other two, but not by much. This makes things difficult since any model that is close to making a sonic transition at the H i-front is also close to sonic at the illuminated face. This often causes pressure failures in the first few cells. Hence, it would be better to have a model where there was only one c-squared maximum. This is difficult to find. I have tried reducing the ionization parameter but this doesn't make much difference.
Many weak-R models fail in similar ways. This seems to be due to the initial velocity being too close to sonic (but "close" can mean as much as a factor of 2). The symptoms are that the density in the first zone gets set to a much lower value than was specified in the init file. Then we get pressure failures a few zones further on. What can be causing this?
The target pressure that we should be aiming at in each zone seems to be pressure.PresTotlInit. This is set in various places, all with the same form:
PresTotCurrent();
...
pressure.PresTotlInit = pressure.<nop>PresTotlCurr;
pressurechange.c:84 iter_startend.c:200, in IterStart() conv_init_solution.c, various places in ConvInitSolution()
PresTotCurrent() simply evaluates the current pressure from all constituents.
Something must be going wrong somewhere, since it ought to be guaranteed that the density at the illuminated face is equal to the hden specified in the init file.
XA020 is an example of this failure mode.
Learnt new things about gdb: automatic display list. Discovered that
problem lies in ConvInitSolution(), and is not at any of the places I
had thought: there are other places where pressure.PresTotlInit is
changed by the '
Even before the target pressure has changed at all, ConvPresTempEdenIoniz is failing to keep the density steady.
Sent mail yesterday to Gary and Robin, querying the following code:
/* Subsonic case: pressure ^ with density ^ => increase density further */
/* Super " case: pressure ^ with density v => decrease density further */
/* printf("Presmode %d flag %g factor %g\n",zonepresmode,(er-erp)*(dense.gas_phase[ipHYDROGEN]-dp),factor); */
if (zonepresmode == SUBSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) < 0)
{
factor = 1+3*width;
}
else if (zonepresmode == SUPERSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) > 0)
{
factor = 1-3*width;
}
I originally thought that the code was in error, since it is doing the opposite of what it advertises in the comment. However, I now think that it is the comment that is wrong, and this code is actually trying to force the solution to the requested branch (supersonic/subsonic) in cases where dP/drho is opposite to what is expected for that branch.
The problem arises when it tries to do this on the first zone, since it ends up forcing down the density to insanely low values.
In principle, it should work:
If we want to be supersonic, then we ought to have dP/drho < 0 at constant mass flux. Hence, it is reasonable to suppose that if dP/drho > 0 then we are really subsonic and we need to decrease the density to push us into the supersonic regime. Could the problem be that the velocity is treated differently in the first zone and is not increasing as the density goes down?
Assuming plane-parallel (dynamics.FluxIndex == 0), then dynamics.FluxScale = wind.windv0, which is set from the initialization commands. The function DynaFlux() returns the value dynamics.FluxScale*struc.DenMass[0], and this seems to be the root of our problems: if the pressure solver pushes down the density in the first zone, then this will decrease DynaFlux(), which is used to set the velocity in PressTotCurrent() as
wind.windv = DynaFlux(radius.depth)/(phycon.xMassDensity);
Actually, this is only done for nzone > 2, so wind.windv in the first zone always stays as wind.windv0, even if the density is going down. Not good. What is needed is to recalculate wind.windv inside the pressure solver, even if it is the first zone. For this to work, we also need that DynaFlux() does not use struc.DenMass[0] but rather the initial value of the density, as specified in the init file.
In PressTotCurrent() we set a new variable
phycon.xMassDensity0 = phycon.xMassDensity
but only if this has not been set before (do this by setting to -1.0 in zero.c)
Then we remove the test (nzone > 1) for resetting wind.windv and change DynaFlux() to use phycon.xMassDensity0 instead of struc.DenMass[0]
This seems to work!
These are moderate ionization parameter models with phi(h)=11, hden=4. There is no well-defined He i-front so there is only one maximum in c-squared (other than the one at the illuminated face).
XA013.in has "set dynamics pressure mode supersonic" and does the typical weak-R thing. On the first iteration, the isothermal Mach number start off around 1.14, increases to 1.4 and then dips to about 1.24 at the i-front, before shooting up on the neutral side (I have "stop efrac -1.5", so it doesn't go very far). Subsequent iterations pull the i-front inwards a bit and the Mach number dips down further as the T-peak at the i-front gets sharper. By the 5th iteration or so, the Mach number hits unity, but since we have the supersonic pressure mode it bounces back again.
XA013A.in has an antishock set at M=1.05. This succeeds at forcing a transsonic solution after a few iterations. However, on the next iteration it switches back to a weak-R solution because the advected terms have changed enormously within one advection length downstream of where the sonic point was, depressing the sound speed there so that the Mach number no longer dips down enough. On further iterations, the solution recovers and goes sonic once again but the whole cycle repeats itself every 4 iterations.
XA013A2.in is similar but with the antishock set at M=1.01 and the convergence criteria tightened up. It behaves similarly.
Next thing to try - it would be nice to reduce the initial velocity so that even the first iteration hits the sonic point in the i-front. This is not easy because c-squared is higher at the illuminated face than it is at the i-front. Perhaps it would work to run a static model through to the minimum in c-squared halfway through the ionized region. Then start a dynamic model with the transmitted continuum from that.
This doesn't work quite how I had hoped: the T in the model using the transmitted continuum (XB013) doesn't match with that from that of the constant pressure model (XA000P). Is this something to do with the diffuse field? This doesn't even work when I try and match one constant-gas-pressure model (XB000P) to another constant-gas-pressure model (XA000P).
Heating is basically H 1 and He 1. Cooling is a mixture of forbidden lines and 'Hp 0 0.0', whatever that is. The heating rate falls off with depth, which is presumably what drives the temperature change. Heating is significantly higher at the illuminated face of XB000P than it was at the back face XA000P, even though these are in theory the same place. The neutral fraction is also somewhat higher there, which would explain this. Maybe this is due to escaping lyman lines?
YES!! I have found the answer. I still can't get the temperatures to match up but it doesn't matter: if I use the "case b no photoionization" option (which stops Ly-a from escaping from the illuminated face and turns off photionization from excited states), then the temperature now increases monotonically towards the i-front. So, I don't need to mess around with the transmitted continuum any more.
On a side note, I don't understand why I don't get an equivalent result from just using the "sphere" command: in a closed geometry, each ly-a photon that escapes from the illuminated face should be exactly balanced by one coming in from the other side of the nebula, shouldn't it? UPDATE: it seems "sphere static" is what I wanted.
Added option "set dynamics pressure solver quadratic" to turn on the quadratic pressure solver via a flag dynamics.lgQuadPresSolver. Changes to dynamics.h, parse_set.c, dynamics.c (DynaZero() and DynaPresChngFactor())
Model XC014A2 goes transsonic on the first iteration. It is trying to go through the sonic point before the c-squared max has been reached, so it has a sequence of pressure failures before getting out the other side. The changes I have made to the zoning logic manage to stop radius.drad from getting too small during this, so we do manage to get out the other side. The quantity degree of pressure failure can be quantified as (P-P0)/P0, where P is the best pressure found (Pcurrent) and P0 is the target pressure (Pcorrect). This behaves quite nicely over the pressure-failure region, rising up to about 1%, then falling back as we come out the other side.
On the next iteration, we don't get a sonic point at all because all the mess around the sonic point has been advected back towards the star and causes a bunch of pressure failures before we even get near the sonic point. The c-squared maximum is all messed up and the solution follows the weak-r path.
What to do?
I had originally thought of trying to converge to Pgas=Pram if we fail to converge to Ptotl=Pcorrect. However, looking at the quadratic solver, it seems that it will look for a minimum in the P(rho) curve if it can't find a root, which is exactly what we need (clever Robin!)
What happened?
Using the quadratic solver, we do indeed get a much better transition through the limbo around the sonic point: the velocity and density now vary more-or-less smoothly. However, we still get the problem that the next iteration goes weak-r on us.
wind velocity -43 index -2 center 10000000000000000 no continuum
Strangely, these seem to work. But I'm not sure mass flux is working correctly. They all converge on the 4th iteration. Not sure that advection is being turned on properly. ah, should be "wind advection ..."
I have a model that is designed to look like a Helix knot - XDglob6.in
The first iteration works, just about. I start out at a radius of 1.e17 cm with a density of 0.1 pcc, a 120,000 K black-body flux of Phi(H)=10**8.5, and a velocity of -64.35 km/s. The ionization front then occurs at a radius of about 1.e15 cm.
There is no csq-maximum but it goes through the sonic point anyway, which it is allowed to in spherical geometry. The sonic point is also at a very low value of H+ fraction - about 1%, although He+ is around 10%. Ne peaks a bit nearer the star, around x(H+) = 0.1, with a value of around 1.e3, but there is also an Ne peak on at the far-neutral boundary, which is worrying. This is where the hydrogen density shoots up to about 1.e6. This is shown in XDglob6-zoom.eps, which shows various physical variables in the region around the i-front.
I also did the optical line structure: XDglob6-lines.eps, although this is a bit meaningless since it is only the first iteration without advection. [O I] is very strong because of the low ionization fraction around the Ne peak. Ha emission is similar to [N II] on the ionized side, albeit 4 times weaker, and extends further in on the neutral side, with a silly sharp peak at the inner edge. [O III] is extremely weak, with a broad peak about 4 times farther out. He I is very concentrated on the neutral side and is still climbing at he inner edge
Gnuplot commands:
plot 'XDglob6.ovr' i 0 u (1.e17-$1):(10**($2-4)) t 'T/1.e4', \ '' i 0 u (1.e17-$1):(10**($5-3)) t 'ne/1.e3', \ '' i 0 u (1.e17-$1):(10**$8) t 'x(H+)', \ '' i 0 u (1.e17-$1):(10**$10) t 'x(He+)', \ '' i 0 u (1.e17-$1):(10**($5-$4)) t 'xe', \ '' i 0 u (1.e17-$1):(10**($4-3)) t 'nh/1.e3', \ 'XDglob6.pre' i 0 u (1.e17-$1):(-$10/30.) t 'v/30', \ '' i 0 u (1.e17-$1):(sqrt(3./5.)*$11/30.) t 'c/30', \ '' i 0 u (1.e17-$1):(sqrt(3./5.)*sqrt($12/$6)*$11/30.) t 'alf/30' plot 'XDglob6.lines' i 0 u (1.e17-$1):(5*$2) t 'ha x 5', \ '' i 0 u (1.e17-$1):(8*$4) t 'he i x 8', \ '' i 0 u (1.e17-$1):(2*$11) t 'n ii x 2', \ '' i 0 u (1.e17-$1):14 t 'o i', \ '' i 0 u (1.e17-$1):(50*$9) t 'o iii x 50'
The second iteration is a complete nonsense. The advection pushes down the temperature while we are still way out in the flow. T drops to zero before we get in to a radius of 5.e16, so the calculation stops. The third iteration is even worse and hardly gets started before T drops to zero. Then the fourth iteration, which has no real previous iteration to get the advection from, comes out like the first iteration again, and so it continues....
Proposed changes to the way advection terms are added (DynaNewStep, DynaSaveLast, DynaIonize)
This could be done in DynaSaveLast, or would it make more sense to have a new function, which was called at the start of a new step?
Two options:
A. Convolution by FFT - order
B. Just iterate over neighbours - order NM
where N is number of zones and M is number of neighbours (i.e., those zones within one smoothing length of each zone).
Option B would be a lot easier to code and will not be too inefficient if we do box smoothing, since smoothing length is a fraction of advection length, which in turn is a small fraction of total depth.
An idea might be to restrict the models to one stopping criterion. As it is, where we have the possibility of stopping on either temperature or electron density then which one is operating is changing from one iteration to the next, which may be contributing to the flip-flops.
With Torsten, I am looking at the globule models again with a view to do