c Main loop goes here
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine runSim()
implicit none
c User-modifiable parameters come first
c Regular variables come later
integer nX,nY,nMax
real*8 eps, tStep,
* xMin, xMax,
* yMin, yMax
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c User-modifiable parameters are eitehr here
c or in the file iniVals.f
integer bufSize !number of past timesteps to buffer
c The buffer of past non-linear terms for
c the adams bashforth method works like a
c "ring buffer" Essentially the current
c head of the buffer is tracked by an
c integer, "bufIdx", which is taken modulo
c bufSize. Previous entries are bufIdx-1 mod bufsize, etc.
c Obviously this works better with an index that starts at zero
parameter (nX=400)
parameter (nY=100)
parameter (nMax=400)
parameter (bufSize=2) ! Size of ring buffer
parameter (xMin=0.d0)
parameter (xMax=4.d0)
parameter (yMin=0.d0)
parameter (yMax=1.d0)
parameter (eps=1.d-20)
c xStep approx 1/100 = 0.01
c yStep approx the same
c xStep*xStep approx 0.0001
c Set tStep to a quarter of that (2.5d-5)
parameter (tStep = 2.5d-5) ! restore to 2.5d-5
c Thats a really small timestep!
integer maxIter, printFreq
c parameter(maxIter=1000000) ! Thats a lotta iters!
parameter(maxIter=9000) ! For testing
parameter(printFreq=5)
c End of user-modifiable stuff
ccccccccccccccccccccccccccccccccccccc
real*8 u, uv, uu, w, psi
real*8 dwdtNL, dwdtLin
real*8 x,y
dimension dwdtNL(0:nMax+1, 0:nMax+1,
* 0:bufSize-1) ! non-linear part of deriv
dimension dwdtLin(0:nMax+1, 0:nMax+1,
* 0:bufSize-1) ! linear part of deriv
dimension x(0:nX+1), y(0:nY+1)
dimension psi(0:nMax+1,0:nMax+1)
dimension w(0:nMax+1, 0:nMax+1)
dimension u(2,0:nMax+1,0:nMax+1)
integer it, bufIdx, nPrintCnt
c bufIdx is the index into past deriv buffers
real*8 currTime
real*8 vMax,vMin
call step(x, xMin, xMax, nx, .TRUE.)
call step(y, yMin, yMax, ny, .FALSE.)
currTime = 0.d0
nPrintCnt = 0
c Set initial conditions
c
call initializeFlow(
* u, nX, nY,
* nMax, x,
* y)
call printUVresult(x,y,u,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
nPrintCnt = nPrintCnt+1
call wFromU(w, u, nX, nY,
* nMax, x(1)-x(0), y(1)-y(0), eps)
call safePsiFromW(psi,
* w, nX, nY, nMax,
* x(1)-x(0), y(1)-y(0), eps)
call printPsi(x,y,psi,nx,ny,nmax,
* currTime,0,
* .FALSE., .FALSE.)
call uFromPsi(u, psi, nX, nY,
* nMax, x(1)-x(0), y(1)-y(0), eps)
call printUVresult(x,y,u,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
call printPsi(x,y,psi,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
nPrintCnt = nPrintCnt+1
do it=0,bufSize
bufIdx = mod(it,bufSize)
write (*,*) "Starting an RK step iter ",it
call DoRkStep(x,y,psi,w,
* dwdtLin, dwdtNL, tStep,
* nX,nY,nMax,
* bufIdx, bufSize, eps)
currTime = currTime+tStep
enddo
call safePsiFromW(psi,
* w, nX, nY, nMax,
* x(1)-x(0), y(1)-y(0), eps)
call uFromPsi(u, psi, nX, nY,
* nMax, x(1)-x(0), y(1)-y(0), eps)
call printUVresult(x,y,u,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
call printPsi(x,y,psi,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
nPrintCnt = nPrintCnt+1
c pause "Done RK Steps"
do it=bufSize+1,maxIter
bufIdx = mod(it,bufSize)
c SIAB stands for
c Semi Implicit Adams Bashforth
c
c I don't know if its a standard acronym
c but DoSemiImplicitAdamsBashforthStep
c is an awfully long function name
call DoSIABStep(x,y,psi,w,
* dwdtLin, dwdtNL, tStep,
* nX,nY,nMax,
* bufIdx, bufSize, eps)
call getUvMinMax(u,nX,nY,nMax,
* vMin, vMax,2)
write (*,*) "on iteration ", it,
* " Min/max are ", vMin, vMax
if(mod(it,printFreq).EQ.0) then
call safePsiFromW(psi,
* w, nX, nY, nMax,
* x(1)-x(0), y(1)-y(0), eps)
call uFromPsi(u, psi, nX, nY,
* nMax, x(1)-x(0), y(1)-y(0), eps)
call printUVresult(x,y,u,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
call printPsi(x,y,psi,nx,ny,nmax,
* currTime,nPrintCnt,
* .FALSE., .FALSE.)
nPrintCnt = nPrintCnt+1
endif
currTime = currTime+tStep
enddo
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c Steps the system once using RK4
c Used to jump-start semi-implicit Adams Bashforth
c
c stands for DoRkStep as in
c Do Runge Kutta step
c not as in "DORKstep"
c
cccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine DoRkStep(x,y,psi,w,
* linDrvs, nlDrvs, dt,
* nX,nY,nMax,
* nBufIdx, nBufSize, eps)
implicit none
c linDrvs means "LINear components of DeRiVativeS"
c nlDrvs means "NonLinear components of DeRiVativeS"
real*8 x,y,psi,w
real*8 linDrvs, nlDrvs
real*8 dt
real*8 tmpArr
real*8 eps
integer nX, nY, nMax,
* nBufIdx, nBufSize
dimension psi(0:nMax+1,0:nMax+1)
dimension w(0:nMax+1,0:nMax+1)
dimension linDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
dimension nlDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
dimension x(0:nX+1), y(0:nY+1)
dimension tmpArr(0:nMax+1,0:nMax+1)
real*8 deltaX, deltaY
deltaX = x(1)-x(0)
deltaY = y(1)-y(0)
call psiFromW(psi,
* w, tmpArr, nX, nY, nMax,
* deltaX, deltaY, eps)
call GetLinTerm(
* psi, w, x,y,
* nX, nY, nMax,
* nMax,
* linDrvs(0,0,nBufIdx))
call GetNonLinTerm(
* psi, w, x,y,
* nX, nY, nMax,
* nMax,
* nlDrvs(0,0,nBufIdx))
c For now we're just going to Euler-step these ones
c Hopefully Euler-stepping for the first couple
c steps is Okay for jump-starting Adams Bashforth
call makeLinComb(
* linDrvs(0,0,nBufIdx),
* nlDrvs(0,0,nBufIdx),tmpArr,
* nX,nY,nMax,nMax,nMax,
* dt,dt, .FALSE.)
call makeLinComb(
* w,
* tmpArr,w,
* nX,nY,nMax,nMax,nMax,
* 1.d0,1.d0, .FALSE.)
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c Steps the system using
c semi-implicit Adams Bashforth
c
cccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine DoSIABStep(x,y,psi,w,
* linDrvs, nlDrvs, dt,
* nX,nY,nMax,
* nBufIdx, nBufSize, eps)
implicit none
c linDrvs means "LINear components of DeRiVativeS"
c nlDrvs means "NonLinear components of DeRiVativeS"
real*8 x,y,psi,w
real*8 linDrvs, nlDrvs
real*8 dt
real*8 abCoeff
real*8 tmpArr
real*8 diag, diagx, diagy
real*8 eps, deltax, deltay
real*8 getViscosity
dimension abCoeff(0:2) !Adams Bashforth coeffs
integer nX, nY, nMax,
* nBufIdx, nBufSize,
* nBTmp ! nBTmp is used in place
c of nBufIdx in places where we want to subtract
c and mod
dimension psi(0:nMax+1,0:nMax+1)
dimension w(0:nMax+1,0:nMax+1)
dimension linDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
dimension nlDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
dimension tmpArr(0:nMax+1, 0:nMax+1)
dimension x(0:nX+1), y(0:nY+1)
deltaX = x(1)-x(0)
deltaY = y(1)-y(0)
call psiFromW(psi,
* w, tmpArr, nX, nY, nMax,
* deltaX, deltaY, eps)
nBTmp = nBufIdx+nBufSize
abCoeff(0) = 23.d0/12.d0
abCoeff(1) = -16.d0/12.d0
abCoeff(2) = 5.d0/12.d0
call GetLinTerm(
* psi, w, x,y,
* nX, nY, nMax,
* nMax,
* linDrvs(0,0,nBufIdx))
call GetNonLinTerm(
* psi, w, x,y,
* nX, nY, nMax,
* nMax,
* nlDrvs(0,0,nBufIdx))
call mul_and_copy(w,
* tmpArr, 1.d0, nX, nY,
* nMax, nMax,
* .TRUE.)
c Now the boundaries on tmpArr are correct
c also tmpArr =w
call makeLinComb(
* tmpArr,
* nlDrvs(0,0,mod(nBTmp,nBufSize)),
* tmpArr,
* nX,nY,nMax,nMax,nMax,
* 1.d0,
* dt*abCoeff(0), .FALSE.)
call makeLinComb(
* tmpArr,
* nlDrvs(0,0,mod(nBTmp-1,nBufSize)),
* tmpArr,
* nX,nY,nMax,nMax,nMax,
* 1.d0,
* dt*abCoeff(1), .FALSE.)
call makeLinComb(
* tmpArr,
* nlDrvs(0,0,mod(nBTmp-2,nBufSize)),
* tmpArr,
* nX,nY,nMax,nMax,nMax,
* 1.d0,
* dt*abCoeff(2), .FALSE.)
c Now we have the non-linear term,
c it is time for the semi-implicit part
c First add half deltaT times
c the linear term in the derivative
call makeLinComb(
* tmpArr,
* linDrvs(0,0,mod(nBTmp,nBufSize)),
* tmpArr,
* nX,nY,nMax,nMax,nMax,
* 1.d0,
* dt*5.d-1, .FALSE.)
c Now we have to use the conjugate-gradient
c solver to find the next w
diagx = -dt*5.d-1/(deltaX*deltaX)
diagy = -dt*5.d-1/(deltaY*deltaY)
diagx = diagx*getViscosity()
diagy = diagy*getViscosity()
diag = 1.d0-2.d0*diagx-2.d0*diagy
write (*,*) "Calling conjgrad "
write (*,*) "dt is ", dt
call conjgrad_2d_lat(diag,diagx,diagy,
* w, tmpArr,
* nX, nY, nMax, eps, .FALSE.)
c subroutine conjgrad_2d_lat(diag,xdiag,ydiag,x0,b,
c * nX, nY, nMax, eps, bYBoundaries)
return
end