diff --git a/WENO.F90 b/WENO.F90
index 18e30b75ba8f4946d322575a8c1040559fb2b166..c07cf33d3edf30a8765265549f417a611d4ac390 100644
--- a/WENO.F90
+++ b/WENO.F90
@@ -23,34 +23,34 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
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 metric_solve_tmp(ndim)
- real R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
- real lambda_pointwise(nvars, -2:3)
- real wavespeeds(nvars)
+ real*8, intent(in) :: states(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
+ real*8, intent(in) :: fluxes(nvars, ndim, -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, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+
+ real*8 flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+ real*8 weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)
+
+ real*8 generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+ real*8 generalized_states_frozen(nvars, -2:3)
+ real*8 generalized_fluxes_frozen(nvars, ndim, -2:3)
+ real*8 metric_solve_tmp(ndim)
+ real*8 R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
+ real*8 lambda_pointwise(nvars, -2:3)
+ real*8 wavespeeds(nvars)
integer i, j, k
- 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)
+ real*8 delta_xi, delta_eta, delta_zeta
+ real*8 characteristic_fluxes_pos(nvars, -2:3)
+ real*8 characteristic_fluxes_neg(nvars, -2:3)
+ real*8 metrics_frozen(ndim, ndim)
+ real*8 combined_frozen_metrics(ndim)
integer v, d
! grid spacing in generalized coordinates
- delta_xi = 1.0
- delta_eta = 1.0
- delta_zeta = 1.0
+ delta_xi = 1.0d0
+ delta_eta = 1.0d0
+ delta_zeta = 1.0d0
!$loopy begin tagged: to_generalized
!call convert_to_generalized(nvars, ndim, nx, ny, nz, &
@@ -61,7 +61,7 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
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), &
+ call mult_mat_vec(ndim, ndim, 1.0d0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
fluxes(v,:,i,j,k), generalized_fluxes(v,:,i,j,k))
end do
end do
@@ -250,21 +250,21 @@ subroutine pointwise_eigenvalues(nvars, d, states, lambda_pointwise)
integer, intent(in) :: nvars
integer, intent(in) :: d
- real, intent(in) :: states(nvars, -2:3)
- real, intent(out) :: lambda_pointwise(nvars, -2:3)
+ real*8, intent(in) :: states(nvars, -2:3)
+ real*8, intent(out) :: lambda_pointwise(nvars, -2:3)
- real u(3), c, p, rho, ke
+ real*8 u(3), c, p, rho, ke
integer v, k, i
do k=-2,3
rho = states(1,k)
- ke = 0.0
+ ke = 0.0d0
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)
+ p = (1.4d0 - 1)*(states(nvars,k) - 0.5d0*rho*ke)
+ c = sqrt(1.4d0*p/rho)
do v=1,nvars-2
lambda_pointwise(v,k) = u(d)
@@ -280,21 +280,21 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
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)
+ 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)
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)
+ 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)
integer i, j
- real b1, b2, alpha, beta
+ real*8 b1, b2, alpha, beta
ik = d
il = d+1
@@ -309,17 +309,17 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
end do
do j=1,2
- ke = 0.0
+ ke = 0.0d0
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)
+ 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
- metric_norm(i) = 0.0
+ metric_norm(i) = 0.0d0
do j=1,ndim
metric_norm(i) = metric_norm(i) + metrics_frozen(i,j)**2
end do
@@ -342,33 +342,33 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
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
+ q = 0.0d0
do i=1,ndim
q = q + u(i)**2
end do
- c = sqrt((1.4 - 1.0)*(H - 0.5*q))
+ c = sqrt((1.4d0 - 1.0d0)*(H - 0.5d0*q))
- b1 = (1.4 - 1.0)/(c**2)
- b2 = 1.0 + b1*q - b1*H
+ b1 = (1.4d0 - 1.0d0)/(c**2)
+ b2 = 1.0d0 + b1*q - b1*H
- u_tilde(1) = 0.0
+ u_tilde(1) = 0.0d0
do i=1,ndim
u_tilde(1) = u_tilde(1) + k(i)*u(i)
end do
- u_tilde(2) = 0.0
+ u_tilde(2) = 0.0d0
do i=1,ndim
u_tilde(2) = u_tilde(2) + l(i)*u(i)
end do
- u_tilde(3) = 0.0
+ u_tilde(3) = 0.0d0
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)
+ 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)
@@ -376,19 +376,19 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
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)
- R(1,1) = 1.0
+ R(1,1) = 1.0d0
R(2,1) = u(1)
R(3,1) = u(2)
R(4,1) = u(3)
- R(5,1) = H - 1.0/b1
+ R(5,1) = H - 1.0d0/b1
- R(1,2) = 0.0
+ R(1,2) = 0.0d0
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(1,3) = 0.0d0
R(2,3) = rho*m(1)
R(3,3) = rho*m(2)
R(4,3) = rho*m(3)
@@ -406,7 +406,7 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
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,1) = 1.0d0 - b2
R_inv(1,2) = b1*u(1)
R_inv(1,3) = b1*u(2)
R_inv(1,4) = b1*u(3)
@@ -416,13 +416,13 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
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(2,5) = 0.0d0
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(3,5) = 0.0d0
R_inv(4,1) = beta*(b2 - u_tilde(1)/c)
R_inv(4,2) = -beta*(b1*u(1) - k(1)/c)
@@ -441,15 +441,15 @@ subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
implicit none
integer, intent(in) :: nvars
- real, intent(in) :: lambda_pointwise(nvars, -2:3)
- real, intent(in) :: lambda_roe(nvars)
- real, intent(out) :: wavespeeds(nvars)
+ real*8, intent(in) :: lambda_pointwise(nvars, -2:3)
+ real*8, intent(in) :: lambda_roe(nvars)
+ real*8, intent(out) :: wavespeeds(nvars)
- real kappa
+ real*8 kappa
integer v
integer k
- kappa = 1.1
+ kappa = 1.1d0
do v=1,nvars
wavespeeds(v) = abs(lambda_roe(v))
@@ -471,32 +471,32 @@ subroutine convert_to_generalized_frozen( &
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*8, intent(in) :: states(nvars, -2:3)
+ real*8, intent(in) :: fluxes(nvars, ndim, -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, ndim, -2:3)
- real jacobian_frozen
+ real*8 jacobian_frozen
integer c, k, v
integer i, j
do j=1,ndim
do i=1,ndim
- metrics_frozen(i,j) = 0.0
+ metrics_frozen(i,j) = 0.0d0
do c=1,2
- metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5*metrics(i,j,c)/metric_jacobians(c)
+ metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5d0*metrics(i,j,c)/metric_jacobians(c)
end do
end do
end do
- jacobian_frozen = 0.5*(metric_jacobians(1) + metric_jacobians(2))
+ jacobian_frozen = 0.5d0*(metric_jacobians(1) + metric_jacobians(2))
do i=1,ndim
- combined_frozen_metrics(i) = 0.0
+ combined_frozen_metrics(i) = 0.0d0
do j=1,ndim
combined_frozen_metrics(i) = combined_frozen_metrics(i) + metrics_frozen(i,j)**2
end do
@@ -510,7 +510,7 @@ subroutine convert_to_generalized_frozen( &
do k=-2,3
do v=1,nvars
- call mult_mat_vec(ndim, ndim, 1.0, metrics_frozen, &
+ call mult_mat_vec(ndim, ndim, 1.0d0, metrics_frozen, &
fluxes(v,:,k), generalized_fluxes_frozen(v,:,k))
end do
end do
@@ -525,10 +525,10 @@ subroutine convert_to_generalized(nvars, ndim, nx, ny, nz, &
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)
+ real*8, intent(in) :: fluxes(nvars, ndim, -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) :: generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
integer i, j, k, v
@@ -536,7 +536,7 @@ subroutine convert_to_generalized(nvars, ndim, nx, ny, nz, &
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), &
+ call mult_mat_vec(ndim, ndim, 1.0d0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
fluxes(v,:,i,j,k), generalized_fluxes(v,:,i,j,k))
end do
end do
@@ -552,10 +552,10 @@ subroutine convert_from_generalized(nvars, ndim, nx, ny, nz, &
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)
+ real*8, intent(in) :: flux_derivatives_generalized(nvars, ndim, -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, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
integer i, j, k, v
integer d
@@ -579,11 +579,11 @@ end subroutine
subroutine solve3x3(a, b, x)
implicit none
- real, intent(in) :: a(3,3), b(3)
- real, intent(out) :: x(3)
+ real*8, intent(in) :: a(3,3), b(3)
+ real*8, intent(out) :: x(3)
- real temp(3,3)
- real det_a(1), det_temp(1)
+ real*8 temp(3,3)
+ real*8 det_a(1), det_temp(1)
integer k
integer i,j
@@ -603,8 +603,8 @@ end subroutine
subroutine determinant3x3(a, det)
implicit none
- real, intent(in) :: a(3,3)
- real, intent(out) :: det(1)
+ real*8, intent(in) :: a(3,3)
+ real*8, 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)) &
@@ -621,26 +621,26 @@ subroutine split_characteristic_fluxes(nvars, &
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) :: wavespeeds(nvars)
- real, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
- real, intent(out) :: characteristic_fluxes_neg(nvars, -2:3)
+ 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)
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
+ 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) &
- + 0.5*R_inv(m,l) &
+ + 0.5d0*R_inv(m,l) &
*(generalized_fluxes_frozen(l,k) + wavespeeds(m)*generalized_states_frozen(l,k))
characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) &
- + 0.5*R_inv(m,l) &
+ + 0.5d0*R_inv(m,l) &
*(generalized_fluxes_frozen(l,k) - wavespeeds(m)*generalized_states_frozen(l,k))
end do
end do
@@ -652,16 +652,16 @@ subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, &
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)
+ 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)
+ real*8, intent(in) :: combined_frozen_metrics
+ 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)
integer v
call consistent_part(nvars, generalized_fluxes, consistent)
@@ -679,8 +679,8 @@ 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)
+ real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
+ real*8, intent(out) :: consistent(nvars)
integer v
@@ -688,7 +688,7 @@ subroutine consistent_part(nvars, generalized_fluxes, consistent)
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.0
+ consistent(v) = consistent(v)/60.0d0
end do
end subroutine
@@ -697,16 +697,16 @@ subroutine dissipation_part_pos(nvars, characteristic_fluxes, &
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*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(in) :: combined_frozen_metrics
+ real*8, intent(in) :: R(nvars, nvars)
+ real*8, intent(out) :: dissipation_pos(nvars)
- real combined_fluxes(nvars)
+ real*8 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)
+ call mult_mat_vec(nvars, nvars, -1.0d0/60, R, combined_fluxes, dissipation_pos)
end subroutine
subroutine weno_combination_pos(nvars, characteristic_fluxes, &
@@ -714,12 +714,12 @@ subroutine weno_combination_pos(nvars, characteristic_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*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(in) :: combined_frozen_metrics
+ real*8, intent(out) :: combined_fluxes(nvars)
- real w(nvars, 3)
- real flux_differences(nvars, 3)
+ real*8 w(nvars, 3)
+ real*8 flux_differences(nvars, 3)
integer v
call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, w)
@@ -736,66 +736,112 @@ subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric
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*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(in) :: combined_frozen_metrics
+ real*8, 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)
+ real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
+ real*8 IS(nvars, 3), alpha(3), eps
integer p, i, j
- real sum_alpha(1)
+ real*8 sum_alpha(1)
- eps = 1e-6*combined_frozen_metrics
+ eps = 1.0d-6!*combined_frozen_metrics
p = 2
- do i=1,nvars
- !w1sum(1) = characteristic_fluxes(i,-2) - 4*characteristic_fluxes(i,-1) &
- ! + 3*characteristic_fluxes(i,0)
- w2sum(1) = -characteristic_fluxes(i,-1) - characteristic_fluxes(i,1) !!! BUG !!!
- !w3sum(1) = -3*characteristic_fluxes(i,0) + 4*characteristic_fluxes(i,1) &
- ! - characteristic_fluxes(i,2)
- !c1sum(1) = characteristic_fluxes(i,-2) &
- ! - 2*characteristic_fluxes(i,-1) + characteristic_fluxes(i,0)
- !c2sum(1) = characteristic_fluxes(i,-1) &
- ! - 2*characteristic_fluxes(i,0) + characteristic_fluxes(i,1)
- c3sum(1) = characteristic_fluxes(i,0) &
- - 2*characteristic_fluxes(i,1) + characteristic_fluxes(i,2)
-
- 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
+ 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(j) + eps)**p
+ 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)
+ !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)
+
+ integer i
+ integer j
+
+ do i=1,nvars
+ w1sum(1) = 0.0d0
+ do j=0,2
+ w1sum(1) = w1sum(1) + weights1(j)*characteristic_fluxes(i,j-2)
+ end do
+ w2sum(1) = 0.0d0
+ do j=0,2
+ w2sum(1) = w2sum(1) + weights2(j)*characteristic_fluxes(i,j-1)
+ end do
+ !w3sum(1) = 0.0d0
+ !do j=0,2
+ ! w3sum(1) = w3sum(1) + weights3(j)*characteristic_fluxes(i,j)
+ !end do
+ !call weighted_sum_pos(3, weights3, characteristic_fluxes(i,0:2), w3sum)
+ call weighted_sum_pos(3, weights3, characteristic_fluxes(i,:), w3sum)
+
+ c1sum(1) = 0.0d0
+ do j=0,2
+ c1sum(1) = c1sum(1) + weightsc(j)*characteristic_fluxes(i,j-2)
+ end do
+ c2sum(1) = 0.0d0
+ do j=0,2
+ c2sum(1) = c2sum(1) + weightsc(j)*characteristic_fluxes(i,j-1)
+ end do
+ c3sum(1) = 0.0d0
+ do j=0,2
+ c3sum(1) = c3sum(1) + weightsc(j)*characteristic_fluxes(i,j)
+ end do
+
+ !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)
+
+ 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
+
+subroutine weighted_sum_pos(n, w, a, a_sum)
+ integer, intent(in) :: n
+ real*8, intent(in) :: w(0:2)
+ !real*8, intent(in) :: a(0:2)
+ real*8, intent(in) :: a(-2:3)
+ real*8, intent(out) :: a_sum(1)
+
+ integer j
+
+ a_sum(1) = 0.0d0
+ do j=0,2
+ a_sum(1) = a_sum(1) + w(j)*a(j)
+ 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)
+ real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(out) :: flux_differences(nvars, 3)
integer i, v
@@ -812,16 +858,16 @@ subroutine dissipation_part_neg(nvars, characteristic_fluxes, &
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*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(in) :: combined_frozen_metrics
+ real*8, intent(in) :: R(nvars, nvars)
+ real*8, intent(out) :: dissipation_neg(nvars)
- real combined_fluxes(nvars)
+ real*8 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)
+ call mult_mat_vec(nvars, nvars, 1.0d0/60, R, combined_fluxes, dissipation_neg)
end subroutine
subroutine weno_combination_neg(nvars, characteristic_fluxes, &
@@ -829,12 +875,12 @@ subroutine weno_combination_neg(nvars, characteristic_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*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(in) :: combined_frozen_metrics
+ real*8, intent(out) :: combined_fluxes(nvars)
- real w(nvars, 3)
- real flux_differences(nvars, 3)
+ real*8 w(nvars, 3)
+ real*8 flux_differences(nvars, 3)
integer v
call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, w)
@@ -851,25 +897,48 @@ subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric
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*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(in) :: combined_frozen_metrics
+ real*8, intent(out) :: w(nvars, 3)
- real IS(3), alpha(3), eps
- real w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
+ real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
+ real*8 IS(nvars, 3), alpha(3), eps
integer p, i, j
- real sum_alpha(1)
+ real*8 sum_alpha(1)
- eps = 1e-6*combined_frozen_metrics
+ eps = 1.0d-6!*combined_frozen_metrics
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/)
+
+ real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
+
+ integer i
+
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)
@@ -878,17 +947,9 @@ subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric
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
+ 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
@@ -896,8 +957,8 @@ 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)
+ real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+ real*8, intent(out) :: flux_differences(nvars, 3)
integer i, v
@@ -913,16 +974,16 @@ 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)
+ real*8, intent(in) :: alpha
+ real*8, intent(in) :: a(m, n)
+ real*8, intent(in) :: b(n)
+ real*8, intent(out) :: c(m)
integer i, j
- real accumulator
+ real*8 accumulator
do i=1,m
- accumulator = 0.0
+ accumulator = 0.0d0
do j=1,n
accumulator = accumulator + alpha*a(i,j)*b(j)
end do
@@ -934,12 +995,12 @@ subroutine sum_array(n, a, a_sum)
implicit none
integer, intent(in) :: n
- real, intent(in) :: a(n)
- real, intent(out) :: a_sum(1)
+ real*8, intent(in) :: a(n)
+ real*8, intent(out) :: a_sum(1)
integer i
- a_sum(1) = 0.0
+ a_sum(1) = 0.0d0
do i=1,n
a_sum(1) = a_sum(1) + a(i)
end do
@@ -949,12 +1010,12 @@ 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)
+ real*8, intent(in) :: a(n), w(n)
+ real*8, intent(out) :: a_sum(1)
integer i
- a_sum(1) = 0.0
+ a_sum(1) = 0.0d0
do i=1,n
a_sum(1) = a_sum(1) + w(i)*a(i)
end do
diff --git a/data_for_test.py b/data_for_test.py
index c83dd53e00786939bcf71b2cf9ec6d38936a3fc4..bde2cd974fb711349107e23ebf51e308b3c05316 100644
--- a/data_for_test.py
+++ b/data_for_test.py
@@ -12,6 +12,7 @@ class FluxDataSingle:
def __init__(self, states_str, fluxes_str,
lam_str, R_str, R_inv_str, wavespeeds_str,
char_fluxes_pos_str, char_fluxes_neg_str,
+ oscillation_pos_str, oscillation_neg_str,
weno_weights_pos_str, weno_weights_neg_str,
consistent_str, dissipation_pos_str, dissipation_neg_str,
weno_flux_str, direction):
@@ -40,6 +41,9 @@ class FluxDataSingle:
self.char_fluxes_neg = u.expand_to_6(
u.transposed_array_from_string(char_fluxes_neg_str))
+ self.oscillation_pos = u.transposed_array_from_string(oscillation_pos_str)
+ self.oscillation_neg = u.transposed_array_from_string(oscillation_neg_str)
+
self.weno_weights_pos = u.transposed_array_from_string(weno_weights_pos_str)
self.weno_weights_neg = u.transposed_array_from_string(weno_weights_neg_str)
@@ -76,6 +80,119 @@ class FluxDataSingle:
single_data = {}
+# {{{ Input data: "flat" case
+
+single_data["Case flat:x"] = FluxDataSingle(
+ states_str="1 1 1 1 5.5,1 1 1 1 5.5",
+ fluxes_str="1 2.6 1 1 7.1,1 2.6 1 1 7.1",
+ lam_str=("1. 1. 1. 2.4966629547095764 -0.49666295470957644,"
+ "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+ R_str=("1 0 0 0.3340765523905305 0.3340765523905305,"
+ "1 0 0 0.8340765523905306 -0.1659234476094695,"
+ "1 1 0 0.3340765523905305 0.3340765523905305,"
+ "1 0 1 0.3340765523905305 0.3340765523905305,"
+ "1.5000000000000009 1 1 2.8719435219727663 1.8719435219727667"),
+ R_inv_str=("0.7321428571428572 0.1785714285714286 0.1785714285714286 "
+ "0.1785714285714286 -0.1785714285714286,"
+ "-1. 0 1. 0 0,"
+ "-1. 0 0 1. 0,"
+ "-0.5991081371313636 0.7327387580875756 -0.2672612419124244 "
+ "-0.2672612419124244 0.2672612419124244,"
+ "1.4008918628686364 -1.2672612419124245 -0.2672612419124244 "
+ "-0.2672612419124244 0.2672612419124244"),
+ wavespeeds_str="1.1 1.1 1.1 2.7463292501805343 0.5463292501805341",
+ char_fluxes_pos_str=(
+ "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326,"
+ "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326"),
+ char_fluxes_neg_str=(
+ "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178,"
+ "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178"),
+ oscillation_pos_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+ oscillation_neg_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+ weno_weights_pos_str=("0.1 0.1 0.1 0.1 0.1,"
+ "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+ weno_weights_neg_str=("0.1 0.1 0.1 0.1 0.1,"
+ "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+ consistent_str="1 2.6 1 1 7.1",
+ dissipation_pos_str=("0 0 0 0 0"),
+ dissipation_neg_str=("0 0 0 0 0"),
+ weno_flux_str=("1 2.6 1 1 7.1"),
+ direction="x")
+single_data["Case flat:y"] = FluxDataSingle(
+ states_str="1 1 1 1 5.5,1 1 1 1 5.5",
+ fluxes_str="1 2.6 1 1 7.1,1 2.6 1 1 7.1",
+ lam_str=("1. 1. 1. 2.4966629547095764 -0.49666295470957644,"
+ "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+ R_str=("1 0 0 0.3340765523905305 0.3340765523905305,"
+ "1 0 0 0.8340765523905306 -0.1659234476094695,"
+ "1 1 0 0.3340765523905305 0.3340765523905305,"
+ "1 0 1 0.3340765523905305 0.3340765523905305,"
+ "1.5000000000000009 1 1 2.8719435219727663 1.8719435219727667"),
+ R_inv_str=("0.7321428571428572 0.1785714285714286 0.1785714285714286 "
+ "0.1785714285714286 -0.1785714285714286,"
+ "-1. 0 1. 0 0,"
+ "-1. 0 0 1. 0,"
+ "-0.5991081371313636 0.7327387580875756 -0.2672612419124244 "
+ "-0.2672612419124244 0.2672612419124244,"
+ "1.4008918628686364 -1.2672612419124245 -0.2672612419124244 "
+ "-0.2672612419124244 0.2672612419124244"),
+ wavespeeds_str="1.1 1.1 1.1 2.7463292501805343 0.5463292501805341",
+ char_fluxes_pos_str=(
+ "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326,"
+ "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326"),
+ char_fluxes_neg_str=(
+ "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178,"
+ "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178"),
+ oscillation_pos_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+ oscillation_neg_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+ weno_weights_pos_str=("0.1 0.1 0.1 0.1 0.1,"
+ "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+ weno_weights_neg_str=("0.1 0.1 0.1 0.1 0.1,"
+ "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+ consistent_str="1 2.6 1 1 7.1",
+ dissipation_pos_str=("0 0 0 0 0"),
+ dissipation_neg_str=("0 0 0 0 0"),
+ weno_flux_str=("1 2.6 1 1 7.1"),
+ direction="y")
+single_data["Case flat:z"] = FluxDataSingle(
+ states_str="1 1 1 1 5.5,1 1 1 1 5.5",
+ fluxes_str="1 2.6 1 1 7.1,1 2.6 1 1 7.1",
+ lam_str=("1. 1. 1. 2.4966629547095764 -0.49666295470957644,"
+ "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+ R_str=("1 0 0 0.3340765523905305 0.3340765523905305,"
+ "1 0 0 0.8340765523905306 -0.1659234476094695,"
+ "1 1 0 0.3340765523905305 0.3340765523905305,"
+ "1 0 1 0.3340765523905305 0.3340765523905305,"
+ "1.5000000000000009 1 1 2.8719435219727663 1.8719435219727667"),
+ R_inv_str=("0.7321428571428572 0.1785714285714286 0.1785714285714286 "
+ "0.1785714285714286 -0.1785714285714286,"
+ "-1. 0 1. 0 0,"
+ "-1. 0 0 1. 0,"
+ "-0.5991081371313636 0.7327387580875756 -0.2672612419124244 "
+ "-0.2672612419124244 0.2672612419124244,"
+ "1.4008918628686364 -1.2672612419124245 -0.2672612419124244 "
+ "-0.2672612419124244 0.2672612419124244"),
+ wavespeeds_str="1.1 1.1 1.1 2.7463292501805343 0.5463292501805341",
+ char_fluxes_pos_str=(
+ "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326,"
+ "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326"),
+ char_fluxes_neg_str=(
+ "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178,"
+ "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178"),
+ oscillation_pos_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+ oscillation_neg_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+ weno_weights_pos_str=("0.1 0.1 0.1 0.1 0.1,"
+ "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+ weno_weights_neg_str=("0.1 0.1 0.1 0.1 0.1,"
+ "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+ consistent_str="1 2.6 1 1 7.1",
+ dissipation_pos_str=("0 0 0 0 0"),
+ dissipation_neg_str=("0 0 0 0 0"),
+ weno_flux_str=("1 2.6 1 1 7.1"),
+ direction="z")
+
+# }}}
+
# {{{ Input data: Case (a)
single_data["Case a:x"] = FluxDataSingle(
@@ -108,6 +225,16 @@ single_data["Case a:x"] = FluxDataSingle(
"-0.05857864376269045 -0.7274934638648571 -0.30602629876675014,"
"-0.06722315769446469 0.24852813742385726 "
"0.24852813742385726 -0.10725061063772967 -0.37456222716537935"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,"
+ "0.5180684032155981 4.777392983773268 4.777392983773268 "
+ "48.69888276854567 0.00844206573002786,"
+ "1.2951710080389953 11.94348245943317 11.94348245943317 "
+ "121.7472069213642 0.02110516432506965"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,"
+ "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+ "0.5129349293057707 0.006262897975282706,"
+ "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+ "1.2823373232644266 0.015657244938206753"),
weno_weights_pos_str=("0.999999997585736 0.9999999999716082 "
"0.9999999999716082 0.9999999999997268 0.9999909282503908,"
"2.2354258683870077e-9 2.6288602366075745e-11 "
@@ -158,6 +285,16 @@ single_data["Case a:y"] = FluxDataSingle(
"-0.05857864376269045 -0.7274934638648571 -0.30602629876675014,"
"-0.06722315769446469 0.24852813742385726 "
"0.24852813742385726 -0.10725061063772967 -0.37456222716537935"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,"
+ "0.5180684032155981 4.777392983773268 4.777392983773268 "
+ "48.69888276854567 0.00844206573002786,"
+ "1.2951710080389953 11.94348245943317 11.94348245943317 "
+ "121.7472069213642 0.02110516432506965"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,"
+ "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+ "0.5129349293057707 0.006262897975282706,"
+ "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+ "1.2823373232644266 0.015657244938206753"),
weno_weights_pos_str=("0.999999997585736 0.9999999999716082 "
"0.9999999999716082 0.9999999999997268 0.9999909282503908,"
"2.2354258683870077e-9 2.6288602366075745e-11 "
@@ -208,6 +345,16 @@ single_data["Case a:z"] = FluxDataSingle(
"-0.05857864376269045 -0.7274934638648571 -0.30602629876675014,"
"-0.06722315769446469 0.24852813742385726 "
"0.24852813742385726 -0.10725061063772967 -0.37456222716537935"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,"
+ "0.5180684032155981 4.777392983773268 4.777392983773268 "
+ "48.69888276854567 0.00844206573002786,"
+ "1.2951710080389953 11.94348245943317 11.94348245943317 "
+ "121.7472069213642 0.02110516432506965"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,"
+ "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+ "0.5129349293057707 0.006262897975282706,"
+ "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+ "1.2823373232644266 0.015657244938206753"),
weno_weights_pos_str=("0.999999997585736 0.9999999999716082 "
"0.9999999999716082 0.9999999999997268 0.9999909282503908,"
"2.2354258683870077e-9 2.6288602366075745e-11 "
@@ -264,6 +411,16 @@ single_data["Case b:x"] = FluxDataSingle(
"-0.6627416997969526 -0.3125273038443238 -1.4795302623674356,"
"-1.0907156303492833 1.2301515190164993 "
"1.2301515190164993 -0.23295627081394965 -7.523052586013577"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,"
+ "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+ "0.006262897975282706 0.5129349293057707,"
+ "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+ "0.015657244938206753 1.2823373232644266"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,"
+ "0.5180684032155981 4.777392983773268 4.777392983773268 "
+ "0.00844206573002786 48.69888276854567,"
+ "1.2951710080389953 11.94348245943317 11.94348245943317 "
+ "0.02110516432506965 121.7472069213642"),
weno_weights_pos_str=("0.9999965203327451 0.999999959029252 "
"0.999999959029252 0.9999835300215946 0.9999999975371708,"
"3.2217041402391125e-6 3.793560951911702e-8 "
@@ -315,6 +472,16 @@ single_data["Case b:y"] = FluxDataSingle(
"-0.6627416997969526 -0.3125273038443238 -1.4795302623674356,"
"-1.0907156303492833 1.2301515190164993 "
"1.2301515190164993 -0.23295627081394965 -7.523052586013577"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,"
+ "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+ "0.006262897975282706 0.5129349293057707,"
+ "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+ "0.015657244938206753 1.2823373232644266"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,"
+ "0.5180684032155981 4.777392983773268 4.777392983773268 "
+ "0.00844206573002786 48.69888276854567,"
+ "1.2951710080389953 11.94348245943317 11.94348245943317 "
+ "0.02110516432506965 121.7472069213642"),
weno_weights_pos_str=("0.9999965203327451 0.999999959029252 "
"0.999999959029252 0.9999835300215946 0.9999999975371708,"
"3.2217041402391125e-6 3.793560951911702e-8 "
@@ -366,6 +533,16 @@ single_data["Case b:z"] = FluxDataSingle(
"-0.6627416997969526 -0.3125273038443238 -1.4795302623674356,"
"-1.0907156303492833 1.2301515190164993 "
"1.2301515190164993 -0.23295627081394965 -7.523052586013577"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,"
+ "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+ "0.006262897975282706 0.5129349293057707,"
+ "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+ "0.015657244938206753 1.2823373232644266"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,"
+ "0.5180684032155981 4.777392983773268 4.777392983773268 "
+ "0.00844206573002786 48.69888276854567,"
+ "1.2951710080389953 11.94348245943317 11.94348245943317 "
+ "0.02110516432506965 121.7472069213642"),
weno_weights_pos_str=("0.9999965203327451 0.999999959029252 "
"0.999999959029252 0.9999835300215946 0.9999999975371708,"
"3.2217041402391125e-6 3.793560951911702e-8 "
@@ -421,6 +598,14 @@ single_data["Case c:x"] = FluxDataSingle(
"-0.1757359312880714 -0.8893252542028551 -0.20004041758408686,"
"-0.009488213714461069 0.49705627484771453 "
"0.7455844122715716 -0.4306378813126712 -0.31431832174320373"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,0.9757858752013722 "
+ "19.109571935093072 42.996536853959384 46.95775189047701 "
+ "0.007062155488717821,2.4394646880034303 47.77392983773268 "
+ "107.49134213489849 117.39437972619254 0.017655388721794552"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,0.025685091003327314 "
+ "0.5030110669373978 1.1317749006091449 0.28052547473186484 "
+ "0.017412585838667068,0.06421272750831829 1.2575276673434945 "
+ "2.8294372515228616 0.7013136868296621 0.04353146459666767"),
weno_weights_pos_str=("0.999999999319454 0.9999999999982254 "
"0.9999999999996494 0.9999999999997061 0.9999870425235151,"
"6.301345524776077e-10 1.6430428067143077e-12 "
@@ -471,6 +656,14 @@ single_data["Case c:y"] = FluxDataSingle(
"-0.1171572875253809 -1.6097440213843948 -0.5689873221894901,"
"-0.018976427428922138 1.4911688245431431 "
"0.49705627484771453 -0.05753157056903602 -1.432380835542713"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,3.9031435008054665 "
+ "171.98614741583754 19.109571935093072 282.0389435835765 "
+ "10.607678229749117,9.757858752013666 429.96536853959395 "
+ "47.77392983773268 705.0973589589414 26.519195574372795"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,0.10274036401330869 "
+ "4.5270996024365795 0.5030110669373978 3.21248465662163 "
+ "0.9939311452005626,0.25685091003327176 11.317749006091447 "
+ "1.2575276673434945 8.031211641554076 2.484827863001407"),
weno_weights_pos_str=("0.9999999999574652 0.999999999999978 "
"0.9999999999982254 0.9999999999999919 0.9999999999942413,"
"3.9384014966385437e-11 2.0284497966079935e-14 "
@@ -521,6 +714,14 @@ single_data["Case c:z"] = FluxDataSingle(
"-0.3514718625761428 -2.447320076091316 -0.8207769392695052,"
"-0.02846464114338254 0.7455844122715716 "
"1.4911688245431431 0.8126310150223119 -3.0474996241899346"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,8.782072876812371 "
+ "42.9965368539594 171.9861474158376 967.6396764339121 "
+ "116.6680636935429,21.955182192030932 107.49134213489847 "
+ "429.9653685395939 2419.099191084781 291.67015923385725"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,0.23116581902994632 "
+ "1.1317749006091449 4.5270996024365795 14.169708155270577 "
+ "6.6110585540523275,0.5779145475748658 2.8294372515228616 "
+ "11.317749006091447 35.42427038817644 16.527646385130815"),
weno_weights_pos_str=("0.999999999991598 0.9999999999996494 "
"0.999999999999978 0.9999999999999993 0.9999999999999524,"
"7.779580658268568e-12 3.2455185423228674e-13 "
@@ -576,6 +777,14 @@ single_data["Case d:x"] = FluxDataSingle(
"-1.9882250993908572 -0.5017887559695788 -2.052422890589809,"
"-1.1162114029572359 2.4603030380329987 "
"3.690454557049497 -0.4290108979594782 -7.986924884679773"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,0.025685091003327314 "
+ "0.5030110669373978 1.1317749006091449 0.017412585838667068 "
+ "0.28052547473186484,0.06421272750831829 1.2575276673434945 "
+ "2.8294372515228616 0.04353146459666767 0.7013136868296621"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,0.9757858752013722 "
+ "19.109571935093072 42.996536853959384 0.007062155488717821 "
+ "46.95775189047701,2.4394646880034303 47.77392983773268 "
+ "107.49134213489849 0.017655388721794552 117.39437972619254"),
weno_weights_pos_str=("0.9999990185022956 0.9999999974390363 "
"0.9999999994941201 0.9999978651320311 0.9999999917661908,"
"9.08762722250494e-7 2.3712585027633853e-9 "
@@ -626,6 +835,14 @@ single_data["Case d:y"] = FluxDataSingle(
"-1.3254833995939053 -3.767119678640133 -1.3413036144786457,"
"-2.2324228059144673 7.380909114098994 "
"2.4603030380329987 -0.9465242322295992 -15.885347333048884"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,0.10274036401330869 "
+ "4.5270996024365795 0.5030110669373978 0.9939311452005626 "
+ "3.21248465662163,0.25685091003327176 11.317749006091447 "
+ "1.2575276673434945 2.484827863001407 8.031211641554076"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,3.9031435008054665 "
+ "171.98614741583754 19.109571935093072 10.607678229749117 "
+ "282.0389435835765,9.757858752013666 429.96536853959395 "
+ "47.77392983773268 26.519195574372795 705.0973589589414"),
weno_weights_pos_str=("0.9999999386221037 0.9999999999683822 "
"0.9999999974390363 0.9999999993440752 0.9999999999372101,"
"5.68308936788955e-8 2.9275831062158654e-11 "
@@ -676,6 +893,14 @@ single_data["Case d:z"] = FluxDataSingle(
"-3.9764501987817162 -8.357934000904585 0.6952990612264269,"
"-3.348634208871708 3.690454557049497 "
"7.380909114098994 0.9962654715332988 -26.24407281945101"),
+ oscillation_pos_str=("0. 0. 0. 0. 0.,0.23116581902994632 "
+ "1.1317749006091449 4.5270996024365795 6.6110585540523275 "
+ "14.169708155270577,0.5779145475748658 2.8294372515228616 "
+ "11.317749006091447 16.527646385130815 35.42427038817644"),
+ oscillation_neg_str=("0. 0. 0. 0. 0.,8.782072876812371 "
+ "42.9965368539594 171.9861474158376 116.6680636935429 "
+ "967.6396764339121,21.955182192030932 107.49134213489847 "
+ "429.9653685395939 291.67015923385725 2419.099191084781"),
weno_weights_pos_str=("0.9999999878747177 0.9999999994941201 "
"0.9999999999683822 0.9999999999851737 0.9999999999967727,"
"1.1227070122726888e-8 4.684070887233427e-10 "
@@ -701,6 +926,7 @@ single_data["Case d:z"] = FluxDataSingle(
@pytest.fixture(scope="session", params=[
+ "Case flat:x", "Case flat:y", "Case flat:z",
"Case a:x", "Case a:y", "Case a:z",
"Case b:x", "Case b:y", "Case b:z",
"Case c:x", "Case c:y", "Case c:z",
diff --git a/test.py b/test.py
index 183d89baf720e19bee8e5082136739f1cceb9401..ed0a4360ff2e3963dd1c1f4913f49ce72fbe08b6 100644
--- a/test.py
+++ b/test.py
@@ -18,7 +18,6 @@ import utilities as u
from data_for_test import flux_test_data_fixture # noqa: F401
-@pytest.mark.slow
def test_weno_flux_uniform_grid(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -38,7 +37,6 @@ def test_weno_flux_uniform_grid(ctx_factory, flux_test_data_fixture):
u.compare_arrays(flux_dev.get(), data.weno_flux)
-@pytest.mark.slow
def test_consistent_part_uniform_grid(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -54,7 +52,6 @@ def test_consistent_part_uniform_grid(ctx_factory, flux_test_data_fixture):
u.compare_arrays(consistent_dev.get(), data.consistent)
-@pytest.mark.slow
def test_dissipation_part_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -72,7 +69,6 @@ def test_dissipation_part_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
u.compare_arrays(dissipation_dev.get(), data.dissipation_pos)
-@pytest.mark.slow
def test_dissipation_part_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -103,10 +99,12 @@ def test_weno_weights_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
combined_frozen_metrics=1.0,
w=weights_dev)
+ sum_weights = np.sum(weights_dev.get(), axis=1)
+ u.compare_arrays(sum_weights, np.ones(data.nvars))
+
u.compare_arrays(weights_dev.get(), data.weno_weights_pos)
-@pytest.mark.slow
def test_weno_weights_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -120,10 +118,42 @@ def test_weno_weights_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
combined_frozen_metrics=1.0,
w=weights_dev)
+ sum_weights = np.sum(weights_dev.get(), axis=1)
+ u.compare_arrays(sum_weights, np.ones(data.nvars))
+
u.compare_arrays(weights_dev.get(), data.weno_weights_neg)
-@pytest.mark.slow
+def test_oscillation_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
+ data = flux_test_data_fixture
+
+ prg = u.get_weno_program_with_root_kernel("oscillation_pos")
+ queue = u.get_queue(ctx_factory)
+
+ oscillation_dev = u.empty_array_on_device(queue, data.nvars, 3)
+
+ prg(queue, nvars=data.nvars,
+ characteristic_fluxes=data.char_fluxes_pos,
+ oscillation=oscillation_dev)
+
+ u.compare_arrays(oscillation_dev.get(), data.oscillation_pos)
+
+
+def test_oscillation_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
+ data = flux_test_data_fixture
+
+ prg = u.get_weno_program_with_root_kernel("oscillation_neg")
+ queue = u.get_queue(ctx_factory)
+
+ oscillation_dev = u.empty_array_on_device(queue, data.nvars, 3)
+
+ prg(queue, nvars=data.nvars,
+ characteristic_fluxes=data.char_fluxes_neg,
+ oscillation=oscillation_dev)
+
+ u.compare_arrays(oscillation_dev.get(), data.oscillation_neg)
+
+
def test_flux_splitting_uniform_grid(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -145,7 +175,6 @@ def test_flux_splitting_uniform_grid(ctx_factory, flux_test_data_fixture):
u.compare_arrays(fluxes_neg_dev.get(), data.char_fluxes_neg)
-@pytest.mark.slow
def test_pointwise_eigenvalues_ideal_gas(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
@@ -160,12 +189,11 @@ def test_pointwise_eigenvalues_ideal_gas(ctx_factory, flux_test_data_fixture):
u.compare_arrays(lam_dev.get(), data.lam_pointwise)
-@pytest.mark.slow
def test_roe_uniform_grid_ideal_gas(ctx_factory, flux_test_data_fixture):
data = flux_test_data_fixture
def identity_matrix(n):
- return np.identity(n).astype(np.float32).copy(order="F")
+ return np.identity(n).astype(np.float64).copy(order="F")
def check_roe_identity(states, R, R_inv):
d_state = states[:, 1] - states[:, 0]
@@ -200,7 +228,6 @@ def test_roe_uniform_grid_ideal_gas(ctx_factory, flux_test_data_fixture):
check_roe_property(data.state_pair, data.flux_pair, R, R_inv, lam)
-@pytest.mark.slow
@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [
("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
@@ -226,7 +253,6 @@ def test_lax_wavespeeds(
u.compare_arrays(lam_dev.get(), lam_expected)
-@pytest.mark.slow
def test_matvec(ctx_factory):
prg = u.get_weno_program_with_root_kernel("mult_mat_vec")
queue = u.get_queue(ctx_factory)
diff --git a/utilities.py b/utilities.py
index 129680a3923cd1a392bb1c8be137d16be0680f4e..c8bde8a3b7baf9e034541188c01b6e73e3bfc00e 100644
--- a/utilities.py
+++ b/utilities.py
@@ -21,7 +21,7 @@ def random_array_on_device(queue, *shape):
def empty_array_on_device(queue, *shape):
- return cl.array.empty(queue, shape, dtype=np.float32, order="F")
+ return cl.array.empty(queue, shape, dtype=np.float64, order="F")
def arrays_from_string(string_arrays):
@@ -38,7 +38,7 @@ def array_from_string(string_array):
return np.array(split_map_to_list(string_array[1:], int, " "))
else:
return np.array(
- split_map_to_list(string_array, float, " "), dtype=np.float32)
+ split_map_to_list(string_array, float, " "), dtype=np.float64)
def array_from_string_2d(string_array):
if string_array[0] == ",":