Commit 481aa4d8 authored by Elia Merzari's avatar Elia Merzari

Update 1/1 Bringin in sync with working version

parent cb870daf
......@@ -69,9 +69,9 @@ c-----------------------------------------------------------------------
ifflow=.false.
if (istep.eq.0) then
call rzero(t,n1)
endif
c if (istep.eq.0) then
c call rzero(t,n1)
c endif
sint1=0.0
sarea1=0.0
......@@ -108,6 +108,8 @@ c-----------------------------------------------------------------------
call gop(sarea1,wtmp,'+ ',1)
flux_moose=sint1/sarea1
call heat_balance(flux_moose)
if (nid.eq.0) then
write(6,*)"*** Temperature: ",tmin," - ",tmax
......@@ -116,9 +118,9 @@ c-----------------------------------------------------------------------
endif
c Will this be overwritten -----------------
write(6,*)"*** Setting polynomial orders ..."
n_legendre=5
m_fourier=2
c write(6,*)"*** Setting polynomial orders ..."
n_legendre=10
m_fourier=5
c-----------------------------------------------
c This is for testing -----------
......@@ -338,6 +340,8 @@ C=======================================================================
enddo
enddo
c For Testing
call nek_testp()
return
end
C=======================================================================
......@@ -407,6 +411,10 @@ C=======================================================================
enddo
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
enddo
sint1=0.0
sarea1=0.0
......@@ -422,7 +430,7 @@ C=======================================================================
call gop(sint1,wtmp,'+ ',1)
call gop(sarea1,wtmp,'+ ',1)
er_avg=sint1/sarea1
er_avg=sqrt(sint1/sarea1)
if (nid.eq.0) write(6,*)"Error: ",er_avg
return
end
......@@ -473,7 +481,7 @@ c-----------------------------------------------------------------------
common/expansion_fcoef/coeff_fij(nl_max,nf_max)
common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
integer e,f
real*8 coeff_base
real*8 coeff_base, flux_base, flux_moose
real*8 fmode(lx1,ly1,lz1,lelt)
ntot=nx1*ny1*nz1*nelt
......@@ -481,6 +489,7 @@ c --------------------------------------
c The flux from MOOSE must have proper sign
c --------------------------------------
coeff_base=-1.0
flux_base=0.25 ! from energy conservation in the problem
c --------------------------------------
call rzero(flux_recon,ntot)
......@@ -497,9 +506,32 @@ c --------------------------------------
c---- Below is for testing
c do i=1,ntot
c flux_recon(i,1,1,1)= 1.0
c flux_recon(i,1,1,1)= 0.0
c enddo
c--------------------------
c---- Renormalization
sint1=0.0
sarea1=0.0
do e=1,lelt
do f=1,6
call surface_int(sint,sarea,flux_recon,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)
flux_moose=sint1/sarea1
do i=1,ntot
flux_recon(i,1,1,1)= flux_recon(i,1,1,1)*flux_base/flux_moose
enddo
c-----------------------
return
end
c-----------------------------------------------------------------------
......@@ -539,4 +571,73 @@ c-----------------------------------------------------------------------
c fl_four=1.0/sqrt(pi)
return
end
c-----------------------------------------------------------------------
subroutine heat_balance(fflux)
include 'SIZE'
include 'TOTAL'
real sint, sint1, sarea, sarea1, wtmp, fflux
real cache(lx1,ly1,lz1,lelt),dtdz(lx1,ly1,lz1,lelt),
& dtdy(lx1,ly1,lz1,lelt),
& dtdx(lx1,ly1,lz1,lelt)
integer e, f
n1=nelt*lx1*ly1*lz1
n2=nelt*lx2*ly2*lz2
nxyz=lx1*ly1*lz1
call gradm1(dtdx,dtdy,dtdz,t(1,1,1,1,1))
sint1=0.0
sarea1=0.0
do e=1,lelt
do i=1,nxyz
cache(i,1,1,e)=t(i,1,1,e,1)*vz(i,1,1,e)
enddo
do f=1,6
call surface_int(sint,sarea,cache,e,f)
if (cbc(f,e,1).eq.'O ') then
sint1=sint1+sint
sarea1=sarea1+sarea
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
test_flux=sint1
sint1=0.0
sarea1=0.0
do e=1,lelt
do f=1,6
call surface_int(sint,sarea,dtdz,e,f)
if (cbc(f,e,1).eq.'v ') then
sint1=sint1+sint
sarea1=sarea1+sarea
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
call gop(sarea1,wtmp,'+ ',1)
flux_inlet=sint1
sint1=0.0
sarea1=0.0
do e=1,lelt
do f=1,6
call surface_int(sint,sarea,vz,e,f)
if (cbc(f,e,1).eq.'W ') then
sarea1=sarea1+sarea
endif
enddo
enddo
call gop(sarea1,wtmp,'+ ',1)
if (nid.eq.0) then
rr_loss=
& flux_inlet*param(8)/(test_flux+flux_inlet*param(8))
write(6,*)"*** Heat balance: ",
& sarea1*fflux,test_flux+flux_inlet*param(8)
write(6,*)"*** Heat loss: ",rr_loss*100.0," %"
endif
return
end
c-----------------------------------------------------------------------
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