Neko  0.9.0
A portable framework for high-order spectral element flow simulations
fast3d.f90
Go to the documentation of this file.
1 ! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
2 !
3 ! The UChicago Argonne, LLC as Operator of Argonne National
4 ! Laboratory holds copyright in the Software. The copyright holder
5 ! reserves all rights except those expressly granted to licensees,
6 ! and U.S. Government license rights.
7 !
8 ! Redistribution and use in source and binary forms, with or without
9 ! modification, are permitted provided that the following conditions
10 ! are met:
11 !
12 ! 1. Redistributions of source code must retain the above copyright
13 ! notice, this list of conditions and the disclaimer below.
14 !
15 ! 2. Redistributions in binary form must reproduce the above copyright
16 ! notice, this list of conditions and the disclaimer (as noted below)
17 ! in the documentation and/or other materials provided with the
18 ! distribution.
19 !
20 ! 3. Neither the name of ANL nor the names of its contributors
21 ! may be used to endorse or promote products derived from this software
22 ! without specific prior written permission.
23 !
24 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
27 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
28 ! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
29 ! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30 ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
31 ! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32 ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33 ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35 ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 !
37 ! Additional BSD Notice
38 ! ---------------------
39 ! 1. This notice is required to be provided under our contract with
40 ! the U.S. Department of Energy (DOE). This work was produced at
41 ! Argonne National Laboratory under Contract
42 ! No. DE-AC02-06CH11357 with the DOE.
43 !
44 ! 2. Neither the United States Government nor UCHICAGO ARGONNE,
45 ! LLC nor any of their employees, makes any warranty,
46 ! express or implied, or assumes any liability or responsibility for the
47 ! accuracy, completeness, or usefulness of any information, apparatus,
48 ! product, or process disclosed, or represents that its use would not
49 ! infringe privately-owned rights.
50 !
51 ! 3. Also, reference herein to any specific commercial products, process,
52 ! or services by trade name, trademark, manufacturer or otherwise does
53 ! not necessarily constitute or imply its endorsement, recommendation,
54 ! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
55 ! The views and opinions of authors expressed
56 ! herein do not necessarily state or reflect those of the United States
57 ! Government or UCHICAGO ARGONNE, LLC, and shall
58 ! not be used for advertising or product endorsement purposes.
59 !
61 module fast3d
62  use num_types, only : rp
63  use speclib
64  use math, only : rzero
65  implicit none
66  private
67 
69 
70 contains
71 
104  subroutine fd_weights_full(xi, x, n, m, c)
105  integer, intent(in) :: n
106  integer, intent(in) :: m
107  real(kind=rp), intent(in) :: x(0:n)
108  real(kind=rp), intent(out) :: c(0:n,0:m)
109  real(kind=rp), intent(in) :: xi
110  real(kind=xp) :: c1, c2, c3, c4, c5
111  integer :: i, j, k, mn
112 
113  c1 = 1d0
114  c4 = x(0) - xi
115 
116  do k = 0, m
117  do j = 0, n
118  c(j,k) = 0d0
119  end do
120  end do
121 
122  c(0,0) = 1d0
123 
124  do i = 1, n
125  mn = min(i,m)
126  c2 = 1d0
127  c5 = c4
128  c4 = x(i) - xi
129  do j = 0, i - 1
130  c3 = x(i) - x(j)
131  c2 = c2 * c3
132  do k = mn, 1, -1
133  c(i,k) = c1 * (k * c(i-1,k-1) - c5 * c(i-1,k)) / c2
134  end do
135  c(i,0) = -c1 * c5 * c(i-1,0) / c2
136  do k = mn, 1, -1
137  c(j,k) = (c4 * c(j,k) - k * c(j,k-1)) / c3
138  end do
139  c(j,0) = c4 * c(j,0) / c3
140  end do
141  c1 = c2
142  end do
143 
144  end subroutine fd_weights_full
145 
146 
167  subroutine semhat(a, b, c, d, z, dgll, jgll, bgl, zgl, dgl, jgl, n, w)
168  integer, intent(in) :: n
169  real(kind=rp), intent(inout) :: a(0:n,0:n)
170  real(kind=rp), intent(inout) :: b(0:n)
171  real(kind=rp), intent(inout) :: c(0:n,0:n)
172  real(kind=rp), intent(inout) :: d(0:n,0:n)
173  real(kind=rp), intent(inout) :: z(0:n)
174  real(kind=rp), intent(inout) :: dgll(0:n,1:n-1),jgll(0:n,1:n-1)
175  real(kind=rp), intent(inout) :: bgl(1:n-1)
176  real(kind=rp), intent(inout) :: zgl(1:n-1)
177  real(kind=rp), intent(inout) :: dgl(1:n-1,0:n)
178  real(kind=rp), intent(inout) :: jgl(1:n-1,0:n)
179  real(kind=rp), intent(inout) :: w(0:2*n+1)
180  integer :: np, nm, n2, i, j, k
181  np = n+1
182  nm = n-1
183  n2 = n-2
184  call zwgll(z, b, np)
185  do i = 0,n
186  call fd_weights_full(z(i), z, n, 1, w)
187  do j = 0,n
188  d(i,j) = w(j+np) ! Derivative matrix
189  end do
190  end do
191 
192  if (n.eq.1) return ! No interpolation for n=1
193 
194  do i = 0,n
195  call fd_weights_full(z(i), z(1), n2, 1, w(1))
196  do j = 1,nm
197  jgll(i,j) = w(j ) ! Interpolation matrix
198  dgll(i,j) = w(j+nm) ! Derivative matrix
199  end do
200  end do
201  call rzero(a, np*np)
202  do j = 0,n
203  do i = 0,n
204  do k = 0,n
205  a(i,j) = a(i,j) + d(k,i)*b(k)*d(k,j)
206  end do
207  c(i,j) = b(i)*d(i,j)
208  end do
209  end do
210  call zwgl(zgl, bgl, nm)
211  do i = 1,n-1
212  call fd_weights_full(zgl(i), z, n, 1, w)
213  do j = 0,n
214  jgl(i,j) = w(j ) ! Interpolation matrix
215  dgl(i,j) = w(j+np) ! Derivative matrix
216  end do
217  end do
218  end subroutine semhat
219 
242  subroutine setup_intp(jh, jht, z_to, z_from, n_to, n_from, derivative)
243  implicit none
244  integer, intent(in) :: n_to, n_from, derivative
245  real(kind=rp), intent(inout) :: jh(n_to, n_from), jht(n_from, n_to)
246  real(kind=rp), intent(inout) :: z_to(n_to), z_from(n_from)
247  real(kind=rp) :: w(n_from, 0:derivative)
248  integer :: i, j
249  do i = 1, n_to
250  ! This will assign w the weights for interpolating to point
251  ! zf(i) values at points zc
252  call fd_weights_full(z_to(i), z_from, n_from-1, derivative, w)
253 
254  ! store each of the weights in the corresponding row/column
255  do j = 1, n_from
256  jh(i,j) = w(j, derivative)
257  jht(j,i) = w(j, derivative)
258  end do
259  end do
260  end subroutine setup_intp
261 end module fast3d
Implements the derivative_t type.
Definition: derivative.f90:36
Fast diagonalization methods from NEKTON.
Definition: fast3d.f90:61
subroutine, public semhat(a, b, c, d, z, dgll, jgll, bgl, zgl, dgl, jgl, n, w)
Generate matrices for single element, 1D operators: a = Laplacian b = diagonal mass matrix c = convec...
Definition: fast3d.f90:168
subroutine, public fd_weights_full(xi, x, n, m, c)
Compute finite-difference stencil weights for evaluating derivatives up to order at a point.
Definition: fast3d.f90:105
subroutine, public setup_intp(jh, jht, z_to, z_from, n_to, n_from, derivative)
Compute interpolation weights for points z_to using values at points z_from.
Definition: fast3d.f90:243
Definition: math.f90:60
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
LIBRARY ROUTINES FOR SPECTRAL METHODS.
Definition: speclib.f90:148
subroutine zwgll(Z, W, NP)
Definition: speclib.f90:169
subroutine dgll(D, DT, Z, NZ, NZD)
Definition: speclib.f90:865
subroutine zwgl(Z, W, NP)
Generate NP Gauss Legendre points Z and weights W associated with Jacobi polynomial ....
Definition: speclib.f90:161