|
PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
|
Numerical utilities. More...
Interfaces and Types | |
| interface | add_arr_mult |
| Add to an array (3) the product of arrays (1) and (2). More... | |
| interface | calc_det |
| Calculate determinant of a matrix. More... | |
| interface | calc_inv |
Calculate inverse of square matrix A. More... | |
| interface | calc_mult |
| Calculate matrix multiplication of two square matrices \(\overline{\text{AB}} = \overline{\text{A}} \ \overline{\text{B}}\). More... | |
| interface | conv_mat |
| Converts a (symmetric) matrix A with the storage convention described in eq_vars.eq_2_type. More... | |
| interface | calc_int |
| Integrates a function using the trapezoidal rule. More... | |
| interface | round_with_tol |
| Rounds an arry of values to limits, with a tolerance \(10^{-5}\) that can optionally be modified. More... | |
| interface | con2dis |
| Convert between points from a continuous grid to a discrete grid. More... | |
| interface | dis2con |
| Convert between points from a discrete grid to a continuous grid. More... | |
| interface | con |
| Either takes the complex conjugate of a square matrix element A, defined on a 3-D grid, or not. More... | |
| interface | bubble_sort |
| Sorting with the bubble sort routine. More... | |
| interface | order_per_fun |
| Order a periodic function to include \(0\ldots 2\pi\) and an overlap. More... | |
| interface | spline |
| Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the main priority. If spline representations are to be reused, manually use the library. More... | |
Functions/Subroutines | |
| subroutine, public | calc_aux_utilities (n_mod) |
| Initialize utilities for fast future reference, depending on program style. | |
| integer function, dimension(:,:), allocatable, public | derivs (order, dims) |
| Returns derivatives of certain order. | |
| integer function, public | check_deriv (deriv, max_deriv, sr_name) |
| checks whether the derivatives requested for a certain subroutine are valid | |
| integer function, public | calc_ext_var (ext_var, var, var_points, ext_point, deriv_in) |
| Extrapolates a function. | |
| integer function, public | calc_coeff_fin_diff (deriv, nr, ind, coeff) |
| Calculate the coefficients for finite differences. | |
| integer function, public | c (ij, sym, n, lim_n) |
| Convert 2-D coordinates (i,j) to the storage convention used in matrices. | |
| recursive integer function, public | fac (n) |
| Calculate factorial. | |
| integer function, public | is_sym (n, nn, sym) |
| Determines whether a matrix making use of the storage convention in eq_vars.eq_2_type is symmetric or not. | |
| recursive integer function, public | lcm (u, v) |
| Returns common multiple using the Euclid's algorithm. | |
| recursive integer function, public | gcd (u, v) |
| Returns least denominator using the GCD. | |
| subroutine, public | shift_f (al, bl, cl, a, b, c) |
| Calculate multiplication through shifting of fourier modes A and B into C. | |
| subroutine, public | solve_vand (n, a, b, x, transp) |
| Solve a Vandermonde system \(\overline{\text{A}} \ \vec{X} = \vec{B}\). | |
| integer function, public | calc_d2_smooth (fil_n, x, y, d2y, style) |
| Calculate second derivative with smoothing formula by Holoborodko, [9]. | |
Variables | |
| integer, dimension(:,:,:), allocatable, public | d |
| 1-D array indices of derivatives | |
| integer, dimension(:,:), allocatable, public | m |
| 1-D array indices of metric indices | |
| integer, dimension(:,:), allocatable, public | f |
| 1-D array indices of Fourier mode combination indices | |
| logical, public | debug_con2dis_reg = .false. |
| plot debug information for con2dis_reg() | |
| logical, public | debug_calc_coeff_fin_diff = .false. |
| plot debug information for calc_coeff_fin_diff() | |
Numerical utilities.
| integer function, public num_utilities::c | ( | integer, dimension(2), intent(in) | ij, |
| logical, intent(in) | sym, | ||
| integer, intent(in), optional | n, | ||
| integer, dimension(2,2), intent(in), optional | lim_n ) |
Convert 2-D coordinates (i,j) to the storage convention used in matrices.
Their size is by default taken to be 3:
\[ \left(\begin{array}{ccc} 1 & 4 & 7 \\ 2 & 5 & 8 \\ 3 & 6 & 9 \end{array}\right) \ \text{or} \ \left(\begin{array}{ccc} 1 & & \\ 2 & 4 & \\ 3 & 5 & 6 \end{array}\right) \ \text{for symmetic matrices} \ . \]
Optionally, this can be changed using n.
The value of \(c\) is then given by
\[c = (j-1) n + i\]
for non-symmetric matrices, and by
\[\begin{aligned} c &= (j-1)n + i - (j-1)\frac{j}{2} & \text{if} \ i > j \\ c &= (i-1)n + j - (i-1)\frac{i}{2} & \text{if} \ j > i \end{aligned}\]
since \(\sum_{k=1}^{j-1} k = \frac{(j-1)j}{2}\).
For local indices in submatrices, the limits in both dimensions have to be passed. The results for the full matrix are then subtracted by amount corresponding to the left, above and below parts with respect to the submatrix:
\[\sum_{i=1}^{m(2)-1} n-i+1 = (m(2)-1) \left(n+1-\frac{m(2)}{2}\right) \]
\[\sum_{i=1}^j \left(m(1)-m(2)+1-i\right) = j \left(m(1)-m(2)+\frac{1}{2} - \frac{j^*}{2}\right) \ , \ \text{with} \ j^* = \text{min}\left(0,j,m(1)-m(2)+1\right) \]
\[\sum_{i=1}^{j-1} n-M(1) = \left(n-M(1)\right) \left(j-1\right) \]
where \(m\) = min and \(M\) = max.
| [in] | ij | 2D coords. (i,j) |
| [in] | sym | whether symmetric |
| [in] | n | max of 2-D coords. |
| [in] | lim_n | min. and max. of 2-D coords. |
Definition at line 2555 of file num_utilities.f90.
| subroutine, public num_utilities::calc_aux_utilities | ( | integer, intent(in), optional | n_mod | ) |
Initialize utilities for fast future reference, depending on program style.
Utilities initialized:
If Fourier modes are also initialized, the quantity n_mod has to be provided as well.
| [in] | n_mod | n_mod for Fourier modes |
Definition at line 2062 of file num_utilities.f90.
| integer function, public num_utilities::calc_coeff_fin_diff | ( | integer, intent(in) | deriv, |
| integer, intent(in) | nr, | ||
| integer, intent(in) | ind, | ||
| real(dp), dimension(:), intent(inout), allocatable | coeff ) |
Calculate the coefficients for finite differences.
These represent a derivative of degree deriv at an index ind, by the weighted sum of a number nr of points with 1<=ind<=nr.
They are found by solving a Vandermonde system, multiplied by the faculty of the derivative.
The result returned in coeff(1:nr) are used in
\[\frac{d f}{d x} = \sum_{j=1}^n c(j) f(i+j) . \]
where \(i\) = ind, \(n\) = nr and \(c\) = coeff
Examples:
deriv of order ord: deriv of order ord: ord, \(d\) = deriv and \(m\) = nr_2.| [in] | deriv | degree of derivative |
| [in] | nr | number of points |
| [in] | ind | position of derivative |
| [in,out] | coeff | output coefficients |
Definition at line 2448 of file num_utilities.f90.
| integer function, public num_utilities::calc_d2_smooth | ( | integer, intent(in) | fil_n, |
| real(dp), dimension(:), intent(in) | x, | ||
| real(dp), dimension(:), intent(in) | y, | ||
| real(dp), dimension(:), intent(inout) | d2y, | ||
| integer, intent(in) | style ) |
Calculate second derivative with smoothing formula by Holoborodko, [9].
The formula used is
\(s\left(k\right) = \left[\left(2N-10\right)s\left(k+1\right) - \left(N+2k+3\right)s\left(k+2\right)\right]/\left(N-2k-1\right)\)
with \(s\left(k>M\right) = 0\) and \(\left(k=N\right) = 1\).
for a given \(N\).
For style 1, central differences are used:
\(\frac{\partial^2 y}{\partial x^2} \approx \frac{1}{2^{N-3}} \left(\sum_{k=1}^M \alpha_k \left(y_k + y_{-k} \right) - 2 y_0 \sum_{k=1}^M \alpha_k \right)\)
and for style 2, backward differences:
\(\frac{\partial^2 y}{\partial x^2} \approx \frac{1}{2^{N-3}} \left(\sum_{k=1}^M \beta_k \left(y_{-M+k} + y_{-M-k} \right) - 2 y_{-M} \sum_{k=1}^M \beta_k \right)\)
with \(\alpha_k\) and \(\beta_k\) given by
\(\begin{aligned} \alpha_k &= \frac{4k^2 s_k}{\left(x_k-x_{-k}\right)^2} \\ \beta_k &= \frac{4k^2 s_k}{\left(x_{-M+k}-x_{-M-k}\right)^2} \end{aligned}\).
| [in] | fil_n | filter length (must be odd) |
| [in] | x | input abscissa |
| [in] | y | input ordinate |
| [in,out] | d2y | second derivative |
Definition at line 2828 of file num_utilities.f90.
| integer function, public num_utilities::calc_ext_var | ( | real(dp), intent(inout) | ext_var, |
| real(dp), dimension(:), intent(in) | var, | ||
| real(dp), dimension(:), intent(in) | var_points, | ||
| real(dp), intent(in) | ext_point, | ||
| integer, intent(in), optional | deriv_in ) |
Extrapolates a function.
This is done using linear or quadratic interpolation, depending on the number of points and values given.
The data should be given sorted, in ascending order, without duplicate points in var_points.
It uses the following solution:
\[\overline{\text{A}} \ \vec{b} = \vec{c}\]
which is written out to
\[\left(\begin{array}{cccc} x_1^0 & x_1^1 & x_1^2 & \ldots \\ x_2^0 & x_2^1 & x_2^2 & \ldots \\ x_3^0 & x_3^1 & x_3^2 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) \left(\begin{array}{c} a_0 \\ a_1 \\ a_2 \\ \vdots \end{array}\right) = \left(\begin{array}{c} v_0 \\ v_1 \\ v_2 \\ \vdots \end{array}\right) \]
where \(v_{i-1}\) = var_i.
This is solved for the coefficients of the interpolating polynomial
ext_var = a_0 + a_1 x + a_2 x^2 + ...
There is an optional flag to look for the \(k^\text{th}\) derivative of the function instead of the function itself, where \(k\) should be lower than the degree of the polynomial.
| [in,out] | ext_var | output |
| [in] | var | abscissa |
| [in] | var_points | ordinate |
| [in] | ext_point | point to which extrapolate |
| [in] | deriv_in | specifies an optional derivative |
Definition at line 2323 of file num_utilities.f90.
| integer function, public num_utilities::check_deriv | ( | integer, dimension(3), intent(in) | deriv, |
| integer, intent(in) | max_deriv, | ||
| character(len=*), intent(in) | sr_name ) |
checks whether the derivatives requested for a certain subroutine are valid
| [in] | deriv | derivative to be checked |
| [in] | max_deriv | maximum derivative allowed |
| [in] | sr_name | name of subroutine |
Definition at line 2259 of file num_utilities.f90.
| integer function, dimension(:,:), allocatable, public num_utilities::derivs | ( | integer, intent(in) | order, |
| integer, intent(in), optional | dims ) |
Returns derivatives of certain order.
| [in] | order | order of derivative |
| [in] | dims | nr. of dimensions |
Definition at line 2245 of file num_utilities.f90.
| recursive integer function, public num_utilities::fac | ( | integer, intent(in) | n | ) |
Calculate factorial.
| [in] | n | input |
Definition at line 2608 of file num_utilities.f90.
| recursive integer function, public num_utilities::gcd | ( | integer, intent(in) | u, |
| integer, intent(in) | v ) |
Returns least denominator using the GCD.
| [in] | u | input |
| [in] | v | input |
Definition at line 2668 of file num_utilities.f90.
| integer function, public num_utilities::is_sym | ( | integer, intent(in) | n, |
| integer, intent(in) | nn, | ||
| logical, intent(inout) | sym ) |
Determines whether a matrix making use of the storage convention in eq_vars.eq_2_type is symmetric or not.
| [in] | n | size of matrix |
| [in] | nn | number of elements in matrix |
| [in,out] | sym | output |
Definition at line 2625 of file num_utilities.f90.
| recursive integer function, public num_utilities::lcm | ( | integer, intent(in) | u, |
| integer, intent(in) | v ) |
Returns common multiple using the Euclid's algorithm.
| [in] | u | input |
| [in] | v | input |
Definition at line 2656 of file num_utilities.f90.
| subroutine, public num_utilities::shift_f | ( | integer, dimension(2), intent(in) | al, |
| integer, dimension(2), intent(in) | bl, | ||
| integer, dimension(2), intent(in) | cl, | ||
| real(dp), dimension(al(1):al(2),2), intent(in) | a, | ||
| real(dp), dimension(bl(1):bl(2),2), intent(in) | b, | ||
| real(dp), dimension(cl(1):cl(2),2), intent(inout) | c ) |
Calculate multiplication through shifting of fourier modes A and B into C.
This works by calculating \(\left[\sum_A \left(\alpha_A \cos m_A \theta + \beta_A \sin m_A \theta \right)\right] \left[\sum_B \left(\alpha_B \cos m_B \theta + \beta_B \sin m_B \theta \right)\right] = \left[\sum_C \left(\alpha_C \cos m_C \theta + \beta_C \sin m_C \theta \right)\right]\)
where the \(\alpha\) and \(\beta\) factors are provide in A, B and C.
This then boils down to finding the four combinations of cosines and sines.
| [in] | al | limits on A mode numbers |
| [in] | bl | limits on B mode numbers |
| [in] | cl | limits on C mode numbers |
| [in] | a | input |
| [in] | b | input |
| [in,out] | c | result |
Definition at line 2695 of file num_utilities.f90.
| subroutine, public num_utilities::solve_vand | ( | integer, intent(in) | n, |
| real(dp), dimension(n), intent(in) | a, | ||
| real(dp), dimension(n), intent(in) | b, | ||
| real(dp), dimension(n), intent(inout) | x, | ||
| logical, intent(in), optional | transp ) |
Solve a Vandermonde system \(\overline{\text{A}} \ \vec{X} = \vec{B}\).
The Vandermonde matrix has the form
\[A = \left(\begin{array}{ccccc} 1 & a_1 & a_1^2 & \cdots & a_1^n \\ 1 & a_2 & a_2^2 & \cdots & a_2^n \\ 1 & a_3 & a_3^2 & \cdots & a_3^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & a_n & a_n^2 & \cdots & a_n^n \end{array}\right)\]
dvand and pvand in https://people.sc.fsu.edu/~jburkardt/f_src/vandermonde/vandermonde.html based on the Björk-Pereyra Algorithm [2]. | [in] | n | matrix size |
| [in] | a | parameters of Vandermonde matrix |
| [in] | b | right-hand side |
| [in,out] | x | solution |
| [in] | transp | transposed Vandermonde matrix |
Definition at line 2750 of file num_utilities.f90.
| integer, dimension(:,:,:), allocatable, public num_utilities::d |
1-D array indices of derivatives
Definition at line 23 of file num_utilities.f90.
| logical, public num_utilities::debug_calc_coeff_fin_diff = .false. |
plot debug information for calc_coeff_fin_diff()
Definition at line 28 of file num_utilities.f90.
| logical, public num_utilities::debug_con2dis_reg = .false. |
plot debug information for con2dis_reg()
Definition at line 27 of file num_utilities.f90.
| integer, dimension(:,:), allocatable, public num_utilities::f |
1-D array indices of Fourier mode combination indices
Definition at line 25 of file num_utilities.f90.
| integer, dimension(:,:), allocatable, public num_utilities::m |
1-D array indices of metric indices
Definition at line 24 of file num_utilities.f90.