Test Program for 1-D Euler Equations



Test Program for 1-D Euler Equations

PROGRAM testFluid

!

! Program to test numerical methods

!

USE Input

USE FluidArrays

USE Trans

USE Matrix

USE Location

IMPLICIT NONE

!

CALL SetIC ! Set initial conditions

CALL SetLocAr ! Set arrays to cross reference data structures

CALL AllocMat ! Allocate Matrices

CALL InitFluid ! Initialize dependent variables

CALL Transient ! Run a transient

!

STOP

!

END

Define Basic Data Types

MODULE IntrType

IMPLICIT NONE

!

! Set parameters that determine the precision of floating point

! arithmetic and number of digits in integers

!

INTEGER, PARAMETER :: sdk = selected_real_kind(13,307)

INTEGER, PARAMETER :: sik = kind(10000000)

END MODULE IntrType

Data Structures

I find that the best way to start a program is to make a list of all the information that you will need. Determine its scope in the program (Global or Local), and establish data structures to carry this information. I usually begin with scalar data.

MODULE ScalarDat

USE IntrType

!

! Global scalar data needed by the program

!

REAL(sdk) :: dt = 0.1_sdk ! time step

REAL(sdk) :: time = 0.0_sdk ! time

REAL(sdk) :: endTime = 1.0_sdk ! end time for the transient

REAL(sdk) :: eps = 1.e-8_sdk ! convergence criterion

REAL(sdk) :: residMax ! maximum factional residual

INTEGER(sik) :: it ! current iteration number

INTEGER(sik) :: itmax = 10 ! number of iterations

INTEGER(sik) :: nsmax = 100000 ! maximum number of time steps

INTEGER(sik) :: nstep = 0 ! current time step number

INTEGER(sik) :: ncell = 10 ! number of computational volumes

INTEGER(sik) :: nvar ! total number of independent variables

INTEGER(sik) :: nbnd = 1 ! number of boundary cells

INTEGER(sik) :: ivLB ! lower bound of velocity indices

INTEGER(sik) :: ivUB ! upper bound of velocity indices

CHARACTER(LEN=1) :: lBC, uBC ! Boundary condition flags

! 'p' for pressure boundary condition

! 'v' for velocity boundary condition

LOGICAL :: converged ! test on iteration convergence

END MODULE ScalarDat

Physical Array Data

MODULE FluidArrays

USE IntrType

!

! Arrays containing state information about the fluid in 1-D flow

! All use SI units

!

! Arrays names ending with "N" are evaluated at the new (n+1) time level

! others are values at the start of the time step.

!

REAL(sdk), ALLOCATABLE :: p(:) ! Start of step pressure

REAL(sdk), ALLOCATABLE :: pN(:) ! End of step pressure

REAL(sdk), ALLOCATABLE :: T(:) ! Start of step temperature

REAL(sdk), ALLOCATABLE :: TN(:) ! End of step temperature

REAL(sdk), ALLOCATABLE :: e(:) ! Start of step specific internal energy

REAL(sdk), ALLOCATABLE :: eN(:) ! End of step specific internal energy

REAL(sdk), ALLOCATABLE :: rho(:) ! Start of step density

REAL(sdk), ALLOCATABLE :: rhoN(:) ! End of step density

REAL(sdk), ALLOCATABLE :: rhoe(:) ! Start of step internal energy per vol

REAL(sdk), ALLOCATABLE :: rhoeN(:) ! End of step internal energy per vol

REAL(sdk), ALLOCATABLE :: v(:) ! Start of step cell edge velocity

REAL(sdk), ALLOCATABLE :: vN(:) ! End of step velocity

REAL(sdk), ALLOCATABLE :: q(:) ! power source per volume

REAL(sdk), ALLOCATABLE :: fric(:) ! wall friction coefficient per volume

!

! Geometry arrays

!

REAL(sdk), ALLOCATABLE :: dx(:) ! length of each computational volume

REAL(sdk), ALLOCATABLE :: vol(:) ! fluid volume of each finite volume

REAL(sdk), ALLOCATABLE :: fa(:) ! fluid flow area at each volume edge

REAL(sdk), ALLOCATABLE :: hd(:) ! hydraulic diameter at each volume edge

!

! Pressure, temperature and velocity are independent variables. The

! TARGET attribute makes it easier to add perturbations from the solution

! arrays at the end of each iteration.

!

TARGET :: pN, TN, vN

CONTAINS

SUBROUTINE StartStep

IMPLICIT NONE

!

! Copy end-of-step variables from the last step into start-of-step

! variables for this time step

!

rho = rhoN

rhoe = rhoeN

p = pN

v = vN

e = eN

T = TN

END SUBROUTINE StartStep

SUBROUTINE InitFluid

USE Eos

IMPLICIT NONE

INTEGER(sik) :: i

!

! Initialize all dependent fluid state variables

!

DO i = LBOUND(p,DIM=1), UBOUND(p,DIM=1)

rhoN(i) = rhoLiq (TN(i), pN(i))

eN(i) = SpEnergy(TN(i), pN(i))

rhoeN(i) = rhoN(i)*eN(i)

ENDDO

!

END SUBROUTINE InitFluid

SUBROUTINE AllocFluidAr

USE ScalarDat

! Allocate state variable arrays, space is reserved at each end for

! storage of boundary information.

ALLOCATE(pN(1-nbnd:ncell+nbnd))

ALLOCATE(p(1-nbnd:ncell+nbnd))

ALLOCATE(TN(1-nbnd:ncell+nbnd))

ALLOCATE(T(1-nbnd:ncell+nbnd))

ALLOCATE(rho(1-nbnd:ncell+nbnd))

ALLOCATE(rhoN(1-nbnd:ncell+nbnd))

ALLOCATE(e(1-nbnd:ncell+nbnd))

ALLOCATE(eN(1-nbnd:ncell+nbnd))

ALLOCATE(rhoe(1-nbnd:ncell+nbnd))

ALLOCATE(rhoeN(1-nbnd:ncell+nbnd))

ALLOCATE(v(1-nbnd:ncell+nbnd+1))

ALLOCATE(vN(1-nbnd:ncell+nbnd+1))

ALLOCATE(q(1:ncell))

ALLOCATE(fric(1:ncell+1))

ALLOCATE(dx(1-nbnd:ncell+nbnd))

ALLOCATE(fa(1-nbnd:ncell+nbnd+1))

ALLOCATE(hd(1-nbnd:ncell+nbnd+1))

ALLOCATE(vol(1-nbnd:ncell+nbnd))

END SUBROUTINE AllocFluidAr

END MODULE FluidArrays

Array Data for Equation Solution

MODULE Matrix

USE IntrType

!

! Array storage for the linear system used in the Newton iteration

!

REAL(sdk), ALLOCATABLE :: a(:,:), x(:), b(:)

INTEGER(sik), ALLOCATABLE :: ipvt(:)

CONTAINS

SUBROUTINE AllocMat

!

! Allocate space for the linear system

!

USE ScalarDat

ALLOCATE ( a(nvar,nvar), x(nvar), b(nvar), ipvt(nvar) )

END SUBROUTINE AllocMat

SUBROUTINE Solve

! Apply a linear solver appropriate to the matrix

USE ScalarDat

USE LUsolve

!

IMPLICIT NONE

INTEGER(sik) :: info

CALL sgefat(a,nvar,nvar,ipvt,info)

CALL sgeslt(a,nvar,nvar,ipvt,b,0)

END SUBROUTINE Solve

END MODULE Matrix

Data Locators

MODULE Location

!

! This contains several types of locaters. The first gives locations and

! associated weighting used to perform averages or differentials of

! quantities. The second provides pointers to move information from

! the data structure associated with the full system of equations and

! unknowns, to the data structure associated with physical state variables

! and spatial locations. The third is a set of translation indices to

! provide the full system variable index for each independent state

! variable. The fourth lists independent variables associated with

USE Intrtype

IMPLICIT NONE

!=======================================================================

! Information for averages (and differences)

TYPE averageT

INTEGER(sik) :: n ! number of points used

INTEGER(sik), POINTER :: i(:) ! cell or face index

REAL(sdk), POINTER :: wf(:) ! weight factor

END TYPE

!

TYPE (averageT), ALLOCATABLE, TARGET:: fluxAv(:), dPdx(:), dV(:)

! fluxAv - information used to construct averages in cell edge

! mass and energy fluxes

! dPdx - information used to evaluate the pressure gradient

! in the momentum equation

! dV - information used to evaluate the velocity derivative

! in momentum transport term

!

!========================================================================

! Translation from system to physical variables

TYPE varPtrT

REAL(sdk), POINTER :: var

INTEGER(sik) :: i, itype

END TYPE

! %var - pointer to the variable storage in the fluid arrays

! %i - index in the fluid array

! %itype - variable type

! 1 - pressure

! 2 - temperature

! 3 - velocity

!

TYPE (varPtrT), ALLOCATABLE :: loc(:)

!

!=======================================================================

! Translation for physical to system variables

!

INTEGER(sik), ALLOCATABLE :: ivarP(:), ivarT(:), ivarV(:)

!

! ivarP - index in the linear solution array "b" for the change

! in the corresponding element of pN

! ivarT - index in the linear solution array "b" for the change

! in the corresponding element of TN

! ivarV - index in the linear solution array "b" for the change

! in the corresponding element of vN

!

!========================================================================

! List of active varibles in each equation

! Also can be considered as a list of column numbers

! in each matrix row with potential for non-zero elements

TYPE varListT

INTEGER(sik) :: n, iloc

INTEGER(sik), POINTER :: i(:)

INTEGER(sik) :: eqnType

END TYPE

! n - number of columns with potential for non-zero elements

! iloc - index of cell or face location for the equation

! i - list of column indices

! eqnType - Equation type associated with this row

! 1 - mass

! 2 - energy

! 3 - momentum

!

TYPE (varListT), ALLOCATABLE :: varList(:)

CONTAINS

SUBROUTINE SetLocAr

!

! Allocate and load indices in the locator arrays

!

! Note that I only record potential nonzero contributions to matricies

! Depending on direction of velocity, some coefficients marked as potentially

! non-zero in structures below may actually be zero.

!

USE ScalarDat

USE FluidArrays

IMPLICIT NONE

INTEGER(sik) :: i, icol, irow, ivar, iLB, iUB, j, jj, nv

!

! Lower and Upper Bounds of velocity (cell edge) arrays

ivLB = 1

ivUB = ncell+1

! Calculate the total number of unknowns

nvar = 3*ncell + 1 ! p & T in each volume, velocity at each edge

IF (lBC.EQ.'v') THEN

ivLB = 2 ! velocity fixed at left edge, do not evaluate

nvar = nvar - 1 ! momentum equation at face 1

ENDIF

IF (uBC.EQ.'v') THEN

ivUB = ncell ! velocity fixed at right edge, do not evaluate

nvar = nvar -1 ! momentum equation at face ncell+1

ENDIF

ALLOCATE(loc(nvar),varList(nvar))

ALLOCATE(fluxAv(ncell+1), dPdx(ivLB:ivUB), dV(ivLB:ivUB))

!

! Load tables used for weight factors

!

DO i = 1,ncell+1

! Number of volumes that might contribute to an edge average

fluxAv(i)%n = 2*nbnd

ALLOCATE (fluxAv(i)%i(fluxAv(i)%n))

ALLOCATE (fluxAv(i)%wf(fluxAv(i)%n))

! Load indices of volumes that might contribute to an edge average

DO j = -nbnd, nbnd-1

fluxAv(i)%i(j+nbnd+1) = i + j

ENDDO

ENDDO

DO i = ivLB,ivUB

! Number of volumes contributing to the pressure derivative

dPdx(i)%n = 2

ALLOCATE (dPdx(i)%i(dPdx(i)%n))

ALLOCATE (dPdx(i)%wf(dPdx(i)%n))

dPdx(i)%i(1) = i - 1 ! Indicies of volumes contributing

dPdx(i)%i(2) = i ! to the derivative

! Weighting factors to generate the derivative.

! Note that calculation here requires a fixed mesh.

dPdx(i)%wf(2) = 2.0_sdk/(dx(i)+dx(i-1))

dPdx(i)%wf(1) = - dPdx(i)%wf(2)

ENDDO

DO i = ivLB, ivUB

! Number of edges contributing to the velocity derivative

dV(i)%n = 1+ 2*nbnd

ALLOCATE (dV(i)%i(dV(i)%n))

ALLOCATE (dV(i)%wf(dV(i)%n))

! Load indices of edges that might contribute

DO j = -nbnd, nbnd

dV(i)%i(j+nbnd+1) = i + j

ENDDO

ENDDO

! ! Translation of state to system variables

iLB = LBOUND(p,DIM=1)

iUB = UBOUND(p,DIM=1)

ALLOCATE ( ivarP(iLB:iUB), ivarT(iLB:iUB))

iLB = LBOUND(v,DIM=1)

iUB = UBOUND(v,DIM=1)

ALLOCATE ( ivarV(iLB:iUB))

ivarP = 0

ivarT = 0

ivarV = 0

ivar = 0 ! variable index in the full system of equations

IF (ivLB.EQ.1) THEN

ivar = ivar+1

loc(ivar)%var => vN(1)

loc(ivar)%i = 1

loc(ivar)%itype = 3

ivarV(1) = ivar

ENDIF

DO i = 1,ncell

ivar = ivar+1 ! Pressure

loc(ivar)%var => pN(i)

loc(ivar)%i = i

loc(ivar)%itype = 1

ivarP(i) = ivar

ivar = ivar+1 ! Temperature

loc(ivar)%var => TN(i)

loc(ivar)%i = i

loc(ivar)%itype = 2

ivarT(i) = ivar

! No system variable at right edge

IF(i+1.GT.ivUB) EXIT ! if it has a constant velocity BC

ivar = ivar+1 ! Velocity

loc(ivar)%var => vN(i+1)

loc(ivar)%i = i+1

loc(ivar)%itype = 3

ivarV(i+1) = ivar

ENDDO

! Calculate and store indices of all variables used

! in each equation. This is a list of columns in

! each row of the matrix with potential for non-zero

! elements.

!

irow = 0

IF(ivarV(1).NE.0) THEN

irow = irow + 1

varList(irow)%n = 1 ! Velocity at this edge appears in the equation

varList(irow)%iloc = 1

varList(irow)%eqnType = 1

! ! Velocity at edge 2 may appear in the equation

IF( ivarV(i+1).NE.0) varList(irow)%n = varList(irow)%n + 1

varList(irow)%n = varList(irow)%n + 2*nbnd ! Variables in Dp & 1/rho terms

ALLOCATE(varList(irow)%i(varList(irow)%n))

j = 1

varList(irow)%i(j) = ivarV(1) ! Local Velocity

varList(irow)%i(j+1) = ivarT(1) ! Temp in 1/rho

varList(irow)%i(j+2) = ivarP(1) ! press in 1/rho and Grad P

IF( ivarV(2).NE.0) varList(irow)%i(j+3) = ivarV(2)

ENDIF

!

DO i = 1,ncell

irow = irow + 1 ! Mass Equation

varList(irow)%iloc = i

varList(irow)%eqnType = 1

varList(irow)%n = 2*( 2*nbnd+1) ! Temperatures and Pressures

IF( ivarV(i).NE.0) varList(irow)%n = varList(irow)%n + 1

IF( ivarV(i+1).NE.0) varList(irow)%n = varList(irow)%n + 1

ALLOCATE(varList(irow)%i(varList(irow)%n))

nv = 0

IF(ivarV(i).NE.0) THEN ! Velocity for left edge flux

nv = nv+1

varList(irow)%i(nv) = ivarV(i)

ENDIF

! Pressure and Temperatures in time derivative and

! Flux terms

DO j = -nbnd,nbnd

nv = nv+1

varList(irow)%i(nv) = ivarP(i+j)

nv = nv+1

varList(irow)%i(nv) = ivarT(i+j)

ENDDO

IF(ivarV(i+1).NE.0) THEN ! Velocity for right edge

nv = nv+1

varList(irow)%i(nv) = ivarV(i+1)

ENDIF

irow = irow + 1 ! Energy Equation

! Same footprint as Mass Equation

ALLOCATE(varList(irow)%i(varList(irow-1)%n))

varList(irow) = varList(irow-1)

varList(irow)%eqnType = 2

!

! Momentum Equation at the Right cell face

IF(ivarV(i+1).EQ.0) CYCLE

irow = irow + 1

varList(irow)%iloc = i+1

varList(irow)%eqnType = 3

varList(irow)%n = 1

IF( ivarV(i+2).NE.0) varList(irow)%n = varList(irow)%n + 1

IF( ivarV(i).NE.0) varList(irow)%n = varList(irow)%n + 1

varList(irow)%n = varList(irow)%n + 2*nbnd ! Variables in 1/rho term

IF(ivarP(i+1).NE.0) THEN ! Dp and 1/rho contributions to right

varList(irow)%n = varList(irow)%n + 2*nbnd

ENDIF

ALLOCATE(varList(irow)%i(varList(irow)%n))

varList(irow)%i = 0

nv = 0

IF( ivarV(i).NE.0) THEN

nv = nv +1

varList(irow)%i(nv) = ivarV(i)

ENDIF

nv = nv + 1

varList(irow)%i(nv) = ivarV(i+1) ! Local Velocity

DO j = -nbnd+1,nbnd ! 1/rho and Grad P terms

IF(ivarP(i+j).EQ.0) CYCLE

nv = nv + 1

varList(irow)%i(nv) = ivarT(i+j)

nv = nv + 1

varList(irow)%i(nv) = ivarP(i+j)

ENDDO

IF( ivarV(i+2).NE.0) varList(irow)%i(nv+1) = ivarV(i+2)

ENDDO

END SUBROUTINE SetLocAr

END MODULE Location

Set Initial Conditions

MODULE Input

USE Intrtype

USE ScalarDat

USE FluidArrays

!

! This module sets initial conditions directly that would normally be

! obtained from an input file

!

CONTAINS

SUBROUTINE SetIC

ncell = 12

dt = 0.1_sdk

CALL AllocFluidAr

DO i = LBOUND(v,DIM=1), UBOUND(v,DIM=1)

vN(i) = 5.0_sdk

fa(i) = 0.4_sdk

ENDDO

DO i = LBOUND(v,DIM=1), UBOUND(p,DIM=1)

pN(i) = 1.0e5_sdk

TN(i) = 300.0_sdk

dx(i) = 0.5_sdk

vol(i) = dx(i)*fa(i)

ENDDO

q = 6.98e8_sdk

fric = 0.1_sdk

lBC = 'v'

uBC = 'p'

END SUBROUTINE SetIC

END MODULE Input

Create a Subprogram for Each Flow Equation

MODULE FlowEqn

!

! This contains functions needed to evaluate flow equation residuals

!

USE IntrType

USE Location

USE FluidArrays

USE ScalarDat

Use Eos

CONTAINS Start With the Mass Conservation Equation

FUNCTION MassEqn(i)

!

! Residual of the mass equation at cell i

!

IMPLICIT NONE

!

! Input

!

! i - cell index at which the equation is evaluated

!

REAL(sdk) :: MassEqn

INTEGER(sik), INTENT(IN) :: i

MassEqn = (rhoN(i) - rho(i))*vol(i)/dt + Flux(rhoN,i+1) - &

Flux(rhoN,i)

END FUNCTION MassEqn

!

Energy Conservation

FUNCTION EnergyEqn(i)

!

! Residual of the energy equation at cell i

!

IMPLICIT NONE

!

! Input

!

! i - cell index at which the equation is evaluated

!

REAL(sdk) :: EnergyEqn

INTEGER(sik), INTENT(IN) :: i

EnergyEqn = (rhoeN(i)- rhoe(i))*vol(i)/dt + &

Flux(rhoeN,i+1) - Flux(rhoeN,i) &

+ pN(i)*(fa(i+1)*vN(i+1)-fa(i)*vN(i)) - q(i)*vol(i)

END FUNCTION EnergyEqn

Momentum Conservation

FUNCTION MomenEqn(i)

!

! Residual of the momentum equation at face i (low edge of cell i)

!

IMPLICIT NONE

!

! Input

!

! i - edge index at which the equation is evaluated

!

REAL(sdk) :: MomenEqn

INTEGER(sik), INTENT(IN) :: i

MomenEqn = (vN(i)- v(i))/dt + Vdv(i) + GradP(i)/EdgeAv(rhoN,i) + &

fric(i)*vN(i)*ABS(vN(i))

END FUNCTION MomenEqn

!

Average to Get Cell Edge Values for State Variables

FUNCTION EdgeAv(x,j)

!

! Calculate the value of scalar field x at cell face j based on

! weighting coefficients in fluxAv

!

IMPLICIT NONE

REAL(sdk) :: EdgeAv

! quantity being fluxed

REAL(sdk), INTENT(IN) :: x(LBOUND(rho,1):UBOUND(rho,1))

INTEGER(sik), INTENT(IN) :: j ! edge index

TYPE (averageT), POINTER :: a

INTEGER(sik) :: i

a => fluxAv(j)

EdgeAv = 0

DO i = 1,a%n

EdgeAv = EdgeAv + a%wf(i)*x(a%i(i))

ENDDO

END FUNCTION EdgeAv

!

Calculate Flux of Mass or Energy at Cell Edge

FUNCTION Flux(x,j)

!

! Calculate the flux of scalar field x at cell face j based on

! weighting coefficients in fluxAv

!

IMPLICIT NONE

REAL(sdk) :: Flux

! Quantity being fluxed

REAL(sdk), INTENT(IN) :: x(LBOUND(rho,1):UBOUND(rho,1))

INTEGER(sik), INTENT(IN) :: j ! edge index

TYPE (averageT), POINTER :: a

INTEGER(sik) :: i

a => fluxAv(j)

Flux = 0 ! First calculate edge value of x

DO i = 1,a%n

Flux = Flux + a%wf(i)*x(a%i(i))

ENDDO

Flux = Flux*vN(j)*fa(j) ! Multiply edge value by edge area and velocity

END FUNCTION Flux

!

Special Function for Momentum Transport

FUNCTION Vdv(j)

!

! Calculate the momentum flux term at face j based upon weighting

! information in derived type array dV

!

IMPLICIT NONE

REAL(sdk) :: Vdv

INTEGER(sik), INTENT(IN) :: j ! Face where momentum Equation is evaluated

INTEGER(sik) :: i

TYPE (averageT), POINTER :: a

a => dV(j)

Vdv = 0

DO i = 1,a%n

Vdv = Vdv + a%wf(i)*vN(a%i(i))

ENDDO

Vdv = vN(j)*Vdv

END FUNCTION Vdv

Special Function for Pressure Gradient

FUNCTION GradP(j)

!

! Calculate an approximation to the pressure gradient at face j

! based on information in the dPdx array

!

IMPLICIT NONE

REAL(sdk) :: GradP

TYPE (averageT), POINTER :: a

INTEGER(sik), INTENT(IN) :: j ! Face where momentum Equation is evaluated

INTEGER(sik) :: i

a => dPdx(j)

GradP = 0

DO i = 1,a%n

GradP = GradP + a%wf(i)*pN(a%i(i))

ENDDO

END FUNCTION GradP

END MODULE FlowEqn

Drive the Transient Solution

MODULE Trans

USE IntrType

!

! Subprograms used in the flow transient

!

CONTAINS

SUBROUTINE Transient

!

! Driver for the transient

!

USE ScalarDat

USE FluidArrays

USE Matrix

IMPLICIT NONE

time = 0.0_sdk

DO nstep = 1,nsmax

CALL StartStep

converged = .FALSE.

CALL SetWf

CALL Residuals(0)

DO it = 1,itmax

CALL SetWf

CALL Jacobian

CALL Solve

CALL EvalVar

CALL Residuals(1)

IF(converged) EXIT

ENDDO

IF(.NOT.converged) THEN

WRITE (*, '(a,f7.2)') 'Iteration failed at time ',time

STOP

ENDIF

time = time + dt

IF( time.GT. endtime-0.01*dt) EXIT

ENDDO

END SUBROUTINE Transient

SUBROUTINE EvalVar

USE FluidArrays

USE Location

USE Matrix

USE ScalarDat

USE Eos

IMPLICIT NONE

INTEGER(sik) :: i

!

! Update the state variable arrays

!

DO i = 1,nvar

loc(i)%var = loc(i)%var + b(i)

ENDDO

DO i = 1,ncell

rhoN(i) = rhoLiq (TN(i), pN(i))

eN(i) = SpEnergy(TN(i), pN(i))

rhoeN(i) = rhoN(i)*eN(i)

ENDDO

END SUBROUTINE EvalVar

SUBROUTINE Residuals(mode)

!

! Evaluate Mass, Energy, and Momentum equation residuals

!

! Input

!

! mode - 0 only evaluate residuals

! 1 find maximum scaled residual

!

USE Matrix

USE FluidArrays

USE Location

USE FlowEqn

USE ScalarDat

!

IMPLICIT NONE

INTEGER(sik), INTENT(IN) :: mode

INTEGER(sik) :: i, ivar

!

DO i = 1,nvar

SELECT CASE (varList(i)%eqnType)

CASE(1) ! Mass Equation

b(i) = MassEqn(varList(i)%iloc)

CASE(2) ! Energy Equation

b(i) = EnergyEqn(varList(i)%iloc)

CASE(3) ! Momentum Equation

b(i) = MomenEqn(varList(i)%iloc)

END SELECT

ENDDO

!

IF(mode.EQ.0) RETURN

!

ivar = 0

residMax = 0.0_sdk

IF(ivarV(1).NE.0) THEN

ivar = ivar+1

residMax = MAX(residMax,b(ivar)/MAX(0.01, v(1), vN(1))*dt)

ENDIF

DO i = 1,ncell

ivar = ivar+1

residMax = MAX(residMax,b(ivar)/(vol(i)*MAX(0.01, rho(i), rhoN(i)))*dt)

ivar = ivar+1

residMax = MAX(residMax,b(ivar)/(vol(i)*MAX(0.01, rhoe(i), rhoeN(i)))*dt)

IF(ivarV(i+1).EQ.0) CYCLE

ivar = ivar+1

residMax = MAX(residMax,b(ivar)/MAX(0.01, v(i+1), vN(i+1))*dt)

ENDDO

!

IF(residMax.LT.eps) converged = .TRUE.

!

END SUBROUTINE Residuals

SUBROUTINE Jacobian

!

! Evaluate the Jacobian Matrix using finite approximations to each

! partial derivative

!

USE Matrix

USE FluidArrays

USE ScalarDat

USE FlowEqn

USE Eos

USE Location

!

IMPLICIT NONE

INTEGER(sik) :: i, ii, ivar, j

REAL(sdk) :: fbase,varBase, rhoBase, rhoeBase, fracDvar=0.001_sdk, dvar

REAL(sdk) :: fpert

! First clear the Jacobian Matrix

a = 0.0_sdk

! Now the set non-zero Jacobian elements

! Loop over all rows in the matrix

DO i = 1,nvar

DO j = 1,varList(i)%n ! Loop over all non-zero columns in the row

ivar = varList(i)%i(j)

IF(ivar.EQ.0) CYCLE

varBase = loc(ivar)%var ! hold the old variable value

!

! Perturb the variable

!

IF(loc(ivar)%itype.EQ.3) THEN ! velocity

IF(loc(ivar)%var.GE.0) THEN

dvar = MAX(fracDvar, fracDvar*loc(ivar)%var)

loc(ivar)%var = loc(ivar)%var + dvar

ELSE

dvar = MIN(-fracDvar, fracDvar*loc(ivar)%var)

loc(ivar)%var = loc(ivar)%var + dvar

ENDIF

ELSE ! pressure or temperature

dvar = fracDvar*loc(ivar)%var

loc(ivar)%var = loc(ivar)%var + dvar

ii = loc(ivar)%i

rhoBase = rhoN(ii)

rhoeBase = rhoeN(ii)

rhoN(ii) = RhoLiq(TN(ii),pN(ii))

rhoeN(ii) = rhoN(ii)*SpEnergy(TN(ii),pN(ii))

ENDIF

! Reevaluate the function

SELECT CASE (varList(i)%eqnType)

CASE(1) ! Mass Equation

fpert = MassEqn(varList(i)%iloc)

CASE(2) ! Energy Equation

fpert = EnergyEqn(varList(i)%iloc)

CASE(3) ! Momentum Equation

fpert = MomenEqn(varList(i)%iloc)

END SELECT

! Jacobian

a(i,ivar) = (b(i) - fpert)/dvar

!

loc(ivar)%var = varBase ! Restore variables

IF(loc(ivar)%itype.NE.3) THEN

ii = loc(ivar)%i

rhoN(ii) = rhoBase

rhoeN(ii) = rhoeBase

ENDIF

ENDDO

ENDDO

END SUBROUTINE Jacobian

SUBROUTINE SetWf

!

! Set weighting factors for edge averages and differences

! Contents of this subroutine are specific to the 1st order upwind

! method used as a sample

!

USE FluidArrays

USE ScalarDat

USE Location

!

IMPLICIT NONE

INTEGER(sik) :: i

DO i = 1,ncell+1

IF(vN(i).LT.0.0_sdk) THEN

fluxAv(i)%wf(1) = 0.0_sdk

fluxAv(i)%wf(2) = 1.0_sdk

IF(i.LT.ivLB.OR.i.GT.ivUB) CYCLE

dV(i)%wf(1) = 0.0_sdk

dV(i)%wf(2) = dPdx(i)%wf(1)

dV(i)%wf(3) = dPdx(i)%wf(2)

ELSE

fluxAv(i)%wf(2) = 0.0_sdk

fluxAv(i)%wf(1) = 1.0_sdk

IF(i.LT.ivLB.OR.i.GT.ivUB) CYCLE

dV(i)%wf(1) = dPdx(i)%wf(1)

dV(i)%wf(2) = dPdx(i)%wf(2)

dV(i)%wf(3) = 0.0_sdk

ENDIF

ENDDO

END SUBROUTINE SetWf

END MODULE Trans

Equation of State

MODULE Eos

USE IntrType

!

! Equation of state information and other physical properties

!

REAL(sdk) :: mu = 0.00086_sdk ! water viscosity (Pa*s)

REAL(sdk) :: cond = 0.62_sdk ! water conductivity (w/m/K)

REAL(sdk) :: cv = 4186.0_sdk ! water specific heat (J/kg/K)

REAL(sdk) :: dRhoDt =-0.2538_sdk ! derivative of density with

! respect to temperature (kg/m**3/K)

CONTAINS

FUNCTION RhoLiq(T,p)

!

! Water density

!

! Input

! T - Temperature (K)

! p - pressure (Pa)

! Output

! RhoLiq - density (kg/m**3)

!

IMPLICIT NONE

REAL(sdk), INTENT(IN) :: T, p

REAL(sdk) :: RhoLiq

!

RhoLiq = 923.86 + dRhoDt*T

!

END FUNCTION RhoLiq

FUNCTION SpEnergy(T,p)

!

! Water Specific Internal Energy

!

! Input

! T - Temperature (K)

! p - pressure (Pa)

! Output

! SpEnergy - Specific Internal Energy (J/kg)

!

IMPLICIT NONE

REAL(sdk), INTENT(IN) :: T, p

REAL(sdk) :: SpEnergy

!

SpEnergy = cv*T

!

END FUNCTION SpEnergy

END MODULE Eos

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download