If not too involved, could someone summarize what would need to be done in order to run non-uniform grids? Would be extremely nice to simulate e.g. the Galactic center region at higher resolution. I am familiar with the internals, but it would be nice to have an expert overview before diving in.

Thanks!

It would be quite straightforward. The grids are defined in Galaxy.h which defines the Galaxy class

double* x; // x grid
double* y; // y grid
double* z; // z grid
double* r; // r grid
There is no assumption of a uniform grid in the propagation code, so these can be set to the required values. The values are set in Galaxy.init() in Galaxy.cc which is invoked at the start of GALPROP. At present they are set to a uniform grid using the parameters x_min, x_max, dx etc, but they can be redefined by modifying Galaxy.init().

Then also the grid dimensions must be set:

int n_rgrid; // number of points in radius (2D)
int n_zgrid; // number of points in z (1D,2D)
int n_xgrid; // number of points in x (3D)
int n_ygrid; // number of points in y (3D)

That should be enough, but we have not tried it so it will be interesting to hear of results.

Make sure to use the latest version at http://sourceforge.net/projects/galprop (currently 528 downloads since the start!)

_________________
Andy Strong, MPE

Very nice! I continue to be impressed by how extensible your team has designed this code!

Seems to be working correctly with your recommended modifications. I will verify that the CR distributions match within reason soon.

Thanks, this is very useful. with your feedback I could introduce this into a future release.
It would be necessary to think about how to specify the variable grid in the parameter file.
I note that DRAGON already has this, so that might also be a useful guide.

_________________
Andy Strong, MPE

I see that particle.dx,dy,dz are used in propel.cc when calculating the Crank-Nicolson coefficients for diffusion and convection. This does impose uniform grid spacing, but they can just be calculated inside propel.cc using the galaxies spacing between grid elements in the relevant directions. At that point, one just needs to modify the assignment of alpha_1,2,3 to use the appropriate dx.

grep of the source folder for particle.dx, particle.dy, particle.dz reveals that this is the only place it is called.

As for the galprop spec, I am currently just adding parameters to specify a portion of each dimension to simulate in high-res via

xmin_inner
xmax_inner
dx_inner
etc..

But of course if we simulate from -2-2 kpc in x, it is for all y and z. Since the grids are currently spaced along each axis only, there are limitations like this, and it would be just a bit more work to actually simulate only a box along each dimension. One would need to specify specify the grid spacing in each dimension as a function of x,y,z.

A more general way of specifying grid spacing might be to assume a default spacing and then add in something like how cr_sources are added, but specifying multiple boxes. If they overlap this could get messy. For my purposes I only care about the GC region...

You are right about propel.cc using dx dy dz, and as you say it is easy fo replace them by the computed x[i+1]-x[i] etc.

For the grid, having regions is ok but a bit restricted. A continuous function could be one solution. Also one should really think about the principles of grid spacing from the numerical point of view, to optimize it.

_________________
Andy Strong, MPE

It is probably not worth going too far with this since GALPROP is quite basic in its approach, while modern methods like in PICARD will eventually take over. Maybe they will eventually use adaptive grids as in MHD codes.

_________________
Andy Strong, MPE

Just for further completeness, one also needs to modify places where galaxy.dx is present. Namely in skymap generation routines:



grep galaxy.dx *.cc
gen_bremss_ionized_skymap.cc: ix=(int)((xx-galaxy.x_min)/galaxy.dx + 0.5);//IMOS20060420
gen_bremss_skymap.cc: ix=(int)((xx-galaxy.x_min)/galaxy.dx);
gen_DM_skymap.cc: ix=(int)((xx-galaxy.x_min)/galaxy.dx);
gen_DM_source.cc: +=pow(DM_profile_av(galaxy.r[ix], galaxy.r[iy], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz),2)
gen_DM_source.cc: pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz)/DMmass,2)
gen_DM_source.cc: *pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz)/DMmass,2)
gen_IC_skymap.cc: ix=(int)((xx - galaxy.x_min)/galaxy.dx);
gen_IC_skymap.cc: ix=(int)((xx-galaxy.x_min)/galaxy.dx + 0.5);
gen_pi0_decay_skymap.cc: ix= ( int ) ( ( xx-galaxy.x_min ) /galaxy.dx );
gen_synch_skymap.cc: ix=(int)((xx-galaxy.x_min)/galaxy.dx + 0.5);//IMOS20060420
gen_synch_skymap.cc: ix=(int)((xx-galaxy.x_min)/galaxy.dx + 0.5);//IMOS20060420
gen_synch_skymap.cc: galaxy.dx, galaxy.dy, galaxy.dz,
gen_synch_skymap.cc: galaxy.dx, galaxy.dy, galaxy.dz,
gen_synch_skymap.cc: galaxy.dx, galaxy.dy, galaxy.dz,
isrf_energy_density.cc: int ix=(int)((rr-galaxy.x_min)/galaxy.dx + 0.5); //IMOS20060412 corrected "0.5"
store_bremss_emiss.cc: if(galaxy.n_spatial_dimensions==3)cdelt1=galaxy.dx;
store_DM_emiss.cc: if(galaxy.n_spatial_dimensions==3) cdelt1=galaxy.dx;
store_IC_emiss.cc: double cdelt1 = (galaxy.n_spatial_dimensions == 2 ? galaxy.dr : galaxy.dx);
store_ionization_rate.cc: if(galaxy.n_spatial_dimensions==3)cdelt1=galaxy.dx;
store_pi0_decay_emiss.cc: cdelt1 = galaxy.dx;
store_synch_emiss.cc: cdelt1 = galaxy.dx;