| real(kind=rk), |
intent(in) |
| | :: |
material_right(nSides,1,4) | integer :: iSide, left, right, iDof
real(kind=rk) :: left_mu, right_mu
real(kind=rk) :: left_epsi, right_epsi
real(kind=rk) :: left_gam, right_gam
real(kind=rk) :: left_chi, right_chi
real(kind=rk) :: left_speedOfLight, right_speedOfLight
real(kind=rk) :: inv_denom_mu, inv_denom_epsi
! --------------------------------------------------------------------------
! flux for D_X, B_Z, D_y and B_z
do iDof = 1, nFaceDofs
do iSide = 1, nSides
! The position of the left and right element in the state
! vector.
left = leftPos(iSide)
right = rightPos(iSide)
! The material parameters of the left and right element
left_mu = material_left(iSide,1,1)
left_epsi = material_left(iSide,1,2)
left_gam = material_left(iSide,1,3)
left_chi = material_left(iSide,1,4)
left_speedOfLight = 1.0_rk / sqrt( left_mu * left_epsi )
right_mu = material_right(iSide,1,1)
right_epsi = material_right(iSide,1,2)
right_gam = material_right(iSide,1,3)
right_chi = material_right(iSide,1,4)
right_speedOfLight = 1.0_rk / sqrt( right_mu * right_epsi )
! The inverse of the denominators
inv_denom_mu = 1.0_rk / ((-1.0_rk)*left_speedOfLight*left_mu - right_speedOfLight*right_mu)
inv_denom_epsi = 1.0_rk / (left_speedOfLight*left_epsi + right_speedOfLight*right_epsi)
! flux for D_x
faceFlux(left,iDof,var(1),2) = inv_denom_mu * ( &
& (-1.0_rk)*left_chi*faceRep(left,iDof,var(1),2)/left_epsi &
& - (-1.0_rk)*right_chi*faceRep(right,iDof,var(1),1)/right_epsi &
& ) &
& + ( &
& left_chi*left_chi*faceRep(left,iDof,var(7),2) &
& + right_chi*right_chi*faceRep(right,iDof,var(7),1) &
& ) / ( left_epsi*left_mu + right_epsi*right_mu )
! flux for B_x
faceFlux(left,iDof,var(4),2) = inv_denom_epsi * ( &
& (-1.0_rk)*left_gam*faceRep(left,iDof,var(4),2)/left_mu &
& - (-1.0_rk)*right_gam*faceRep(right,iDof,var(4),1)/right_mu) &
& + ( &
& left_gam*left_gam*faceRep(left,iDof,var(8),2) &
& + right_gam*right_gam*faceRep(right,iDof,var(8),1) &
& ) / ( left_epsi*left_mu + right_epsi*right_mu )
! flux for phi (electric correction)
faceFlux(left,iDof,var(7),2) = ( &
& left_speedOfLight*faceRep(left,iDof,var(1),2) &
& + right_speedOfLight*faceRep(right,iDof,var(1),1) &
& + left_speedOfLight*left_chi*faceRep(left,iDof,var(7),2) &
& - right_speedOfLight*right_chi*faceRep(right,iDof,var(7),1) &
& ) / (left_speedOfLight + right_speedOfLight)
! flux for psi (magnetic correction)
faceFlux(left,iDof,var(8),2) = ( &
& left_speedOfLight*faceRep(left,iDof,var(4),2) &
& + right_speedOfLight*faceRep(right,iDof,var(4),1) &
& + left_speedOfLight*left_gam*faceRep(left,iDof,var(8),2) &
& - right_speedOfLight*right_gam*faceRep(right,iDof,var(8),1) &
& ) / (left_speedOfLight + right_speedOfLight)
! the flux for D_y
faceFlux(left,iDof,var(2),2) = ( &
& ( ((-1.0_rk)*faceRep(left,iDof,var(2),2) / left_epsi) &
& - ((-1.0_rk)*faceRep(right,iDof,var(2),1) / right_epsi) ) &
& - ( left_speedOfLight * faceRep(left,iDof,var(6),2) &
& + right_speedOfLight * faceRep(right,iDof,var(6),1) ))
! the flux for B_z
faceFlux(left,iDof,var(6),2) = ( &
& ( left_speedOfLight * faceRep(left,iDof,var(2),2) &
& + right_speedOfLight * faceRep(right,iDof,var(2),1) )&
& + ( ( faceRep(left,iDof,var(6),2) / left_mu ) &
& - ( faceRep(right,iDof,var(6),1) / right_mu ) ))
! the flux for D_z
faceFlux(left,iDof,var(3),2) = ( &
& ( ( (-1.0_rk) * faceRep(left,iDof,var(3),2) / left_epsi ) &
& - ( (-1.0_rk) * faceRep(right,iDof,var(3),1) / right_epsi ) ) &
& + ( left_speedOfLight * faceRep(left,iDof,var(5),2) &
& + right_speedOfLight * faceRep(right,iDof,var(5),1) ) &
& )
! the flux for B_y
faceFlux(left,iDof,var(5),2) = ( &
& ( (-1.0_rk) * left_speedOfLight * faceRep(left,iDof,var(3),2) &
& - right_speedOfLight * faceRep(right,iDof,var(3),1) ) &
& + ( ( faceRep(left,iDof,var(5),2) / left_mu ) &
& - ( faceRep(right,iDof,var(5),1) / right_mu) ) &
& )
! Normalize the calculated fluxes
faceFlux(left,iDof,var(2),2) = inv_denom_mu * faceFlux(left,iDof,var(2),2)
faceFlux(left,iDof,var(6),2) = inv_denom_epsi * faceFlux(left,iDof,var(6),2)
faceFlux(left,iDof,var(3),2) = inv_denom_mu * faceFlux(left,iDof,var(3),2)
faceFlux(left,iDof,var(5),2) = inv_denom_epsi * faceFlux(left,iDof,var(5),2)
! Assign the same flux for both adjacent elements
faceFlux(right,iDof,:,1) = faceFlux(left,iDof,:,2)
end do
end do |