A generalized scheme to updat the solution explicitly using any RK method and even first order euler explicit.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(inout), | dimension(-2:imx+2,-2:jmx+2,-2:kmx+2,1:n_var), target | :: | qp | Store primitive variable at cell center |
|
| real(kind=wp), | intent(in), | dimension(:, :, :, :) | :: | residue | Store residue at each cell-center |
|
| real(kind=wp), | intent(in), | dimension(1:imx-1,1:jmx-1,1:kmx-1) | :: | delta_t | Local time increment value at each cell center |
|
| type(celltype), | intent(in), | dimension(-2:imx+2,-2:jmx+2,-2:kmx+2) | :: | cells | Cell center quantities: volume |
|
| type(schemetype), | intent(in) | :: | scheme | finite-volume Schemes |
||
| type(flowtype), | intent(in) | :: | flow | Information about fluid flow: freestream-speed, ref-viscosity,etc. |
||
| character(len=*), | intent(in) | :: | type | |||
| real(kind=wp), | intent(in), | optional | :: | time_factor | ||
| real(kind=wp), | intent(in), | optional | :: | store_factor | ||
| logical, | intent(in), | optional | :: | use | ||
| real(kind=wp), | intent(inout), | optional | dimension(1:imx-1,1:jmx-1,1:kmx-1,1:n_var) | :: | Rn | |
| real(kind=wp), | intent(in), | optional | dimension(-2:imx+2,-2:jmx+2,-2:kmx+2,1:n_var), target | :: | un |
subroutine update_with(qp, residue, delta_t, cells, scheme, flow, type, time_factor, store_factor, use, Rn, un)
!< A generalized scheme to updat the solution explicitly using
!< any RK method and even first order euler explicit.
implicit none
real(wp), dimension(-2:imx+2,-2:jmx+2,-2:kmx+2,1:n_var), intent(inout), target:: qp
!< Store primitive variable at cell center
type(celltype), dimension(-2:imx+2,-2:jmx+2,-2:kmx+2), intent(in) :: cells
!< Cell center quantities: volume
real(wp), dimension(1:imx-1,1:jmx-1,1:kmx-1), intent(in) :: delta_t
!< Local time increment value at each cell center
real(wp), dimension(:, :, :, :), intent(in) :: residue
!< Store residue at each cell-center
type(schemetype), intent(in) :: scheme
!< finite-volume Schemes
type(flowtype), intent(in) :: flow
!< Information about fluid flow: freestream-speed, ref-viscosity,etc.
character(len=*), intent(in) :: type
real(wp), intent(in), optional :: time_factor ! time factor
real(wp), intent(in), optional :: store_factor
logical, intent(in), optional :: use
real(wp), dimension(-2:imx+2,-2:jmx+2,-2:kmx+2,1:n_var), intent(in), optional, target :: un
real(wp), dimension(1:imx-1,1:jmx-1,1:kmx-1,1:n_var), intent(inout), optional :: Rn
real(wp) :: TF = 1.0 !time factor
real(wp) :: SF = 1.0!store factor
Logical :: TU = .FALSE. !to use or not
real(wp), dimension(:,:,:,:), pointer :: Quse
integer :: i,j,k
real(wp) :: KE=0.
real(wp) :: beta
!sa variables
real(wp) :: vort
real(wp) :: fv1
real(wp) :: fv2
real(wp) :: fw
real(wp) :: g
real(wp) :: scap
real(wp) :: rsa
real(wp) :: kd2
real(wp) :: xi
real(wp) :: mass_residue
real(wp) :: x_mom_residue, y_mom_residue, z_mom_residue
real(wp) :: energy_residue
real(wp) :: TKE_residue, Omega_residue, kl_residue
if(present(time_factor)) TF=time_factor
if(present(store_factor)) SF=store_factor
if(present(use)) TU=use
!check if user want to update from particular solution
if(present(un))then
Quse(-2:imx+2,-2:jmx+2,-2:kmx+2,1:n_var)=>un(:,:,:,:)
else
Quse(-2:imx+2,-2:jmx+2,-2:kmx+2,1:n_var)=>qp(:,:,:,:)
end if
select case(type)
case('primitive')
!update primitive variable
do k = 1,kmx-1
do j = 1,jmx-1
do i = 1,imx-1
mass_residue = residue(i,j,k,1)
x_mom_residue = residue(i,j,k,2)
y_mom_residue = residue(i,j,k,3)
z_mom_residue = residue(i,j,k,4)
energy_residue = residue(i,j,k,5)
u1(1:n_var) = Quse(i,j,k,1:n_var)
! finding primitive residue
R(1) = mass_residue
R(2) = -1*(u1(2)/u1(1))*mass_residue + x_mom_residue/u1(1)
R(3) = -1*(u1(3)/u1(1))*mass_residue + y_mom_residue/u1(1)
R(4) = -1*(u1(4)/u1(1))*mass_residue + z_mom_residue/u1(1)
R(5) = 0.5*(flow%gm-1.)*(sum(u1(2:4)**2)*mass_residue) &
-(flow%gm-1.)*u1(2)*x_mom_residue &
-(flow%gm-1.)*u1(3)*y_mom_residue &
-(flow%gm-1.)*u1(4)*z_mom_residue &
+(flow%gm-1.)*energy_residue
select case(scheme%turbulence)
case('none')
!do nothing
continue
case('sst', 'sst2003')
TKE_residue = residue(i,j,k,6)
omega_residue = residue(i,j,k,7)
beta = beta1*sst_F1(i,j,k) + (1. - sst_F1(i,j,k))*beta2
R(5) = R(5) - (flow%gm-1.)*TKE_residue
R(6) = -(u1(6)/u1(1))*mass_residue&
+(1./(1.+bstar*u1(6)*delta_t(i,j,k)))*TKE_residue/u1(1)
R(7) = -(u1(7)/u1(1))*mass_residue&
+(1./(1.+2.*beta*u1(6)*delta_t(i,j,k)))*omega_residue/u1(1)
case('kkl')
TKE_residue = residue(i,j,k,6)
kl_residue = residue(i,j,k,7)
eta = u1(1)*dist(i,j,k)*(sqrt(0.3*u1(6))/(20*mu(i,j,k)))
fphi = (1+cd1*eta)/(1+eta**4)
R(5) = R(5) - (flow%gm-1.)*TKE_residue
R(6) = -(u1(6)/u1(1))*mass_residue&
+ (1./(1.+((2.5*((cmu**0.75)*u1(1)*(u1(6)**1.5)/max(u1(7),1.e-20))&
+(2*mu(i,j,k)/dist(i,j,k)**2))*delta_t(i,j,k))))*TKE_residue/u1(1)
R(7) = -(u1(7)/u1(1))*mass_residue&
+(1./(1.+(6*mu(i,j,k)*fphi/dist(i,j,k)**2)*delta_t(i,j,k)))*kl_residue/u1(1)
case DEFAULT
Fatal_error
end select
!check if user want to store residue
if(present(Rn)) then
Rn(i,j,k,1:n_var) = Rn(i,j,k,1:n_var) + SF*R(1:n_var)
if(TU) R(:) = Rn(i,j,k,:)
end if
!update
u2(:) = u1(:) - R(:)*(TF*delta_t(i,j,k)/cells(i,j,k)%volume)
!check solution for non pyhysical results
if((u2(1) < 0.) .or. (u2(5)) < 0.)then
Fatal_error
else !update
qp(i,j,k,1:5) = u2(1:5)
select case(trim(scheme%turbulence))
case('sst', 'sst2003', 'kkl')
if(u2(6)>0.) qp(i,j,k,6) = u2(6)
if(u2(7)>0.) qp(i,j,k,7) = u2(7)
case DEFAULT
! do nothing
continue
end select
end if
end do
end do
end do
case('conservative')
!include "update_conservative.inc"
!update conservative variable
do k = 1,kmx-1
do j = 1,jmx-1
do i = 1,imx-1
! getting conservative variable
u1(1) = Quse(i,j,k,1)
u1(2:) = Quse(i,j,k,2:)*u1(1)
select case(scheme%turbulence)
case('sst', 'sst2003', 'kkl')
KE = 0.0!u1(6)
case('sa','saBC')
KE=0.0
case DEFAULT
KE = 0.
end select
u1(5) = (u1(5)/(flow%gm-1.) + 0.5*sum(u1(2:4)**2))/u1(1) + KE
! get R
R(1:n_var) = residue(i,j,k,1:n_var)
! point implicit destruction term
select case(trim(scheme%turbulence))
case('none')
!do nothing
continue
case('sst', 'sst2003')
beta = beta1*sst_F1(i,j,k) + (1. - sst_F1(i,j,k))*beta2
R(6) = R(6)/(1+(beta*qp(i,j,k,7)*delta_t(i,j,k)))
R(7) = R(7)/(1+(2*beta*qp(i,j,k,7)*delta_t(i,j,k)))
case('kkl')
eta = u1(1)*dist(i,j,k)*(sqrt(0.3*u1(6))/(20*mu(i,j,k)))
fphi = (1+cd1*eta)/(1+eta**4)
R(6) = R(6)/(1.+((2.5*((cmu**0.75)*sqrt(u1(1))*(u1(6)**1.5)/max(u1(7),1.e-20))&
+(2*mu(i,j,k)/(dist(i,j,k)**2)))*delta_t(i,j,k)))
R(7) = R(7)/(1.+(6*mu(i,j,k)*fphi/(dist(i,j,k)**2))*delta_t(i,j,k))
case('sa', 'saBC')
vort = sqrt( ((gradw_y(i,j,k)- gradv_z(i,j,k))**2 &
+ (gradu_z(i,j,k)- gradw_x(i,j,k))**2 &
+ (gradv_x(i,j,k)- gradu_y(i,j,k))**2 &
)&
)
kd2 = (kappa_sa*dist(i,j,k))**2
xi = U1(6)*qp(i,j,k,1)/mu(i,j,k)
fv1 = xi**3/(xi**3 + cv1**3)
fv2 = 1.0 - xi/(1 + xi*fv1)
scap = vort + U1(6)*fv2/(kd2)
rsa = min(U1(6)/(Scap*kd2), 10.0)
g = rsa + cw2*(rsa**6 - rsa)
fw = g*( (1.0+cw3**6)/(g**6+cw3**6) )**(1.0/6.0)
R(6) = R(6)/(1.+((-1.0*u1(1)*cb1*scap)+(2.0*u1(1)*cw1*fw*u1(6)/(dist(i,j,k)**2)))*delta_t(i,j,k))
case DEFAULT
Fatal_error
end select
!check if user want to store residue
if(present(Rn)) then
Rn(i,j,k,1:n_var) = Rn(i,j,k,1:n_var) + SF*R(1:n_var)
if(TU) R(:) = Rn(i,j,k,:)
end if
!update
u2(1:n_var) = u1(1:n_var) - R(1:n_var)*(TF*delta_t(i,j,k)/cells(i,j,k)%volume)
! getting primitve variable back variable
u2(1) = u2(1)
u2(2:) = u2(2:)/u2(1)
select case(scheme%turbulence)
case('sst', 'sst2003', 'kkl')
KE = 0.0!u2(6)
case('sa', 'saBC')
!u2(6) = u2(6)*u2(1)
KE=0.0
case DEFAULT
KE = 0.
end select
u2(5) = (flow%gm-1.)*u2(1)*(u2(5) - (0.5*sum(u2(2:4)**2)) - KE)
!check solution for non pyhysical results
if((u2(1) < 0.) .or. (u2(5)) < 0. .or. any(isnan(u2)))then
print*, u2(:)
print*, "R: ", R
print*, "old ", U1
Fatal_error
else !update
qp(i,j,k,1:5) = u2(1:5)
select case(trim(scheme%turbulence))
case('sst', 'sst2003', 'kkl')
if(u2(6)>=0.) then
qp(i,j,k,6) = u2(6)
else
! qp(i,j,k,6) = tk_inf
! qp(i,j,k,6) = (max(qp(i-1,j,k,6),0.) + max(qp(i+1,j,k,6),0.) &
! +max(qp(i,j-1,k,6),0.) + max(qp(i,j+1,k,6),0.) &
! )/4
! qp(i,j,k,6) = 1.e-3*maxval(qp(i-1:i+1,j-1:j+1,k-1:k+1,6))
end if
if(u2(7)>=0.) then
qp(i,j,k,7) = u2(7)
else
! qp(i,j,k,7) = tkl_inf
! qp(i,j,k,7) = (max(qp(i-1,j,k,7),0.) + max(qp(i+1,j,k,7),0.) &
! +max(qp(i,j-1,k,7),0.) + max(qp(i,j+1,k,7),0.) &
! )/4
end if
case('sa', 'saBC')
qp(i,j,k,6) = max(u2(6), 1.e-12)
case DEFAULT
! do nothing
continue
end select
end if
!print*, i,j, R(1:n_var)
end do
end do
end do
case DEFAULT
Fatal_error
end select
end subroutine update_with