Commit bebe10d5 authored by Elia Merzari's avatar Elia Merzari

Adding fix for additional helmoltz solve ... due to update to Nek

parent 1a1f1163
......@@ -34,7 +34,7 @@ C=======================================================================
if (lastep .eq. 1) then
pstep=2
else
call nek_advance
call nek_advance_moose
pstep=2
endif
......@@ -51,7 +51,7 @@ C=======================================================================
integer pstep
pstep=pstep+1
call heat(pstep)
call heat_mod(pstep)
return
end
......@@ -349,3 +349,248 @@ c-----------------------------------------------------------------------
return
end
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
subroutine cdscal_mod (igeom)
C
C Solve the convection-diffusion equation for passive scalar IPSCAL
C
include 'SIZE'
include 'INPUT'
include 'GEOM'
include 'MVGEOM'
include 'SOLN'
include 'MASS'
include 'TSTEP'
COMMON /CPRINT/ IFPRINT
LOGICAL IFPRINT
LOGICAL IFCONV
C
COMMON /SCRNS/ TA(LX1,LY1,LZ1,LELT)
$ ,TB(LX1,LY1,LZ1,LELT)
COMMON /SCRVH/ H1(LX1,LY1,LZ1,LELT)
$ ,H2(LX1,LY1,LZ1,LELT)
include 'ORTHOT'
if (ifdg) then
call cdscal_dg(igeom)
return
endif
napprox(1) = laxt ! Fix this... pff 10/10/15
nel = nelfld(ifield)
n = nx1*ny1*nz1*nel
if (igeom.eq.1) then ! geometry at t^{n-1}
call makeq
call lagscal
else ! geometry at t^n
IF (IFPRINT) THEN
IF (IFMODEL .AND. IFKEPS) THEN
NFLDT = NFIELD - 1
IF (IFIELD.EQ.NFLDT.AND.NID.EQ.0) THEN
WRITE (6,*) ' Turbulence Model - k/epsilon solution'
ENDIF
ELSE
IF (IFIELD.EQ.2.AND.NID.EQ.0) THEN
WRITE (6,*) ' Temperature/Passive scalar solution'
ENDIF
ENDIF
ENDIF
if1=ifield-1
write(name4,1) if1-1
1 format('PS',i2)
if(ifield.eq.2) write(name4,'(A4)') 'TEMP'
C
C New geometry
C
isd = 1
if (ifaxis.and.ifaziv.and.ifield.eq.2) isd = 2
c if (ifaxis.and.ifmhd) isd = 2 !This is a problem if T is to be T!
do 1000 iter=1,nmxnl ! iterate for nonlin. prob. (e.g. radiation b.c.)
INTYPE = 0
IF (IFTRAN) INTYPE = -1
CALL SETHLM (H1,H2,INTYPE)
CALL BCNEUSC (TA,-1)
CALL ADD2 (H2,TA,N)
CALL BCDIRSC (T(1,1,1,1,IFIELD-1))
CALL AXHELM (TA,T(1,1,1,1,IFIELD-1),H1,H2,IMESH,isd)
CALL SUB3 (TB,BQ(1,1,1,1,IFIELD-1),TA,N)
CALL BCNEUSC (TA,1)
CALL ADD2 (TB,TA,N)
c CALL HMHOLTZ (name4,TA,TB,H1,H2
c $ ,TMASK(1,1,1,1,IFIELD-1)
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
$ ,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
$ ,tmask(1,1,1,1,ifield-1)
$ ,tmult(1,1,1,1,ifield-1)
$ ,imesh,tolht(ifield),nmxh,1
$ ,approx,napprox,binvm1)
endif
call add2 (t(1,1,1,1,ifield-1),ta,n)
call cvgnlps (ifconv) ! Check convergence for nonlinear problem
if (ifconv) goto 2000
C Radiation case, smooth convergence, avoid flip-flop (ER).
CALL CMULT (TA,0.5,N)
CALL SUB2 (T(1,1,1,1,IFIELD-1),TA,N)
1000 CONTINUE
2000 CONTINUE
CALL BCNEUSC (TA,1)
endif
return
end
c-----------------------------------------------------------------------
subroutine heat_mod (igeom)
C
C Driver for temperature or passive scalar.
C
C Current version:
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 temperatur or the velocity mesh.
C (5) A passive scalar has its own multiplicity (B.C.).
C
include 'SIZE'
include 'INPUT'
include 'TSTEP'
include 'TURBO'
include 'DEALIAS'
real*8 ts, dnekclock
ts = dnekclock()
if (nio.eq.0 .and. igeom.eq.2)
& write(*,'(13x,a)') 'Solving for Hmholtz scalars'
do ifield = 2,nfield
if (idpss(ifield-1).eq.0) then ! helmholtz
intype = -1
if (.not.iftmsh(ifield)) imesh = 1
if ( iftmsh(ifield)) imesh = 2
call unorm
call settolt
call cdscal_mod(igeom)
endif
enddo
if (nio.eq.0 .and. igeom.eq.2)
& write(*,'(4x,i7,a,1p2e12.4)')
& istep,' Scalars done',time,dnekclock()-ts
return
end
c-----------------------------------------------------------------------
subroutine nek_advance_moose
include 'SIZE'
include 'TOTAL'
include 'CTIMER'
common /cgeom/ igeom
call nekgsync
if (iftran) call settime
if (ifmhd ) call cfl_check
call setsolv
call comment
if (ifcmt) then
if (nio.eq.0.and.istep.le.1) write(6,*) 'CMT branch active'
call cmt_nek_advance
return
endif
if (ifsplit) then ! PN/PN formulation
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)
if (ifgeom) then
call gengeom (igeom)
call geneig (igeom)
endif
if (ifheat) call heat_mod (igeom)
if (igeom.eq.2) then
call setprop
if (iflomach) call qthermal(.true.,.false.,dummy)
endif
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
enddo
else ! PN-2/PN-2 formulation
call setprop
do igeom=1,ngeom
if (igeom.gt.2) call userchk_set_xfer
if (ifgeom) then
call gengeom (igeom)
call geneig (igeom)
endif
if (ifneknekm.and.igeom.eq.2) call multimesh_create
if (ifmhd) then
if (ifheat) call heat (igeom)
call induct (igeom)
elseif (ifpert) then
if (ifbase.and.ifheat) call heat (igeom)
if (ifbase.and.ifflow) call fluid (igeom)
if (ifflow) call fluidp (igeom)
if (ifheat) call heatp (igeom)
else ! std. nek case
if (ifheat) call heat_mod (igeom)
if (ifflow) call fluid (igeom)
if (ifmvbd) call meshv (igeom)
endif
if (igeom.eq.ngeom.and.param(103).gt.0)
$ call q_filter(param(103))
call setup_convect (igeom) ! Save convective velocity _after_ filter
enddo
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