Neko 1.99.1
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 :: temp_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, temp_indices(1))
49 call neko_scratch_registry%request_field(k_rho_2, temp_indices(2))
50 call neko_scratch_registry%request_field(k_rho_3, temp_indices(3))
51 call neko_scratch_registry%request_field(k_rho_4, temp_indices(4))
52 call neko_scratch_registry%request_field(k_m_x_1, temp_indices(5))
53 call neko_scratch_registry%request_field(k_m_x_2, temp_indices(6))
54 call neko_scratch_registry%request_field(k_m_x_3, temp_indices(7))
55 call neko_scratch_registry%request_field(k_m_x_4, temp_indices(8))
56 call neko_scratch_registry%request_field(k_m_y_1, temp_indices(9))
57 call neko_scratch_registry%request_field(k_m_y_2, temp_indices(10))
58 call neko_scratch_registry%request_field(k_m_y_3, temp_indices(11))
59 call neko_scratch_registry%request_field(k_m_y_4, temp_indices(12))
60 call neko_scratch_registry%request_field(k_m_z_1, temp_indices(13))
61 call neko_scratch_registry%request_field(k_m_z_2, temp_indices(14))
62 call neko_scratch_registry%request_field(k_m_z_3, temp_indices(15))
63 call neko_scratch_registry%request_field(k_m_z_4, temp_indices(16))
64 call neko_scratch_registry%request_field(k_e_1, temp_indices(17))
65 call neko_scratch_registry%request_field(k_e_2, temp_indices(18))
66 call neko_scratch_registry%request_field(k_e_3, temp_indices(19))
67 call neko_scratch_registry%request_field(k_e_4, temp_indices(20))
68 call neko_scratch_registry%request_field(temp_rho, temp_indices(21))
69 call neko_scratch_registry%request_field(temp_m_x, temp_indices(22))
70 call neko_scratch_registry%request_field(temp_m_y, temp_indices(23))
71 call neko_scratch_registry%request_field(temp_m_z, temp_indices(24))
72 call neko_scratch_registry%request_field(temp_e, temp_indices(25))
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(temp_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 :: temp, f_x, f_y, f_z, &
164 visc_rho, visc_m_x, visc_m_y, visc_m_z, visc_E
165 integer :: temp_indices(9)
166
167 n = coef%dof%size()
168 call neko_scratch_registry%request_field(temp, temp_indices(1))
169 call neko_scratch_registry%request_field(f_x, temp_indices(2))
170 call neko_scratch_registry%request_field(f_y, temp_indices(3))
171 call neko_scratch_registry%request_field(f_z, temp_indices(4))
172
174 call div(rhs_rho_field%x, m_x%x, m_y%x, m_z%x, coef)
175
177 ! m_x
178 do concurrent(i = 1:n)
179 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) &
180 + p%x(i,1,1,1)
181 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)
182 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)
183 end do
184 call div(rhs_m_x%x, f_x%x, f_y%x, f_z%x, coef)
185 ! m_y
186 do concurrent(i = 1:n)
187 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)
188 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) &
189 + p%x(i,1,1,1)
190 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)
191 end do
192 call div(rhs_m_y%x, f_x%x, f_y%x, f_z%x, coef)
193 ! m_z
194 do concurrent(i = 1:n)
195 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)
196 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)
197 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) &
198 + p%x(i,1,1,1)
199 end do
200 call div(rhs_m_z%x, f_x%x, f_y%x, f_z%x, coef)
201
203 do concurrent(i = 1:n)
204 f_x%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) &
205 * u%x(i,1,1,1)
206 f_y%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) &
207 * v%x(i,1,1,1)
208 f_z%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) &
209 * w%x(i,1,1,1)
210 end do
211 call div(rhs_e%x, f_x%x, f_y%x, f_z%x, coef)
212
213 call gs%op(rhs_rho_field, gs_op_add)
214 call gs%op(rhs_m_x, gs_op_add)
215 call gs%op(rhs_m_y, gs_op_add)
216 call gs%op(rhs_m_z, gs_op_add)
217 call gs%op(rhs_e, gs_op_add)
218
219 do concurrent(i = 1:rhs_e%dof%size())
220 rhs_rho_field%x(i,1,1,1) = rhs_rho_field%x(i,1,1,1) * coef%mult(i,1,1,1)
221 rhs_m_x%x(i,1,1,1) = rhs_m_x%x(i,1,1,1) * coef%mult(i,1,1,1)
222 rhs_m_y%x(i,1,1,1) = rhs_m_y%x(i,1,1,1) * coef%mult(i,1,1,1)
223 rhs_m_z%x(i,1,1,1) = rhs_m_z%x(i,1,1,1) * coef%mult(i,1,1,1)
224 rhs_e%x(i,1,1,1) = rhs_e%x(i,1,1,1) * coef%mult(i,1,1,1)
225 end do
226
227 call neko_scratch_registry%request_field(visc_rho, temp_indices(5))
228 call neko_scratch_registry%request_field(visc_m_x, temp_indices(6))
229 call neko_scratch_registry%request_field(visc_m_y, temp_indices(7))
230 call neko_scratch_registry%request_field(visc_m_z, temp_indices(8))
231 call neko_scratch_registry%request_field(visc_e, temp_indices(9))
232
233 ! Calculate artificial diffusion
234 call ax%compute(visc_rho%x, rho_field%x, coef, p%msh, p%Xh)
235 call ax%compute(visc_m_x%x, m_x%x, coef, p%msh, p%Xh)
236 call ax%compute(visc_m_y%x, m_y%x, coef, p%msh, p%Xh)
237 call ax%compute(visc_m_z%x, m_z%x, coef, p%msh, p%Xh)
238 call ax%compute(visc_e%x, e%x, coef, p%msh, p%Xh)
239
240 call gs%op(visc_rho, gs_op_add)
241 call gs%op(visc_m_x, gs_op_add)
242 call gs%op(visc_m_y, gs_op_add)
243 call gs%op(visc_m_z, gs_op_add)
244 call gs%op(visc_e, gs_op_add)
245
246 ! Move div to the rhs and apply the artificial viscosity
247 do concurrent(i = 1:n)
248 rhs_rho_field%x(i,1,1,1) = -rhs_rho_field%x(i,1,1,1) &
249 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_rho%x(i,1,1,1)
250 rhs_m_x%x(i,1,1,1) = -rhs_m_x%x(i,1,1,1) &
251 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_x%x(i,1,1,1)
252 rhs_m_y%x(i,1,1,1) = -rhs_m_y%x(i,1,1,1) &
253 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_y%x(i,1,1,1)
254 rhs_m_z%x(i,1,1,1) = -rhs_m_z%x(i,1,1,1) &
255 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_z%x(i,1,1,1)
256 rhs_e%x(i,1,1,1) = -rhs_e%x(i,1,1,1) &
257 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_e%x(i,1,1,1)
258 end do
259
260 call neko_scratch_registry%relinquish_field(temp_indices)
261 end subroutine evaluate_rhs_sx
262
263end 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:107
runge_kutta_time_scheme
Definition runge_kutta_scheme.f90:33
scratch_registry
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
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:84
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:55
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