Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
fluid_pnpn.f90
Go to the documentation of this file.
1! Copyright (c) 2022-2025, 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!
35 use, intrinsic :: iso_fortran_env, only: error_unit
36 use coefs, only : coef_t
37 use symmetry, only : symmetry_t
39 use logger, only: neko_log, log_size
40 use num_types, only : rp
41 use krylov, only : ksp_monitor_t
43 pnpn_prs_res_factory, pnpn_vel_res_factory, &
44 pnpn_prs_res_stress_factory, pnpn_vel_res_stress_factory
46 rhs_maker_oifs_t, rhs_maker_sumab_fctry, rhs_maker_bdf_fctry, &
47 rhs_maker_ext_fctry, rhs_maker_oifs_fctry
51 use fluid_aux, only : fluid_step_info
52 use projection, only : projection_t
56 use advection, only : advection_t, advection_factory
58 use json_module, only : json_file, json_core, json_value
61 use json_module, only : json_file
62 use ax_product, only : ax_t, ax_helm_factory
63 use field, only : field_t
64 use dirichlet, only : dirichlet_t
68 use non_normal, only : non_normal_t
69 use checkpoint, only : chkp_t
70 use mesh, only : mesh_t
71 use user_intf, only : user_t
73 use gs_ops, only : gs_op_add
75 use mathops, only : opadd2cm, opcolv
76 use bc_list, only: bc_list_t
80 use bc, only : bc_t
81 use file, only : file_t
82 use operators, only : ortho
83 use time_state, only : time_state_t
84 use comm, only : neko_comm
85 use mpi_f08, only : mpi_allreduce, mpi_in_place, mpi_max, mpi_lor, &
86 mpi_integer, mpi_logical
87 implicit none
88 private
89
90
91
93
95 type(field_t) :: p_res, u_res, v_res, w_res
96
99 type(field_t) :: dp, du, dv, dw
100
101 !
102 ! Implicit operators, i.e. the left-hand-side of the Helmholz problem.
103 !
104
105 ! Coupled Helmholz operator for velocity
106 class(ax_t), allocatable :: ax_vel
107 ! Helmholz operator for pressure
108 class(ax_t), allocatable :: ax_prs
109
110 !
111 ! Projections for solver speed-up
112 !
113
115 type(projection_t) :: proj_prs
116 type(projection_vel_t) :: proj_vel
117
118 !
119 ! Special Karniadakis scheme boundary conditions in the pressure equation
120 !
121
123 type(facet_normal_t) :: bc_prs_surface
124
126 type(facet_normal_t) :: bc_sym_surface
127
128 !
129 ! Boundary conditions and lists for residuals and solution increments
130 !
131
133 type(zero_dirichlet_t) :: bc_vel_res
135 type(zero_dirichlet_t) :: bc_du
137 type(zero_dirichlet_t) :: bc_dv
139 type(zero_dirichlet_t) :: bc_dw
141 type(zero_dirichlet_t) :: bc_dp
142
144 type(bc_list_t) :: bclst_vel_res
145 type(bc_list_t) :: bclst_du
146 type(bc_list_t) :: bclst_dv
147 type(bc_list_t) :: bclst_dw
148 type(bc_list_t) :: bclst_dp
149
150
151 ! Checker for wether we have a strong pressure bc. If not, the pressure
152 ! is demeaned at every time step.
153 logical :: prs_dirichlet = .false.
154
155
156 ! The advection operator.
157 class(advection_t), allocatable :: adv
158
159 ! Time OIFS interpolation scheme for advection.
160 logical :: oifs
161
162 ! Time variables
163 type(field_t) :: abx1, aby1, abz1
164 type(field_t) :: abx2, aby2, abz2
165
166 ! Advection terms for the oifs method
167 type(field_t) :: advx, advy, advz
168
170 class(pnpn_prs_res_t), allocatable :: prs_res
171
173 class(pnpn_vel_res_t), allocatable :: vel_res
174
176 class(rhs_maker_sumab_t), allocatable :: sumab
177
179 class(rhs_maker_ext_t), allocatable :: makeabf
180
182 class(rhs_maker_bdf_t), allocatable :: makebdf
183
185 class(rhs_maker_oifs_t), allocatable :: makeoifs
186
188 type(fluid_volflow_t) :: vol_flow
189
191 logical :: full_stress_formulation = .false.
192
193 contains
195 procedure, pass(this) :: init => fluid_pnpn_init
197 procedure, pass(this) :: free => fluid_pnpn_free
199 procedure, pass(this) :: step => fluid_pnpn_step
201 procedure, pass(this) :: restart => fluid_pnpn_restart
203 procedure, pass(this) :: setup_bcs => fluid_pnpn_setup_bcs
205 procedure, pass(this) :: write_boundary_conditions => &
207 end type fluid_pnpn_t
208
209 interface
210
217 module subroutine pressure_bc_factory(object, scheme, json, coef, user)
218 class(bc_t), pointer, intent(inout) :: object
219 type(fluid_pnpn_t), intent(in) :: scheme
220 type(json_file), intent(inout) :: json
221 type(coef_t), target, intent(in) :: coef
222 type(user_t), intent(in) :: user
223 end subroutine pressure_bc_factory
224 end interface
225
226 interface
227
234 module subroutine velocity_bc_factory(object, scheme, json, coef, user)
235 class(bc_t), pointer, intent(inout) :: object
236 type(fluid_pnpn_t), intent(inout) :: scheme
237 type(json_file), intent(inout) :: json
238 type(coef_t), target, intent(in) :: coef
239 type(user_t), intent(in) :: user
240 end subroutine velocity_bc_factory
241 end interface
242
243contains
244
245 subroutine fluid_pnpn_init(this, msh, lx, params, user, chkp)
246 class(fluid_pnpn_t), target, intent(inout) :: this
247 type(mesh_t), target, intent(inout) :: msh
248 integer, intent(in) :: lx
249 type(json_file), target, intent(inout) :: params
250 type(user_t), target, intent(in) :: user
251 type(chkp_t), target, intent(inout) :: chkp
252 character(len=15), parameter :: scheme = 'Modular (Pn/Pn)'
253 integer :: i
254 class(bc_t), pointer :: bc_i, vel_bc
255 real(kind=rp) :: abs_tol
256 character(len=LOG_SIZE) :: log_buf
257 integer :: ierr, integer_val, solver_maxiter
258 character(len=:), allocatable :: solver_type, precon_type
259 logical :: monitor, found
260 logical :: advection
261 type(json_file) :: numerics_params, precon_params
262
263 call this%free()
264
265 ! Initialize base class
266 call this%init_base(msh, lx, params, scheme, user, .true.)
267
268 ! Add pressure field to the registry. For this scheme it is in the same
269 ! Xh as the velocity
270 call neko_field_registry%add_field(this%dm_Xh, 'p')
271 this%p => neko_field_registry%get_field('p')
272
273 !
274 ! Select governing equations via associated residual and Ax types
275 !
276
277 call json_get(params, 'case.numerics.time_order', integer_val)
278 allocate(this%ext_bdf)
279 call this%ext_bdf%init(integer_val)
280
281 call json_get_or_default(params, "case.fluid.full_stress_formulation", &
282 this%full_stress_formulation, .false.)
283
284 if (this%full_stress_formulation .eqv. .true.) then
285 ! Setup backend dependent Ax routines
286 call ax_helm_factory(this%Ax_vel, full_formulation = .true.)
287
288 ! Setup backend dependent prs residual routines
289 call pnpn_prs_res_stress_factory(this%prs_res)
290
291 ! Setup backend dependent vel residual routines
292 call pnpn_vel_res_stress_factory(this%vel_res)
293 else
294 ! Setup backend dependent Ax routines
295 call ax_helm_factory(this%Ax_vel, full_formulation = .false.)
296
297 ! Setup backend dependent prs residual routines
298 call pnpn_prs_res_factory(this%prs_res)
299
300 ! Setup backend dependent vel residual routines
301 call pnpn_vel_res_factory(this%vel_res)
302 end if
303
304
305
306 if (params%valid_path('case.fluid.nut_field')) then
307 if (this%full_stress_formulation .eqv. .false.) then
308 call neko_error("You need to set full_stress_formulation to " // &
309 "true for the fluid to have a spatially varying " // &
310 "viscocity field.")
311 end if
312 call json_get(params, 'case.fluid.nut_field', this%nut_field_name)
313 else
314 this%nut_field_name = ""
315 end if
316
317 ! Setup Ax for the pressure
318 call ax_helm_factory(this%Ax_prs, full_formulation = .false.)
319
320
321 ! Setup backend dependent summation of AB/BDF
322 call rhs_maker_sumab_fctry(this%sumab)
323
324 ! Setup backend dependent summation of extrapolation scheme
325 call rhs_maker_ext_fctry(this%makeabf)
326
327 ! Setup backend depenent contributions to F from lagged BD terms
328 call rhs_maker_bdf_fctry(this%makebdf)
329
330 ! Setup backend dependent summations of the OIFS method
331 call rhs_maker_oifs_fctry(this%makeoifs)
332
333 ! Initialize variables specific to this plan
334 associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
335 dm_xh => this%dm_Xh, nelv => this%msh%nelv)
336
337 call this%p_res%init(dm_xh, "p_res")
338 call this%u_res%init(dm_xh, "u_res")
339 call this%v_res%init(dm_xh, "v_res")
340 call this%w_res%init(dm_xh, "w_res")
341 call this%abx1%init(dm_xh, "abx1")
342 call this%aby1%init(dm_xh, "aby1")
343 call this%abz1%init(dm_xh, "abz1")
344 call this%abx2%init(dm_xh, "abx2")
345 call this%aby2%init(dm_xh, "aby2")
346 call this%abz2%init(dm_xh, "abz2")
347 call this%advx%init(dm_xh, "advx")
348 call this%advy%init(dm_xh, "advy")
349 call this%advz%init(dm_xh, "advz")
350 end associate
351
352 call this%du%init(this%dm_Xh, 'du')
353 call this%dv%init(this%dm_Xh, 'dv')
354 call this%dw%init(this%dm_Xh, 'dw')
355 call this%dp%init(this%dm_Xh, 'dp')
356
357 ! Set up boundary conditions
358 call this%setup_bcs(user, params)
359
360 ! Check if we need to output boundaries
361 call json_get_or_default(params, 'case.output_boundary', found, .false.)
362 if (found) call this%write_boundary_conditions()
363
364 call this%proj_prs%init(this%dm_Xh%size(), this%pr_projection_dim, &
365 this%pr_projection_activ_step)
366
367 call this%proj_vel%init(this%dm_Xh%size(), this%vel_projection_dim, &
368 this%vel_projection_activ_step)
369
370
371
372 ! Determine the time-interpolation scheme
373 call json_get_or_default(params, 'case.numerics.oifs', this%oifs, .false.)
374 if (params%valid_path('case.fluid.flow_rate_force')) then
375 call this%vol_flow%init(this%dm_Xh, params)
376 end if
377
378 ! Setup pressure solver
379 call neko_log%section("Pressure solver")
380
381 call json_get_or_default(params, &
382 'case.fluid.pressure_solver.max_iterations', &
383 solver_maxiter, 800)
384 call json_get(params, 'case.fluid.pressure_solver.type', solver_type)
385 call json_get(params, 'case.fluid.pressure_solver.preconditioner.type', &
386 precon_type)
387 call json_extract_object(params, &
388 'case.fluid.pressure_solver.preconditioner', precon_params)
389 call json_get(params, 'case.fluid.pressure_solver.absolute_tolerance', &
390 abs_tol)
391 call json_get_or_default(params, 'case.fluid.pressure_solver.monitor', &
392 monitor, .false.)
393 call neko_log%message('Type : ('// trim(solver_type) // &
394 ', ' // trim(precon_type) // ')')
395 write(log_buf, '(A,ES13.6)') 'Abs tol :', abs_tol
396 call neko_log%message(log_buf)
397
398 call this%solver_factory(this%ksp_prs, this%dm_Xh%size(), &
399 solver_type, solver_maxiter, abs_tol, monitor)
400 call this%precon_factory_(this%pc_prs, this%ksp_prs, &
401 this%c_Xh, this%dm_Xh, this%gs_Xh, this%bcs_prs, &
402 precon_type, precon_params)
403 call neko_log%end_section()
404
405 ! Initialize the advection factory
406 call json_get_or_default(params, 'case.fluid.advection', advection, .true.)
407 call json_extract_object(params, 'case.numerics', numerics_params)
408 call advection_factory(this%adv, numerics_params, this%c_Xh, &
409 this%ulag, this%vlag, this%wlag, &
410 chkp%dtlag, chkp%tlag, this%ext_bdf, &
411 .not. advection)
412 ! Should be in init_base maybe?
413 this%chkp => chkp
414 ! This is probably scheme specific
415 ! Should not be init really, but more like, add fluid or something...
416 call this%chkp%init(this%u, this%v, this%w, this%p)
417
418 this%chkp%abx1 => this%abx1
419 this%chkp%abx2 => this%abx2
420 this%chkp%aby1 => this%aby1
421 this%chkp%aby2 => this%aby2
422 this%chkp%abz1 => this%abz1
423 this%chkp%abz2 => this%abz2
424 call this%chkp%add_lag(this%ulag, this%vlag, this%wlag)
425
426
427
428 call neko_log%end_section()
429
430 end subroutine fluid_pnpn_init
431
432 subroutine fluid_pnpn_restart(this, chkp)
433 class(fluid_pnpn_t), target, intent(inout) :: this
434 type(chkp_t), intent(inout) :: chkp
435 real(kind=rp) :: dtlag(10), tlag(10)
436 type(field_t) :: u_temp, v_temp, w_temp
437 integer :: i, j, n
438
439 dtlag = chkp%dtlag
440 tlag = chkp%tlag
441
442 n = this%u%dof%size()
443 if (allocated(chkp%previous_mesh%elements) .or. &
444 chkp%previous_Xh%lx .ne. this%Xh%lx) then
445 associate(u => this%u, v => this%v, w => this%w, p => this%p, &
446 c_xh => this%c_Xh, ulag => this%ulag, vlag => this%vlag, &
447 wlag => this%wlag)
448 do concurrent(j = 1:n)
449 u%x(j,1,1,1) = u%x(j,1,1,1) * c_xh%mult(j,1,1,1)
450 v%x(j,1,1,1) = v%x(j,1,1,1) * c_xh%mult(j,1,1,1)
451 w%x(j,1,1,1) = w%x(j,1,1,1) * c_xh%mult(j,1,1,1)
452 p%x(j,1,1,1) = p%x(j,1,1,1) * c_xh%mult(j,1,1,1)
453 end do
454 do i = 1, this%ulag%size()
455 do concurrent(j = 1:n)
456 ulag%lf(i)%x(j,1,1,1) = ulag%lf(i)%x(j,1,1,1) &
457 * c_xh%mult(j,1,1,1)
458 vlag%lf(i)%x(j,1,1,1) = vlag%lf(i)%x(j,1,1,1) &
459 * c_xh%mult(j,1,1,1)
460 wlag%lf(i)%x(j,1,1,1) = wlag%lf(i)%x(j,1,1,1) &
461 * c_xh%mult(j,1,1,1)
462 end do
463 end do
464 end associate
465 end if
466
467 if (neko_bcknd_device .eq. 1) then
468 associate(u => this%u, v => this%v, w => this%w, &
469 ulag => this%ulag, vlag => this%vlag, wlag => this%wlag,&
470 p => this%p)
471 call device_memcpy(u%x, u%x_d, u%dof%size(), &
472 host_to_device, sync = .false.)
473 call device_memcpy(v%x, v%x_d, v%dof%size(), &
474 host_to_device, sync = .false.)
475 call device_memcpy(w%x, w%x_d, w%dof%size(), &
476 host_to_device, sync = .false.)
477 call device_memcpy(p%x, p%x_d, p%dof%size(), &
478 host_to_device, sync = .false.)
479 call device_memcpy(ulag%lf(1)%x, ulag%lf(1)%x_d, &
480 u%dof%size(), host_to_device, sync = .false.)
481 call device_memcpy(ulag%lf(2)%x, ulag%lf(2)%x_d, &
482 u%dof%size(), host_to_device, sync = .false.)
483
484 call device_memcpy(vlag%lf(1)%x, vlag%lf(1)%x_d, &
485 v%dof%size(), host_to_device, sync = .false.)
486 call device_memcpy(vlag%lf(2)%x, vlag%lf(2)%x_d, &
487 v%dof%size(), host_to_device, sync = .false.)
488
489 call device_memcpy(wlag%lf(1)%x, wlag%lf(1)%x_d, &
490 w%dof%size(), host_to_device, sync = .false.)
491 call device_memcpy(wlag%lf(2)%x, wlag%lf(2)%x_d, &
492 w%dof%size(), host_to_device, sync = .false.)
493 call device_memcpy(this%abx1%x, this%abx1%x_d, &
494 w%dof%size(), host_to_device, sync = .false.)
495 call device_memcpy(this%abx2%x, this%abx2%x_d, &
496 w%dof%size(), host_to_device, sync = .false.)
497 call device_memcpy(this%aby1%x, this%aby1%x_d, &
498 w%dof%size(), host_to_device, sync = .false.)
499 call device_memcpy(this%aby2%x, this%aby2%x_d, &
500 w%dof%size(), host_to_device, sync = .false.)
501 call device_memcpy(this%abz1%x, this%abz1%x_d, &
502 w%dof%size(), host_to_device, sync = .false.)
503 call device_memcpy(this%abz2%x, this%abz2%x_d, &
504 w%dof%size(), host_to_device, sync = .false.)
505 call device_memcpy(this%advx%x, this%advx%x_d, &
506 w%dof%size(), host_to_device, sync = .false.)
507 call device_memcpy(this%advy%x, this%advy%x_d, &
508 w%dof%size(), host_to_device, sync = .false.)
509 call device_memcpy(this%advz%x, this%advz%x_d, &
510 w%dof%size(), host_to_device, sync = .false.)
511 end associate
512 end if
513 ! Make sure that continuity is maintained (important for interpolation)
514 ! Do not do this for lagged rhs
515 ! (derivatives are not necessairly coninous across elements)
516
517 if (allocated(chkp%previous_mesh%elements) &
518 .or. chkp%previous_Xh%lx .ne. this%Xh%lx) then
519 call this%gs_Xh%op(this%u, gs_op_add)
520 call this%gs_Xh%op(this%v, gs_op_add)
521 call this%gs_Xh%op(this%w, gs_op_add)
522 call this%gs_Xh%op(this%p, gs_op_add)
523
524 do i = 1, this%ulag%size()
525 call this%gs_Xh%op(this%ulag%lf(i), gs_op_add)
526 call this%gs_Xh%op(this%vlag%lf(i), gs_op_add)
527 call this%gs_Xh%op(this%wlag%lf(i), gs_op_add)
528 end do
529 end if
530
531 end subroutine fluid_pnpn_restart
532
533 subroutine fluid_pnpn_free(this)
534 class(fluid_pnpn_t), intent(inout) :: this
535
536 !Deallocate velocity and pressure fields
537 call this%scheme_free()
538
539 call this%bc_prs_surface%free()
540 call this%bc_sym_surface%free()
541 call this%bclst_vel_res%free()
542 call this%bclst_dp%free()
543 call this%proj_prs%free()
544 call this%proj_vel%free()
545
546 call this%p_res%free()
547 call this%u_res%free()
548 call this%v_res%free()
549 call this%w_res%free()
550
551 call this%du%free()
552 call this%dv%free()
553 call this%dw%free()
554 call this%dp%free()
555
556 call this%abx1%free()
557 call this%aby1%free()
558 call this%abz1%free()
559
560 call this%abx2%free()
561 call this%aby2%free()
562 call this%abz2%free()
563
564 call this%advx%free()
565 call this%advy%free()
566 call this%advz%free()
567
568 if (allocated(this%Ax_vel)) then
569 deallocate(this%Ax_vel)
570 end if
571
572 if (allocated(this%Ax_prs)) then
573 deallocate(this%Ax_prs)
574 end if
575
576 if (allocated(this%prs_res)) then
577 deallocate(this%prs_res)
578 end if
579
580 if (allocated(this%vel_res)) then
581 deallocate(this%vel_res)
582 end if
583
584 if (allocated(this%sumab)) then
585 deallocate(this%sumab)
586 end if
587
588 if (allocated(this%makeabf)) then
589 deallocate(this%makeabf)
590 end if
591
592 if (allocated(this%makebdf)) then
593 deallocate(this%makebdf)
594 end if
595
596 if (allocated(this%makeoifs)) then
597 deallocate(this%makeoifs)
598 end if
599
600 if (allocated(this%ext_bdf)) then
601 deallocate(this%ext_bdf)
602 end if
603
604 call this%vol_flow%free()
605
606 end subroutine fluid_pnpn_free
607
614 subroutine fluid_pnpn_step(this, time, dt_controller)
615 class(fluid_pnpn_t), target, intent(inout) :: this
616 type(time_state_t), intent(in) :: time
617 type(time_step_controller_t), intent(in) :: dt_controller
618 ! number of degrees of freedom
619 integer :: n
620 ! Solver results monitors (pressure + 3 velocity)
621 type(ksp_monitor_t) :: ksp_results(4)
622
623 type(file_t) :: dump_file
624 class(bc_t), pointer :: bc_i
625 type(non_normal_t), pointer :: bc_j
626
627 if (this%freeze) return
628
629 n = this%dm_Xh%size()
630
631 call profiler_start_region('Fluid', 1)
632 associate(u => this%u, v => this%v, w => this%w, p => this%p, &
633 u_e => this%u_e, v_e => this%v_e, w_e => this%w_e, &
634 du => this%du, dv => this%dv, dw => this%dw, dp => this%dp, &
635 u_res => this%u_res, v_res => this%v_res, w_res => this%w_res, &
636 p_res => this%p_res, ax_vel => this%Ax_vel, ax_prs => this%Ax_prs, &
637 xh => this%Xh, &
638 c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
639 ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
640 msh => this%msh, prs_res => this%prs_res, &
641 source_term => this%source_term, vel_res => this%vel_res, &
642 sumab => this%sumab, makeoifs => this%makeoifs, &
643 makeabf => this%makeabf, makebdf => this%makebdf, &
644 vel_projection_dim => this%vel_projection_dim, &
645 pr_projection_dim => this%pr_projection_dim, &
646 oifs => this%oifs, &
647 rho => this%rho, mu_tot => this%mu_tot, &
648 f_x => this%f_x, f_y => this%f_y, f_z => this%f_z, &
649 t => time%t, tstep => time%tstep, dt => time%dt, &
650 ext_bdf => this%ext_bdf, event => glb_cmd_event)
651
652 ! Extrapolate the velocity if it's not done in nut_field estimation
653 call sumab%compute_fluid(u_e, v_e, w_e, u, v, w, &
654 ulag, vlag, wlag, ext_bdf%advection_coeffs, ext_bdf%nadv)
655
656 ! Compute the source terms
657 call this%source_term%compute(time)
658
659 ! Add Neumann bc contributions to the RHS
660 call this%bcs_vel%apply_vector(f_x%x, f_y%x, f_z%x, &
661 this%dm_Xh%size(), time, strong = .false.)
662
663 if (oifs) then
664 ! Add the advection operators to the right-hand-side.
665 call this%adv%compute(u, v, w, &
666 this%advx, this%advy, this%advz, &
667 xh, this%c_Xh, dm_xh%size(), dt)
668
669 ! At this point the RHS contains the sum of the advection operator and
670 ! additional source terms, evaluated using the velocity field from the
671 ! previous time-step. Now, this value is used in the explicit time
672 ! scheme to advance both terms in time.
673
674 call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
675 this%abx2, this%aby2, this%abz2, &
676 f_x%x, f_y%x, f_z%x, &
677 rho%x(1,1,1,1), ext_bdf%advection_coeffs, n)
678
679 ! Now, the source terms from the previous time step are added to the RHS.
680 call makeoifs%compute_fluid(this%advx%x, this%advy%x, this%advz%x, &
681 f_x%x, f_y%x, f_z%x, &
682 rho%x(1,1,1,1), dt, n)
683 else
684 ! Add the advection operators to the right-hand-side.
685 call this%adv%compute(u, v, w, &
686 f_x, f_y, f_z, &
687 xh, this%c_Xh, dm_xh%size())
688
689 ! At this point the RHS contains the sum of the advection operator and
690 ! additional source terms, evaluated using the velocity field from the
691 ! previous time-step. Now, this value is used in the explicit time
692 ! scheme to advance both terms in time.
693 call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
694 this%abx2, this%aby2, this%abz2, &
695 f_x%x, f_y%x, f_z%x, &
696 rho%x(1,1,1,1), ext_bdf%advection_coeffs, n)
697
698 ! Add the RHS contributions coming from the BDF scheme.
699 call makebdf%compute_fluid(ulag, vlag, wlag, f_x%x, f_y%x, f_z%x, &
700 u, v, w, c_xh%B, rho%x(1,1,1,1), dt, &
701 ext_bdf%diffusion_coeffs, ext_bdf%ndiff, n)
702 end if
703
704 call ulag%update()
705 call vlag%update()
706 call wlag%update()
707
708 call this%bc_apply_vel(time, strong = .true.)
709 call this%bc_apply_prs(time)
710
711 ! Update material properties if necessary
712 call this%update_material_properties(time)
713
714 ! Compute pressure residual.
715 call profiler_start_region('Pressure_residual', 18)
716
717 call prs_res%compute(p, p_res,&
718 u, v, w, &
719 u_e, v_e, w_e, &
720 f_x, f_y, f_z, &
721 c_xh, gs_xh, &
722 this%bc_prs_surface, this%bc_sym_surface,&
723 ax_prs, ext_bdf%diffusion_coeffs(1), dt, &
724 mu_tot, rho, event)
725
726 ! De-mean the pressure residual when no strong pressure boundaries present
727 if (.not. this%prs_dirichlet) call ortho(p_res%x, this%glb_n_points, n)
728
729 call gs_xh%op(p_res, gs_op_add, event)
730 call device_event_sync(event)
731
732 ! Set the residual to zero at strong pressure boundaries.
733 call this%bclst_dp%apply_scalar(p_res%x, p%dof%size(), time)
734
735
736 call profiler_end_region('Pressure_residual', 18)
737
738
739 call this%proj_prs%pre_solving(p_res%x, tstep, c_xh, n, dt_controller, &
740 'Pressure')
741
742 call this%pc_prs%update()
743
744 call profiler_start_region('Pressure_solve', 3)
745
746 ! Solve for the pressure increment.
747 ksp_results(1) = &
748 this%ksp_prs%solve(ax_prs, dp, p_res%x, n, c_xh, &
749 this%bclst_dp, gs_xh)
750 ksp_results(1)%name = 'Pressure'
751
752
753 call profiler_end_region('Pressure_solve', 3)
754
755 call this%proj_prs%post_solving(dp%x, ax_prs, c_xh, &
756 this%bclst_dp, gs_xh, n, tstep, dt_controller)
757
758 ! Update the pressure with the increment. Demean if necessary.
759 call field_add2(p, dp, n)
760 if (.not. this%prs_dirichlet) call ortho(p%x, this%glb_n_points, n)
761
762 ! Compute velocity residual.
763 call profiler_start_region('Velocity_residual', 19)
764 call vel_res%compute(ax_vel, u, v, w, &
765 u_res, v_res, w_res, &
766 p, &
767 f_x, f_y, f_z, &
768 c_xh, msh, xh, &
769 mu_tot, rho, ext_bdf%diffusion_coeffs(1), &
770 dt, dm_xh%size())
771
772 call gs_xh%op(u_res, gs_op_add, event)
773 call device_event_sync(event)
774 call gs_xh%op(v_res, gs_op_add, event)
775 call device_event_sync(event)
776 call gs_xh%op(w_res, gs_op_add, event)
777 call device_event_sync(event)
778
779 ! Set residual to zero at strong velocity boundaries.
780 call this%bclst_vel_res%apply(u_res, v_res, w_res, time)
781
782
783 call profiler_end_region('Velocity_residual', 19)
784
785 call this%proj_vel%pre_solving(u_res%x, v_res%x, w_res%x, &
786 tstep, c_xh, n, dt_controller, 'Velocity')
787
788 call this%pc_vel%update()
789
790 call profiler_start_region("Velocity_solve", 4)
791 ksp_results(2:4) = this%ksp_vel%solve_coupled(ax_vel, du, dv, dw, &
792 u_res%x, v_res%x, w_res%x, n, c_xh, &
793 this%bclst_du, this%bclst_dv, this%bclst_dw, gs_xh, &
794 this%ksp_vel%max_iter)
795 call profiler_end_region("Velocity_solve", 4)
796 if (this%full_stress_formulation) then
797 ksp_results(2)%name = 'Momentum'
798 else
799 ksp_results(2)%name = 'X-Velocity'
800 ksp_results(3)%name = 'Y-Velocity'
801 ksp_results(4)%name = 'Z-Velocity'
802 end if
803
804 call this%proj_vel%post_solving(du%x, dv%x, dw%x, ax_vel, c_xh, &
805 this%bclst_du, this%bclst_dv, this%bclst_dw, gs_xh, n, tstep, &
806 dt_controller)
807
808 if (neko_bcknd_device .eq. 1) then
809 call device_opadd2cm(u%x_d, v%x_d, w%x_d, &
810 du%x_d, dv%x_d, dw%x_d, 1.0_rp, n, msh%gdim)
811 else
812 call opadd2cm(u%x, v%x, w%x, du%x, dv%x, dw%x, 1.0_rp, n, msh%gdim)
813 end if
814
815 if (this%forced_flow_rate) then
816 ! Horrible mu hack?!
817 call this%vol_flow%adjust( u, v, w, p, u_res, v_res, w_res, p_res, &
818 c_xh, gs_xh, ext_bdf, rho%x(1,1,1,1), mu_tot, &
819 dt, this%bclst_dp, this%bclst_du, this%bclst_dv, &
820 this%bclst_dw, this%bclst_vel_res, ax_vel, ax_prs, this%ksp_prs, &
821 this%ksp_vel, this%pc_prs, this%pc_vel, this%ksp_prs%max_iter, &
822 this%ksp_vel%max_iter)
823 end if
824
825 call fluid_step_info(time, ksp_results, &
826 this%full_stress_formulation, this%strict_convergence)
827
828 end associate
829 call profiler_end_region('Fluid', 1)
830 end subroutine fluid_pnpn_step
831
834 subroutine fluid_pnpn_setup_bcs(this, user, params)
835 class(fluid_pnpn_t), target, intent(inout) :: this
836 type(user_t), target, intent(in) :: user
837 type(json_file), intent(inout) :: params
838 integer :: i, n_bcs, zone_index, j, zone_size, global_zone_size, ierr
839 class(bc_t), pointer :: bc_i
840 type(json_core) :: core
841 type(json_value), pointer :: bc_object
842 type(json_file) :: bc_subdict
843 logical :: found
844 ! Monitor which boundary zones have been marked
845 logical, allocatable :: marked_zones(:)
846 integer, allocatable :: zone_indices(:)
847
848 ! Lists for the residuals and solution increments
849 call this%bclst_vel_res%init()
850 call this%bclst_du%init()
851 call this%bclst_dv%init()
852 call this%bclst_dw%init()
853 call this%bclst_dp%init()
854
855 call this%bc_vel_res%init_from_components(this%c_Xh)
856 call this%bc_du%init_from_components(this%c_Xh)
857 call this%bc_dv%init_from_components(this%c_Xh)
858 call this%bc_dw%init_from_components(this%c_Xh)
859 call this%bc_dp%init_from_components(this%c_Xh)
860
861 ! Special PnPn boundary conditions for pressure
862 call this%bc_prs_surface%init_from_components(this%c_Xh)
863 call this%bc_sym_surface%init_from_components(this%c_Xh)
864
865 ! Populate bcs_vel and bcs_prs based on the case file
866 if (params%valid_path('case.fluid.boundary_conditions')) then
867 call params%info('case.fluid.boundary_conditions', n_children = n_bcs)
868 call params%get_core(core)
869 call params%get('case.fluid.boundary_conditions', bc_object, found)
870
871 !
872 ! Velocity bcs
873 !
874 call this%bcs_vel%init(n_bcs)
875
876 allocate(marked_zones(size(this%msh%labeled_zones)))
877 marked_zones = .false.
878
879 do i = 1, n_bcs
880 ! Create a new json containing just the subdict for this bc
881 call json_extract_item(core, bc_object, i, bc_subdict)
882
883 call json_get(bc_subdict, "zone_indices", zone_indices)
884
885 ! Check that we are not trying to assing a bc to zone, for which one
886 ! has already been assigned and that the zone has more than 0 size
887 ! in the mesh.
888 do j = 1, size(zone_indices)
889 zone_size = this%msh%labeled_zones(zone_indices(j))%size
890 call mpi_allreduce(zone_size, global_zone_size, 1, &
891 mpi_integer, mpi_max, neko_comm, ierr)
892
893 if (global_zone_size .eq. 0) then
894 write(error_unit, '(A, A, I0, A, A, I0, A)') "*** ERROR ***: ",&
895 "Zone index ", zone_indices(j), &
896 " is invalid as this zone has 0 size, meaning it ", &
897 "is not in the mesh. Check fluid boundary condition ", &
898 i, "."
899 error stop
900 end if
901
902 if (marked_zones(zone_indices(j)) .eqv. .true.) then
903 write(error_unit, '(A, A, I0, A, A, A, A)') "*** ERROR ***: ", &
904 "Zone with index ", zone_indices(j), &
905 " has already been assigned a boundary condition. ", &
906 "Please check your boundary_conditions entry for the ", &
907 "fluid and make sure that each zone index appears only ", &
908 "in a single boundary condition."
909 error stop
910 else
911 marked_zones(zone_indices(j)) = .true.
912 end if
913 end do
914
915 bc_i => null()
916 call velocity_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
917
918 ! Not all bcs require an allocation for velocity in particular,
919 ! so we check.
920 if (associated(bc_i)) then
921
922 ! We need to treat mixed bcs separately because they are by
923 ! convention marked weak and currently contain nested
924 ! bcs, some of which are strong.
925 select type (bc_i)
926 type is (symmetry_t)
927 ! Symmetry has 3 internal bcs, but only one actually contains
928 ! markings.
929 ! Symmetry's apply_scalar doesn't do anything, so we need to mark
930 ! individual nested bcs to the du,dv,dw, whereas the vel_res can
931 ! just get symmetry as a whole, because on this list we call
932 ! apply_vector.
933 ! Additionally we have to mark the special surface bc for p.
934 call this%bclst_vel_res%append(bc_i)
935 call this%bc_du%mark_facets(bc_i%bc_x%marked_facet)
936 call this%bc_dv%mark_facets(bc_i%bc_y%marked_facet)
937 call this%bc_dw%mark_facets(bc_i%bc_z%marked_facet)
938
939 call this%bcs_vel%append(bc_i)
940
941 call this%bc_sym_surface%mark_facets(bc_i%marked_facet)
942 type is (non_normal_t)
943 ! This is a bc for the residuals and increments, not the
944 ! velocity itself. So, don't append to bcs_vel
945 call this%bclst_vel_res%append(bc_i)
946 call this%bc_du%mark_facets(bc_i%bc_x%marked_facet)
947 call this%bc_dv%mark_facets(bc_i%bc_y%marked_facet)
948 call this%bc_dw%mark_facets(bc_i%bc_z%marked_facet)
949 type is (shear_stress_t)
950 ! Same as symmetry
951 call this%bclst_vel_res%append(bc_i%symmetry)
952 call this%bclst_du%append(bc_i%symmetry%bc_x)
953 call this%bclst_dv%append(bc_i%symmetry%bc_y)
954 call this%bclst_dw%append(bc_i%symmetry%bc_z)
955
956 call this%bcs_vel%append(bc_i)
957 type is (wall_model_bc_t)
958 ! Same as symmetry
959 call this%bclst_vel_res%append(bc_i%symmetry)
960 call this%bclst_du%append(bc_i%symmetry%bc_x)
961 call this%bclst_dv%append(bc_i%symmetry%bc_y)
962 call this%bclst_dw%append(bc_i%symmetry%bc_z)
963
964 call this%bcs_vel%append(bc_i)
965 class default
966
967 ! For the default case we use our dummy zero_dirichlet bcs to
968 ! mark the same faces as in ordinary velocity dirichlet
969 ! conditions.
970 ! Additionally we mark the special PnPn pressure bc.
971 if (bc_i%strong .eqv. .true.) then
972 call this%bc_vel_res%mark_facets(bc_i%marked_facet)
973 call this%bc_du%mark_facets(bc_i%marked_facet)
974 call this%bc_dv%mark_facets(bc_i%marked_facet)
975 call this%bc_dw%mark_facets(bc_i%marked_facet)
976
977 call this%bc_prs_surface%mark_facets(bc_i%marked_facet)
978 end if
979
980 call this%bcs_vel%append(bc_i)
981 end select
982 end if
983 end do
984
985 ! Make sure all labeled zones with non-zero size have been marked
986 do i = 1, size(this%msh%labeled_zones)
987 if ((this%msh%labeled_zones(i)%size .gt. 0) .and. &
988 (marked_zones(i) .eqv. .false.)) then
989 write(error_unit, '(A, A, I0)') "*** ERROR ***: ", &
990 "No fluid boundary condition assigned to zone ", i
991 error stop
992 end if
993 end do
994
995 !
996 ! Pressure bcs
997 !
998 call this%bcs_prs%init(n_bcs)
999
1000 do i = 1, n_bcs
1001 ! Create a new json containing just the subdict for this bc
1002 call json_extract_item(core, bc_object, i, bc_subdict)
1003 bc_i => null()
1004 call pressure_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
1005
1006 ! Not all bcs require an allocation for pressure in particular,
1007 ! so we check.
1008 if (associated(bc_i)) then
1009 call this%bcs_prs%append(bc_i)
1010
1011 ! Mark strong bcs in the dummy dp bc to force zero change.
1012 if (bc_i%strong .eqv. .true.) then
1013 call this%bc_dp%mark_facets(bc_i%marked_facet)
1014 end if
1015
1016 end if
1017
1018 end do
1019 else
1020 ! Check that there are no labeled zones, i.e. all are periodic.
1021 do i = 1, size(this%msh%labeled_zones)
1022 if (this%msh%labeled_zones(i)%size .gt. 0) then
1023 call neko_error("No boundary_conditions entry in the case file!")
1024 end if
1025 end do
1026
1027 ! For a pure periodic case, we still need to initilise the bc lists
1028 ! to a zero size to avoid issues with apply() in step()
1029 call this%bcs_vel%init()
1030 call this%bcs_prs%init()
1031
1032 end if
1033
1034 call this%bc_prs_surface%finalize()
1035 call this%bc_sym_surface%finalize()
1036
1037 call this%bc_vel_res%finalize()
1038 call this%bc_du%finalize()
1039 call this%bc_dv%finalize()
1040 call this%bc_dw%finalize()
1041 call this%bc_dp%finalize()
1042
1043 call this%bclst_vel_res%append(this%bc_vel_res)
1044 call this%bclst_du%append(this%bc_du)
1045 call this%bclst_dv%append(this%bc_dv)
1046 call this%bclst_dw%append(this%bc_dw)
1047 call this%bclst_dp%append(this%bc_dp)
1048
1049 ! If we have no strong pressure bcs, we will demean the pressure
1050 this%prs_dirichlet = .not. this%bclst_dp%is_empty()
1051 call mpi_allreduce(mpi_in_place, this%prs_dirichlet, 1, &
1052 mpi_logical, mpi_lor, neko_comm)
1053
1054 end subroutine fluid_pnpn_setup_bcs
1055
1057 subroutine fluid_pnpn_write_boundary_conditions(this)
1058 use inflow, only : inflow_t
1060 use blasius, only : blasius_t
1062 use dong_outflow, only : dong_outflow_t
1063 class(fluid_pnpn_t), target, intent(inout) :: this
1064 type(dirichlet_t) :: bdry_mask
1065 type(field_t), pointer :: bdry_field
1066 type(file_t) :: bdry_file
1067 integer :: temp_index, i
1068 class(bc_t), pointer :: bci
1069 character(len=LOG_SIZE) :: log_buf
1070
1071 call neko_log%section("Fluid boundary conditions")
1072 write(log_buf, '(A)') 'Marking using integer keys in bdry0.f00000'
1073 call neko_log%message(log_buf)
1074 write(log_buf, '(A)') 'Condition-value pairs: '
1075 call neko_log%message(log_buf)
1076 write(log_buf, '(A)') ' no_slip = 1'
1077 call neko_log%message(log_buf)
1078 write(log_buf, '(A)') ' velocity_value = 2'
1079 call neko_log%message(log_buf)
1080 write(log_buf, '(A)') ' outflow, normal_outflow (+dong) = 3'
1081 call neko_log%message(log_buf)
1082 write(log_buf, '(A)') ' symmetry = 4'
1083 call neko_log%message(log_buf)
1084 write(log_buf, '(A)') ' periodic = 6'
1085 call neko_log%message(log_buf)
1086 write(log_buf, '(A)') ' user_velocity = 7'
1087 call neko_log%message(log_buf)
1088 write(log_buf, '(A)') ' user_pressure = 8'
1089 call neko_log%message(log_buf)
1090 write(log_buf, '(A)') ' shear_stress = 9'
1091 call neko_log%message(log_buf)
1092 write(log_buf, '(A)') ' wall_modelling = 10'
1093 call neko_log%message(log_buf)
1094 write(log_buf, '(A)') ' blasius_profile = 11'
1095 call neko_log%message(log_buf)
1096 call neko_log%end_section()
1097
1098 call this%scratch%request_field(bdry_field, temp_index)
1099 bdry_field = 0.0_rp
1100
1101
1102
1103 call bdry_mask%init_from_components(this%c_Xh, 6.0_rp)
1104 call bdry_mask%mark_zone(this%msh%periodic)
1105 call bdry_mask%finalize()
1106 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1107 call bdry_mask%free()
1108
1109 do i = 1, this%bcs_prs%size()
1110 bci => this%bcs_prs%get(i)
1111 select type (bc => bci)
1112 type is (zero_dirichlet_t)
1113 call bdry_mask%init_from_components(this%c_Xh, 3.0_rp)
1114 call bdry_mask%mark_facets(bci%marked_facet)
1115 call bdry_mask%finalize()
1116 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1117 call bdry_mask%free()
1118 type is (dong_outflow_t)
1119 call bdry_mask%init_from_components(this%c_Xh, 3.0_rp)
1120 call bdry_mask%mark_facets(bci%marked_facet)
1121 call bdry_mask%finalize()
1122 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1123 call bdry_mask%free()
1124 type is (field_dirichlet_t)
1125 call bdry_mask%init_from_components(this%c_Xh, 8.0_rp)
1126 call bdry_mask%mark_facets(bci%marked_facet)
1127 call bdry_mask%finalize()
1128 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1129 call bdry_mask%free()
1130 end select
1131 end do
1132
1133 do i = 1, this%bcs_vel%size()
1134 bci => this%bcs_vel%get(i)
1135 select type (bc => bci)
1136 type is (zero_dirichlet_t)
1137 call bdry_mask%init_from_components(this%c_Xh, 1.0_rp)
1138 call bdry_mask%mark_facets(bci%marked_facet)
1139 call bdry_mask%finalize()
1140 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1141 call bdry_mask%free()
1142 type is (inflow_t)
1143 call bdry_mask%init_from_components(this%c_Xh, 2.0_rp)
1144 call bdry_mask%mark_facets(bci%marked_facet)
1145 call bdry_mask%finalize()
1146 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1147 call bdry_mask%free()
1148 type is (symmetry_t)
1149 call bdry_mask%init_from_components(this%c_Xh, 4.0_rp)
1150 call bdry_mask%mark_facets(bci%marked_facet)
1151 call bdry_mask%finalize()
1152 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1153 call bdry_mask%free()
1154 type is (field_dirichlet_vector_t)
1155 call bdry_mask%init_from_components(this%c_Xh, 7.0_rp)
1156 call bdry_mask%mark_facets(bci%marked_facet)
1157 call bdry_mask%finalize()
1158 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1159 call bdry_mask%free()
1160 type is (shear_stress_t)
1161 call bdry_mask%init_from_components(this%c_Xh, 9.0_rp)
1162 call bdry_mask%mark_facets(bci%marked_facet)
1163 call bdry_mask%finalize()
1164 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1165 call bdry_mask%free()
1166 type is (wall_model_bc_t)
1167 call bdry_mask%init_from_components(this%c_Xh, 10.0_rp)
1168 call bdry_mask%mark_facets(bci%marked_facet)
1169 call bdry_mask%finalize()
1170 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1171 call bdry_mask%free()
1172 type is (blasius_t)
1173 call bdry_mask%init_from_components(this%c_Xh, 11.0_rp)
1174 call bdry_mask%mark_facets(bci%marked_facet)
1175 call bdry_mask%finalize()
1176 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1177 call bdry_mask%free()
1178 end select
1179 end do
1180
1181
1182 call bdry_file%init('bdry.fld')
1183 call bdry_file%write(bdry_field)
1184
1185 call this%scratch%relinquish_field(temp_index)
1186 end subroutine fluid_pnpn_write_boundary_conditions
1187
1188end module fluid_pnpn
Copy data between host and device (or device and device)
Definition device.F90:66
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Retrieves a parameter by name or throws an error.
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 list of bc_t.
Definition bc_list.f90:34
Defines a boundary condition.
Definition bc.f90:34
Defines a Blasius profile dirichlet condition.
Definition blasius.f90:34
Defines a checkpoint.
Coefficients.
Definition coef.f90:34
Definition comm.F90:1
type(mpi_comm), public neko_comm
MPI communicator.
Definition comm.F90:42
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
subroutine, public device_event_sync(event)
Synchronize an event.
Definition device.F90:1309
integer, parameter, public host_to_device
Definition device.F90:47
type(c_ptr), bind(C), public glb_cmd_event
Event for the global command queue.
Definition device.F90:57
Defines a dirichlet boundary condition.
Definition dirichlet.f90:34
Defines a dong outflow condition.
Dirichlet condition applied in the facet normal direction.
Defines inflow dirichlet conditions.
Defines user dirichlet condition for a scalar field.
subroutine, public field_add2(a, b, n)
Vector addition .
subroutine, public field_copy(a, b, n)
Copy a vector .
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines a field.
Definition field.f90:34
Module for file I/O operations.
Definition file.f90:34
Auxiliary routines for fluid solvers.
Definition fluid_aux.f90:2
subroutine, public fluid_step_info(time, ksp_results, full_stress_formulation, strict_convergence)
Prints for prs, velx, vely, velz the following: Number of iterations, start residual,...
Definition fluid_aux.f90:19
Modular version of the Classic Nek5000 Pn/Pn formulation for fluids.
subroutine fluid_pnpn_setup_bcs(this, user, params)
Sets up the boundary condition for the scheme.
subroutine fluid_pnpn_restart(this, chkp)
subroutine fluid_pnpn_write_boundary_conditions(this)
Write a field with boundary condition specifications.
subroutine fluid_pnpn_step(this, time, dt_controller)
Advance fluid simulation in time.
subroutine fluid_pnpn_free(this)
subroutine fluid_pnpn_init(this, msh, lx, params, user, chkp)
Boundary condition factory for pressure.
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Defines inflow dirichlet conditions.
Definition inflow.f90:34
Utilities for retrieving parameters from the case files.
subroutine, public json_extract_object(json, name, object)
Extract object as a separate JSON dictionary.
Implements the base abstract type for Krylov solvers plus helper types.
Definition krylov.f90:34
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
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.
integer, parameter neko_bcknd_device
Dirichlet condition on axis aligned plane in the non normal direction.
integer, parameter, public dp
Definition num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
subroutine, public ortho(x, glb_n_points, n)
Othogonalize with regard to vector (1,1,1,1,1,1...,1)^T.
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 Couple projections for velocity.
Project x onto X, the space of old solutions and back again.
Routines to generate the right-hand sides for the convection-diffusion equation. Employs the EXT/BDF ...
Definition rhs_maker.f90:38
Defines a shear stress boundary condition for a vector field. Maintainer: Timofey Mukha.
Implements the source_term_t type and a wrapper source_term_wrapper_t.
Mixed Dirichlet-Neumann axis aligned symmetry plane.
Definition symmetry.f90:34
Module with things related to the simulation time.
Implements type time_step_controller.
Interfaces for user interaction with NEKO.
Definition user_intf.f90:34
Utilities.
Definition utils.f90:35
subroutine, public neko_type_error(base_type, wrong_type, known_types)
Reports an error allocating a type for a particular base pointer class.
Definition utils.f90:251
Defines the wall_model_bc_t type. Maintainer: Timofey Mukha.
Defines a zero-valued Dirichlet boundary condition.
Base abstract type for computing the advection operator.
Definition advection.f90:46
Base type for a matrix-vector product providing .
Definition ax.f90:43
Base type for a boundary condition.
Definition bc.f90:61
A list of allocatable `bc_t`. Follows the standard interface of lists.
Definition bc_list.f90:48
Blasius profile for inlet (vector valued).
Definition blasius.f90:53
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Generic Dirichlet boundary condition on .
Definition dirichlet.f90:48
Dong outflow condition Follows "A Convective-like Energy-Stable Open Boundary Condition for Simulati...
Dirichlet condition in facet normal direction.
User defined dirichlet condition, for which the user can work with an entire field....
Extension of the user defined dirichlet condition field_dirichlet
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Definition file.f90:55
Dirichlet condition for inlet (vector valued)
Definition inflow.f90:47
Type for storing initial and final residuals in a Krylov solver.
Definition krylov.f90:56
Dirichlet condition in non normal direction of a plane.
Abstract type to compute pressure residual.
Definition pnpn_res.f90:48
Abstract type to compute velocity residual.
Definition pnpn_res.f90:54
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
A shear stress boundary condition.
Mixed Dirichlet-Neumann symmetry plane condition.
Definition symmetry.f90:50
A struct that contains all info about the time, expand as needed.
A type collecting all the overridable user routines and flag to suppress type injection from custom m...
A shear stress boundary condition, computing the stress values using a wall model.
Zero-valued Dirichlet boundary condition. Used for no-slip walls, but also for various auxillary cond...