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] == ",":