Skip to content

GitLab

Commit 3993263c authored by Elia Merzari's avatar Elia Merzari
Browse files

Adding projection to Legendre-Fourier, and reconstruction

parent 21bd7fcb
Hide whitespace changes
Inline Side-by-side
examples/integration_example/integration_example.usr
View file @ 3993263c
......@@ -47,8 +47,13 @@ c-----------------------------------------------------------------------
subroutine userchk
include 'SIZE'
include 'TOTAL'
common/test_passing/ flux_moose, temp_nek
parameter (nl_max=100)
parameter (nf_max=100)
common/expansion_tdata/n_legendre, m_fourier
common/expansion_tcoef/coeff_tij(nl_max,nf_max)
common/expansion_fcoef/coeff_fij(nl_max,nf_max)
common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
common/test_passing/ flux_moose, temp_nek
real sint, sint1, sarea, sarea1, wtmp
integer e, f
......@@ -57,10 +62,8 @@ c-----------------------------------------------------------------------
pmin=glmin(pr,n2)
pmax=glmax(pr,n2)
wmin=glmin(vz,n1)
wmax=glmax(vz,n1)
tmin=glmin(t,n1)
tmax=glmax(t,n1)
......@@ -86,9 +89,24 @@ c-----------------------------------------------------------------------
if (nid.eq.0) then
write(6,*)"*** Temperature: ",tmin," - ",tmax
write(6,*)"*** Av. Temperature: ",temp_nek
write(6,*)"*** Average Temperature: ",temp_nek
endif
c Will this be overwritten -----------------
write(6,*)"*** Setting polynomial orders ..."
n_legendre=10
m_fourier=10
c-----------------------------------------------
c This is for testing -----------
c call flux_reconstruction()
c call nek_expansion()
c call nek_diag()
c if (nid.eq.0) write(6,*)" *** Expansion done"
c call nek_testp()
c call exitt()
c--------------------------------
return
end
c-----------------------------------------------------------------------
......@@ -98,6 +116,7 @@ c-----------------------------------------------------------------------
include 'NEKUSE'
common/test_passing/ flux_moose, temp_nek
common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
integer e,ieg
real ucx, ucy, ucz, ucy_e, yy
......@@ -113,7 +132,7 @@ c-----------------------------------------------------------------------
uy = 0.0
uz = 1.0
temp = 0.0
flux = 1.0 !flux_moose
flux = flux_recon(ix,iy,iz,e) !flux_moose
return
end
......@@ -221,7 +240,6 @@ C=======================================================================
end
C=======================================================================
subroutine nek_finalize_step()
include 'SIZE'
include 'TSTEP'
include 'INPUT'
......@@ -244,6 +262,242 @@ C=======================================================================
call nek_comm_settings(isyc,0)
call comment
return
end
C=======================================================================
subroutine nek_expansion()
parameter (nl_max=100)
parameter (nf_max=100)
include 'SIZE'
include 'TOTAL'
common/expansion_tdata/n_legendre, m_fourier
common/expansion_tcoef/coeff_tij(nl_max,nf_max)
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
real*8 sint, sint1
ntot=nx1*ny1*nz1*nelt
zmin=glmin(zm1,ntot)
zmax=glmax(zm1,ntot)
call rzero(fmode,ntot)
do i0=1,n_legendre
do j0=1,m_fourier
call nek_mode(fmode,i0,j0)
sint1=0.0
do e=1,nelt
do f=1,6
sint=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
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
!
coeff_tij(i0,j0)=sint1*2.0/(0.5*(zmax-zmin))
! Note that R=0.5 here!!!!
enddo
enddo
return
end
C=======================================================================
subroutine nek_diag()
parameter (nl_max=100)
parameter (nf_max=100)
include 'SIZE'
include 'TOTAL'
common/expansion_tdata/n_legendre, m_fourier
common/expansion_tcoef/coeff_tij(nl_max,nf_max)
common/diags_coeff/diag_c(nl_max,nf_max)
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt)
real*8 cache(lx1,ly1,lz1,lelt)
real*8 sint, sint1
ntot=nx1*ny1*nz1*nelt
zmin=glmin(zm1,ntot)
zmax=glmax(zm1,ntot)
call rzero(fmode,ntot)
do i0=1,n_legendre
do j0=1,m_fourier
call nek_mode(fmode,i0,j0)
sint1=0.0
do e=1,nelt
do f=1,6
sint=0.0
if (cbc(f,e,1).eq.'W ') then
call col3(cache,fmode,fmode,ntot)
call surface_int(sint,sarea,cache,e,f)
sint1=sint1+sint
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
diag_c(i0,j0)=sint1*2.0/(0.5*(zmax-zmin))
if (nid.eq.0) write(6,*)i0,j0,diag_c(i0,j0)
enddo
enddo
return
end
C=======================================================================
subroutine nek_testp()
parameter (nl_max=100)
parameter (nf_max=100)
include 'SIZE'
include 'TOTAL'
common/expansion_tdata/n_legendre, m_fourier
common/expansion_tcoef/coeff_tij(nl_max,nf_max)
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
real*8 fun(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
call rzero(fmode,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
enddo
enddo
call rzero(cache,ntot)
call sub3(cache,fun,t,ntot)
sint1=0.0
sarea1=0.0
do e=1,lelt
do f=1,6
call surface_int(sint,sarea,cache,e,f)
if (cbc(f,e,1).eq.'W ') then
sint1=sint1+sint
sarea1=sarea1+sarea
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
call gop(sarea1,wtmp,'+ ',1)
er_avg=sint1/sarea1
if (nid.eq.0) write(6,*)"Error: ",er_avg
return
end
C=======================================================================
subroutine nek_mode(fmode,im,jm)
include 'SIZE'
include 'TOTAL'
integer e,f
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)
y=ym1(j1,j2,1,e)
z=zm1(j1,j2,1,e)
z_leg=2*((z-zmin)/zmax)-1
theta=atan2(y,x)
fmode(j1,j2,1,e)=
& pl_leg(z_leg,im-1)*fl_four(theta,jm-1)
enddo
enddo
endif
enddo
enddo
return
end
c-----------------------------------------------------------------------
subroutine flux_reconstruction()
include 'SIZE'
include 'TOTAL'
parameter (nl_max=100)
parameter (nf_max=100)
common/expansion_tdata/n_legendre, m_fourier
common/expansion_tcoef/coeff_tij(nl_max,nf_max)
common/expansion_fcoef/coeff_fij(nl_max,nf_max)
common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
integer e,f
real*8 fmode(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
call rzero(flux_recon,ntot)
call rzero(fmode,ntot)
do i0=1,n_legendre
do j0=1,m_fourier
call nek_mode(fmode,i0,j0)
do i=1,ntot
flux_recon(i,1,1,1)= flux_recon(i,1,1,1)
& + fmode(i,1,1,1)*coeff_fij(i0,j0)
enddo
enddo
enddo
c---- Below is for testing
c do i=1,ntot
c flux_recon(i,1,1,1)= 1.0
c enddo
c--------------------------
return
end
c-----------------------------------------------------------------------
function pl_leg(x,n)
!======================================
! calculates Legendre polynomials Pn(x)
! using the recurrence relation
! if n > 100 the function retuns 0.0
!======================================
real*8 pl,pl_leg
real*8 x
real*8 pln(0:n)
integer n, k
pln(0) = 1.0
pln(1) = x
if (n.le.1) then
pl = pln(n)
else
do k=1,n-1
pln(k+1)=
& ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
end do
pl = pln(n)
end if
pl_leg=pl*sqrt(dble(2*n+1)/2.0)
return
end
!======================================
function fl_four(x,n)
real*8 fl_four,pi
real*8 x, A
integer n, k
pi=4.0*atan(1.0)
A=1.0/sqrt(2*pi)
if (n.gt.0) A=1.0/sqrt(pi)
fl_four=A*cos(n*x)
c fl_four=1.0/sqrt(pi)
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