Newer
Older
subroutine compute_flux_derivatives(states, fluxes, flux_derivatives)
real, intent(in) :: states(nvars, nx, ny, nz)
real, intent(in) :: fluxes(nvars, nx, ny, nz)
real, intent(in) :: metrics(ndim, ndim, nx, ny, nz)
real, intent(in) :: metric_jacobians(nx, ny, nz)
real, intent(out) :: flux_derivatives(ndim, nx, ny, nz)
! grid spacing in generalized coordinates
delta_xi = 1.0
delta_eta = 1.0
delta_zeta = 1.0
flux_derivatives_generalized = 0.0
do k=1,nz
do j=1,ny
do i=0,nx
call weno_flux(states(:, i-2:i+3, j, k),
fluxes(:, i-2:i+3, j, k),
metrics(:, :, i-2:i+3, j, k),
metric_jacobians(i-2:i+3, j, k),
flux)
flux_derivatives_generalized(1, i, j, k) += flux / delta_xi
flux_derivatives_generalized(1, i+1, j, k) -= flux / delta_xi
end do
end do
end do
! two more loops:
! j is inner index, filling flux_derivatives_generalized(2,:,:,:), using delta_eta
! k is inner index, filling flux_derivatives_generalized(3,:,:,:), using delta_zeta
call convert_from_generalized(flux_derivatives_generalized,
metrics,
metric_jacobians,
flux_derivatives)
end subroutine
subroutine weno_flux(states, fluxes, metrics, metric_jacobians, flux)
real, intent(in) :: states(nvars, -2:3)
real, intent(in) :: fluxes(nvars, -2:3)
real, intent(in) :: metrics(ndim, ndim, -2:3)
real, intent(in) :: metric_jacobians(-2:3)
real, intent(out) :: flux
! compute frozen metrics
call roe_eigenvectors(..., R)
! compute characteristic fluxes
call split_characteristic_fluxes(..., characteristic_fluxes_pos, characteristic_fluxes_neg)
call consistent_part(fluxes, metrics, metric_jacobians, consistent)
call dissipation_part_pos(characteristic_fluxes_pos, combined_frozen_metrics, R, dissipation_pos)
call dissipation_part_neg(characteristic_fluxes_neg, combined_frozen_metrics, R, dissipation_neg)
flux = consistent + dissipation_pos + dissipation_neg
end subroutine
subroutine roe_eigenvectors(..., R)
! returns eigenvectors of Roe matrix
end subroutine
subroutine consistent_part(fluxes, metrics, metric_jacobians, consistent)
call convert_to_generalized(fluxes, metrics, metric_jacobians, generalized_fluxes)
consistent = sum([1, -8, 37, 37, -8, 1]*generalized_fluxes)/60
end subroutine
subroutine convert_to_generalized(fluxes,
metrics,
metric_jacobians,
generalized_fluxes)
real, intent(in) :: fluxes(nvars, -2:3)
real, intent(in) :: metrics(ndim, ndim, -2:3)
real, intent(in) :: metric_jacobians(-2:3)
real, intent(out) :: generalized_fluxes(-2:3)
! transforms fluxes from physical coordinates to generalized coordinates
! uses metrics at each point
end subroutine
subroutine dissipation_part_pos(states, fluxes, metrics, metric_jacobians, dissipation_pos)
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
real, intent(in) :: characteristic_fluxes(-2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(in) :: R(nvars, nvars)
call weno_combination_pos(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
dissipation_pos = -matmul(R, combined_fluxes)/60
end subroutine
subroutine weno_combination_pos(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: combined_fluxes(nvars)
real :: w(nvars, 3)
real :: flux_differences(nvars, 3)
call weno_weights_pos(characteristic_fluxes, combined_frozen_metrics, w)
call flux_differences_pos(characteristic_fluxes, flux_differences)
combined_fluxes = (20*w(:,1) - 1)*flux_differences(:,1)
- (10*(w(:,1) + w(:,2)) - 5)*flux_differences(:,2)
+ flux_differences(:,3)
end subroutine
subroutine weno_weights_pos(characteristic_fluxes, combined_frozen_metrics, w)
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: w(nvars, 3)
C = [0.1, 0.6, 0.3]
eps = 1e-6*combined_frozen_metrics
p = 2
do i=1,nvars
IS(1) = (1/4)*(sum([1, -4, 3]*characteristic_fluxes(i,-2:0)))**2
+ (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i,-2:0)))**2
IS(2) = (1/4)*(sum([-1, 0, 1]*characteristic_fluxes(i,-1:1)))**2
+ (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i,-1:1)))**2
IS(3) = (1/4)*(sum([-3, 4, -1]*characteristic_fluxes(i,0:2)))**2
+ (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i,0:2)))**2
alpha = C/(IS + eps)**p
w(i,:) = alpha/sum(alpha)
end do
end subroutine
subroutine flux_differences_pos(characteristic_fluxes, flux_differences)
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(out) :: flux_differences(nvars, 3)
do i=1,3
flux_differences(:,i) = sum([-1, 3, -3, 1]*characteristic_fluxes(:,-3+i:i))
end do
end subroutine
subroutine dissipation_part_neg(states, fluxes, metrics, metric_jacobians, dissipation_neg)
real, intent(in) :: characteristic_fluxes(-2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(in) :: R(nvars, nvars)
call weno_combination_neg(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
dissipation_neg = matmul(R, combined_fluxes)/60
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
subroutine weno_combination_neg(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: combined_fluxes(nvars)
real :: w(nvars, 3)
real :: flux_differences(nvars, 3)
call weno_weights_neg(characteristic_fluxes, combined_frozen_metrics, w)
call flux_differences_neg(characteristic_fluxes, flux_differences)
combined_fluxes = (20*w(:,1) - 1)*flux_differences(:,1)
- (10*(w(:,1) + w(:,2)) - 5)*flux_differences(:,2)
+ flux_differences(:,3)
end subroutine
subroutine weno_weights_neg(characteristic_fluxes, combined_frozen_metrics, w)
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(in) :: combined_frozen_metrics
real, intent(out) :: w(nvars, 3)
C = [0.1, 0.6, 0.3]
eps = 1e-6*combined_frozen_metrics
p = 2
do i=1,nvars
IS(1) = (1/4)*(sum([1, -4, 3]*characteristic_fluxes(i, 3:1)))**2
+ (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 3:1)))**2
IS(2) = (1/4)*(sum([-1, 0, 1]*characteristic_fluxes(i, 2:0)))**2
+ (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 2:0)))**2
IS(3) = (1/4)*(sum([-3, 4, -1]*characteristic_fluxes(i, 1:-1)))**2
+ (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 1:-1)))**2
alpha = C/(IS + eps)**p
w(i,:) = alpha/sum(alpha)
end do
end subroutine
subroutine flux_differences_neg(characteristic_fluxes, flux_differences)
real, intent(in) :: characteristic_fluxes(nvars, -2:3)
real, intent(out) :: flux_differences(nvars, 3)
do i=1,3
flux_differences(:,i) = sum([-1, 3, -3, 1]*characteristic_fluxes(:,1-i:4-i))
end do
subroutine convert_from_generalized(flux_derivatives_generalized,
metrics,
metric_jacobians,
flux_derivatives)
! transforms flux derivatives from generalized coordinates to physical coordinates
end subroutine