Poisson solver via spectral methods within the shearingsheet. More...
Data Types | |
| type | gravity_sboxspectral |
Functions/Subroutines | |
| subroutine | updategravity_single (this, Mesh, Physics, Fluxes, pvar, time, dt) |
| Calculates the acceleration from potential. More... | |
| subroutine | fft_forward (this, Mesh, Physics) |
| Calculates the FFT forward. More... | |
Variables | |
| real, parameter | sqrttwopi = 2.50662827463100050241576528481104525300698674 |
spectral poisson solver shearing box | |
| subroutine | initgravity_sboxspectral (this, Mesh, Physics, config, IO) |
| Constructor of gravity sboxspectral module. More... | |
| subroutine | calcpotential (this, Mesh, Physics, time, pvar) |
| Computes the potential with FFT method within a shearingsheet. More... | |
| subroutine | fft_backward (this, Mesh, Physics) |
| Calculates the FFT backward transformation. More... | |
| subroutine | calcdiskheight_single (this, Mesh, Physics, pvar, bccsound, h_ext, height) |
| Compute disk pressure scale height for geometrically thin self-gravitating shearingsheet. More... | |
| subroutine | infogravity (this, Mesh) |
| Prints out information. More... | |
| subroutine | fieldshift (this, Mesh, Physics, delt, field, mass2D) |
| Shifts the whole field to the next periodic point. More... | |
| subroutine | finalize (this) |
| Closes the gravity term of the shearingsheet spectral solver. More... | |
Detailed Description
Poisson solver via spectral methods within the shearingsheet.
Program and data initialization for a shearing-box simulation
There are currently two different implementations. Please be aware that this module requires a certain domain decomposition in pencils and cannot be chosen arbitrarily.
- FFTW on one core
- FFTW distributed parallel
Function/Subroutine Documentation
◆ calcdiskheight_single()
| subroutine gravity_sboxspectral_mod::calcdiskheight_single | ( | class(gravity_sboxspectral), 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 shearingsheet.
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 & (2) \\ \rho_c &=& 1/\sqrt{2\pi} \Sigma / h & (3) \\ 1/h_{ext}^2 &=& \Omega_{K}^2 / c_{s}^2 \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 (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. \( \Phi_{sg} \) is the gravitational potential in the shearing sheet. Replacing \( \rho_c \) in (2) by (3) and inserting the result in (1) together with (4) yields a quadratic equation with the two solutions
\[ h = -p \pm \sqrt{q+p^2} = p (-1 \pm \sqrt{1 + q/p^2}) \]
where \( p = \sqrt{2\pi} G \Sigma / \Omega_{K}^2 \) and \( q = c_{s}^2 / \Omega_{K}^2 \). The sign of the root has to be positive, because in the non self- gravitating limit \( p=0, q=1 \) and \( h \) needs to be \( c_{s}/\Omega_{K} \) not \( -c_{s}/\Omega_{K} \).
Definition at line 839 of file gravity_sboxspectral.f90.
◆ calcpotential()
| subroutine gravity_sboxspectral_mod::calcpotential | ( | class(gravity_sboxspectral), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| real, intent(in) | time, | ||
| class(marray_compound), intent(inout) | pvar | ||
| ) |
Computes the potential with FFT method within a shearingsheet.
Calculates Poisson's equation:
\[ \Delta \Phi = 4 \pi G \Sigma \delta (z), \]
where \( \delta (z) \) is the \( \delta \)-distribution, thus for razor-thin disks.
The implemenation is done by the following steps:
- Shift surface density field to next periodic point
\[ \Sigma(x,y) \longrightarrow \Sigma(x,y'). \]
See gravity_sboxspectral.fieldshift for more informations. - FFT of shifted surface density field
\[ \Sigma(x,y') \longrightarrow \Sigma(\mathbf{k}) \]
- Calculate
\[ \Phi(\mathbf{k}) = -\frac{2\pi G}{|\mathbf{k}|}\Sigma(\mathbf{k}) e^{i\mathbf{k\cdot x - |kz|}}. \]
- FFT \(^{-1}\) of Potential
\[ \Phi(\mathbf{k}) \longrightarrow \Phi(x,y') \]
- Backshift of Potential
\[ \Phi(x,y') \longrightarrow \Phi(x,y). \]
The acceleration is eventually calculated in gravity_sboxspectral.updategravity_single.
- Todo:
- Violation of best practices. It should take Fmass2D and mass2D as dummy argument, but this is different for every combination NEC$ IEXPAND
- Todo:
- Violation of best practices. It should take Fmass2D and mass2D as dummy argument, but this is different for every combination. NEC$ EXPAND
Definition at line 433 of file gravity_sboxspectral.f90.
◆ fft_backward()
|
private |
Calculates the FFT backward transformation.
Definition at line 786 of file gravity_sboxspectral.f90.
◆ fft_forward()
|
private |
Calculates the FFT forward.
Definition at line 737 of file gravity_sboxspectral.f90.
◆ fieldshift()
| subroutine gravity_sboxspectral_mod::fieldshift | ( | class(gravity_sboxspectral), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| real | delt, | ||
| real, dimension(mesh%igmin:mesh%igmax,mesh%jgmin:mesh%jgmax,mesh%kgmin:mesh%kgmax), intent(in) | field, | ||
| real, dimension(mesh%imin:mesh%imax,mesh%jmin:mesh%jmax), intent(out) | mass2D | ||
| ) |
Shifts the whole field to the next periodic point.
Implementation is done by
\[ y' = y - q \Omega x (t-t_p), \]
with \( t_p = \text{NINT}(q\Omega t) / (q\Omega) \). In order to map the continuous shift at the discrete field linear interpolation is used, and assumes periodic behaviour along the y-direction.
Definition at line 899 of file gravity_sboxspectral.f90.
◆ finalize()
| subroutine gravity_sboxspectral_mod::finalize | ( | class(gravity_sboxspectral), intent(inout) | this | ) |
Closes the gravity term of the shearingsheet spectral solver.
Definition at line 960 of file gravity_sboxspectral.f90.
◆ infogravity()
|
private |
Prints out information.
Definition at line 872 of file gravity_sboxspectral.f90.
◆ initgravity_sboxspectral()
| subroutine gravity_sboxspectral_mod::initgravity_sboxspectral | ( | class(gravity_sboxspectral), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| type(dict_typ), pointer | config, | ||
| type(dict_typ), pointer | IO | ||
| ) |
Constructor of gravity sboxspectral module.
This subroutine reads the necessary config data for the the specral gravity solver within the shearingbox. It initializes the gravity type and various mesh data arrays. Some of those are marked for output. The possible combinations are
- FFTW - serial, parallel
Definition at line 148 of file gravity_sboxspectral.f90.
◆ updategravity_single()
| subroutine gravity_sboxspectral_mod::updategravity_single | ( | class(gravity_sboxspectral), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| class(fluxes_base), intent(in) | Fluxes, | ||
| class(marray_compound), intent(inout) | pvar, | ||
| real, intent(in) | time, | ||
| real, intent(in) | dt | ||
| ) |
Calculates the acceleration from potential.
Calculates
\[ \mathbf{a} = -\nabla \Phi. \]
Uses second order symmetric difference quotient.
Definition at line 339 of file gravity_sboxspectral.f90.
Variable Documentation
◆ sqrttwopi
|
private |
Definition at line 91 of file gravity_sboxspectral.f90.
