-
Andreas Klöckner authoredAndreas Klöckner authored
WENO.F90 28.26 KiB
! 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, &
states, fluxes, metrics, metric_jacobians, flux_derivatives)
implicit none
integer, intent(in) :: nvars, ndim, nx, ny, nz
real, intent(in) :: states(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
real, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)
real generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real generalized_states_frozen(nvars, -2:3)
real generalized_fluxes_frozen(nvars, ndim, -2:3)
real R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
real lambda_pointwise(nvars, -2:3)
real lambda(nvars)
integer i, j, k, l, m
real delta_xi, delta_eta, delta_zeta
real characteristic_fluxes_pos(nvars, -2:3)
real characteristic_fluxes_neg(nvars, -2:3)
real metrics_frozen(ndim, ndim)
real combined_frozen_metrics(ndim)
integer v
! grid spacing in generalized coordinates
delta_xi = 1.0
delta_eta = 1.0
delta_zeta = 1.0
!$loopy begin tagged: to_generalized
call convert_to_generalized(nvars, ndim, nx, ny, nz, &
fluxes, metrics, metric_jacobians, generalized_fluxes)
!$loopy end tagged: to_generalized
! for the next three loops, note carefully that the inner loop indices have a different range
!$loopy begin tagged: flux_x_compute
do k=1,nz
do j=1,ny
do i=0,nx
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), &
metrics_frozen, & combined_frozen_metrics, &
generalized_states_frozen, &
generalized_fluxes_frozen)
call pointwise_eigenvalues(nvars, 1, states(:, 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)
call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda)
call split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen(:, 1, :), &
R_inv, &
lambda, &
characteristic_fluxes_pos, &
characteristic_fluxes_neg)
call weno_flux(nvars, &
generalized_fluxes(:, 1, 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
do v=1,nvars
flux_derivatives_generalized(v, 1, i, j, k) = &
(weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k)) / delta_xi
end do
end do
end do
end do
!$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, 2, states(:, 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)
call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda)
call split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen(:, 2, :), &
R_inv, &
lambda, &
characteristic_fluxes_pos, &
characteristic_fluxes_neg)
call weno_flux(nvars, &
generalized_fluxes(:, 2, 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
do v=1,nvars
flux_derivatives_generalized(v, 2, i, j, k) = &
flux_derivatives_generalized(v, 2, i, j, k) &
+ (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k)) / delta_eta
end do
end do
end do
end do
!$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, 3, states(:, 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)
call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda)
call split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen(:, 3, :), &
R_inv, &
lambda, &
characteristic_fluxes_pos, &
characteristic_fluxes_neg)
call weno_flux(nvars, &
generalized_fluxes(:, 3, 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
do v=1,nvars
flux_derivatives_generalized(v, 3, i, j, k) = &
flux_derivatives_generalized(v, 3, i, j, k) + &
(weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) / delta_zeta
end do
end do
end do end do
!$loopy end tagged: flux_z_diff
!$loopy begin tagged: from_generalized
call convert_from_generalized(nvars, ndim, nx, ny, nz, &
flux_derivatives_generalized, &
metrics, &
metric_jacobians, &
flux_derivatives)
!$loopy end tagged: from_generalized
end subroutine
subroutine pointwise_eigenvalues(nvars, d, states, lambda_pointwise)
implicit none
integer, intent(in) :: nvars
integer, intent(in) :: d
real, intent(in) :: states(nvars, -2:3)
real, intent(out) :: lambda_pointwise(nvars, -2:3)
real u(3), c, p, rho, ke
integer v, k, i
do k=-2,3
rho = states(1,k)
ke = 0.0
do i=1,3
u(i) = states(i+1,k)/rho
ke = ke + u(i)**2
end do
p = (1.4 - 1)*(states(nvars,k) - 0.5*rho*ke)
c = sqrt(1.4*p/rho)
do v=1,nvars-2
lambda_pointwise(v,k) = u(d)
end do
lambda_pointwise(nvars-1,k) = u(d) + c
lambda_pointwise(nvars,k) = u(d) - c
end do
end subroutine
subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lambda_roe)
implicit none
integer, intent(in) :: nvars
integer, intent(in) :: ndim
integer, intent(in) :: d
real, intent(in) :: states(nvars, 2)
real, intent(in) :: metrics_frozen(ndim, ndim)
real, intent(out) :: R(nvars, nvars)
real, intent(out) :: R_inv(nvars, nvars)
real, intent(out) :: lambda_roe(nvars)
integer ik, il, im
real metric_norm(ndim)
real p, ke
real u_orig(ndim, 2), H_orig(2)
real sqrt_rho(2), sqrt_rho_sum
real u(ndim), rho, c, H, q
real u_tilde(ndim), k(ndim), l(ndim), m(ndim)
integer i, j
real 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 u_orig(1,j) = states(ik+1,j)/states(1,j)
u_orig(2,j) = states(il+1,j)/states(1,j)
u_orig(3,j) = states(im+1,j)/states(1,j)
end do
do j=1,2
ke = 0.0
do i=1,ndim
ke = ke + u_orig(i,j)**2
end do
p = (1.4 - 1.0)*(states(nvars,j) - 0.5*states(1,j)*ke)
H_orig(j) = (states(nvars,j) + p)/states(1,j)
end do
do i=1,ndim
metric_norm(i) = 0.0
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
q = 0.0
do i=1,ndim
q = q + u(i)**2
end do
c = sqrt((1.4 - 1.0)*(H - 0.5*q))
b1 = (1.4 - 1.0)/(c**2)
b2 = 1.0 + b1*q**2 - b1*H
u_tilde(1) = 0.0
do i=1,ndim
u_tilde(1) = u_tilde(1) + k(i)*u(i)
end do
u_tilde(2) = 0.0
do i=1,ndim
u_tilde(2) = u_tilde(2) + l(i)*u(i)
end do
u_tilde(3) = 0.0
do i=1,ndim
u_tilde(3) = u_tilde(3) + m(i)*u(i)
end do
alpha = rho/(2.0*c)
beta = 1.0/(2.0*alpha)
lambda_roe(1) = u(1)
lambda_roe(2) = u(1)
lambda_roe(3) = u(1) lambda_roe(4) = u(1) + c
lambda_roe(5) = u(1) - c
R(1,1) = 1.0
R(2,1) = u(1)
R(3,1) = u(2)
R(4,1) = u(3)
R(5,1) = H - 1.0/b1
R(1,2) = 0.0
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)
R(1,3) = 0.0
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,1) = 1.0 - b2
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(2,5) = 0.0
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(3,5) = 0.0
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
end subroutine
subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: lambda_pointwise(nvars, -2:3)
real, intent(in) :: lambda_roe(nvars) real, intent(out) :: lambda(nvars)
real kappa
integer v
integer k
kappa = 1.1
do v=1,nvars
lambda(v) = abs(lambda_roe(v))
do k=-2,3
lambda(v) = max(lambda(v), abs(lambda_pointwise(v,k)))
end do
lambda(v) = kappa*lambda(v)
end do
end subroutine
subroutine convert_to_generalized_frozen( &
nvars, ndim, &
states, &
fluxes, &
metrics, &
metric_jacobians, &
metrics_frozen, &
combined_frozen_metrics, &
generalized_states_frozen, &
generalized_fluxes_frozen)
integer, intent(in) :: nvars, ndim
real, intent(in) :: states(nvars, -2:3)
real, intent(in) :: fluxes(nvars, ndim, -2:3)
real, intent(in) :: metrics(ndim, ndim, 2)
real, intent(in) :: metric_jacobians(2)
real, intent(out) :: metrics_frozen(ndim, ndim)
real, intent(out) :: combined_frozen_metrics(ndim)
real, intent(out) :: generalized_states_frozen(nvars, -2:3)
real, intent(out) :: generalized_fluxes_frozen(nvars, ndim, -2:3)
real jacobian_frozen
integer c, k, v
integer i, j
do j=1,ndim
do i=1,ndim
metrics_frozen(i,j) = 0.0
do c=1,2
metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5*metrics(i,j,c)/metric_jacobians(c)
end do
end do
end do
jacobian_frozen = 0.5*(metric_jacobians(1) + metric_jacobians(2))
do i=1,ndim
combined_frozen_metrics(i) = 0.0
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
do k=-2,3
do v=1,nvars
call mult_mat_vec(ndim, ndim, 1.0, metrics_frozen, &
fluxes(v,:,k), generalized_fluxes_frozen(v,:,k))
end do end do
end subroutine
subroutine convert_to_generalized(nvars, ndim, nx, ny, nz, &
fluxes, &
metrics, &
metric_jacobians, &
generalized_fluxes)
implicit none
integer, intent(in) :: nvars, ndim, nx, ny, nz
real, intent(in) :: fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
real, intent(out) :: generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
integer i, j, k, v
do k=-2,nz+3
do j=-2,ny+3
do i=-2,nx+3
do v=1,nvars
call mult_mat_vec(ndim, ndim, 1.0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
fluxes(v,:,i,j,k), generalized_fluxes(v,:,i,j,k))
end do
end do
end do
end do
end subroutine
subroutine convert_from_generalized(nvars, ndim, nx, ny, nz, &
flux_derivatives_generalized, &
metrics, &
metric_jacobians, &
flux_derivatives)
implicit none
integer, intent(in) :: nvars, ndim, nx, ny, nz
real, intent(in) :: flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
real, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
real, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
integer i, j, k, v
integer d
do k=-2,nz+3
do j=-2,ny+3
do i=-2,nx+3
do v=1,nvars
call solve3x3(metrics(:,:,i,j,k), flux_derivatives_generalized(v,:,i,j,k), &
flux_derivatives(v,:,i,j,k))
do d=1,ndim
flux_derivatives(v,d,i,j,k) = flux_derivatives(v,d,i,j,k)*metric_jacobians(i,j,k)
end do
end do
end do
end do
end do
end subroutine
subroutine solve3x3(a, b, x)
implicit none
real, intent(in) :: a(3,3), b(3)
real, intent(out) :: x(3)
real temp(3,3)
real det_a(1), det_temp(1)
integer k integer i,j
call determinant3x3(a, det_a)
do k=1,3
do j=1,3
do i=1,3
temp(i,j) = a(i,j)
temp(i,k) = b(k)
end do
end do
call determinant3x3(temp, det_temp)
x(k) = det_temp(1)/det_a(1)
end do
end subroutine
subroutine determinant3x3(a, det)
implicit none
real, intent(in) :: a(3,3)
real, intent(out) :: det(1)
det(1) = a(1,1)*(a(2,2)*a(3,3) - a(2,3)*a(3,2)) &
+ a(1,2)*(a(2,3)*a(3,1) - a(2,1)*a(3,3)) &
+ a(1,3)*(a(2,1)*a(3,2) - a(2,2)*a(3,1))
end subroutine
subroutine split_characteristic_fluxes(nvars, &
generalized_states_frozen, &
generalized_fluxes_frozen, &
R_inv, &
lambda, &
characteristic_fluxes_pos, &
characteristic_fluxes_neg)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: generalized_states_frozen(nvars, -2:3)
real, intent(in) :: generalized_fluxes_frozen(nvars, -2:3)
real, intent(in) :: R_inv(nvars, nvars)
real, intent(in) :: lambda(nvars)
real, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
real, intent(out) :: characteristic_fluxes_neg(nvars, -2:3)
integer k, m
integer l
do k=-2,3
do m=1,nvars
characteristic_fluxes_pos(m,k) = 0.0
characteristic_fluxes_neg(m,k) = 0.0
do l=1,nvars
characteristic_fluxes_pos(m,k) = characteristic_fluxes_pos(m,k) &
+ 0.5*R_inv(m,l) &
*(generalized_fluxes_frozen(l,k) + lambda(m)*generalized_states_frozen(l,k))
characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) &
+ 0.5*R_inv(m,l) &
*(generalized_fluxes_frozen(l,k) - lambda(m)*generalized_states_frozen(l,k))
end do
end do
end do
end subroutine
subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, &
characteristic_fluxes_neg, combined_frozen_metrics, R, flux)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: generalized_fluxes(nvars, -2:3)
real, intent(in) :: characteristic_fluxes_pos(nvars, -2:3)
real, intent(in) :: characteristic_fluxes_neg(nvars, -2:3) real, intent(in) :: combined_frozen_metrics
real, intent(in) :: R(nvars, nvars)
real, intent(out) :: flux(nvars)
real consistent(nvars)
real dissipation_pos(nvars)
real dissipation_neg(nvars)
integer v
call consistent_part(nvars, generalized_fluxes, consistent)
call dissipation_part_pos(nvars, characteristic_fluxes_pos, &
combined_frozen_metrics, R, dissipation_pos)
call dissipation_part_neg(nvars, characteristic_fluxes_neg, &
combined_frozen_metrics, R, dissipation_neg)
do v=1,nvars
flux(v) = consistent(v) + dissipation_pos(v) + dissipation_neg(v)
end do
end subroutine
subroutine consistent_part(nvars, generalized_fluxes, consistent)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: generalized_fluxes(nvars, -2:3)
real, intent(out) :: consistent(nvars)
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)
end do
end subroutine
subroutine dissipation_part_pos(nvars, characteristic_fluxes, &
combined_frozen_metrics, R, dissipation_pos)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(in) :: R(nvars, nvars)
real, intent(out) :: dissipation_pos(nvars)
real combined_fluxes(nvars)
call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
call mult_mat_vec(nvars, nvars, -1.0/60, R, combined_fluxes, dissipation_pos)
end subroutine
subroutine weno_combination_pos(nvars, characteristic_fluxes, &
combined_frozen_metrics, combined_fluxes)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: combined_fluxes(nvars)
real w(nvars, 3)
real flux_differences(nvars, 3)
integer v
call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, w)
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
end subroutine
subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, w)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: w(nvars, 3)
real :: C(3) = (/0.1, 0.6, 0.3/)
real :: weights1(0:2) = (/1, -4, 3/)
real :: weights2(0:2) = (/-1, 0, 1/)
real :: weights3(0:2) = (/-3, 4, 1/)
real :: weightsc(0:2) = (/1, -2, 1/)
real IS(3), alpha(3), eps
real w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
integer p, i, j
real sum_alpha(1)
eps = 1e-6*combined_frozen_metrics
p = 2
do i=1,nvars
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)
IS(1) = (1.0/4)*w1sum(1)**2 + (13.0/12)*c1sum(1)**2
IS(2) = (1.0/4)*w2sum(1)**2 + (13.0/12)*c2sum(1)**2
IS(3) = (1.0/4)*w3sum(1)**2 + (13.0/12)*c3sum(1)**2
do j=1,3
alpha(j) = C(j)/(IS(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 flux_differences_pos(nvars, characteristic_fluxes, flux_differences)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(out) :: flux_differences(nvars, 3)
integer i, v
do i=1,3
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)
end do
end do
end subroutine
subroutine dissipation_part_neg(nvars, characteristic_fluxes, & combined_frozen_metrics, R, dissipation_neg)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(in) :: R(nvars, nvars)
real, intent(out) :: dissipation_neg(nvars)
real combined_fluxes(nvars)
call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
call mult_mat_vec(nvars, nvars, 1.0/60, R, combined_fluxes, dissipation_neg)
end subroutine
subroutine weno_combination_neg(nvars, characteristic_fluxes, &
combined_frozen_metrics, combined_fluxes)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: combined_fluxes(nvars)
real w(nvars, 3)
real flux_differences(nvars, 3)
integer v
call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, w)
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
end subroutine
subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, w)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: w(nvars, 3)
real :: C(3) = (/0.1, 0.6, 0.3/)
real :: weights1(0:2) = (/1, -4, 3/)
real :: weights2(0:2) = (/-1, 0, 1/)
real :: weights3(0:2) = (/-3, 4, 1/)
real :: weightsc(0:2) = (/1, -2, 1/)
real IS(3), alpha(3), eps
real w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
integer p, i, j
real sum_alpha(1)
eps = 1e-6*combined_frozen_metrics
p = 2
do i=1,nvars
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)
IS(1) = (1.0/4)*w1sum(1)**2 + (13.0/12)*c1sum(1)**2
IS(2) = (1.0/4)*w2sum(1)**2 + (13.0/12)*c2sum(1)**2
IS(3) = (1.0/4)*w3sum(1)**2 + (13.0/12)*c3sum(1)**2
do j=1,3
alpha(j) = C(j)/(IS(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 flux_differences_neg(nvars, characteristic_fluxes, flux_differences)
implicit none
integer, intent(in) :: nvars
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(out) :: flux_differences(nvars, 3)
integer i, v
do i=1,3
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)
end do
end do
end subroutine
subroutine mult_mat_vec(m, n, alpha, a, b, c)
implicit none
integer, intent(in) :: m, n
real, intent(in) :: alpha
real, intent(in) :: a(m, n)
real, intent(in) :: b(n)
real, intent(out) :: c(m)
integer i, j
do i=1,m
c(i) = 0.0
do j=1,n
c(i) = c(i) + alpha*a(i,j)*b(j)
end do
end do
end subroutine
subroutine sum_array(n, a, a_sum)
implicit none
integer, intent(in) :: n
real, intent(in) :: a(n)
real, intent(out) :: a_sum(1)
integer i
a_sum(1) = 0.0
do i=1,n
a_sum(1) = a_sum(1) + a(i)
end do
end subroutine
subroutine weighted_sum(n, w, a, a_sum)
implicit none
integer, intent(in) :: n
real, intent(in) :: a(n), w(n) real, intent(out) :: a_sum(1)
integer i
a_sum(1) = 0.0
do i=1,n
a_sum(1) = a_sum(1) + w(i)*a(i)
end do
end subroutine
!$loopy begin
!
! prg = lp.parse_fortran(lp.c_preprocess(SOURCE), FILENAME)
! prg = lp.fix_parameters(prg, ndim=3, nvars=5, _remove=False)
! RESULT = prg
!
!$loopy end