It is currently Wed Dec 11, 2024 9:11 am

All times are UTC - 8 hours [ DST ]




Post new topic Reply to topic  [ 9 posts ] 
Author Message
 Post subject: Non-uniform spatial grid
PostPosted: Thu Jan 21, 2016 6:13 pm 
Offline

Joined: Thu Apr 05, 2012 10:40 pm
Posts: 14
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!


Top
 Profile  
 
PostPosted: Thu Jan 21, 2016 11:40 pm 
Offline

Joined: Wed May 17, 2006 7:02 am
Posts: 285
Location: MPE Garching
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


Top
 Profile  
 
PostPosted: Fri Jan 22, 2016 1:18 pm 
Offline

Joined: Thu Apr 05, 2012 10:40 pm
Posts: 14
Very nice! I continue to be impressed by how extensible your team has designed this code!


Top
 Profile  
 
PostPosted: Fri Jan 22, 2016 7:44 pm 
Offline

Joined: Thu Apr 05, 2012 10:40 pm
Posts: 14
Seems to be working correctly with your recommended modifications. I will verify that the CR distributions match within reason soon.


Top
 Profile  
 
PostPosted: Fri Jan 22, 2016 11:31 pm 
Offline

Joined: Wed May 17, 2006 7:02 am
Posts: 285
Location: MPE Garching
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


Top
 Profile  
 
PostPosted: Fri Jan 29, 2016 1:46 pm 
Offline

Joined: Thu Apr 05, 2012 10:40 pm
Posts: 14
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...


Top
 Profile  
 
PostPosted: Fri Jan 29, 2016 11:26 pm 
Offline

Joined: Wed May 17, 2006 7:02 am
Posts: 285
Location: MPE Garching
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


Top
 Profile  
 
PostPosted: Sat Jan 30, 2016 5:52 am 
Offline

Joined: Wed May 17, 2006 7:02 am
Posts: 285
Location: MPE Garching
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


Top
 Profile  
 
PostPosted: Fri Feb 12, 2016 3:53 pm 
Offline

Joined: Thu Apr 05, 2012 10:40 pm
Posts: 14
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;


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 9 posts ] 

All times are UTC - 8 hours [ DST ]


Who is online

Users browsing this forum: No registered users and 0 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
Powered by phpBB® Forum Software © phpBB Group