Neko 1.99.2
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
euler_res_sx.f90
Go to the documentation of this file.
1module euler_res_sx
2 use euler_residual, only : euler_rhs_t
3 use field, only : field_t
4 use ax_product, only : ax_t
5 use coefs, only : coef_t
6 use gather_scatter, only : gs_t
7 use num_types, only : rp
8 use operators, only: div
9 use math, only: subcol3, copy, sub2, add2, add3, &
10 col2, col3, addcol3, cmult, cfill, invcol3
11 use gs_ops, only : gs_op_add
12 use scratch_registry, only: neko_scratch_registry
13 use runge_kutta_time_scheme, only : runge_kutta_time_scheme_t
14 use field_list, only : field_list_t
15 implicit none
16 private
17
18 type, public, extends(euler_rhs_t) :: euler_res_sx_t
19 contains
20 procedure, nopass :: step => advance_primitive_variables_sx
21 procedure, nopass :: evaluate_rhs_sx
22 end type euler_res_sx_t
23
24contains
25 subroutine advance_primitive_variables_sx(rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, &
26 coef, gs, h, c_avisc_low, rk_scheme, dt)
27 type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E
28 type(field_t), intent(in) :: p, u, v, w, h
29 class(ax_t), intent(inout) :: Ax
30 type(coef_t), intent(inout) :: coef
31 type(gs_t), intent(inout) :: gs
32 real(kind=rp) :: c_avisc_low
33 class(runge_kutta_time_scheme_t), intent(in) :: rk_scheme
34 real(kind=rp), intent(in) :: dt
35 integer :: n, s, i, j, k
36 real(kind=rp) :: t, c
37 type(field_t), pointer :: k_rho_1, k_rho_2, k_rho_3, k_rho_4, &
38 k_m_x_1, k_m_x_2, k_m_x_3, k_m_x_4, &
39 k_m_y_1, k_m_y_2, k_m_y_3, k_m_y_4, &
40 k_m_z_1, k_m_z_2, k_m_z_3, k_m_z_4, &
41 k_e_1, k_e_2, k_e_3, k_e_4, &
42 temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_e
43 integer :: tmp_indices(25)
44 type(field_list_t) :: k_rho, k_m_x, k_m_y, k_m_z, k_E
45
46 n = p%dof%size()
47 s = rk_scheme%order
48 call neko_scratch_registry%request_field(k_rho_1, tmp_indices(1), .true.)
49 call neko_scratch_registry%request_field(k_rho_2, tmp_indices(2), .true.)
50 call neko_scratch_registry%request_field(k_rho_3, tmp_indices(3), .true.)
51 call neko_scratch_registry%request_field(k_rho_4, tmp_indices(4), .true.)
52 call neko_scratch_registry%request_field(k_m_x_1, tmp_indices(5), .true.)
53 call neko_scratch_registry%request_field(k_m_x_2, tmp_indices(6), .true.)
54 call neko_scratch_registry%request_field(k_m_x_3, tmp_indices(7), .true.)
55 call neko_scratch_registry%request_field(k_m_x_4, tmp_indices(8), .true.)
56 call neko_scratch_registry%request_field(k_m_y_1, tmp_indices(9), .true.)
57 call neko_scratch_registry%request_field(k_m_y_2, tmp_indices(10), .true.)
58 call neko_scratch_registry%request_field(k_m_y_3, tmp_indices(11), .true.)
59 call neko_scratch_registry%request_field(k_m_y_4, tmp_indices(12), .true.)
60 call neko_scratch_registry%request_field(k_m_z_1, tmp_indices(13), .true.)
61 call neko_scratch_registry%request_field(k_m_z_2, tmp_indices(14), .true.)
62 call neko_scratch_registry%request_field(k_m_z_3, tmp_indices(15), .true.)
63 call neko_scratch_registry%request_field(k_m_z_4, tmp_indices(16), .true.)
64 call neko_scratch_registry%request_field(k_e_1, tmp_indices(17), .true.)
65 call neko_scratch_registry%request_field(k_e_2, tmp_indices(18), .true.)
66 call neko_scratch_registry%request_field(k_e_3, tmp_indices(19), .true.)
67 call neko_scratch_registry%request_field(k_e_4, tmp_indices(20), .true.)
68 call neko_scratch_registry%request_field(temp_rho, tmp_indices(21), .false.)
69 call neko_scratch_registry%request_field(temp_m_x, tmp_indices(22), .false.)
70 call neko_scratch_registry%request_field(temp_m_y, tmp_indices(23), .false.)
71 call neko_scratch_registry%request_field(temp_m_z, tmp_indices(24), .false.)
72 call neko_scratch_registry%request_field(temp_e, tmp_indices(25), .false.)
73
74 call k_rho%init(4)
75 call k_rho%assign(1, k_rho_1)
76 call k_rho%assign(2, k_rho_2)
77 call k_rho%assign(3, k_rho_3)
78 call k_rho%assign(4, k_rho_4)
79 call k_m_x%init(4)
80 call k_m_x%assign(1, k_m_x_1)
81 call k_m_x%assign(2, k_m_x_2)
82 call k_m_x%assign(3, k_m_x_3)
83 call k_m_x%assign(4, k_m_x_4)
84 call k_m_y%init(4)
85 call k_m_y%assign(1, k_m_y_1)
86 call k_m_y%assign(2, k_m_y_2)
87 call k_m_y%assign(3, k_m_y_3)
88 call k_m_y%assign(4, k_m_y_4)
89 call k_m_z%init(4)
90 call k_m_z%assign(1, k_m_z_1)
91 call k_m_z%assign(2, k_m_z_2)
92 call k_m_z%assign(3, k_m_z_3)
93 call k_m_z%assign(4, k_m_z_4)
94 call k_e%init(4)
95 call k_e%assign(1, k_e_1)
96 call k_e%assign(2, k_e_2)
97 call k_e%assign(3, k_e_3)
98 call k_e%assign(4, k_e_4)
99
100 ! Runge-Kutta stages
101 do i = 1, s
102 call copy(temp_rho%x, rho_field%x, n)
103 call copy(temp_m_x%x, m_x%x, n)
104 call copy(temp_m_y%x, m_y%x, n)
105 call copy(temp_m_z%x, m_z%x, n)
106 call copy(temp_e%x, e%x, n)
107
108 do j = 1, i-1
109 do concurrent(k = 1:n)
110 temp_rho%x(k,1,1,1) = temp_rho%x(k,1,1,1) &
111 + dt * rk_scheme%coeffs_A(i, j) * k_rho%items(j)%ptr%x(k,1,1,1)
112 temp_m_x%x(k,1,1,1) = temp_m_x%x(k,1,1,1) &
113 + dt * rk_scheme%coeffs_A(i, j) * k_m_x%items(j)%ptr%x(k,1,1,1)
114 temp_m_y%x(k,1,1,1) = temp_m_y%x(k,1,1,1) &
115 + dt * rk_scheme%coeffs_A(i, j) * k_m_y%items(j)%ptr%x(k,1,1,1)
116 temp_m_z%x(k,1,1,1) = temp_m_z%x(k,1,1,1) &
117 + dt * rk_scheme%coeffs_A(i, j) * k_m_z%items(j)%ptr%x(k,1,1,1)
118 temp_e%x(k,1,1,1) = temp_e%x(k,1,1,1) &
119 + dt * rk_scheme%coeffs_A(i, j) * k_e%items(j)%ptr%x(k,1,1,1)
120 end do
121 end do
122
123 ! Compute f(U) = rhs(U) for the intermediate values
124 call evaluate_rhs_sx(k_rho%items(i)%ptr, k_m_x%items(i)%ptr, &
125 k_m_y%items(i)%ptr, k_m_z%items(i)%ptr, &
126 k_e%items(i)%ptr, &
127 temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_e, &
128 p, u, v, w, ax, &
129 coef, gs, h, c_avisc_low)
130 end do
131
132 ! Update the solution
133 do i = 1, s
134 do concurrent(k = 1:n)
135 rho_field%x(k,1,1,1) = rho_field%x(k,1,1,1) &
136 + dt * rk_scheme%coeffs_b(i) * k_rho%items(i)%ptr%x(k,1,1,1)
137 m_x%x(k,1,1,1) = m_x%x(k,1,1,1) &
138 + dt * rk_scheme%coeffs_b(i) * k_m_x%items(i)%ptr%x(k,1,1,1)
139 m_y%x(k,1,1,1) = m_y%x(k,1,1,1) &
140 + dt * rk_scheme%coeffs_b(i) * k_m_y%items(i)%ptr%x(k,1,1,1)
141 m_z%x(k,1,1,1) = m_z%x(k,1,1,1) &
142 + dt * rk_scheme%coeffs_b(i) * k_m_z%items(i)%ptr%x(k,1,1,1)
143 e%x(k,1,1,1) = e%x(k,1,1,1) &
144 + dt * rk_scheme%coeffs_b(i) * k_e%items(i)%ptr%x(k,1,1,1)
145 end do
146 end do
147
148 call neko_scratch_registry%relinquish_field(tmp_indices)
149 end subroutine advance_primitive_variables_sx
150
151 subroutine evaluate_rhs_sx(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_E, &
152 rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, &
153 coef, gs, h, c_avisc_low)
154 type(field_t), intent(inout) :: rhs_rho_field, rhs_m_x, &
155 rhs_m_y, rhs_m_z, rhs_E
156 type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E
157 type(field_t), intent(in) :: p, u, v, w, h
158 class(ax_t), intent(inout) :: Ax
159 type(coef_t), intent(inout) :: coef
160 type(gs_t), intent(inout) :: gs
161 real(kind=rp) :: c_avisc_low
162 integer :: i, n
163 type(field_t), pointer :: f_x, f_y, f_z, &
164 visc_rho, visc_m_x, visc_m_y, visc_m_z, visc_E
165 integer :: tmp_indices(8)
166
167 n = coef%dof%size()
168 call neko_scratch_registry%request_field(f_x, tmp_indices(1), .false.)
169 call neko_scratch_registry%request_field(f_y, tmp_indices(2), .false.)
170 call neko_scratch_registry%request_field(f_z, tmp_indices(3), .false.)
171
173 call div(rhs_rho_field%x, m_x%x, m_y%x, m_z%x, coef)
174
176 ! m_x
177 do concurrent(i = 1:n)
178 f_x%x(i,1,1,1) = m_x%x(i,1,1,1) * m_x%x(i,1,1,1) / rho_field%x(i, 1, 1, 1) &
179 + p%x(i,1,1,1)
180 f_y%x(i,1,1,1) = m_x%x(i,1,1,1) * m_y%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
181 f_z%x(i,1,1,1) = m_x%x(i,1,1,1) * m_z%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
182 end do
183 call div(rhs_m_x%x, f_x%x, f_y%x, f_z%x, coef)
184 ! m_y
185 do concurrent(i = 1:n)
186 f_x%x(i,1,1,1) = m_y%x(i,1,1,1) * m_x%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
187 f_y%x(i,1,1,1) = m_y%x(i,1,1,1) * m_y%x(i,1,1,1) / rho_field%x(i, 1, 1, 1) &
188 + p%x(i,1,1,1)
189 f_z%x(i,1,1,1) = m_y%x(i,1,1,1) * m_z%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
190 end do
191 call div(rhs_m_y%x, f_x%x, f_y%x, f_z%x, coef)
192 ! m_z
193 do concurrent(i = 1:n)
194 f_x%x(i,1,1,1) = m_z%x(i,1,1,1) * m_x%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
195 f_y%x(i,1,1,1) = m_z%x(i,1,1,1) * m_y%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
196 f_z%x(i,1,1,1) = m_z%x(i,1,1,1) * m_z%x(i,1,1,1) / rho_field%x(i, 1, 1, 1) &
197 + p%x(i,1,1,1)
198 end do
199 call div(rhs_m_z%x, f_x%x, f_y%x, f_z%x, coef)
200
202 do concurrent(i = 1:n)
203 f_x%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) &
204 * u%x(i,1,1,1)
205 f_y%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) &
206 * v%x(i,1,1,1)
207 f_z%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) &
208 * w%x(i,1,1,1)
209 end do
210 call div(rhs_e%x, f_x%x, f_y%x, f_z%x, coef)
211
212 call gs%op(rhs_rho_field, gs_op_add)
213 call gs%op(rhs_m_x, gs_op_add)
214 call gs%op(rhs_m_y, gs_op_add)
215 call gs%op(rhs_m_z, gs_op_add)
216 call gs%op(rhs_e, gs_op_add)
217
218 do concurrent(i = 1:rhs_e%dof%size())
219 rhs_rho_field%x(i,1,1,1) = rhs_rho_field%x(i,1,1,1) * coef%mult(i,1,1,1)
220 rhs_m_x%x(i,1,1,1) = rhs_m_x%x(i,1,1,1) * coef%mult(i,1,1,1)
221 rhs_m_y%x(i,1,1,1) = rhs_m_y%x(i,1,1,1) * coef%mult(i,1,1,1)
222 rhs_m_z%x(i,1,1,1) = rhs_m_z%x(i,1,1,1) * coef%mult(i,1,1,1)
223 rhs_e%x(i,1,1,1) = rhs_e%x(i,1,1,1) * coef%mult(i,1,1,1)
224 end do
225
226 call neko_scratch_registry%request_field(visc_rho, tmp_indices(4), .false.)
227 call neko_scratch_registry%request_field(visc_m_x, tmp_indices(5), .false.)
228 call neko_scratch_registry%request_field(visc_m_y, tmp_indices(6), .false.)
229 call neko_scratch_registry%request_field(visc_m_z, tmp_indices(7), .false.)
230 call neko_scratch_registry%request_field(visc_e, tmp_indices(8), .false.)
231
232 ! Calculate artificial diffusion
233 call ax%compute(visc_rho%x, rho_field%x, coef, p%msh, p%Xh)
234 call ax%compute(visc_m_x%x, m_x%x, coef, p%msh, p%Xh)
235 call ax%compute(visc_m_y%x, m_y%x, coef, p%msh, p%Xh)
236 call ax%compute(visc_m_z%x, m_z%x, coef, p%msh, p%Xh)
237 call ax%compute(visc_e%x, e%x, coef, p%msh, p%Xh)
238
239 call gs%op(visc_rho, gs_op_add)
240 call gs%op(visc_m_x, gs_op_add)
241 call gs%op(visc_m_y, gs_op_add)
242 call gs%op(visc_m_z, gs_op_add)
243 call gs%op(visc_e, gs_op_add)
244
245 ! Move div to the rhs and apply the artificial viscosity
246 do concurrent(i = 1:n)
247 rhs_rho_field%x(i,1,1,1) = -rhs_rho_field%x(i,1,1,1) &
248 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_rho%x(i,1,1,1)
249 rhs_m_x%x(i,1,1,1) = -rhs_m_x%x(i,1,1,1) &
250 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_x%x(i,1,1,1)
251 rhs_m_y%x(i,1,1,1) = -rhs_m_y%x(i,1,1,1) &
252 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_y%x(i,1,1,1)
253 rhs_m_z%x(i,1,1,1) = -rhs_m_z%x(i,1,1,1) &
254 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_z%x(i,1,1,1)
255 rhs_e%x(i,1,1,1) = -rhs_e%x(i,1,1,1) &
256 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_e%x(i,1,1,1)
257 end do
258
259 call neko_scratch_registry%relinquish_field(tmp_indices)
260 end subroutine evaluate_rhs_sx
261
262end module euler_res_sx
ax_product
Defines a Matrix-vector product.
Definition ax.f90:34
coefs
Coefficients.
Definition coef.f90:34
euler_res_sx
Definition euler_res_sx.f90:1
euler_res_sx::advance_primitive_variables_sx
subroutine advance_primitive_variables_sx(rho_field, m_x, m_y, m_z, e, p, u, v, w, ax, coef, gs, h, c_avisc_low, rk_scheme, dt)
Definition euler_res_sx.f90:27
euler_res_sx::evaluate_rhs_sx
subroutine evaluate_rhs_sx(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_e, rho_field, m_x, m_y, m_z, e, p, u, v, w, ax, coef, gs, h, c_avisc_low)
Definition euler_res_sx.f90:154
euler_residual
Definition euler_res.f90:33
field_list
Definition field_list.f90:1
field
Defines a field.
Definition field.f90:34
gather_scatter
Gather-scatter.
Definition gather_scatter.f90:34
gs_ops
Defines Gather-scatter operations.
Definition gs_ops.f90:34
gs_ops::gs_op_add
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
math
Definition math.f90:60
math::cmult
subroutine, public cmult(a, c, n)
Multiplication by constant c .
Definition math.f90:411
math::subcol3
subroutine, public subcol3(a, b, c, n)
Returns .
Definition math.f90:881
math::add3
subroutine, public add3(a, b, c, n)
Vector addition .
Definition math.f90:739
math::add2
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:726
math::cfill
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition math.f90:487
math::invcol3
subroutine, public invcol3(a, b, c, n)
Invert a vector .
Definition math.f90:626
math::addcol3
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:958
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:854
math::copy
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:249
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:867
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:768
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
operators
Operators.
Definition operators.f90:34
operators::div
subroutine, public div(res, ux, uy, uz, coef)
Compute the divergence of a vector field.
Definition operators.f90:115
runge_kutta_time_scheme
Definition runge_kutta_scheme.f90:33
scratch_registry
Defines a registry for storing and requesting temporary objects This can be used when you have a func...
Definition scratch_registry.f90:37
scratch_registry::neko_scratch_registry
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Definition scratch_registry.f90:117
ax_product::ax_t
Base type for a matrix-vector product providing .
Definition ax.f90:43
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:56
euler_res_sx::euler_res_sx_t
Definition euler_res_sx.f90:18
euler_residual::euler_rhs_t
Abstract type to compute rhs.
Definition euler_res.f90:44
field::field_t
Definition field.f90:47
field_list::field_list_t
field_list_t, To be able to group fields together
Definition field_list.f90:15
gather_scatter::gs_t
Definition gather_scatter.f90:68
runge_kutta_time_scheme::runge_kutta_time_scheme_t
Definition runge_kutta_scheme.f90:44