2D poisson solver using spectral methods for direct integration More...
Data Types | |
| type | gravity_spectral |
Functions/Subroutines | |
| elemental real function | greenfunction (dr2, green, sigma) |
| Green's function. More... | |
Variables | |
| real, parameter | sqrttwopi = 2.50662827463100050241576528481104525300698674 |
Variables in Parallel Mode | |
| subroutine | initgravity_spectral (this, Mesh, Physics, config, IO) |
| integer function | calcmcut (this, Mesh) |
| subroutine | calcpotential (this, Mesh, Physics, time, pvar) |
| subroutine | precomputei (this, Mesh, Physics) |
| Precomputes the fourier transform. More... | |
| subroutine | updategravity_single (this, Mesh, Physics, Fluxes, pvar, time, dt) |
| subroutine | calcdiskheight_single (this, Mesh, Physics, pvar, bccsound, h_ext, height) |
| compute disk pressure scale height for geometrically thin self-gravitating disks More... | |
| subroutine | finalize (this) |
Detailed Description
2D poisson solver using spectral methods for direct integration
This code implements the methods described in [chan2006] and [li2009] . It is a native 2D solver and works only in flat, polar geometries. This is why it should be paid attention to the dimension of the arrays. They where only allocated in 3D if necessary.
Function/Subroutine Documentation
◆ calcdiskheight_single()
| subroutine gravity_spectral_mod::calcdiskheight_single | ( | class(gravity_spectral), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| class(marray_compound), intent(inout) | pvar, | ||
| type(marray_base), intent(inout) | bccsound, | ||
| type(marray_base), intent(inout) | h_ext, | ||
| type(marray_base), intent(inout) | height | ||
| ) |
compute disk pressure scale height for geometrically thin self-gravitating disks
The algorithm solves the equations
\begin{eqnarray} 1/h^2 &=& 1/h_{sg}^2 + 1/h_{ext}^2 & (1)\\ \left.\frac{\partial^2\Phi_{sg}}{\partial z^2}\right|_{z=0} = (c_s/h_{sg})^2 &=& 4\pi G \rho_c - \left.\Delta_{\xi\eta} \Phi_{sg}\right|_{z=0} & (2) \\ \rho_c &=& 1/\sqrt{2\pi} \Sigma / h & (3) \end{eqnarray}
for the disk scale height \( h \) where \( h_{sg} \) is the self-gravitating pressure scale height due to the material in the disk and \( h_{ext} \) the scale height caused by an external (probably point mass) potential. \( \Sigma \) is the surface density, \( c_s \) is the speed of sound and \( \rho_c \) is the density in the equatorial plane. \( \Delta_{\xi\eta} \) denotes the 2D Laplacian in curvilinear coordinates and \( \Phi_{sg} \) is the gravitational potential of the disk. One can rewrite the last term in equation (2) replacing the Laplacian of the potential with the divergence of the 2D gravitational acceleration according to
\[ \left.\Delta_{\xi\eta} \Phi_{sg}\right|_{z=0} = -\nabla_{\xi\eta}\cdot\vec{g}_{sg}. \]
Replacing \( \rho_c \) in (2) by (3) and inserting the result in (1) yields a quadratic equation in \( h_{ext}/h \) with the two solutions
\[ h_{ext}/h = p \pm \sqrt{q+p^2} \]
where
\[ p = \sqrt{2\pi} G \Sigma h_{ext}/c_s^2 \]
and
\[ q = 1 + \left(h_{ext}/c_s\right)^2 + \nabla_{\xi\eta}\cdot\vec{g}_{sg}. \]
The sign of the root must be positive, because in the non self-gravitating limit \( p=0, q=1 \) and \( h_{ext}/h \) must become \( 1 \) not \( -1 \) .
- Attention
- This routine only works in 2D
Definition at line 688 of file gravity_spectral.f90.
◆ calcmcut()
|
private |
Definition at line 300 of file gravity_spectral.f90.
◆ calcpotential()
|
private |
Definition at line 335 of file gravity_spectral.f90.
◆ finalize()
|
private |
Definition at line 757 of file gravity_spectral.f90.
◆ greenfunction()
| elemental real function gravity_spectral_mod::greenfunction | ( | real, intent(in) | dr2, |
| integer, intent(in) | green, | ||
| real, intent(in) | sigma | ||
| ) |
Green's function.
The Green's function for an arbitrary vertical disk profile \( Z(r,z) \) is given by the integral
\[ G(r,r',\varphi) = -\int_{-\infty}^{\infty} \frac{Z(r',z')}{\sqrt{r^2+r'^2 - 2\,r\,r'\cos{(\varphi)} + \varepsilon^2 + z'^2}} dz'. \]
where \( \varepsilon \) is a smoothing parameter to avoid division by zero. This integral can be evaluated analytically, e. g., for a vertical Gaussian density distribution
\[ Z(r,z) = \frac{1}{\sqrt{2\pi H^2}} \exp{\left(-\tfrac{z^2}{2 H^2}\right)} \]
with pressure scale height \( H(r) \) depending only on the radial coordinate. In this case the Green's function is given by
\[ G(r,r',\varphi) = - \frac{\exp{\left(\tfrac{R^2}{4}\right)} K_0\left(\tfrac{R^2}{4}\right)}{\sqrt{2\pi} H(r')}, \quad\textsf{with}\quad R^2 = \frac{r^2 + r'^2 - 2\,r\,r'\cos{(\varphi)} + \varepsilon^2}{H(r')^2} \]
where \( K_0(x) \) is the modified Bessel function of the second kind of order 0.
Definition at line 471 of file gravity_spectral.f90.
◆ initgravity_spectral()
|
private |
Definition at line 122 of file gravity_spectral.f90.
◆ precomputei()
|
private |
Precomputes the fourier transform.
Precompute the fourier transform with respect to \( \varphi-\varphi' \) of
\[ I(r,r',\varphi-\varphi') = 2 \pi r' G(r,r',\varphi-\varphi'), \]
where \( G(r,r',\varphi-\varphi') \) is the softened Green's function (see greenfunction )
Definition at line 502 of file gravity_spectral.f90.
◆ updategravity_single()
|
private |
- Attention
- This routine only works in 2D
Definition at line 608 of file gravity_spectral.f90.
Variable Documentation
◆ sqrttwopi
|
private |
Definition at line 62 of file gravity_spectral.f90.
