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, 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
163 class(
ksp_t),
intent(inout) :: ksp_prs, ksp_vel
164 class(
pc_t),
intent(inout) :: pc_prs, pc_vel
165 real(kind=
rp),
intent(inout) :: bd
166 real(kind=
rp),
intent(in) :: rho, mu, dt
167 integer,
intent(in) :: vel_max_iter, prs_max_iter
169 real(kind=
rp) :: xlmin, xlmax
170 real(kind=
rp) :: ylmin, ylmax
171 real(kind=
rp) :: zlmin, zlmax
173 type(
field_t),
pointer :: ta1, ta2, ta3
174 integer :: temp_indices(3)
176 call this%scratch%request_field(ta1, temp_indices(1))
177 call this%scratch%request_field(ta2, temp_indices(2))
178 call this%scratch%request_field(ta3, temp_indices(3))
181 associate(msh => c_xh%msh, p_vol => this%p_vol, &
182 u_vol => this%u_vol, v_vol => this%v_vol, w_vol => this%w_vol)
185 xlmin =
glmin(c_xh%dof%x, n)
186 xlmax =
glmax(c_xh%dof%x, n)
187 ylmin =
glmin(c_xh%dof%y, n)
188 ylmax =
glmax(c_xh%dof%y, n)
189 zlmin =
glmin(c_xh%dof%z, n)
190 zlmax =
glmax(c_xh%dof%z, n)
191 if (this%flow_dir .eq. 1)
then
192 this%domain_length = xlmax - xlmin
194 if (this%flow_dir .eq. 2)
then
195 this%domain_length = ylmax - ylmin
197 if (this%flow_dir .eq. 3)
then
198 this%domain_length = zlmax - zlmin
206 c_xh%h1(i,1,1,1) = 1.0_rp / rho
207 c_xh%h2(i,1,1,1) = 0.0_rp
214 if (this%flow_dir .eq. 1)
then
215 call cdtp(p_res%x, c_xh%h1, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
218 if (this%flow_dir .eq. 2)
then
219 call cdtp(p_res%x, c_xh%h1, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
222 if (this%flow_dir .eq. 3)
then
223 call cdtp(p_res%x, c_xh%h1, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
226 call gs_xh%op(p_res, gs_op_add)
229 ksp_result = ksp_prs%solve(ax, p_vol, p_res%x, n, &
230 c_xh, bclst_dp, gs_xh, prs_max_iter)
234 call opgrad(u_res%x, v_res%x, w_res%x, p_vol%x, c_xh)
236 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
237 (neko_bcknd_opencl .eq. 1))
then
243 call opchsign(u_res%x, v_res%x, w_res%x, msh%gdim, n)
244 call copy(ta1%x, c_xh%B, n)
245 call copy(ta2%x, c_xh%B, n)
246 call copy(ta3%x, c_xh%B, n)
249 ta1%x, ta2%x, ta3%x, n)
254 if (this%flow_dir .eq. 1)
then
256 else if (this%flow_dir .eq. 2)
then
258 else if (this%flow_dir .eq. 3)
then
262 if (this%flow_dir .eq. 1)
then
263 call add2(u_res%x, ta1%x, n)
264 else if (this%flow_dir .eq. 2)
then
265 call add2(v_res%x, ta2%x, n)
266 else if (this%flow_dir .eq. 3)
then
267 call add2(w_res%x, ta3%x, n)
276 c_xh%h1(i,1,1,1) = mu
277 c_xh%h2(i,1,1,1) = rho * (bd / dt)
282 call gs_xh%op(u_res, gs_op_add)
283 call gs_xh%op(v_res, gs_op_add)
284 call gs_xh%op(w_res, gs_op_add)
287 u_res%x, v_res%x, w_res%x, n)
290 ksp_result = ksp_vel%solve(ax, u_vol, u_res%x, n, &
291 c_xh, bclst_du, gs_xh, vel_max_iter)
292 ksp_result = ksp_vel%solve(ax, v_vol, v_res%x, n, &
293 c_xh, bclst_dv, gs_xh, vel_max_iter)
294 ksp_result = ksp_vel%solve(ax, w_vol, w_res%x, n, &
295 c_xh, bclst_dw, gs_xh, vel_max_iter)
298 if (this%flow_dir .eq. 1)
then
300 device_glsc2(u_vol%x_d, c_xh%B_d, n) / this%domain_length
303 if (this%flow_dir .eq. 2)
then
305 device_glsc2(v_vol%x_d, c_xh%B_d, n) / this%domain_length
308 if (this%flow_dir .eq. 3)
then
310 device_glsc2(w_vol%x_d, c_xh%B_d, n) / this%domain_length
313 if (this%flow_dir .eq. 1)
then
314 this%base_flow =
glsc2(u_vol%x, c_xh%B, n) / this%domain_length
317 if (this%flow_dir .eq. 2)
then
318 this%base_flow =
glsc2(v_vol%x, c_xh%B, n) / this%domain_length
321 if (this%flow_dir .eq. 3)
then
322 this%base_flow =
glsc2(w_vol%x, c_xh%B, n) / this%domain_length
327 call this%scratch%relinquish_field(temp_indices)
340 c_Xh, gs_Xh, ext_bdf, rho, mu, dt, &
341 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
342 Ax, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
345 type(field_t),
intent(inout) :: u, v, w, p
346 type(field_t),
intent(inout) :: u_res, v_res, w_res, p_res
347 type(coef_t),
intent(inout) :: c_Xh
348 type(gs_t),
intent(inout) :: gs_Xh
349 type(time_scheme_controller_t),
intent(inout) :: ext_bdf
350 real(kind=rp),
intent(in) :: rho, mu, dt
351 type(bc_list_t),
intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
352 type(bc_list_t),
intent(inout) :: bclst_vel_res
353 class(ax_t),
intent(inout) :: Ax
354 class(ksp_t),
intent(inout) :: ksp_prs, ksp_vel
355 class(pc_t),
intent(inout) :: pc_prs, pc_vel
356 integer,
intent(in) :: prs_max_iter, vel_max_iter
357 real(kind=rp) :: ifcomp, flow_rate, xsec
358 real(kind=rp) :: current_flow, delta_flow, scale
359 integer :: n, ierr, i
360 type(field_t),
pointer :: ta1, ta2, ta3
362 associate(u_vol => this%u_vol, v_vol => this%v_vol, &
363 w_vol => this%w_vol, p_vol => this%p_vol)
372 if (dt .ne. this%dtlag .or. &
373 ext_bdf%diffusion_coeffs(1) .ne. this%bdlag)
then
378 this%bdlag = ext_bdf%diffusion_coeffs(1)
380 call mpi_allreduce(mpi_in_place, ifcomp, 1, &
381 mpi_real_precision, mpi_sum, neko_comm, ierr)
383 if (ifcomp .gt. 0d0)
then
384 call this%compute(u_res, v_res, w_res, p_res, &
385 ext_bdf, gs_xh, c_xh, rho, mu, ext_bdf%diffusion_coeffs(1), dt, &
386 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
387 ax, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
390 if (neko_bcknd_device .eq. 1)
then
391 if (this%flow_dir .eq. 1)
then
393 device_glsc2(u%x_d, c_xh%B_d, n) / this%domain_length
394 else if (this%flow_dir .eq. 2)
then
396 device_glsc2(v%x_d, c_xh%B_d, n) / this%domain_length
397 else if (this%flow_dir .eq. 3)
then
399 device_glsc2(w%x_d, c_xh%B_d, n) / this%domain_length
402 if (this%flow_dir .eq. 1)
then
403 current_flow = glsc2(u%x, c_xh%B, n) / this%domain_length
404 else if (this%flow_dir .eq. 2)
then
405 current_flow = glsc2(v%x, c_xh%B, n) / this%domain_length
406 else if (this%flow_dir .eq. 3)
then
407 current_flow = glsc2(w%x, c_xh%B, n) / this%domain_length
411 if (this%avflow)
then
412 xsec = c_xh%volume / this%domain_length
413 flow_rate = this%flow_rate*xsec
415 flow_rate = this%flow_rate
418 delta_flow = flow_rate - current_flow
419 scale = delta_flow / this%base_flow
421 if (neko_bcknd_device .eq. 1)
then
422 call device_add2s2(u%x_d, u_vol%x_d, scale, n)
423 call device_add2s2(v%x_d, v_vol%x_d, scale, n)
424 call device_add2s2(w%x_d, w_vol%x_d, scale, n)
425 call device_add2s2(p%x_d, p_vol%x_d, scale, n)
427 do concurrent(i = 1: n)
428 u%x(i,1,1,1) = u%x(i,1,1,1) + scale * u_vol%x(i,1,1,1)
429 v%x(i,1,1,1) = v%x(i,1,1,1) + scale * v_vol%x(i,1,1,1)
430 w%x(i,1,1,1) = w%x(i,1,1,1) + scale * w_vol%x(i,1,1,1)
431 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, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
Compute flow adjustment.
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, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
Adjust flow volume.
subroutine fluid_vol_flow_free(this)
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 ...