-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathtimestep.f90
More file actions
130 lines (128 loc) · 6.32 KB
/
timestep.f90
File metadata and controls
130 lines (128 loc) · 6.32 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
! TODO: add fail flag
subroutine timeIntegration(q,p,fail,I2E,B2E,In,Bn,qnrm,Jinv,Jinv2,detJ,detJ2,Minv,w,phi,gphi,w1,rBC,resids,&
phiL,phiR,qlist,resnorm,gamma,Rgas,CFL,convtol,min_iter,max_iter,nelem,niface,nbface,nqelem,Ng,Ng1)
! -----------------------------------------------------------------------
! Purpose: use forward Euler to timestep the governing equations
!
! Inputs:
! q[nelem,4] = the initial state
! I2E[niface,4]
! B2E[nbface,3]
! In[niface,3] = normal vector for interior faces, and length
! Bn[nbface,3] = normal vector for boundary faces, and length
! rBC[5] = information on freestream/boundary values
! gamma = heat ratio
! Rgas = specific gas constant
! dt = the delta t for each timestep
! ndt = number of time steps
!
! Outs:
! q[nelem,4] = the state after ndt timesteps
! resids[nelem,4] = the residual after ndt timesteps
!
! -----------------------------------------------------------------------
implicit none
integer, intent(in) :: nelem,niface,nbface,min_iter,max_iter,p,nqelem,Ng,Ng1
real(8), intent(inout), dimension(nelem,(p+1)*(p+2)/2,4) :: q
integer, intent(in), dimension(niface,4) :: I2E
integer, intent(in), dimension(nbface,3) :: B2E
real(8), intent(in), dimension(niface,3) :: In ! includes length as 3rd component
real(8), intent(in), dimension(nbface,3) :: Bn ! includes length as 3rd component
real(8), intent(in), dimension(nqelem,Ng1,2) :: qnrm ! nrm for high-q elems
real(8), intent(in), dimension(nelem) :: detJ
real(8), intent(in), dimension(nqelem,Ng) :: detJ2
real(8), intent(in), dimension(nelem,2,2) :: Jinv
real(8), intent(in), dimension(nqelem,Ng,2,2) :: Jinv2
real(8), intent(in), dimension(Ng,(p+1)*(p+2)/2) :: phi
real(8), intent(in), dimension(Ng,(p+1)*(p+2)/2,2) :: gphi
real(8), intent(in), dimension(Ng) :: w
real(8), intent(in), dimension(Ng1) :: w1
real(8), intent(in), dimension(3,Ng1,(p+1)*(p+2)/2) :: phiL,phiR
integer, intent(in), dimension(nqelem) :: qlist
real(8), intent(in), dimension(nelem,(p+1)*(p+2)/2,(p+1)*(p+2)/2) :: Minv
real(8), intent(out),dimension(nelem,(p+1)*(p+2)/2,4) :: resids
real(8), intent(in), dimension(5) :: rBC
real(8), intent(in) :: gamma, Rgas, convtol, CFL
real(8), intent(out), dimension(max_iter) :: resnorm
logical, intent(out) :: fail
!f2py intent(in) node,E2N,I2E,B2E,In,Bn,gamma,rBC,Jinv,Minv
!f2py intent(out) resids, resnorm
!f2py intent(in,out) q
integer :: ielem, iter
integer, dimension(3) :: loc
real(8), dimension(nelem) :: wavespeed
real(8), dimension(nelem,(p+1)*(p+2)/2,4) :: q1,q2
real(8), dimension((p+1)*(p+2)/2,4) :: new_resids
real(8) :: dt
resnorm(:) = -1 ! set to high value to allow first pass of while loop
fail = .false.
do iter = 1,max_iter
call getResidual(q,p,I2E,B2E,In,Bn,qnrm,rBC,resids,Jinv,Jinv2,detJ,detJ2,w,phi,gphi,w1,&
phiL,phiR,qlist,wavespeed,gamma,Rgas,nelem,niface,nbface,nqelem,Ng,Ng1)
resnorm(iter) = maxval(resids)
loc = maxloc(resids)
if (mod(iter,2000) == 0 .or. iter == 1) then
print*, " ---------------------------------------------------------------------"
print*, " iter maximum residual element basis q"
print*, " ---------------------------------------------------------------------"
endif
if (mod(iter,100) == 0) then
print*, iter, resnorm(iter), loc(1), loc(2), loc(3)
endif
do ielem = 1,nelem
call applyMatRes(Minv(ielem,:,:),resids(ielem,:,:),new_resids,p)
dt = CFL*detJ(ielem)/wavespeed(ielem) ! we let detJ = 2*area
q1(ielem,:,:) = q(ielem,:,:) - dt*new_resids ! the area term cancels with the dt expression
enddo
if (p == 0) then
q = q1
elseif (p == 1) then
call getResidual(q1,p,I2E,B2E,In,Bn,qnrm,rBC,resids,Jinv,Jinv2,detJ,detJ2,w,phi,gphi,w1,&
phiL,phiR,qlist,wavespeed,gamma,Rgas,nelem,niface,nbface,nqelem,Ng,Ng1)
resnorm(iter) = maxval(resids)
loc = maxloc(resids)
do ielem = 1,nelem
call applyMatRes(Minv(ielem,:,:),resids(ielem,:,:),new_resids,p)
dt = CFL*detJ(ielem)/wavespeed(ielem)
q(ielem,:,:) = 0.5d0*(q(ielem,:,:) + q1(ielem,:,:) - dt*new_resids)
enddo
elseif (p == 2) then
call getResidual(q1,p,I2E,B2E,In,Bn,qnrm,rBC,resids,Jinv,Jinv2,detJ,detJ2,w,phi,gphi,w1,&
phiL,phiR,qlist,wavespeed,gamma,Rgas,nelem,niface,nbface,nqelem,Ng,Ng1)
do ielem = 1,nelem
call applyMatRes(Minv(ielem,:,:),resids(ielem,:,:),new_resids,p)
dt = CFL*detJ(ielem)/wavespeed(ielem)
q2(ielem,:,:) = 0.25d0*(3.d0*q(ielem,:,:) + q1(ielem,:,:) - dt*new_resids)
enddo
call getResidual(q2,p,I2E,B2E,In,Bn,qnrm,rBC,resids,Jinv,Jinv2,detJ,detJ2,w,phi,gphi,w1,&
phiL,phiR,qlist,wavespeed,gamma,Rgas,nelem,niface,nbface,nqelem,Ng,Ng1)
resnorm(iter) = maxval(resids)
loc = maxloc(resids)
do ielem = 1,nelem
call applyMatRes(Minv(ielem,:,:),resids(ielem,:,:),new_resids,p)
dt = CFL*detJ(ielem)/wavespeed(ielem)
q(ielem,:,:) = 1.d0/3.d0*(q(ielem,:,:) + 2.d0*q2(ielem,:,:) - 2.d0*dt*new_resids)
enddo
endif
if (iter > min_iter .and. resnorm(iter) < convtol .and. resnorm(iter) == resnorm(iter))then
fail = .false.
exit
elseif (resnorm(iter) /= resnorm(iter)) then
fail = .true.
exit
endif
enddo
end subroutine timeIntegration
subroutine applyMatRes(Minv,resids,new_resids,p)
! we use matmul instead of DGEMM because the sizes are small enough
! that matmul is actually faster
implicit none
integer, intent(in) :: p
real(8), intent(in), dimension((p+1)*(p+2)/2,(p+1)*(p+2)/2) :: Minv
real(8), intent(in),dimension((p+1)*(p+2)/2,4) :: resids
real(8), intent(out),dimension((p+1)*(p+2)/2,4) :: new_resids
integer :: Nb
Nb = (p+1)*(p+2)/2
! call dgemm('N','N',Nb,4,Nb,1.d0,Minv,Nb,resids,Nb,0.d0,new_resids,Nb)
new_resids = matmul(Minv,resids)
end subroutine applyMatRes