Poisson solver using spectral methods within the shearingbox. More...
Data Types | |
| type | field_typ |
| type | gravity_sboxspectral |
Functions/Subroutines | |
| subroutine | initgravity_sboxspectral (this, Mesh, Physics, config, IO) |
| Constructor of gravity sboxspectral module. More... | |
| subroutine | setoutput (this, Mesh, Physics, config, IO) |
| subroutine | updategravity_single (this, Mesh, Physics, Fluxes, pvar, time, dt) |
| Calculates the acceleration from potential. More... | |
| subroutine | calcpotential (this, Mesh, Physics, time, pvar) |
| Computes the potential with FFT method within a shearingsheet. More... | |
| subroutine | setboundarydata (this, Mesh, Physics, delt) |
| Set ghost cell data for the potential. More... | |
| subroutine | fft_forward (this, Mesh, Physics, field) |
| Calculates the FFT forward. More... | |
| subroutine | fft_backward (this, Mesh, Physics, field) |
| 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 | fieldshift (this, Mesh, Physics, delt, field, shifted_field) |
| Shifts the whole field to the next periodic point. More... | |
| subroutine | finalize (this) |
| Closes the gravity term of the shearingsheet spectral solver. More... | |
Variables | |
| real, parameter | sqrttwopi = 2.50662827463100050241576528481104525300698674 |
Detailed Description
Poisson solver using spectral methods within the shearingbox.
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
- Attention
- The parallel execution of this spectral solver requires pencil decomposition of the computation domain, i.e., splitting is only allowed in 2nd and/or 3rd dimension. Check the key "decomposition" when setting the mesh properties (see Mesh ).
References
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.
For razor-thin 2D 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 & (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-selfgravitating limit \( p=0, q=1 \) and \( h \) needs to be \( c_s/\Omega_{K} \) not \( -c_s/\Omega_{K} \).
In the 3D case the midplane density \( \rho_c \) can be obtained directly from the data. Thus there is no need for equation (3) and the scale height is simply given by
\[ h = \frac{c_s}{\sqrt{4\pi G \rho_c + \Omega_{K}^2}}. \]
Definition at line 1020 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
Definition at line 522 of file gravity_sboxspectral.f90.
◆ fft_backward()
|
private |
Calculates the FFT backward transformation.
Definition at line 961 of file gravity_sboxspectral.f90.
◆ fft_forward()
|
private |
Calculates the FFT forward.
Definition at line 914 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, intent(in) | delt, | ||
| real, dimension(mesh%igmin:mesh%igmax,mesh%jgmin:mesh%jgmax,mesh%kgmin:mesh%kgmax), intent(in) | field, | ||
| type(field_typ), dimension(mesh%kmin-mesh%kp1:mesh%kmax+mesh%kp1), intent(out) | shifted_field | ||
| ) |
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 1084 of file gravity_sboxspectral.f90.
◆ finalize()
| subroutine gravity_sboxspectral_mod::finalize | ( | type(gravity_sboxspectral), intent(inout) | this | ) |
Closes the gravity term of the shearingsheet spectral solver.
Definition at line 1141 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 156 of file gravity_sboxspectral.f90.
◆ setboundarydata()
| subroutine gravity_sboxspectral_mod::setboundarydata | ( | class(gravity_sboxspectral), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| real, intent(in) | delt | ||
| ) |
Set ghost cell data for the potential.
Parallelization in x-y-plane is only supported for shear along 1st dimension, i.e., in west-east direction, so we need MPI communication only in this case.
Definition at line 738 of file gravity_sboxspectral.f90.
◆ setoutput()
| subroutine gravity_sboxspectral_mod::setoutput | ( | 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 | ||
| ) |
Definition at line 385 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 424 of file gravity_sboxspectral.f90.
Variable Documentation
◆ sqrttwopi
|
private |
Definition at line 83 of file gravity_sboxspectral.f90.
