77 use json_module,
only : json_file
89 real(kind=
rp) :: flow_rate
90 real(kind=
rp) :: dtlag = 0d0
91 real(kind=
rp) :: bdlag = 0d0
92 type(
field_t) :: u_vol, v_vol, w_vol, p_vol
93 real(kind=
rp) :: domain_length, base_flow
107 type(
dofmap_t),
target,
intent(inout) :: dm_Xh
108 type(json_file),
intent(inout) :: params
111 real(kind=
rp) :: rate
116 call json_get(params,
'case.fluid.flow_rate_force.direction', direction)
117 call json_get(params,
'case.fluid.flow_rate_force.value', rate)
118 call json_get(params,
'case.fluid.flow_rate_force.use_averaged_flow',&
121 this%flow_dir = direction
122 this%avflow = average
123 this%flow_rate = rate
125 if (this%flow_dir .ne. 0)
then
126 call this%u_vol%init(dm_xh,
'u_vol')
127 call this%v_vol%init(dm_xh,
'v_vol')
128 call this%w_vol%init(dm_xh,
'w_vol')
129 call this%p_vol%init(dm_xh,
'p_vol')
139 call this%u_vol%free()
140 call this%v_vol%free()
141 call this%w_vol%free()
142 call this%p_vol%free()
144 call this%scratch%free()
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)
237 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
238 (neko_bcknd_opencl .eq. 1))
then
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)
Retrieves a parameter by name or throws an error.
Defines a Matrix-vector product.
Defines a boundary condition.
subroutine, public bc_list_apply_vector(bclst, x, y, z, n, t, tstep)
Apply a list of boundary conditions to a vector field.
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
subroutine, public device_add2(a_d, b_d, n)
Vector addition .
subroutine, public device_rzero(a_d, n)
Zero a real vector.
subroutine, public device_add2s2(a_d, b_d, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
real(kind=rp) function, public device_glsc2(a_d, b_d, n)
Weighted inner product .
subroutine, public device_copy(a_d, b_d, n)
Copy a vector .
subroutine, public device_cfill(a_d, c, n)
Set all elements to a constant c .
subroutine, public device_opchsign(a1_d, a2_d, a3_d, gdim, n)
Defines a mapping of the degrees of freedom.
subroutine fluid_vol_flow_init(this, dm_Xh, params)
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.
subroutine fluid_vol_flow_free(this)
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.
Utilities for retrieving parameters from the case files.
Implements the base abstract type for Krylov solvers plus helper types.
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
subroutine, public add2(a, b, n)
Vector addition .
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Collection of vector field operations operating on and . Note that in general the indices and ....
subroutine, public opchsign(a1, a2, a3, gdim, n)
for and .
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
subroutine, public opgrad(ux, uy, uz, u, coef, es, ee)
Compute the weak gradient of a scalar field, i.e. the gradient multiplied by the mass matrix.
subroutine, public cdtp(dtx, x, dr, ds, dt, coef, es, ee)
Apply D^T to a scalar field, where D is the derivative matrix.
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
Compound scheme for the advection and diffusion operators in a transport equation.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Type for storing initial and final residuals in a Krylov solver.
Base abstract type for a canonical Krylov method, solving .
Defines a canonical Krylov preconditioner.
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...