poisson_cyl
Module
Description
$Id$
Quick access
- Variables:
apply_boundcond,b,bessel_grid,bsmooth,calc_potential,check_setup,compute_acceleration,construct_large_grid,copy_density_large_grid,copy_potential_small_grid,coslt,decide_fourier_routine,dkr,dphi,dr,dr1,dth,dth1,du,dxc1,dxlt,dyc1,dz1,fft_fouriergrid,fourier_cosine_terms,fourier_to_physical_proc,fourier_to_physical_procs,gauss_seidel_iterate,generate_coordinates,generate_fourier_density,generate_kernals,generate_massfields,get_communication_matrix,gnewton,gphi,gr,innerradius,inverse_laplacian_expandgrid,inverse_laplacian_fft,inverse_laplacian_isoz,inverse_laplacian_logradial_fft,ipoisson_method,iproc_fourier,ipx_fourier,ipy_fourier,iregion,irl,iroot,iru,iteration_threshold,ixlower_fourier,ixupper_fourier,iylower_fourier,iyupper_fourier,jrl,jru,kkx_fft,kky_fft,kphi,kphi_fft_imag,kphi_fft_real,kr_fft,kr_fft_imag,kr_fft_real,krl,kru,legendre_qmod,lexpand_grid,lisoz,lklimit_shear,lkmax,lkmin,lmakecartoon,lnorepeatsumming,lprecalcdists,lrazor_thin,lreadoctree,lsemispectral,lshowtime,lsolve_bessel,lsolve_cyl2cart,lsolve_relax_sor,lsquareregions,ltreestatus,luseprevioussum,lwriteoctree,m_fft,mmax,ngroup,niter_poisson,nkhgrid,nkr,nkt,nktgrid,nlt,nnx,nny,nprecalc,nsingle,nslice,nth,nthgrid,octree_maxdist,octree_smoothdist,octree_theta,phi2d,phi2d_global,phi_previous_step,physical_to_fourier_proc,physical_to_fourier_procs,r2_ext,r2_int,rad1,refine,regdist1_group,regdist1_single,regsmooth_group,regsmooth_single,remap_to_pencil_fouriergrid,residual,restrict,rhs_previous_step,rjac,rn,sinlt,smooth_full_weight,sphi,sphi_fft_imag,sphi_fft_real,sqrtrad_1,sr,sr_fft_imag,sr_fft_real,stop_fatal,themap_group,themap_single,theta0,theta1,tht,transp_pencil_fouriergrid,trilinear_interpolate,u2d,u2d_global,unmap_from_pencil_fouriergrid,v_cycle,vols,xc0,xlt,xmesh,xrecv,yc0,ymesh,yrecv,zed,zmesh,zrecv- Routines:
calculate_cross_bessel_functions(),calculate_cross_legendre_functions(),do_barneshut(),five_point_solver(),get_acceleration(),get_border_values(),get_dist(),get_serial_array(),initialize_poisson(),integrate_border(),inverse_laplacian(),inverse_laplacian_bessel(),inverse_laplacian_cyl2cart(),inverse_laplacian_fft_z(),inverse_laplacian_semispectral(),inverse_laplacian_sor(),inverse_laplacian_z_2nd_neumann(),mkmap(),read_octree(),read_poisson_init_pars(),read_poisson_run_pars(),roundtwo(),sincoslf(),write_octree(),write_poisson_init_pars(),write_poisson_run_pars()
Needed modules
Variables
- poisson/dkr [real]
- poisson/dr [real]
- poisson/dr1 [real]
- poisson/dth [real]
- poisson/dth1 [real]
- poisson/dz1 [real]
- poisson/ipoisson_method [character,optional/default='nothing']
- poisson/iteration_threshold [real,optional/default=0.001]
- poisson/kmax [real,optional/default=0.0]
- poisson/lrazor_thin [logical,optional/default=.true.]
- poisson/lsolve_bessel [logical,optional/default=.false.]
- poisson/lsolve_cyl2cart [logical,optional/default=.false.]
- poisson/lsolve_relax_sor [logical,optional/default=.false.]
- poisson/mmax [integer,parameter=8]
- poisson/nkhgrid [integer]
- poisson/nkr [integer]
- poisson/nkt [integer]
- poisson/nktgrid [integer]
- poisson/nr [integer]
- poisson/nth [integer]
- poisson/nthgrid [integer]
- poisson/r0 [real]
- poisson/rjac [real,optional/default=1.0]
- poisson/rn [real]
- poisson/theta0 [real]
- poisson/theta1 [real]
Subroutines and functions
- subroutine poisson/initialize_poisson()
Perform any post-parameter-read initialization i.e. calculate derived parameters.
18-oct-07/anders: adapted
- Call to:
fatal_error(),calculate_cross_bessel_functions(),calculate_cross_legendre_functions(),inverse_laplacian_bessel(),inverse_laplacian_cyl2cart(),inverse_laplacian_sor(),fourier_transform_xy_xy_other(),fourier_transform_y(),get_border_values(),five_point_solver(),integrate_border(),get_serial_array(),stop_it()
- subroutine poisson/inverse_laplacian(f, phi)
Dispatch solving the Poisson equation to inverse_laplacian_fft or inverse_laplacian_semispectral, based on the boundary conditions
17-jul-2007/wolf: coded wrapper
- subroutine poisson/inverse_laplacian_cyl2cart(phi)
Solve the 2D Poisson equation in cylindrical coordinates by transforming to a periodic cartesian grid before Fourier transforming.
This subroutine is faster than inverse_laplacian_bessel for low resolutions and low number of processors. But it gets slower when increasing any of them, due to the great amount of communication used.
The frequent broadcast of a big array gave problems at the PIA cluster in Heidelberg after some thousands of time-steps. The problem was probably due to memory leaking. No problem occured at the UPPMAX cluster in Uppsala. So beware that using this broadcasting extravaganza subroutine might not work in all clusters.
01-12-07/wlad: coded 28-02-08/wlad: merged the serial and mpi versions
- Parameters:
- Use :
- Call to:
fatal_error(),fourier_transform_xy_xy_other(),fourier_transform_y(),get_border_values(),five_point_solver(),integrate_border(),get_serial_array(),stop_it()
- subroutine poisson/inverse_laplacian_bessel(phi)
Solve the 2D Poisson equation in cylindrical coordinates
This beautiful and elegant theory for calculating the potential of thin disks using bessel functions and hankel transforms is not so useful numerically because of the ammount of integrations to perform.
In any case, we managed to optimize the computations so that only 30% of the computational time is spent in this routine. The number is only slightly dependant on parallelization and resolution, as opposed to inverse_laplacian_cyl2cart.
For example, having this
- do ikr=1,nkr
- SS(ikr)=sum(bessel_grid(2:nr-1,ikr,ikt)*sigma_tilde_rad(2:nr-1))+&
- .5*(bessel_grid( 1,ikr,ikt)*sigma_tilde_rad(1) +&
bessel_grid(nr,ikr,ikt)*sigma_tilde_rad(nr))
enddo
instead of
- do ikr=1,nkr
tmp=bessel_grid(:,ikr,ikt)*sigma_tilde*rad SS(ikr)=sum(tmp(2:nr-1))+.5*(tmp(1)+tmp(nr))
enddo
lead to a factor 2 speed up. That’s a lot for a subroutine that usually takes most of the computing time. So, think (twice) before you modify it.
At every wavelength, the density-potential pair in fourier space is
Phi_tilde_k = exp(-k|z|) Jm(k*r) ; Sigma_tilde_k =-k/(2piG) Jm(k*r)
So, the total potential and the total spectral density is
Phi_tilde = Int [ S(k) Jm(k*r) exp(-k|z|) ]dk
Sigma_tilde=-1/(2piG) Int [ S(k) Jm(k*r) k] dk
The function above is the Hankel transform of S(k), so S(k) is the inverse Hankel transform of Sigma_tilde
S(k)=-2piG Int[ Jm(k*r) Sigma_tilde(r,m) r] dr
06-03-08/wlad: coded
- Parameters:
- Use :
- Call to:
fatal_error(),fourier_transform_y(),get_border_values(),five_point_solver(),integrate_border(),get_serial_array(),stop_it()
- subroutine poisson/inverse_laplacian_sor(f, phi)
Solve the 3D Poisson equation in cylindrical coordinates by using mesh-relaxation with the SOR algorithm. The borders have to be pre-computed. For this, we use the analytical expression of Cohl and Tohline (1999) using harmonic expansions in Legendre functions.
This is very expensive and does not need to be done every time. Just when the potential has changed by some small threshold
23-04-08/wlad: coded
- Parameters:
- Use :
- Call to:
fatal_error(),get_border_values(),fourier_transform_y(),five_point_solver(),integrate_border(),get_serial_array(),stop_it()
- subroutine poisson/five_point_solver(lhs, rhs, a_band, b_band, b_band1, c_band, d_band, e_band, lupdate, lverbose)
Invert a five point matrix
- u(i,j) = 1/b*(a*u(i-1,j) + c*u(i+1,j) + d*u(i,j+1) +
e*u(i,j-1) + b*u(i,j) - rhs(i,j))
using Chebychev (checkboard, red and black) acceleration
- Parameters:
- Called from:
inverse_laplacian_cyl2cart(),inverse_laplacian_bessel(),inverse_laplacian_sor()- Call to:
fatal_error(),integrate_border(),get_serial_array(),stop_it()
- subroutine poisson/get_border_values(rhs, phi, lupdate_grid, lverbose)
Calculate the value of the potential in the borders of the grid
28-apr-08/wlad: coded
- Parameters:
- Use :
- Called from:
inverse_laplacian_cyl2cart(),inverse_laplacian_bessel(),inverse_laplacian_sor()- Call to:
- subroutine poisson/integrate_border(rhs_serial, ir, ith, iz, potential)
Get the potential in the borders by direct integration (see routine calculate_cross_legendre_functions below)
- Phi(r,tht,z) = -G/(pi*sqrt(r)) * &
Int r’dr’dz’dth’/sqrt(r’) * rho’ * & Sum_m=0^infty eps(m)*Q_(m-1/2)*cos(m*(tht-tht’))
_ infty
G / d3x’ rho(x’) ___
- Phi(x) = _ _______ | ___________ e(m) Q (chi) cos[m(t-t’)]
pi /r _/ /r’ m=0
This is time consuming and only done for the border, if they need updating. In the grid, we iteratively solve the discrete Poisson equation via Chebychev acceleration.
28-apr-08/wlad: coded
- Parameters:
- Called from:
inverse_laplacian_cyl2cart(),inverse_laplacian_bessel(),inverse_laplacian_sor(),five_point_solver(),get_border_values()- Call to:
- subroutine poisson/calculate_cross_bessel_functions()
Calculate the Bessel functions related to the cylindrical grid
bessel_grid(ir,ikr,m)=J_m(kr(ikr)*rad(ir))
06-03-08/wlad: coded
- Use :
- Call to:
- subroutine poisson/calculate_cross_legendre_functions()
Calculate the Legendre functions of second kind related to the cylindrical grid. They are needed to solve the potential by the analytical expression (Cohl & Tohline 1999)
- Phi(r,tht,z) = -G/(pi*sqrt(r)) * &
Int r’dr’dz’dth’/sqrt(r’) * rho’ * & Sum_m=0^infty eps(m)*Q_(m-1/2)*cos(m*(tht-tht’))
_ infty
G / d3x’ rho(x’) ___
- Phi(x) = _ _______ | ___________ e(m) Q (chi) cos[m(t-t’)]
pi /r _/ /r’ m=0
The Neumann factor e(m) is 1 for m=0 and 2 for other harmonics, chi is a function of the grid
chi= [r^2 + r’^2 + (z-z’)^2]/(2rr’)
And the Legendre functions can be found by a recurrency relation
Q_(-1/2) = mu*K(mu) Q_( 1/2) = chi*mu*K(mu) - (1+chi)*mu*E(mu)
4*(m-1) 2m-3
Q_(m-1/2) = _______ chi*Q_(m-3/2) - ______ Q_(m-5/2)
2m-1 2m-1
where mu=sqrt(2./(1+chi)) and K and E are the first and second kind complete elliptical integrals. The summation is truncated at a maximum harmonic mmax
28-apr-08/wlad: coded
- Use :
- Call to:
- subroutine poisson/get_serial_array(array, array_serial, var)
Feed in the array present on the processor, and return (broadcast) a serial array with the elements of that
- Parameters:
array (*) [real] ::
in all processors. Used for getting all y or z elements.
28-apr-08/wlad: coded
array_serial (*) [real]
var [character]
- Use :
- Called from:
inverse_laplacian_cyl2cart(),inverse_laplacian_bessel(),inverse_laplacian_sor(),five_point_solver(),get_border_values(),integrate_border(),calculate_cross_bessel_functions(),calculate_cross_legendre_functions()- Call to:
- subroutine poisson/read_poisson_init_pars(iomsg)
- Parameters:
iomsg [character,out]
- Use :
- subroutine poisson/write_poisson_init_pars(unit)
- Parameters:
unit [integer,in]
- subroutine poisson/read_poisson_run_pars(iomsg)
- Parameters:
iomsg [character,out]
- Use :
- subroutine poisson/write_poisson_run_pars(unit)
- Parameters:
unit [integer,in]