76 use json_module,
only : json_file
88 real(kind=
rp) :: flow_rate
89 real(kind=
rp) :: dtlag = 0d0
90 real(kind=
rp) :: bdlag = 0d0
91 type(
field_t) :: u_vol, v_vol, w_vol, p_vol
92 real(kind=
rp) :: domain_length, base_flow
106 type(
dofmap_t),
intent(inout) :: dm_Xh
107 type(json_file),
intent(inout) :: params
108 logical average, found
110 real(kind=
rp) :: rate
115 call json_get(params,
'case.fluid.flow_rate_force.direction', direction)
116 call json_get(params,
'case.fluid.flow_rate_force.value', rate)
117 call json_get(params,
'case.fluid.flow_rate_force.use_averaged_flow',&
120 this%flow_dir = direction
121 this%avflow = average
122 this%flow_rate = rate
124 if (this%flow_dir .ne. 0)
then
125 call this%u_vol%init(dm_xh,
'u_vol')
126 call this%v_vol%init(dm_xh,
'v_vol')
127 call this%w_vol%init(dm_xh,
'w_vol')
128 call this%p_vol%init(dm_xh,
'p_vol')
138 call this%u_vol%free()
139 call this%v_vol%free()
140 call this%w_vol%free()
141 call this%p_vol%free()
143 call this%scratch%free()
151 ext_bdf, gs_Xh, c_Xh, rho, mu, bd, dt, &
152 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
153 Ax, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
155 type(
field_t),
intent(inout) :: u_res, v_res, w_res, p_res
156 type(
coef_t),
intent(inout) :: c_Xh
157 type(
gs_t),
intent(inout) :: gs_Xh
159 type(
bc_list_t),
intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
160 type(
bc_list_t),
intent(inout) :: bclst_vel_res
161 class(
ax_t),
intent(inout) :: Ax
162 class(
ksp_t),
intent(inout) :: ksp_prs, ksp_vel
163 class(
pc_t),
intent(inout) :: pc_prs, pc_vel
164 real(kind=
rp),
intent(inout) :: bd
165 real(kind=
rp),
intent(in) :: rho, mu, dt
166 integer,
intent(in) :: vel_max_iter, prs_max_iter
168 real(kind=
rp) :: xlmin, xlmax
169 real(kind=
rp) :: ylmin, ylmax
170 real(kind=
rp) :: zlmin, zlmax
172 type(
field_t),
pointer :: ta1, ta2, ta3
173 integer :: temp_indices(3)
175 call this%scratch%request_field(ta1, temp_indices(1))
176 call this%scratch%request_field(ta2, temp_indices(2))
177 call this%scratch%request_field(ta3, temp_indices(3))
180 associate(msh => c_xh%msh, p_vol => this%p_vol, &
181 u_vol => this%u_vol, v_vol => this%v_vol, w_vol => this%w_vol)
184 xlmin =
glmin(c_xh%dof%x, n)
185 xlmax =
glmax(c_xh%dof%x, n)
186 ylmin =
glmin(c_xh%dof%y, n)
187 ylmax =
glmax(c_xh%dof%y, n)
188 zlmin =
glmin(c_xh%dof%z, n)
189 zlmax =
glmax(c_xh%dof%z, n)
190 if (this%flow_dir.eq.1) this%domain_length = xlmax - xlmin
191 if (this%flow_dir.eq.2) this%domain_length = ylmax - ylmin
192 if (this%flow_dir.eq.3) this%domain_length = zlmax - zlmin
199 c_xh%h1(i,1,1,1) = 1.0_rp / rho
200 c_xh%h2(i,1,1,1) = 0.0_rp
207 if (this%flow_dir .eq. 1)
then
208 call cdtp(p_res%x, c_xh%h1, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
211 if (this%flow_dir .eq. 2)
then
212 call cdtp(p_res%x, c_xh%h1, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
215 if (this%flow_dir .eq. 3)
then
216 call cdtp(p_res%x, c_xh%h1, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
219 call gs_xh%op(p_res, gs_op_add)
222 ksp_result = ksp_prs%solve(ax, p_vol, p_res%x, n, &
223 c_xh, bclst_dp, gs_xh, prs_max_iter)
227 call opgrad(u_res%x, v_res%x, w_res%x, p_vol%x, c_xh)
236 call opchsign(u_res%x, v_res%x, w_res%x, msh%gdim, n)
237 call copy(ta1%x, c_xh%B, n)
238 call copy(ta2%x, c_xh%B, n)
239 call copy(ta3%x, c_xh%B, n)
242 ta1%x, ta2%x, ta3%x, n)
247 if (this%flow_dir .eq. 1)
then
249 else if (this%flow_dir .eq. 2)
then
251 else if (this%flow_dir .eq. 3)
then
255 if (this%flow_dir .eq. 1)
then
256 call add2(u_res%x, ta1%x, n)
257 else if (this%flow_dir .eq. 2)
then
258 call add2(v_res%x, ta2%x, n)
259 else if (this%flow_dir .eq. 3)
then
260 call add2(w_res%x, ta3%x, n)
269 c_xh%h1(i,1,1,1) = mu
270 c_xh%h2(i,1,1,1) = rho * (bd / dt)
275 call gs_xh%op(u_res, gs_op_add)
276 call gs_xh%op(v_res, gs_op_add)
277 call gs_xh%op(w_res, gs_op_add)
280 u_res%x, v_res%x, w_res%x, n)
283 ksp_result = ksp_vel%solve(ax, u_vol, u_res%x, n, &
284 c_xh, bclst_du, gs_xh, vel_max_iter)
285 ksp_result = ksp_vel%solve(ax, v_vol, v_res%x, n, &
286 c_xh, bclst_dv, gs_xh, vel_max_iter)
287 ksp_result = ksp_vel%solve(ax, w_vol, w_res%x, n, &
288 c_xh, bclst_dw, gs_xh, vel_max_iter)
291 if (this%flow_dir .eq. 1)
then
293 device_glsc2(u_vol%x_d, c_xh%B_d, n) / this%domain_length
296 if (this%flow_dir .eq. 2)
then
298 device_glsc2(v_vol%x_d, c_xh%B_d, n) / this%domain_length
301 if (this%flow_dir .eq. 3)
then
303 device_glsc2(w_vol%x_d, c_xh%B_d, n) / this%domain_length
306 if (this%flow_dir .eq. 1)
then
307 this%base_flow =
glsc2(u_vol%x, c_xh%B, n) / this%domain_length
310 if (this%flow_dir .eq. 2)
then
311 this%base_flow =
glsc2(v_vol%x, c_xh%B, n) / this%domain_length
314 if (this%flow_dir .eq. 3)
then
315 this%base_flow =
glsc2(w_vol%x, c_xh%B, n) / this%domain_length
320 call this%scratch%relinquish_field(temp_indices)
333 c_Xh, gs_Xh, ext_bdf, rho, mu, dt, &
334 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
335 Ax, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
338 type(field_t),
intent(inout) :: u, v, w, p
339 type(field_t),
intent(inout) :: u_res, v_res, w_res, p_res
340 type(coef_t),
intent(inout) :: c_Xh
341 type(gs_t),
intent(inout) :: gs_Xh
342 type(time_scheme_controller_t),
intent(inout) :: ext_bdf
343 real(kind=rp),
intent(in) :: rho, mu, dt
344 type(bc_list_t),
intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
345 type(bc_list_t),
intent(inout) :: bclst_vel_res
346 class(ax_t),
intent(inout) :: Ax
347 class(ksp_t),
intent(inout) :: ksp_prs, ksp_vel
348 class(pc_t),
intent(inout) :: pc_prs, pc_vel
349 integer,
intent(in) :: prs_max_iter, vel_max_iter
350 real(kind=rp) :: ifcomp, flow_rate, xsec
351 real(kind=rp) :: current_flow, delta_flow, base_flow, scale
353 type(field_t),
pointer :: ta1, ta2, ta3
354 integer :: temp_indices(3)
356 associate(u_vol => this%u_vol, v_vol => this%v_vol, &
357 w_vol => this%w_vol, p_vol => this%p_vol)
366 if (dt .ne. this%dtlag .or. ext_bdf%diffusion_coeffs(1) .ne. this%bdlag)
then
371 this%bdlag = ext_bdf%diffusion_coeffs(1)
373 call mpi_allreduce(mpi_in_place, ifcomp, 1, &
374 mpi_real_precision, mpi_sum, neko_comm, ierr)
376 if (ifcomp .gt. 0d0)
then
377 call this%compute(u_res, v_res, w_res, p_res, &
378 ext_bdf, gs_xh, c_xh, rho, mu, ext_bdf%diffusion_coeffs(1), dt, &
379 bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
380 ax, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
383 if (neko_bcknd_device .eq. 1)
then
384 if (this%flow_dir .eq. 1)
then
386 device_glsc2(u%x_d, c_xh%B_d, n) / this%domain_length
387 else if (this%flow_dir .eq. 2)
then
389 device_glsc2(v%x_d, c_xh%B_d, n) / this%domain_length
390 else if (this%flow_dir .eq. 3)
then
392 device_glsc2(w%x_d, c_xh%B_d, n) / this%domain_length
395 if (this%flow_dir .eq. 1)
then
396 current_flow = glsc2(u%x, c_xh%B, n) / this%domain_length
397 else if (this%flow_dir .eq. 2)
then
398 current_flow = glsc2(v%x, c_xh%B, n) / this%domain_length
399 else if (this%flow_dir .eq. 3)
then
400 current_flow = glsc2(w%x, c_xh%B, n) / this%domain_length
404 if (this%avflow)
then
405 xsec = c_xh%volume / this%domain_length
406 flow_rate = this%flow_rate*xsec
408 flow_rate = this%flow_rate
411 delta_flow = flow_rate - current_flow
412 scale = delta_flow / this%base_flow
414 if (neko_bcknd_device .eq. 1)
then
415 call device_add2s2(u%x_d, u_vol%x_d, scale, n)
416 call device_add2s2(v%x_d, v_vol%x_d, scale, n)
417 call device_add2s2(w%x_d, w_vol%x_d, scale, n)
418 call device_add2s2(p%x_d, p_vol%x_d, scale, n)
420 call add2s2(u%x, u_vol%x, scale, n)
421 call add2s2(v%x, v_vol%x, scale, n)
422 call add2s2(w%x, w_vol%x, scale, n)
423 call add2s2(p%x, p_vol%x, scale, n)
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)
subroutine, public device_rzero(a_d, n)
real(kind=rp) function, public device_glsc2(a_d, b_d, n)
subroutine, public device_copy(a_d, b_d, n)
subroutine, public device_cfill(a_d, c, n)
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 .
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 opchsign(a1, a2, a3, gdim, n)
for and .
integer, parameter neko_bcknd_hip
integer, parameter neko_bcknd_device
integer, parameter neko_bcknd_opencl
integer, parameter neko_bcknd_cuda
integer, parameter, public rp
Global precision used in computations.
subroutine, public cdtp(dtx, x, dr, ds, dt, coef)
Apply D^T to a scalar field, where D is the derivative matrix.
subroutine, public opgrad(ux, uy, uz, u, coef, es, ee)
Compute the gradient of a scalar field.
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 ...