Commit 1c33df58 authored by Elia Merzari's avatar Elia Merzari
Browse files

Replace nek_moose.f with addition of spherical harmonics

parent 51647a0e
Pipeline #3447 failed with stages
......@@ -9,7 +9,7 @@ C=======================================================================
integer icall, kstep, i, pstep
common /cht_coupler/ pstep
save icall
if (icall.eq.0) then
call nekgsync()
if (instep.eq.0) then
......@@ -62,9 +62,9 @@ C=======================================================================
include 'INPUT'
include 'CTIMER'
common /cht_coupler/ pstep
integer pstep
integer pstep
logical ifdoin
character*3 prefix
character*3 prefix
real*4 papi_mflops
integer*8 papi_flops
integer icall, kstep, knstep, i
......@@ -76,7 +76,7 @@ C=======================================================================
call userchk
prefix='his'
ifdoin=.false.
ifdoin=.false.
call prepost (ifdoin,prefix)
call in_situ_check()
......@@ -91,7 +91,7 @@ C=======================================================================
include 'SIZE'
include 'TOTAL'
include 'NEKMOOSE'
include 'NEKMOOSE'
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
......@@ -131,11 +131,60 @@ c For Testing
c call nek_testp()
return
end
C=======================================================================
subroutine nek_expansion_sh()
include 'SIZE'
include 'TOTAL'
include 'NEKMOOSE'
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
real*8 sint, sint1, sarea, sarea1
real*8 pi
pi=4.0*atan(1.0)
ntot=nx1*ny1*nz1*nelt
call rzero(fmode,ntot)
do i0=1,n_legendre
i1=i0-1
n_four1=2*(i0-1)+1
do j0=1,n_four1
j1=j0-i1-1
call nek_mode_sh(fmode,i1,j1,0.0,0.0,0.0)
sarea1=0.0
sint1= 0.0
do e=1,nelt
do f=1,6
sint=0.0
sarea=0.0
if (cbc(f,e,1).eq.'W ') then
call col3(cache,fmode,t,ntot)
call surface_int(sint,sarea,cache,e,f)
sint1=sint1+sint
sarea1=sarea1+sarea
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
call gop(sarea1,wtmp,'+ ',1)
!
coeff_tij(i0,j0)=sint1*4.0*pi/sarea1
!
enddo
enddo
c For Testing
c call nek_testp()
return
end
C=======================================================================
subroutine nek_diag()
include 'SIZE'
include 'TOTAL'
include 'NEKMOOSE'
include 'NEKMOOSE'
common/diags_coeff/diag_c(nl_max,nf_max)
integer e,f
......@@ -151,7 +200,7 @@ C=======================================================================
call rzero(fmode,ntot)
do i0=1,n_legendre
do j0=1,m_fourier
call nek_mode(fmode,i0,j0)
call nek_mode(fmode,i0,j0)
sint1=0.0
do e=1,nelt
do f=1,6
......@@ -175,28 +224,28 @@ C=======================================================================
subroutine nek_testp()
include 'SIZE'
include 'TOTAL'
include 'NEKMOOSE'
include 'NEKMOOSE'
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
real*8 fun(lx1,ly1,lz1,lelt)
real*8 fun(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
call rzero(fmode,ntot)
call rzero(fun,ntot)
call rzero(fun,ntot)
do i0=1,n_legendre
do j0=1,m_fourier
call nek_mode(fmode,i0,j0)
do i=1,ntot
fun(i,1,1,1)=fun(i,1,1,1)+fmode(i,1,1,1)*coeff_tij(i0,j0)
enddo
fun(i,1,1,1)=fun(i,1,1,1)+fmode(i,1,1,1)*coeff_tij(i0,j0)
enddo
enddo
enddo
call rzero(cache,ntot)
call rzero(cache,ntot)
call sub3(cache,fun,t,ntot)
do i=1,ntot
c1=cache(i,1,1,1)**2
cache(i,1,1,1)=c1
cache(i,1,1,1)=c1
enddo
sint1=0.0
sarea1=0.0
do e=1,lelt
......@@ -220,7 +269,7 @@ C=======================================================================
include 'SIZE'
include 'TOTAL'
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt)
real*8 fmode(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
zmin=glmin(zm1,ntot)
zmax=glmax(zm1,ntot)
......@@ -250,21 +299,66 @@ C=======================================================================
enddo
enddo
return
end
end
C=======================================================================
subroutine nek_mode_sh(fmode,im,jm,x0,y0,z0)
include 'SIZE'
include 'TOTAL'
integer e,f
real*8 x0,y0,z0,r0
real*8 fmode(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
zmin=glmin(zm1,ntot)
zmax=glmax(zm1,ntot)
do e=1,nelt
do f=1,6
if (cbc(f,e,1).eq.'W ') then
call dsset(nx1,ny1,nz1)
iface = eface1(f)
js1 = skpdat(1,iface)
jf1 = skpdat(2,iface)
jskip1 = skpdat(3,iface)
js2 = skpdat(4,iface)
jf2 = skpdat(5,iface)
jskip2 = skpdat(6,iface)
do j2=js2,jf2,jskip2
do j1=js1,jf1,jskip1
x = xm1(j1,j2,1,e)-x0
y = ym1(j1,j2,1,e)-y0
z = zm1(j1,j2,1,e)-z0
c phi is the latitude
phi = z/(x**2 + y**2 + z**2)
rr2=(xx**2+yy**2+zz**2)
c theta is the longitude
theta=atan2(y,x)
call pmns_polynomial_value(1, im, jm, phi, cx )
fmode(j1,j2,1,e)=cx*fl_four(theta,jm)/rr2
enddo
enddo
endif
enddo
enddo
return
end
c-----------------------------------------------------------------------
subroutine flux_reconstruction()
include 'SIZE'
include 'TOTAL'
include 'NEKMOOSE'
include 'NEKMOOSE'
integer e,f
real*8 coeff_base, flux_base, flux_moose
real*8 coeff_base, flux_base, flux_moose
real*8 fmode(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
c --------------------------------------
c The flux from MOOSE must have proper sign
c --------------------------------------
coeff_base=-1.0
coeff_base=-1.0
flux_base=0.25 ! from energy conservation in the problem
c --------------------------------------
......@@ -438,22 +532,22 @@ c $ ,TMULT(1,1,1,1,IFIELD-1)
c $ ,IMESH,TOLHT(IFIELD),NMXH,isd)
if(iftmsh(ifield)) then
call hsolve (name4,TA,TB,H1,H2
call hsolve (name4,TA,TB,H1,H2
$ ,tmask(1,1,1,1,ifield-1)
$ ,tmult(1,1,1,1,ifield-1)
$ ,imesh,tolht(ifield),nmxh,1
$ ,approx,napprox,bintm1)
else
call hsolve (name4,TA,TB,H1,H2
call hsolve (name4,TA,TB,H1,H2
$ ,tmask(1,1,1,1,ifield-1)
$ ,tmult(1,1,1,1,ifield-1)
$ ,imesh,tolht(ifield),nmxh,1
$ ,approx,napprox,binvm1)
endif
endif
call add2 (t(1,1,1,1,ifield-1),ta,n)
call cvgnlps (ifconv) ! Check convergence for nonlinear problem
call cvgnlps (ifconv) ! Check convergence for nonlinear problem
if (ifconv) goto 2000
C Radiation case, smooth convergence, avoid flip-flop (ER).
......@@ -478,21 +572,21 @@ C (1) Varaiable properties.
C (2) Implicit time stepping.
C (3) User specified tolerance for the Helmholtz solver
C (not based on eigenvalues).
C (4) A passive scalar can be defined on either the
C (4) A passive scalar can be defined on either the
C temperatur or the velocity mesh.
C (5) A passive scalar has its own multiplicity (B.C.).
C (5) A passive scalar has its own multiplicity (B.C.).
C
include 'SIZE'
include 'INPUT'
include 'TSTEP'
include 'TURBO'
include 'TURBO'
include 'DEALIAS'
real*8 ts, dnekclock
ts = dnekclock()
if (nio.eq.0 .and. igeom.eq.2)
if (nio.eq.0 .and. igeom.eq.2)
& write(*,'(13x,a)') 'Solving for Hmholtz scalars'
do ifield = 2,nfield
......@@ -507,7 +601,7 @@ C
enddo
if (nio.eq.0 .and. igeom.eq.2)
& write(*,'(4x,i7,a,1p2e12.4)')
& write(*,'(4x,i7,a,1p2e12.4)')
& istep,' Scalars done',time,dnekclock()-ts
return
......@@ -537,9 +631,9 @@ c-----------------------------------------------------------------------
do igeom=1,ngeom
! within cvode we use the lagged wx for
! extrapolation, that's why we have to call it before gengeom
if (ifheat .and. ifcvode) call heat_cvode (igeom)
! within cvode we use the lagged wx for
! extrapolation, that's why we have to call it before gengeom
if (ifheat .and. ifcvode) call heat_cvode (igeom)
if (ifgeom) then
call gengeom (igeom)
......@@ -548,7 +642,7 @@ c-----------------------------------------------------------------------
if (ifheat) call heat_mod (igeom)
if (igeom.eq.2) then
if (igeom.eq.2) then
call setprop
if (iflomach) call qthermal(.true.,.false.,dummy)
endif
......@@ -556,7 +650,7 @@ c-----------------------------------------------------------------------
if (ifflow) call fluid (igeom)
if (ifmvbd) call meshv (igeom)
if (param(103).gt.0) call q_filter (param(103))
call setup_convect (igeom) ! Save convective velocity _after_ filter
call setup_convect (igeom) ! Save convective velocity _after_ filter
enddo
else ! PN-2/PN-2 formulation
......@@ -587,7 +681,7 @@ c-----------------------------------------------------------------------
if (ifmvbd) call meshv (igeom)
endif
if (igeom.eq.ngeom.and.param(103).gt.0)
if (igeom.eq.ngeom.and.param(103).gt.0)
$ call q_filter(param(103))
call setup_convect (igeom) ! Save convective velocity _after_ filter
......@@ -597,5 +691,163 @@ c-----------------------------------------------------------------------
return
end
c-----------------------------------------------------------------------
subroutine pmns_polynomial_value ( mm, n, m, x, cx )
c*********************************************************************72
c
cc PMNS_POLYNOMIAL_VALUE: sphere-normalized Legendre polynomial Pmns(n,m,x).
c
c Discussion:
c
c The unnormalized associated Legendre functions P_N^M(X) have
c the property that
c
c Integral ( -1 .le. X .le. 1 ) ( P_N^M(X) )^2 dX
c = 2 * ( N + M )! / ( ( 2 * N + 1 ) * ( N - M )! )
c
c By dividing the function by the square root of this term,
c the normalized associated Legendre functions have norm 1.
c
c However, we plan to use these functions to build spherical
c harmonics, so we use a slightly different normalization factor of
c
c sqrt ( ( ( 2 * N + 1 ) * ( N - M )! ) / ( 4 * pi * ( N + M )! ) )
c
c Licensing:
c
c This code is distributed under the GNU LGPL license.
c
c Modified:
c
c 07 August 2013
c
c Author:
c
c John Burkardt
c
c Reference:
c
c Milton Abramowitz, Irene Stegun,
c Handbook of Mathematical Functions,
c National Bureau of Standards, 1964,
c ISBN: 0-486-61272-4,
c LC: QA47.A34.
c
c Parameters:
c
c Input, integer MM, the number of evaluation points.
c
c Input, integer N, the maximum first index of the Legendre
c function, which must be at least 0.
c
c Input, integer M, the second index of the Legendre function,
c which must be at least 0, and no greater than N.
c
c Input, double precision X(MM), the evaluation points.
c
c Output, double precision CX(MM,0:N), the function values.
c
implicit none
integer mm
integer n
double precision cx(mm,0:n)
double precision factor
integer i
integer j
integer m
double precision pi
parameter ( pi = 3.141592653589793D+00 )
double precision r8_factorial
double precision x(mm)
do j = 0, n
do i = 1, mm
cx(i,j) = 0.0D+00
end do
end do
if ( m .le. n ) then
do i = 1, mm
cx(1:mm,m) = 1.0D+00
end do
factor = 1.0D+00
do j = 1, m
do i = 1, mm
cx(i,m) = - cx(i,m) * factor * sqrt ( 1.0D+00 - x(i)**2 )
end do
factor = factor + 2.0D+00
end do
end if
if ( m + 1 .le. n ) then
do i = 1, mm
cx(i,m+1) = x(i) * dble ( 2 * m + 1 ) * cx(i,m)
end do
end if
do j = m + 2, n
do i = 1, mm
cx(i,j) = ( dble ( 2 * j - 1 ) * x(i) * cx(i,j-1)
& + dble ( - j - m + 1 ) * cx(i,j-2) )
& / dble ( j - m )
end do
end do
c
c Normalization.
c
do j = m, n
factor = sqrt ( ( dble ( 2 * j + 1 ) * r8_factorial ( j - m ) )
& / ( 4.0D+00 * pi * r8_factorial ( j + m ) ) )
do i = 1, mm
cx(i,j) = cx(i,j) * factor
end do
end do
return
end
function r8_factorial ( n )
c*********************************************************************72
c
cc R8_FACTORIAL computes the factorial of N.
c
c Discussion:
c
c factorial ( N ) = product ( 1 .le. I .le. N ) I
c
c Licensing:
c
c This code is distributed under the GNU LGPL license.
c
c Modified:
c
c 07 June 2008
c
c Author:
c
c John Burkardt
c
c Parameters:
c
c Input, integer N, the argument of the factorial function.
c If N is less than 1, the function value is returned as 1.
c
c Output, double precision R8_FACTORIAL, the factorial of N.
c
implicit none
integer i
integer n
double precision r8_factorial
r8_factorial = 1.0D+00
do i = 1, n
r8_factorial = r8_factorial * dble ( i )
end do
return
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment