Generic capability to read in wind models

[kslong]

Currently, it is possible to read in a cylindrical model to Python, but we cannot read in a simple spherical model. It would be very useful to do this, if for example we wanted to make a more direct comparison with a Cloudy model or for calculating the spectrum that we would get from a pre-calculated stellar atmosphere, such as provided by Tlusty.

[Higginbottom]

The machinery to read in a hydro file is all about spherical coordinates, so it should be possible to use this…

On 30 Oct 2017, at 21:25, Knox S. Long notifications@github.com wrote:

Currently, it is possible to read in a cylindrical model to Python, but we cannot read in a simple spherical model. It would be very useful to do this, if for example we wanted to make a more direct comparison with a Cloudy model or for calculating the spectrum that we would get from a pre-calculated stellar atmosphere, such as provided by Tlusty.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/agnwinds/python/issues/309, or mute the thread https://github.com/notifications/unsubscribe-auth/ADLMvSVolGi_GQU9mJ_PtpAu4ENYIVICks5sxj7YgaJpZM4QL4K0.

[kslong]

I should have been clearer . Cloudy and stellar models are 1D. The hydro models are 2D spherical polar. We cannot reard in 1D.

[kslong]

The standard fort.7 file that Tlusty produces looks like this:

   50 1130
 3.208000D-07 5.065000D-07 8.008000D-07 1.266000D-06 2.003000D-06 3.165000D-06
 5.000000D-06 7.896000D-06 1.246000D-05 1.967000D-05 3.102000D-05 4.892000D-05
 7.711000D-05 1.215000D-04 1.912000D-04 3.008000D-04 4.725000D-04 7.413000D-04
 1.161000D-03 1.816000D-03 2.837000D-03 4.424000D-03 6.887000D-03 1.070000D-02
 1.657000D-02 2.554000D-02 3.912000D-02 5.948000D-02 8.964000D-02 1.339000D-01
 1.984000D-01 2.919000D-01 4.278000D-01 6.263000D-01 9.178000D-01 1.349000D+00
 1.988000D+00 2.936000D+00 4.329000D+00 6.341000D+00 9.195000D+00 1.351000D+01
 2.092000D+01 3.357000D+01 5.500000D+01 9.069000D+01 1.496000D+02 2.444000D+02
 3.946000D+02 3.959000D+02
   1.475795D+04   7.015288D+07   1.523895D-16   6.435971D-02   1.491812D-05
   9.424856D-06   1.271046D-05   1.875043D-05   2.648736D-05   3.550609D-05
   4.599836D-05   1.771616D-01   6.361591D+07   7.733330D+00   3.179224D-06
   1.966796D-07   3.448150D-06   1.395219D-07   2.050797D-07   3.091883D-08
   6.854136D-07   5.628826D-07   6.350994D-08   5.109919D-08   1.156610D-07
   2.133967D-08   3.965577D-07   4.448020D-07   7.497742D-08   4.099334D-07

This was just part of the file but that is what the first part looks like, and the remainder is similar.

The format is described in Volume 3, of the new Tlusty guides, in section 8.1.

The first line has the number of depth points and number of parameters. Then comes the depth grid which is in grams/cm**2. Then you get the information for all of the depth grids, these are wrapped. The first 3 parameters are T, ne, and rho (in cgs coordinates). There is no explicitly distance for stars. Disks are different; there you do get a distance from the disk plane. The remainder of the values involve level populations, which probably do not need to be importated into Python directly.

So to create something reasonable for Python one would either need to write a routine to interpret this within Python, or more reasonable, write a translation routine that would produce parameters more like what Python would use.

[kslong]

The Kurucz model atomospheres, seem to look like the following:

TEFF   5000.  GRAVITY 2.00000 LTE 
TITLE  [+0.5] VTURB=2  L/H=1.25 NOVER NEW ODF                                   
 OPACITY IFOP 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 0
 CONVECTION ON   1.25 TURBULENCE OFF  0.00  0.00  0.00  0.00
ABUNDANCE SCALE   3.16228 ABUNDANCE CHANGE 1 0.91790 2 0.07810
 ABUNDANCE CHANGE  3 -10.94  4 -10.64  5  -9.49  6  -3.52  7  -4.12  8  -3.21
 ABUNDANCE CHANGE  9  -7.48 10  -3.96 11  -5.71 12  -4.46 13  -5.57 14  -4.49
 ABUNDANCE CHANGE 15  -6.59 16  -4.71 17  -6.54 18  -5.64 19  -6.92 20  -5.68
 ABUNDANCE CHANGE 21  -8.87 22  -7.02 23  -8.04 24  -6.37 25  -6.65 26  -4.54
 ABUNDANCE CHANGE 27  -7.12 28  -5.79 29  -7.83 30  -7.44 31  -9.16 32  -8.63
 ABUNDANCE CHANGE 33  -9.67 34  -8.63 35  -9.41 36  -8.73 37  -9.44 38  -9.07
 ABUNDANCE CHANGE 39  -9.80 40  -9.44 41 -10.62 42 -10.12 43 -20.00 44 -10.20
 ABUNDANCE CHANGE 45 -10.92 46 -10.35 47 -11.10 48 -10.27 49 -10.38 50 -10.04
 ABUNDANCE CHANGE 51 -11.04 52  -9.80 53 -10.53 54  -9.87 55 -10.91 56  -9.91
 ABUNDANCE CHANGE 57 -10.87 58 -10.46 59 -11.33 60 -10.54 61 -20.00 62 -11.03
 ABUNDANCE CHANGE 63 -11.53 64 -10.92 65 -11.69 66 -10.90 67 -11.78 68 -11.11
 ABUNDANCE CHANGE 69 -12.04 70 -10.96 71 -11.98 72 -11.16 73 -12.17 74 -10.93
 ABUNDANCE CHANGE 75 -11.76 76 -10.59 77 -10.69 78 -10.24 79 -11.03 80 -10.91
 ABUNDANCE CHANGE 81 -11.14 82 -10.09 83 -11.33 84 -20.00 85 -20.00 86 -20.00
 ABUNDANCE CHANGE 87 -20.00 88 -20.00 89 -20.00 90 -11.95 91 -20.00 92 -12.54
 ABUNDANCE CHANGE 93 -20.00 94 -20.00 95 -20.00 96 -20.00 97 -20.00 98 -20.00
 ABUNDANCE CHANGE 99 -20.00
READ DECK6 72 RHOX,T,P,XNE,ABROSS,ACCRAD,VTURB, FLXCNV,VCONV,VELSND
 2.91268412E-03   2660.1 2.912E-01 3.057E+07 4.578E-05 2.093E-02 2.000E+05 0.000E+00 0.000E+00 1.225E+06
 3.83035748E-03   2706.2 3.830E-01 4.377E+07 5.115E-05 2.007E-02 2.000E+05 0.000E+00 0.000E+00 1.110E+06
 4.93448092E-03   2748.1 4.934E-01 6.078E+07 5.628E-05 1.927E-02 2.000E+05 0.000E+00 0.000E+00 1.018E+06
 6.25256704E-03   2791.3 6.252E-01 8.375E+07 6.372E-05 1.861E-02 2.000E+05 0.000E+00 0.000E+00 9.454E+05
 7.80746936E-03   2833.2 7.806E-01 1.134E+08 7.179E-05 1.800E-02 2.000E+05 0.000E+00 0.000E+00 8.864E+05
 9.65812494E-03   2874.1 9.657E-01 1.516E+08 8.019E-05 1.745E-02 2.000E+05 0.000E+00 0.000E+00 8.378E+05
 1.18372528E-02   2915.3 1.184E+00 2.010E+08 9.193E-05 1.700E-02 2.000E+05 0.000E+00 0.000E+00 7.980E+05
 1.43731207E-02   2955.0 1.437E+00 2.624E+08 1.052E-04 1.658E-02 2.000E+05 0.000E+00 0.000E+00 7.654E+05
 1.73426699E-02   2993.4 1.734E+00 3.383E+08 1.192E-04 1.621E-02 2.000E+05 0.000E+00 0.000E+00 7.383E+05
 2.08441024E-02   3030.9 2.084E+00 4.322E+08 1.348E-04 1.586E-02 2.000E+05 0.000E+00 0.000E+00 7.155E+05
 2.49554221E-02   3068.1 2.495E+00 5.480E+08 1.537E-04 1.555E-02 2.000E+05 0.000E+00 0.000E+00 6.967E+05
 2.97709316E-02   3104.3 2.977E+00 6.888E+08 1.746E-04 1.534E-02 2.000E+05 0.000E+00 0.000E+00 6.810E+05

The format looks fairly self explanatory, but there is a descriptin of it here

[kslong]

I've changed the title of this issue because we need to deal more generically with reading in wind models (given the work we are starting with Lee Hartman & Kelley Milliner). They are going to need to read in models in cylindrical coordinates; I have a need to be able to read in models in spherical (1d) coordinates.

We currently have the capability to read in the zeus model outputs. These models are in polar (r-theta) coordinates. The velocities are given as v_r, v_theta, and v_phi. Care is taken to make sure that the velocities and positions are at the edges of cells, and rho (and T) are at the centers of cells. The input format is tailored a bit to the zeus interaction because the files have some header type information. Correct? @Higginbottom

We do not have the ability to read in models on a cylindrical grid or on a simple 1d spherical grid. We need to create these, and in the process decide whether we need a more generic way to read in polar coordinates.

Various questions arise which need comments by everyone interested.

1. What should the format of the input files be, and should it have header information that effective tells you the coordinate system. The alternative is to select a coordinate system and then read in a file. Currently, we ask for the wind type before the coordinate system, so in principle, the file that contains the read-in wind could also contain the information about what coordinate system it is. One could image a file that begins say with

   # CoordSystem POLAR

   which would describe the coordinate system.

2. Should we assume that all of the variables that are read in are at the same position, and that we have to do some kind of interpolation to get the ones that are say cell_centered, or should we assume the user understands that she has to give us rho at the center of the cell but velocities at the edges, or should we provide flexibility in this regard, e.g add another header line allows us to deal with this internally, e.g a choice?
3. In the current hydro scheme, we read positions, i,j, velocities, rho, and a temperature. If we think this is all we ever want to transmit, then we can use sscanf and let the number of valid entries tell us that we have a temperature, but if there were other optional things to transmit we cannot do this. (Life would be easier if c programs could read astopy tables).
4. I hate to come back to this question because @smangham will likely come back with some modern answer that we are Luddites unless we use some kind of standard file transfer that would allow us to read keywords and tables. Still, I think this is a good time to revisit this issue. Note: I assume we are going to get the original wind models in some format that is not immediately readable into Python, and that rather than incorporating all of these into Python directly, we should be expecting to write a little translator routine in the programming language named for our radiative transfer program that outputs a file Python will read.

I've probably not gotten all of the questions down here, and we may not want to go all the way to a final solution to meet our immediate needs, but I definitely think we should decide where we would like to go now and why.

[Higginbottom]

On 2 Nov 2017, at 20:55, Knox S. Long notifications@github.com wrote:

I've changed the title of this issue because we need to deal more generically with reading in wind models (given the work we are starting with Lee Hartman & Kelley Milliner). They are going to need to read in models in cylindrical coordinates; I have a need to be able to read in models in spherical (1d) coordinates.

We currently have the capability to read in the zeus model outputs. These models are in polar (r-theta) coordinates. The velocities are given as v_r, v_theta, and v_phi. Care is taken to make sure that the velocities and positions are at the edges of cells, and rho (and T) are at the centers of cells. The input format is tailored a bit to the zeus interaction because the files have some header type information. Correct? @Higginbottom https://github.com/higginbottom This is all correct. Python does the interpolating to ensure that the relevant quantities are worked out at the centre or edge of cell, but it is all based upon knowing where the zeus quantities are from. Clearly if we wanted to make it general, there would have to be some rule, or keyword...

We do not have the ability to read in models on a cylindrical grid or on a simple 1d spherical grid. We need to create these, and in the process decide whether we need a more generic way to read in polar coordinates.

Various questions arise which need comments by everyone interested.

What should the format of the input files be, and should it have header information that effective tells you the coordinate system. The alternative is to select a coordinate system and then read in a file. Currently, we ask for the wind type before the coordinate system, so in principle, the file that contains the read-in wind could also contain the information about what coordinate system it is. One could image a file that begins say with CoordSystem POLAR

which would describe the coordinate system.

2. Should we assume that all of the variables that are read in are at the same position, and that we have to do some kind of interpolation to get the ones that are say cell_centered, or should we assume the user understands that she has to give us rho at the center of the cell but velocities at the edges, or should we provide flexibility in this regard, e.g add another header line allows us to deal with this internally, e.g a choice?

I guess that as long as you give the coordinates at which each parameter is defined, one can interpolate….

3. In the current hydro scheme, we read positions, i,j, velocities, rho, and a temperature. If we think this is all we ever want to transmit, then we can use sscanf and let the number of valid entries tell us that we have a temperature, but if there were other optional things to transmit we cannot do this. (Life would be easier if c programs could read astopy tables).
4. I hate to come back to this question because @smangham https://github.com/smangham will likely come back with some modern answer that we are Luddites unless we use some kind of standard file transfer that would allow us to read keywords and tables. Still, I think this is a good time to revisit this issue. Note: I assume we are going to get the original wind models in some format that is not immediately readable into Python, and that rather than incorporating all of these into Python directly, we should be expecting to write a little translator routine in the programming language named for our radiative transfer program that outputs a file Python will read.

I've probably not gotten all of the questions down here, and we may not want to go all the way to a final solution to meet our immediate needs, but I definitely think we should decide where we would like to go now and why.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/agnwinds/python/issues/309#issuecomment-341555147, or mute the thread https://github.com/notifications/unsubscribe-auth/ADLMvSzM8IV2QQUuDwfnXhNDphfuWu5cks5syiwugaJpZM4QL4K0.

[kslong]

As of now, the code sort of works for a spherical model. Doing all of this is tricky, and is going to be even more tricky for the other coordinate systems.

One of my more general thoughts about this, is that in general the file we want to read in, if given to us by others will often require some little python script (as Nick does for the zeus interface). This little script should be used to standardize the inputs. This could simplify what Python needs to do.

As an aside: I believe the way the test this code is to use our various models, e.g SV, and take the windsave2table files, "munge" them slightly, and read them back in for read in mode.

One more point. The hydro model has not been integrated completely into domains, and it should be. But it is also tricky, or at least the most general case is tricky, because what one does is read in the model first, and then sometime later one interpolates the model into the wind and plasma arrays. But we store the read in models into a set of fixed arrays, and if one tried to do this twice, one would have to be quite careful, not to write one over the other. This should seldom be a practical problem, but it is the sort of thing that if we could solve it easily would simplify the code. Note that this could be as simple as adding a counter that would say one has already read in, say a cylindrical model, and cannot do it twice.

[kslong]

OK, the current version on the import branch has the beginning of an attempt to put cylindrical models in. It still does not work, but the philosophy if @jhmatthews wishes to take up the cudgel, is to read as much as one can into wmain at the start, e.g the positions and the velocities. If one does this, then one can use interpolation to calculate the velocities at the centers of the cells, which one needs for the gradients. The biggest issue with all of this is the padding cells. An easy way out would be to require the padding cells in the input grid, and if it is absent from the data we receive from, say, Lee, to add the padding before reading in the data. The issue about what to do about the wind cones remains.

Note that I have added to little routines to py_progs that are intended to be used for creating input grids for the python import model option. One runs a normal model, say for a cv, and then windsave2table on the completed model. This produces various text files. The routines take these text files and modify them to the specific format that python import expects.

[kslong]

I finally have a version of the cylindrical import model that does something other than fall over immediately. Along the way, I split out the routines for the various coordinate system import.

For the cylindrical model, I have adopted the philosophy that the input model should have all of the buffer cells included. I do not try to construct them. I assume also that the quantities that should be cell centered, e.g rho, are cell-centered, and those that are supposed to be located on the edges of cells are located on the edges of cells.

There are still issues as of today, with knowing which cells are in the wind or not. I will assume that this is in the input file but I have not fixed this yet. Note that all this means a good deal of the burden is going to be on whatever script is used to create the input file. In particular, there are some issues about how to extrapolate cell-centered variables, that have not been worked out.

[jhmatthews]

I've got a CV model fundamentally working, in that I can run Knox's python script on a windsave2table output from cv_standard.pf, and it produces a reasonable looking spectrum:

The differences are probably due to the fact that, at the moment, the region outside the wind cone appears to be populated with densities and velocities, which needs fixing.

[jhmatthews]

One of the to do items here is to deal with the inwind variable. I've first added something that copies the inwind variable from the file, and checked if it gets overwritten by the other machinery. It currently gets overwritten in two places:

• the "not in wind" cells get modified by "cylind_volumes()" - which is fixable by adding an additional flag such as "INWIND_NOT_ASSIGNED". Fairly trivial.
• the cells get modified by check_corners_in_wind(), which incorrectly overwrites some cells.

The issue that pops up is related to the use of where_in_wind. where_in_wind takes a position vector, x, and works out if that position is all in the wind or not. To do so, it uses the minimum and maximum angles associated with the wind cone. For an imported model, you don't necessarily have an angle - you have a series of cells with densities and so on, and the angles are set to zero. Fundamentally, one has already decided which cells are in the wind and should be allowed to have arbitrary geometries (in theory).

The reason this breaks the code I've written is because it finds that cells without a volume have four corners in wind, because angles are not properly set. But, more generally it will clearly fail for photon propagation, because we won't be able to establish if the photon is in the wind or not.

I see two options:

• replace the current where_in_wind method with one that loops over cells to find if we are in the wind or not, based on the inwind flag in the wind pointer. The problem with this is the "partially in wind" cells, which are dealt with ok with the current method.
• have a flag to choose which method to use (the above loop over cells or do the trigonometry as we do at the moment)

I favour the latter. It preserves current behaviour with existing kinematic models, is presumably faster in those models, and allows one to deal with "part in wind" cells.

What do you think folks? @kslong @Higginbottom ?

[kslong]

I'm not quite understanding the problem. Are we talking about in the initial setup of the model, or later on, e.g when photons are flying.

In the setup of the model, for imported models, I was assuming we would jump_around any steps that involved windcones, using the wind_type(==IMPORT), and that \in cylind_volumes we would assign the total volumens.

Partly, I'm confused by the idea that we would ever need to cycle through the cells to find the relevant. So far, all of our grids are regular, so we can always use the method that goes from x to n, and then we once the grid is set up, we could then avoid any checks of wind-cones etc. In that sense this just seems like an ordering of how we check things.

But as I say I am probably confused.

[jhmatthews]

Both, I suppose. The problem is that currently we have a function where_in_wind which is called in a number of places:

• to check how many corners each cell has in the wind (but as you say, this could be skipped)
• to decide if a photon is in the wind or not (which is always needed)

At the moment, this function works by taking a position, and comparing it to the geometry of the model -- it doesn't actually use the grid at all, it just compares the positions to the inner and outer edges of the flow with trig.

So the superficial issue is that the current import model breaks because it doesn't have angles associated with it and so where_in_wind doesn't work. But the more general issue is whether where_in_wind() should be able to operate without having angles defined, and I believe it should, and thus it needs to use the grid information and cycle over them instead of doing geometry. The second general issue is that that means we can't really deal with "partially in wind" cells, because we would have no way of defining where in the cell the wind stops.

[kslong]

I have looked at where_in_wind, now. As you say it uses wind_cones etc to decide if a positions is in the wind region or not, and also importantly which wind domain the object is in. It searches through the domains in reverse order, consistent with the philosophy that one might want to put a corona above a wind model.

I believe what is done for the zeus hydro model, and what I assumed we were doing for this, was to set the wind_cones (and rho_min, rho_max, etc) so that all but the buffer cells are in the wind. I had assumed we would do this for the imported models. (I can imagine in future that there would be models where we might actually want to read in information about the wind cones.)

I guess the question here is what do you do about the possibility that some of the cells in a typical import model are actually not in the wind. I could imagine therefore adding a check, that involves adding a line that asks where_in_grid one is to get the index into wmain for that position, and then since using the value of inwind to continue.

Regarding the imported models, I think we need to set all of the cells on import to either in the wind or not, and so this additional check would always give a clean result.

Now what I am referring to above is for when one makes a call to during the photon flight portions of the program. I do think one needs to judiciously step around some of the logic when getting things up, by using the fact that the model wind_type is IMPORT.

[jhmatthews]

Great Knox - that is exactly what I had started doing. I have made some progress on this with inwind and so on, but there's still something up with the volumes that I need to figure out. Will update tomorrow.

[jhmatthews]

My recent commit to import contains the following changes (https://github.com/agnwinds/python/commit/e6355efb84c4a0385f589c965e4baaa1a42bb242 w/ additional comments in https://github.com/agnwinds/python/commit/c3c443fbddbf5b6e594449627bbae41a61eb71da)

• added a step to copy the inwind value from the imported model to the wind pointer
• an additional flag called W_NOT_ASSIGNED which is initially given to wind cells. This way, the cylind_volumes machinery knows whether the inwind flag has already been set by an imported model or needs to be calculated.
• relevant if statements in cylind_volumes to deal with the above and not assign partial volumes in the case of an imported volume. Volumes are either set to all in wind volumes or zero.
• Currently, part in wind cells are set to not in wind. This may not be the best choice.
• added an if statement in where_in_wind. If we are using an imported model for that domain, then we loop over the cells and use where_in_grid to find the cell n, and the wmain[n].inwind value to find if we are in the wind. Current dev behaviour is preserved for non-imported models
• Skip "check_corners_in_wind" section for imported models

Current status: It sort of works - the wind printout now looks reasonable and the out of wind portions of the standard cv model are no longer filled with plasma. However:

• the spectra doesn't match the non-imported model - there's no CIV P-Cyg at all now (see below)
• The volumes and densities don't appear to match the standard model (non imported)
• there are a number of errors, see below:

Error summary:

Error summary: End of program, Thread 0 only
Recurrences --  Description
        7 -- getatomic_data: line input incomplete: %s
       94 -- vwind_xyz: Cannot determine an xyz velocity on z axis. Returnin 0,0,v[2]
        1 -- check_grid: velocity changes by >1,000 km/s in %i cells
        1 -- check_grid: some cells have large changes. Consider modifying zlog_scale or grid dims
        5 -- trans_phot: %ld photons were lost due to DFUDGE (=%8.4e) pushing them outside of the wind after scatter
        6 -- trans_phot: Trying to scatter a photon in a cell with no wind volume
        6 -- trans_phot: grid %3d x %8.2e %8.2e %8.2e
        6 -- trans_phot: This photon is effectively lost!
        6 -- walls: distance %g<0. Position %g %g %g 
        2 -- disk_init: Got to ftot %e at r %e < rmax %e. OK if freqs are high
      285 -- sane_check: %d %e
      285 -- trans_phot: sane_check photon %d has weight %e before extract
        2 -- error_count: This error will no longer be logged: %s
        1 -- pillbox %d interfaces to pillbox is impossible
Run py_error.py for full error report.
Thread 0 Finalized. All done
xxx Fri Nov 17 14:26:03 2017    248.1 COMPLETE             cv_import
Completed entire program.  The elapsed TIME was 248.114481

Spectrum comparison:

[kslong]

At present for the cv_test model, there as @jhmatthews notes, no partial cells. For the original cv model, there are 137 in the wind (same as test) but 102 partially in the wind. Consistent with this, many more photons just escape without scattering in the import version, 98000, vs 87000 in cycle.

At present, the import version gives a few hundred of vwind_xyz on the z axis errors in cycle 0, whereas the standard versions gives none of these errors. The import version also gets more rate matrix determination problems, indicating there are two photons in a cell.

Surprisingly, x and z are not the same for the cv model and the imported version. Why is this? It's because windsave2table prints out xcen, zcen, in the 2d case. Sign, for in the normal world we probably do want the center positions but for diagnosing our import models, we do not. Either, we need to write out a new file or agree we want the edge positions. Velocities are printed out at the edges in windsave2table.

[jhmatthews]

Perhaps the solution for the moment is to include a flag in windsave2table such as --edge?

It seems also that it might work better in Lee & Kelly's case, so I will see if I can get their model working tomorrow when i work on this again.

[kslong]

That sounds reasonable, or maybe just modify it to add another two columns to the master table.

[kslong]

windsave2table has now been modified so that it prints out the edge positions as well as the center position into the so-called master.txt table. (The other tables remain the same.). After doing this, I have verified that one can round-trip everything, that is create a master table, run the 2d conversion routine to create an import file, run python, and then run windsave2table again, and the coordinates for the cv model are the same as before.

We are still getting some errors that do not appear in the cv runs (which, if I had to guess are associated with the absence of windcones).

The issue that remains for imported models is what to do about cells which are partially in the wind. As currently implemented these are excluded from the wind. We could, instead include them in python. We could use the scripts which create an import file to fix things up how ever one wants. We could use, add a rdpar parameter, which appears if the import file has partially filled cells, to ask the user what she wants to do.

[kslong]

Sigh, somewhere along the way we, probably Knox, have broken py for rtheta coordinates, never mind trying to import it.

rtheta.pf.txt

This input file works fine under dev branch, but fails on the import branch. I found this when I was beginning to work on importing models in rtheta coordinates

This is the offending error:

100001 -- translate_in_wind: Photon is in cell %d with negligible volume, moving photon %.2e Occurrences %d

[kslong]

Tsk, Tsk, @jhmatthews This one was on you. You created a new initialization for wind_cells, W_NOT_ASSIGNED, but did not check what this implied for creating volumes in rtheta.c. Should be fixed now.

[jhmatthews]

This looks much better, and thanks for sorting my oversight.

One issue I've spotted is that there are on the order of 1% errors in the volume, and this is simply because we print out to 2 significant figures in windsave2table with %8.2e. I'd suggest we make this more accurate (4sf?)

Once that's done, I think this initial work can probably be merged into dev but with additional testing going forward. Do you agree?

I'm going to try and convert Kelly's model into a readable format now.

[kslong]

No worries, James. I am going to try to finish up the rtheta import today. I agree that we should merge back to dev soon. Let’s plan on your doing this on Wednesday, if that is OK.

I had also noticed that the volumes did not match, but did not understand why. I will fix the windsave2table formats

[kslong]

The rtheta import model is now working at least in principle. Much more testing is required for this, and, given that I found an error in the cylindrical import along the way, in all of the various import models. I updated the minor version number so we can track this addition.

Nevertheless @jhmatthews I think it is good enough that am now happy for you to merge all of this into dev, and any further development can continue there.

I still need to worry about the number of significant figures in windsave2table. The only real problem with adding significant figures is that it becomes hard to read the tables because they are too wide.

[kslong]

Although we thought we had this working, the new input files we got from Kelly, which are cylindrical models do not work at all. No photons are found in the wind. The models we tested on created by python for a cylindrical grid still do work.

What appears to be happening in the FU Ori model is that once it gets into the routine ds_to_wind the photon is translated for directly out of the wind region.

Inspection shows that the difference between the input models for FU Ori and the CV model is that the cv model has a data point in the grid at 0,9 whereas the FUori model starts at a large positive value. This means photons in the cv case start out in the wind, and have the possibility of never encountering ds_to_wind (which gives a infinity result)

Inspection also shows that the wind cones that are created in the fuori and cv model case are both identical. It does not look as if this was intended. We do calculate the wind_rho_max and min in the import routines, and in the fuori case these rho_min is defined.

What is happening it appears that setup_wind_cones is being called (from python.c) before this is done however. The reason for this I believe is that in the normal kinematic model case we want the windcones before we set up the wind.

So there are two possible solutions. (a) Run setup_wind cones again for the import cases, (b) provide the information earlier. (a) is a kluge, but nevertheless I tried this, and it did not quite work. That is because when theta min and theta max are 0, it is treated as an odd special case.

Note that the wind_cones are defined by the intercept of the cone with the z axis and an angle. When the angle is zero, the intercept is infinite.

Note that ds_to_wind does not really check for anything but the wind_cones and the star, except for the special case of a corona, where it looks for the pillbox.

[kslong]

So this is quite complicated, @jhmatthews I've created a new branch import for solving this problem. I added code to ds_to_wind that stops the photon at the edge of a cylindrical region to ds_to_wind, which gets translate_in_space to stop at the edge of the cylindrical region. But this is not enough because on the second pass through where_in_wind says the photon is not in the wind, because Kelly's import file has only a part of the cylindrical region filled.

These means, see below, that where_in_wind returns not in wind, even though it is in one of the cells.

58 if (where_in_wind (pp->x, &ndomain) < 0) 59 { 60 istat = translate_in_space (pp); 61 } 62 else if ((pp->grid = where_in_grid (ndomain, pp->x)) >= 0) 63 { 64 istat = translate_in_wind (w, pp, tau_scat, tau, nres); 65 } 66 else 67 { 68 istat = pp->istat = -1; / It's not in the wind and it's not in the grid. Bummer! / 69 Error ("translate: Found photon that was not in wind or grid, istat %i\n", where_in_wind (pp->x, &ndomain)); 70 }

So either we have to add additional code to translate_in_wind, so that it moves the photon to a cell that is in the wind, or we have to fix the logic here so that the the program understands that to treat all cells as if they are in the wind for translating from place to place.

These problems did not arise with our kinematic models, because we use the wind cones to populate the inwind variable of wmain.

Aside: the moving of photons in cells that are not in the wind may be a performance hit. To reduce this one might try to find the "best" wind_cones to use. Exactly how to do that thoug wouild require real thought.

[kslong]

@jhmatthews I've taken the point of view that we want translate_in_wind to translate the photon until it finds a cell that is in the wind. It works, more or less for the self-generated models, but not for Kelly's example because only some of the cells are actually in the wind.

However, this is proving really complicated for the general case, because some of our routines cannot handle the situation where a cell is in the grid but not in the wind. We never needed to move photons through such cells previously.

Routines like where_in_grid to take a domain number, which is good

but routines like cylind_ds_in_cell, which just take a photon_ptr, and get ndom indirectly through wmain ... , e.g

ndom = wmain[p->grid].ndom;

But if you are not in the wind region, this number is a set to a negative value.

Aside: An alternative approach here although it would mean starting over would be in imported models to consider all cells to be in the wind, but prevent them from doing anything except pass through a cell.

[kslong]

@jhmatthews I have gotten a version of the program to actually get photons to go through the wind, but I had to suppress some errors to make this happen. If the program is working, which is questionable, then the reason the errors have to be suppressed is that there are rays that get into the cylindrical grid but never actually encounter a grid cell that is in the wind.

[kslong]

@jhmatthews I spent more time on this today. I fixed at least one problem. Photons now go through all of the wind regions one expects; previously they only went through a part of the wind. I modified the command line print our of nphotons so that the wind region is clearly indicated.

But I am sill getting lots of errors having to do with photons not in the grid. It's possible that these photons are now simply in "space" and that the program is doing what it is supposed to be doing, and one just needs to capture this a little better. The location where photons are throwing errors, looks to be at the outside edges where we expect photons to be leaving, which suggests that may be the case. Green here is Kelly's grid. Blue are the location of "problem photons"

but actually verifying this has eluded me, so I am reluctant to simply suppress the errors.

One possible way to test that everything is working is to try two domain one which is a normal say sv model which is bigger than the imported model. I tried this though and if failed or reasons that are not obviously connected to the problem at hand.

[kslong]

One way that occurred to me to test whether things were working as far as getting things throught the plasma was concerned, was to take a model replacing Kelly's densities with low values, and to compare this to a model with the same radiation sources, but, say a KWD model, also with very low densities. I tried this and the results were very similar. This suggests we may be OK.

[kslong]

I've replaced the errors that when the import model gave out_of_grid warnings with the ability to track photons as they go through the wind, in advanced move with the @Diag.save_photons option. (Note that this produces a very large file for even a small run, and I had to edit it down to be able to plot it in a finite period of time.)

Photons do not seem to be stopping where they should not. The next figure shows a sample of the final place photons exit the calculation.

There are two populations one at large distance and one at a small distance. Honing in on the smaller distance one finds:

There is not indication that anything is getting stuck at the edge of the grid in these plots. To first order, I believe the transfer of photons through the grid is actually working.

There is some oddness if you look at the standard log file. You see more photons at the edge of the wind by a factor of 2 or so than in adjacent cells. This may be cause for some residual concern.

Also, it should be noted that all of this only works for cylindrical grids (but that's what we are concerned about right now).

Note that we are getting a fair number of other errors for the sample run, and it is unclear whether any of them are related to the simple problem of importing the model, and transporting photons:

Error summary: End of program
Recurrences --  Description
        7 -- getatomic_data: line input incomplete: %s
       16 -- wind_div_v: div v %e negative in cell %d Domain %d. Major problem if inwind (%d) == 0
       42 -- lucy_mazzali1: fudge>10 www %8.2e t_e %8.2e  t_r %8.2e 
       13 -- wind_cooling: xtotal emission %8.4e is < 0!
        1 -- disk_init: Got to ftot %e at r %e < rmax %e. OK if freqs are high
     1579 -- sane_check: %d %e
      553 -- trans_phot: sane_check photon %d has weight %e before extract
     1026 -- Radiation:sane_check photon weight is %e for tau %e
        4 -- error_count: This error will no longer be logged: %s
     4333 -- nnscat is %i for photon %i in scatter mode %i! nres %i NLINES %i

[kslong]

349 is a continuation of this.

[kslong]

Since work on this has moved to #349, ksl is closing this issue. @saultyevil when you get around to working on this you will want to read through this as well, since it illustrates the problems one has transporting photons through an arbitrary wind.