Generalized Pseudospectral Method Module written by Son, Sang-Kil in Dec. 2003 - Feb. 2004 rewritten by Son, Sang-Kil in Feb. 2005 - May 2005, Aug. 2005 revised by Son, Sang-Kil in Sep. 2005
This module makes 1D grid points with the following methods: - Fourier Grid Hamiltonian method (FGH) : [1...N] - Legendre-Gauss collocation method (LG) : [1...N] - Legendre-Gauss-Lobatto collocation method (LGL) : [1...N] - extended Legendre-Gauss-Lobatto collocation method (ext-LGL) : [0...N+1] * Legendre-Gauss-Radau collocation method (LGR) : [1...N]
Also, we can choose any mapping functions like: - linear mapping: r = a x + b - algebraic (reciprocal) mapping: r = L (1+x)/(1-x) - inverse-reciprocal mapping: r = Rmax - L(1-x)/(1+x+alpha) - arbitrary mapping: it should be defined by mapping_f and mapping_fp
USAGE: type(Grid) :: G G = generate_Grid( "Legendre-Gauss-Lobatto", N, Rmin, Rmax, L ) -- or -- G = generate_Grid( "Legendre-Gauss", N, Rmin, Rmax, L, "algebraic" ) -- or -- if you want an arbitrary mapping function, G = generate_Grid( "LGL", N ) / set r(:) and rp(:) as you want / call prepare_GPS_matrix( G )
USEFUL SUBROUTINES - SR apply_mapping( G, Rmin, Rmax, L, mapping ): reset a mapping - SR prepare_GPS_matrix( G ): re-construct D1 & D2 after changing a mapping - SR impose_LGL_boundary( G, boundary="both" ): impose Neumann B.C. - F GPS_interpolate( G, f(:), r ): get a function value at r - F GPS_grid_type( G ): return grid type - F GPS_mapping_r( G, x ): return r corresponding to x - F GPS_mapping_x( G, r ): return x corresponding to r - F GPS_mapping_rp( G, x ): return r' corresponding to x
Ref) Yao & Chu, CPL 204, 381 (1993): for GPS Tong & Chu, CP 217, 119 (1997) Telnov & Chu, PRA 59, 2864 (1999) Marston & Balint-Kurti, JCP 91, 3571 (1989) Canuto, et al., Spectral Methods in Fluid Dynamics (1988) Fornberg, A Practical Guide to Pseudospectral Methods (1996)
NOTE: 1) All quantities include r'_i for convenience. "D1" is defined by d1_ij / sqrt( r'_i * r'_j ). "D2" is defined by d2_ij / ( r'_i * r'_j ). "weight" is defined by w_i * r'_i. originally, integral := SUM( f(:) * r'(:) * w(:) ) but here, integral := SUM( f(:) * weight(:) ) And "Deriv" is defined by (original_d1)_ij / r'_i. originally, f'(i) := SUM( orig_D1(i,:) * f(:) ) / r'(i) but here, f'(i) := SUM( Deriv(i,:) * f(:) )
2) This module is designed to generate N points in (-1,1) for any method, i.e. x1, x2, ..., xN in (-1,1) range. -1 and 1 are excluded. For Legendre-Gauss-Lobatto, N points come from zeros of P'_N+1(x)=0. For Legendre-Gauss, N points come from zeros of P'_N(x)=0. Thus, some expressions should be modified from Canuto's book. For Legendre-Gauss-Lobatto, N in Canuto's notation --> N+1 For Legendre-Gauss, N+1 in Canuto's notation --> N
Grid for 1D N : # of grid points, 1...N Rmin, Rmax : the range in real space x(:) : grid points in the computational space (-1:1) r(:) : grid points in the real space (Rmin:Rmax), not [ ] rp(:) : dr/dx = 1 / (dx/dr) weight(:) : weights for integration D1(:,:) : first derivatives matrix D2(:,:) : second derivatives matrix Deriv(:,:) : original 1st derivatives matrix divided by r'(i) we can compute f' directly by f'(i) = sum( Deriv(i,:)*f(:) ) L : a mapping parameter for the reciprocal mapping mapping : mapping type ( 0: arbitrary, 1: linear, 2: algeb., 3: inv.rec. ) grid_type : grid type ( 0: FGH, 1: LG, 2: LGL, 3: ext-LGL )
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| integer, | public | :: | n | ||||
| real(kind=long), | public | :: | rmin | ||||
| real(kind=long), | public | :: | rmax | ||||
| real(kind=long), | public | :: | l | ||||
| real(kind=long), | public, | allocatable | :: | x(:) | |||
| real(kind=long), | public, | allocatable | :: | r(:) | |||
| real(kind=long), | public, | allocatable | :: | rp(:) | |||
| real(kind=long), | public, | allocatable | :: | weight(:) | |||
| real(kind=long), | public, | allocatable | :: | wf_coeff(:) | |||
| real(kind=long), | public, | allocatable | :: | deriv(:,:) | |||
| real(kind=long), | public, | allocatable | :: | d1(:,:) | |||
| real(kind=long), | public, | allocatable | :: | d2(:,:) | |||
| integer, | public | :: | mapping | ||||
| integer, | public | :: | grid_type |
general function to generate grids grid_type= LGL, LG, or FGH mapping_type= linear or algebraic(reciprocal) This function does 1) set parameters of N and grid_type 2) allocate memory for r, rp, D1, D2, and so on 2) make computational grids 3) make real grids by means of a mapping function 4) construct D1 and D2
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| character(len=*) | :: | grid_type | ||||
| integer, | intent(in) | :: | n | |||
| type(grid), | intent(out) | :: | g | |||
| real(kind=long), | intent(in), | optional | :: | rmin | ||
| real(kind=long), | intent(in), | optional | :: | rmax | ||
| real(kind=long), | intent(in), | optional | :: | l | ||
| character(len=*), | intent(in), | optional | :: | mapping |
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(grid), | intent(inout) | :: | g |