refactoring, put rk substeps into rk4
This commit is contained in:
ignis
parent
0849b31e7b
commit
fea70aa4cb
1 changed files with 47 additions and 60 deletions
107
m_fdm_calc.f90
107
m_fdm_calc.f90
|
|
@ -693,6 +693,27 @@ module m_fdm_calc
|
|||
END SUBROUTINE fns
|
||||
|
||||
|
||||
SUBROUTINE solve(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
integer :: i,j,k,xx,yy,zz
|
||||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
|
||||
|
||||
! advance reacting scalars - either Euler or Adams-Bashforth
|
||||
if (fors) then
|
||||
call EE1(xx,yy,zz,uu_,vv_,ww_,yy1,yy2)
|
||||
fors = .false.
|
||||
else
|
||||
call AB2(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
||||
yy2 = yyf
|
||||
end if
|
||||
|
||||
return
|
||||
END SUBROUTINE solve
|
||||
|
||||
|
||||
subroutine EE1(xx,yy,zz,uu_,vv_,ww_,yy1,rhs1)
|
||||
|
||||
implicit none
|
||||
|
|
@ -774,48 +795,7 @@ module m_fdm_calc
|
|||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
|
||||
|
||||
|
||||
CALL rotarysubstep(1, xx, yy, zz, uu_, vv_, ww_, yy1, yy1, yy2)
|
||||
CALL rotarysubstep(2, xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
|
||||
CALL rotarysubstep(3, xx, yy, zz, uu_, vv_, ww_, yy2, yyf, yy1)
|
||||
CALL rotarysubstep(4, xx, yy, zz, uu_, vv_, ww_, yyf, yy1, yy2)
|
||||
CALL lastsubstep (xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
|
||||
|
||||
return
|
||||
END SUBROUTINE RK4
|
||||
|
||||
SUBROUTINE solve(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
integer :: i,j,k,xx,yy,zz
|
||||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
|
||||
|
||||
! advance reacting scalars - either Euler or Adams-Bashforth
|
||||
if (fors) then
|
||||
call EE1(xx,yy,zz,uu_,vv_,ww_,yy1,yy2)
|
||||
fors = .false.
|
||||
else
|
||||
call AB2(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
||||
yy2 = yyf
|
||||
end if
|
||||
|
||||
return
|
||||
END SUBROUTINE solve
|
||||
|
||||
subroutine rotarysubstep(istage,xx,yy,zz,uu_,vv_,ww_,r0,r1,r2)
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: xx,yy,zz,istage
|
||||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
real*8 :: r0(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq)
|
||||
|
||||
real*8 :: a(5),b(5)
|
||||
real*8 :: at, bt
|
||||
|
||||
integer :: i,j,k,l
|
||||
|
||||
a(1)= 970286171893./4311952581923.
|
||||
a(2)= 6584761158862./12103376702013.
|
||||
|
|
@ -829,6 +809,29 @@ module m_fdm_calc
|
|||
b(4)= 2114624349019./3568978502595.
|
||||
b(5)= 5198255086312./14908931495163.
|
||||
|
||||
|
||||
CALL rotarysubstep(1, xx, yy, zz, uu_, vv_, ww_, yy1, yy1, yy2)
|
||||
CALL rotarysubstep(2, xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
|
||||
CALL rotarysubstep(3, xx, yy, zz, uu_, vv_, ww_, yy2, yyf, yy1)
|
||||
CALL rotarysubstep(4, xx, yy, zz, uu_, vv_, ww_, yyf, yy1, yy2)
|
||||
CALL lastsubstep (xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
|
||||
|
||||
return
|
||||
|
||||
contains
|
||||
|
||||
subroutine rotarysubstep(istage,xx,yy,zz,uu_,vv_,ww_,r0,r1,r2)
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: xx,yy,zz,istage
|
||||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
real*8 :: r0(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq)
|
||||
|
||||
real*8 :: at, bt
|
||||
|
||||
integer :: i,j,k,l
|
||||
|
||||
CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_)
|
||||
|
||||
at=a(istage)*fdmdt
|
||||
|
|
@ -860,15 +863,11 @@ module m_fdm_calc
|
|||
real*8 :: r0(xx,yy,zz,neq), r1(xx,yy,zz,neq), r2(xx,yy,zz,neq)
|
||||
real*8 :: uu_(xx,yy,zz), vv_(xx,yy,zz), ww_(xx,yy,zz)
|
||||
|
||||
real*8 :: b
|
||||
|
||||
integer :: i, j, k
|
||||
|
||||
b = 5198255086312./14908931495163.
|
||||
|
||||
CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_)
|
||||
|
||||
bt = b*fdmdt
|
||||
bt = b(5)*fdmdt
|
||||
|
||||
r0 = r1 + bt*r2
|
||||
|
||||
|
|
@ -901,23 +900,9 @@ module m_fdm_calc
|
|||
integer :: i,j,k,xx,yy,zz,istage
|
||||
real*8 :: at,bt , wrate , yr
|
||||
real*8 :: ri(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq),f(xx,yy,zz,neq)
|
||||
real*8 :: a(5),b(5)
|
||||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
integer :: nfinal, iscr, mspec, mpict, msave, nmindt, nv
|
||||
|
||||
|
||||
a(1)= 970286171893./4311952581923.
|
||||
a(2)= 6584761158862./12103376702013.
|
||||
a(3)= 2251764453980./15575788980749.
|
||||
a(4)= 26877169314380./34165994151039.
|
||||
a(5)=0.
|
||||
|
||||
b(1)= 1153189308089./22510343858157.
|
||||
b(2)= 1772645290293./4653164025191.
|
||||
b(3)= -1672844663538./4480602732383.
|
||||
b(4)= 2114624349019./3568978502595.
|
||||
b(5)= 5198255086312./14908931495163.
|
||||
|
||||
CALL fns(ri,f,xx,yy,zz,uu_,vv_,ww_)
|
||||
|
||||
IF(istage<5) THEN
|
||||
|
|
@ -957,6 +942,8 @@ module m_fdm_calc
|
|||
return
|
||||
END SUBROUTINE substep
|
||||
|
||||
END SUBROUTINE RK4
|
||||
|
||||
subroutine tp2 (a, b, n1, n2)
|
||||
! a = transpose(b)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue