152 ext_bdf, gs_Xh, c_Xh, rho, mu, bd, dt, &
153 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
154 Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
156 type(
field_t),
intent(inout) :: u_res, v_res, w_res, p_res
157 type(
coef_t),
intent(inout) :: c_Xh
158 type(
gs_t),
intent(inout) :: gs_Xh
160 type(
bc_list_t),
intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
161 type(
bc_list_t),
intent(inout) :: bclst_vel_res
162 class(
ax_t),
intent(inout) :: Ax_vel
163 class(
ax_t),
intent(inout) :: Ax_prs
164 class(
ksp_t),
intent(inout) :: ksp_prs, ksp_vel
165 class(
pc_t),
intent(inout) :: pc_prs, pc_vel
166 real(kind=
rp),
intent(inout) :: bd
167 real(kind=
rp),
intent(in) :: rho, mu, dt
168 integer,
intent(in) :: vel_max_iter, prs_max_iter
170 real(kind=
rp) :: xlmin, xlmax
171 real(kind=
rp) :: ylmin, ylmax
172 real(kind=
rp) :: zlmin, zlmax
174 type(
field_t),
pointer :: ta1, ta2, ta3
175 integer :: temp_indices(3)
177 call this%scratch%request_field(ta1, temp_indices(1))
178 call this%scratch%request_field(ta2, temp_indices(2))
179 call this%scratch%request_field(ta3, temp_indices(3))
182 associate(msh => c_xh%msh, p_vol => this%p_vol, &
183 u_vol => this%u_vol, v_vol => this%v_vol, w_vol => this%w_vol)
186 xlmin =
glmin(c_xh%dof%x, n)
187 xlmax =
glmax(c_xh%dof%x, n)
188 ylmin =
glmin(c_xh%dof%y, n)
189 ylmax =
glmax(c_xh%dof%y, n)
190 zlmin =
glmin(c_xh%dof%z, n)
191 zlmax =
glmax(c_xh%dof%z, n)
192 if (this%flow_dir .eq. 1)
then
193 this%domain_length = xlmax - xlmin
195 if (this%flow_dir .eq. 2)
then
196 this%domain_length = ylmax - ylmin
198 if (this%flow_dir .eq. 3)
then
199 this%domain_length = zlmax - zlmin
207 c_xh%h1(i,1,1,1) = 1.0_rp / rho
208 c_xh%h2(i,1,1,1) = 0.0_rp
215 if (this%flow_dir .eq. 1)
then
216 call cdtp(p_res%x, c_xh%h1, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
219 if (this%flow_dir .eq. 2)
then
220 call cdtp(p_res%x, c_xh%h1, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
223 if (this%flow_dir .eq. 3)
then
224 call cdtp(p_res%x, c_xh%h1, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
227 call gs_xh%op(p_res, gs_op_add)
230 ksp_results(1) = ksp_prs%solve(ax_prs, p_vol, p_res%x, n, &
231 c_xh, bclst_dp, gs_xh, prs_max_iter)
235 call opgrad(u_res%x, v_res%x, w_res%x, p_vol%x, c_xh)
244 call opchsign(u_res%x, v_res%x, w_res%x, msh%gdim, n)
245 call copy(ta1%x, c_xh%B, n)
246 call copy(ta2%x, c_xh%B, n)
247 call copy(ta3%x, c_xh%B, n)
250 ta1%x, ta2%x, ta3%x, n)
255 if (this%flow_dir .eq. 1)
then
257 else if (this%flow_dir .eq. 2)
then
259 else if (this%flow_dir .eq. 3)
then
263 if (this%flow_dir .eq. 1)
then
264 call add2(u_res%x, ta1%x, n)
265 else if (this%flow_dir .eq. 2)
then
266 call add2(v_res%x, ta2%x, n)
267 else if (this%flow_dir .eq. 3)
then
268 call add2(w_res%x, ta3%x, n)
277 c_xh%h1(i,1,1,1) = mu
278 c_xh%h2(i,1,1,1) = rho * (bd / dt)
283 call gs_xh%op(u_res, gs_op_add)
284 call gs_xh%op(v_res, gs_op_add)
285 call gs_xh%op(w_res, gs_op_add)
288 u_res%x, v_res%x, w_res%x, n)
291 ksp_results(2:4) = ksp_vel%solve_coupled(ax_vel, &
292 u_vol, v_vol, w_vol, &
293 u_res%x, v_res%x, w_res%x, &
295 bclst_du, bclst_dv, bclst_dw, &
299 if (this%flow_dir .eq. 1)
then
301 device_glsc2(u_vol%x_d, c_xh%B_d, n) / this%domain_length
304 if (this%flow_dir .eq. 2)
then
306 device_glsc2(v_vol%x_d, c_xh%B_d, n) / this%domain_length
309 if (this%flow_dir .eq. 3)
then
311 device_glsc2(w_vol%x_d, c_xh%B_d, n) / this%domain_length
314 if (this%flow_dir .eq. 1)
then
315 this%base_flow =
glsc2(u_vol%x, c_xh%B, n) / this%domain_length
318 if (this%flow_dir .eq. 2)
then
319 this%base_flow =
glsc2(v_vol%x, c_xh%B, n) / this%domain_length
322 if (this%flow_dir .eq. 3)
then
323 this%base_flow =
glsc2(w_vol%x, c_xh%B, n) / this%domain_length
328 call this%scratch%relinquish_field(temp_indices)
341 c_Xh, gs_Xh, ext_bdf, rho, mu, dt, &
342 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
343 Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
346 type(field_t),
intent(inout) :: u, v, w, p
347 type(field_t),
intent(inout) :: u_res, v_res, w_res, p_res
348 type(coef_t),
intent(inout) :: c_Xh
349 type(gs_t),
intent(inout) :: gs_Xh
350 type(time_scheme_controller_t),
intent(inout) :: ext_bdf
351 real(kind=rp),
intent(in) :: rho, mu, dt
352 type(bc_list_t),
intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
353 type(bc_list_t),
intent(inout) :: bclst_vel_res
354 class(ax_t),
intent(inout) :: Ax_vel
355 class(ax_t),
intent(inout) :: Ax_prs
356 class(ksp_t),
intent(inout) :: ksp_prs, ksp_vel
357 class(pc_t),
intent(inout) :: pc_prs, pc_vel
358 integer,
intent(in) :: prs_max_iter, vel_max_iter
359 real(kind=rp) :: ifcomp, flow_rate, xsec
360 real(kind=rp) :: current_flow, delta_flow, scale
361 integer :: n, ierr, i
362 type(field_t),
pointer :: ta1, ta2, ta3
364 associate(u_vol => this%u_vol, v_vol => this%v_vol, &
365 w_vol => this%w_vol, p_vol => this%p_vol)
374 if (dt .ne. this%dtlag .or. &
375 ext_bdf%diffusion_coeffs(1) .ne. this%bdlag)
then
380 this%bdlag = ext_bdf%diffusion_coeffs(1)
382 call mpi_allreduce(mpi_in_place, ifcomp, 1, &
383 mpi_real_precision, mpi_sum, neko_comm, ierr)
385 if (ifcomp .gt. 0d0)
then
386 call this%compute(u_res, v_res, w_res, p_res, &
387 ext_bdf, gs_xh, c_xh, rho, mu, ext_bdf%diffusion_coeffs(1), dt, &
388 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
389 ax_vel, ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, &
393 if (neko_bcknd_device .eq. 1)
then
394 if (this%flow_dir .eq. 1)
then
396 device_glsc2(u%x_d, c_xh%B_d, n) / this%domain_length
397 else if (this%flow_dir .eq. 2)
then
399 device_glsc2(v%x_d, c_xh%B_d, n) / this%domain_length
400 else if (this%flow_dir .eq. 3)
then
402 device_glsc2(w%x_d, c_xh%B_d, n) / this%domain_length
405 if (this%flow_dir .eq. 1)
then
406 current_flow = glsc2(u%x, c_xh%B, n) / this%domain_length
407 else if (this%flow_dir .eq. 2)
then
408 current_flow = glsc2(v%x, c_xh%B, n) / this%domain_length
409 else if (this%flow_dir .eq. 3)
then
410 current_flow = glsc2(w%x, c_xh%B, n) / this%domain_length
414 if (this%avflow)
then
415 xsec = c_xh%volume / this%domain_length
416 flow_rate = this%flow_rate*xsec
418 flow_rate = this%flow_rate
421 delta_flow = flow_rate - current_flow
422 scale = delta_flow / this%base_flow
424 if (neko_bcknd_device .eq. 1)
then
425 call device_add2s2(u%x_d, u_vol%x_d, scale, n)
426 call device_add2s2(v%x_d, v_vol%x_d, scale, n)
427 call device_add2s2(w%x_d, w_vol%x_d, scale, n)
428 call device_add2s2(p%x_d, p_vol%x_d, scale, n)
430 do concurrent(i = 1: n)
431 u%x(i,1,1,1) = u%x(i,1,1,1) + scale * u_vol%x(i,1,1,1)
432 v%x(i,1,1,1) = v%x(i,1,1,1) + scale * v_vol%x(i,1,1,1)
433 w%x(i,1,1,1) = w%x(i,1,1,1) + scale * w_vol%x(i,1,1,1)
434 p%x(i,1,1,1) = p%x(i,1,1,1) + scale * p_vol%x(i,1,1,1)
subroutine fluid_vol_flow(this, u, v, w, p, u_res, v_res, w_res, p_res, c_xh, gs_xh, ext_bdf, rho, mu, dt, bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, ax_vel, ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
Adjust flow volume.
subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, ext_bdf, gs_xh, c_xh, rho, mu, bd, dt, bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, ax_vel, ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
Compute flow adjustment.