Newer
Older
! Copyright (C) 2019 Timothy A. Smith
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, &
Timothy A. Smith
committed
states, fluxes, metrics, metric_jacobians, flux_derivatives)
integer, intent(in) :: nvars, ndim, nx, ny, nz, d
real*8, intent(in) :: states(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
real*8, intent(in) :: fluxes(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
real*8, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real*8, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
real*8, intent(out) :: flux_derivatives(nvars, nx, ny, nz)
real*8 weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)
real*8 metrics_frozen(ndim, ndim)
real*8 combined_frozen_metrics(ndim)
real*8 generalized_states_frozen(nvars, -2:3)
real*8 generalized_fluxes_frozen(nvars, -2:3)
real*8 R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
real*8 lambda_pointwise(nvars, -2:3)
real*8 wavespeeds(nvars)
real*8 characteristic_fluxes_pos(nvars, -2:3)
real*8 characteristic_fluxes_neg(nvars, -2:3)
integer i, j, k, v
! for the next three loops, note carefully that the inner loop indices have a different range
!$loopy begin tagged: flux_x_compute
Timothy A. Smith
committed
call convert_to_generalized_frozen(nvars, ndim, &
states(:, i-2:i+3, j, k), &
fluxes(:, i-2:i+3, j, k), &
metrics(:, :, i:i+1, j, k), &
metric_jacobians(i:i+1, j, k), &
Timothy A. Smith
committed
combined_frozen_metrics, &
generalized_states_frozen, &
call pointwise_eigenvalues(nvars, ndim, 1, states(:, i-2:i+3, j, k), &
metrics(:, :, i-2:i+3, j, k), lambda_pointwise)
call roe_eigensystem(nvars, ndim, 1, states(:, i:i+1, j, k), &
metrics_frozen, R, R_inv, lambda_roe)
Timothy A. Smith
committed
call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
Timothy A. Smith
committed
call split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen, &
Timothy A. Smith
committed
wavespeeds, &
characteristic_fluxes_pos, &
Timothy A. Smith
committed
Timothy A. Smith
committed
call weno_flux(nvars, &
fluxes(:, i-2:i+3, j, k), &
characteristic_fluxes_pos, characteristic_fluxes_neg, &
combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k))
end do
end do
end do
!$loopy end tagged: flux_x_compute
!$loopy begin tagged: flux_x_diff
do k=1,nz
do j=1,ny
do i=1,nx
flux_derivatives(v, i, j, k) = &
(weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k))
!$loopy end tagged: flux_x_diff
!$loopy begin tagged: flux_y_compute
do k=1,nz
do i=1,nx
do j=0,ny
call convert_to_generalized_frozen(nvars, ndim, &
states(:, i, j-2:j+3, k), &
fluxes(:, i, j-2:j+3, k), &
metrics(:, :, i, j:j+1, k), &
metric_jacobians(i, j:j+1, k), &
metrics_frozen, &
combined_frozen_metrics, &
generalized_states_frozen, &
generalized_fluxes_frozen)
call pointwise_eigenvalues(nvars, ndim, 2, states(:, i, j-2:j+3, k), &
metrics(:, :, i, j-2:j+3, k), lambda_pointwise)
call roe_eigensystem(nvars, ndim, 2, states(:, i, j:j+1, k), &
metrics_frozen, R, R_inv, lambda_roe)
Timothy A. Smith
committed
call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
call split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen, &
Timothy A. Smith
committed
wavespeeds, &
characteristic_fluxes_pos, &
characteristic_fluxes_neg)
call weno_flux(nvars, &
fluxes(:, i, j-2:j+3, k), &
characteristic_fluxes_pos, characteristic_fluxes_neg, &
combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k))
end do
end do
end do
!$loopy end tagged: flux_y_compute
!$loopy begin tagged: flux_y_diff
do k=1,nz
do j=1,ny
do i=1,nx
flux_derivatives(v, i, j, k) = &
(weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k))
!$loopy end tagged: flux_y_diff
!$loopy begin tagged: flux_z_compute
do j=1,ny
do i=1,nx
do k=0,nz
call convert_to_generalized_frozen(nvars, ndim, &
states(:, i, j, k-2:k+3), &
fluxes(:, i, j, k-2:k+3), &
metrics(:, :, i, j, k:k+1), &
metric_jacobians(i, j, k:k+1), &
metrics_frozen, &
combined_frozen_metrics, &
generalized_states_frozen, &
generalized_fluxes_frozen)
call pointwise_eigenvalues(nvars, ndim, 3, states(:, i, j, k-2:k+3), &
metrics(:, :, i, j, k-2:k+3), lambda_pointwise)
call roe_eigensystem(nvars, ndim, 3, states(:, i, j, k:k+1), &
metrics_frozen, R, R_inv, lambda_roe)
Timothy A. Smith
committed
call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
call split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen, &
Timothy A. Smith
committed
wavespeeds, &
characteristic_fluxes_pos, &
characteristic_fluxes_neg)
call weno_flux(nvars, &
fluxes(:, i, j, k-2:k+3), &
characteristic_fluxes_pos, characteristic_fluxes_neg, &
combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k))
end do
end do
end do
!$loopy end tagged: flux_z_compute
!$loopy begin tagged: flux_z_diff
do k=1,nz
do j=1,ny
do i=1,nx
flux_derivatives(v, i, j, k) = &
(weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1))
flux_derivatives(v, i, j, k) = 0.0d0
!$loopy end tagged: flux_z_diff
subroutine pointwise_eigenvalues(nvars, ndim, d, states, metrics, lambda_pointwise)
Timothy A. Smith
committed
implicit none
integer, intent(in) :: nvars
Timothy A. Smith
committed
integer, intent(in) :: d
real*8, intent(in) :: states(nvars, -2:3)
real*8, intent(in) :: metrics(ndim, ndim, -2:3)
real*8, intent(out) :: lambda_pointwise(nvars, -2:3)
real*8 u(3), c, p, rho, ke, metric_norm
Timothy A. Smith
committed
do k=-2,3
rho = states(1,k)
Timothy A. Smith
committed
do i=1,3
u(i) = states(i+1,k)/rho
ke = ke + u(i)**2
metric_norm = metric_norm + metrics(d,i,k)**2
Timothy A. Smith
committed
end do
p = (1.4d0 - 1)*(states(nvars,k) - 0.5d0*rho*ke)
c = sqrt(1.4d0*p/rho)
Timothy A. Smith
committed
Timothy A. Smith
committed
do v=1,nvars-2
Timothy A. Smith
committed
lambda_pointwise(v,k) = u(d)
end do
lambda_pointwise(nvars-1,k) = u(d) + c*metric_norm
lambda_pointwise(nvars,k) = u(d) - c*metric_norm
Timothy A. Smith
committed
end do
subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lambda_roe)
Timothy A. Smith
committed
implicit none
integer, intent(in) :: nvars
Timothy A. Smith
committed
integer, intent(in) :: d
real*8, intent(in) :: states(nvars, 2)
real*8, intent(in) :: metrics_frozen(ndim, ndim)
real*8, intent(out) :: R(nvars, nvars)
real*8, intent(out) :: R_inv(nvars, nvars)
real*8, intent(out) :: lambda_roe(nvars)
real*8 metric_norm(ndim)
real*8 p, ke
real*8 u_orig(ndim, 2), H_orig(2)
real*8 sqrt_rho(2), sqrt_rho_sum
real*8 u(ndim), rho, c, H, q
real*8 u_tilde(ndim), k(ndim), l(ndim), m(ndim)
real*8 b1, b2, alpha, beta
ik = d
il = d+1
im = d+2
if (il > 3) il = il - 3
if (im > 3) im = im - 3
do j=1,2
do i=1,ndim
u_orig(i,j) = states(i+1,j)/states(1,j)
end do
do i=1,ndim
ke = ke + u_orig(i,j)**2
end do
p = (1.4d0 - 1)*(states(nvars,j) - 0.5d0*states(1,j)*ke)
H_orig(j) = (states(nvars,j) + p)/states(1,j)
end do
do i=1,ndim
do j=1,ndim
metric_norm(i) = metric_norm(i) + metrics_frozen(i,j)**2
end do
metric_norm(i) = sqrt(metric_norm(i))
end do
do i=1,ndim
k(i) = metrics_frozen(ik,i)/metric_norm(ik)
l(i) = metrics_frozen(il,i)/metric_norm(il)
m(i) = metrics_frozen(im,i)/metric_norm(im)
end do
sqrt_rho(1) = sqrt(states(1,1))
sqrt_rho(2) = sqrt(states(1,2))
sqrt_rho_sum = sqrt_rho(1) + sqrt_rho(2)
do i=1,ndim
u(i) = (u_orig(i,1)*sqrt_rho(1) + u_orig(i,2)*sqrt_rho(2))/sqrt_rho_sum
end do
rho = sqrt_rho(1)*sqrt_rho(2)
H = (H_orig(1)*sqrt_rho(1) + H_orig(2)*sqrt_rho(2))/sqrt_rho_sum
do i=1,ndim
q = q + u(i)**2
end do
c = sqrt((1.4d0 - 1.0d0)*(H - 0.5d0*q))
b1 = (1.4d0 - 1.0d0)/(c**2)
b2 = 1.0d0 + b1*q - b1*H
do i=1,ndim
u_tilde(1) = u_tilde(1) + k(i)*u(i)
end do
do i=1,ndim
u_tilde(2) = u_tilde(2) + l(i)*u(i)
end do
do i=1,ndim
u_tilde(3) = u_tilde(3) + m(i)*u(i)
end do
alpha = rho/(2.0d0*c)
beta = 1.0d0/(2.0d0*alpha)
lambda_roe(1) = u_tilde(1)*metric_norm(ik)
lambda_roe(2) = u_tilde(1)*metric_norm(ik)
lambda_roe(3) = u_tilde(1)*metric_norm(ik)
lambda_roe(4) = u_tilde(1)*metric_norm(ik) + c*metric_norm(ik)
lambda_roe(5) = u_tilde(1)*metric_norm(ik) - c*metric_norm(ik)
Timothy A. Smith
committed
Timothy A. Smith
committed
R(2,1) = u(1)
R(3,1) = u(2)
R(4,1) = u(3)
Timothy A. Smith
committed
Timothy A. Smith
committed
R(2,2) = rho*l(1)
R(3,2) = rho*l(2)
R(4,2) = rho*l(3)
R(5,2) = rho*u_tilde(2)
Timothy A. Smith
committed
R(2,3) = rho*m(1)
R(3,3) = rho*m(2)
R(4,3) = rho*m(3)
R(5,3) = rho*u_tilde(3)
R(1,4) = alpha
R(2,4) = alpha*(u(1) + c*k(1))
R(3,4) = alpha*(u(2) + c*k(2))
R(4,4) = alpha*(u(3) + c*k(3))
R(5,4) = alpha*(H + c*u_tilde(1))
R(1,5) = alpha
R(2,5) = alpha*(u(1) - c*k(1))
R(3,5) = alpha*(u(2) - c*k(2))
R(4,5) = alpha*(u(3) - c*k(3))
R(5,5) = alpha*(H - c*u_tilde(1))
R_inv(1,2) = b1*u(1)
R_inv(1,3) = b1*u(2)
R_inv(1,4) = b1*u(3)
R_inv(1,5) = -b1
R_inv(2,1) = -u_tilde(2)/rho
R_inv(2,2) = l(1)/rho
R_inv(2,3) = l(2)/rho
R_inv(2,4) = l(3)/rho
R_inv(3,1) = -u_tilde(3)/rho
R_inv(3,2) = m(1)/rho
R_inv(3,3) = m(2)/rho
R_inv(3,4) = m(3)/rho
R_inv(4,1) = beta*(b2 - u_tilde(1)/c)
R_inv(4,2) = -beta*(b1*u(1) - k(1)/c)
R_inv(4,3) = -beta*(b1*u(2) - k(2)/c)
R_inv(4,4) = -beta*(b1*u(3) - k(3)/c)
R_inv(4,5) = beta*b1
R_inv(5,1) = beta*(b2 + u_tilde(1)/c)
R_inv(5,2) = -beta*(b1*u(1) + k(1)/c)
R_inv(5,3) = -beta*(b1*u(2) + k(2)/c)
R_inv(5,4) = -beta*(b1*u(3) + k(3)/c)
R_inv(5,5) = beta*b1
Timothy A. Smith
committed
subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
Timothy A. Smith
committed
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: lambda_pointwise(nvars, -2:3)
real*8, intent(in) :: lambda_roe(nvars)
real*8, intent(out) :: wavespeeds(nvars)
Timothy A. Smith
committed
integer v
Timothy A. Smith
committed
kappa = 1.1d0
Timothy A. Smith
committed
wavespeeds(v) = abs(lambda_roe(v))
Timothy A. Smith
committed
wavespeeds(v) = max(wavespeeds(v), abs(lambda_pointwise(v,k)))
Timothy A. Smith
committed
wavespeeds(v) = kappa*wavespeeds(v)
subroutine convert_to_generalized_frozen( &
Timothy A. Smith
committed
nvars, ndim, &
states, &
fluxes, &
metrics, &
metric_jacobians, &
Timothy A. Smith
committed
combined_frozen_metrics, &
generalized_states_frozen, &
Timothy A. Smith
committed
integer, intent(in) :: nvars, ndim
real*8, intent(in) :: states(nvars, -2:3)
real*8, intent(in) :: fluxes(nvars, -2:3)
real*8, intent(in) :: metrics(ndim, ndim, 2)
real*8, intent(in) :: metric_jacobians(2)
real*8, intent(out) :: metrics_frozen(ndim, ndim)
real*8, intent(out) :: combined_frozen_metrics(ndim)
real*8, intent(out) :: generalized_states_frozen(nvars, -2:3)
real*8, intent(out) :: generalized_fluxes_frozen(nvars, -2:3)
Timothy A. Smith
committed
Timothy A. Smith
committed
do j=1,ndim
do i=1,ndim
Timothy A. Smith
committed
do c=1,2
metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5d0*metrics(i,j,c)/metric_jacobians(c)
jacobian_frozen = 0.5d0*(metric_jacobians(1) + metric_jacobians(2))
Timothy A. Smith
committed
do i=1,ndim
combined_frozen_metrics(i) = 0.0d0
Timothy A. Smith
committed
do j=1,ndim
combined_frozen_metrics(i) = combined_frozen_metrics(i) + metrics_frozen(i,j)**2
end do
end do
do k=-2,3
do v=1,nvars
generalized_states_frozen(v,k) = states(v,k)/jacobian_frozen
end do
end do
! FIXME: this is a bug
! ... we actually need to recompute the fluxes using frozen metrics
! ... but this works when the metric matrix is identity
generalized_fluxes_frozen(v,k) = fluxes(v,k)
end do
end do
end subroutine
Timothy A. Smith
committed
subroutine split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen, &
R_inv, &
Timothy A. Smith
committed
wavespeeds, &
characteristic_fluxes_pos, &
Timothy A. Smith
committed
integer, intent(in) :: nvars
real*8, intent(in) :: generalized_states_frozen(nvars, -2:3)
real*8, intent(in) :: generalized_fluxes_frozen(nvars, -2:3)
real*8, intent(in) :: R_inv(nvars, nvars)
real*8, intent(in) :: wavespeeds(nvars)
real*8, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
real*8, intent(out) :: characteristic_fluxes_neg(nvars, -2:3)
Timothy A. Smith
committed
integer k, m
Timothy A. Smith
committed
characteristic_fluxes_pos(m,k) = 0.0d0
characteristic_fluxes_neg(m,k) = 0.0d0
do l=1,nvars
characteristic_fluxes_pos(m,k) = characteristic_fluxes_pos(m,k) &
Timothy A. Smith
committed
*(generalized_fluxes_frozen(l,k) + wavespeeds(m)*generalized_states_frozen(l,k))
characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) &
Timothy A. Smith
committed
*(generalized_fluxes_frozen(l,k) - wavespeeds(m)*generalized_states_frozen(l,k))
Timothy A. Smith
committed
end subroutine
Timothy A. Smith
committed
subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, &
Timothy Smith
committed
characteristic_fluxes_neg, combined_frozen_metric, R, flux)
Timothy A. Smith
committed
integer, intent(in) :: nvars
real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
real*8, intent(in) :: characteristic_fluxes_pos(nvars, -2:3)
real*8, intent(in) :: characteristic_fluxes_neg(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(in) :: R(nvars, nvars)
real*8, intent(out) :: flux(nvars)
real*8 consistent(nvars)
real*8 dissipation_pos(nvars)
real*8 dissipation_neg(nvars)
Timothy A. Smith
committed
call consistent_part(nvars, generalized_fluxes, consistent)
call dissipation_part_pos(nvars, characteristic_fluxes_pos, &
Timothy Smith
committed
combined_frozen_metric, R, dissipation_pos)
Timothy A. Smith
committed
call dissipation_part_neg(nvars, characteristic_fluxes_neg, &
Timothy Smith
committed
combined_frozen_metric, R, dissipation_neg)
Timothy A. Smith
committed
do v=1,nvars
flux(v) = consistent(v) + dissipation_pos(v) + dissipation_neg(v)
end do
Timothy A. Smith
committed
end subroutine
Timothy A. Smith
committed
subroutine consistent_part(nvars, generalized_fluxes, consistent)
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
real*8, intent(out) :: consistent(nvars)
Timothy A. Smith
committed
integer v
do v=1,nvars
consistent(v) = generalized_fluxes(v,-2) - 8*generalized_fluxes(v,-1) &
+ 37*generalized_fluxes(v,0) + 37*generalized_fluxes(v,1) &
- 8*generalized_fluxes(v,2) + generalized_fluxes(v,3)
consistent(v) = consistent(v)/60.0d0
Timothy A. Smith
committed
end do
Timothy A. Smith
committed
end subroutine
Timothy A. Smith
committed
subroutine dissipation_part_pos(nvars, characteristic_fluxes, &
Timothy Smith
committed
combined_frozen_metric, R, dissipation_pos)
Timothy A. Smith
committed
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(in) :: R(nvars, nvars)
real*8, intent(out) :: dissipation_pos(nvars)
Timothy A. Smith
committed
real*8 combined_fluxes(nvars)
Timothy Smith
committed
call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metric, combined_fluxes)
call mult_mat_vec(nvars, nvars, -1.0d0/60, R, combined_fluxes, dissipation_pos)
Timothy A. Smith
committed
subroutine weno_combination_pos(nvars, characteristic_fluxes, &
Timothy Smith
committed
combined_frozen_metric, combined_fluxes)
Timothy A. Smith
committed
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(out) :: combined_fluxes(nvars)
real*8 w(nvars, 3)
real*8 flux_differences(nvars, 3)
Timothy Smith
committed
call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric, w)
Timothy A. Smith
committed
call flux_differences_pos(nvars, characteristic_fluxes, flux_differences)
do v=1,nvars
combined_fluxes(v) = (20*w(v,1) - 1)*flux_differences(v,1) &
- (10*(w(v,1) + w(v,2)) - 5)*flux_differences(v,2) &
+ flux_differences(v,3)
end do
Timothy Smith
committed
subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric, w)
Timothy A. Smith
committed
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(out) :: w(nvars, 3)
real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
real*8 IS(nvars, 3), alpha(3), eps
integer p, i, j
real*8 sum_alpha(1)
Timothy Smith
committed
eps = 1.0d-6*combined_frozen_metric
p = 2
call oscillation_pos(nvars, characteristic_fluxes, IS)
do i=1,nvars
sum_alpha(1) = 0.0d0
do j=1,3
alpha(j) = C(j)/(IS(i, j) + eps)**p
sum_alpha(1) = sum_alpha(1) + alpha(j)
end do
!call sum_array(3, alpha, sum_alpha)
do j=1,3
w(i,j) = alpha(j)/sum_alpha(1)
end do
end do
end subroutine
subroutine oscillation_pos(nvars, characteristic_fluxes, oscillation)
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
real*8, intent(out) :: oscillation(nvars, 3)
real*8 :: weights1(0:2) = (/ 1.0d0, -4.0d0, 3.0d0/)
real*8 :: weights2(0:2) = (/-1.0d0, 0.0d0, 1.0d0/)
real*8 :: weights3(0:2) = (/-3.0d0, 4.0d0, -1.0d0/)
real*8 :: weightsc(0:2) = (/ 1.0d0, -2.0d0, 1.0d0/)
real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
Timothy A. Smith
committed
integer i
call weighted_sum(3, weights1, characteristic_fluxes(i,-2:0), w1sum)
call weighted_sum(3, weightsc, characteristic_fluxes(i,-2:0), c1sum)
call weighted_sum(3, weights2, characteristic_fluxes(i,-1:1), w2sum)
call weighted_sum(3, weightsc, characteristic_fluxes(i,-1:1), c2sum)
call weighted_sum(3, weights3, characteristic_fluxes(i,0:2), w3sum)
call weighted_sum(3, weightsc, characteristic_fluxes(i,0:2), c3sum)
Timothy A. Smith
committed
oscillation(i, 1) = (1.0d0/4)*w1sum(1)**2 + (13.0d0/12)*c1sum(1)**2
oscillation(i, 2) = (1.0d0/4)*w2sum(1)**2 + (13.0d0/12)*c2sum(1)**2
oscillation(i, 3) = (1.0d0/4)*w3sum(1)**2 + (13.0d0/12)*c3sum(1)**2
end do
end subroutine
Timothy A. Smith
committed
subroutine flux_differences_pos(nvars, characteristic_fluxes, flux_differences)
Timothy A. Smith
committed
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
real*8, intent(out) :: flux_differences(nvars, 3)
Timothy A. Smith
committed
integer i, v
Timothy A. Smith
committed
do v=1,nvars
flux_differences(v,i) = -characteristic_fluxes(v,-3+i) + 3*characteristic_fluxes(v,-2+i) &
- 3*characteristic_fluxes(v,-1+i) + characteristic_fluxes(v,i)
Timothy A. Smith
committed
end do
Timothy A. Smith
committed
subroutine dissipation_part_neg(nvars, characteristic_fluxes, &
Timothy Smith
committed
combined_frozen_metric, R, dissipation_neg)
Timothy A. Smith
committed
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(in) :: R(nvars, nvars)
real*8, intent(out) :: dissipation_neg(nvars)
Timothy A. Smith
committed
real*8 combined_fluxes(nvars)
Timothy Smith
committed
call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metric, combined_fluxes)
call mult_mat_vec(nvars, nvars, 1.0d0/60, R, combined_fluxes, dissipation_neg)
Timothy A. Smith
committed
subroutine weno_combination_neg(nvars, characteristic_fluxes, &
Timothy Smith
committed
combined_frozen_metric, combined_fluxes)
Timothy A. Smith
committed
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(out) :: combined_fluxes(nvars)
real*8 w(nvars, 3)
real*8 flux_differences(nvars, 3)
Timothy Smith
committed
call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric, w)
Timothy A. Smith
committed
call flux_differences_neg(nvars, characteristic_fluxes, flux_differences)
do v=1,nvars
combined_fluxes(v) = (20*w(v,1) - 1)*flux_differences(v,1) &
- (10*(w(v,1) + w(v,2)) - 5)*flux_differences(v,2) &
+ flux_differences(v,3)
end do
Timothy Smith
committed
subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric, w)
Timothy A. Smith
committed
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
Timothy Smith
committed
real*8, intent(in) :: combined_frozen_metric
real*8, intent(out) :: w(nvars, 3)
real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
real*8 IS(nvars, 3), alpha(3), eps
integer p, i, j
real*8 sum_alpha(1)
Timothy Smith
committed
eps = 1.0d-6*combined_frozen_metric
p = 2
call oscillation_neg(nvars, characteristic_fluxes, IS)
do i=1,nvars
do j=1,3
alpha(j) = C(j)/(IS(i, j) + eps)**p
end do
call sum_array(3, alpha, sum_alpha)
do j=1,3
w(i,j) = alpha(j)/sum_alpha(1)
end do
end do
end subroutine
subroutine oscillation_neg(nvars, characteristic_fluxes, oscillation)
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
real*8, intent(out) :: oscillation(nvars, 3)
real*8 :: weights1(0:2) = (/ 1.0d0, -4.0d0, 3.0d0/)
real*8 :: weights2(0:2) = (/-1.0d0, 0.0d0, 1.0d0/)
real*8 :: weights3(0:2) = (/-3.0d0, 4.0d0, -1.0d0/)
real*8 :: weightsc(0:2) = (/ 1.0d0, -2.0d0, 1.0d0/)
Timothy A. Smith
committed
real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
integer i
call weighted_sum(3, weights1, characteristic_fluxes(i,3:1:-1), w1sum)
call weighted_sum(3, weightsc, characteristic_fluxes(i,3:1:-1), c1sum)
call weighted_sum(3, weights2, characteristic_fluxes(i,2:0:-1), w2sum)
call weighted_sum(3, weightsc, characteristic_fluxes(i,2:0:-1), c2sum)
call weighted_sum(3, weights3, characteristic_fluxes(i,1:-1:-1), w3sum)
call weighted_sum(3, weightsc, characteristic_fluxes(i,1:-1:-1), c3sum)
oscillation(i, 1) = (1.0d0/4)*w1sum(1)**2 + (13.0d0/12)*c1sum(1)**2
oscillation(i, 2) = (1.0d0/4)*w2sum(1)**2 + (13.0d0/12)*c2sum(1)**2
oscillation(i, 3) = (1.0d0/4)*w3sum(1)**2 + (13.0d0/12)*c3sum(1)**2
Timothy A. Smith
committed
subroutine flux_differences_neg(nvars, characteristic_fluxes, flux_differences)
implicit none
integer, intent(in) :: nvars
real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
real*8, intent(out) :: flux_differences(nvars, 3)
Timothy A. Smith
committed
integer i, v
Timothy A. Smith
committed
do v=1,nvars
flux_differences(v,i) = -characteristic_fluxes(v,1-i) + 3*characteristic_fluxes(v,2-i) &
- 3*characteristic_fluxes(v,3-i) + characteristic_fluxes(v,4-i)
Timothy A. Smith
committed
end do
subroutine mult_mat_vec(m, n, alpha, a, b, c)
real*8, intent(in) :: alpha
real*8, intent(in) :: a(m, n)
real*8, intent(in) :: b(n)
real*8, intent(out) :: c(m)
Andreas Klöckner
committed
accumulator = accumulator + alpha*a(i,j)*b(j)
c(i) = accumulator
end do
end subroutine
subroutine sum_array(n, a, a_sum)
implicit none
integer, intent(in) :: n
real*8, intent(in) :: a(n)
real*8, intent(out) :: a_sum(1)
a_sum(1) = a_sum(1) + a(i)
subroutine weighted_sum(n, w, a, a_sum)
implicit none
integer, intent(in) :: n
real*8, intent(in) :: a(n), w(n)
real*8, intent(out) :: a_sum(1)
integer i
do i=1,n
a_sum(1) = a_sum(1) + w(i)*a(i)
end do
end subroutine
! prg = lp.parse_fortran(lp.c_preprocess(SOURCE), FILENAME)
Kaushik Kulkarni
committed
! prg = lp.fix_parameters(prg, ndim=3, nvars=5, _remove=False)
Timothy A. Smith
committed
!
! knl = prg["compute_flux_derivatives"]
Timothy A. Smith
committed
!
! knl = lp.assume(knl, "nx > 0 and ny > 0 and nz > 0")
Timothy A. Smith
committed
!
! knl = lp.set_temporary_scope(knl, "weno_flux_tmp",
Timothy A. Smith
committed
! lp.AddressSpace.GLOBAL)
!
Timothy A. Smith
committed
!
! RESULT = prg
!
!$loopy end