Neko  0.9.0
A portable framework for high-order spectral element flow simulations
fluid_pnpn.f90
Go to the documentation of this file.
1 ! Copyright (c) 2022-2024, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
34 module fluid_pnpn
35  use num_types, only : rp
36  use krylov, only : ksp_monitor_t
38  pnpn_prs_res_factory, pnpn_vel_res_factory, &
39  pnpn_prs_res_stress_factory, pnpn_vel_res_stress_factory
41  rhs_maker_oifs_t, rhs_maker_sumab_fctry, rhs_maker_bdf_fctry, &
42  rhs_maker_ext_fctry, rhs_maker_oifs_fctry
43  use fluid_volflow, only : fluid_volflow_t
44  use fluid_scheme, only : fluid_scheme_t
46  use fluid_aux, only : fluid_step_info
48  use projection, only : projection_t
50  use advection, only : advection_t, advection_factory
53  use json_module, only : json_file
54  use ax_product, only : ax_t, ax_helm_factory
55  use field, only : field_t
56  use dirichlet, only : dirichlet_t
57  use facet_normal, only : facet_normal_t
58  use non_normal, only : non_normal_t
59  use mesh, only : mesh_t
60  use user_intf, only : user_t
62  use gs_ops, only : gs_op_add
63  use neko_config, only : neko_bcknd_device
64  use math, only : col2, glsum
65  use mathops, only : opadd2cm, opcolv
68  use utils, only : neko_error
69  use field_math, only : field_add2, field_copy
70  implicit none
71  private
72 
73 
74  type, public, extends(fluid_scheme_t) :: fluid_pnpn_t
75  type(field_t) :: p_res, u_res, v_res, w_res
76 
77  type(field_t) :: dp, du, dv, dw
78 
79  ! Coupled Helmholz operator for velocity
80  class(ax_t), allocatable :: ax_vel
81  ! Helmholz operator for pressure
82  class(ax_t), allocatable :: ax_prs
83 
84  type(projection_t) :: proj_prs
85  type(projection_t) :: proj_u
86  type(projection_t) :: proj_v
87  type(projection_t) :: proj_w
88 
89  type(facet_normal_t) :: bc_prs_surface
90  type(facet_normal_t) :: bc_sym_surface
91  type(dirichlet_t) :: bc_vel_res
92  type(dirichlet_t) :: bc_field_dirichlet_p
93  type(dirichlet_t) :: bc_field_dirichlet_u
94  type(dirichlet_t) :: bc_field_dirichlet_v
95  type(dirichlet_t) :: bc_field_dirichlet_w
96  type(non_normal_t) :: bc_vel_res_non_normal
97  type(bc_list_t) :: bclst_vel_res
98  type(bc_list_t) :: bclst_du
99  type(bc_list_t) :: bclst_dv
100  type(bc_list_t) :: bclst_dw
101  type(bc_list_t) :: bclst_dp
102 
103  class(advection_t), allocatable :: adv
104 
105  ! Time interpolation scheme
106  logical :: oifs
107 
108  ! Time variables
109  type(field_t) :: abx1, aby1, abz1
110  type(field_t) :: abx2, aby2, abz2
111  ! Advection terms for the oifs method
112  type(field_t) :: advx, advy, advz
113 
115  class(pnpn_prs_res_t), allocatable :: prs_res
116 
118  class(pnpn_vel_res_t), allocatable :: vel_res
119 
121  class(rhs_maker_sumab_t), allocatable :: sumab
122 
124  class(rhs_maker_ext_t), allocatable :: makeabf
125 
127  class(rhs_maker_bdf_t), allocatable :: makebdf
128 
130  class(rhs_maker_oifs_t), allocatable :: makeoifs
131 
133  type(fluid_volflow_t) :: vol_flow
134 
135  contains
136  procedure, pass(this) :: init => fluid_pnpn_init
137  procedure, pass(this) :: free => fluid_pnpn_free
138  procedure, pass(this) :: step => fluid_pnpn_step
139  procedure, pass(this) :: restart => fluid_pnpn_restart
140  end type fluid_pnpn_t
141 
142 contains
143 
144  subroutine fluid_pnpn_init(this, msh, lx, params, user, time_scheme)
145  class(fluid_pnpn_t), target, intent(inout) :: this
146  type(mesh_t), target, intent(inout) :: msh
147  integer, intent(inout) :: lx
148  type(json_file), target, intent(inout) :: params
149  type(user_t), target, intent(in) :: user
150  type(time_scheme_controller_t), target, intent(in) :: time_scheme
151  character(len=15), parameter :: scheme = 'Modular (Pn/Pn)'
152 
153  call this%free()
154 
155  ! Initialize base class
156  call this%scheme_init(msh, lx, params, .true., .true., scheme, user)
157 
158  if (this%variable_material_properties .eqv. .true.) then
159  ! Setup backend dependent Ax routines
160  call ax_helm_factory(this%Ax_vel, full_formulation = .true.)
161 
162  ! Setup backend dependent prs residual routines
163  call pnpn_prs_res_stress_factory(this%prs_res)
164 
165  ! Setup backend dependent vel residual routines
166  call pnpn_vel_res_stress_factory(this%vel_res)
167  else
168  ! Setup backend dependent Ax routines
169  call ax_helm_factory(this%Ax_vel, full_formulation = .false.)
170 
171  ! Setup backend dependent prs residual routines
172  call pnpn_prs_res_factory(this%prs_res)
173 
174  ! Setup backend dependent vel residual routines
175  call pnpn_vel_res_factory(this%vel_res)
176  end if
177 
178  ! Setup Ax for the pressure
179  call ax_helm_factory(this%Ax_prs, full_formulation = .false.)
180 
181 
182  ! Setup backend dependent summation of AB/BDF
183  call rhs_maker_sumab_fctry(this%sumab)
184 
185  ! Setup backend dependent summation of extrapolation scheme
186  call rhs_maker_ext_fctry(this%makeabf)
187 
188  ! Setup backend depenent contributions to F from lagged BD terms
189  call rhs_maker_bdf_fctry(this%makebdf)
190 
191  ! Setup backend dependent summations of the OIFS method
192  call rhs_maker_oifs_fctry(this%makeoifs)
193 
194  ! Initialize variables specific to this plan
195  associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
196  dm_xh => this%dm_Xh, nelv => this%msh%nelv)
197 
198  call this%p_res%init(dm_xh, "p_res")
199  call this%u_res%init(dm_xh, "u_res")
200  call this%v_res%init(dm_xh, "v_res")
201  call this%w_res%init(dm_xh, "w_res")
202  call this%abx1%init(dm_xh, "abx1")
203  call this%aby1%init(dm_xh, "aby1")
204  call this%abz1%init(dm_xh, "abz1")
205  call this%abx2%init(dm_xh, "abx2")
206  call this%aby2%init(dm_xh, "aby2")
207  call this%abz2%init(dm_xh, "abz2")
208  call this%advx%init(dm_xh, "advx")
209  call this%advy%init(dm_xh, "advy")
210  call this%advz%init(dm_xh, "advz")
211  this%abx1 = 0.0_rp
212  this%aby1 = 0.0_rp
213  this%abz1 = 0.0_rp
214  this%abx2 = 0.0_rp
215  this%aby2 = 0.0_rp
216  this%abz2 = 0.0_rp
217  this%advx = 0.0_rp
218  this%advy = 0.0_rp
219  this%advz = 0.0_rp
220 
221  call this%du%init(dm_xh, 'du')
222  call this%dv%init(dm_xh, 'dv')
223  call this%dw%init(dm_xh, 'dw')
224  call this%dp%init(dm_xh, 'dp')
225 
226  end associate
227 
228  ! Initialize velocity surface terms in pressure rhs
229  call this%bc_prs_surface%init_base(this%c_Xh)
230  call this%bc_prs_surface%mark_zone(msh%inlet)
231  call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,&
232  'v', this%bc_labels)
233  ! This impacts the rhs of the pressure,
234  ! need to check what is correct to add here
235  call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,&
236  'd_vel_u', this%bc_labels)
237  call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,&
238  'd_vel_v', this%bc_labels)
239  call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,&
240  'd_vel_w', this%bc_labels)
241  call this%bc_prs_surface%finalize()
242  ! Initialize symmetry surface terms in pressure rhs
243  call this%bc_sym_surface%init_base(this%c_Xh)
244  call this%bc_sym_surface%mark_zone(msh%sympln)
245  call this%bc_sym_surface%mark_zones_from_list(msh%labeled_zones,&
246  'sym', this%bc_labels)
247  ! Same here, should du, dv, dw be marked here?
248  call this%bc_sym_surface%finalize()
249  ! Initialize dirichlet bcs for velocity residual
250  call this%bc_vel_res_non_normal%init_base(this%c_Xh)
251  call this%bc_vel_res_non_normal%mark_zone(msh%outlet_normal)
252  call this%bc_vel_res_non_normal%mark_zones_from_list(msh%labeled_zones,&
253  'on', this%bc_labels)
254  call this%bc_vel_res_non_normal%mark_zones_from_list(msh%labeled_zones,&
255  'on+dong', &
256  this%bc_labels)
257  call this%bc_vel_res_non_normal%finalize()
258  call this%bc_vel_res_non_normal%init(this%c_Xh)
259 
260  call this%bc_field_dirichlet_p%init_base(this%c_Xh)
261  call this%bc_field_dirichlet_p%mark_zones_from_list(msh%labeled_zones, &
262  'on+dong', this%bc_labels)
263  call this%bc_field_dirichlet_p%mark_zones_from_list(msh%labeled_zones, &
264  'o+dong', this%bc_labels)
265  call this%bc_field_dirichlet_p%mark_zones_from_list(msh%labeled_zones, &
266  'd_pres', this%bc_labels)
267  call this%bc_field_dirichlet_p%finalize()
268  call this%bc_field_dirichlet_p%set_g(0.0_rp)
269  call bc_list_init(this%bclst_dp)
270  call bc_list_add(this%bclst_dp, this%bc_field_dirichlet_p)
271  !Add 0 prs bcs
272  call bc_list_add(this%bclst_dp, this%bc_prs)
273 
274  call this%bc_field_dirichlet_u%init_base(this%c_Xh)
275  call this%bc_field_dirichlet_u%mark_zones_from_list( &
276  msh%labeled_zones, 'd_vel_u', this%bc_labels)
277  call this%bc_field_dirichlet_u%finalize()
278  call this%bc_field_dirichlet_u%set_g(0.0_rp)
279 
280  call this%bc_field_dirichlet_v%init_base(this%c_Xh)
281  call this%bc_field_dirichlet_v%mark_zones_from_list(msh%labeled_zones, &
282  'd_vel_v', &
283  this%bc_labels)
284  call this%bc_field_dirichlet_v%finalize()
285  call this%bc_field_dirichlet_v%set_g(0.0_rp)
286 
287  call this%bc_field_dirichlet_w%init_base(this%c_Xh)
288  call this%bc_field_dirichlet_w%mark_zones_from_list(msh%labeled_zones, &
289  'd_vel_w', &
290  this%bc_labels)
291  call this%bc_field_dirichlet_w%finalize()
292  call this%bc_field_dirichlet_w%set_g(0.0_rp)
293 
294  call this%bc_vel_res%init_base(this%c_Xh)
295  call this%bc_vel_res%mark_zone(msh%inlet)
296  call this%bc_vel_res%mark_zone(msh%wall)
297  call this%bc_vel_res%mark_zones_from_list(msh%labeled_zones, &
298  'v', this%bc_labels)
299  call this%bc_vel_res%mark_zones_from_list(msh%labeled_zones, &
300  'w', this%bc_labels)
301  call this%bc_vel_res%finalize()
302  call this%bc_vel_res%set_g(0.0_rp)
303  call bc_list_init(this%bclst_vel_res)
304  call bc_list_add(this%bclst_vel_res, this%bc_vel_res)
305  call bc_list_add(this%bclst_vel_res, this%bc_vel_res_non_normal)
306  call bc_list_add(this%bclst_vel_res, this%bc_sym)
307  call bc_list_add(this%bclst_vel_res, this%bc_sh%symmetry)
308  call bc_list_add(this%bclst_vel_res, this%bc_wallmodel%symmetry)
309 
310  !Initialize bcs for u, v, w velocity components
311  call bc_list_init(this%bclst_du)
312  call bc_list_add(this%bclst_du, this%bc_sym%bc_x)
313  call bc_list_add(this%bclst_du, this%bc_sh%symmetry%bc_x)
314  call bc_list_add(this%bclst_du, this%bc_wallmodel%symmetry%bc_x)
315  call bc_list_add(this%bclst_du, this%bc_vel_res_non_normal%bc_x)
316  call bc_list_add(this%bclst_du, this%bc_vel_res)
317  call bc_list_add(this%bclst_du, this%bc_field_dirichlet_u)
318 
319  call bc_list_init(this%bclst_dv)
320  call bc_list_add(this%bclst_dv, this%bc_sym%bc_y)
321  call bc_list_add(this%bclst_dv, this%bc_sh%symmetry%bc_y)
322  call bc_list_add(this%bclst_dv, this%bc_wallmodel%symmetry%bc_y)
323  call bc_list_add(this%bclst_dv, this%bc_vel_res_non_normal%bc_y)
324  call bc_list_add(this%bclst_dv, this%bc_vel_res)
325  call bc_list_add(this%bclst_dv, this%bc_field_dirichlet_v)
326 
327  call bc_list_init(this%bclst_dw)
328  call bc_list_add(this%bclst_dw, this%bc_sym%bc_z)
329  call bc_list_add(this%bclst_dw, this%bc_sh%symmetry%bc_z)
330  call bc_list_add(this%bclst_dw, this%bc_wallmodel%symmetry%bc_z)
331  call bc_list_add(this%bclst_dw, this%bc_vel_res_non_normal%bc_z)
332  call bc_list_add(this%bclst_dw, this%bc_vel_res)
333  call bc_list_add(this%bclst_dw, this%bc_field_dirichlet_w)
334 
335  !Intialize projection space thingy
336 
337  if (this%variable_material_properties .and. &
338  this%vel_projection_dim .gt. 0) then
339  call neko_error("Velocity projection not available for full stress &
340  &formulation")
341  end if
342 
343 
344  call this%proj_prs%init(this%dm_Xh%size(), this%pr_projection_dim, &
345  this%pr_projection_activ_step)
346 
347  call this%proj_u%init(this%dm_Xh%size(), this%vel_projection_dim, &
348  this%vel_projection_activ_step)
349  call this%proj_v%init(this%dm_Xh%size(), this%vel_projection_dim, &
350  this%vel_projection_activ_step)
351  call this%proj_w%init(this%dm_Xh%size(), this%vel_projection_dim, &
352  this%vel_projection_activ_step)
353 
354 
355  ! Add lagged term to checkpoint
356  call this%chkp%add_lag(this%ulag, this%vlag, this%wlag)
357 
358  ! Determine the time-interpolation scheme
359  call json_get_or_default(params, 'case.numerics.oifs', this%oifs, .false.)
360 
361  ! Initialize the advection factory
362  call advection_factory(this%adv, params, this%c_Xh, &
363  this%ulag, this%vlag, this%wlag, &
364  this%chkp%dtlag, this%chkp%tlag, time_scheme)
365 
366  if (params%valid_path('case.fluid.flow_rate_force')) then
367  call this%vol_flow%init(this%dm_Xh, params)
368  end if
369 
370  end subroutine fluid_pnpn_init
371 
372  subroutine fluid_pnpn_restart(this, dtlag, tlag)
373  class(fluid_pnpn_t), target, intent(inout) :: this
374  real(kind=rp) :: dtlag(10), tlag(10)
375  type(field_t) :: u_temp, v_temp, w_temp
376  integer :: i, n
377 
378  n = this%u%dof%size()
379  if (allocated(this%chkp%previous_mesh%elements) .or. &
380  this%chkp%previous_Xh%lx .ne. this%Xh%lx) then
381  call col2(this%u%x, this%c_Xh%mult, this%u%dof%size())
382  call col2(this%v%x, this%c_Xh%mult, this%u%dof%size())
383  call col2(this%w%x, this%c_Xh%mult, this%u%dof%size())
384  call col2(this%p%x, this%c_Xh%mult, this%u%dof%size())
385  do i = 1, this%ulag%size()
386  call col2(this%ulag%lf(i)%x, this%c_Xh%mult, this%u%dof%size())
387  call col2(this%vlag%lf(i)%x, this%c_Xh%mult, this%u%dof%size())
388  call col2(this%wlag%lf(i)%x, this%c_Xh%mult, this%u%dof%size())
389  end do
390  end if
391 
392  if (neko_bcknd_device .eq. 1) then
393  associate(u => this%u, v => this%v, w => this%w, &
394  ulag => this%ulag, vlag => this%vlag, wlag => this%wlag,&
395  p => this%p)
396  call device_memcpy(u%x, u%x_d, u%dof%size(), &
397  host_to_device, sync = .false.)
398  call device_memcpy(v%x, v%x_d, v%dof%size(), &
399  host_to_device, sync = .false.)
400  call device_memcpy(w%x, w%x_d, w%dof%size(), &
401  host_to_device, sync = .false.)
402  call device_memcpy(p%x, p%x_d, p%dof%size(), &
403  host_to_device, sync = .false.)
404  call device_memcpy(ulag%lf(1)%x, ulag%lf(1)%x_d, &
405  u%dof%size(), host_to_device, sync = .false.)
406  call device_memcpy(ulag%lf(2)%x, ulag%lf(2)%x_d, &
407  u%dof%size(), host_to_device, sync = .false.)
408 
409  call device_memcpy(vlag%lf(1)%x, vlag%lf(1)%x_d, &
410  v%dof%size(), host_to_device, sync = .false.)
411  call device_memcpy(vlag%lf(2)%x, vlag%lf(2)%x_d, &
412  v%dof%size(), host_to_device, sync = .false.)
413 
414  call device_memcpy(wlag%lf(1)%x, wlag%lf(1)%x_d, &
415  w%dof%size(), host_to_device, sync = .false.)
416  call device_memcpy(wlag%lf(2)%x, wlag%lf(2)%x_d, &
417  w%dof%size(), host_to_device, sync = .false.)
418  call device_memcpy(this%abx1%x, this%abx1%x_d, &
419  w%dof%size(), host_to_device, sync = .false.)
420  call device_memcpy(this%abx2%x, this%abx2%x_d, &
421  w%dof%size(), host_to_device, sync = .false.)
422  call device_memcpy(this%aby1%x, this%aby1%x_d, &
423  w%dof%size(), host_to_device, sync = .false.)
424  call device_memcpy(this%aby2%x, this%aby2%x_d, &
425  w%dof%size(), host_to_device, sync = .false.)
426  call device_memcpy(this%abz1%x, this%abz1%x_d, &
427  w%dof%size(), host_to_device, sync = .false.)
428  call device_memcpy(this%abz2%x, this%abz2%x_d, &
429  w%dof%size(), host_to_device, sync = .false.)
430  call device_memcpy(this%advx%x, this%advx%x_d, &
431  w%dof%size(), host_to_device, sync = .false.)
432  call device_memcpy(this%advy%x, this%advy%x_d, &
433  w%dof%size(), host_to_device, sync = .false.)
434  call device_memcpy(this%advz%x, this%advz%x_d, &
435  w%dof%size(), host_to_device, sync = .false.)
436  end associate
437  end if
438  ! Make sure that continuity is maintained (important for interpolation)
439  ! Do not do this for lagged rhs
440  ! (derivatives are not necessairly coninous across elements)
441 
442  if (allocated(this%chkp%previous_mesh%elements) &
443  .or. this%chkp%previous_Xh%lx .ne. this%Xh%lx) then
444  call this%gs_Xh%op(this%u, gs_op_add)
445  call this%gs_Xh%op(this%v, gs_op_add)
446  call this%gs_Xh%op(this%w, gs_op_add)
447  call this%gs_Xh%op(this%p, gs_op_add)
448 
449  do i = 1, this%ulag%size()
450  call this%gs_Xh%op(this%ulag%lf(i), gs_op_add)
451  call this%gs_Xh%op(this%vlag%lf(i), gs_op_add)
452  call this%gs_Xh%op(this%wlag%lf(i), gs_op_add)
453  end do
454  end if
455  !! If we would decide to only restart from lagged fields instead of saving
456  !! abx1, aby1 etc.
457  !! Observe that one also needs to recompute the focing at the old time steps
458  !u_temp = this%ulag%lf(2)
459  !v_temp = this%vlag%lf(2)
460  !w_temp = this%wlag%lf(2)
461  !! Compute the source terms
462  !call this%source_term%compute(tlag(2), -1)
463  !
464  !! Pre-multiply the source terms with the mass matrix.
465  !if (NEKO_BCKND_DEVICE .eq. 1) then
466  ! call device_opcolv(this%f_x%x_d, this%f_y%x_d, this%f_z%x_d, &
467  ! this%c_Xh%B_d, this%msh%gdim, n)
468  !else
469  ! call opcolv(this%f_x%x, this%f_y%x, this%f_z%x, &
470  ! this%c_Xh%B, this%msh%gdim, n)
471  !end if
472 
473  !! Add the advection operators to the right-hand-side.
474  !call this%adv%compute(u_temp, v_temp, w_temp, &
475  ! this%f_x%x, this%f_y%x, this%f_z%x, &
476  ! this%Xh, this%c_Xh, this%dm_Xh%size())
477  !this%abx2 = this%f_x
478  !this%aby2 = this%f_y
479  !this%abz2 = this%f_z
480  !
481  !u_temp = this%ulag%lf(1)
482  !v_temp = this%vlag%lf(1)
483  !w_temp = this%wlag%lf(1)
484  !call this%source_term%compute(tlag(1), 0)
485 
486  !! Pre-multiply the source terms with the mass matrix.
487  !if (NEKO_BCKND_DEVICE .eq. 1) then
488  ! call device_opcolv(this%f_x%x_d, this%f_y%x_d, this%f_z%x_d, &
489  ! this%c_Xh%B_d, this%msh%gdim, n)
490  !else
491  ! call opcolv(this%f_x%x, this%f_y%x, this%f_z%x, &
492  ! this%c_Xh%B, this%msh%gdim, n)
493  !end if
494 
495  !! Pre-multiply the source terms with the mass matrix.
496  !if (NEKO_BCKND_DEVICE .eq. 1) then
497  ! call device_opcolv(this%f_x%x_d, this%f_y%x_d, this%f_z%x_d, &
498  ! this%c_Xh%B_d, this%msh%gdim, n)
499  !else
500  ! call opcolv(this%f_x%x, this%f_y%x, this%f_z%x, &
501  ! this%c_Xh%B, this%msh%gdim, n)
502  !end if
503 
504  !call this%adv%compute(u_temp, v_temp, w_temp, &
505  ! this%f_x%x, this%f_y%x, this%f_z%x, &
506  ! this%Xh, this%c_Xh, this%dm_Xh%size())
507  !this%abx1 = this%f_x
508  !this%aby1 = this%f_y
509  !this%abz1 = this%f_z
510 
511  end subroutine fluid_pnpn_restart
512 
513  subroutine fluid_pnpn_free(this)
514  class(fluid_pnpn_t), intent(inout) :: this
515 
516  !Deallocate velocity and pressure fields
517  call this%scheme_free()
518 
519  call this%bc_prs_surface%free()
520  call this%bc_sym_surface%free()
521  call bc_list_free(this%bclst_vel_res)
522  call bc_list_free(this%bclst_dp)
523  call this%proj_prs%free()
524  call this%proj_u%free()
525  call this%proj_v%free()
526  call this%proj_w%free()
527 
528  call this%p_res%free()
529  call this%u_res%free()
530  call this%v_res%free()
531  call this%w_res%free()
532 
533  call this%du%free()
534  call this%dv%free()
535  call this%dw%free()
536  call this%dp%free()
537 
538  call this%abx1%free()
539  call this%aby1%free()
540  call this%abz1%free()
541 
542  call this%abx2%free()
543  call this%aby2%free()
544  call this%abz2%free()
545 
546  call this%advx%free()
547  call this%advy%free()
548  call this%advz%free()
549 
550  if (allocated(this%Ax_vel)) then
551  deallocate(this%Ax_vel)
552  end if
553 
554  if (allocated(this%Ax_prs)) then
555  deallocate(this%Ax_prs)
556  end if
557 
558  if (allocated(this%prs_res)) then
559  deallocate(this%prs_res)
560  end if
561 
562  if (allocated(this%vel_res)) then
563  deallocate(this%vel_res)
564  end if
565 
566  if (allocated(this%sumab)) then
567  deallocate(this%sumab)
568  end if
569 
570  if (allocated(this%makeabf)) then
571  deallocate(this%makeabf)
572  end if
573 
574  if (allocated(this%makebdf)) then
575  deallocate(this%makebdf)
576  end if
577 
578  if (allocated(this%makeoifs)) then
579  deallocate(this%makeoifs)
580  end if
581 
582  call this%vol_flow%free()
583 
584  end subroutine fluid_pnpn_free
585 
592  subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
593  class(fluid_pnpn_t), target, intent(inout) :: this
594  real(kind=rp), intent(inout) :: t
595  integer, intent(inout) :: tstep
596  real(kind=rp), intent(in) :: dt
597  type(time_scheme_controller_t), intent(inout) :: ext_bdf
598  type(time_step_controller_t), intent(in) :: dt_controller
599  ! number of degrees of freedom
600  integer :: n
601  ! Solver results monitors (pressure + 3 velocity)
602  type(ksp_monitor_t) :: ksp_results(4)
603  ! Extrapolated velocity for the pressure residual
604  type(field_t), pointer :: u_e, v_e, w_e
605  ! Indices for tracking temporary fields
606  integer :: temp_indices(3)
607 
608  if (this%freeze) return
609 
610  n = this%dm_Xh%size()
611 
612  call profiler_start_region('Fluid', 1)
613  associate(u => this%u, v => this%v, w => this%w, p => this%p, &
614  du => this%du, dv => this%dv, dw => this%dw, dp => this%dp, &
615  u_res => this%u_res, v_res => this%v_res, w_res => this%w_res, &
616  p_res => this%p_res, ax_vel => this%Ax_vel, ax_prs => this%Ax_prs, &
617  xh => this%Xh, &
618  c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
619  ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
620  msh => this%msh, prs_res => this%prs_res, &
621  source_term => this%source_term, vel_res => this%vel_res, &
622  sumab => this%sumab, makeoifs => this%makeoifs, &
623  makeabf => this%makeabf, makebdf => this%makebdf, &
624  vel_projection_dim => this%vel_projection_dim, &
625  pr_projection_dim => this%pr_projection_dim, &
626  rho => this%rho, mu => this%mu, oifs => this%oifs, &
627  rho_field => this%rho_field, mu_field => this%mu_field, &
628  f_x => this%f_x, f_y => this%f_y, f_z => this%f_z, &
629  if_variable_dt => dt_controller%if_variable_dt, &
630  dt_last_change => dt_controller%dt_last_change)
631 
632  ! Get temporary arrays
633  call this%scratch%request_field(u_e, temp_indices(1))
634  call this%scratch%request_field(v_e, temp_indices(2))
635  call this%scratch%request_field(w_e, temp_indices(3))
636  call sumab%compute_fluid(u_e, v_e, w_e, u, v, w, &
637  ulag, vlag, wlag, ext_bdf%advection_coeffs, ext_bdf%nadv)
638 
639  ! Compute the source terms
640  call this%source_term%compute(t, tstep)
641 
642  ! Add Neumann bc contributions to the RHS
643  call bc_list_apply_vector(this%bclst_vel_neumann, f_x%x, f_y%x, f_z%x, &
644  this%dm_Xh%size(), t, tstep)
645 
646  ! Compute the grandient jump penalty term
647  if (this%if_gradient_jump_penalty .eqv. .true.) then
648  call this%gradient_jump_penalty_u%compute(u, v, w, u)
649  call this%gradient_jump_penalty_v%compute(u, v, w, v)
650  call this%gradient_jump_penalty_w%compute(u, v, w, w)
651  call this%gradient_jump_penalty_u%perform(f_x)
652  call this%gradient_jump_penalty_v%perform(f_y)
653  call this%gradient_jump_penalty_w%perform(f_z)
654  end if
655 
656  if (oifs) then
657  ! Add the advection operators to the right-hand-side.
658  call this%adv%compute(u, v, w, &
659  this%advx, this%advy, this%advz, &
660  xh, this%c_Xh, dm_xh%size(), dt)
661 
662  ! At this point the RHS contains the sum of the advection operator and
663  ! additional source terms, evaluated using the velocity field from the
664  ! previous time-step. Now, this value is used in the explicit time
665  ! scheme to advance both terms in time.
666  call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
667  this%abx2, this%aby2, this%abz2, &
668  f_x%x, f_y%x, f_z%x, &
669  rho, ext_bdf%advection_coeffs, n)
670 
671  ! Now, the source terms from the previous time step are added to the RHS.
672  call makeoifs%compute_fluid(this%advx%x, this%advy%x, this%advz%x, &
673  f_x%x, f_y%x, f_z%x, &
674  rho, dt, n)
675  else
676  ! Add the advection operators to the right-hand-side.
677  call this%adv%compute(u, v, w, &
678  f_x, f_y, f_z, &
679  xh, this%c_Xh, dm_xh%size())
680 
681  ! At this point the RHS contains the sum of the advection operator and
682  ! additional source terms, evaluated using the velocity field from the
683  ! previous time-step. Now, this value is used in the explicit time
684  ! scheme to advance both terms in time.
685  call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
686  this%abx2, this%aby2, this%abz2, &
687  f_x%x, f_y%x, f_z%x, &
688  rho, ext_bdf%advection_coeffs, n)
689 
690  ! Add the RHS contributions coming from the BDF scheme.
691  call makebdf%compute_fluid(ulag, vlag, wlag, f_x%x, f_y%x, f_z%x, &
692  u, v, w, c_xh%B, rho, dt, &
693  ext_bdf%diffusion_coeffs, ext_bdf%ndiff, n)
694  end if
695 
696  call ulag%update()
697  call vlag%update()
698  call wlag%update()
699 
703  call this%user_field_bc_vel%update(this%user_field_bc_vel%field_list, &
704  this%user_field_bc_vel%bc_list, this%c_Xh, t, tstep, "fluid")
705 
706  call this%bc_apply_vel(t, tstep)
707  call this%bc_apply_prs(t, tstep)
708 
709  ! Update material properties if necessary
710  call this%update_material_properties()
711 
712  ! Compute pressure.
713  call profiler_start_region('Pressure_residual', 18)
714  call prs_res%compute(p, p_res,&
715  u, v, w, &
716  u_e, v_e, w_e, &
717  f_x, f_y, f_z, &
718  c_xh, gs_xh, &
719  this%bc_prs_surface, this%bc_sym_surface,&
720  ax_prs, ext_bdf%diffusion_coeffs(1), dt, &
721  mu_field, rho_field)
722 
723  call gs_xh%op(p_res, gs_op_add)
724  call bc_list_apply_scalar(this%bclst_dp, p_res%x, p%dof%size(), t, tstep)
725  call profiler_end_region('Pressure_residual', 18)
726 
727  call this%proj_prs%pre_solving(p_res%x, tstep, c_xh, n, dt_controller, &
728  'Pressure')
729 
730  call this%pc_prs%update()
731  call profiler_start_region('Pressure_solve', 3)
732  ksp_results(1) = &
733  this%ksp_prs%solve(ax_prs, dp, p_res%x, n, c_xh, this%bclst_dp, gs_xh)
734 
735  call profiler_end_region('Pressure_solve', 3)
736 
737  call this%proj_prs%post_solving(dp%x, ax_prs, c_xh, &
738  this%bclst_dp, gs_xh, n, tstep, dt_controller)
739 
740  call field_add2(p, dp, n)
741 
742  ! Compute velocity.
743  call profiler_start_region('Velocity_residual', 19)
744  call vel_res%compute(ax_vel, u, v, w, &
745  u_res, v_res, w_res, &
746  p, &
747  f_x, f_y, f_z, &
748  c_xh, msh, xh, &
749  mu_field, rho_field, ext_bdf%diffusion_coeffs(1), &
750  dt, dm_xh%size())
751 
752  call gs_xh%op(u_res, gs_op_add)
753  call gs_xh%op(v_res, gs_op_add)
754  call gs_xh%op(w_res, gs_op_add)
755 
756  call bc_list_apply_vector(this%bclst_vel_res,&
757  u_res%x, v_res%x, w_res%x, dm_xh%size(),&
758  t, tstep)
759 
760  ! We should implement a bc that takes three field_bcs and implements
761  ! vector_apply
762  if (neko_bcknd_device .eq. 1) then
763  call this%bc_field_dirichlet_u%apply_scalar_dev(u_res%x_d, t, tstep)
764  call this%bc_field_dirichlet_v%apply_scalar_dev(v_res%x_d, t, tstep)
765  call this%bc_field_dirichlet_w%apply_scalar_dev(w_res%x_d, t, tstep)
766  else
767  call this%bc_field_dirichlet_u%apply_scalar(u_res%x, n, t, tstep)
768  call this%bc_field_dirichlet_v%apply_scalar(v_res%x, n, t, tstep)
769  call this%bc_field_dirichlet_w%apply_scalar(w_res%x, n, t, tstep)
770  end if
771 
772  call profiler_end_region('Velocity_residual', 19)
773 
774  call this%proj_u%pre_solving(u_res%x, tstep, c_xh, n, dt_controller)
775  call this%proj_v%pre_solving(v_res%x, tstep, c_xh, n, dt_controller)
776  call this%proj_w%pre_solving(w_res%x, tstep, c_xh, n, dt_controller)
777 
778  call this%pc_vel%update()
779 
780  call profiler_start_region("Velocity_solve", 4)
781  ksp_results(2:4) = this%ksp_vel%solve_coupled(ax_vel, du, dv, dw, &
782  u_res%x, v_res%x, w_res%x, n, c_xh, &
783  this%bclst_du, this%bclst_dv, this%bclst_dw, gs_xh, &
784  this%ksp_vel%max_iter)
785  call profiler_end_region("Velocity_solve", 4)
786 
787  call this%proj_u%post_solving(du%x, ax_vel, c_xh, &
788  this%bclst_du, gs_xh, n, tstep, dt_controller)
789  call this%proj_v%post_solving(dv%x, ax_vel, c_xh, &
790  this%bclst_dv, gs_xh, n, tstep, dt_controller)
791  call this%proj_w%post_solving(dw%x, ax_vel, c_xh, &
792  this%bclst_dw, gs_xh, n, tstep, dt_controller)
793 
794  if (neko_bcknd_device .eq. 1) then
795  call device_opadd2cm(u%x_d, v%x_d, w%x_d, &
796  du%x_d, dv%x_d, dw%x_d, 1.0_rp, n, msh%gdim)
797  else
798  call opadd2cm(u%x, v%x, w%x, du%x, dv%x, dw%x, 1.0_rp, n, msh%gdim)
799  end if
800 
801  if (this%forced_flow_rate) then
802  call this%vol_flow%adjust( u, v, w, p, u_res, v_res, w_res, p_res, &
803  c_xh, gs_xh, ext_bdf, rho, mu, dt, &
804  this%bclst_dp, this%bclst_du, this%bclst_dv, &
805  this%bclst_dw, this%bclst_vel_res, ax_vel, ax_prs, this%ksp_prs, &
806  this%ksp_vel, this%pc_prs, this%pc_vel, this%ksp_prs%max_iter, &
807  this%ksp_vel%max_iter)
808  end if
809 
810  call fluid_step_info(tstep, t, dt, ksp_results)
811 
812  call this%scratch%relinquish_field(temp_indices)
813 
814  end associate
815  call profiler_end_region('Fluid', 1)
816  end subroutine fluid_pnpn_step
817 
818 
819 
820 end module fluid_pnpn
Copy data between host and device (or device and device)
Definition: device.F90:51
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition: json_utils.f90:54
Retrieves a parameter by name or throws an error.
Definition: json_utils.f90:45
Subroutines to add advection terms to the RHS of a transport equation.
Definition: advection.f90:34
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.f90:34
subroutine, public bc_list_add(bclst, bc)
Add a condition to a list of boundary conditions.
Definition: bc.f90:505
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_init(bclst, size)
Constructor for a list of boundary conditions.
Definition: bc.f90:464
subroutine, public bc_list_free(bclst)
Destructor for a list of boundary conditions.
Definition: bc.f90:491
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
subroutine, public device_opadd2cm(a1_d, a2_d, a3_d, b1_d, b2_d, b3_d, c, n, gdim)
subroutine, public device_opcolv(a1_d, a2_d, a3_d, c_d, gdim, n)
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
Defines a dirichlet boundary condition.
Definition: dirichlet.f90:34
Dirichlet condition applied in the facet normal direction.
subroutine, public field_add2(a, b, n)
Vector addition .
Definition: field_math.f90:255
subroutine, public field_copy(a, b, n)
Copy a vector .
Definition: field_math.f90:126
Defines a field.
Definition: field.f90:34
Auxiliary routines for fluid solvers.
Definition: fluid_aux.f90:2
subroutine, public fluid_step_info(step, t, dt, ksp_results)
Prints for prs, velx, vely, velz the following: Number of iterations, start residual,...
Definition: fluid_aux.f90:17
Modular version of the Classic Nek5000 Pn/Pn formulation for fluids.
Definition: fluid_pnpn.f90:34
subroutine fluid_pnpn_init(this, msh, lx, params, user, time_scheme)
Definition: fluid_pnpn.f90:145
subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
Advance fluid simulation in time.
Definition: fluid_pnpn.f90:593
subroutine fluid_pnpn_free(this)
Definition: fluid_pnpn.f90:514
subroutine fluid_pnpn_restart(this, dtlag, tlag)
Definition: fluid_pnpn.f90:373
Fluid formulations.
Defines Gather-scatter operations.
Definition: gs_ops.f90:34
integer, parameter, public gs_op_add
Definition: gs_ops.f90:36
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
Implements the base abstract type for Krylov solvers plus helper types.
Definition: krylov.f90:34
Definition: math.f90:60
real(kind=rp) function, public glsum(a, n)
Sum a vector of length n.
Definition: math.f90:360
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:729
Collection of vector field operations operating on and . Note that in general the indices and ....
Definition: mathops.f90:65
subroutine, public opcolv(a1, a2, a3, c, gdim, n)
Definition: mathops.f90:97
subroutine, public opadd2cm(a1, a2, a3, b1, b2, b3, c, n, gdim)
Definition: mathops.f90:142
Defines a mesh.
Definition: mesh.f90:34
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
Dirichlet condition on axis aligned plane in the non normal direction.
Definition: non_normal.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Defines Pressure and velocity residuals in the Pn-Pn formulation.
Definition: pnpn_res.f90:34
Profiling interface.
Definition: profiler.F90:34
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
Definition: profiler.F90:78
subroutine, public profiler_end_region(name, region_id)
End the most recently started profiler region.
Definition: profiler.F90:109
Project x onto X, the space of old solutions and back again.
Definition: projection.f90:63
Routines to generate the right-hand sides for the convection-diffusion equation. Employs the EXT/BDF ...
Definition: rhs_maker.f90:38
Implements the source_term_t type and a wrapper source_term_wrapper_t.
Definition: source_term.f90:34
Compound scheme for the advection and diffusion operators in a transport equation.
Base class for time integration schemes.
Definition: time_scheme.f90:61
Implements type time_step_controller.
Interfaces for user interaction with NEKO.
Definition: user_intf.f90:34
Utilities.
Definition: utils.f90:35
Base abstract type for computing the advection operator.
Definition: advection.f90:46
Base type for a matrix-vector product providing .
Definition: ax.f90:43
A list of boundary conditions.
Definition: bc.f90:104
Generic Dirichlet boundary condition on .
Definition: dirichlet.f90:44
Dirichlet condition in facet normal direction.
Base type of all fluid formulations.
Type for storing initial and final residuals in a Krylov solver.
Definition: krylov.f90:56
Dirichlet condition in non normal direction of a plane.
Definition: non_normal.f90:44
Abstract type to compute pressure residual.
Definition: pnpn_res.f90:47
Abstract type to compute velocity residual.
Definition: pnpn_res.f90:53
Abstract type to add contributions to F from lagged BD terms.
Definition: rhs_maker.f90:59
Abstract type to sum up contributions to kth order extrapolation scheme.
Definition: rhs_maker.f90:52
Abstract type to add contributions of kth order OIFS scheme.
Definition: rhs_maker.f90:66
Abstract type to compute extrapolated velocity field for the pressure equation.
Definition: rhs_maker.f90:46
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
Provides a tool to set time step dt.