Neko 1.99.5
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
speclib.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!
60!==============================================================================
61!
75!
76!------------------------------------------------------------------------------
77!
78! ABBRIVIATIONS:
79!
80! M - Set of mesh points
81! Z - Set of collocation/quadrature points
82! W - Set of quadrature weights
83! H - Lagrangian interpolant
84! D - Derivative operator
85! I - Interpolation operator
86! GL - Gauss Legendre
87! GLL - Gauss-Lobatto Legendre
88! GJ - Gauss Jacobi
89! GLJ - Gauss-Lobatto Jacobi
90!
91!
92! MAIN ROUTINES:
93!
94! Points and weights:
95!
96! ZWGL Compute Gauss Legendre points and weights
97! ZWGLL Compute Gauss-Lobatto Legendre points and weights
98! ZWGJ Compute Gauss Jacobi points and weights (general)
99! ZWGLJ Compute Gauss-Lobatto Jacobi points and weights (general)
100!
101! Lagrangian interpolants:
102!
103! HGL Compute Gauss Legendre Lagrangian interpolant
104! HGLL Compute Gauss-Lobatto Legendre Lagrangian interpolant
105! HGJ Compute Gauss Jacobi Lagrangian interpolant (general)
106! HGLJ Compute Gauss-Lobatto Jacobi Lagrangian interpolant (general)
107!
108! Derivative operators:
109!
110! DGLL Compute Gauss-Lobatto Legendre derivative matrix
111! DGLLGL Compute derivative matrix for a staggered mesh (GLL->GL)
112! DGJ Compute Gauss Jacobi derivative matrix (general)
113! DGLJ Compute Gauss-Lobatto Jacobi derivative matrix (general)
114! DGLJGJ Compute derivative matrix for a staggered mesh
115! (GLJ->GJ) (general)
116!
117! Interpolation operators:
118!
119! IGLM Compute interpolation operator GL -> M
120! IGLLM Compute interpolation operator GLL -> M
121! IGJM Compute interpolation operator GJ -> M (general)
122! IGLJM Compute interpolation operator GLJ -> M (general)
123!
124! Other:
125!
126! PNLEG Compute Legendre polynomial of degree N
127! legendre_poly Compute Legendre polynomial of degree 0-N
128! PNDLEG Compute derivative of Legendre polynomial of degree N
129!
130! Comments:
131!
132! Note that many of the above routines exist in both single and
133! double precision. If the name of the single precision routine is
134! SUB, the double precision version is called SUBD. In most cases
135! all the "low-level" arithmetic is done in double precision, even
136! for the single precsion versions.
137!
138! Useful references:
139!
140! [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society,
141! Providence, Rhode Island, 1939.
142! [2] Abramowitz & Stegun: Handbook of Mathematical Functions,
143! Dover, New York, 1972.
144! [3] Canuto, Hussaini, Quarteroni & Zang: Spectral Methods in Fluid
145! Dynamics, Springer-Verlag, 1988.
146!
147!
148!==============================================================================
149module speclib
150 use num_types, only : rp, xp
151 use math, only: abscmp
152 use utils, only : neko_error
153 use, intrinsic :: iso_fortran_env, only : stderr => error_unit
154 implicit none
155
156contains
157
164 subroutine zwgl(Z, W, NP)
165 integer, intent(in) :: NP
166 real(kind=rp), intent(inout) :: z(np), w(np)
167 real(kind=rp) alpha, beta
168 alpha = 0.0_rp
169 beta = 0.0_rp
170 call zwgj(z, w, np, alpha, beta)
171 end subroutine zwgl
172
173
179 subroutine zwgll(Z, W, NP)
180 integer, intent(in) :: NP
181 real(kind=rp), intent(inout) :: z(np), w(np)
182 real(kind=rp) alpha, beta
183 alpha = 0.0_rp
184 beta = 0.0_rp
185 call zwglj(z, w, np, alpha, beta)
186 end subroutine zwgll
187
192 subroutine zwgj(Z, W, NP, ALPHA, BETA)
193 integer, intent(in) :: NP
194 real(kind=rp), intent(inout) :: z(np), w(np)
195 real(kind=rp), intent(in) :: alpha, beta
196
197 integer, parameter :: NMAX = 84
198 integer, parameter :: NZD = nmax
199
200 real(kind=xp) zd(nzd), wd(nzd), alphad, betad
201 integer :: I, NPMAX
202
203 npmax = nzd
204 if (np .gt. npmax) then
205 write (stderr, *) 'Too large polynomial degree in ZWGJ'
206 write (stderr, *) 'Maximum polynomial degree is', nmax
207 write (stderr, *) 'Here NP=', np
208 call neko_error
209 end if
210
211 alphad = real(alpha, kind=xp)
212 betad = real(beta, kind=xp)
213 call zwgjd(zd, wd, np, alphad, betad)
214 do i = 1, np
215 z(i) = real(zd(i), kind=rp)
216 w(i) = real(wd(i), kind=rp)
217 end do
218 end subroutine zwgj
219
224 subroutine zwgjd(Z, W, NP, ALPHA, BETA)
225 integer, intent(in) :: NP
226 real(kind=xp), intent(inout) :: z(np), w(np)
227 real(kind=xp), intent(in) :: alpha, beta
228
229 real(kind=xp) :: dn, apb
230 real(kind=xp) :: fac1, fac2, fac3, fnorm
231 real(kind=xp) :: rcoef, p, pd, pm1, pdm1, pm2, pdm2
232 real(kind=xp) :: dnp1, dnp2
233 integer :: N, NP1, NP2, I
234
235 n = np - 1
236 dn = real(n, kind=xp)
237
238 apb = alpha + beta
239
240 if (np .le. 0) then
241 write (stderr, *) 'ZWGJD: Minimum number of Gauss points is 1', np
242 call neko_error
243 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
244 call neko_error('ZWGJD: Alpha and Beta must be greater than -1')
245 end if
246
247 if (np .eq. 1) then
248 z(1) = (beta - alpha) / (apb + 2.0_xp)
249 w(1) = gammaf(alpha + 1.0_xp) * gammaf(beta + 1.0_xp) / &
250 gammaf(apb + 2.0_xp) * 2.0_xp**(apb + 1.0_xp)
251 return
252 end if
253
254 call jacg(z, np, alpha, beta)
255
256 np1 = n + 1
257 np2 = n + 2
258 dnp1 = real(np1, kind=xp)
259 dnp2 = real(np2, kind=xp)
260 fac1 = dnp1 + alpha + beta + 1.0_xp
261 fac2 = fac1 + dnp1
262 fac3 = fac2 + 1.0_xp
263 fnorm = pnormj(np1, alpha, beta)
264 rcoef = (fnorm*fac2*fac3) / (2.0_xp*fac1*dnp2)
265 do i = 1, np
266 call jacobf(p, pd, pm1, pdm1, pm2, pdm2, np2, alpha, beta, z(i))
267 w(i) = -rcoef/(p*pdm1)
268 end do
269 end subroutine zwgjd
270
275 subroutine zwglj(Z, W, NP, ALPHA, BETA)
276 integer, intent(in) :: NP
277 real(kind=rp), intent(inout) :: z(np), w(np)
278 real(kind=rp), intent(in) :: alpha, beta
279
280 integer, parameter :: NMAX = 84
281 integer, parameter :: NZD = nmax
282
283 real(kind=xp) zd(nzd), wd(nzd), alphad, betad
284 integer :: I, NPMAX
285
286 npmax = nzd
287 if (np .gt. npmax) then
288 write (stderr, *) 'Too large polynomial degree in ZWGLJ'
289 write (stderr, *) 'Maximum polynomial degree is', nmax
290 write (stderr, *) 'Here NP=', np
291 call neko_error
292 end if
293 alphad = real(alpha, kind=xp)
294 betad = real(beta, kind=xp)
295 call zwgljd(zd, wd, np, alphad, betad)
296 do i = 1, np
297 z(i) = real(zd(i), kind=rp)
298 w(i) = real(wd(i), kind=rp)
299 end do
300 end subroutine zwglj
301
302
307 subroutine zwgljd(Z, W, NP, ALPHA, BETA)
308
309 integer, intent(in) :: NP
310 real(kind=xp), intent(inout) :: z(np), w(np)
311 real(kind=xp), intent(in) :: alpha, beta
312
313 real(kind=xp) :: alpg, betg
314 real(kind=xp) :: p, pd, pm1, pdm1, pm2, pdm2
315 integer :: N, NM1, I
316
317 n = np - 1
318 nm1 = n - 1
319
320 if (np .le. 1) then
321 write (stderr, *) 'ZWGLJD: Minimum number of Gauss-Lobatto points is 2'
322 write (stderr, *) 'ZWGLJD: alpha, beta:', alpha, beta, np
323 call neko_error
324 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
325 call neko_error('ZWGLJD: Alpha and Beta must be greater than -1')
326 end if
327
328 if (nm1 .gt. 0) then
329 alpg = alpha + 1.0_xp
330 betg = beta + 1.0_xp
331 call zwgjd(z(2), w(2), nm1, alpg, betg)
332 end if
333
334 z(1) = -1.0_xp
335 z(np) = 1.0_xp
336 do i = 2, np - 1
337 w(i) = w(i) / (1.0_xp-z(i)**2)
338 end do
339 call jacobf(p, pd, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(1))
340 w(1) = endw1(n, alpha, beta) / (2.0_xp*pd)
341 call jacobf(p, pd, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(np))
342 w(np) = endw2(n, alpha, beta) / (2.0_xp*pd)
343
344 end subroutine zwgljd
345
347 real(kind=xp) function endw1(N, ALPHA, BETA)
348
349 real(kind=xp), intent(in) :: alpha, beta
350 integer, intent(in) :: n
351
352 real(kind=xp) :: apb
353 real(kind=xp) :: f1, f2, f3, fint1, fint2
354 real(kind=xp) :: a1, a2, a3, di, abn, abnn
355 integer :: i
356
357 if (n .eq. 0) then
358 endw1 = 0.0_xp
359 return
360 end if
361
362 apb = alpha + beta
363 f1 = gammaf(alpha + 2.0_xp)*gammaf(beta + 1.0_xp) / gammaf(apb + 3.0_xp)
364 f1 = f1*(apb + 2.0_xp)*2.0_xp**(apb + 2.0_xp)/2.0_xp
365 if (n .eq. 1) then
366 endw1 = f1
367 return
368 end if
369
370 fint1 = gammaf(alpha + 2.0_xp)*gammaf(beta + 1.0_xp) / gammaf(apb + 3.0_xp)
371 fint1 = fint1*2.0_xp**(apb + 2.0_xp)
372 fint2 = gammaf(alpha + 2.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 4.0_xp)
373 fint2 = fint2*2.0_xp**(apb + 3.0_xp)
374 f2 = (-2.0_xp*(beta + 2.0_xp)*fint1 + (apb + 4.0_xp)*fint2) * &
375 (apb + 3.0_xp) / 4.0_xp
376 if (n .eq. 2) then
377 endw1 = f2
378 return
379 end if
380
381 do i = 3, n
382 di = real(i - 1, kind=xp)
383 abn = alpha + beta + di
384 abnn = abn + di
385 a1 = -(2.0_xp*(di + alpha) * (di + beta)) / (abn*abnn*(abnn + 1.0_xp))
386 a2 = (2.0_xp*(alpha - beta)) / (abnn*(abnn + 2.0_xp))
387 a3 = (2.0_xp*(abn + 1.0_xp)) / ((abnn + 2.0_xp) * (abnn + 1.0_xp))
388 f3 = -(a2*f2 + a1*f1) / a3
389 f1 = f2
390 f2 = f3
391 end do
392 endw1 = f3
393 end function endw1
394
396 real(kind=xp) function endw2(N, ALPHA, BETA)
397
398 real(kind=xp), intent(in) :: alpha, beta
399 integer, intent(in) :: n
400
401 real(kind=xp) :: apb
402 real(kind=xp) :: f1, f2, f3, fint1, fint2
403 real(kind=xp) :: a1, a2, a3, di, abn, abnn
404 integer :: i
405
406
407 if (n .eq. 0) then
408 endw2 = 0.0_xp
409 return
410 end if
411
412 apb = alpha + beta
413 f1 = gammaf(alpha + 1.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 3.0_xp)
414 f1 = f1*(apb + 2.0_xp)*2.0_xp**(apb + 2.0_xp)/2.0_xp
415 if (n .eq. 1) then
416 endw2 = f1
417 return
418 end if
419
420 fint1 = gammaf(alpha + 1.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 3.0_xp)
421 fint1 = fint1*2.0_xp**(apb + 2.0_xp)
422 fint2 = gammaf(alpha + 2.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 4.0_xp)
423 fint2 = fint2*2.0_xp**(apb + 3.0_xp)
424 f2 = (2.0_xp*(alpha + 2.0_xp)*fint1 - (apb + 4.0_xp)*fint2) * &
425 (apb + 3.0_xp) / 4.0_xp
426 if (n .eq. 2) then
427 endw2 = f2
428 return
429 end if
430
431 do i = 3, n
432 di = ((i-1))
433 abn = alpha + beta + di
434 abnn = abn + di
435 a1 = -(2.0_xp*(di + alpha) * (di + beta)) / (abn*abnn*(abnn + 1.0_xp))
436 a2 = (2.0_xp*(alpha - beta)) / (abnn*(abnn + 2.0_xp))
437 a3 = (2.0_xp*(abn + 1.0_xp)) / ((abnn + 2.0_xp) * (abnn + 1.0_xp))
438 f3 = -(a2*f2 + a1*f1)/a3
439 f1 = f2
440 f2 = f3
441 end do
442 endw2 = f3
443 end function endw2
444
446 real(kind=xp) function gammaf(X)
447 real(kind=xp), intent(in) :: x
448 real(kind=xp), parameter :: pi = 4.0_xp*atan(1.0_xp)
449
450 gammaf = 1.0_xp
451 if (abscmp(x, -0.5_xp)) gammaf = -2.0_xp*sqrt(pi)
452 if (abscmp(x, 0.5_xp)) gammaf = sqrt(pi)
453 if (abscmp(x, 1.0_xp)) gammaf = 1.0_xp
454 if (abscmp(x, 2.0_xp)) gammaf = 1.0_xp
455 if (abscmp(x, 1.5_xp)) gammaf = sqrt(pi) / 2.0_xp
456 if (abscmp(x, 2.5_xp)) gammaf = 1.5_xp * sqrt(pi) / 2.0_xp
457 if (abscmp(x, 3.5_xp)) gammaf = 0.5_xp * (2.5_xp * (1.5_xp * sqrt(pi)))
458 if (abscmp(x, 3.0_xp)) gammaf = 2.0_xp
459 if (abscmp(x, 4.0_xp)) gammaf = 6.0_xp
460 if (abscmp(x, 5.0_xp)) gammaf = 24.0_xp
461 if (abscmp(x, 6.0_xp)) gammaf = 120.0_xp
462 end function gammaf
463
465 real(kind=xp) function pnormj(N, ALPHA, BETA)
466 real(kind=xp), intent(in) :: alpha, beta
467 integer, intent(in) :: n
468
469 real(kind=xp) :: dn, dindx
470 real(kind=xp) :: const, prod, frac
471 integer :: i
472
473 dn = real(n, kind=xp)
474 const = alpha + beta + 1.0_xp
475 if (n .le. 1) then
476 prod = gammaf(dn + alpha)*gammaf(dn + beta)
477 prod = prod / (gammaf(dn)*gammaf(dn + alpha + beta))
478 pnormj = prod * 2.0_xp**const / (2.0_xp*dn + const)
479 return
480 end if
481
482 prod = gammaf(alpha + 1.0_xp)*gammaf(beta + 1.0_xp)
483 prod = prod/(2.0_xp*(1.0_xp + const)*gammaf(const + 1.0_xp))
484 prod = prod*(1.0_xp + alpha) * (2.0_xp + alpha)
485 prod = prod*(1.0_xp + beta) * (2.0_xp + beta)
486 do i = 3, n
487 dindx = real(i, kind=xp)
488 frac = (dindx + alpha) * (dindx + beta) / (dindx*(dindx + alpha + beta))
489 prod = prod*frac
490 end do
491 pnormj = prod*2.0_xp**const / (2.0_xp*dn + const)
492 end function pnormj
493
500 subroutine jacg(XJAC, NP, ALPHA, BETA)
501 integer, intent(in) :: NP
502 real(kind=xp), intent(inout) :: xjac(np)
503 real(kind=xp), intent(in) :: alpha, beta
504
505 integer, parameter :: KSTOP = 10
506 real(kind=rp), parameter :: eps = 1.0e-12_rp
507 real(kind=xp), parameter :: pi = 4.0_xp*atan(1.0_xp)
508
509 real(kind=xp) :: dth, x, x1, x2, xlast, delx, xmin
510 real(kind=xp) :: p, pd, pm1, pdm1, pm2, pdm2
511 real(kind=xp) :: recsum, swap
512 integer :: I, J, K, N, JM, JMIN
513
514 n = np - 1
515 dth = pi / (2.0_xp*real(n, kind=xp) + 2.0_xp)
516 do j = 1, np
517 if (j .eq. 1) then
518 x = cos((2.0_xp*(real(j, kind=xp) - 1.0_xp) + 1.0_xp)*dth)
519 else
520 x1 = cos((2.0_xp*(real(j, kind=xp) - 1.0_xp) + 1.0_xp)*dth)
521 x2 = xlast
522 x = (x1 + x2) / 2.0_xp
523 end if
524
525 do k = 1, kstop
526 call jacobf(p, pd, pm1, pdm1, pm2, pdm2, np, alpha, beta, x)
527 recsum = 0.0_xp
528 jm = j - 1
529 do i = 1, jm
530 recsum = recsum + 1.0_xp / (x-xjac(np - i + 1))
531 end do
532 delx = -p / (pd - recsum*p)
533 x = x + delx
534 if (abs(delx) .lt. eps) exit
535 end do
536
537 xjac(np-j + 1) = x
538 xlast = x
539 end do
540
541 do i = 1, np
542 xmin = 2.
543 do j = i, np
544 if (xjac(j) .lt. xmin) then
545 xmin = xjac(j)
546 jmin = j
547 end if
548 end do
549 if (jmin .ne. i) then
550 swap = xjac(i)
551 xjac(i) = xjac(jmin)
552 xjac(jmin) = swap
553 end if
554 end do
555 end subroutine jacg
556
559 subroutine jacobf(POLY, PDER, POLYM1, PDERM1, POLYM2, PDERM2, N, ALP, BET, X)
560
561 real(kind=xp), intent(inout) :: poly, pder, polym1, pderm1, polym2, pderm2
562 real(kind=xp), intent(in) :: alp, bet, x
563 integer, intent(in) :: N
564
565 real(kind=xp) :: apb, polyl, pderl, polyn, pdern
566 real(kind=xp) :: psave, pdsave
567 real(kind=xp) :: a1, a2, a3, a4, b3
568 real(kind=xp) :: dk
569 integer :: K
570
571 apb = alp + bet
572 poly = 1.0_xp
573 pder = 0.0_xp
574 if (n .eq. 0) return
575
576 polyl = poly
577 pderl = pder
578 poly = (alp - bet + (apb + 2.0_xp)*x) / 2.0_xp
579 pder = (apb + 2.0_xp) / 2.0_xp
580 if (n .eq. 1) return
581
582 do k = 2, n
583 dk = real(k, kind=xp)
584 a1 = 2.0_xp*dk*(dk + apb) * (2.0_xp*dk + apb - 2.0_xp)
585 a2 = (2.0_xp*dk + apb - 1.0_xp) * (alp**2 - bet**2)
586 b3 = (2.0_xp*dk + apb - 2.0_xp)
587 a3 = b3*(b3 + 1.0_xp) * (b3 + 2.0_xp)
588 a4 = 2.0_xp*(dk + alp - 1.0_xp) * (dk + bet - 1.0_xp) * (2.0_xp*dk + apb)
589 polyn = ((a2 + a3*x)*poly - a4*polyl) / a1
590 pdern = ((a2 + a3*x)*pder - a4*pderl + a3*poly) / a1
591 psave = polyl
592 pdsave = pderl
593 polyl = poly
594 poly = polyn
595 pderl = pder
596 pder = pdern
597 end do
598 polym1 = polyl
599 pderm1 = pderl
600 polym2 = psave
601 pderm2 = pdsave
602 end subroutine jacobf
603
607 real(kind=xp) function hgj(II, Z, ZGJ, NP, ALPHA, BETA)
608 integer, intent(in) :: np, ii
609 real(kind=xp), intent(in) :: z, zgj(np), alpha, beta
610
611 integer, parameter :: nmax = 84
612 integer, parameter :: nzd = nmax
613
614 real(kind=xp) zd, zgjd(nzd)
615 integer :: i, npmax
616
617 npmax = nzd
618 if (np .gt. npmax) then
619 write (stderr, *) 'Too large polynomial degree in HGJ'
620 write (stderr, *) 'Maximum polynomial degree is', nmax
621 write (stderr, *) 'Here NP=', np
622 call neko_error
623 end if
624
625 zd = z
626 do i = 1, np
627 zgjd(i) = zgj(i)
628 end do
629 hgj = hgjd(ii, zd, zgjd, np, alpha, beta)
630 end function hgj
631
635 real(kind=xp) function hgjd(II, Z, ZGJ, NP, ALPHA, BETA)
636 integer, intent(in) :: np, ii
637 real(kind=xp), intent(in) :: z, zgj(np), alpha, beta
638
639 real(kind=xp) :: eps, zi, dz
640 real(kind=xp) :: pz, pdz, pzi, pdzi, pm1, pdm1, pm2, pdm2
641
642 eps = 1.0e-5_xp
643 zi = zgj(ii)
644 dz = z - zi
645 if (abs(dz) .lt. eps) then
646 hgjd = 1.0_xp
647 return
648 end if
649 call jacobf(pzi, pdzi, pm1, pdm1, pm2, pdm2, np, alpha, beta, zi)
650 call jacobf(pz, pdz, pm1, pdm1, pm2, pdm2, np, alpha, beta, z)
651 hgjd = pz / (pdzi*(z-zi))
652 end function hgjd
653
657 real(kind=xp) function hglj(II, Z, ZGLJ, NP, ALPHA, BETA)
658 integer, intent(in) :: np, ii
659 real(kind=xp), intent(in) :: z, zglj(np), alpha, beta
660
661 integer, parameter :: nmax = 84
662 integer, parameter :: nzd = nmax
663
664 real(kind=xp) zd, zgljd(nzd)
665 integer :: i, npmax
666
667 npmax = nzd
668 if (np .gt. npmax) then
669 write (stderr, *) 'Too large polynomial degree in HGLJ'
670 write (stderr, *) 'Maximum polynomial degree is', nmax
671 write (stderr, *) 'Here NP=', np
672 call neko_error
673 end if
674 zd = z
675 do i = 1, np
676 zgljd(i) = zglj(i)
677 end do
678 hglj = hgljd(ii, zd, zgljd, np, alpha, beta)
679 end function hglj
680
684 real(kind=xp) function hgljd(I, Z, ZGLJ, NP, ALPHA, BETA)
685 integer, intent(in) :: np, i
686 real(kind=xp), intent(in) :: z, zglj(np), alpha, beta
687
688 real(kind=xp) :: eps, zi, dz, dn
689 real(kind=xp) :: p, pd, pi, pdi, pm1, pdm1, pm2, pdm2
690 real(kind=xp) :: eigval, const
691 integer :: n
692
693 eps = 1.0e-5_xp
694 zi = zglj(i)
695 dz = z-zi
696 if (abs(dz) .lt. eps) then
697 hgljd = 1.0_xp
698 return
699 end if
700
701 n = np - 1
702 dn = real(n, kind=xp)
703 eigval = -dn*(dn + alpha + beta + 1.0_xp)
704 call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, n, alpha, beta, zi)
705 const = eigval*pi + alpha*(1.0_xp + zi)*pdi - beta*(1.0_xp - zi)*pdi
706 call jacobf(p, pd, pm1, pdm1, pm2, pdm2, n, alpha, beta, z)
707 hgljd = (1.0_xp - z**2)*pd / (const*(z - zi))
708 end function hgljd
709
715 subroutine dgj(D, DT, Z, NZ, NZD, ALPHA, BETA)
716 integer, intent(in) :: NZ, NZD
717 real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
718 real(kind=xp), intent(in) :: z(nz), alpha, beta
719
720 integer, parameter :: NMAX = 84
721 integer, parameter :: NZDD = nmax
722
723 real(kind=xp) :: dd(nzdd, nzdd), dtd(nzdd, nzdd), zd(nzdd)
724 integer :: I, J
725
726 if (nz .le. 0) then
727 call neko_error('DGJ: Minimum number of Gauss points is 1')
728 else if (nz .gt. nmax) then
729 write (stderr, *) 'Too large polynomial degree in DGJ'
730 write (stderr, *) 'Maximum polynomial degree is', nmax
731 write (stderr, *) 'Here Nz=', nz
732 call neko_error
733 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
734 call neko_error('DGJ: Alpha and Beta must be greater than -1')
735 end if
736
737 do i = 1, nz
738 zd(i) = z(i)
739 end do
740 call dgjd(dd, dtd, zd, nz, nzdd, alpha, beta)
741 do i = 1, nz
742 do j = 1, nz
743 d(i, j) = dd(i, j)
744 dt(i, j) = dtd(i, j)
745 end do
746 end do
747 end subroutine dgj
748
754 subroutine dgjd(D, DT, Z, NZ, NZD, ALPHA, BETA)
755 integer, intent(in) :: NZ, NZD
756 real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
757 real(kind=xp), intent(in) :: z(nz), alpha, beta
758
759 real(kind=xp) :: dn
760 real(kind=xp) :: pdi, pdj, pi, pj, pm1, pdm1, pm2, pdm2
761 integer :: I, J, N
762
763 n = nz - 1
764 dn = real(n, kind=xp)
765
766
767 if (nz .le. 1) then
768 call neko_error('DGJD: Minimum number of Gauss-Lobatto points is 2')
769 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
770 call neko_error('DGJD: Alpha and Beta must be greater than -1')
771 end if
772
773 do i = 1, nz
774 do j = 1, nz
775 call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, nz, alpha, beta, z(i))
776 call jacobf(pj, pdj, pm1, pdm1, pm2, pdm2, nz, alpha, beta, z(j))
777 if (i .ne. j) then
778 d(i, j) = pdi / (pdj*(z(i) - z(j)))
779 else
780 d(i, j) = ((alpha + beta + 2.0_xp)*z(i) + alpha - beta) / &
781 (2.0_xp*(1.0_xp - z(i)**2))
782 end if
783 dt(j, i) = d(i, j)
784 end do
785 end do
786 end subroutine dgjd
787
793 subroutine dglj(D, DT, Z, NZ, NZD, ALPHA, BETA)
794 integer, parameter :: NMAX = 84
795 integer, parameter :: NZDD = nmax
796 integer, intent(in) :: NZ, NZD
797 real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
798 real(kind=xp), intent(in) :: z(nz), alpha, beta
799
800 real(kind=xp) :: dd(nzdd, nzdd), dtd(nzdd, nzdd), zd(nzdd)
801 integer :: I, J
802
803 if (nz .le. 1) then
804 call neko_error('DGLJ: Minimum number of Gauss-Lobatto points is 2')
805 else if (nz .gt. nmax) then
806 write (stderr, *) 'Too large polynomial degree in DGLJ'
807 write (stderr, *) 'Maximum polynomial degree is', nmax
808 write (stderr, *) 'Here NZ=', nz
809 call neko_error
810 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
811 call neko_error('DGLJ: Alpha and Beta must be greater than -1')
812 end if
813
814 do i = 1, nz
815 zd(i) = z(i)
816 end do
817 call dgljd(dd, dtd, zd, nz, nzdd, alpha, beta)
818 do i = 1, nz
819 do j = 1, nz
820 d(i, j) = dd(i, j)
821 dt(i, j) = dtd(i, j)
822 end do
823 end do
824 end subroutine dglj
825
826
832 subroutine dgljd(D, DT, Z, NZ, NZD, ALPHA, BETA)
833 integer, intent(in) :: NZ, NZD
834 real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
835 real(kind=xp), intent(in) :: z(nz), alpha, beta
836
837 real(kind=xp) :: dn, eigval
838 real(kind=xp) :: pdi, pdj, pi, pj, pm1, pdm1, pm2, pdm2
839 real(kind=xp) :: ci, cj
840 integer :: I, J, N
841
842 n = nz - 1
843 dn = real(n, kind=xp)
844
845 eigval = -dn*(dn + alpha + beta + 1.0_xp)
846
847 if (nz .le. 1) then
848 call neko_error('DGLJD: Minimum number of Gauss-Lobatto points is 2')
849 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
850 call neko_error('DGLJD: Alpha and Beta must be greater than -1')
851 end if
852
853 do i = 1, nz
854 do j = 1, nz
855 call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(i))
856 call jacobf(pj, pdj, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(j))
857 ci = eigval*pi - (beta*(1.0_xp - z(i)) - alpha*(1.0_xp + z(i)))*pdi
858 cj = eigval*pj - (beta*(1.0_xp - z(j)) - alpha*(1.0_xp + z(j)))*pdj
859
860 ! Todo: This should have some elses in there
861 if (i .ne. j) then
862 d(i, j) = ci / (cj*(z(i) - z(j)))
863 else if (i .eq. 1) then
864 d(i, j) = (eigval + alpha) / (2.0_xp*(beta + 2.0_xp))
865 else if (i .eq. nz) then
866 d(i, j) = -(eigval + beta) / (2.0_xp*(alpha + 2.0_xp))
867 else
868 d(i, j) = (alpha*(1.0_xp + z(i)) - beta*(1.0_xp - z(i))) / &
869 (2.0_xp*(1.0_xp - z(i)**2))
870 end if
871 dt(j, i) = d(i, j)
872 end do
873 end do
874 end subroutine dgljd
875
880 subroutine dgll(D, DT, Z, NZ, NZD)
881
882 integer, intent(in) :: NZ, NZD
883 real(kind=rp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
884 real(kind=rp), intent(in) :: z(nz)
885
886 integer, parameter :: NMAX = 84
887
888 real(kind=xp) :: d0, fn
889 integer :: I, J, N
890
891 n = nz - 1
892 if (nz .gt. nmax) then
893 write (stderr, *) 'Subroutine DGLL'
894 write (stderr, *) 'Maximum polynomial degree =', nmax
895 write (stderr, *) 'Polynomial degree =', nz
896 call neko_error
897 else if (nz .eq. 1) then
898 d(1, 1) = 0.0_rp
899 return
900 end if
901
902 fn = real(n, kind=xp)
903 d0 = fn*(fn + 1.0_xp)/4.0_xp
904 do i = 1, nz
905 do j = 1, nz
906 if (i .ne. j) then
907 d(i, j) = pnleg(real(z(i), xp), n)/ &
908 (pnleg(real(z(j), xp), n) * (z(i) - z(j)))
909 else if (i .eq. 1) then
910 d(i, j) = -d0
911 else if (i .eq. nz) then
912 d(i, j) = d0
913 else
914 d(i, j) = 0.0_rp
915 end if
916 dt(j, i) = d(i, j)
917 end do
918 end do
919 end subroutine dgll
920
923 real(kind=xp) function hgll(I, Z, ZGLL, NZ)
924 integer, intent(in) :: i, nz
925 real(kind=xp), intent(in) :: zgll(nz), z
926
927 real(kind=xp) :: eps, dz
928 real(kind=xp) :: alfan
929 integer :: n
930
931 eps = 1.0e-5_xp
932 dz = z - zgll(i)
933 if (abs(dz) .lt. eps) then
934 hgll = 1.0_xp
935 return
936 end if
937
938 n = nz - 1
939 alfan = real(n, kind=xp) * (real(n, kind=xp) + 1.0_xp)
940 hgll = -(1.0_xp - z*z)*pndleg(z, n) / (alfan*pnleg(zgll(i), n) * &
941 (z - zgll(i)))
942 end function hgll
943
946 real(kind=xp) function hgl (I, Z, ZGL, NZ)
947 integer, intent(in) :: i, nz
948 real(kind=xp), intent(in) :: zgl(nz), z
949 real(kind=xp) :: eps, dz
950
951 integer :: n
952
953 eps = 1.0e-5_xp
954 dz = z - zgl(i)
955 if (abs(dz) .lt. eps) then
956 hgl = 1.0_xp
957 return
958 end if
959
960 n = nz - 1
961 hgl = pnleg(z, nz) / (pndleg(zgl(i), nz) * (z - zgl(i)))
962 end function hgl
963
967 real(kind=xp) function pnleg(Z, N)
968
969!---------------------------------------------------------------------
970!
971! This next statement is to overcome the underflow bug in the i860.
972! It can be removed at a later date. 11 Aug 1990 pff.
973!
974! IMPLICIT REAL(KIND=XP) (A-H,O-Z)
975! REAL(KIND=XP) Z, P1, P2, P3
976! !IF(ABS(Z) .LT. 1.0E-25) Z = 0.0
977
978 real(kind=xp), intent(in) :: z
979 integer, intent(in) :: n
980
981 real(kind=xp) :: p1, p2, p3, fk
982 integer :: k
983
984 p1 = 1.0_xp
985 if (n .eq. 0) then
986 pnleg = p1
987 return
988 end if
989
990 p2 = z
991 p3 = p2
992 do k = 1, n-1
993 fk = real(k, kind=xp)
994 p3 = ((2.0_xp*fk + 1.0_xp)*z*p2 - fk*p1) / (fk + 1.0_xp)
995 p1 = p2
996 p2 = p3
997 end do
998 pnleg = p3
999 end function pnleg
1000
1003 subroutine legendre_poly(L, x, N)
1004 integer, intent(in) :: N
1005 real(kind=rp), intent(inout) :: l(0:n)
1006 real(kind=rp), intent(in) :: x
1007
1008 real(kind=rp) :: dj
1009 integer :: j
1010
1011 l(0) = 1.0_rp
1012 if (n .eq. 0) return
1013 l(1) = x
1014
1015 do j = 1, n-1
1016 dj = real(j, kind=rp)
1017 l(j + 1) = ((2.0_rp*dj + 1.0_rp)*x*l(j) - dj*l(j-1)) / (dj + 1.0_rp)
1018 end do
1019 end subroutine legendre_poly
1020
1024 real(kind=xp) function pndleg(Z, N)
1025 real(kind=xp), intent(in) :: z
1026 integer, intent(in) :: n
1027
1028 real(kind=xp) :: p1, p2, p3, p1d, p2d, p3d, fk
1029 integer :: k
1030
1031 if (n .eq. 0) then
1032 pndleg = 0.0_xp
1033 return
1034 end if
1035
1036 p1 = 1.0_xp
1037 p2 = z
1038 p1d = 0.0_xp
1039 p2d = 1.0_xp
1040 p3d = 1.0_xp
1041 do k = 1, n-1
1042 fk = real(k, kind=xp)
1043 p3 = ((2.0_xp*fk + 1.0_xp)*z*p2 - fk*p1) / (fk + 1.0_xp)
1044 p3d = ((2.0_xp*fk + 1.0_xp)*p2 + (2.0_xp*fk + 1.0_xp)*z*p2d - fk*p1d) / &
1045 (fk + 1.0_xp)
1046 p1 = p2
1047 p2 = p3
1048 p1d = p2d
1049 p2d = p3d
1050 end do
1051 pndleg = p3d
1052 end function pndleg
1053
1061 subroutine dgllgl(D, DT, ZM1, ZM2, IM12, NZM1, NZM2, ND1, ND2)
1062 integer, intent(in) :: NZM1, NZM2, ND1, ND2
1063 real(kind=xp), intent(inout) :: d(nd2, nd1), dt(nd1, nd2)
1064 real(kind=xp), intent(in) :: zm1(nd1), zm2(nd2), im12(nd2, nd1)
1065
1066 real(kind=xp) eps, zp, zq
1067 integer :: IP, JQ, NM1
1068
1069 if (nzm1 .eq. 1) then
1070 d(1, 1) = 0.0_xp
1071 dt(1, 1) = 0.0_xp
1072 return
1073 end if
1074 eps = 1.0e-6_xp
1075 nm1 = nzm1 - 1
1076 do ip = 1, nzm2
1077 do jq = 1, nzm1
1078 zp = zm2(ip)
1079 zq = zm1(jq)
1080 if ((abs(zp) .lt. eps) .and. (abs(zq) .lt. eps)) then
1081 d(ip, jq) = 0.0_xp
1082 else
1083 d(ip, jq) = (pnleg(zp, nm1) / pnleg(zq, nm1) - im12(ip, jq)) / &
1084 (zp - zq)
1085 end if
1086 dt(jq, ip) = d(ip, jq)
1087 end do
1088 end do
1089 end subroutine dgllgl
1090
1099 subroutine dgljgj(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
1100 integer, intent(in) :: NPGL, NPG, ND1, ND2
1101 real(kind=xp), intent(inout) :: d(nd2, nd1), dt(nd1, nd2)
1102 real(kind=xp), intent(in) :: zgl(nd1), zg(nd2), iglg(nd2, nd1), alpha, beta
1103
1104 integer, parameter :: NMAX = 84
1105 integer, parameter :: NDD = nmax
1106
1107 real(kind=xp) dd(ndd, ndd), dtd(ndd, ndd)
1108 real(kind=xp) zgd(ndd), zgld(ndd), iglgd(ndd, ndd)
1109 integer :: I, J
1110
1111 if (npgl .le. 1) then
1112 call neko_error('DGLJGJ: Minimum number of Gauss-Lobatto points is 2')
1113 else if (npgl .gt. nmax) then
1114 write(stderr, *) 'Polynomial degree too high in DGLJGJ'
1115 write(stderr, *) 'Maximum polynomial degree is', nmax
1116 write(stderr, *) 'Here NPGL=', npgl
1117 call neko_error
1118 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
1119 call neko_error('DGLJGJ: Alpha and Beta must be greater than -1')
1120 end if
1121
1122 do i = 1, npg
1123 zgd(i) = zg(i)
1124 do j = 1, npgl
1125 iglgd(i, j) = iglg(i, j)
1126 end do
1127 end do
1128 do i = 1, npgl
1129 zgld(i) = zgl(i)
1130 end do
1131 call dgljgjd(dd, dtd, zgld, zgd, iglgd, npgl, npg, ndd, ndd, alpha, beta)
1132 do i = 1, npg
1133 do j = 1, npgl
1134 d(i, j) = dd(i, j)
1135 dt(j, i) = dtd(j, i)
1136 end do
1137 end do
1138 end subroutine dgljgj
1139
1148 subroutine dgljgjd(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
1149 integer, intent(in) :: NPGL, NPG, ND1, ND2
1150 real(kind=xp), intent(inout) :: d(nd2, nd1), dt(nd1, nd2)
1151 real(kind=xp), intent(in) :: zgl(nd1), zg(nd2), iglg(nd2, nd1), alpha, beta
1152
1153 real(kind=xp) :: eps, eigval, dn
1154 real(kind=xp) :: pdi, pdj, pi, pj, pm1, pdm1, pm2, pdm2
1155 real(kind=xp) :: dz, faci, facj, const
1156 integer :: I, J, NGL
1157
1158 if (npgl .le. 1) then
1159 call neko_error('DGLJGJD: Minimum number of Gauss-Lobatto points is 2')
1160 else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
1161 call neko_error('DGLJGJD: Alpha and Beta must be greater than -1')
1162 end if
1163
1164 eps = 1.0e-6_xp
1165
1166 ngl = npgl-1
1167 dn = real(ngl, kind=xp)
1168 eigval = -dn*(dn + alpha + beta + 1.0_xp)
1169
1170 do i = 1, npg
1171 do j = 1, npgl
1172 dz = abs(zg(i)-zgl(j))
1173 if (dz .lt. eps) then
1174 d(i, j) = (alpha*(1.0_xp + zg(i)) - beta*(1.0_xp - zg(i))) / &
1175 (2.0_xp*(1.0_xp - zg(i)**2))
1176 else
1177 call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, ngl, alpha, beta, zg(i))
1178 call jacobf(pj, pdj, pm1, pdm1, pm2, pdm2, ngl, alpha, &
1179 beta, zgl(j))
1180 faci = alpha*(1.0_xp + zg(i)) - beta*(1.0_xp - zg(i))
1181 facj = alpha*(1.0_xp + zgl(j)) - beta*(1.0_xp - zgl(j))
1182 const = eigval*pj + facj*pdj
1183 d(i, j) = ((eigval*pi + faci*pdi) * (zg(i) - zgl(j)) - &
1184 (1.0_xp - zg(i)**2)*pdi) / (const*(zg(i) - zgl(j))**2)
1185 end if
1186 dt(j, i) = d(i, j)
1187 end do
1188 end do
1189 end subroutine dgljgjd
1190
1196 subroutine iglm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
1197 integer, intent(in) :: NZ1, NZ2, ND1, ND2
1198 real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
1199 real(kind=xp), intent(in) :: z1(nd1), z2(nd2)
1200 real(kind=xp) :: zi
1201 integer :: I, J
1202
1203 if (nz1 .eq. 1) then
1204 i12(1, 1) = 1.0_xp
1205 it12(1, 1) = 1.0_xp
1206 return
1207 end if
1208
1209 do i = 1, nz2
1210 zi = z2(i)
1211 do j = 1, nz1
1212 i12(i, j) = hgl(j, zi, z1, nz1)
1213 it12(j, i) = i12(i, j)
1214 end do
1215 end do
1216 end subroutine iglm
1217
1223 subroutine igllm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
1224 integer, intent(in) :: NZ1, NZ2, ND1, ND2
1225 real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
1226 real(kind=xp), intent(in) :: z1(nd1), z2(nd2)
1227 real(kind=xp) :: zi
1228 integer :: I, J
1229
1230 if (nz1 .eq. 1) then
1231 i12(1, 1) = 1.0_xp
1232 it12(1, 1) = 1.0_xp
1233 return
1234 end if
1235
1236 do i = 1, nz2
1237 zi = z2(i)
1238 do j = 1, nz1
1239 i12(i, j) = hgll(j, zi, z1, nz1)
1240 it12(j, i) = i12(i, j)
1241 end do
1242 end do
1243 end subroutine igllm
1244
1251 subroutine igjm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
1252 integer, intent(in) :: NZ1, NZ2, ND1, ND2
1253 real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
1254 real(kind=xp), intent(in) :: z1(nd1), z2(nd2), alpha, beta
1255 real(kind=xp) :: zi
1256 integer :: I, J
1257
1258 if (nz1 .eq. 1) then
1259 i12(1, 1) = 1.0_xp
1260 it12(1, 1) = 1.0_xp
1261 return
1262 end if
1263
1264 do i = 1, nz2
1265 zi = z2(i)
1266 do j = 1, nz1
1267 i12(i, j) = hgj(j, zi, z1, nz1, alpha, beta)
1268 it12(j, i) = i12(i, j)
1269 end do
1270 end do
1271 end subroutine igjm
1272
1279 subroutine igljm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
1280 integer, intent(in) :: NZ1, NZ2, ND1, ND2
1281 real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
1282 real(kind=xp), intent(in) :: z1(nd1), z2(nd2), alpha, beta
1283 real(kind=xp) :: zi
1284 integer :: I, J
1285
1286 if (nz1 .eq. 1) then
1287 i12(1, 1) = 1.0_xp
1288 it12(1, 1) = 1.0_xp
1289 return
1290 end if
1291
1292 do i = 1, nz2
1293 zi = z2(i)
1294 do j = 1, nz1
1295 i12(i, j) = hglj(j, zi, z1, nz1, alpha, beta)
1296 it12(j, i) = i12(i, j)
1297 end do
1298 end do
1299 end subroutine igljm
1300end module speclib
real
double real
Definition device_config.h:12
utils::neko_error
Definition utils.f90:45
math
Definition math.f90:60
math::pi
real(kind=rp), parameter, public pi
Definition math.f90:76
num_types
Definition num_types.f90:1
num_types::xp
integer, parameter, public xp
Definition num_types.f90:14
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
speclib
LIBRARY ROUTINES FOR SPECTRAL METHODS.
Definition speclib.f90:149
speclib::dgj
subroutine dgj(d, dt, z, nz, nzd, alpha, beta)
Compute the derivative matrix D and its transpose DT associated with the Nth order Lagrangian interpo...
Definition speclib.f90:716
speclib::pndleg
real(kind=xp) function pndleg(z, n)
Compute the derivative of the Nth order Legendre polynomial at Z. (Simpler than JACOBF) Based on the ...
Definition speclib.f90:1025
speclib::dgll
subroutine dgll(d, dt, z, nz, nzd)
Compute the derivative matrix D and its transpose DT associated with the Nth order Lagrangian interpo...
Definition speclib.f90:881
speclib::igjm
subroutine igjm(i12, it12, z1, z2, nz1, nz2, nd1, nd2, alpha, beta)
Compute the one-dimensional interpolation operator (matrix) I12 ands its transpose IT12 for interpola...
Definition speclib.f90:1252
speclib::hgj
real(kind=xp) function hgj(ii, z, zgj, np, alpha, beta)
Compute the value of the Lagrangian interpolant HGJ through the NP Gauss Jacobi points ZGJ at the poi...
Definition speclib.f90:608
speclib::dgllgl
subroutine dgllgl(d, dt, zm1, zm2, im12, nzm1, nzm2, nd1, nd2)
Compute the (one-dimensional) derivative matrix D and its transpose DT associated with taking the der...
Definition speclib.f90:1062
speclib::zwgll
subroutine zwgll(z, w, np)
Generate NP Gauss-Lobatto Legendre points (Z) and weights (W) associated with Jacobi polynomial P(N)(...
Definition speclib.f90:180
speclib::pnormj
real(kind=xp) function pnormj(n, alpha, beta)
Definition speclib.f90:466
speclib::dgjd
subroutine dgjd(d, dt, z, nz, nzd, alpha, beta)
Compute the derivative matrix D and its transpose DT associated with the Nth order Lagrangian interpo...
Definition speclib.f90:755
speclib::gammaf
real(kind=xp) function gammaf(x)
Definition speclib.f90:447
speclib::jacg
subroutine jacg(xjac, np, alpha, beta)
Compute NP Gauss points XJAC, which are the zeros of the Jacobi polynomial J(NP) with parameters ALPH...
Definition speclib.f90:501
speclib::zwgljd
subroutine zwgljd(z, w, np, alpha, beta)
Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W) associated with Jacobi polynomial P(N)(al...
Definition speclib.f90:308
speclib::hglj
real(kind=xp) function hglj(ii, z, zglj, np, alpha, beta)
Compute the value of the Lagrangian interpolant HGLJ through the NZ Gauss-Lobatto Jacobi points ZGLJ ...
Definition speclib.f90:658
speclib::hgljd
real(kind=xp) function hgljd(i, z, zglj, np, alpha, beta)
Compute the value of the Lagrangian interpolant HGLJD through the NZ Gauss-Lobatto Jacobi points ZJAC...
Definition speclib.f90:685
speclib::zwglj
subroutine zwglj(z, w, np, alpha, beta)
Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W) associated with Jacobi polynomial P(N)(al...
Definition speclib.f90:276
speclib::zwgj
subroutine zwgj(z, w, np, alpha, beta)
Generate NP GAUSS JACOBI points (Z) and weights (W) associated with Jacobi polynomial P(N)(alpha>-1,...
Definition speclib.f90:193
speclib::iglm
subroutine iglm(i12, it12, z1, z2, nz1, nz2, nd1, nd2)
Compute the one-dimensional interpolation operator (matrix) I12 ands its transpose IT12 for interpola...
Definition speclib.f90:1197
speclib::jacobf
subroutine jacobf(poly, pder, polym1, pderm1, polym2, pderm2, n, alp, bet, x)
Computes the Jacobi polynomial (POLY) and its derivative (PDER) of degree N at X.
Definition speclib.f90:560
speclib::zwgjd
subroutine zwgjd(z, w, np, alpha, beta)
Generate NP GAUSS JACOBI points (Z) and weights (W) associated with Jacobi polynomial P(N)(alpha>-1,...
Definition speclib.f90:225
speclib::dglj
subroutine dglj(d, dt, z, nz, nzd, alpha, beta)
Compute the derivative matrix D and its transpose DT associated with the Nth order Lagrangian interpo...
Definition speclib.f90:794
speclib::dgljgj
subroutine dgljgj(d, dt, zgl, zg, iglg, npgl, npg, nd1, nd2, alpha, beta)
Compute the (one-dimensional) derivative matrix D and its transpose DT associated with taking the der...
Definition speclib.f90:1100
speclib::pnleg
real(kind=xp) function pnleg(z, n)
Compute the value of the Nth order Legendre polynomial at Z. (Simpler than JACOBF) Based on the recur...
Definition speclib.f90:968
speclib::hgl
real(kind=xp) function hgl(i, z, zgl, nz)
Compute the value of the Lagrangian interpolant HGL through the NZ Gauss Legendre points ZGL at the p...
Definition speclib.f90:947
speclib::dgljd
subroutine dgljd(d, dt, z, nz, nzd, alpha, beta)
Compute the derivative matrix D and its transpose DT associated with the Nth order Lagrangian interpo...
Definition speclib.f90:833
speclib::legendre_poly
subroutine legendre_poly(l, x, n)
Evaluate Legendre polynomials of degrees 0-N at point x and store in array L.
Definition speclib.f90:1004
speclib::endw1
real(kind=xp) function endw1(n, alpha, beta)
Definition speclib.f90:348
speclib::igllm
subroutine igllm(i12, it12, z1, z2, nz1, nz2, nd1, nd2)
Compute the one-dimensional interpolation operator (matrix) I12 ands its transpose IT12 for interpola...
Definition speclib.f90:1224
speclib::dgljgjd
subroutine dgljgjd(d, dt, zgl, zg, iglg, npgl, npg, nd1, nd2, alpha, beta)
Compute the (one-dimensional) derivative matrix D and its transpose DT associated with taking the der...
Definition speclib.f90:1149
speclib::zwgl
subroutine zwgl(z, w, np)
Generate NP Gauss Legendre points Z and weights W associated with Jacobi polynomial ....
Definition speclib.f90:165
speclib::hgll
real(kind=xp) function hgll(i, z, zgll, nz)
Compute the value of the Lagrangian interpolant L through the NZ Gauss-Lobatto Legendre points ZGLL a...
Definition speclib.f90:924
speclib::hgjd
real(kind=xp) function hgjd(ii, z, zgj, np, alpha, beta)
Compute the value of the Lagrangian interpolant HGJD through the NZ Gauss-Lobatto Jacobi points ZGJ a...
Definition speclib.f90:636
speclib::endw2
real(kind=xp) function endw2(n, alpha, beta)
Definition speclib.f90:397
speclib::igljm
subroutine igljm(i12, it12, z1, z2, nz1, nz2, nd1, nd2, alpha, beta)
Compute the one-dimensional interpolation operator (matrix) I12 ands its transpose IT12 for interpola...
Definition speclib.f90:1280
utils
Utilities.
Definition utils.f90:35