GALPROP https://galprop.stanford.edu/forum/ |
|
Non-uniform spatial grid https://galprop.stanford.edu/forum/viewtopic.php?f=19&t=225 |
Page 1 of 1 |
Author: | erccarls [ Thu Jan 21, 2016 6:13 pm ] |
Post subject: | Non-uniform spatial grid |
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! |
Author: | strong [ Thu Jan 21, 2016 11:40 pm ] |
Post subject: | Re: Non-uniform spatial grid |
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!) |
Author: | erccarls [ Fri Jan 22, 2016 1:18 pm ] |
Post subject: | Re: Non-uniform spatial grid |
Very nice! I continue to be impressed by how extensible your team has designed this code! |
Author: | erccarls [ Fri Jan 22, 2016 7:44 pm ] |
Post subject: | Re: Non-uniform spatial grid |
Seems to be working correctly with your recommended modifications. I will verify that the CR distributions match within reason soon. |
Author: | strong [ Fri Jan 22, 2016 11:31 pm ] |
Post subject: | Re: Non-uniform spatial grid |
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. |
Author: | erccarls [ Fri Jan 29, 2016 1:46 pm ] |
Post subject: | Re: Non-uniform spatial grid |
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... |
Author: | strong [ Fri Jan 29, 2016 11:26 pm ] |
Post subject: | Re: Non-uniform spatial grid |
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. |
Author: | strong [ Sat Jan 30, 2016 5:52 am ] |
Post subject: | Re: Non-uniform spatial grid |
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. |
Author: | erccarls [ Fri Feb 12, 2016 3:53 pm ] |
Post subject: | Re: Non-uniform spatial grid |
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; |
Page 1 of 1 | All times are UTC - 8 hours [ DST ] |
Powered by phpBB® Forum Software © phpBB Group https://www.phpbb.com/ |