Neko 0.9.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
fluid_volflow.f90
Go to the documentation of this file.
1! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
2!
3! The UChicago Argonne, LLC as Operator of Argonne National
4! Laboratory holds copyright in the Software. The copyright holder
5! reserves all rights except those expressly granted to licensees,
6! and U.S. Government license rights.
7!
8! Redistribution and use in source and binary forms, with or without
9! modification, are permitted provided that the following conditions
10! are met:
11!
12! 1. Redistributions of source code must retain the above copyright
13! notice, this list of conditions and the disclaimer below.
14!
15! 2. Redistributions in binary form must reproduce the above copyright
16! notice, this list of conditions and the disclaimer (as noted below)
17! in the documentation and/or other materials provided with the
18! distribution.
19!
20! 3. Neither the name of ANL nor the names of its contributors
21! may be used to endorse or promote products derived from this software
22! without specific prior written permission.
23!
24! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
27! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
28! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
29! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
31! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36!
37! Additional BSD Notice
38! ---------------------
39! 1. This notice is required to be provided under our contract with
40! the U.S. Department of Energy (DOE). This work was produced at
41! Argonne National Laboratory under Contract
42! No. DE-AC02-06CH11357 with the DOE.
43!
44! 2. Neither the United States Government nor UCHICAGO ARGONNE,
45! LLC nor any of their employees, makes any warranty,
46! express or implied, or assumes any liability or responsibility for the
47! accuracy, completeness, or usefulness of any information, apparatus,
48! product, or process disclosed, or represents that its use would not
49! infringe privately-owned rights.
50!
51! 3. Also, reference herein to any specific commercial products, process,
52! or services by trade name, trademark, manufacturer or otherwise does
53! not necessarily constitute or imply its endorsement, recommendation,
54! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
55! The views and opinions of authors expressed
56! herein do not necessarily state or reflect those of the United States
57! Government or UCHICAGO ARGONNE, LLC, and shall
58! not be used for advertising or product endorsement purposes.
59!
61 use operators, only : opgrad, cdtp
62 use num_types, only : rp
63 use mathops, only : opchsign
64 use krylov, only : ksp_t, ksp_monitor_t
65 use precon, only : pc_t
66 use dofmap, only : dofmap_t
67 use field, only : field_t
68 use coefs, only : coef_t
70 use math, only : copy, glsc2, glmin, glmax, add2, add2s2
71 use comm
76 use gather_scatter, only : gs_t, gs_op_add
77 use json_module, only : json_file
78 use json_utils, only: json_get
81 use ax_product, only : ax_t
82 implicit none
83 private
84
86 type, public :: fluid_volflow_t
87 integer :: flow_dir
88 logical :: avflow
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
95 type(scratch_registry_t) :: scratch
96 contains
97 procedure, pass(this) :: init => fluid_vol_flow_init
98 procedure, pass(this) :: free => fluid_vol_flow_free
99 procedure, pass(this) :: adjust => fluid_vol_flow
100 procedure, private, pass(this) :: compute => fluid_vol_flow_compute
101 end type fluid_volflow_t
102
103contains
104
105 subroutine fluid_vol_flow_init(this, dm_Xh, params)
106 class(fluid_volflow_t), intent(inout) :: this
107 type(dofmap_t), target, intent(inout) :: dm_Xh
108 type(json_file), intent(inout) :: params
109 logical average
110 integer :: direction
111 real(kind=rp) :: rate
112
113 call this%free()
114
115 !Initialize vol_flow (if there is a forced volume flow)
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',&
119 average)
120
121 this%flow_dir = direction
122 this%avflow = average
123 this%flow_rate = rate
124
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')
130 end if
131
132 this%scratch = scratch_registry_t(dm_xh, 3, 1)
133
134 end subroutine fluid_vol_flow_init
135
136 subroutine fluid_vol_flow_free(this)
137 class(fluid_volflow_t), intent(inout) :: this
138
139 call this%u_vol%free()
140 call this%v_vol%free()
141 call this%w_vol%free()
142 call this%p_vol%free()
143
144 call this%scratch%free()
145
146 end subroutine fluid_vol_flow_free
147
151 subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, &
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)
155 class(fluid_volflow_t), intent(inout) :: this
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
159 type(time_scheme_controller_t), intent(inout) :: ext_bdf
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
169 integer :: n, i
170 real(kind=rp) :: xlmin, xlmax
171 real(kind=rp) :: ylmin, ylmax
172 real(kind=rp) :: zlmin, zlmax
173 type(ksp_monitor_t) :: ksp_results(4)
174 type(field_t), pointer :: ta1, ta2, ta3
175 integer :: temp_indices(3)
176
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))
180
181
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)
184
185 n = c_xh%dof%size()
186 xlmin = glmin(c_xh%dof%x, n)
187 xlmax = glmax(c_xh%dof%x, n)
188 ylmin = glmin(c_xh%dof%y, n) ! for Y!
189 ylmax = glmax(c_xh%dof%y, n)
190 zlmin = glmin(c_xh%dof%z, n) ! for Z!
191 zlmax = glmax(c_xh%dof%z, n)
192 if (this%flow_dir .eq. 1) then
193 this%domain_length = xlmax - xlmin
194 end if
195 if (this%flow_dir .eq. 2) then
196 this%domain_length = ylmax - ylmin
197 end if
198 if (this%flow_dir .eq. 3) then
199 this%domain_length = zlmax - zlmin
200 end if
201
202 if (neko_bcknd_device .eq. 1) then
203 call device_cfill(c_xh%h1_d, 1.0_rp/rho, n)
204 call device_rzero(c_xh%h2_d, n)
205 else
206 do i = 1, n
207 c_xh%h1(i,1,1,1) = 1.0_rp / rho
208 c_xh%h2(i,1,1,1) = 0.0_rp
209 end do
210 end if
211 c_xh%ifh2 = .false.
212
213 ! Compute pressure
214
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)
217 end if
218
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)
221 end if
222
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)
225 end if
226
227 call gs_xh%op(p_res, gs_op_add)
228 call bc_list_apply_scalar(bclst_dp, p_res%x, n)
229 call pc_prs%update()
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)
232
233 ! Compute velocity
234
235 call opgrad(u_res%x, v_res%x, w_res%x, p_vol%x, c_xh)
236
237 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
238 (neko_bcknd_opencl .eq. 1)) then
239 call device_opchsign(u_res%x_d, v_res%x_d, w_res%x_d, msh%gdim, n)
240 call device_copy(ta1%x_d, c_xh%B_d, n)
241 call device_copy(ta2%x_d, c_xh%B_d, n)
242 call device_copy(ta3%x_d, c_xh%B_d, n)
243 else
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)
248 end if
249 call bc_list_apply_vector(bclst_vel_res,&
250 ta1%x, ta2%x, ta3%x, n)
251
252 ! add forcing
253
254 if (neko_bcknd_device .eq. 1) then
255 if (this%flow_dir .eq. 1) then
256 call device_add2(u_res%x_d, ta1%x_d, n)
257 else if (this%flow_dir .eq. 2) then
258 call device_add2(v_res%x_d, ta2%x_d, n)
259 else if (this%flow_dir .eq. 3) then
260 call device_add2(w_res%x_d, ta3%x_d, n)
261 end if
262 else
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)
269 end if
270 end if
271
272 if (neko_bcknd_device .eq. 1) then
273 call device_cfill(c_xh%h1_d, mu, n)
274 call device_cfill(c_xh%h2_d, rho * (bd / dt), n)
275 else
276 do i = 1, n
277 c_xh%h1(i,1,1,1) = mu
278 c_xh%h2(i,1,1,1) = rho * (bd / dt)
279 end do
280 end if
281 c_xh%ifh2 = .true.
282
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)
286
287 call bc_list_apply_vector(bclst_vel_res,&
288 u_res%x, v_res%x, w_res%x, n)
289 call pc_vel%update()
290
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, &
294 n, c_xh, &
295 bclst_du, bclst_dv, bclst_dw, &
296 gs_xh, vel_max_iter)
297
298 if (neko_bcknd_device .eq. 1) then
299 if (this%flow_dir .eq. 1) then
300 this%base_flow = &
301 device_glsc2(u_vol%x_d, c_xh%B_d, n) / this%domain_length
302 end if
303
304 if (this%flow_dir .eq. 2) then
305 this%base_flow = &
306 device_glsc2(v_vol%x_d, c_xh%B_d, n) / this%domain_length
307 end if
308
309 if (this%flow_dir .eq. 3) then
310 this%base_flow = &
311 device_glsc2(w_vol%x_d, c_xh%B_d, n) / this%domain_length
312 end if
313 else
314 if (this%flow_dir .eq. 1) then
315 this%base_flow = glsc2(u_vol%x, c_xh%B, n) / this%domain_length
316 end if
317
318 if (this%flow_dir .eq. 2) then
319 this%base_flow = glsc2(v_vol%x, c_xh%B, n) / this%domain_length
320 end if
321
322 if (this%flow_dir .eq. 3) then
323 this%base_flow = glsc2(w_vol%x, c_xh%B, n) / this%domain_length
324 end if
325 end if
326 end associate
327
328 call this%scratch%relinquish_field(temp_indices)
329 end subroutine fluid_vol_flow_compute
330
340 subroutine fluid_vol_flow(this, u, v, w, p, u_res, v_res, w_res, p_res, &
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)
344
345 class(fluid_volflow_t), intent(inout) :: this
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
363
364 associate(u_vol => this%u_vol, v_vol => this%v_vol, &
365 w_vol => this%w_vol, p_vol => this%p_vol)
366
367 n = c_xh%dof%size()
368
369 ! If either dt or the backwards difference coefficient change,
370 ! then recompute base flow solution corresponding to unit forcing:
371
372 ifcomp = 0.0_rp
373
374 if (dt .ne. this%dtlag .or. &
375 ext_bdf%diffusion_coeffs(1) .ne. this%bdlag) then
376 ifcomp = 1.0_rp
377 end if
378
379 this%dtlag = dt
380 this%bdlag = ext_bdf%diffusion_coeffs(1)
381
382 call mpi_allreduce(mpi_in_place, ifcomp, 1, &
383 mpi_real_precision, mpi_sum, neko_comm, ierr)
384
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, &
390 vel_max_iter)
391 end if
392
393 if (neko_bcknd_device .eq. 1) then
394 if (this%flow_dir .eq. 1) then
395 current_flow = &
396 device_glsc2(u%x_d, c_xh%B_d, n) / this%domain_length ! for X
397 else if (this%flow_dir .eq. 2) then
398 current_flow = &
399 device_glsc2(v%x_d, c_xh%B_d, n) / this%domain_length ! for Y
400 else if (this%flow_dir .eq. 3) then
401 current_flow = &
402 device_glsc2(w%x_d, c_xh%B_d, n) / this%domain_length ! for Z
403 end if
404 else
405 if (this%flow_dir .eq. 1) then
406 current_flow = glsc2(u%x, c_xh%B, n) / this%domain_length ! for X
407 else if (this%flow_dir .eq. 2) then
408 current_flow = glsc2(v%x, c_xh%B, n) / this%domain_length ! for Y
409 else if (this%flow_dir .eq. 3) then
410 current_flow = glsc2(w%x, c_xh%B, n) / this%domain_length ! for Z
411 end if
412 end if
413
414 if (this%avflow) then
415 xsec = c_xh%volume / this%domain_length
416 flow_rate = this%flow_rate*xsec
417 else
418 flow_rate = this%flow_rate
419 end if
420
421 delta_flow = flow_rate - current_flow
422 scale = delta_flow / this%base_flow
423
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)
429 else
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)
435 end do
436 end if
437 end associate
438
439 end subroutine fluid_vol_flow
440
441
442end module fluid_volflow
Retrieves a parameter by name or throws an error.
Defines a Matrix-vector product.
Definition ax.f90:34
Defines a boundary condition.
Definition bc.f90:34
subroutine, public bc_list_apply_vector(bclst, x, y, z, n, t, tstep)
Apply a list of boundary conditions to a vector field.
Definition bc.f90:587
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
Definition bc.f90:530
Coefficients.
Definition coef.f90:34
Definition comm.F90:1
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.
Definition dofmap.f90:35
Defines a field.
Definition field.f90:34
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_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)
Gather-scatter.
Utilities for retrieving parameters from the case files.
Implements the base abstract type for Krylov solvers plus helper types.
Definition krylov.f90:34
Definition math.f90:60
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
Definition math.f90:876
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:587
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition math.f90:377
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:239
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition math.f90:673
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition math.f90:407
Collection of vector field operations operating on and . Note that in general the indices and ....
Definition mathops.f90:65
subroutine, public opchsign(a1, a2, a3, gdim, n)
for and .
Definition mathops.f90:76
Build configurations.
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.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
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.
Krylov preconditioner.
Definition precon.f90:34
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
type(scratch_registry_t) function init(dof, size, expansion_size)
Constructor, optionally taking initial registry and expansion size as argument.
Compound scheme for the advection and diffusion operators in a transport equation.
Base type for a matrix-vector product providing .
Definition ax.f90:43
A list of boundary conditions.
Definition bc.f90:104
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Type for storing initial and final residuals in a Krylov solver.
Definition krylov.f90:56
Base abstract type for a canonical Krylov method, solving .
Definition krylov.f90:66
Defines a canonical Krylov preconditioner.
Definition precon.f90:40
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...