diff --git a/examples/poisson/driver.f90 b/examples/poisson/driver.f90 index 400f21350f1..70cbe9ea1fe 100644 --- a/examples/poisson/driver.f90 +++ b/examples/poisson/driver.f90 @@ -30,16 +30,16 @@ program poisson end if stop end if - - call neko_init + + call neko_init call get_command_argument(1, fname) call get_command_argument(2, lxchar) call get_command_argument(3, iterchar) read(lxchar, *) lx read(iterchar, *) niter - + nmsh_file = file_t(fname) - call nmsh_file%read(msh) + call nmsh_file%read(msh) call Xh%init(GLL, lx, lx, lx) @@ -47,20 +47,19 @@ program poisson call gs_h%init(dm) call coef%init(gs_h) - + call x%init(dm, "x") n = Xh%lx * Xh%ly * Xh%lz * msh%nelv - call dir_bc%init(coef) - call dir_bc%set_g(real(0.0d0,rp)) - + call dir_bc%init_from_components(coef, 0.0_rp) + !user specified call set_bc(dir_bc, msh) - + call dir_bc%finalize() - call bc_list_init(bclst) - call bc_list_add(bclst,dir_bc) + call bclst%init() + call bclst%append(dir_bc) call solver%init(n, niter, abs_tol = tol) allocate(f(n)) @@ -71,23 +70,23 @@ program poisson call bc_list_apply(bclst,f,n) ksp_mon = solver%solve(ax, x, f, n, coef, bclst, gs_h, niter) n_glb = Xh%lx * Xh%ly * Xh%lz * msh%glb_nelv - + call MPI_Barrier(NEKO_COMM, ierr) call set_timer_flop_cnt(0, msh%glb_nelv, x%Xh%lx, niter, n_glb, ksp_mon) ksp_mon = solver%solve(ax, x, f, n, coef, bclst, gs_h, niter) call set_timer_flop_cnt(1, msh%glb_nelv, x%Xh%lx, niter, n_glb, ksp_mon) - + fname = 'out.fld' mf = file_t(fname) call mf%write(x) deallocate(f) call solver%free() call dir_bc%free() - call bc_list_free(bclst) + call bclst%free() call Xh%free() call x%free() - call msh%free() + call msh%free() call neko_finalize end program poisson diff --git a/examples/rayleigh-benard-cylinder/rayleigh.case b/examples/rayleigh-benard-cylinder/rayleigh.case index f60ea65040c..48b86abfc88 100644 --- a/examples/rayleigh-benard-cylinder/rayleigh.case +++ b/examples/rayleigh-benard-cylinder/rayleigh.case @@ -59,7 +59,12 @@ "Pr": 1.0, "initial_condition": { "type": "user" - } + }, + "boundary_conditions": [ + { + "type": "user_pointwise" + } + ], }, "simulation_components": [ { diff --git a/examples/rayleigh-benard/rayleigh.case b/examples/rayleigh-benard/rayleigh.case index 0681be28e8b..5d9a9d3ea45 100644 --- a/examples/rayleigh-benard/rayleigh.case +++ b/examples/rayleigh-benard/rayleigh.case @@ -55,17 +55,24 @@ "scalar": { "enabled": true, "Pr": 0.71, - "boundary_types": [ - "", - "", - "", - "", - "d=1", - "d=0" + "boundary_conditions": [ + { + "type": "dirichlet", + "zone_index": 5, + "value": 1.0 + }, + { + "type": "dirichlet", + "zone_index": 6, + "value": 0.0 + }, + { + "type": "user_pointwise" + } ], "initial_condition": { "type": "user" } } } -} \ No newline at end of file +} diff --git a/src/.depends b/src/.depends index 0e7666db337..b6e846692f0 100644 --- a/src/.depends +++ b/src/.depends @@ -96,7 +96,7 @@ math/fdm.o : math/fdm.f90 comm/comm.o common/utils.o device/device.o math/bcknd/ math/bcknd/cpu/fdm_cpu.o : math/bcknd/cpu/fdm_cpu.f90 math/bcknd/cpu/tensor_cpu.o config/num_types.o math/bcknd/sx/fdm_sx.o : math/bcknd/sx/fdm_sx.f90 math/bcknd/sx/tensor_sx.o config/num_types.o math/bcknd/xsmm/fdm_xsmm.o : math/bcknd/xsmm/fdm_xsmm.f90 math/bcknd/xsmm/tensor_xsmm.o config/num_types.o -math/schwarz.o : math/schwarz.f90 config/neko_config.o device/device.o math/fdm.o math/bcknd/device/device_math.o math/bcknd/device/device_schwarz.o gs/gather_scatter.o bc/bc.o sem/dofmap.o sem/space.o mesh/mesh.o math/math.o config/num_types.o +math/schwarz.o : math/schwarz.f90 bc/bc_list.o config/neko_config.o device/device.o math/fdm.o math/bcknd/device/device_math.o math/bcknd/device/device_schwarz.o gs/gather_scatter.o bc/bc.o sem/dofmap.o sem/space.o mesh/mesh.o math/math.o config/num_types.o math/vector.o : math/vector.f90 common/utils.o math/bcknd/device/device_math.o device/device.o config/num_types.o config/neko_config.o math/matrix.o : math/matrix.f90 common/utils.o math/bcknd/device/device_math.o device/device.o config/num_types.o config/neko_config.o math/signed_distance.o : math/signed_distance.f90 adt/stack.o mesh/search_tree/aabb.o mesh/point.o common/utils.o mesh/aabb_tree.o mesh/tri_mesh.o mesh/tri.o field/field.o config/num_types.o @@ -127,41 +127,43 @@ common/global_interpolation.o : common/global_interpolation.F90 comm/mpi_types.o common/profiler.o : common/profiler.F90 common/craypat.o device/hip/roctx.o device/cuda/nvtx.o device/device.o config/neko_config.o common/craypat.o : common/craypat.F90 common/utils.o adt/stack.o bc/bc.o : bc/bc.f90 common/utils.o adt/tuple.o adt/stack.o mesh/facet_zone.o mesh/mesh.o sem/space.o sem/coef.o sem/dofmap.o device/device.o config/num_types.o config/neko_config.o -bc/dirichlet.o : bc/dirichlet.f90 bc/bc.o config/num_types.o bc/bcknd/device/device_dirichlet.o -bc/neumann.o : bc/neumann.f90 math/math.o sem/coef.o common/utils.o bc/bc.o config/num_types.o +bc/bc_list.o : bc/bc_list.f90 bc/bc.o device/device.o config/num_types.o config/neko_config.o +bc/bc_fctry.o : bc/bc_fctry.f90 common/user_intf.o bc/usr_scalar.o mesh/mesh.o mesh/facet_zone.o sem/coef.o bc/neumann.o bc/dirichlet.o common/log.o common/json_utils.o bc/bc.o +bc/dirichlet.o : bc/dirichlet.f90 common/json_utils.o sem/coef.o bc/bc.o config/num_types.o bc/bcknd/device/device_dirichlet.o +bc/neumann.o : bc/neumann.f90 math/math.o common/json_utils.o sem/coef.o common/utils.o bc/bc.o config/num_types.o bc/dong_outflow.o : bc/dong_outflow.f90 common/json_utils.o field/field_registry.o bc/bcknd/device/device_dong_outflow.o common/utils.o sem/coef.o sem/dofmap.o field/field.o bc/bc.o config/num_types.o device/device.o bc/dirichlet.o config/neko_config.o -bc/wall.o : bc/wall.f90 bc/bc.o config/num_types.o bc/bcknd/device/device_wall.o -bc/inflow.o : bc/inflow.f90 bc/bc.o config/num_types.o bc/bcknd/device/device_inflow.o -bc/field_dirichlet.o : bc/field_dirichlet.f90 sem/dofmap.o math/bcknd/device/device_math.o math/math.o field/field_list.o field/field.o common/utils.o device/device.o bc/bc.o bc/dirichlet.o sem/coef.o config/num_types.o -bc/field_dirichlet_vector.o : bc/field_dirichlet_vector.f90 bc/field_dirichlet.o sem/dofmap.o math/bcknd/device/device_math.o math/math.o field/field_list.o field/field.o common/utils.o device/device.o bc/bc.o bc/dirichlet.o sem/coef.o config/num_types.o +bc/zero_dirichlet.o : bc/zero_dirichlet.f90 sem/coef.o bc/bc.o config/num_types.o bc/bcknd/device/device_zero_dirichlet.o +bc/inflow.o : bc/inflow.f90 common/json_utils.o sem/coef.o bc/bc.o config/num_types.o bc/bcknd/device/device_inflow.o +bc/field_dirichlet.o : bc/field_dirichlet.f90 sem/dofmap.o math/bcknd/device/device_math.o math/math.o field/field_list.o field/field.o common/utils.o device/device.o bc/bc_list.o bc/bc.o bc/dirichlet.o sem/coef.o config/num_types.o +bc/field_dirichlet_vector.o : bc/field_dirichlet_vector.f90 bc/field_dirichlet.o sem/dofmap.o math/bcknd/device/device_math.o math/math.o field/field_list.o field/field.o common/utils.o device/device.o bc/bc_list.o bc/bc.o bc/dirichlet.o sem/coef.o config/num_types.o bc/usr_inflow.o : bc/usr_inflow.f90 bc/bc.o common/utils.o bc/bcknd/device/device_inhom_dirichlet.o device/device.o bc/inflow.o sem/coef.o config/num_types.o bc/usr_scalar.o : bc/usr_scalar.f90 common/utils.o bc/bcknd/device/device_inhom_dirichlet.o device/device.o bc/bc.o sem/coef.o config/num_types.o bc/facet_normal.o : bc/facet_normal.f90 common/utils.o bc/bc.o sem/coef.o math/math.o config/num_types.o bc/bcknd/device/device_facet_normal.o -bc/symmetry.o : bc/symmetry.f90 sem/coef.o adt/tuple.o adt/stack.o common/utils.o math/math.o bc/bc.o bc/dirichlet.o config/num_types.o config/neko_config.o bc/bcknd/device/device_symmetry.o +bc/symmetry.o : bc/symmetry.f90 bc/zero_dirichlet.o sem/coef.o adt/tuple.o adt/stack.o common/utils.o math/math.o bc/bc.o bc/dirichlet.o config/num_types.o config/neko_config.o bc/bcknd/device/device_symmetry.o bc/non_normal.o : bc/non_normal.f90 adt/stack.o common/utils.o math/math.o sem/coef.o device/device.o adt/tuple.o bc/dirichlet.o config/num_types.o config/neko_config.o bc/symmetry.o -bc/blasius.o : bc/blasius.f90 bc/bc.o fluid/flow_profile.o bc/bcknd/device/device_inhom_dirichlet.o device/device.o common/utils.o sem/coef.o config/num_types.o -bc/shear_stress.o : bc/shear_stress.f90 math/bcknd/device/device_math.o math/math.o bc/neumann.o bc/dirichlet.o sem/coef.o common/utils.o bc/bc.o config/num_types.o +bc/blasius.o : bc/blasius.f90 common/json_utils.o bc/bc.o fluid/flow_profile.o bc/bcknd/device/device_inhom_dirichlet.o device/device.o common/utils.o sem/coef.o config/num_types.o +bc/shear_stress.o : bc/shear_stress.f90 common/json_utils.o math/bcknd/device/device_math.o math/math.o bc/neumann.o bc/dirichlet.o sem/coef.o common/utils.o bc/bc.o config/num_types.o bc/wall_model_bc.o : bc/wall_model_bc.f90 bc/shear_stress.o wall_models/spalding.o wall_models/rough_log_law.o wall_models/wall_model.o sem/coef.o common/utils.o bc/bc.o config/num_types.o krylov/precon.o : krylov/precon.f90 config/num_types.o -krylov/krylov.o : krylov/krylov.f90 config/neko_config.o krylov/bcknd/device/pc_identity_device.o krylov/pc_identity.o bc/bc.o common/utils.o field/field.o mesh/mesh.o sem/coef.o krylov/precon.o config/num_types.o math/ax.o gs/gather_scatter.o +krylov/krylov.o : krylov/krylov.f90 config/neko_config.o krylov/bcknd/device/pc_identity_device.o krylov/pc_identity.o bc/bc_list.o common/utils.o field/field.o mesh/mesh.o sem/coef.o krylov/precon.o config/num_types.o math/ax.o gs/gather_scatter.o krylov/pc_identity.o : krylov/pc_identity.f90 config/num_types.o krylov/precon.o math/math.o krylov/precon_fctry.o : krylov/precon_fctry.f90 config/neko_config.o common/utils.o krylov/pc_hsmg.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/cpu/pc_jacobi.o krylov/bcknd/device/pc_identity_device.o krylov/pc_identity.o krylov/precon.o krylov/krylov_fctry.o : krylov/krylov_fctry.f90 config/neko_config.o common/utils.o krylov/precon.o krylov/krylov.o config/num_types.o krylov/bcknd/device/gmres_device.o krylov/bcknd/sx/gmres_sx.o krylov/bcknd/cpu/gmres.o krylov/bcknd/cpu/bicgstab.o krylov/bcknd/device/fusedcg_cpld_device.o krylov/bcknd/device/fusedcg_device.o krylov/bcknd/device/pipecg_device.o krylov/bcknd/sx/pipecg_sx.o krylov/bcknd/cpu/pipecg.o krylov/bcknd/cpu/cacg.o krylov/bcknd/device/cg_device.o krylov/bcknd/sx/cg_sx.o krylov/bcknd/cpu/cg.o -krylov/bcknd/cpu/cg.o : krylov/bcknd/cpu/cg.f90 comm/comm.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o -krylov/bcknd/cpu/cacg.o : krylov/bcknd/cpu/cacg.f90 math/mxm_wrapper.o comm/comm.o common/utils.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o -krylov/bcknd/cpu/pipecg.o : krylov/bcknd/cpu/pipecg.f90 comm/comm.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o -krylov/bcknd/cpu/bicgstab.o : krylov/bcknd/cpu/bicgstab.f90 common/utils.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o -krylov/bcknd/cpu/gmres.o : krylov/bcknd/cpu/gmres.f90 comm/comm.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/cpu/cg.o : krylov/bcknd/cpu/cg.f90 comm/comm.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o +krylov/bcknd/cpu/cacg.o : krylov/bcknd/cpu/cacg.f90 math/mxm_wrapper.o comm/comm.o common/utils.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o +krylov/bcknd/cpu/pipecg.o : krylov/bcknd/cpu/pipecg.f90 comm/comm.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/cpu/bicgstab.o : krylov/bcknd/cpu/bicgstab.f90 common/utils.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o +krylov/bcknd/cpu/gmres.o : krylov/bcknd/cpu/gmres.f90 comm/comm.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o krylov/bcknd/cpu/pc_jacobi.o : krylov/bcknd/cpu/pc_jacobi.f90 gs/gather_scatter.o sem/dofmap.o config/num_types.o sem/coef.o krylov/precon.o math/math.o -krylov/bcknd/sx/cg_sx.o : krylov/bcknd/sx/cg_sx.f90 math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o -krylov/bcknd/sx/pipecg_sx.o : krylov/bcknd/sx/pipecg_sx.f90 comm/comm.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o -krylov/bcknd/sx/gmres_sx.o : krylov/bcknd/sx/gmres_sx.f90 comm/comm.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/sx/cg_sx.o : krylov/bcknd/sx/cg_sx.f90 math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o +krylov/bcknd/sx/pipecg_sx.o : krylov/bcknd/sx/pipecg_sx.f90 comm/comm.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/sx/gmres_sx.o : krylov/bcknd/sx/gmres_sx.f90 comm/comm.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o krylov/bcknd/sx/pc_jacobi_sx.o : krylov/bcknd/sx/pc_jacobi_sx.f90 gs/gather_scatter.o config/num_types.o sem/dofmap.o sem/coef.o krylov/precon.o math/math.o -krylov/bcknd/device/cg_device.o : krylov/bcknd/device/cg_device.f90 math/bcknd/device/device_math.o device/device.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o -krylov/bcknd/device/pipecg_device.o : krylov/bcknd/device/pipecg_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o -krylov/bcknd/device/fusedcg_device.o : krylov/bcknd/device/fusedcg_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o -krylov/bcknd/device/fusedcg_cpld_device.o : krylov/bcknd/device/fusedcg_cpld_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o -krylov/bcknd/device/gmres_device.o : krylov/bcknd/device/gmres_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o krylov/bcknd/device/pc_identity_device.o bc/bc.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/device/cg_device.o : krylov/bcknd/device/cg_device.f90 math/bcknd/device/device_math.o device/device.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o math/ax.o krylov/precon.o krylov/krylov.o config/num_types.o +krylov/bcknd/device/pipecg_device.o : krylov/bcknd/device/pipecg_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/device/fusedcg_device.o : krylov/bcknd/device/fusedcg_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/device/fusedcg_cpld_device.o : krylov/bcknd/device/fusedcg_cpld_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o +krylov/bcknd/device/gmres_device.o : krylov/bcknd/device/gmres_device.F90 comm/comm.o device/device.o math/bcknd/device/device_math.o math/math.o krylov/bcknd/device/pc_identity_device.o bc/bc_list.o gs/gather_scatter.o sem/coef.o field/field.o config/num_types.o math/ax.o krylov/precon.o krylov/krylov.o krylov/bcknd/device/pc_jacobi_device.o : krylov/bcknd/device/pc_jacobi_device.F90 gs/gather_scatter.o device/device.o math/bcknd/device/device_math.o config/num_types.o sem/dofmap.o sem/coef.o krylov/precon.o krylov/bcknd/device/pc_identity_device.o : krylov/bcknd/device/pc_identity_device.f90 config/num_types.o krylov/precon.o math/bcknd/device/device_math.o device/device.o time_schemes/time_scheme.o : time_schemes/time_scheme.f90 config/num_types.o config/neko_config.o @@ -182,12 +184,12 @@ common/bcknd/device/rhs_maker_device.o : common/bcknd/device/rhs_maker_device.F9 common/material_properties.o : common/material_properties.f90 comm/comm.o common/utils.o common/user_intf.o common/log.o common/json_utils.o config/num_types.o config/neko_config.o : config/neko_config.f90 case.o : case.f90 common/material_properties.o mesh/point_zone_registry.o field/scratch_registry.o common/json_utils.o scalar/scalar_pnpn.o common/user_intf.o common/jobctrl.o common/log.o time_schemes/time_scheme_controller.o comm/comm.o mesh/mesh.o common/utils.o io/file.o common/statistics.o scalar/scalar_ic.o fluid/flow_ic.o common/sampler.o comm/redist.o comm/parmetis.o field/mesh_field.o io/fluid_stats_output.o io/mean_flow_output.o io/mean_sqr_flow_output.o io/chkp_output.o io/fluid_output.o fluid/fluid_scheme.o fluid/fluid_pnpn.o fluid/fluid_fctry.o config/num_types.o -common/user_intf.o : common/user_intf.f90 common/log.o common/utils.o common/json_utils.o config/num_types.o bc/field_dirichlet.o bc/usr_scalar.o bc/usr_inflow.o mesh/mesh.o bc/bc.o sem/coef.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o field/field.o -fluid/fluid_scheme.o : fluid/fluid_scheme.f90 common/time_step_controller.o field/field_series.o common/material_properties.o common/utils.o common/user_intf.o field/scratch_registry.o common/json_utils.o field/field_registry.o common/log.o math/operators.o math/mathops.o time_schemes/time_scheme_controller.o math/math.o mesh/mesh.o bc/bc.o fluid/fluid_stats.o krylov/precon_fctry.o krylov/krylov_fctry.o krylov/precon.o krylov/pc_hsmg.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/cpu/pc_jacobi.o bc/field_dirichlet_vector.o bc/field_dirichlet.o bc/non_normal.o bc/symmetry.o bc/dong_outflow.o bc/dirichlet.o bc/blasius.o bc/usr_inflow.o bc/inflow.o bc/wall.o sem/coef.o krylov/krylov.o sem/dofmap.o sem/space.o field/field.o field/field_list.o fluid/fluid_source_term.o fluid/fluid_user_source_term.o comm/comm.o config/num_types.o fluid/mean_flow.o common/checkpoint.o config/neko_config.o fluid/mean_sqr_flow.o gs/gather_scatter.o +common/user_intf.o : common/user_intf.f90 common/log.o common/utils.o common/json_utils.o config/num_types.o bc/field_dirichlet.o bc/usr_scalar.o bc/usr_inflow.o mesh/mesh.o bc/bc_list.o sem/coef.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o field/field.o +fluid/fluid_scheme.o : fluid/fluid_scheme.f90 bc/bc_list.o common/time_step_controller.o field/field_series.o common/material_properties.o common/utils.o common/user_intf.o field/scratch_registry.o common/json_utils.o field/field_registry.o common/log.o math/operators.o math/mathops.o time_schemes/time_scheme_controller.o math/math.o mesh/mesh.o bc/bc.o fluid/fluid_stats.o krylov/precon_fctry.o krylov/krylov_fctry.o krylov/precon.o krylov/pc_hsmg.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/cpu/pc_jacobi.o bc/field_dirichlet_vector.o bc/field_dirichlet.o bc/non_normal.o bc/symmetry.o bc/dong_outflow.o bc/dirichlet.o bc/blasius.o bc/usr_inflow.o bc/inflow.o sem/coef.o bc/zero_dirichlet.o krylov/krylov.o sem/dofmap.o sem/space.o field/field.o field/field_list.o fluid/fluid_source_term.o fluid/fluid_user_source_term.o comm/comm.o config/num_types.o fluid/mean_flow.o common/checkpoint.o config/neko_config.o fluid/mean_sqr_flow.o gs/gather_scatter.o fluid/fluid_aux.o : fluid/fluid_aux.f90 krylov/krylov.o config/num_types.o common/log.o -fluid/fluid_pnpn.o : fluid/fluid_pnpn.f90 bc/bc.o math/mathops.o math/math.o config/neko_config.o gs/gather_scatter.o common/time_step_controller.o sem/coef.o common/user_intf.o mesh/mesh.o bc/non_normal.o bc/facet_normal.o bc/dirichlet.o field/field.o math/ax.o fluid/advection_fctry.o common/material_properties.o common/json_utils.o common/profiler.o fluid/advection.o common/log.o device/device.o common/projection.o time_schemes/time_scheme_controller.o fluid/fluid_aux.o math/bcknd/device/device_mathops.o math/bcknd/device/device_math.o field/field_series.o fluid/fluid_scheme.o fluid/fluid_volflow.o common/rhs_maker.o common/rhs_maker_fctry.o math/ax_helm_fctry.o fluid/pnpn_res.o fluid/pnpn_res_fctry.o krylov/krylov.o config/num_types.o +fluid/fluid_pnpn.o : fluid/fluid_pnpn.f90 bc/zero_dirichlet.o bc/bc_list.o math/mathops.o math/math.o config/neko_config.o gs/gather_scatter.o common/time_step_controller.o sem/coef.o common/user_intf.o mesh/mesh.o bc/non_normal.o bc/facet_normal.o bc/dirichlet.o field/field.o math/ax.o fluid/advection_fctry.o common/material_properties.o common/json_utils.o common/profiler.o fluid/advection.o common/log.o device/device.o common/projection.o time_schemes/time_scheme_controller.o fluid/fluid_aux.o math/bcknd/device/device_mathops.o math/bcknd/device/device_math.o field/field_series.o fluid/fluid_scheme.o fluid/fluid_volflow.o common/rhs_maker.o common/rhs_maker_fctry.o math/ax_helm_fctry.o fluid/pnpn_res.o fluid/pnpn_res_fctry.o krylov/krylov.o config/num_types.o fluid/fluid_fctry.o : fluid/fluid_fctry.f90 common/utils.o fluid/fluid_pnpn.o fluid/fluid_scheme.o -fluid/fluid_volflow.o : fluid/fluid_volflow.f90 math/ax.o bc/bc.o field/scratch_registry.o common/json_utils.o gs/gather_scatter.o math/bcknd/device/device_mathops.o math/bcknd/device/device_math.o config/neko_config.o comm/comm.o math/math.o time_schemes/time_scheme_controller.o sem/coef.o field/field.o sem/dofmap.o krylov/precon.o krylov/krylov.o math/mathops.o config/num_types.o math/operators.o +fluid/fluid_volflow.o : fluid/fluid_volflow.f90 math/ax.o bc/bc_list.o field/scratch_registry.o common/json_utils.o gs/gather_scatter.o math/bcknd/device/device_mathops.o math/bcknd/device/device_math.o config/neko_config.o comm/comm.o math/math.o time_schemes/time_scheme_controller.o sem/coef.o field/field.o sem/dofmap.o krylov/precon.o krylov/krylov.o math/mathops.o config/num_types.o math/operators.o fluid/pnpn_res.o : fluid/pnpn_res.f90 config/num_types.o mesh/mesh.o sem/space.o bc/facet_normal.o sem/coef.o field/field.o math/ax.o gs/gather_scatter.o fluid/pnpn_res_fctry.o : fluid/pnpn_res_fctry.f90 fluid/bcknd/sx/pnpn_res_sx.o fluid/bcknd/cpu/pnpn_res_cpu.o fluid/bcknd/device/pnpn_res_device.o fluid/pnpn_res.o common/utils.o config/neko_config.o fluid/mean_flow.o : fluid/mean_flow.f90 field/field.o field/mean_field.o @@ -218,11 +220,11 @@ math/bcknd/cpu/ax_helm_full_cpu.o : math/bcknd/cpu/ax_helm_full_cpu.f90 common/u math/bcknd/sx/ax_helm_sx.o : math/bcknd/sx/ax_helm_sx.f90 math/math.o mesh/mesh.o sem/space.o sem/coef.o config/num_types.o math/ax_helm.o math/bcknd/xsmm/ax_helm_xsmm.o : math/bcknd/xsmm/ax_helm_xsmm.F90 math/mxm_wrapper.o mesh/mesh.o sem/space.o sem/coef.o config/num_types.o math/ax_helm.o math/bcknd/device/ax_helm_device.o : math/bcknd/device/ax_helm_device.F90 device/device.o math/bcknd/device/device_math.o mesh/mesh.o sem/space.o sem/coef.o config/num_types.o math/ax_helm.o -common/projection.o : common/projection.f90 common/time_step_controller.o common/log.o common/profiler.o common/bcknd/device/device_projection.o math/bcknd/device/device_math.o device/device.o config/neko_config.o gs/gather_scatter.o comm/comm.o bc/bc.o math/ax.o sem/coef.o math/math.o config/num_types.o +common/projection.o : common/projection.f90 common/time_step_controller.o common/log.o common/profiler.o common/bcknd/device/device_projection.o math/bcknd/device/device_math.o device/device.o config/neko_config.o gs/gather_scatter.o comm/comm.o bc/bc_list.o math/ax.o sem/coef.o math/math.o config/num_types.o common/bcknd/device/device_projection.o : common/bcknd/device/device_projection.F90 common/utils.o config/num_types.o comm/parmetis.o : comm/parmetis.F90 mesh/mesh.o field/mesh_field.o config/num_types.o common/utils.o mesh/point.o comm/comm.o comm/redist.o : comm/redist.f90 mesh/element.o mesh/facet_zone.o io/format/nmsh.o mesh/mesh.o comm/comm.o mesh/curve.o adt/stack.o mesh/point.o adt/htable.o comm/mpi_types.o field/mesh_field.o -krylov/pc_hsmg.o : krylov/pc_hsmg.f90 krylov/krylov_fctry.o krylov/krylov.o mesh/mesh.o sem/coef.o field/field.o sem/dofmap.o sem/space.o common/profiler.o math/bcknd/device/device_math.o device/device.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/cpu/pc_jacobi.o math/schwarz.o bc/dirichlet.o bc/bc.o sem/interpolation.o gs/gather_scatter.o math/ax_helm_fctry.o math/ax.o krylov/precon.o common/utils.o math/math.o config/num_types.o config/neko_config.o +krylov/pc_hsmg.o : krylov/pc_hsmg.f90 bc/bc_list.o bc/zero_dirichlet.o krylov/krylov_fctry.o krylov/krylov.o mesh/mesh.o sem/coef.o field/field.o sem/dofmap.o sem/space.o common/profiler.o math/bcknd/device/device_math.o device/device.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/cpu/pc_jacobi.o math/schwarz.o bc/dirichlet.o bc/bc.o sem/interpolation.o gs/gather_scatter.o math/ax_helm_fctry.o math/ax.o krylov/precon.o common/utils.o math/math.o config/num_types.o config/neko_config.o common/signal.o : common/signal.f90 common/utils.o common/jobctrl.o : common/jobctrl.f90 common/log.o comm/comm.o common/utils.o common/signal.o config/num_types.o device/cuda_intf.o : device/cuda_intf.F90 common/utils.o @@ -241,14 +243,14 @@ math/bcknd/device/fdm_device.o : math/bcknd/device/fdm_device.F90 device/device. math/bcknd/device/device_mathops.o : math/bcknd/device/device_mathops.F90 common/utils.o config/num_types.o bc/bcknd/device/device_dirichlet.o : bc/bcknd/device/device_dirichlet.F90 common/utils.o config/num_types.o bc/bcknd/device/device_inflow.o : bc/bcknd/device/device_inflow.F90 common/utils.o config/num_types.o -bc/bcknd/device/device_wall.o : bc/bcknd/device/device_wall.F90 common/utils.o +bc/bcknd/device/device_zero_dirichlet.o : bc/bcknd/device/device_zero_dirichlet.F90 common/utils.o bc/bcknd/device/device_symmetry.o : bc/bcknd/device/device_symmetry.F90 common/utils.o config/num_types.o bc/bcknd/device/device_facet_normal.o : bc/bcknd/device/device_facet_normal.F90 common/utils.o bc/bcknd/device/device_inhom_dirichlet.o : bc/bcknd/device/device_inhom_dirichlet.F90 common/utils.o config/num_types.o bc/bcknd/device/device_dong_outflow.o : bc/bcknd/device/device_dong_outflow.F90 common/utils.o config/num_types.o scalar/bcknd/device/scalar_residual_device.o : scalar/bcknd/device/scalar_residual_device.F90 math/bcknd/device/device_mathops.o math/bcknd/device/device_math.o math/operators.o gs/gather_scatter.o scalar/scalar_residual.o -scalar/scalar_scheme.o : scalar/scalar_scheme.f90 common/time_step_controller.o field/field_series.o config/neko_config.o math/bcknd/device/device_math.o math/math.o scalar/scalar_source_term.o comm/comm.o common/utils.o common/material_properties.o common/user_intf.o common/json_utils.o bc/usr_scalar.o field/field_registry.o common/log.o time_schemes/time_scheme_controller.o mesh/facet_zone.o mesh/mesh.o bc/field_dirichlet.o bc/bc.o krylov/precon.o krylov/precon_fctry.o krylov/pc_hsmg.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/cpu/pc_jacobi.o krylov/krylov_fctry.o bc/neumann.o bc/dirichlet.o sem/coef.o krylov/krylov.o sem/dofmap.o sem/space.o field/field_list.o field/field.o config/num_types.o common/checkpoint.o gs/gather_scatter.o -scalar/scalar_pnpn.o : scalar/scalar_pnpn.f90 field/scratch_registry.o common/time_step_controller.o config/neko_config.o common/material_properties.o common/user_intf.o common/json_utils.o common/profiler.o fluid/advection_fctry.o fluid/advection.o common/log.o math/math.o common/projection.o time_schemes/time_scheme_controller.o scalar/scalar_aux.o math/bcknd/device/device_math.o krylov/krylov.o bc/facet_normal.o field/field_series.o math/ax.o scalar/scalar_residual.o gs/gather_scatter.o device/device.o sem/coef.o common/checkpoint.o mesh/mesh.o bc/bc.o field/field.o bc/neumann.o bc/dirichlet.o scalar/scalar_scheme.o common/rhs_maker.o common/rhs_maker_fctry.o math/ax_helm_fctry.o scalar/scalar_residual_fctry.o config/num_types.o +scalar/scalar_scheme.o : scalar/scalar_scheme.f90 common/time_step_controller.o config/neko_config.o math/bcknd/device/device_math.o math/math.o bc/bc_fctry.o field/field_series.o scalar/scalar_source_term.o comm/comm.o common/utils.o common/material_properties.o common/user_intf.o common/json_utils.o bc/usr_scalar.o field/field_registry.o common/log.o time_schemes/time_scheme_controller.o mesh/facet_zone.o mesh/mesh.o bc/field_dirichlet.o krylov/precon.o krylov/precon_fctry.o bc/bc_list.o bc/bc.o krylov/pc_hsmg.o krylov/bcknd/sx/pc_jacobi_sx.o krylov/bcknd/device/pc_jacobi_device.o krylov/bcknd/cpu/pc_jacobi.o krylov/krylov_fctry.o bc/neumann.o bc/dirichlet.o sem/coef.o krylov/krylov.o sem/dofmap.o sem/space.o field/field_list.o field/field.o config/num_types.o common/checkpoint.o gs/gather_scatter.o +scalar/scalar_pnpn.o : scalar/scalar_pnpn.f90 field/scratch_registry.o common/time_step_controller.o bc/zero_dirichlet.o config/neko_config.o common/material_properties.o common/user_intf.o common/json_utils.o common/profiler.o fluid/advection_fctry.o fluid/advection.o common/log.o math/math.o common/projection.o time_schemes/time_scheme_controller.o scalar/scalar_aux.o math/bcknd/device/device_math.o krylov/krylov.o bc/facet_normal.o field/field_series.o math/ax.o scalar/scalar_residual.o gs/gather_scatter.o device/device.o sem/coef.o common/checkpoint.o mesh/mesh.o bc/bc_list.o field/field.o bc/neumann.o bc/dirichlet.o scalar/scalar_scheme.o common/rhs_maker.o common/rhs_maker_fctry.o math/ax_helm_fctry.o scalar/scalar_residual_fctry.o config/num_types.o scalar/scalar_aux.o : scalar/scalar_aux.f90 krylov/krylov.o config/num_types.o common/log.o scalar/scalar_residual.o : scalar/scalar_residual.f90 config/num_types.o mesh/mesh.o sem/space.o bc/facet_normal.o scalar/source_scalar.o sem/coef.o field/field.o math/ax.o gs/gather_scatter.o scalar/scalar_residual_fctry.o : scalar/scalar_residual_fctry.f90 scalar/bcknd/sx/scalar_residual_sx.o scalar/bcknd/cpu/scalar_residual_cpu.o scalar/bcknd/device/scalar_residual_device.o scalar/scalar_residual.o config/neko_config.o @@ -293,5 +295,5 @@ wall_models/rough_log_law.o : wall_models/rough_log_law.f90 common/utils.o commo wall_models/spalding.o : wall_models/spalding.f90 common/utils.o common/log.o common/json_utils.o field/field_registry.o wall_models/wall_model.o config/neko_config.o sem/coef.o sem/dofmap.o config/num_types.o field/field.o wall_models/wall_model.o : wall_models/wall_model.f90 common/log.o comm/comm.o math/math.o common/utils.o math/vector.o device/device.o config/neko_config.o sem/coef.o sem/dofmap.o field/field_registry.o field/field.o config/num_types.o wall_models/wall_model_fctry.o : wall_models/wall_model_fctry.f90 common/json_utils.o common/utils.o wall_models/rough_log_law.o wall_models/spalding.o sem/coef.o sem/dofmap.o les/vreman.o wall_models/wall_model.o config/num_types.o -neko.o : neko.f90 bc/field_dirichlet_vector.o bc/field_dirichlet.o mesh/point_zone_registry.o mesh/point_zones/sphere_point_zone.o mesh/point_zones/box_point_zone.o mesh/point_zone.o sem/point_interpolator.o common/time_interpolator.o io/data_streamer.o simulation_components/simcomp_executor.o field/scratch_registry.o field/field_registry.o qoi/drag_torque.o common/system.o sem/spectral_error_indicator.o simulation_components/probes.o simulation_components/simulation_component.o math/tensor.o math/matrix.o math/vector.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o fluid/fluid_stats.o sem/cpr.o sem/map_1d.o math/bcknd/device/device_math.o device/device.o common/jobctrl.o common/signal.o comm/parmetis.o common/user_intf.o common/projection.o math/mathops.o math/operators.o simulation.o io/output.o common/sampler.o case.o config/neko_config.o math/ax.o math/ax_helm_fctry.o krylov/precon_fctry.o krylov/krylov_fctry.o bc/dirichlet.o bc/wall.o bc/bc.o sem/coef.o gs/gather_scatter.o comm/mpi_types.o field/field.o io/file.o common/global_interpolation.o math/mxm_wrapper.o io/format/map.o field/mesh_field.o mesh/point.o mesh/mesh.o adt/tuple.o adt/stack.o adt/uset.o adt/htable.o sem/space.o sem/dofmap.o sem/speclib.o math/math.o common/log.o common/utils.o comm/comm.o config/num_types.o +neko.o : neko.f90 bc/field_dirichlet_vector.o bc/field_dirichlet.o mesh/point_zone_registry.o mesh/point_zones/sphere_point_zone.o mesh/point_zones/box_point_zone.o mesh/point_zone.o sem/point_interpolator.o common/time_interpolator.o io/data_streamer.o simulation_components/simcomp_executor.o field/scratch_registry.o field/field_registry.o qoi/drag_torque.o common/system.o sem/spectral_error_indicator.o simulation_components/probes.o simulation_components/simulation_component.o math/tensor.o math/matrix.o math/vector.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o fluid/fluid_stats.o sem/cpr.o sem/map_1d.o math/bcknd/device/device_math.o device/device.o common/jobctrl.o common/signal.o comm/parmetis.o common/user_intf.o common/projection.o math/mathops.o math/operators.o simulation.o io/output.o common/sampler.o case.o config/neko_config.o math/ax.o math/ax_helm_fctry.o krylov/precon_fctry.o krylov/krylov_fctry.o bc/dirichlet.o bc/zero_dirichlet.o bc/bc.o sem/coef.o gs/gather_scatter.o comm/mpi_types.o field/field.o io/file.o common/global_interpolation.o math/mxm_wrapper.o io/format/map.o field/mesh_field.o mesh/point.o mesh/mesh.o adt/tuple.o adt/stack.o adt/uset.o adt/htable.o sem/space.o sem/dofmap.o sem/speclib.o math/math.o common/log.o common/utils.o comm/comm.o config/num_types.o driver.o : driver.f90 neko.o diff --git a/src/.depends_device b/src/.depends_device index efef5d7c4cb..987c611b531 100644 --- a/src/.depends_device +++ b/src/.depends_device @@ -18,14 +18,14 @@ krylov/bcknd/device/hip/gmres_aux.o : krylov/bcknd/device/hip/gmres_aux.hip kryl gs/bcknd/device/hip/gs.o : gs/bcknd/device/hip/gs.hip gs/bcknd/device/hip/gs_kernels.h bc/bcknd/device/hip/dirichlet.o : bc/bcknd/device/hip/dirichlet.hip bc/bcknd/device/hip/dirichlet_kernel.h bc/bcknd/device/hip/inflow.o : bc/bcknd/device/hip/inflow.hip bc/bcknd/device/hip/inflow_kernel.h -bc/bcknd/device/hip/no_slip_wall.o : bc/bcknd/device/hip/no_slip_wall.hip bc/bcknd/device/hip/no_slip_wall_kernel.h +bc/bcknd/device/hip/zero_dirichlet.o : bc/bcknd/device/hip/zero_dirichlet.hip bc/bcknd/device/hip/zero_dirichlet_kernel.h bc/bcknd/device/hip/symmetry.o : bc/bcknd/device/hip/symmetry.hip bc/bcknd/device/hip/symmetry_kernel.h bc/bcknd/device/hip/facet_normal.o : bc/bcknd/device/hip/facet_normal.hip bc/bcknd/device/hip/facet_normal_kernel.h bc/bcknd/device/hip/inhom_dirichlet.o : bc/bcknd/device/hip/inhom_dirichlet.hip bc/bcknd/device/hip/inhom_dirichlet_kernel.h bc/bcknd/device/hip/dong_outflow.o : bc/bcknd/device/hip/dong_outflow.hip bc/bcknd/device/hip/dong_outflow_kernel.h common/bcknd/device/hip/rhs_maker.o : common/bcknd/device/hip/rhs_maker.hip common/bcknd/device/hip/sumab_kernel.h common/bcknd/device/hip/makeext_kernel.h common/bcknd/device/hip/makebdf_kernel.h common/bcknd/device/hip/projection.o : common/bcknd/device/hip/projection.hip common/bcknd/device/hip/projection_kernel.h -fluid/bcknd/device/hip/pnpn_res.o : fluid/bcknd/device/hip/pnpn_res.hip fluid/bcknd/device/hip/vel_res_update_kernel.h fluid/bcknd/device/hip/prs_res_kernel.h +fluid/bcknd/device/hip/pnpn_res.o : fluid/bcknd/device/hip/pnpn_res.hip fluid/bcknd/device/hip/vel_res_update_kernel.h fluid/bcknd/device/hip/prs_res_kernel.h scalar/bcknd/device/hip/scalar_residual.o : scalar/bcknd/device/hip/scalar_residual.hip scalar/bcknd/device/hip/scalar_residual_update_kernel.h sem/bcknd/device/hip/coef.o : sem/bcknd/device/hip/coef.hip sem/bcknd/device/hip/coef_kernel.h device/cuda/check.o : device/cuda/check.cu device/cuda/check.h @@ -46,22 +46,22 @@ krylov/bcknd/device/cuda/pipecg_aux.o : krylov/bcknd/device/cuda/pipecg_aux.cu k krylov/bcknd/device/cuda/fusedcg_aux.o : krylov/bcknd/device/cuda/fusedcg_aux.cu krylov/bcknd/device/cuda/fusedcg_kernel.h krylov/bcknd/device/cuda/fusedcg_cpld_aux.o : krylov/bcknd/device/cuda/fusedcg_cpld_aux.cu krylov/bcknd/device/cuda/fusedcg_cpld_kernel.h krylov/bcknd/device/cuda/gmres_aux.o : krylov/bcknd/device/cuda/gmres_aux.cu krylov/bcknd/device/cuda/gmres_kernel.h -gs/bcknd/device/cuda/gs.o : gs/bcknd/device/cuda/gs.cu gs/bcknd/device/cuda/gs_kernels.h +gs/bcknd/device/cuda/gs.o : gs/bcknd/device/cuda/gs.cu gs/bcknd/device/cuda/gs_kernels.h bc/bcknd/device/cuda/dirichlet.o : bc/bcknd/device/cuda/dirichlet.cu bc/bcknd/device/cuda/dirichlet_kernel.h bc/bcknd/device/cuda/inflow.o : bc/bcknd/device/cuda/inflow.cu bc/bcknd/device/cuda/inflow_kernel.h -bc/bcknd/device/cuda/no_slip_wall.o : bc/bcknd/device/cuda/no_slip_wall.cu bc/bcknd/device/cuda/no_slip_wall_kernel.h +bc/bcknd/device/cuda/zero_dirichlet.o : bc/bcknd/device/cuda/zero_dirichlet.cu bc/bcknd/device/cuda/zero_dirichlet_kernel.h bc/bcknd/device/cuda/symmetry.o : bc/bcknd/device/cuda/symmetry.cu bc/bcknd/device/cuda/symmetry_kernel.h bc/bcknd/device/cuda/facet_normal.o : bc/bcknd/device/cuda/facet_normal.cu bc/bcknd/device/cuda/facet_normal_kernel.h bc/bcknd/device/cuda/inhom_dirichlet.o : bc/bcknd/device/cuda/inhom_dirichlet.cu bc/bcknd/device/cuda/inhom_dirichlet_kernel.h bc/bcknd/device/cuda/dong_outflow.o : bc/bcknd/device/cuda/dong_outflow.cu bc/bcknd/device/cuda/dong_outflow_kernel.h common/bcknd/device/cuda/rhs_maker.o : common/bcknd/device/cuda/rhs_maker.cu common/bcknd/device/cuda/sumab_kernel.h common/bcknd/device/cuda/makeext_kernel.h common/bcknd/device/cuda/makebdf_kernel.h common/bcknd/device/cuda/projection.o : common/bcknd/device/cuda/projection.cu common/bcknd/device/cuda/projection_kernel.h -fluid/bcknd/device/cuda/pnpn_res.o : fluid/bcknd/device/cuda/pnpn_res.cu fluid/bcknd/device/cuda/vel_res_update_kernel.h fluid/bcknd/device/cuda/prs_res_kernel.h +fluid/bcknd/device/cuda/pnpn_res.o : fluid/bcknd/device/cuda/pnpn_res.cu fluid/bcknd/device/cuda/vel_res_update_kernel.h fluid/bcknd/device/cuda/prs_res_kernel.h scalar/bcknd/device/cuda/scalar_residual.o : scalar/bcknd/device/cuda/scalar_residual.cu scalar/bcknd/device/cuda/scalar_residual_update_kernel.h sem/bcknd/device/cuda/coef.o : sem/bcknd/device/cuda/coef.cu sem/bcknd/device/cuda/coef_kernel.h device/opencl/check.o : device/opencl/check.c device/opencl/check.h device/opencl/jit.o : device/opencl/jit.c device/opencl/jit.h -math/bcknd/device/opencl/math.o : math/bcknd/device/opencl/math.c math/bcknd/device/opencl/math_kernel.cl.h +math/bcknd/device/opencl/math.o : math/bcknd/device/opencl/math.c math/bcknd/device/opencl/math_kernel.cl.h math/bcknd/device/opencl/schwarz.o : math/bcknd/device/opencl/schwarz.c math/bcknd/device/opencl/schwarz_kernel.cl.h math/bcknd/device/opencl/tensor.o : math/bcknd/device/opencl/tensor.c math/bcknd/device/opencl/tensor_kernel.cl.h math/bcknd/device/opencl/fdm.o : math/bcknd/device/opencl/fdm.c math/bcknd/device/opencl/fdm_kernel.cl.h @@ -74,10 +74,10 @@ math/bcknd/device/opencl/opr_lambda2.o : math/bcknd/device/opencl/opr_lambda2.c math/bcknd/device/opencl/opr_cfl.o : math/bcknd/device/opencl/opr_cfl.c math/bcknd/device/opencl/cfl_kernel.cl.h math/bcknd/device/opencl/ax_helm.o : math/bcknd/device/opencl/ax_helm.c math/bcknd/device/opencl/ax_helm_kernel.cl.h krylov/bcknd/device/opencl/pc_jacobi.o : krylov/bcknd/device/opencl/pc_jacobi.c krylov/bcknd/device/opencl/jacobi_kernel.cl.h -gs/bcknd/device/opencl/gs.o : gs/bcknd/device/opencl/gs.c gs/bcknd/device/opencl/gs_kernels.cl.h +gs/bcknd/device/opencl/gs.o : gs/bcknd/device/opencl/gs.c gs/bcknd/device/opencl/gs_kernels.cl.h bc/bcknd/device/opencl/dirichlet.o : bc/bcknd/device/opencl/dirichlet.c bc/bcknd/device/opencl/dirichlet_kernel.cl.h bc/bcknd/device/opencl/inflow.o : bc/bcknd/device/opencl/inflow.c bc/bcknd/device/opencl/inflow_kernel.cl.h -bc/bcknd/device/opencl/no_slip_wall.o : bc/bcknd/device/opencl/no_slip_wall.c bc/bcknd/device/opencl/no_slip_wall_kernel.cl.h +bc/bcknd/device/opencl/zero_dirichlet.o : bc/bcknd/device/opencl/zero_dirichlet.c bc/bcknd/device/opencl/zero_dirichlet_kernel.cl.h bc/bcknd/device/opencl/symmetry.o : bc/bcknd/device/opencl/symmetry.c bc/bcknd/device/opencl/symmetry_kernel.cl.h bc/bcknd/device/opencl/facet_normal.o : bc/bcknd/device/opencl/facet_normal.c bc/bcknd/device/opencl/facet_normal_kernel.cl.h bc/bcknd/device/opencl/inhom_dirichlet.o : bc/bcknd/device/opencl/inhom_dirichlet.c bc/bcknd/device/opencl/inhom_dirichlet_kernel.cl.h diff --git a/src/Makefile.am b/src/Makefile.am index 9ea00d73834..7cf73f5ae22 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -130,10 +130,12 @@ neko_fortran_SOURCES = \ common/profiler.F90\ common/craypat.F90\ bc/bc.f90\ + bc/bc_list.f90\ + bc/bc_fctry.f90\ bc/dirichlet.f90\ bc/neumann.f90\ bc/dong_outflow.f90\ - bc/wall.f90\ + bc/zero_dirichlet.f90\ bc/inflow.f90\ bc/field_dirichlet.f90\ bc/field_dirichlet_vector.f90\ @@ -244,7 +246,7 @@ neko_fortran_SOURCES = \ math/bcknd/device/device_mathops.F90\ bc/bcknd/device/device_dirichlet.F90\ bc/bcknd/device/device_inflow.F90\ - bc/bcknd/device/device_wall.F90\ + bc/bcknd/device/device_zero_dirichlet.F90\ bc/bcknd/device/device_symmetry.F90\ bc/bcknd/device/device_facet_normal.F90\ bc/bcknd/device/device_inhom_dirichlet.F90\ @@ -335,7 +337,7 @@ libneko_a_SOURCES += \ gs/bcknd/device/hip/gs.hip\ bc/bcknd/device/hip/dirichlet.hip\ bc/bcknd/device/hip/inflow.hip\ - bc/bcknd/device/hip/no_slip_wall.hip\ + bc/bcknd/device/hip/zero_dirichlet.hip\ bc/bcknd/device/hip/symmetry.hip\ bc/bcknd/device/hip/facet_normal.hip\ bc/bcknd/device/hip/inhom_dirichlet.hip\ @@ -370,7 +372,7 @@ libneko_a_SOURCES += \ gs/bcknd/device/cuda/gs.cu\ bc/bcknd/device/cuda/dirichlet.cu\ bc/bcknd/device/cuda/inflow.cu\ - bc/bcknd/device/cuda/no_slip_wall.cu\ + bc/bcknd/device/cuda/zero_dirichlet.cu\ bc/bcknd/device/cuda/symmetry.cu\ bc/bcknd/device/cuda/facet_normal.cu\ bc/bcknd/device/cuda/inhom_dirichlet.cu\ @@ -403,7 +405,7 @@ libneko_a_SOURCES += \ gs/bcknd/device/opencl/gs.c\ bc/bcknd/device/opencl/dirichlet.c\ bc/bcknd/device/opencl/inflow.c\ - bc/bcknd/device/opencl/no_slip_wall.c\ + bc/bcknd/device/opencl/zero_dirichlet.c\ bc/bcknd/device/opencl/symmetry.c\ bc/bcknd/device/opencl/facet_normal.c\ bc/bcknd/device/opencl/inhom_dirichlet.c\ @@ -491,7 +493,7 @@ endif EXTRA_DIST = \ bc/bcknd/device/cuda/dirichlet_kernel.h\ bc/bcknd/device/cuda/inflow_kernel.h\ - bc/bcknd/device/cuda/no_slip_wall_kernel.h\ + bc/bcknd/device/cuda/zero_dirichlet_kernel.h\ bc/bcknd/device/cuda/symmetry_kernel.h\ bc/bcknd/device/cuda/facet_normal_kernel.h\ bc/bcknd/device/cuda/inhom_dirichlet_kernel.h\ @@ -520,7 +522,7 @@ EXTRA_DIST = \ sem/bcknd/device/cuda/coef_kernel.h\ bc/bcknd/device/hip/dirichlet_kernel.h\ bc/bcknd/device/hip/inflow_kernel.h\ - bc/bcknd/device/hip/no_slip_wall_kernel.h\ + bc/bcknd/device/hip/zero_dirichlet_kernel.h\ bc/bcknd/device/hip/symmetry_kernel.h\ bc/bcknd/device/hip/facet_normal_kernel.h\ bc/bcknd/device/hip/inhom_dirichlet_kernel.h\ @@ -555,7 +557,7 @@ EXTRA_DIST = \ sem/bcknd/device/hip/coef_kernel.h\ bc/bcknd/device/opencl/dirichlet_kernel.cl\ bc/bcknd/device/opencl/inflow_kernel.cl\ - bc/bcknd/device/opencl/no_slip_wall_kernel.cl\ + bc/bcknd/device/opencl/zero_dirichlet_kernel.cl\ bc/bcknd/device/opencl/symmetry_kernel.cl\ bc/bcknd/device/opencl/facet_normal_kernel.cl\ bc/bcknd/device/opencl/inhom_dirichlet_kernel.cl\ diff --git a/src/bc/bc.f90 b/src/bc/bc.f90 index 1b4a11958eb..0a8fbc74a79 100644 --- a/src/bc/bc.f90 +++ b/src/bc/bc.f90 @@ -44,6 +44,7 @@ module bc use tuple, only : tuple_i4_t use utils, only : neko_error, linear_index, split_string use, intrinsic :: iso_c_binding, only : c_ptr, C_NULL_PTR + use json_module, only : json_file implicit none private @@ -67,6 +68,10 @@ module bc type(c_ptr) :: msk_d = C_NULL_PTR !> Device pointer for facet type(c_ptr) :: facet_d = C_NULL_PTR + !> Wether the bc is strongly enforced. Essentially valid for all Dirichlet + !! types of bcs. These need to be masked out for solvers etc, so that + !! values are not affected. + logical :: strong = .true. contains !> Constructor procedure, pass(this) :: init_base => bc_init_base @@ -82,32 +87,34 @@ module bc procedure, pass(this) :: mark_zone => bc_mark_zone !> Finalize the construction of the bc by populting the msk and facet !! arrays - procedure, pass(this) :: finalize => bc_finalize - !> Apply the boundary condition to a scalar field + procedure, pass(this) :: finalize_base => bc_finalize_base + + !> Apply the boundary condition to a scalar field. Dispatches to the CPU + !! or the device version. + procedure, pass(this) :: apply_scalar_generic => bc_apply_scalar_generic + !> Apply the boundary condition to a vector field. Dispatches to the CPU + !! or the device version. + procedure, pass(this) :: apply_vector_generic => bc_apply_vector_generic + !> Apply the boundary condition to a scalar field on the CPU. procedure(bc_apply_scalar), pass(this), deferred :: apply_scalar - !> Apply the boundary condition to a vector field + !> Apply the boundary condition to a vector field on the CPU. procedure(bc_apply_vector), pass(this), deferred :: apply_vector - !> Device version of \ref apply_scalar + !> Device version of \ref apply_scalar on the device. procedure(bc_apply_scalar_dev), pass(this), deferred :: apply_scalar_dev - !> Device version of \ref apply_vector + !> Device version of \ref apply_vector on the device. procedure(bc_apply_vector_dev), pass(this), deferred :: apply_vector_dev - !> Destructor + !> Deferred destructor. procedure(bc_destructor), pass(this), deferred :: free + !> Deferred constructor. + procedure(bc_constructor), pass(this), deferred :: init + !> Deferred finalizer. + procedure(bc_finalize), pass(this), deferred :: finalize end type bc_t - !> Pointer to boundary condtiion - type, private :: bcp_t - class(bc_t), pointer :: bcp - end type bcp_t - - !> A list of boundary conditions - type, public :: bc_list_t - type(bcp_t), allocatable :: bc(:) - !> Number of items. - integer :: n - !> Capacity. - integer :: size - end type bc_list_t + !> Pointer to a @ref `bc_t`. + type, public :: bc_ptr_t + class(bc_t), pointer :: ptr + end type bc_ptr_t abstract interface !> Apply the boundary condition to a scalar field @@ -147,6 +154,16 @@ subroutine bc_apply_vector(this, x, y, z, n, t, tstep) end subroutine bc_apply_vector end interface + abstract interface + !> Constructor + subroutine bc_constructor(this, coef, json) + import :: bc_t, coef_t, json_file + class(bc_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + end subroutine bc_constructor + end interface + abstract interface !> Destructor subroutine bc_destructor(this) @@ -155,6 +172,14 @@ subroutine bc_destructor(this) end subroutine bc_destructor end interface + abstract interface + !> Finalize by building the mask and facet arrays. + subroutine bc_finalize(this) + import :: bc_t + class(bc_t), intent(inout), target :: this + end subroutine bc_finalize + end interface + abstract interface !> Apply the boundary condition to a scalar field on the device !! @param x_d Device pointer to the field. @@ -187,13 +212,6 @@ subroutine bc_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) end subroutine bc_apply_vector_dev end interface - interface bc_list_apply - module procedure bc_list_apply_scalar, bc_list_apply_vector - end interface bc_list_apply - - public :: bc_list_init, bc_list_free, bc_list_add, & - bc_list_apply_scalar, bc_list_apply_vector, bc_list_apply - contains !> Constructor @@ -244,6 +262,98 @@ subroutine bc_free_base(this) end subroutine bc_free_base + !> Apply the boundary condition to a vector field. Dispatches to the CPU + !! or the device version. + !! @param x The x comp of the field for which to apply the bc. + !! @param y The y comp of the field for which to apply the bc. + !! @param z The z comp of the field for which to apply the bc. + !! @param n The size of x, y, and z. + !! @param t Current time. + !! @param tstep The current time iteration. + subroutine bc_apply_vector_generic(this, x, y, z, n, t, tstep) + class(bc_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(inout), dimension(n) :: y + real(kind=rp), intent(inout), dimension(n) :: z + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + type(c_ptr) :: x_d + type(c_ptr) :: y_d + type(c_ptr) :: z_d + integer :: i + + + if (NEKO_BCKND_DEVICE .eq. 1) then + x_d = device_get_ptr(x) + y_d = device_get_ptr(y) + z_d = device_get_ptr(z) + if (present(t) .and. present(tstep)) then + call this%apply_vector_dev(x_d, y_d, z_d, t=t, tstep=tstep) + else if (present(t)) then + call this%apply_vector_dev(x_d, y_d, z_d, t=t) + else if (present(tstep)) then + call this%apply_vector_dev(x_d, y_d, z_d, tstep=tstep) + else + call this%apply_vector_dev(x_d, y_d, z_d) + end if + else + if (present(t) .and. present(tstep)) then + call this%apply_vector(x, y, z, n, t=t, tstep=tstep) + else if (present(t)) then + call this%apply_vector(x, y, z, n, t=t) + else if (present(tstep)) then + call this%apply_vector(x, y, z, n, tstep=tstep) + else + call this%apply_vector(x, y, z, n) + end if + end if + + end subroutine bc_apply_vector_generic + + !> Apply the boundary condition to a scalar field. Dispatches to the CPU + !! or the device version. + !! @param x The x comp of the field for which to apply the bc. + !! @param y The y comp of the field for which to apply the bc. + !! @param z The z comp of the field for which to apply the bc. + !! @param n The size of x, y, and z. + !! @param t Current time. + !! @param tstep The current time iteration. + subroutine bc_apply_scalar_generic(this, x, n, t, tstep) + class(bc_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + type(c_ptr) :: x_d + integer :: i + + + if (NEKO_BCKND_DEVICE .eq. 1) then + x_d = device_get_ptr(x) + if (present(t) .and. present(tstep)) then + call this%apply_scalar_dev(x_d, t=t, tstep=tstep) + else if (present(t)) then + call this%apply_scalar_dev(x_d, t=t) + else if (present(tstep)) then + call this%apply_scalar_dev(x_d, tstep=tstep) + else + call this%apply_scalar_dev(x_d) + end if + else + if (present(t) .and. present(tstep)) then + call this%apply_scalar(x, n, t=t, tstep=tstep) + else if (present(t)) then + call this%apply_scalar(x, n, t=t) + else if (present(tstep)) then + call this%apply_scalar(x, n, tstep=tstep) + else + call this%apply_scalar(x, n) + end if + end if + + end subroutine bc_apply_scalar_generic + !> Mark @a facet on element @a el as part of the boundary condition !! @param facet The index of the facet. !! @param el The index of the element. @@ -284,15 +394,14 @@ subroutine bc_mark_zone(this, bc_zone) end do end subroutine bc_mark_zone - !> Mark all facets from a list of zones, also marks type of bc in the mesh. + !> Mark all facets from the list of zones in th mesh. + !! Also marks type of bc in the mesh. !! The facet_type in mesh is because of the fdm from Nek5000... !! That is a hack that should be removed at some point... - !! @param bc_zone Array of boundary zones. !! @param bc_key Boundary condition label, e.g. 'w' for wall. !! @param bc_label List of boundary condition labels. - subroutine bc_mark_zones_from_list(this, bc_zones, bc_key, bc_labels) + subroutine bc_mark_zones_from_list(this, bc_key, bc_labels) class(bc_t), intent(inout) :: this - class(facet_zone_t), intent(inout) :: bc_zones(:) character(len=*) :: bc_key character(len=100), allocatable :: split_key(:) character(len=NEKO_MSH_MAX_ZLBL_LEN) :: bc_labels(NEKO_MSH_MAX_ZLBLS) @@ -318,45 +427,50 @@ subroutine bc_mark_zones_from_list(this, bc_zones, bc_key, bc_labels) msh_bc_type = 2 end if - do i = 1, NEKO_MSH_MAX_ZLBLS - !Check if several bcs are defined for this zone - !bcs are seperated by /, but we could use something else - if (index(trim(bc_labels(i)), '/') .eq. 0) then - if (trim(bc_key) .eq. trim(bc_labels(i))) then - call bc_mark_zone(this, bc_zones(i)) - ! Loop across all faces in the mesh - do j = 1,this%msh%nelv - do k = 1, 2 * this%msh%gdim - if (this%msh%facet_type(k,j) .eq. -i) then - this%msh%facet_type(k,j) = msh_bc_type - end if - end do - end do - end if - else - split_key = split_string(trim(bc_labels(i)),'/') - do l = 1, size(split_key) - if (trim(split_key(l)) .eq. trim(bc_key)) then - call bc_mark_zone(this, bc_zones(i)) - ! Loop across all faces in the mesh - do j = 1,this%msh%nelv - do k = 1, 2 * this%msh%gdim - if (this%msh%facet_type(k,j) .eq. -i) then - this%msh%facet_type(k,j) = msh_bc_type - end if - end do - end do - end if - end do - end if - end do +! do i = 1, NEKO_MSH_MAX_ZLBLS +! !Check if several bcs are defined for this zone +! !bcs are seperated by /, but we could use something else +! +! split_key = split_string(trim(bc_labels(i)),'/') +! do l = 1, size(split_key) +! if (trim(split_key(l)) .eq. trim(bc_key)) then +! call bc_mark_zone(this, this%msh%labeled_zones(i)) +! ! if (index(trim(bc_labels(i)), '/') .eq. 0) then +! ! if (trim(bc_key) .eq. trim(bc_labels(i))) then +! ! call bc_mark_zone(this, bc_zones(i)) +! ! Loop across all faces in the mesh +! do j = 1,this%msh%nelv +! do k = 1, 2 * this%msh%gdim +! if (this%msh%facet_type(k,j) .eq. -i) then +! this%msh%facet_type(k,j) = msh_bc_type +! end if +! end do +! end do +! ! end if +! else +! split_key = split_string(trim(bc_labels(i)),'/') +! do l = 1, size(split_key) +! if (trim(split_key(l)) .eq. trim(bc_key)) then +! call bc_mark_zone(this, bc_zones(i)) +! ! Loop across all faces in the mesh +! do j = 1,this%msh%nelv +! do k = 1, 2 * this%msh%gdim +! if (this%msh%facet_type(k,j) .eq. -i) then +! this%msh%facet_type(k,j) = msh_bc_type +! end if +! end do +! end do +! end if +! end do +! end if +! end do end subroutine bc_mark_zones_from_list !> Finalize the construction of the bc by populting the `msk` and `facet` !! arrays. !! @details This will linearize the marked facet's indicies in the msk array. - subroutine bc_finalize(this) + subroutine bc_finalize_base(this) class(bc_t), target, intent(inout) :: this type(tuple_i4_t), pointer :: bfp(:) type(tuple_i4_t) :: bc_facet @@ -452,188 +566,6 @@ subroutine bc_finalize(this) HOST_TO_DEVICE, sync=.false.) end if - end subroutine bc_finalize - - !> Constructor for a list of boundary conditions - !! @param size The size of the list to allocate. - subroutine bc_list_init(bclst, size) - type(bc_list_t), intent(inout), target :: bclst - integer, optional :: size - integer :: n, i - - call bc_list_free(bclst) - - if (present(size)) then - n = size - else - n = 1 - end if - - allocate(bclst%bc(n)) - - do i = 1, n - bclst%bc(i)%bcp => null() - end do - - bclst%n = 0 - bclst%size = n - - end subroutine bc_list_init - - !> Destructor for a list of boundary conditions - !! @note This will only nullify all pointers, not deallocate any - !! conditions pointed to by the list - subroutine bc_list_free(bclst) - type(bc_list_t), intent(inout) :: bclst - - if (allocated(bclst%bc)) then - deallocate(bclst%bc) - end if - - bclst%n = 0 - bclst%size = 0 - - end subroutine bc_list_free - - !> Add a condition to a list of boundary conditions - !! @param bc The boundary condition to add. - subroutine bc_list_add(bclst, bc) - type(bc_list_t), intent(inout) :: bclst - class(bc_t), intent(inout), target :: bc - type(bcp_t), allocatable :: tmp(:) - - !> Do not add if bc is empty - if(bc%marked_facet%size() .eq. 0) return - - if (bclst%n .ge. bclst%size) then - bclst%size = bclst%size * 2 - allocate(tmp(bclst%size)) - tmp(1:bclst%n) = bclst%bc - call move_alloc(tmp, bclst%bc) - end if - - bclst%n = bclst%n + 1 - bclst%bc(bclst%n)%bcp => bc - - end subroutine bc_list_add - - !> Apply a list of boundary conditions to a scalar field - !! @param x The field to apply the boundary conditions to. - !! @param n The size of x. - !! @param t Current time. - !! @param tstep Current time-step. - subroutine bc_list_apply_scalar(bclst, x, n, t, tstep) - type(bc_list_t), intent(inout) :: bclst - integer, intent(in) :: n - real(kind=rp), intent(inout), dimension(n) :: x - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - type(c_ptr) :: x_d - integer :: i - - if (NEKO_BCKND_DEVICE .eq. 1) then - x_d = device_get_ptr(x) - if (present(t) .and. present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar_dev(x_d, t=t, tstep=tstep) - end do - else if (present(t)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar_dev(x_d, t=t) - end do - else if (present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar_dev(x_d, tstep=tstep) - end do - else - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar_dev(x_d) - end do - end if - else - if (present(t) .and. present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar(x, n, t, tstep) - end do - else if (present(t)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar(x, n, t=t) - end do - else if (present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar(x, n, tstep=tstep) - end do - else - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_scalar(x, n) - end do - end if - end if - - end subroutine bc_list_apply_scalar - - !> Apply a list of boundary conditions to a vector field. - !! @param x The x comp of the field for which to apply the bcs. - !! @param y The y comp of the field for which to apply the bcs. - !! @param z The z comp of the field for which to apply the bcs. - !! @param n The size of x, y, z. - !! @param t Current time. - !! @param tstep Current time-step. - subroutine bc_list_apply_vector(bclst, x, y, z, n, t, tstep) - type(bc_list_t), intent(inout) :: bclst - integer, intent(in) :: n - real(kind=rp), intent(inout), dimension(n) :: x - real(kind=rp), intent(inout), dimension(n) :: y - real(kind=rp), intent(inout), dimension(n) :: z - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - type(c_ptr) :: x_d - type(c_ptr) :: y_d - type(c_ptr) :: z_d - integer :: i - - if (NEKO_BCKND_DEVICE .eq. 1) then - x_d = device_get_ptr(x) - y_d = device_get_ptr(y) - z_d = device_get_ptr(z) - if (present(t) .and. present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d, t, tstep) - end do - else if (present(t)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d, t=t) - end do - else if (present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d, tstep=tstep) - end do - else - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d) - end do - end if - else - if (present(t) .and. present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector(x, y, z, n, t, tstep) - end do - else if (present(t)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector(x, y, z, n, t=t) - end do - else if (present(tstep)) then - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector(x, y, z, n, tstep=tstep) - end do - else - do i = 1, bclst%n - call bclst%bc(i)%bcp%apply_vector(x, y, z, n) - end do - end if - end if - - end subroutine bc_list_apply_vector - + end subroutine bc_finalize_base end module bc diff --git a/src/bc/bc_fctry.f90 b/src/bc/bc_fctry.f90 new file mode 100644 index 00000000000..1a4322d9de9 --- /dev/null +++ b/src/bc/bc_fctry.f90 @@ -0,0 +1,93 @@ +! Copyright (c) 2024, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +!> Defines a factory subroutine for boundary conditions. +module bc_fctry + use bc, only : bc_t + use json_module, only : json_file + use json_utils, only : json_get + use logger, only : neko_log + use dirichlet, only : dirichlet_t + use neumann, only : neumann_t + use coefs, only : coef_t + use facet_zone, only : facet_zone_t + use mesh, only : NEKO_MSH_MAX_ZLBLS + use usr_scalar, only : usr_scalar_t + use user_intf, only : user_t + implicit none + private + + public :: bc_factory + +contains + + !> Boundary condition factory. Both constructs and initializes the object. + !! Will mark a mesh zone for the bc and finalize. + !! @param[in] coef SEM coefficients. + !! @param[inout] json JSON object for initializing the bc. + subroutine bc_factory(object, json, coef, user) + !class(bc_t), allocatable, intent(inout) :: object + class(bc_t), pointer, intent(inout) :: object + type(json_file), intent(inout) :: json + type(coef_t), intent(in) :: coef + type(user_t), intent(in) :: user + character(len=:), allocatable :: type + integer :: zone_index + + call json_get(json, "type", type) + + if (trim(type) .eq. "user_pointwise") then + ! Note, the bc is now in the list even if the mask is zero. + allocate(usr_scalar_t::object) + call object%init(coef, json) + call object%finalize() + + select type(obj => object) + type is(usr_scalar_t) + call obj%set_eval(user%scalar_user_bc) + end select + return + else if (trim(type) .eq. "dirichlet") then + allocate(dirichlet_t::object) + else if (trim(type) .eq. "neumann") then + allocate(neumann_t::object) + end if + + call json_get(json, "zone_index", zone_index) + call object%init(coef, json) + call object%mark_zone(coef%msh%labeled_zones(zone_index)) + call object%finalize() + + end subroutine bc_factory + +end module bc_fctry diff --git a/src/bc/bc_list.f90 b/src/bc/bc_list.f90 new file mode 100644 index 00000000000..46789c66770 --- /dev/null +++ b/src/bc/bc_list.f90 @@ -0,0 +1,304 @@ +! Copyright (c) 2020-2024, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +!> Defines a list of `bc_t`. +module bc_list + use neko_config + use num_types + use device + use, intrinsic :: iso_c_binding, only : c_ptr + use bc, only : bc_t, bc_ptr_t + implicit none + private + + !> A list of allocatable @ref `bc_t`. + !! Follows the standard interface of lists. + type, public :: bc_list_t + type(bc_ptr_t), allocatable :: items(:) + !> Number of items. + integer :: size_ + !> Capacity. + integer :: capacity + contains + !> Constructor. + procedure, pass(this) :: init => bc_list_init + !> Destructor. + procedure, pass(this) :: free => bc_list_free + !> Append an item to the end of the list. + procedure, pass(this) :: append => bc_list_append + !> Apply all boundary conditions in the list. + generic :: apply => apply_scalar, apply_vector + !> Appply the boundary conditions to a scalar field. + procedure, pass(this) :: apply_scalar => bc_list_apply_scalar + !> Appply the boundary conditions to a vector field. + procedure, pass(this) :: apply_vector => bc_list_apply_vector + !> Return the number of items in the list. + procedure, pass(this) :: size => bc_list_size + !> Return wether a given item is a strong bc + procedure, pass(this) :: strong => bc_list_strong + end type bc_list_t + +contains + + !> Constructor. + !! @param size The size of the list to allocate. + subroutine bc_list_init(this, size) + class(bc_list_t), intent(inout), target :: this + integer, optional :: size + integer :: n, i + + call this%free() + + if (present(size)) then + n = size + else + n = 1 + end if + + allocate(this%items(n)) + + do i = 1, n + this%items(i)%ptr => null() + end do + + this%size_ = 0 + this%capacity = n + + end subroutine bc_list_init + + !> Destructor. + !! @note This will only nullify all pointers, not deallocate any + !! conditions pointed to by the list + subroutine bc_list_free(this) + class(bc_list_t), intent(inout) :: this + integer :: i + + if (allocated(this%items)) then + do i =1, this%size() +! call this%items(i)%ptr%free() + nullify(this%items(i)%ptr) + end do + + deallocate(this%items) + end if + + this%size_ = 0 + this%capacity = 0 + end subroutine bc_list_free + + !> Append a condition to the end of the list. + !! @param bc The boundary condition to add. + subroutine bc_list_append(this, bc) + class(bc_list_t), intent(inout) :: this + class(bc_t), intent(inout), target :: bc + type(bc_ptr_t), allocatable :: tmp(:) + + !> Do not add if bc is empty + if(bc%marked_facet%size() .eq. 0) return + + if (this%size_ .ge. this%capacity) then + this%capacity = this%capacity * 2 + allocate(tmp(this%capacity)) + tmp(1:this%size_) = this%items + call move_alloc(tmp, this%items) + end if + + this%size_ = this%size_ + 1 + this%items(this%size_)%ptr => bc + + end subroutine bc_list_append + + !> Apply a list of boundary conditions to a scalar field + !! @param x The field to apply the boundary conditions to. + !! @param n The size of x. + !! @param t Current time. + !! @param tstep Current time-step. + subroutine bc_list_apply_scalar(this, x, n, t, tstep, strong) + class(bc_list_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + logical, intent(in), optional :: strong + type(c_ptr) :: x_d + integer :: i + logical, allocatable :: execute(:) + + allocate(execute(this%size())) + + execute = .true. + if (present(strong)) then + do i=1, this%size() + if (.not. (this%strong(i) .eqv. strong)) then + execute(i) = .false. + end if + end do + end if + + if (NEKO_BCKND_DEVICE .eq. 1) then + x_d = device_get_ptr(x) + if (present(t) .and. present(tstep)) then + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar_dev(x_d, t=t, tstep=tstep) + end if + end do + else if (present(t)) then + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar_dev(x_d, t=t) + end if + end do + else if (present(tstep)) then + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar_dev(x_d, tstep=tstep) + end if + end do + else + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar_dev(x_d) + end if + end do + end if + else + if (present(t) .and. present(tstep)) then + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar(x, n, t, tstep) + end if + end do + else if (present(t)) then + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar(x, n, t=t) + end if + end do + else if (present(tstep)) then + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar(x, n, tstep=tstep) + end if + end do + else + do i = 1, this%size() + if (execute(i)) then + call this%items(i)%ptr%apply_scalar(x, n) + end if + end do + end if + end if + end subroutine bc_list_apply_scalar + + !> Apply a list of boundary conditions to a vector field. + !! @param x The x comp of the field for which to apply the bcs. + !! @param y The y comp of the field for which to apply the bcs. + !! @param z The z comp of the field for which to apply the bcs. + !! @param n The size of x, y, z. + !! @param t Current time. + !! @param tstep Current time-step. + subroutine bc_list_apply_vector(this, x, y, z, n, t, tstep) + class(bc_list_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(inout), dimension(n) :: y + real(kind=rp), intent(inout), dimension(n) :: z + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + type(c_ptr) :: x_d + type(c_ptr) :: y_d + type(c_ptr) :: z_d + integer :: i + + if (NEKO_BCKND_DEVICE .eq. 1) then + x_d = device_get_ptr(x) + y_d = device_get_ptr(y) + z_d = device_get_ptr(z) + if (present(t) .and. present(tstep)) then + do i = 1, this%size() + call this%items(i)%ptr%apply_vector_dev(x_d, y_d, z_d, t, tstep) + end do + else if (present(t)) then + do i = 1, this%size() + call this%items(i)%ptr%apply_vector_dev(x_d, y_d, z_d, t=t) + end do + else if (present(tstep)) then + do i = 1, this%size() + call this%items(i)%ptr%apply_vector_dev(x_d, y_d, z_d, tstep=tstep) + end do + else + do i = 1, this%size() + call this%items(i)%ptr%apply_vector_dev(x_d, y_d, z_d) + end do + end if + else + if (present(t) .and. present(tstep)) then + do i = 1, this%size() + call this%items(i)%ptr%apply_vector(x, y, z, n, t, tstep) + end do + else if (present(t)) then + do i = 1, this%size() + call this%items(i)%ptr%apply_vector(x, y, z, n, t=t) + end do + else if (present(tstep)) then + do i = 1, this%size() + call this%items(i)%ptr%apply_vector(x, y, z, n, tstep=tstep) + end do + else + do i = 1, this%size() + call this%items(i)%ptr%apply_vector(x, y, z, n) + end do + end if + end if + + end subroutine bc_list_apply_vector + + !> Return the number of items in the list. + pure function bc_list_size(this) result(size) + class(bc_list_t), intent(in), target :: this + integer :: size + + size = this%size_ + end function bc_list_size + + !> Return the number of items in the list. + pure function bc_list_strong(this, i) result(strong) + class(bc_list_t), intent(in), target :: this + integer, intent(in) :: i + logical :: strong + + strong = this%items(i)%ptr%strong + end function bc_list_strong + + +end module bc_list diff --git a/src/bc/bcknd/device/cuda/no_slip_wall.cu b/src/bc/bcknd/device/cuda/zero_dirichlet.cu similarity index 84% rename from src/bc/bcknd/device/cuda/no_slip_wall.cu rename to src/bc/bcknd/device/cuda/zero_dirichlet.cu index 318bb126348..3feb7b180d6 100644 --- a/src/bc/bcknd/device/cuda/no_slip_wall.cu +++ b/src/bc/bcknd/device/cuda/zero_dirichlet.cu @@ -1,5 +1,5 @@ /* - Copyright (c) 2021-2022, The Neko Authors + Copyright (c) 2021-2024, The Neko Authors All rights reserved. Redistribution and use in source and binary forms, with or without @@ -32,37 +32,37 @@ POSSIBILITY OF SUCH DAMAGE. */ -#include "no_slip_wall_kernel.h" +#include "zero_dirichlet_kernel.h" #include #include extern "C" { - /** - * Fortran wrapper for device no-slip wall apply vector + /** + * Fortran wrapper for device zero_dirichlet apply vector */ - void cuda_no_slip_wall_apply_scalar(void *msk, void *x, int *m) { + void cuda_zero_dirichlet_apply_scalar(void *msk, void *x, int *m) { const dim3 nthrds(1024, 1, 1); const dim3 nblcks(((*m)+1024 - 1)/ 1024, 1, 1); - no_slip_wall_apply_scalar_kernel + zero_dirichlet_apply_scalar_kernel <<>>((int *) msk, (real *) x, *m); CUDA_CHECK(cudaGetLastError()); } - - /** - * Fortran wrapper for device no-slip wall apply vector + + /** + * Fortran wrapper for device zero_dirichlet wall apply vector */ - void cuda_no_slip_wall_apply_vector(void *msk, void *x, void *y, + void cuda_zero_dirichlet_apply_vector(void *msk, void *x, void *y, void *z, int *m) { const dim3 nthrds(1024, 1, 1); const dim3 nblcks(((*m)+1024 - 1)/ 1024, 1, 1); - no_slip_wall_apply_vector_kernel + zero_dirichlet_apply_vector_kernel <<>>((int *) msk, (real *) x, (real *) y, @@ -70,5 +70,5 @@ extern "C" { *m); CUDA_CHECK(cudaGetLastError()); } - + } diff --git a/src/bc/bcknd/device/cuda/no_slip_wall_kernel.h b/src/bc/bcknd/device/cuda/zero_dirichlet_kernel.h similarity index 69% rename from src/bc/bcknd/device/cuda/no_slip_wall_kernel.h rename to src/bc/bcknd/device/cuda/zero_dirichlet_kernel.h index 25f56e7280d..2390eab0ac4 100644 --- a/src/bc/bcknd/device/cuda/no_slip_wall_kernel.h +++ b/src/bc/bcknd/device/cuda/zero_dirichlet_kernel.h @@ -1,5 +1,5 @@ /* - Copyright (c) 2021-2022, The Neko Authors + Copyright (c) 2021-2024, The Neko Authors All rights reserved. Redistribution and use in source and binary forms, with or without @@ -32,17 +32,17 @@ POSSIBILITY OF SUCH DAMAGE. */ -#ifndef __BC_NO_SLIP_WALL_KERNEL__ -#define __BC_NO_SLIP_WALL_KERNEL__ +#ifndef __BC_ZERO_DIRICHLET_KERNEL__ +#define __BC_ZERO_DIRICHLET_KERNEL__ /** - * Device kernel for scalar apply for a no-slip wall conditon + * Device kernel for scalar apply for a zero-dirichlet conditon */ template< typename T > -__global__ void no_slip_wall_apply_scalar_kernel(const int * __restrict__ msk, - T * __restrict__ x, - const int m) { - +__global__ void zero_dirichlet_apply_scalar_kernel(const int * __restrict__ msk, + T * __restrict__ x, + const int m) { + const int idx = blockIdx.x * blockDim.x + threadIdx.x; const int str = blockDim.x * gridDim.x; @@ -53,15 +53,15 @@ __global__ void no_slip_wall_apply_scalar_kernel(const int * __restrict__ msk, } /** - * Device kernel for vector apply for a no-slip wall conditon + * Device kernel for vector apply for a zero-dirichlet conditon */ template< typename T > -__global__ void no_slip_wall_apply_vector_kernel(const int * __restrict__ msk, - T * __restrict__ x, - T * __restrict__ y, - T * __restrict__ z, - const int m) { - +__global__ void zero_dirichlet_apply_vector_kernel(const int * __restrict__ msk, + T * __restrict__ x, + T * __restrict__ y, + T * __restrict__ z, + const int m) { + const int idx = blockIdx.x * blockDim.x + threadIdx.x; const int str = blockDim.x * gridDim.x; @@ -73,4 +73,4 @@ __global__ void no_slip_wall_apply_vector_kernel(const int * __restrict__ msk, } } -#endif // __BC_NO_SLIP_WALL_KERNEL__ \ No newline at end of file +#endif // __BC_ZERO_DIRICHLET_KERNEL__ \ No newline at end of file diff --git a/src/bc/bcknd/device/device_wall.F90 b/src/bc/bcknd/device/device_zero_dirichlet.F90 similarity index 62% rename from src/bc/bcknd/device/device_wall.F90 rename to src/bc/bcknd/device/device_zero_dirichlet.F90 index d5d7ed78a67..7c9f13f996e 100644 --- a/src/bc/bcknd/device/device_wall.F90 +++ b/src/bc/bcknd/device/device_zero_dirichlet.F90 @@ -30,107 +30,107 @@ ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! -module device_wall +module device_zero_dirichlet use utils, only : neko_error use, intrinsic :: iso_c_binding, only : c_ptr, c_int private #ifdef HAVE_HIP interface - subroutine hip_no_slip_wall_apply_scalar(msk, x, m) & - bind(c, name='hip_no_slip_wall_apply_scalar') + subroutine hip_zero_dirichlet_apply_scalar(msk, x, m) & + bind(c, name='hip_zero_dirichlet_apply_scalar') use, intrinsic :: iso_c_binding implicit none integer(c_int) :: m type(c_ptr), value :: msk, x - end subroutine hip_no_slip_wall_apply_scalar + end subroutine hip_zero_dirichlet_apply_scalar end interface interface - subroutine hip_no_slip_wall_apply_vector(msk, x, y, z, m) & - bind(c, name='hip_no_slip_wall_apply_vector') + subroutine hip_zero_dirichlet_apply_vector(msk, x, y, z, m) & + bind(c, name='hip_zero_dirichlet_apply_vector') use, intrinsic :: iso_c_binding implicit none integer(c_int) :: m type(c_ptr), value :: msk, x, y, z - end subroutine hip_no_slip_wall_apply_vector + end subroutine hip_zero_dirichlet_apply_vector end interface #elif HAVE_CUDA interface - subroutine cuda_no_slip_wall_apply_scalar(msk, x, m) & - bind(c, name='cuda_no_slip_wall_apply_scalar') + subroutine cuda_zero_dirichlet_apply_scalar(msk, x, m) & + bind(c, name='cuda_zero_dirichlet_apply_scalar') use, intrinsic :: iso_c_binding implicit none integer(c_int) :: m type(c_ptr), value :: msk, x - end subroutine cuda_no_slip_wall_apply_scalar + end subroutine cuda_zero_dirichlet_apply_scalar end interface interface - subroutine cuda_no_slip_wall_apply_vector(msk, x, y, z, m) & - bind(c, name='cuda_no_slip_wall_apply_vector') + subroutine cuda_zero_dirichlet_apply_vector(msk, x, y, z, m) & + bind(c, name='cuda_zero_dirichlet_apply_vector') use, intrinsic :: iso_c_binding implicit none integer(c_int) :: m type(c_ptr), value :: msk, x, y, z - end subroutine cuda_no_slip_wall_apply_vector + end subroutine cuda_zero_dirichlet_apply_vector end interface #elif HAVE_OPENCL interface - subroutine opencl_no_slip_wall_apply_scalar(msk, x, m) & - bind(c, name='opencl_no_slip_wall_apply_scalar') + subroutine opencl_zero_dirichlet_apply_scalar(msk, x, m) & + bind(c, name='opencl_zero_dirichlet_apply_scalar') use, intrinsic :: iso_c_binding implicit none integer(c_int) :: m type(c_ptr), value :: msk, x - end subroutine opencl_no_slip_wall_apply_scalar + end subroutine opencl_zero_dirichlet_apply_scalar end interface interface - subroutine opencl_no_slip_wall_apply_vector(msk, x, y, z, m) & - bind(c, name='opencl_no_slip_wall_apply_vector') + subroutine opencl_zero_dirichlet_apply_vector(msk, x, y, z, m) & + bind(c, name='opencl_zero_dirichlet_apply_vector') use, intrinsic :: iso_c_binding implicit none integer(c_int) :: m type(c_ptr), value :: msk, x, y, z - end subroutine opencl_no_slip_wall_apply_vector + end subroutine opencl_zero_dirichlet_apply_vector end interface #endif - public :: device_no_slip_wall_apply_scalar, device_no_slip_wall_apply_vector + public :: device_zero_dirichlet_apply_scalar, device_zero_dirichlet_apply_vector contains - subroutine device_no_slip_wall_apply_scalar(msk, x, m) + subroutine device_zero_dirichlet_apply_scalar(msk, x, m) integer, intent(in) :: m type(c_ptr) :: msk, x #ifdef HAVE_HIP - call hip_no_slip_wall_apply_scalar(msk, x, m) + call hip_zero_dirichlet_apply_scalar(msk, x, m) #elif HAVE_CUDA - call cuda_no_slip_wall_apply_scalar(msk, x, m) + call cuda_zero_dirichlet_apply_scalar(msk, x, m) #elif HAVE_OPENCL - call opencl_no_slip_wall_apply_scalar(msk, x, m) + call opencl_zero_dirichlet_apply_scalar(msk, x, m) #else call neko_error('No device backend configured') #endif - end subroutine device_no_slip_wall_apply_scalar + end subroutine device_zero_dirichlet_apply_scalar - subroutine device_no_slip_wall_apply_vector(msk, x, y, z, m) + subroutine device_zero_dirichlet_apply_vector(msk, x, y, z, m) integer, intent(in) :: m type(c_ptr) :: msk, x, y, z #ifdef HAVE_HIP - call hip_no_slip_wall_apply_vector(msk, x, y, z, m) + call hip_zero_dirichlet_apply_vector(msk, x, y, z, m) #elif HAVE_CUDA - call cuda_no_slip_wall_apply_vector(msk, x, y, z, m) + call cuda_zero_dirichlet_apply_vector(msk, x, y, z, m) #elif HAVE_OPENCL - call opencl_no_slip_wall_apply_vector(msk, x, y, z, m) + call opencl_zero_dirichlet_apply_vector(msk, x, y, z, m) #else call neko_error('No device backend configured') #endif - end subroutine device_no_slip_wall_apply_vector + end subroutine device_zero_dirichlet_apply_vector -end module device_wall +end module device_zero_dirichlet diff --git a/src/bc/bcknd/device/hip/no_slip_wall.hip b/src/bc/bcknd/device/hip/zero_dirichlet.hip similarity index 80% rename from src/bc/bcknd/device/hip/no_slip_wall.hip rename to src/bc/bcknd/device/hip/zero_dirichlet.hip index ddb7741080a..22dfdad41cc 100644 --- a/src/bc/bcknd/device/hip/no_slip_wall.hip +++ b/src/bc/bcknd/device/hip/zero_dirichlet.hip @@ -1,5 +1,5 @@ /* - Copyright (c) 2021-2022, The Neko Authors + Copyright (c) 2021-2024, The Neko Authors All rights reserved. Redistribution and use in source and binary forms, with or without @@ -35,37 +35,37 @@ #include #include #include -#include "no_slip_wall_kernel.h" +#include "zero_dirichlet_kernel.h" extern "C" { - /** - * Fortran wrapper for device no-slip wall apply vector + /** + * Fortran wrapper for device zero-dirichlet apply vector */ - void hip_no_slip_wall_apply_scalar(void *msk, void *x, int *m) { + void hip_zero_dirichlet_apply_scalar(void *msk, void *x, int *m) { const dim3 nthrds(1024, 1, 1); const dim3 nblcks(((*m)+1024 - 1)/ 1024, 1, 1); - hipLaunchKernelGGL(HIP_KERNEL_NAME(no_slip_wall_apply_scalar_kernel), + hipLaunchKernelGGL(HIP_KERNEL_NAME(zero_dirichlet_apply_scalar_kernel), nblcks, nthrds, 0, (hipStream_t) glb_cmd_queue, (int *) msk, (real *) x, *m); HIP_CHECK(hipGetLastError()); } - - /** - * Fortran wrapper for device no-slip wall apply vector + + /** + * Fortran wrapper for device zero-dirichlet apply vector */ - void hip_no_slip_wall_apply_vector(void *msk, void *x, void *y, + void hip_zero_dirichlet_apply_vector(void *msk, void *x, void *y, void *z, int *m) { const dim3 nthrds(1024, 1, 1); const dim3 nblcks(((*m)+1024 - 1)/ 1024, 1, 1); - hipLaunchKernelGGL(HIP_KERNEL_NAME(no_slip_wall_apply_vector_kernel), + hipLaunchKernelGGL(HIP_KERNEL_NAME(zero_dirichlet_apply_vector_kernel), nblcks, nthrds, 0, (hipStream_t) glb_cmd_queue, (int *) msk, (real *) x, (real *) y, (real *) z, *m); HIP_CHECK(hipGetLastError()); } - + } diff --git a/src/bc/bcknd/device/hip/no_slip_wall_kernel.h b/src/bc/bcknd/device/hip/zero_dirichlet_kernel.h similarity index 82% rename from src/bc/bcknd/device/hip/no_slip_wall_kernel.h rename to src/bc/bcknd/device/hip/zero_dirichlet_kernel.h index b83e811ce45..db38c691332 100644 --- a/src/bc/bcknd/device/hip/no_slip_wall_kernel.h +++ b/src/bc/bcknd/device/hip/zero_dirichlet_kernel.h @@ -1,5 +1,5 @@ /* - Copyright (c) 2021-2022, The Neko Authors + Copyright (c) 2021-2024, The Neko Authors All rights reserved. Redistribution and use in source and binary forms, with or without @@ -32,19 +32,19 @@ POSSIBILITY OF SUCH DAMAGE. */ -#ifndef __BC_NO_SLIP_WALL_KERNEL__ -#define __BC_NO_SLIP_WALL_KERNEL__ +#ifndef __BC_ZERO_DIRICHLET_KERNEL__ +#define __BC_ZERO_DIRICHLET_KERNEL__ #include /** - * Device kernel for scalar apply for a no-slip wall conditon + * Device kernel for scalar apply for a zero-dirichlet conditon */ template< typename T > -__global__ void no_slip_wall_apply_scalar_kernel(const int * __restrict__ msk, +__global__ void zero_dirichlet_apply_scalar_kernel(const int * __restrict__ msk, T * __restrict__ x, const int m) { - + const int idx = blockIdx.x * blockDim.x + threadIdx.x; const int str = blockDim.x * gridDim.x; @@ -55,15 +55,15 @@ __global__ void no_slip_wall_apply_scalar_kernel(const int * __restrict__ msk, } /** - * Device kernel for vector apply for a no-slip wall conditon + * Device kernel for vector apply for a zero-dirichlet conditon */ template< typename T > -__global__ void no_slip_wall_apply_vector_kernel(const int * __restrict__ msk, +__global__ void zero_dirichlet_apply_vector_kernel(const int * __restrict__ msk, T * __restrict__ x, T * __restrict__ y, T * __restrict__ z, const int m) { - + const int idx = blockIdx.x * blockDim.x + threadIdx.x; const int str = blockDim.x * gridDim.x; @@ -75,4 +75,4 @@ __global__ void no_slip_wall_apply_vector_kernel(const int * __restrict__ msk, } } -#endif // __BC_NO_SLIP_WALL_KERNEL__ \ No newline at end of file +#endif // __BC_ZERO_DIRICHLET_KERNEL__ \ No newline at end of file diff --git a/src/bc/bcknd/device/opencl/no_slip_wall.c b/src/bc/bcknd/device/opencl/zero_dirichlet.c similarity index 94% rename from src/bc/bcknd/device/opencl/no_slip_wall.c rename to src/bc/bcknd/device/opencl/zero_dirichlet.c index 394f6a6ae4c..42a6db26016 100644 --- a/src/bc/bcknd/device/opencl/no_slip_wall.c +++ b/src/bc/bcknd/device/opencl/zero_dirichlet.c @@ -1,5 +1,5 @@ /* - Copyright (c) 2021-2022, The Neko Authors + Copyright (c) 2021-2024, The Neko Authors All rights reserved. Redistribution and use in source and binary forms, with or without @@ -44,18 +44,18 @@ #include #include -#include "no_slip_wall_kernel.cl.h" +#include "zero_dirichlet_kernel.cl.h" -/** - * Fortran wrapper for device no-slip wall apply vector +/** + * Fortran wrapper for device zero-dirichlet apply vector */ void opencl_no_slip_wall_apply_scalar(void *msk, void *x, int *m) { - + cl_int err; - + if (no_slip_wall_program == NULL) opencl_kernel_jit(no_slip_wall_kernel, (cl_program *) &no_slip_wall_program); - + cl_kernel kernel = clCreateKernel(no_slip_wall_program, "no_slip_wall_apply_scalar_kernel", &err); CL_CHECK(err); @@ -63,7 +63,7 @@ void opencl_no_slip_wall_apply_scalar(void *msk, void *x, int *m) { CL_CHECK(clSetKernelArg(kernel, 0, sizeof(cl_mem), (void *) &msk)); CL_CHECK(clSetKernelArg(kernel, 1, sizeof(cl_mem), (void *) &x)); CL_CHECK(clSetKernelArg(kernel, 2, sizeof(int), m)); - + const int nb = ((*m) + 256 - 1) / 256; const size_t global_item_size = 256 * nb; const size_t local_item_size = 256; @@ -73,16 +73,16 @@ void opencl_no_slip_wall_apply_scalar(void *msk, void *x, int *m) { 0, NULL, NULL)); } -/** - * Fortran wrapper for device no-slip wall apply vector +/** + * Fortran wrapper for device zero-dirichlet apply vector */ void opencl_no_slip_wall_apply_vector(void *msk, void *x, void *y, void *z, int *m) { cl_int err; - + if (no_slip_wall_program == NULL) opencl_kernel_jit(no_slip_wall_kernel, (cl_program *) &no_slip_wall_program); - + cl_kernel kernel = clCreateKernel(no_slip_wall_program, "no_slip_wall_apply_vector_kernel", &err); CL_CHECK(err); @@ -92,7 +92,7 @@ void opencl_no_slip_wall_apply_vector(void *msk, void *x, void *y, CL_CHECK(clSetKernelArg(kernel, 2, sizeof(cl_mem), (void *) &y)); CL_CHECK(clSetKernelArg(kernel, 3, sizeof(cl_mem), (void *) &z)); CL_CHECK(clSetKernelArg(kernel, 4, sizeof(int), m)); - + const int nb = ((*m) + 256 - 1) / 256; const size_t global_item_size = 256 * nb; const size_t local_item_size = 256; @@ -100,5 +100,5 @@ void opencl_no_slip_wall_apply_vector(void *msk, void *x, void *y, CL_CHECK(clEnqueueNDRangeKernel((cl_command_queue) glb_cmd_queue, kernel, 1, NULL, &global_item_size, &local_item_size, 0, NULL, NULL)); - + } diff --git a/src/bc/bcknd/device/opencl/no_slip_wall_kernel.cl b/src/bc/bcknd/device/opencl/zero_dirichlet_kernel.cl similarity index 70% rename from src/bc/bcknd/device/opencl/no_slip_wall_kernel.cl rename to src/bc/bcknd/device/opencl/zero_dirichlet_kernel.cl index 9241165b9cf..c871313ce65 100644 --- a/src/bc/bcknd/device/opencl/no_slip_wall_kernel.cl +++ b/src/bc/bcknd/device/opencl/zero_dirichlet_kernel.cl @@ -32,16 +32,16 @@ POSSIBILITY OF SUCH DAMAGE. */ -#ifndef __BC_NO_SLIP_WALL_KERNEL__ -#define __BC_NO_SLIP_WALL_KERNEL__ +#ifndef __BC_ZERO_DIRICHLET_KERNEL__ +#define __BC_ZERO_DIRICHLET_KERNEL__ /** - * Device kernel for scalar apply for a no-slip wall conditon + * Device kernel for scalar apply for a zero-valued Dirichlet conditon */ -__kernel void no_slip_wall_apply_scalar_kernel(__global const int *msk, - __global real *x, - const int m) { - +__kernel void zero_dirichlet_apply_scalar_kernel(__global const int *msk, + __global real *x, + const int m) { + const int idx = get_global_id(0); const int str = get_global_size(0); @@ -52,14 +52,14 @@ __kernel void no_slip_wall_apply_scalar_kernel(__global const int *msk, } /** - * Device kernel for vector apply for a no-slip wall conditon + * Device kernel for vector apply for a zero-valued Dirichlet conditon */ -__kernel void no_slip_wall_apply_vector_kernel(__global const int *msk, - __global real *x, - __global real *y, - __global real *z, - const int m) { - +__kernel void zero_dirichlet_apply_vector_kernel(__global const int *msk, + __global real *x, + __global real *y, + __global real *z, + const int m) { + const int idx = get_global_id(0); const int str = get_global_size(0); @@ -71,4 +71,4 @@ __kernel void no_slip_wall_apply_vector_kernel(__global const int *msk, } } -#endif // __BC_NO_SLIP_WALL_KERNEL__ \ No newline at end of file +#endif // __BC_ZERO_DIRICHLET_KERNEL__ \ No newline at end of file diff --git a/src/bc/blasius.f90 b/src/bc/blasius.f90 index 81d04481c10..a70d155819a 100644 --- a/src/bc/blasius.f90 +++ b/src/bc/blasius.f90 @@ -41,10 +41,14 @@ module blasius use, intrinsic :: iso_fortran_env use, intrinsic :: iso_c_binding use bc, only : bc_t + use json_module, only : json_file + use json_utils, only : json_get implicit none private - !> Blasius profile for inlet (vector valued) + !> Blasius profile for inlet (vector valued). + !! @warning Works only with axis-aligned sugar-cube elements and assumes + !! the boundary is alinged with zOy. type, public, extends(bc_t) :: blasius_t real(kind=rp), dimension(3) :: uinf = (/0d0, 0d0, 0d0 /) real(kind=rp) :: delta @@ -58,12 +62,54 @@ module blasius procedure, pass(this) :: apply_scalar_dev => blasius_apply_scalar_dev procedure, pass(this) :: apply_vector_dev => blasius_apply_vector_dev procedure, pass(this) :: set_params => blasius_set_params + !> Constructor + procedure, pass(this) :: init => blasius_init !> Destructor. procedure, pass(this) :: free => blasius_free + !> Finalize. + procedure, pass(this) :: finalize => blasius_finalize end type blasius_t contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine blasius_init(this, coef, json) + class(blasius_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) :: json + real(kind=rp) :: delta + real(kind=rp), allocatable :: uinf(:) + character(len=:), allocatable :: approximation + + call this%init_base(coef) + + call json_get(json, 'case.fluid.blasius.delta', delta) + call json_get(json, 'case.fluid.blasius.approximation', approximation) + call json_get(json, 'case.fluid.blasius.freestream_velocity', uinf) + + this%delta = delta + this%uinf = uinf + + select case(trim(approximation)) + case('linear') + this%bla => blasius_linear + case('quadratic') + this%bla => blasius_quadratic + case('cubic') + this%bla => blasius_cubic + case('quartic') + this%bla => blasius_quartic + case('sin') + this%bla => blasius_sin + case default + call neko_error('Invalid Blasius approximation') + end select + + + end subroutine blasius_init + subroutine blasius_free(this) class(blasius_t), target, intent(inout) :: this @@ -237,4 +283,10 @@ subroutine blasius_set_params(this, uinf, delta, type) end select end subroutine blasius_set_params + !> Finalize + subroutine blasius_finalize(this) + class(blasius_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine blasius_finalize end module blasius diff --git a/src/bc/dirichlet.f90 b/src/bc/dirichlet.f90 index 8cf6a53da33..52df385b944 100644 --- a/src/bc/dirichlet.f90 +++ b/src/bc/dirichlet.f90 @@ -35,6 +35,9 @@ module dirichlet use device_dirichlet use num_types, only : rp use bc, only : bc_t + use coefs, only : coef_t + use json_module, only : json_file + use json_utils, only : json_get use, intrinsic :: iso_c_binding, only : c_ptr implicit none private @@ -49,12 +52,46 @@ module dirichlet procedure, pass(this) :: apply_scalar_dev => dirichlet_apply_scalar_dev procedure, pass(this) :: apply_vector_dev => dirichlet_apply_vector_dev procedure, pass(this) :: set_g => dirichlet_set_g + !> Constructor from JSON. + procedure, pass(this) :: init => dirichlet_init + !> Constructor from components. + procedure, pass(this) :: init_from_components => & + dirichlet_init_from_components !> Destructor. procedure, pass(this) :: free => dirichlet_free + !> Finalize. + procedure, pass(this) :: finalize => dirichlet_finalize end type dirichlet_t contains + !> Constructor from JSON. + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine dirichlet_init(this, coef, json) + class(dirichlet_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + real(kind=rp) :: g + + call this%init_base(coef) + call json_get(json , "value", g) + + this%g = g + end subroutine dirichlet_init + + !> Constructor from components. + !! @param[in] coef The SEM coefficients. + !! @param[in] g The value to apply at the boundary. + subroutine dirichlet_init_from_components(this, coef, g) + class(dirichlet_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + real(kind=rp), intent(in) :: g + + call this%init_base(coef) + this%g = g + end subroutine dirichlet_init_from_components + !> Boundary condition apply for a generic Dirichlet condition !! to a vector @a x subroutine dirichlet_apply_scalar(this, x, n, t, tstep) @@ -139,4 +176,11 @@ subroutine dirichlet_free(this) end subroutine dirichlet_free + !> Finalize + subroutine dirichlet_finalize(this) + class(dirichlet_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine dirichlet_finalize + end module dirichlet diff --git a/src/bc/dong_outflow.f90 b/src/bc/dong_outflow.f90 index 1a478b1778d..dc665e77c5d 100644 --- a/src/bc/dong_outflow.f90 +++ b/src/bc/dong_outflow.f90 @@ -68,14 +68,20 @@ module dong_outflow procedure, pass(this) :: apply_vector => dong_outflow_apply_vector procedure, pass(this) :: apply_scalar_dev => dong_outflow_apply_scalar_dev procedure, pass(this) :: apply_vector_dev => dong_outflow_apply_vector_dev + !> Constructor procedure, pass(this) :: init => dong_outflow_init !> Destructor. procedure, pass(this) :: free => dong_outflow_free + !> Finalize. + procedure, pass(this) :: finalize => dong_outflow_finalize end type dong_outflow_t contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. subroutine dong_outflow_init(this, coef, json) - class(dong_outflow_t), intent(inout) :: this + class(dong_outflow_t), target, intent(inout) :: this type(coef_t), intent(in) :: coef type(json_file), intent(inout) :: json real(kind=rp), allocatable :: temp_x(:) @@ -85,7 +91,8 @@ subroutine dong_outflow_init(this, coef, json) integer :: i, m, k, facet, idx(4) real(kind=rp) :: normal_xyz(3) -! call this%dirichlet_t%init + call this%free() + call this%init_base(coef) call json_get_or_default(json, 'case.fluid.outflow_condition.delta', & this%delta, 0.01_rp) @@ -203,4 +210,11 @@ subroutine dong_outflow_free(this) end subroutine dong_outflow_free + !> Finalize + subroutine dong_outflow_finalize(this) + class(dong_outflow_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine dong_outflow_finalize + end module dong_outflow diff --git a/src/bc/facet_normal.f90 b/src/bc/facet_normal.f90 index 96ac1f204e6..fa9133efeee 100644 --- a/src/bc/facet_normal.f90 +++ b/src/bc/facet_normal.f90 @@ -38,6 +38,7 @@ module facet_normal use coefs, only : coef_t use bc, only : bc_t use utils + use json_module, only : json_file use, intrinsic :: iso_c_binding, only : c_ptr implicit none private @@ -51,12 +52,27 @@ module facet_normal procedure, pass(this) :: apply_vector_dev => facet_normal_apply_vector_dev procedure, pass(this) :: apply_surfvec => facet_normal_apply_surfvec procedure, pass(this) :: apply_surfvec_dev => facet_normal_apply_surfvec_dev + !> Constructor + procedure, pass(this) :: init => facet_normal_init !> Destructor. procedure, pass(this) :: free => facet_normal_free + !> Finalize. + procedure, pass(this) :: finalize => facet_normal_finalize end type facet_normal_t contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine facet_normal_init(this, coef, json) + class(facet_normal_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + + call this%init_base(coef) + end subroutine facet_normal_init + !> No-op scalar apply subroutine facet_normal_apply_scalar(this, x, n, t, tstep) class(facet_normal_t), intent(inout) :: this @@ -170,4 +186,11 @@ subroutine facet_normal_free(this) end subroutine facet_normal_free + !> Finalize + subroutine facet_normal_finalize(this) + class(facet_normal_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine facet_normal_finalize + end module facet_normal diff --git a/src/bc/field_dirichlet.f90 b/src/bc/field_dirichlet.f90 index e7154f67ada..0c4cd3ecda4 100644 --- a/src/bc/field_dirichlet.f90 +++ b/src/bc/field_dirichlet.f90 @@ -30,12 +30,13 @@ ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! -!> Defines inflow dirichlet conditions +!> Defines user dirichlet condition for a scalar field. module field_dirichlet use num_types, only: rp use coefs, only: coef_t use dirichlet, only: dirichlet_t - use bc, only: bc_list_t, bc_t + use bc, only: bc_t + use bc_list, only : bc_list_t use device, only: c_ptr, c_size_t use utils, only: split_string use field, only : field_t @@ -44,6 +45,7 @@ module field_dirichlet use device_math, only: device_masked_copy use dofmap, only : dofmap_t use utils, only: neko_error + use json_module, only : json_file use field_list, only : field_list_t implicit none private @@ -68,8 +70,8 @@ module field_dirichlet !! of the boundary fields. procedure(field_dirichlet_update), nopass, pointer :: update => null() contains - !> Constructor. - procedure, pass(this) :: init_field => field_dirichlet_init + !> Constructor for the field. + procedure, pass(this) :: init_field => field_dirichlet_init_field !> Apply scalar by performing a masked copy. procedure, pass(this) :: apply_scalar => field_dirichlet_apply_scalar !> (No-op) Apply vector. @@ -78,8 +80,12 @@ module field_dirichlet procedure, pass(this) :: apply_vector_dev => field_dirichlet_apply_vector_dev !> Apply scalar (device). procedure, pass(this) :: apply_scalar_dev => field_dirichlet_apply_scalar_dev - !> Destructor + !> Constructor. + procedure, pass(this) :: init => field_dirichlet_init + !> Destructor. procedure, pass(this) :: free => field_dirichlet_free + !> Finalize. + procedure, pass(this) :: finalize => field_dirichlet_finalize end type field_dirichlet_t @@ -110,10 +116,20 @@ end subroutine field_dirichlet_update public :: field_dirichlet_update contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine field_dirichlet_init(this, coef, json) + class(field_dirichlet_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + + call this%init_base(coef) + end subroutine field_dirichlet_init !> Initializes this%field_bc. !! @param bc_name Name of this%field_bc - subroutine field_dirichlet_init(this, bc_name) + subroutine field_dirichlet_init_field(this, bc_name) class(field_dirichlet_t), target, intent(inout) :: this character(len=*), intent(in) :: bc_name @@ -121,7 +137,7 @@ subroutine field_dirichlet_init(this, bc_name) call this%field_list%init(1) call this%field_list%assign_to_field(1, this%field_bc) - end subroutine field_dirichlet_init + end subroutine field_dirichlet_init_field !> Destructor. Currently this%field_bc is being freed in `fluid_scheme::free` subroutine field_dirichlet_free(this) @@ -212,4 +228,10 @@ subroutine field_dirichlet_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) end subroutine field_dirichlet_apply_vector_dev + !> Finalize + subroutine field_dirichlet_finalize(this) + class(field_dirichlet_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine field_dirichlet_finalize end module field_dirichlet diff --git a/src/bc/field_dirichlet_vector.f90 b/src/bc/field_dirichlet_vector.f90 index 2c736b55cc4..8d68cb40f6d 100644 --- a/src/bc/field_dirichlet_vector.f90 +++ b/src/bc/field_dirichlet_vector.f90 @@ -35,7 +35,8 @@ module field_dirichlet_vector use num_types, only: rp use coefs, only: coef_t use dirichlet, only: dirichlet_t - use bc, only: bc_list_t, bc_t, bc_list_free + use bc, only: bc_t + use bc_list, only: bc_list_t use device, only: c_ptr, c_size_t use utils, only: split_string use field, only : field_t @@ -45,6 +46,7 @@ module field_dirichlet_vector use dofmap, only : dofmap_t use field_dirichlet, only: field_dirichlet_t, field_dirichlet_update use utils, only: neko_error + use json_module, only : json_file use field_list, only : field_list_t implicit none private @@ -67,9 +69,13 @@ module field_dirichlet_vector procedure(field_dirichlet_update), nopass, pointer :: update => null() contains !> Initializes this%field_bc. - procedure, pass(this) :: init_field => field_dirichlet_vector_init + procedure, pass(this) :: init_field => field_dirichlet_vector_init_field + !> Constructor + procedure, pass(this) :: init => field_dirichlet_vector_init !> Destructor procedure, pass(this) :: free => field_dirichlet_vector_free + !> Finalize. + procedure, pass(this) :: finalize => field_dirichlet_vector_finalize !> Apply scalar by performing a masked copy. procedure, pass(this) :: apply_scalar => field_dirichlet_vector_apply_scalar !> (No-op) Apply vector. @@ -84,14 +90,25 @@ module field_dirichlet_vector contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine field_dirichlet_vector_init(this, coef, json) + class(field_dirichlet_vector_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + + call this%init_base(coef) + end subroutine field_dirichlet_vector_init + !> Initializes this%field_bc. - subroutine field_dirichlet_vector_init(this, bc_name) + subroutine field_dirichlet_vector_init_field(this, bc_name) class(field_dirichlet_vector_t), intent(inout) :: this character(len=*), intent(in) :: bc_name call neko_error("Fields must be initialized individually!") - end subroutine field_dirichlet_vector_init + end subroutine field_dirichlet_vector_init_field !> Destructor. Currently unused as is, all field_dirichlet attributes !! are freed in `fluid_scheme::free`. @@ -103,7 +120,7 @@ subroutine field_dirichlet_vector_free(this) call this%bc_w%free() call this%field_list%free() - call bc_list_free(this%bc_list) + call this%bc_list%free() if (associated(this%update)) then nullify(this%update) @@ -206,4 +223,11 @@ subroutine field_dirichlet_vector_apply_vector_dev(this, x_d, y_d, z_d, t, tstep end subroutine field_dirichlet_vector_apply_vector_dev + !> Finalize + subroutine field_dirichlet_vector_finalize(this) + class(field_dirichlet_vector_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine field_dirichlet_vector_finalize + end module field_dirichlet_vector diff --git a/src/bc/inflow.f90 b/src/bc/inflow.f90 index b564bdec40f..284cfbc046a 100644 --- a/src/bc/inflow.f90 +++ b/src/bc/inflow.f90 @@ -36,6 +36,9 @@ module inflow use num_types, only : rp use bc, only : bc_t use, intrinsic :: iso_c_binding, only : c_ptr, c_loc + use coefs, only : coef_t + use json_module, only : json_file + use json_utils, only : json_get implicit none private @@ -47,13 +50,30 @@ module inflow procedure, pass(this) :: apply_vector => inflow_apply_vector procedure, pass(this) :: apply_scalar_dev => inflow_apply_scalar_dev procedure, pass(this) :: apply_vector_dev => inflow_apply_vector_dev - procedure, pass(this) :: set_inflow => inflow_set_vector + !> Constructor + procedure, pass(this) :: init => inflow_init !> Destructor. procedure, pass(this) :: free => inflow_free + !> Finalize. + procedure, pass(this) :: finalize => inflow_finalize end type inflow_t contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine inflow_init(this, coef, json) + class(inflow_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + real(kind=rp), allocatable :: x(:) + + call this%init_base(coef) + call json_get(json, 'case.fluid.inflow_condition.value', x) + this%x = x + end subroutine inflow_init + !> No-op scalar apply subroutine inflow_apply_scalar(this, x, n, t, tstep) class(inflow_t), intent(inout) :: this @@ -105,20 +125,19 @@ subroutine inflow_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) end subroutine inflow_apply_vector_dev - !> Set inflow vector - subroutine inflow_set_vector(this, x) - class(inflow_t), intent(inout) :: this - real(kind=rp), dimension(3), intent(inout) :: x - this%x = x - end subroutine inflow_set_vector - !> Destructor subroutine inflow_free(this) class(inflow_t), target, intent(inout) :: this call this%free_base() - end subroutine inflow_free + !> Finalize + subroutine inflow_finalize(this) + class(inflow_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine inflow_finalize + end module inflow diff --git a/src/bc/neumann.f90 b/src/bc/neumann.f90 index ad261d05584..9a33d9ccbe1 100644 --- a/src/bc/neumann.f90 +++ b/src/bc/neumann.f90 @@ -37,6 +37,8 @@ module neumann use, intrinsic :: iso_c_binding, only : c_ptr use utils, only : neko_error, nonlinear_index use coefs, only : coef_t + use json_module, only : json_file + use json_utils, only : json_get use math, only : cfill implicit none private @@ -47,19 +49,54 @@ module neumann !! to the right-hand-side. type, public, extends(bc_t) :: neumann_t real(kind=rp), allocatable, private :: flux_(:) + real(kind=rp), private :: init_flux_ contains procedure, pass(this) :: apply_scalar => neumann_apply_scalar procedure, pass(this) :: apply_vector => neumann_apply_vector procedure, pass(this) :: apply_scalar_dev => neumann_apply_scalar_dev procedure, pass(this) :: apply_vector_dev => neumann_apply_vector_dev - procedure, pass(this) :: finalize_neumann => neumann_finalize_neumann procedure, pass(this) :: flux => neumann_flux + !> Constructor + procedure, pass(this) :: init => neumann_init + !> Constructor from components + procedure, pass(this) :: init_from_components => & + neumann_init_from_components !> Destructor. procedure, pass(this) :: free => neumann_free + !> Finalize. + procedure, pass(this) :: finalize => neumann_finalize end type neumann_t contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine neumann_init(this, coef, json) + class(neumann_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) :: json + real(kind=rp) :: flux + + call this%init_base(coef) + this%strong = .false. + + call json_get(json, "value", flux) + this%init_flux_ = flux + end subroutine neumann_init + + !> Constructor from components. + !! @param[in] coef The SEM coefficients. + !! @param[in] g The value to apply at the boundary. + subroutine neumann_init_from_components(this, coef, flux) + class(neumann_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + real(kind=rp), intent(in) :: flux + + call this%init_base(coef) + this%init_flux_ = flux + end subroutine neumann_init_from_components + !> Boundary condition apply for a generic Neumann condition !! to a vector @a x subroutine neumann_apply_scalar(this, x, n, t, tstep) @@ -138,6 +175,16 @@ subroutine neumann_free(this) end subroutine neumann_free + !> Finalize by setting the flux. + subroutine neumann_finalize(this) + class(neumann_t), target, intent(inout) :: this + + call this%finalize_base() + allocate(this%flux_(this%msk(0))) + + call cfill(this%flux_, this%init_flux_, this%msk(0)) + end subroutine neumann_finalize + !> Get the flux. pure function neumann_flux(this) result(flux) class(neumann_t), intent(in) :: this @@ -146,16 +193,4 @@ pure function neumann_flux(this) result(flux) flux = this%flux_ end function neumann_flux - !> Finalize by setting the flux - !> @param flux The desired flux. - subroutine neumann_finalize_neumann(this, flux) - class(neumann_t), intent(inout) :: this - real(kind=rp), intent(in) :: flux - - call this%finalize() - allocate(this%flux_(this%msk(0))) - - call cfill(this%flux_, flux, this%msk(0)) - end subroutine neumann_finalize_neumann - end module neumann diff --git a/src/bc/non_normal.f90 b/src/bc/non_normal.f90 index c2d6488ebb3..767baee1515 100644 --- a/src/bc/non_normal.f90 +++ b/src/bc/non_normal.f90 @@ -42,39 +42,52 @@ module non_normal use math use utils use stack + use json_module, only : json_file use, intrinsic :: iso_c_binding implicit none private - !> Dirichlet condition in non normal direction of a plane + !> Dirichlet condition in non normal direction of a plane. + !! @warning Only works for axis-aligned plane boundaries. type, public, extends(symmetry_t) :: non_normal_t contains !> Constructor. procedure, pass(this) :: init => non_normal_init !> Destructor. procedure, pass(this) :: free => non_normal_free + !> Finalize. + procedure, pass(this) :: finalize => non_normal_finalize end type non_normal_t contains - !> Constructor. - subroutine non_normal_init(this, coef) - class(non_normal_t), intent(inout) :: this + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine non_normal_init(this, coef, json) + class(non_normal_t), target, intent(inout) :: this type(coef_t), intent(in) :: coef - integer :: i, j, l + type(json_file), intent(inout) ::json + + call this%free() + + call this%init_base(coef) + call this%bc_x%init(this%coef, json) + call this%bc_y%init(this%coef, json) + call this%bc_z%init(this%coef, json) + + end subroutine non_normal_init + + !> Finalize + subroutine non_normal_finalize(this) + class(non_normal_t), target, intent(inout) :: this + integer :: i, j, k, l type(tuple_i4_t), pointer :: bfp(:) real(kind=rp) :: sx,sy,sz real(kind=rp), parameter :: TOL = 1d-3 type(tuple_i4_t) :: bc_facet integer :: facet, el - call this%free() - - call this%init_base(coef) - call this%bc_x%init_base(this%coef) - call this%bc_y%init_base(this%coef) - call this%bc_z%init_base(this%coef) - associate(c=>this%coef, nx => this%coef%nx, ny => this%coef%ny, & nz => this%coef%nz) bfp => this%marked_facet%array() @@ -82,38 +95,7 @@ subroutine non_normal_init(this, coef) bc_facet = bfp(i) facet = bc_facet%x(1) el = bc_facet%x(2) - sx = 0d0 - sy = 0d0 - sz = 0d0 - select case (facet) - case(1,2) - do l = 2, c%Xh%lx - 1 - do j = 2, c%Xh%lx -1 - sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) - sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) - sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) - end do - end do - case(3,4) - do l = 2, c%Xh%lx - 1 - do j = 2, c%Xh%lx - 1 - sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) - sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) - sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) - end do - end do - case(5,6) - do l = 2, c%Xh%lx - 1 - do j = 2, c%Xh%lx - 1 - sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) - sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) - sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) - end do - end do - end select - sx = sx / (c%Xh%lx - 2)**2 - sy = sy / (c%Xh%lx - 2)**2 - sz = sz / (c%Xh%lx - 2)**2 + call this%get_normal_axis(sx, sy, sz, facet, el) if (sx .lt. TOL) then call this%bc_y%mark_facet(facet, el) @@ -132,12 +114,11 @@ subroutine non_normal_init(this, coef) end do end associate call this%bc_x%finalize() - call this%bc_x%set_g(0.0_rp) call this%bc_y%finalize() - call this%bc_y%set_g(0.0_rp) call this%bc_z%finalize() - call this%bc_z%set_g(0.0_rp) - end subroutine non_normal_init + + call this%finalize_base() + end subroutine non_normal_finalize !> Destructor subroutine non_normal_free(this) diff --git a/src/bc/shear_stress.f90 b/src/bc/shear_stress.f90 index 5250688e5e7..8a3b990ae50 100644 --- a/src/bc/shear_stress.f90 +++ b/src/bc/shear_stress.f90 @@ -32,176 +32,177 @@ ! !> Defines a shear stress boundary condition for a vector field. module shear_stress - use num_types - use bc, only : bc_t - use, intrinsic :: iso_c_binding, only : c_ptr - use utils, only : neko_error, nonlinear_index - use coefs, only : coef_t - use dirichlet, only : dirichlet_t - use neumann, only : neumann_t - use math, only : copy - use device_math, only : device_copy - implicit none - private - - !> A shear stress boundary condition. - !! @warning Currently hard coded to prescribe a given stress in x, 0 stress - !! in z, and 0 velocity in y. - type, public, extends(bc_t) :: shear_stress_t - !> The stress in the 1st wall-parallel direction. - real(kind=rp), allocatable, private :: tau1_(:) - !> The stress in the 2nd wall-parallel direction. - real(kind=rp), allocatable, private :: tau2_(:) - !> Dirchlet bc for the wall-normal component. - type(dirichlet_t) :: dirichlet - !> Neumann condition for the 1st wall-parallel direction. - type(neumann_t) :: neumann1 - !> Neumann condition for the 2nd wall-parallel direction. - type(neumann_t) :: neumann2 - contains - procedure, pass(this) :: apply_scalar => shear_stress_apply_scalar - procedure, pass(this) :: apply_vector => shear_stress_apply_vector - procedure, pass(this) :: apply_scalar_dev => shear_stress_apply_scalar_dev - procedure, pass(this) :: apply_vector_dev => shear_stress_apply_vector_dev - procedure, pass(this) :: init_shear_stress => & - shear_stress_init_shear_stress - procedure, pass(this) :: tau1 => shear_stress_tau1 - procedure, pass(this) :: tau2 => shear_stress_tau2 - procedure, pass(this) :: set_stress => shear_stress_set_stress - procedure, pass(this) :: shear_stress_finalize - procedure, pass(this) :: free => shear_stress_free - end type shear_stress_t - - contains - - !> Apply shear stress for a scalar field @a x. - subroutine shear_stress_apply_scalar(this, x, n, t, tstep) - class(shear_stress_t), intent(inout) :: this - integer, intent(in) :: n - real(kind=rp), intent(inout), dimension(n) :: x - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - integer :: i, m, k, facet - ! Store non-linear index - integer :: idx(4) - - call neko_error("The shear stress bc is not applicable to scalar fields.") - - end subroutine shear_stress_apply_scalar - - !> Boundary condition apply for a generic shear_stress condition - !! to vectors @a x, @a y and @a z - subroutine shear_stress_apply_vector(this, x, y, z, n, t, tstep) - class(shear_stress_t), intent(inout) :: this - integer, intent(in) :: n - real(kind=rp), intent(inout), dimension(n) :: x - real(kind=rp), intent(inout), dimension(n) :: y - real(kind=rp), intent(inout), dimension(n) :: z - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - - call this%neumann1%apply_scalar(x, n, t, tstep) - call this%neumann2%apply_scalar(z, n, t, tstep) - call this%dirichlet%apply_scalar(y, n, t, tstep) - - end subroutine shear_stress_apply_vector - - !> Boundary condition apply for a generic shear_stress condition - !! to a vector @a x (device version) - subroutine shear_stress_apply_scalar_dev(this, x_d, t, tstep) - class(shear_stress_t), intent(inout), target :: this - type(c_ptr) :: x_d - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - - call neko_error("shear_stress bc not implemented on the device") - - end subroutine shear_stress_apply_scalar_dev - - !> Boundary condition apply for a generic shear_stress condition - !! to vectors @a x, @a y and @a z (device version) - subroutine shear_stress_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) - class(shear_stress_t), intent(inout), target :: this - type(c_ptr) :: x_d - type(c_ptr) :: y_d - type(c_ptr) :: z_d - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - - call neko_error("shear_stress bc not implemented on the device") - - end subroutine shear_stress_apply_vector_dev - - !> Constructor. - !> @param tau1 The desired stress in the 1st wall-parallel direction. - !> @param tau2 The desired stress in the 2nd wall-parallel direction. - !> @param coef The SEM coefficients. - subroutine shear_stress_init_shear_stress(this, coef) - class(shear_stress_t), intent(inout) :: this - type(coef_t), target, intent(in) :: coef - - call this%init_base(coef) - - call this%dirichlet%init_base(coef) - call this%neumann1%init_base(coef) - call this%neumann2%init_base(coef) - - call this%dirichlet%set_g(0.0_rp) - - this%coef => coef - end subroutine shear_stress_init_shear_stress - - !> Finalize by allocating the stress arrays and marking the facets for - !! the bc components. - subroutine shear_stress_finalize(this, tau1, tau2) - class(shear_stress_t), intent(inout) :: this - real(kind=rp), intent(in) :: tau1 - real(kind=rp), intent(in) :: tau2 - - call this%neumann1%mark_facets(this%marked_facet) - call this%neumann2%mark_facets(this%marked_facet) - call this%dirichlet%mark_facets(this%marked_facet) - - call this%dirichlet%finalize() - call this%neumann1%finalize_neumann(tau1) - call this%neumann2%finalize_neumann(tau2) - - end subroutine shear_stress_finalize - - !> Get the stress in the 1st wall-parallel direction. - pure function shear_stress_tau1(this) result(tau1) - class(shear_stress_t), intent(in) :: this - real(kind=rp) :: tau1(this%msk(0)) - - tau1 = this%tau1_ - end function shear_stress_tau1 - - !> Get the stress in the 2nd wall-parallel direction. - pure function shear_stress_tau2(this) result(tau2) - class(shear_stress_t), intent(in) :: this - real(kind=rp) :: tau2(this%msk(0)) - - tau2 = this%tau2_ - end function shear_stress_tau2 - - !> Set the shear stress components. - subroutine shear_stress_set_stress(this, tau1, tau2) - class(shear_stress_t), intent(inout) :: this - real(kind=rp), intent(in) :: tau1(this%msk(0)) - real(kind=rp), intent(in) :: tau2(this%msk(0)) - - call copy(this%tau1_, tau1, this%msk(0)) - call copy(this%tau2_, tau2, this%msk(0)) - - end subroutine shear_stress_set_stress - - !> Destructor. - subroutine shear_stress_free(this) - class(shear_stress_t), target, intent(inout) :: this - call this%free_base - call this%dirichlet%free - call this%neumann1%free - call this%neumann2%free + use num_types + use bc, only : bc_t + use, intrinsic :: iso_c_binding, only : c_ptr + use utils, only : neko_error, nonlinear_index + use coefs, only : coef_t + use dirichlet, only : dirichlet_t + use neumann, only : neumann_t + use math, only : copy + use device_math, only : device_copy + use json_module, only : json_file + use json_utils, only : json_get + implicit none + private + + !> A shear stress boundary condition. + !! @warning Currently hard coded to prescribe a given stress in x, 0 stress + !! in z, and 0 velocity in y. + type, public, extends(bc_t) :: shear_stress_t + !> The stress in the 1st wall-parallel direction. + real(kind=rp), allocatable, private :: tau1_(:) + !> The stress in the 2nd wall-parallel direction. + real(kind=rp), allocatable, private :: tau2_(:) + !> Dirchlet bc for the wall-normal component. + type(dirichlet_t) :: dirichlet + !> Neumann condition for the 1st wall-parallel direction. + type(neumann_t) :: neumann1 + !> Neumann condition for the 2nd wall-parallel direction. + type(neumann_t) :: neumann2 + contains + procedure, pass(this) :: apply_scalar => shear_stress_apply_scalar + procedure, pass(this) :: apply_vector => shear_stress_apply_vector + procedure, pass(this) :: apply_scalar_dev => shear_stress_apply_scalar_dev + procedure, pass(this) :: apply_vector_dev => shear_stress_apply_vector_dev + procedure, pass(this) :: tau1 => shear_stress_tau1 + procedure, pass(this) :: tau2 => shear_stress_tau2 + procedure, pass(this) :: set_stress => shear_stress_set_stress + procedure, pass(this) :: free => shear_stress_free + !> Constructor + procedure, pass(this) :: init => shear_stress_init + !> Finalize. + procedure, pass(this) :: finalize => shear_stress_finalize + end type shear_stress_t + +contains + + !> Apply shear stress for a scalar field @a x. + subroutine shear_stress_apply_scalar(this, x, n, t, tstep) + class(shear_stress_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + integer :: i, m, k, facet + ! Store non-linear index + integer :: idx(4) + + call neko_error("The shear stress bc is not applicable to scalar fields.") + + end subroutine shear_stress_apply_scalar + + !> Boundary condition apply for a generic shear_stress condition + !! to vectors @a x, @a y and @a z + subroutine shear_stress_apply_vector(this, x, y, z, n, t, tstep) + class(shear_stress_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(inout), dimension(n) :: y + real(kind=rp), intent(inout), dimension(n) :: z + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + + call this%neumann1%apply_scalar(x, n, t, tstep) + call this%neumann2%apply_scalar(z, n, t, tstep) + call this%dirichlet%apply_scalar(y, n, t, tstep) + + end subroutine shear_stress_apply_vector + + !> Boundary condition apply for a generic shear_stress condition + !! to a vector @a x (device version) + subroutine shear_stress_apply_scalar_dev(this, x_d, t, tstep) + class(shear_stress_t), intent(inout), target :: this + type(c_ptr) :: x_d + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + + call neko_error("shear_stress bc not implemented on the device") + + end subroutine shear_stress_apply_scalar_dev + + !> Boundary condition apply for a generic shear_stress condition + !! to vectors @a x, @a y and @a z (device version) + subroutine shear_stress_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) + class(shear_stress_t), intent(inout), target :: this + type(c_ptr) :: x_d + type(c_ptr) :: y_d + type(c_ptr) :: z_d + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + + call neko_error("shear_stress bc not implemented on the device") + + end subroutine shear_stress_apply_vector_dev + + !> Constructor. + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine shear_stress_init(this, coef, json) + class(shear_stress_t), target, intent(inout) :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) :: json + real(kind=rp), allocatable :: tau(:) + + call this%init_base(coef) + + call json_get(json, "value", tau) + + call this%dirichlet%init_from_components(coef, 0.0_rp) + call this%neumann1%init_from_components(coef, tau(1)) + call this%neumann2%init_from_components(coef, tau(2)) + end subroutine shear_stress_init + + !> Finalize by allocating the stress arrays and marking the facets for + !! the bc components. + subroutine shear_stress_finalize(this) + class(shear_stress_t), target, intent(inout) :: this + + call this%neumann1%mark_facets(this%marked_facet) + call this%neumann2%mark_facets(this%marked_facet) + call this%dirichlet%mark_facets(this%marked_facet) + + call this%dirichlet%finalize() + call this%neumann1%finalize() + call this%neumann2%finalize() + + end subroutine shear_stress_finalize + + !> Get the stress in the 1st wall-parallel direction. + pure function shear_stress_tau1(this) result(tau1) + class(shear_stress_t), intent(in) :: this + real(kind=rp) :: tau1(this%msk(0)) + + tau1 = this%tau1_ + end function shear_stress_tau1 + + !> Get the stress in the 2nd wall-parallel direction. + pure function shear_stress_tau2(this) result(tau2) + class(shear_stress_t), intent(in) :: this + real(kind=rp) :: tau2(this%msk(0)) + + tau2 = this%tau2_ + end function shear_stress_tau2 + + !> Set the shear stress components. + subroutine shear_stress_set_stress(this, tau1, tau2) + class(shear_stress_t), intent(inout) :: this + real(kind=rp), intent(in) :: tau1(this%msk(0)) + real(kind=rp), intent(in) :: tau2(this%msk(0)) + + call copy(this%tau1_, tau1, this%msk(0)) + call copy(this%tau2_, tau2, this%msk(0)) + + end subroutine shear_stress_set_stress + + !> Destructor. + subroutine shear_stress_free(this) + class(shear_stress_t), target, intent(inout) :: this + call this%free_base + call this%dirichlet%free + call this%neumann1%free + call this%neumann2%free - end subroutine shear_stress_free - end module shear_stress + end subroutine shear_stress_free + +end module shear_stress diff --git a/src/bc/symmetry.f90 b/src/bc/symmetry.f90 index 9d7b9fcff4f..a27efcea55c 100644 --- a/src/bc/symmetry.f90 +++ b/src/bc/symmetry.f90 @@ -42,44 +42,64 @@ module symmetry use stack use tuple use coefs, only : coef_t + use json_module, only : json_file + use zero_dirichlet, only : zero_dirichlet_t use, intrinsic :: iso_c_binding, only : c_ptr implicit none private - !> Mixed Dirichlet-Neumann symmetry plane condition + !> Mixed Dirichlet-Neumann symmetry plane condition. + !! @warning Only works for axis-aligned plane boundaries. type, public, extends(bc_t) :: symmetry_t - type(dirichlet_t) :: bc_x - type(dirichlet_t) :: bc_y - type(dirichlet_t) :: bc_z + type(zero_dirichlet_t) :: bc_x + type(zero_dirichlet_t) :: bc_y + type(zero_dirichlet_t) :: bc_z contains - procedure, pass(this) :: init => symmetry_init procedure, pass(this) :: apply_scalar => symmetry_apply_scalar procedure, pass(this) :: apply_vector => symmetry_apply_vector procedure, pass(this) :: apply_scalar_dev => symmetry_apply_scalar_dev procedure, pass(this) :: apply_vector_dev => symmetry_apply_vector_dev + !> Constructor + procedure, pass(this) :: init => symmetry_init !> Destructor. procedure, pass(this) :: free => symmetry_free + !> Get the axis coressponding to the direction of the normal. + procedure, pass(this) :: get_normal_axis => symmetry_get_normal_axis + !> Finalize. + procedure, pass(this) :: finalize => symmetry_finalize end type symmetry_t contains - - !> Initialize symmetry mask for each axis - subroutine symmetry_init(this, coef) - class(symmetry_t), intent(inout) :: this + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine symmetry_init(this, coef, json) + class(symmetry_t), intent(inout), target :: this type(coef_t), intent(in) :: coef - integer :: i, j, l + type(json_file), intent(inout) ::json + + call this%free() + + call this%init_base(coef) + call this%bc_x%init(this%coef, json) + call this%bc_y%init(this%coef, json) + call this%bc_z%init(this%coef, json) + + end subroutine symmetry_init + + !> Finalize. + !! Marks the appropriate faces for application of a homogeneous Dirchlet + !! condition based on the direction of the axis normal. + subroutine symmetry_finalize(this) + class(symmetry_t), target, intent(inout) :: this + integer :: i, m, j, l type(tuple_i4_t), pointer :: bfp(:) real(kind=rp) :: sx,sy,sz real(kind=rp), parameter :: TOL = 1d-3 type(tuple_i4_t) :: bc_facet integer :: facet, el - call this%free() - - call this%init_base(coef) - call this%bc_x%init_base(this%coef) - call this%bc_y%init_base(this%coef) - call this%bc_z%init_base(this%coef) + call this%finalize_base() associate(c=>this%coef, nx => this%coef%nx, ny => this%coef%ny, & nz => this%coef%nz) @@ -88,38 +108,7 @@ subroutine symmetry_init(this, coef) bc_facet = bfp(i) facet = bc_facet%x(1) el = bc_facet%x(2) - sx = 0d0 - sy = 0d0 - sz = 0d0 - select case (facet) - case(1,2) - do l = 2, c%Xh%lx - 1 - do j = 2, c%Xh%lx -1 - sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) - sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) - sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) - end do - end do - case(3,4) - do l = 2, c%Xh%lx - 1 - do j = 2, c%Xh%lx - 1 - sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) - sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) - sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) - end do - end do - case(5,6) - do l = 2, c%Xh%lx - 1 - do j = 2, c%Xh%lx - 1 - sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) - sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) - sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) - end do - end do - end select - sx = sx / (c%Xh%lx - 2)**2 - sy = sy / (c%Xh%lx - 2)**2 - sz = sz / (c%Xh%lx - 2)**2 + call this%get_normal_axis(sx, sy, sz, facet, el) if (sx .lt. TOL) then call this%bc_x%mark_facet(facet, el) @@ -135,13 +124,56 @@ subroutine symmetry_init(this, coef) end do end associate call this%bc_x%finalize() - call this%bc_x%set_g(0.0_rp) call this%bc_y%finalize() - call this%bc_y%set_g(0.0_rp) call this%bc_z%finalize() - call this%bc_z%set_g(0.0_rp) + end subroutine symmetry_finalize - end subroutine symmetry_init + subroutine symmetry_get_normal_axis(this, sx, sy, sz, facet, el) + class(symmetry_t), target, intent(inout) :: this + real(kind=rp), intent(out) :: sx, sy, sz + integer, intent(in) :: facet + integer, intent(in) :: el + integer :: i, m, j, l + type(tuple_i4_t), pointer :: bfp(:) + real(kind=rp), parameter :: TOL = 1d-3 + type(tuple_i4_t) :: bc_facet + + associate(c=>this%coef, nx => this%coef%nx, ny => this%coef%ny, & + nz => this%coef%nz) + sx = 0d0 + sy = 0d0 + sz = 0d0 + select case (facet) + case(1,2) + do l = 2, c%Xh%lx - 1 + do j = 2, c%Xh%lx -1 + sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) + sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) + sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) + end do + end do + case(3,4) + do l = 2, c%Xh%lx - 1 + do j = 2, c%Xh%lx - 1 + sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) + sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) + sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) + end do + end do + case(5,6) + do l = 2, c%Xh%lx - 1 + do j = 2, c%Xh%lx - 1 + sx = sx + abs(abs(nx(l, j, facet, el)) - 1d0) + sy = sy + abs(abs(ny(l, j, facet, el)) - 1d0) + sz = sz + abs(abs(nz(l, j, facet, el)) - 1d0) + end do + end do + end select + sx = sx / (c%Xh%lx - 2)**2 + sy = sy / (c%Xh%lx - 2)**2 + sz = sz / (c%Xh%lx - 2)**2 + end associate + end subroutine symmetry_get_normal_axis !> No-op scalar apply subroutine symmetry_apply_scalar(this, x, n, t, tstep) @@ -162,9 +194,9 @@ subroutine symmetry_apply_vector(this, x, y, z, n, t, tstep) real(kind=rp), intent(in), optional :: t integer, intent(in), optional :: tstep - call this%bc_x%apply_scalar(x,n) - call this%bc_y%apply_scalar(y,n) - call this%bc_z%apply_scalar(z,n) + call this%bc_x%apply_scalar(x, n) + call this%bc_y%apply_scalar(y, n) + call this%bc_z%apply_scalar(z, n) end subroutine symmetry_apply_vector diff --git a/src/bc/usr_inflow.f90 b/src/bc/usr_inflow.f90 index 3e78d902b2a..c0c00353668 100644 --- a/src/bc/usr_inflow.f90 +++ b/src/bc/usr_inflow.f90 @@ -39,6 +39,7 @@ module usr_inflow use device_inhom_dirichlet use utils, only : neko_error, nonlinear_index, neko_warning use bc, only : bc_t + use json_module, only : json_file implicit none private @@ -57,8 +58,12 @@ module usr_inflow procedure, pass(this) :: set_eval => usr_inflow_set_eval procedure, pass(this) :: apply_vector_dev => usr_inflow_apply_vector_dev procedure, pass(this) :: apply_scalar_dev => usr_inflow_apply_scalar_dev + !> Constructor + procedure, pass(this) :: init => usr_inflow_init !> Destructor. procedure, pass(this) :: free => usr_inflow_free + !> Finalize. + procedure, pass(this) :: finalize => usr_inflow_finalize end type usr_inflow_t abstract interface @@ -104,6 +109,17 @@ end subroutine usr_inflow_eval contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine usr_inflow_init(this, coef, json) + class(usr_inflow_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + + call this%init_base(coef) + end subroutine usr_inflow_init + subroutine usr_inflow_free(this) class(usr_inflow_t), target, intent(inout) :: this @@ -346,4 +362,10 @@ subroutine usr_inflow_validate(this) end subroutine usr_inflow_validate + !> Finalize + subroutine usr_inflow_finalize(this) + class(usr_inflow_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine usr_inflow_finalize end module usr_inflow diff --git a/src/bc/usr_scalar.f90 b/src/bc/usr_scalar.f90 index b3f2cbd481f..058eff65259 100644 --- a/src/bc/usr_scalar.f90 +++ b/src/bc/usr_scalar.f90 @@ -38,6 +38,7 @@ module usr_scalar use device use device_inhom_dirichlet use utils, only : neko_error, nonlinear_index, neko_warning + use json_module, only : json_file implicit none private @@ -52,8 +53,12 @@ module usr_scalar procedure, pass(this) :: set_eval => usr_scalar_set_eval procedure, pass(this) :: apply_vector_dev => usr_scalar_apply_vector_dev procedure, pass(this) :: apply_scalar_dev => usr_scalar_apply_scalar_dev + !> Constructor + procedure, pass(this) :: init => usr_scalar_init !> Destructor. procedure, pass(this) :: free => usr_scalar_free + !> Finalize. + procedure, pass(this) :: finalize => usr_scalar_finalize end type usr_scalar_t abstract interface @@ -98,6 +103,17 @@ end subroutine usr_scalar_bc_eval contains + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine usr_scalar_init(this, coef, json) + class(usr_scalar_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + + call this%init_base(coef) + end subroutine usr_scalar_init + subroutine usr_scalar_free(this) class(usr_scalar_t), target, intent(inout) :: this @@ -322,4 +338,11 @@ subroutine usr_scalar_validate(this) end subroutine usr_scalar_validate + !> Finalize + subroutine usr_scalar_finalize(this) + class(usr_scalar_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine usr_scalar_finalize + end module usr_scalar diff --git a/src/bc/wall_model_bc.f90 b/src/bc/wall_model_bc.f90 index 8e3fc214783..f921c6831dd 100644 --- a/src/bc/wall_model_bc.f90 +++ b/src/bc/wall_model_bc.f90 @@ -32,113 +32,114 @@ ! !> Defines the `wall_model_bc_t` type. module wall_model_bc - use num_types - use bc, only : bc_t - use, intrinsic :: iso_c_binding, only : c_ptr - use utils, only : neko_error, nonlinear_index - use coefs, only : coef_t - use wall_model, only : wall_model_t - use rough_log_law, only : rough_log_law_t - use spalding, only : spalding_t - use shear_stress, only : shear_stress_t - implicit none - private - - !> A shear stress boundary condition, computing the stress values using a - !! wall model. - type, public, extends(shear_stress_t) :: wall_model_bc_t - !> The wall model to compute the stress. - class(wall_model_t), allocatable :: wall_model - contains - procedure, pass(this) :: apply_scalar => wall_model_bc_apply_scalar - procedure, pass(this) :: apply_vector => wall_model_bc_apply_vector - procedure, pass(this) :: apply_scalar_dev => wall_model_bc_apply_scalar_dev - procedure, pass(this) :: apply_vector_dev => wall_model_bc_apply_vector_dev - procedure, pass(this) :: init_wall_model_bc => & - wall_model_bc_init_wall_model_bc - end type wall_model_bc_t - - contains - - !> Apply shear stress for a scalar field @a x. - subroutine wall_model_bc_apply_scalar(this, x, n, t, tstep) - class(wall_model_bc_t), intent(inout) :: this - integer, intent(in) :: n - real(kind=rp), intent(inout), dimension(n) :: x - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - - call neko_error("The wall model bc is not applicable to scalar fields.") - - end subroutine wall_model_bc_apply_scalar - - !> Boundary condition apply for a generic wall_model_bc condition - !! to vectors @a x, @a y and @a z - subroutine wall_model_bc_apply_vector(this, x, y, z, n, t, tstep) - class(wall_model_bc_t), intent(inout) :: this - integer, intent(in) :: n - real(kind=rp), intent(inout), dimension(n) :: x - real(kind=rp), intent(inout), dimension(n) :: y - real(kind=rp), intent(inout), dimension(n) :: z - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - integer :: i, m, k, fid - real(kind=rp) :: magtau - - call this%wall_model%compute(t, tstep) - - do i=1, this%msk(0) - magtau = sqrt(this%wall_model%tau_x(i)**2 + this%wall_model%tau_y(i)**2& - + this%wall_model%tau_z(i)**2) - - ! Mark sampling nodes with a -1 for debugging - this%wall_model%tau_field%x(this%wall_model%ind_r(i), & - this%wall_model%ind_s(i), & - this%wall_model%ind_t(i), & - this%wall_model%ind_e(i)) = -1.0_rp - this%wall_model%tau_field%x(this%msk(i),1,1,1) = magtau - end do - - call this%shear_stress_t%set_stress(this%wall_model%tau_x, & - this%wall_model%tau_z) - call this%shear_stress_t%apply_vector(x, y, z, n, t, tstep) - - end subroutine wall_model_bc_apply_vector - - !> Boundary condition apply for a generic wall_model_bc condition - !! to a vector @a x (device version) - subroutine wall_model_bc_apply_scalar_dev(this, x_d, t, tstep) - class(wall_model_bc_t), intent(inout), target :: this - type(c_ptr) :: x_d - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - - call neko_error("wall_model_bc bc not implemented on the device") - - end subroutine wall_model_bc_apply_scalar_dev - - !> Boundary condition apply for a generic wall_model_bc condition - !! to vectors @a x, @a y and @a z (device version) - subroutine wall_model_bc_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) - class(wall_model_bc_t), intent(inout), target :: this - type(c_ptr) :: x_d - type(c_ptr) :: y_d - type(c_ptr) :: z_d - real(kind=rp), intent(in), optional :: t - integer, intent(in), optional :: tstep - - call neko_error("wall_model_bc bc not implemented on the device") - - end subroutine wall_model_bc_apply_vector_dev - - !> Constructor. - !> @param coef The SEM coefficients. - subroutine wall_model_bc_init_wall_model_bc(this, coef) - class(wall_model_bc_t), intent(inout) :: this - type(coef_t), target, intent(in) :: coef - - call this%shear_stress_t%init_shear_stress(coef) - - end subroutine wall_model_bc_init_wall_model_bc - - end module wall_model_bc \ No newline at end of file + use num_types + use bc, only : bc_t + use, intrinsic :: iso_c_binding, only : c_ptr + use utils, only : neko_error, nonlinear_index + use coefs, only : coef_t + use wall_model, only : wall_model_t + use rough_log_law, only : rough_log_law_t + use spalding, only : spalding_t + use shear_stress, only : shear_stress_t + use json_module, only : json_file + implicit none + private + + !> A shear stress boundary condition, computing the stress values using a + !! wall model. + type, public, extends(shear_stress_t) :: wall_model_bc_t + !> The wall model to compute the stress. + class(wall_model_t), allocatable :: wall_model + contains + procedure, pass(this) :: apply_scalar => wall_model_bc_apply_scalar + procedure, pass(this) :: apply_vector => wall_model_bc_apply_vector + procedure, pass(this) :: apply_scalar_dev => wall_model_bc_apply_scalar_dev + procedure, pass(this) :: apply_vector_dev => wall_model_bc_apply_vector_dev + procedure, pass(this) :: init => wall_model_bc_init + end type wall_model_bc_t + +contains + + !> Apply shear stress for a scalar field @a x. + subroutine wall_model_bc_apply_scalar(this, x, n, t, tstep) + class(wall_model_bc_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + + call neko_error("The wall model bc is not applicable to scalar fields.") + + end subroutine wall_model_bc_apply_scalar + + !> Boundary condition apply for a generic wall_model_bc condition + !! to vectors @a x, @a y and @a z + subroutine wall_model_bc_apply_vector(this, x, y, z, n, t, tstep) + class(wall_model_bc_t), intent(inout) :: this + integer, intent(in) :: n + real(kind=rp), intent(inout), dimension(n) :: x + real(kind=rp), intent(inout), dimension(n) :: y + real(kind=rp), intent(inout), dimension(n) :: z + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + integer :: i, m, k, fid + real(kind=rp) :: magtau + + call this%wall_model%compute(t, tstep) + + do i=1, this%msk(0) + magtau = sqrt(this%wall_model%tau_x(i)**2 + this%wall_model%tau_y(i)**2& + + this%wall_model%tau_z(i)**2) + + ! Mark sampling nodes with a -1 for debugging + this%wall_model%tau_field%x(this%wall_model%ind_r(i), & + this%wall_model%ind_s(i), & + this%wall_model%ind_t(i), & + this%wall_model%ind_e(i)) = -1.0_rp + this%wall_model%tau_field%x(this%msk(i),1,1,1) = magtau + end do + + call this%shear_stress_t%set_stress(this%wall_model%tau_x, & + this%wall_model%tau_z) + call this%shear_stress_t%apply_vector(x, y, z, n, t, tstep) + + end subroutine wall_model_bc_apply_vector + + !> Boundary condition apply for a generic wall_model_bc condition + !! to a vector @a x (device version) + subroutine wall_model_bc_apply_scalar_dev(this, x_d, t, tstep) + class(wall_model_bc_t), intent(inout), target :: this + type(c_ptr) :: x_d + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + + call neko_error("wall_model_bc bc not implemented on the device") + + end subroutine wall_model_bc_apply_scalar_dev + + !> Boundary condition apply for a generic wall_model_bc condition + !! to vectors @a x, @a y and @a z (device version) + subroutine wall_model_bc_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) + class(wall_model_bc_t), intent(inout), target :: this + type(c_ptr) :: x_d + type(c_ptr) :: y_d + type(c_ptr) :: z_d + real(kind=rp), intent(in), optional :: t + integer, intent(in), optional :: tstep + + call neko_error("wall_model_bc bc not implemented on the device") + + end subroutine wall_model_bc_apply_vector_dev + + !> Constructor. + !> @param coef The SEM coefficients. + subroutine wall_model_bc_init(this, coef, json) + class(wall_model_bc_t), target, intent(inout) :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) :: json + + call this%shear_stress_t%init(coef, json) + + end subroutine wall_model_bc_init + +end module wall_model_bc \ No newline at end of file diff --git a/src/bc/wall.f90 b/src/bc/zero_dirichlet.f90 similarity index 50% rename from src/bc/wall.f90 rename to src/bc/zero_dirichlet.f90 index d9eb57bb213..27707abc99d 100644 --- a/src/bc/wall.f90 +++ b/src/bc/zero_dirichlet.f90 @@ -1,4 +1,4 @@ -! Copyright (c) 2020-2021, The Neko Authors +! Copyright (c) 2020-2024, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -30,32 +30,51 @@ ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! -!> Defines wall boundary conditions -module wall - use device_wall +!> Defines a zero-valued Dirichlet boundary condition. +module zero_dirichlet + use device_zero_dirichlet use num_types, only : rp use bc, only : bc_t use, intrinsic :: iso_c_binding, only : c_ptr + use coefs, only : coef_t + use json_module, only : json_file implicit none private - !> No-slip Wall boundary condition - type, public, extends(bc_t) :: no_slip_wall_t + !> Zero-valued Dirichlet boundary condition. + !! Used for no-slip walls, but also for various auxillary conditions, + !! such as for residuals. + type, public, extends(bc_t) :: zero_dirichlet_t contains - procedure, pass(this) :: apply_scalar => no_slip_wall_apply_scalar - procedure, pass(this) :: apply_vector => no_slip_wall_apply_vector - procedure, pass(this) :: apply_scalar_dev => no_slip_wall_apply_scalar_dev - procedure, pass(this) :: apply_vector_dev => no_slip_wall_apply_vector_dev + procedure, pass(this) :: apply_scalar => zero_dirichlet_apply_scalar + procedure, pass(this) :: apply_vector => zero_dirichlet_apply_vector + procedure, pass(this) :: apply_scalar_dev => zero_dirichlet_apply_scalar_dev + procedure, pass(this) :: apply_vector_dev => zero_dirichlet_apply_vector_dev + !> Constructor. + procedure, pass(this) :: init => zero_dirichlet_init !> Destructor. - procedure, pass(this) :: free => no_slip_wall_free - end type no_slip_wall_t + procedure, pass(this) :: free => zero_dirichlet_free + !> Finalize. + procedure, pass(this) :: finalize => zero_dirichlet_finalize + end type zero_dirichlet_t contains - !> Boundary condition apply for a no-slip wall condition + !> Constructor + !! @param[in] coef The SEM coefficients. + !! @param[inout] json The JSON object configuring the boundary condition. + subroutine zero_dirichlet_init(this, coef, json) + class(zero_dirichlet_t), intent(inout), target :: this + type(coef_t), intent(in) :: coef + type(json_file), intent(inout) ::json + + call this%init_base(coef) + end subroutine zero_dirichlet_init + + !> Apply boundary condition to a scalar field. !! to a vector @a x - subroutine no_slip_wall_apply_scalar(this, x, n, t, tstep) - class(no_slip_wall_t), intent(inout) :: this + subroutine zero_dirichlet_apply_scalar(this, x, n, t, tstep) + class(zero_dirichlet_t), intent(inout) :: this integer, intent(in) :: n real(kind=rp), intent(inout), dimension(n) :: x real(kind=rp), intent(in), optional :: t @@ -68,12 +87,11 @@ subroutine no_slip_wall_apply_scalar(this, x, n, t, tstep) x(k) = 0d0 end do - end subroutine no_slip_wall_apply_scalar + end subroutine zero_dirichlet_apply_scalar - !> Boundary condition apply for a no-slip wall condition - !! to vectors @a x, @a y and @a z - subroutine no_slip_wall_apply_vector(this, x, y, z, n, t, tstep) - class(no_slip_wall_t), intent(inout) :: this + !> Apply boundary condition to a vector field. + subroutine zero_dirichlet_apply_vector(this, x, y, z, n, t, tstep) + class(zero_dirichlet_t), intent(inout) :: this integer, intent(in) :: n real(kind=rp), intent(inout), dimension(n) :: x real(kind=rp), intent(inout), dimension(n) :: y @@ -90,41 +108,46 @@ subroutine no_slip_wall_apply_vector(this, x, y, z, n, t, tstep) z(k) = 0d0 end do - end subroutine no_slip_wall_apply_vector + end subroutine zero_dirichlet_apply_vector - !> Boundary condition apply for a no-slip wall condition - !! to a vector @a x (device version) - subroutine no_slip_wall_apply_scalar_dev(this, x_d, t, tstep) - class(no_slip_wall_t), intent(inout), target :: this + !> Apply boundary condition to a scalar field, device version. + subroutine zero_dirichlet_apply_scalar_dev(this, x_d, t, tstep) + class(zero_dirichlet_t), intent(inout), target :: this type(c_ptr) :: x_d real(kind=rp), intent(in), optional :: t integer, intent(in), optional :: tstep - call device_no_slip_wall_apply_scalar(this%msk_d, x_d, size(this%msk)) + call device_zero_dirichlet_apply_scalar(this%msk_d, x_d, size(this%msk)) - end subroutine no_slip_wall_apply_scalar_dev + end subroutine zero_dirichlet_apply_scalar_dev - !> Boundary condition apply for a no-slip wall condition - !! to vectors @a x, @a y and @a z (device version) - subroutine no_slip_wall_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) - class(no_slip_wall_t), intent(inout), target :: this + !> Apply boundary condition to a vector field, device version. + subroutine zero_dirichlet_apply_vector_dev(this, x_d, y_d, z_d, t, tstep) + class(zero_dirichlet_t), intent(inout), target :: this type(c_ptr) :: x_d type(c_ptr) :: y_d type(c_ptr) :: z_d real(kind=rp), intent(in), optional :: t integer, intent(in), optional :: tstep - call device_no_slip_wall_apply_vector(this%msk_d, x_d, y_d, z_d, & + call device_zero_dirichlet_apply_vector(this%msk_d, x_d, y_d, z_d, & size(this%msk)) - end subroutine no_slip_wall_apply_vector_dev + end subroutine zero_dirichlet_apply_vector_dev !> Destructor - subroutine no_slip_wall_free(this) - class(no_slip_wall_t), target, intent(inout) :: this + subroutine zero_dirichlet_free(this) + class(zero_dirichlet_t), target, intent(inout) :: this call this%free_base() - end subroutine no_slip_wall_free + end subroutine zero_dirichlet_free + + !> Finalize + subroutine zero_dirichlet_finalize(this) + class(zero_dirichlet_t), target, intent(inout) :: this + + call this%finalize_base() + end subroutine zero_dirichlet_finalize -end module wall +end module zero_dirichlet diff --git a/src/case.f90 b/src/case.f90 index 2ad77dd71da..f23caab762c 100644 --- a/src/case.f90 +++ b/src/case.f90 @@ -264,11 +264,6 @@ subroutine case_init_common(C) end if end if - ! Setup user boundary conditions for the scalar. - if (scalar) then - call C%scalar%set_user_bc(C%usr%scalar_user_bc) - end if - ! ! Setup initial conditions ! diff --git a/src/common/projection.f90 b/src/common/projection.f90 index 31bd8b706cd..c2069f44fa7 100644 --- a/src/common/projection.f90 +++ b/src/common/projection.f90 @@ -65,7 +65,7 @@ module projection use math use coefs, only : coef_t use ax_product, only : ax_t - use bc, only : bc_list_t, bc_list_apply_scalar + use bc_list, only : bc_list_t use comm use gather_scatter, only : gs_t, GS_OP_ADD use neko_config, only : NEKO_BCKND_DEVICE @@ -75,6 +75,7 @@ module projection use profiler, only : profiler_start_region, profiler_end_region use logger, only : LOG_SIZE, neko_log use, intrinsic :: iso_c_binding + use bc_list, only : bc_list_t use time_step_controller, only : time_step_controller_t implicit none @@ -317,7 +318,7 @@ subroutine bcknd_project_back(this,x,Ax,coef, bclst, gs_h, n) call Ax%compute(this%bb(1,this%m), x, coef, coef%msh, coef%Xh) call gs_h%gs_op_vector(this%bb(1,this%m), n, GS_OP_ADD) - call bc_list_apply_scalar(bclst, this%bb(1,this%m), n) + call bclst%apply_scalar(this%bb(1, this%m), n) if (NEKO_BCKND_DEVICE .eq. 1) then call device_proj_ortho(this, this%xx_d, this%bb_d, coef%mult_d, n) diff --git a/src/common/user_intf.f90 b/src/common/user_intf.f90 index a8193208e70..4aefdc270b6 100644 --- a/src/common/user_intf.f90 +++ b/src/common/user_intf.f90 @@ -39,7 +39,7 @@ module user_intf use scalar_user_source_term, only : scalar_user_source_term_t, & scalar_source_compute_pointwise, scalar_source_compute_vector use coefs, only : coef_t - use bc, only: bc_list_t + use bc_list, only: bc_list_t use mesh, only : mesh_t use usr_inflow, only : usr_inflow_t, usr_inflow_eval use usr_scalar, only : usr_scalar_t, usr_scalar_bc_eval diff --git a/src/fluid/fluid_pnpn.f90 b/src/fluid/fluid_pnpn.f90 index 9aa9a3760ea..85ec97148fd 100644 --- a/src/fluid/fluid_pnpn.f90 +++ b/src/fluid/fluid_pnpn.f90 @@ -69,8 +69,8 @@ module fluid_pnpn use neko_config, only : NEKO_BCKND_DEVICE use math, only : col2, add2 use mathops, only : opadd2cm, opcolv - use bc, only: bc_list_t, bc_list_init, bc_list_add, bc_list_free, & - bc_list_apply_scalar, bc_list_apply_vector + use bc_list, only: bc_list_t + use zero_dirichlet, only : zero_dirichlet_t implicit none private @@ -89,11 +89,11 @@ module fluid_pnpn type(facet_normal_t) :: bc_prs_surface !< Surface term in pressure rhs type(facet_normal_t) :: bc_sym_surface !< Surface term in pressure rhs - type(dirichlet_t) :: bc_vel_res !< Dirichlet condition vel. res. - type(dirichlet_t) :: bc_field_dirichlet_p !< Dirichlet condition vel. res. - type(dirichlet_t) :: bc_field_dirichlet_u !< Dirichlet condition vel. res. - type(dirichlet_t) :: bc_field_dirichlet_v !< Dirichlet condition vel. res. - type(dirichlet_t) :: bc_field_dirichlet_w !< Dirichlet condition vel. res. + type(zero_dirichlet_t) :: bc_vel_res !< Dirichlet condition vel. res. + type(zero_dirichlet_t) :: bc_field_dirichlet_p !< Dirichlet condition vel. res. + type(zero_dirichlet_t) :: bc_field_dirichlet_u !< Dirichlet condition vel. res. + type(zero_dirichlet_t) :: bc_field_dirichlet_v !< Dirichlet condition vel. res. + type(zero_dirichlet_t) :: bc_field_dirichlet_w !< Dirichlet condition vel. res. type(non_normal_t) :: bc_vel_res_non_normal !< Dirichlet condition vel. res. type(bc_list_t) :: bclst_vel_res type(bc_list_t) :: bclst_du @@ -196,100 +196,80 @@ subroutine fluid_pnpn_init(this, msh, lx, params, user, material_properties) end associate ! Initialize velocity surface terms in pressure rhs - call this%bc_prs_surface%init_base(this%c_Xh) - call this%bc_prs_surface%mark_zone(msh%inlet) - call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,& - 'v', this%bc_labels) + call this%bc_prs_surface%init(this%c_Xh, params) + call this%bc_prs_surface%mark_zones_from_list('v', this%bc_labels) !This impacts the rhs of the pressure, need to check what is correct to add here - call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_u', this%bc_labels) - call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_v', this%bc_labels) - call this%bc_prs_surface%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_w', this%bc_labels) + call this%bc_prs_surface%mark_zones_from_list('d_vel_u', this%bc_labels) + call this%bc_prs_surface%mark_zones_from_list('d_vel_v', this%bc_labels) + call this%bc_prs_surface%mark_zones_from_list('d_vel_w', this%bc_labels) call this%bc_prs_surface%finalize() ! Initialize symmetry surface terms in pressure rhs - call this%bc_sym_surface%init_base(this%c_Xh) - call this%bc_sym_surface%mark_zone(msh%sympln) - call this%bc_sym_surface%mark_zones_from_list(msh%labeled_zones,& - 'sym', this%bc_labels) + call this%bc_sym_surface%init(this%c_Xh, params) + call this%bc_sym_surface%mark_zones_from_list('sym', this%bc_labels) ! Same here, should du, dv, dw be marked here? call this%bc_sym_surface%finalize() ! Initialize dirichlet bcs for velocity residual - call this%bc_vel_res_non_normal%init_base(this%c_Xh) - call this%bc_vel_res_non_normal%mark_zone(msh%outlet_normal) - call this%bc_vel_res_non_normal%mark_zones_from_list(msh%labeled_zones,& - 'on', this%bc_labels) - call this%bc_vel_res_non_normal%mark_zones_from_list(msh%labeled_zones,& - 'on+dong', & + call this%bc_vel_res_non_normal%init(this%c_Xh, params) + call this%bc_vel_res_non_normal%mark_zones_from_list('on', this%bc_labels) + call this%bc_vel_res_non_normal%mark_zones_from_list('on+dong', & this%bc_labels) call this%bc_vel_res_non_normal%finalize() - call this%bc_vel_res_non_normal%init(this%c_Xh) - call this%bc_field_dirichlet_p%init_base(this%c_Xh) - call this%bc_field_dirichlet_p%mark_zones_from_list(msh%labeled_zones, 'on+dong', & + call this%bc_field_dirichlet_p%init(this%c_Xh, params) + call this%bc_field_dirichlet_p%mark_zones_from_list('on+dong', & this%bc_labels) - call this%bc_field_dirichlet_p%mark_zones_from_list(msh%labeled_zones, & - 'o+dong', this%bc_labels) - call this%bc_field_dirichlet_p%mark_zones_from_list(msh%labeled_zones, 'd_pres', & + call this%bc_field_dirichlet_p%mark_zones_from_list('o+dong', & + this%bc_labels) + call this%bc_field_dirichlet_p%mark_zones_from_list('d_pres', & this%bc_labels) call this%bc_field_dirichlet_p%finalize() - call this%bc_field_dirichlet_p%set_g(0.0_rp) - call bc_list_init(this%bclst_dp) - call bc_list_add(this%bclst_dp, this%bc_field_dirichlet_p) + call this%bclst_dp%init() + call this%bclst_dp%append(this%bc_field_dirichlet_p) !Add 0 prs bcs - call bc_list_add(this%bclst_dp, this%bc_prs) + call this%bclst_dp%append(this%bc_prs) - call this%bc_field_dirichlet_u%init_base(this%c_Xh) - call this%bc_field_dirichlet_u%mark_zones_from_list(msh%labeled_zones, 'd_vel_u', & + call this%bc_field_dirichlet_u%init(this%c_Xh, params) + call this%bc_field_dirichlet_u%mark_zones_from_list('d_vel_u', & this%bc_labels) call this%bc_field_dirichlet_u%finalize() - call this%bc_field_dirichlet_u%set_g(0.0_rp) - call this%bc_field_dirichlet_v%init_base(this%c_Xh) - call this%bc_field_dirichlet_v%mark_zones_from_list(msh%labeled_zones, 'd_vel_v', & + call this%bc_field_dirichlet_v%init(this%c_Xh, params) + call this%bc_field_dirichlet_v%mark_zones_from_list('d_vel_v', & this%bc_labels) call this%bc_field_dirichlet_v%finalize() - call this%bc_field_dirichlet_v%set_g(0.0_rp) - call this%bc_field_dirichlet_w%init_base(this%c_Xh) - call this%bc_field_dirichlet_w%mark_zones_from_list(msh%labeled_zones, 'd_vel_w', & + call this%bc_field_dirichlet_w%init(this%c_Xh, params) + call this%bc_field_dirichlet_w%mark_zones_from_list('d_vel_w', & this%bc_labels) call this%bc_field_dirichlet_w%finalize() - call this%bc_field_dirichlet_w%set_g(0.0_rp) - - call this%bc_vel_res%init_base(this%c_Xh) - call this%bc_vel_res%mark_zone(msh%inlet) - call this%bc_vel_res%mark_zone(msh%wall) - call this%bc_vel_res%mark_zones_from_list(msh%labeled_zones, & - 'v', this%bc_labels) - call this%bc_vel_res%mark_zones_from_list(msh%labeled_zones, & - 'w', this%bc_labels) + + call this%bc_vel_res%init(this%c_Xh, params) + call this%bc_vel_res%mark_zones_from_list('v', this%bc_labels) + call this%bc_vel_res%mark_zones_from_list('w', this%bc_labels) call this%bc_vel_res%finalize() - call this%bc_vel_res%set_g(0.0_rp) - call bc_list_init(this%bclst_vel_res) - call bc_list_add(this%bclst_vel_res, this%bc_vel_res) - call bc_list_add(this%bclst_vel_res, this%bc_vel_res_non_normal) - call bc_list_add(this%bclst_vel_res, this%bc_sym) + call this%bclst_vel_res%init() + call this%bclst_vel_res%append(this%bc_vel_res) + call this%bclst_vel_res%append(this%bc_vel_res_non_normal) + call this%bclst_vel_res%append(this%bc_sym) !Initialize bcs for u, v, w velocity components - call bc_list_init(this%bclst_du) - call bc_list_add(this%bclst_du,this%bc_sym%bc_x) - call bc_list_add(this%bclst_du,this%bc_vel_res_non_normal%bc_x) - call bc_list_add(this%bclst_du, this%bc_vel_res) - call bc_list_add(this%bclst_du, this%bc_field_dirichlet_u) - - call bc_list_init(this%bclst_dv) - call bc_list_add(this%bclst_dv,this%bc_sym%bc_y) - call bc_list_add(this%bclst_dv,this%bc_vel_res_non_normal%bc_y) - call bc_list_add(this%bclst_dv, this%bc_vel_res) - call bc_list_add(this%bclst_dv, this%bc_field_dirichlet_v) - - call bc_list_init(this%bclst_dw) - call bc_list_add(this%bclst_dw,this%bc_sym%bc_z) - call bc_list_add(this%bclst_dw,this%bc_vel_res_non_normal%bc_z) - call bc_list_add(this%bclst_dw, this%bc_vel_res) - call bc_list_add(this%bclst_dw, this%bc_field_dirichlet_w) + call this%bclst_du%init() + call this%bclst_du%append(this%bc_sym%bc_x) + call this%bclst_du%append(this%bc_vel_res_non_normal%bc_x) + call this%bclst_du%append(this%bc_vel_res) + call this%bclst_du%append(this%bc_field_dirichlet_u) + + call this%bclst_dv%init() + call this%bclst_dv%append(this%bc_sym%bc_y) + call this%bclst_dv%append(this%bc_vel_res_non_normal%bc_y) + call this%bclst_dv%append(this%bc_vel_res) + call this%bclst_dv%append(this%bc_field_dirichlet_v) + + call this%bclst_dw%init() + call this%bclst_dw%append(this%bc_sym%bc_z) + call this%bclst_dw%append(this%bc_vel_res_non_normal%bc_z) + call this%bclst_dw%append(this%bc_vel_res) + call this%bclst_dw%append(this%bc_field_dirichlet_w) !Intialize projection space thingy call this%proj_prs%init(this%dm_Xh%size(), this%pr_projection_dim, & @@ -446,8 +426,8 @@ subroutine fluid_pnpn_free(this) call this%bc_prs_surface%free() call this%bc_sym_surface%free() - call bc_list_free(this%bclst_vel_res) - call bc_list_free(this%bclst_dp) + call this%bclst_vel_res%free() + call this%bclst_dp%free() call this%proj_prs%free() call this%proj_u%free() call this%proj_v%free() @@ -600,7 +580,7 @@ subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) dt, mu, rho) call gs_Xh%op(p_res, GS_OP_ADD) - call bc_list_apply_scalar(this%bclst_dp, p_res%x, p%dof%size(), t, tstep) + call this%bclst_dp%apply_scalar(p_res%x, p%dof%size(), t, tstep) call profiler_end_region call this%proj_prs%pre_solving(p_res%x, tstep, c_Xh, n, dt_controller, 'Pressure') @@ -636,7 +616,7 @@ subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) call gs_Xh%op(v_res, GS_OP_ADD) call gs_Xh%op(w_res, GS_OP_ADD) - call bc_list_apply_vector(this%bclst_vel_res,& + call this%bclst_vel_res%apply_vector(& u_res%x, v_res%x, w_res%x, dm_Xh%size(),& t, tstep) diff --git a/src/fluid/fluid_scheme.f90 b/src/fluid/fluid_scheme.f90 index 343e899102a..0da6853a62f 100644 --- a/src/fluid/fluid_scheme.f90 +++ b/src/fluid/fluid_scheme.f90 @@ -46,8 +46,8 @@ module fluid_scheme use space use dofmap, only : dofmap_t use krylov, only : ksp_t + use zero_dirichlet, only : zero_dirichlet_t use coefs, only: coef_t - use wall, only : no_slip_wall_t use inflow, only : inflow_t use usr_inflow, only : usr_inflow_t, usr_inflow_eval use blasius, only : blasius_t @@ -81,6 +81,7 @@ module fluid_scheme use material_properties, only : material_properties_t use field_series use time_step_controller + use bc_list, only : bc_list_t implicit none !> Base type of all fluid formulations @@ -110,7 +111,7 @@ module fluid_scheme integer :: pr_projection_dim !< Size of the projection space for ksp_pr integer :: vel_projection_activ_step !< Steps to activate projection for ksp_vel integer :: pr_projection_activ_step !< Steps to activate projection for ksp_pr - type(no_slip_wall_t) :: bc_wall !< No-slip wall for velocity + type(zero_dirichlet_t) :: bc_wall !< No-slip wall for velocity class(bc_t), allocatable :: bc_inflow !< Dirichlet inflow for velocity ! Attributes for field dirichlet BCs @@ -322,15 +323,12 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, & this%bc_labels) end if - call bc_list_init(this%bclst_vel) + call this%bclst_vel%init() - call this%bc_sym%init_base(this%c_Xh) - call this%bc_sym%mark_zone(msh%sympln) - call this%bc_sym%mark_zones_from_list(msh%labeled_zones,& - 'sym', this%bc_labels) + call this%bc_sym%init(this%c_Xh, params) + call this%bc_sym%mark_zones_from_list('sym', this%bc_labels) call this%bc_sym%finalize() - call this%bc_sym%init(this%c_Xh) - call bc_list_add(this%bclst_vel, this%bc_sym) + call this%bclst_vel%append(this%bc_sym) ! ! Inflow @@ -347,45 +345,20 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, & call neko_error('Invalid inflow condition '//string_val1) end if - call this%bc_inflow%init_base(this%c_Xh) - call this%bc_inflow%mark_zone(msh%inlet) - call this%bc_inflow%mark_zones_from_list(msh%labeled_zones,& - 'v', this%bc_labels) + call this%bc_inflow%init(this%c_Xh, params) + call this%bc_inflow%mark_zones_from_list('v', this%bc_labels) call this%bc_inflow%finalize() - call bc_list_add(this%bclst_vel, this%bc_inflow) - - if (trim(string_val1) .eq. "uniform") then - call json_get(params, 'case.fluid.inflow_condition.value', real_vec) - select type(bc_if => this%bc_inflow) - type is(inflow_t) - call bc_if%set_inflow(real_vec) - end select - else if (trim(string_val1) .eq. "blasius") then - select type(bc_if => this%bc_inflow) - type is(blasius_t) - call json_get(params, 'case.fluid.blasius.delta', real_val) - call json_get(params, 'case.fluid.blasius.approximation',& - string_val2) - call json_get(params, 'case.fluid.blasius.freestream_velocity',& - real_vec) - - call bc_if%set_params(real_vec, real_val, string_val2) - - end select - else if (trim(string_val1) .eq. "user") then - end if + call this%bclst_vel%append(this%bc_inflow) end if - call this%bc_wall%init_base(this%c_Xh) - call this%bc_wall%mark_zone(msh%wall) - call this%bc_wall%mark_zones_from_list(msh%labeled_zones,& - 'w', this%bc_labels) + call this%bc_wall%init(this%c_Xh, params) + call this%bc_wall%mark_zones_from_list('w', this%bc_labels) call this%bc_wall%finalize() - call bc_list_add(this%bclst_vel, this%bc_wall) + call this%bclst_vel%append(this%bc_wall) ! Setup field dirichlet bc for u-velocity call this%user_field_bc_vel%bc_u%init_base(this%c_Xh) - call this%user_field_bc_vel%bc_u%mark_zones_from_list(msh%labeled_zones,& + call this%user_field_bc_vel%bc_u%mark_zones_from_list(& 'd_vel_u', this%bc_labels) call this%user_field_bc_vel%bc_u%finalize() @@ -395,7 +368,7 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, & ! Setup field dirichlet bc for v-velocity call this%user_field_bc_vel%bc_v%init_base(this%c_Xh) - call this%user_field_bc_vel%bc_v%mark_zones_from_list(msh%labeled_zones,& + call this%user_field_bc_vel%bc_v%mark_zones_from_list(& 'd_vel_v', this%bc_labels) call this%user_field_bc_vel%bc_v%finalize() @@ -405,8 +378,8 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, & ! Setup field dirichlet bc for w-velocity call this%user_field_bc_vel%bc_w%init_base(this%c_Xh) - call this%user_field_bc_vel%bc_w%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_w', this%bc_labels) + call this%user_field_bc_vel%bc_w%mark_zones_from_list('d_vel_w',& + this%bc_labels) call this%user_field_bc_vel%bc_w%finalize() call MPI_Allreduce(this%user_field_bc_vel%bc_w%msk(0), integer_val, 1, & @@ -415,16 +388,13 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, & ! Setup our global field dirichlet bc call this%user_field_bc_vel%init_base(this%c_Xh) - call this%user_field_bc_vel%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_u', this%bc_labels) - call this%user_field_bc_vel%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_v', this%bc_labels) - call this%user_field_bc_vel%mark_zones_from_list(msh%labeled_zones,& - 'd_vel_w', this%bc_labels) + call this%user_field_bc_vel%mark_zones_from_list('d_vel_u', this%bc_labels) + call this%user_field_bc_vel%mark_zones_from_list('d_vel_v', this%bc_labels) + call this%user_field_bc_vel%mark_zones_from_list('d_vel_w', this%bc_labels) call this%user_field_bc_vel%finalize() ! Add the field bc to velocity bcs - call bc_list_add(this%bclst_vel, this%user_field_bc_vel) + call this%bclst_vel%append(this%user_field_bc_vel) ! ! Associate our field dirichlet update to the user one. @@ -446,11 +416,11 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, & call this%user_field_bc_vel%field_list%assign_to_field(4, & this%user_field_bc_prs%field_bc) - call bc_list_init(this%user_field_bc_vel%bc_list, size=4) + call this%user_field_bc_vel%bc_list%init(size=4) ! Note, bc_list_add only adds if the bc is not empty - call bc_list_add(this%user_field_bc_vel%bc_list, this%user_field_bc_vel%bc_u) - call bc_list_add(this%user_field_bc_vel%bc_list, this%user_field_bc_vel%bc_v) - call bc_list_add(this%user_field_bc_vel%bc_list, this%user_field_bc_vel%bc_w) + call this%user_field_bc_vel%bc_list%append(this%user_field_bc_vel%bc_u) + call this%user_field_bc_vel%bc_list%append(this%user_field_bc_vel%bc_v) + call this%user_field_bc_vel%bc_list%append(this%user_field_bc_vel%bc_w) ! ! Check if we need to output boundary types to a separate field @@ -541,46 +511,33 @@ subroutine fluid_scheme_init_all(this, msh, lx, params, kspv_init, kspp_init,& ! ! Setup pressure boundary conditions ! - call bc_list_init(this%bclst_prs) - call this%bc_prs%init_base(this%c_Xh) - call this%bc_prs%mark_zones_from_list(msh%labeled_zones,& - 'o', this%bc_labels) - call this%bc_prs%mark_zones_from_list(msh%labeled_zones,& - 'on', this%bc_labels) + call this%bclst_prs%init() + call this%bc_prs%init(this%c_Xh, params) + call this%bc_prs%mark_zones_from_list('o', this%bc_labels) + call this%bc_prs%mark_zones_from_list('on', this%bc_labels) ! Field dirichlet pressure bc call this%user_field_bc_prs%init_base(this%c_Xh) - call this%user_field_bc_prs%mark_zones_from_list(msh%labeled_zones,& - 'd_pres', this%bc_labels) + call this%user_field_bc_prs%mark_zones_from_list('d_pres', this%bc_labels) call this%user_field_bc_prs%finalize() call MPI_Allreduce(this%user_field_bc_prs%msk(0), integer_val, 1, & MPI_INTEGER, MPI_SUM, NEKO_COMM, ierr) if (integer_val .gt. 0) call this%user_field_bc_prs%init_field('d_pres') - call bc_list_add(this%bclst_prs, this%user_field_bc_prs) - call bc_list_add(this%user_field_bc_vel%bc_list, this%user_field_bc_prs) - - if (msh%outlet%size .gt. 0) then - call this%bc_prs%mark_zone(msh%outlet) - end if - if (msh%outlet_normal%size .gt. 0) then - call this%bc_prs%mark_zone(msh%outlet_normal) - end if + call this%bclst_prs%append(this%user_field_bc_prs) + call this%user_field_bc_vel%bc_list%append(this%user_field_bc_prs) call this%bc_prs%finalize() - call this%bc_prs%set_g(0.0_rp) - call bc_list_add(this%bclst_prs, this%bc_prs) + call this%bclst_prs%append(this%bc_prs) call this%bc_dong%init_base(this%c_Xh) - call this%bc_dong%mark_zones_from_list(msh%labeled_zones,& - 'o+dong', this%bc_labels) - call this%bc_dong%mark_zones_from_list(msh%labeled_zones,& - 'on+dong', this%bc_labels) + call this%bc_dong%mark_zones_from_list('o+dong', this%bc_labels) + call this%bc_dong%mark_zones_from_list('on+dong', this%bc_labels) call this%bc_dong%finalize() call this%bc_dong%init(this%c_Xh, params) - call bc_list_add(this%bclst_prs, this%bc_dong) + call this%bclst_prs%append(this%bc_dong) ! Pressure solver if (kspp_init) then @@ -668,7 +625,7 @@ subroutine fluid_scheme_free(this) call this%c_Xh%free() - call bc_list_free(this%bclst_vel) + call this%bclst_vel%free() call this%scratch%free() @@ -772,7 +729,7 @@ subroutine fluid_scheme_bc_apply_vel(this, t, tstep) real(kind=rp), intent(in) :: t integer, intent(in) :: tstep - call bc_list_apply_vector(this%bclst_vel,& + call this%bclst_vel%apply_vector(& this%u%x, this%v%x, this%w%x, this%dm_Xh%size(), t, tstep) end subroutine fluid_scheme_bc_apply_vel @@ -784,8 +741,7 @@ subroutine fluid_scheme_bc_apply_prs(this, t, tstep) real(kind=rp), intent(in) :: t integer, intent(in) :: tstep - call bc_list_apply_scalar(this%bclst_prs, this%p%x, & - this%p%dof%size(), t, tstep) + call this%bclst_prs%apply_scalar(this%p%x, this%p%dof%size(), t, tstep) end subroutine fluid_scheme_bc_apply_prs @@ -881,52 +837,41 @@ subroutine fluid_scheme_set_bc_type_output(this, params) this%bdry = 0.0_rp call bdry_mask%init_base(this%c_Xh) - call bdry_mask%mark_zone(this%msh%wall) - call bdry_mask%mark_zones_from_list(this%msh%labeled_zones,& - 'w', this%bc_labels) + call bdry_mask%mark_zones_from_list('w', this%bc_labels) call bdry_mask%finalize() call bdry_mask%set_g(1.0_rp) call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size()) call bdry_mask%free() call bdry_mask%init_base(this%c_Xh) - call bdry_mask%mark_zone(this%msh%inlet) - call bdry_mask%mark_zones_from_list(this%msh%labeled_zones,& - 'v', this%bc_labels) + call bdry_mask%mark_zones_from_list('v', this%bc_labels) call bdry_mask%finalize() call bdry_mask%set_g(2.0_rp) call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size()) call bdry_mask%free() call bdry_mask%init_base(this%c_Xh) - call bdry_mask%mark_zone(this%msh%outlet) - call bdry_mask%mark_zones_from_list(this%msh%labeled_zones,& - 'o', this%bc_labels) + call bdry_mask%mark_zones_from_list('o', this%bc_labels) call bdry_mask%finalize() call bdry_mask%set_g(3.0_rp) call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size()) call bdry_mask%free() call bdry_mask%init_base(this%c_Xh) - call bdry_mask%mark_zone(this%msh%sympln) - call bdry_mask%mark_zones_from_list(this%msh%labeled_zones,& - 'sym', this%bc_labels) + call bdry_mask%mark_zones_from_list('sym', this%bc_labels) call bdry_mask%finalize() call bdry_mask%set_g(4.0_rp) call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size()) call bdry_mask%free() call bdry_mask%init_base(this%c_Xh) - call bdry_mask%mark_zone(this%msh%outlet_normal) - call bdry_mask%mark_zones_from_list(this%msh%labeled_zones,& - 'on', this%bc_labels) + call bdry_mask%mark_zones_from_list('on', this%bc_labels) call bdry_mask%finalize() call bdry_mask%set_g(5.0_rp) call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size()) call bdry_mask%free() call bdry_mask%init_base(this%c_Xh) - call bdry_mask%mark_zone(this%msh%periodic) call bdry_mask%finalize() call bdry_mask%set_g(6.0_rp) call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size()) diff --git a/src/fluid/fluid_volflow.f90 b/src/fluid/fluid_volflow.f90 index 5a3577bf22e..bac7cc3a677 100644 --- a/src/fluid/fluid_volflow.f90 +++ b/src/fluid/fluid_volflow.f90 @@ -76,7 +76,7 @@ module fluid_volflow use json_module, only : json_file use json_utils, only: json_get use scratch_registry, only : scratch_registry_t - use bc, only : bc_list_t, bc_list_apply_vector, bc_list_apply_scalar + use bc_list, only : bc_list_t use ax_product, only : ax_t implicit none private @@ -217,7 +217,7 @@ subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, & end if call gs_Xh%op(p_res, GS_OP_ADD) - call bc_list_apply_scalar(bclst_dp, p_res%x, n) + call bclst_dp%apply_scalar(p_res%x, n) call pc_prs%update() ksp_result = ksp_prs%solve(Ax, p_vol, p_res%x, n, & c_Xh, bclst_dp, gs_Xh, prs_max_iter) @@ -238,7 +238,7 @@ subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, & call copy(ta2%x, c_Xh%B, n) call copy(ta3%x, c_Xh%B, n) end if - call bc_list_apply_vector(bclst_vel_res,& + call bclst_vel_res%apply_vector(& ta1%x, ta2%x, ta3%x, n) ! add forcing @@ -276,8 +276,7 @@ subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, & call gs_Xh%op(v_res, GS_OP_ADD) call gs_Xh%op(w_res, GS_OP_ADD) - call bc_list_apply_vector(bclst_vel_res,& - u_res%x, v_res%x, w_res%x, n) + call bclst_vel_res%apply_vector(u_res%x, v_res%x, w_res%x, n) call pc_vel%update() ksp_result = ksp_vel%solve(Ax, u_vol, u_res%x, n, & diff --git a/src/krylov/bcknd/cpu/bicgstab.f90 b/src/krylov/bcknd/cpu/bicgstab.f90 index 9d77ae760cd..d9a4ab28bc8 100644 --- a/src/krylov/bcknd/cpu/bicgstab.f90 +++ b/src/krylov/bcknd/cpu/bicgstab.f90 @@ -39,7 +39,7 @@ module bicgstab use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, NEKO_EPS, add2s2, x_update, & p_update, abscmp use utils, only : neko_error @@ -196,7 +196,7 @@ function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r call this%M%solve(p_hat, p, n) call Ax%compute(v, p_hat, coef, x%msh, x%Xh) call gs_h%op(v, n, GS_OP_ADD) - call bc_list_apply(blst, v, n) + call blst%apply(v, n) alpha = rho_1 / glsc3(f, coef%mult, v, n) call copy(s, r, n) call add2s2(s, v, -alpha, n) @@ -210,7 +210,7 @@ function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r call this%M%solve(s_hat, s, n) call Ax%compute(t, s_hat, coef, x%msh, x%Xh) call gs_h%op(t, n, GS_OP_ADD) - call bc_list_apply(blst, t, n) + call blst%apply(t, n) omega = glsc3(t, coef%mult, s, n) & / glsc3(t, coef%mult, t, n) call x_update(x%x, p_hat, s_hat, alpha, omega, n) diff --git a/src/krylov/bcknd/cpu/cacg.f90 b/src/krylov/bcknd/cpu/cacg.f90 index 9f2f1df6f08..86e2e216aa4 100644 --- a/src/krylov/bcknd/cpu/cacg.f90 +++ b/src/krylov/bcknd/cpu/cacg.f90 @@ -39,7 +39,7 @@ module cacg use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply, bc_list_apply_scalar + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, x_update, abscmp use utils, only : neko_warning use comm @@ -176,7 +176,7 @@ function cacg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resul if (mod(i,2) .eq. 0) then call Ax%compute(PR(1,i), PR(1,i-1), coef, x%msh, x%Xh) call gs_h%gs_op_vector(PR(1,i), n, GS_OP_ADD) - call bc_list_apply_scalar(blst, PR(1,i), n) + call blst%apply_scalar(PR(1,i), n) else call this%M%solve(PR(1,i), PR(1,i-1), n) end if @@ -188,7 +188,7 @@ function cacg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resul else call Ax%compute(PR(1,i+1), PR(1,i), coef, x%msh, x%Xh) call gs_h%gs_op_vector(PR(1,i+1), n, GS_OP_ADD) - call bc_list_apply_scalar(blst, PR(1,1+i), n) + call blst%apply_scalar(PR(1,1+i), n) end if end do diff --git a/src/krylov/bcknd/cpu/cg.f90 b/src/krylov/bcknd/cpu/cg.f90 index b2c3fa3200d..626245e2499 100644 --- a/src/krylov/bcknd/cpu/cg.f90 +++ b/src/krylov/bcknd/cpu/cg.f90 @@ -39,7 +39,7 @@ module cg use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, abscmp use comm implicit none @@ -178,7 +178,7 @@ function cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results call Ax%compute(w, p(1,p_cur), coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) pap = glsc3(w, coef%mult, p(1,p_cur), n) diff --git a/src/krylov/bcknd/cpu/gmres.f90 b/src/krylov/bcknd/cpu/gmres.f90 index 40eb600faee..a3f046ea6dc 100644 --- a/src/krylov/bcknd/cpu/gmres.f90 +++ b/src/krylov/bcknd/cpu/gmres.f90 @@ -39,7 +39,7 @@ module gmres use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, rone, copy, sub2, cmult2, abscmp use comm implicit none @@ -206,7 +206,7 @@ function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resu call copy(r, f, n) call Ax%compute(w, x%x, coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) call sub2(r, w, n) end if @@ -227,7 +227,7 @@ function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resu call Ax%compute(w, z(1,j), coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) do l = 1, j h(l,j) = 0.0_rp diff --git a/src/krylov/bcknd/cpu/pipecg.f90 b/src/krylov/bcknd/cpu/pipecg.f90 index 3092ec14e96..01d359d9ea7 100644 --- a/src/krylov/bcknd/cpu/pipecg.f90 +++ b/src/krylov/bcknd/cpu/pipecg.f90 @@ -39,7 +39,7 @@ module pipecg use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, abscmp use comm implicit none @@ -188,7 +188,7 @@ function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_res call this%M%solve(u(1,u_prev), r, n) call Ax%compute(w, u(1,u_prev), coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) rtr = glsc3(r, coef%mult, r, n) rnorm = sqrt(rtr)*norm_fac @@ -217,7 +217,7 @@ function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_res call this%M%solve(mi, w, n) call Ax%compute(ni, mi, coef, x%msh, x%Xh) call gs_h%op(ni, n, GS_OP_ADD) - call bc_list_apply(blst, ni, n) + call blst%apply(ni, n) call MPI_Wait(request, status, ierr) gamma2 = gamma1 diff --git a/src/krylov/bcknd/device/cg_device.f90 b/src/krylov/bcknd/device/cg_device.f90 index 8d977b5633c..923f4b69cf1 100644 --- a/src/krylov/bcknd/device/cg_device.f90 +++ b/src/krylov/bcknd/device/cg_device.f90 @@ -39,7 +39,7 @@ module cg_device use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : abscmp use device use device_math, only : device_rzero, device_copy, device_glsc3, & @@ -201,7 +201,7 @@ function cg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_ call Ax%compute(this%w, this%p, coef, x%msh, x%Xh) call gs_h%op(this%w, n, GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply(blst, this%w, n) + call blst%apply(this%w, n) pap = device_glsc3(this%w_d, coef%mult_d, this%p_d, n) diff --git a/src/krylov/bcknd/device/fusedcg_cpld_device.F90 b/src/krylov/bcknd/device/fusedcg_cpld_device.F90 index c4e1d499037..15ece0c104a 100644 --- a/src/krylov/bcknd/device/fusedcg_cpld_device.F90 +++ b/src/krylov/bcknd/device/fusedcg_cpld_device.F90 @@ -39,7 +39,7 @@ module fusedcg_cpld_device use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, abscmp use device_math, only : device_rzero, device_copy, device_glsc3, device_glsc2 use device @@ -67,7 +67,7 @@ module fusedcg_cpld_device real(kind=rp), allocatable :: alpha(:) type(c_ptr) :: w1_d = C_NULL_PTR type(c_ptr) :: w2_d = C_NULL_PTR - type(c_ptr) :: w3_d = C_NULL_PTR + type(c_ptr) :: w3_d = C_NULL_PTR type(c_ptr) :: r1_d = C_NULL_PTR type(c_ptr) :: r2_d = C_NULL_PTR type(c_ptr) :: r3_d = C_NULL_PTR @@ -261,10 +261,10 @@ subroutine fusedcg_cpld_device_init(this, n, max_iter, M, rel_tol, abs_tol) do i = 1, DEVICE_FUSEDCG_CPLD_P_SPACE+1 this%p1_d(i) = C_NULL_PTR call device_map(this%p1(:,i), this%p1_d(i), n) - + this%p2_d(i) = C_NULL_PTR call device_map(this%p2(:,i), this%p2_d(i), n) - + this%p3_d(i) = C_NULL_PTR call device_map(this%p3(:,i), this%p3_d(i), n) end do @@ -281,7 +281,7 @@ subroutine fusedcg_cpld_device_init(this, n, max_iter, M, rel_tol, abs_tol) HOST_TO_DEVICE, sync=.false.) ptr = c_loc(this%p3_d) call device_memcpy(ptr, this%p3_d_d, p_size, & - HOST_TO_DEVICE, sync=.false.) + HOST_TO_DEVICE, sync=.false.) if (present(rel_tol) .and. present(abs_tol)) then call this%ksp_init(max_iter, rel_tol, abs_tol) else if (present(rel_tol)) then @@ -308,11 +308,11 @@ subroutine fusedcg_cpld_device_free(this) if (allocated(this%w1)) then deallocate(this%w1) end if - + if (allocated(this%w2)) then deallocate(this%w2) end if - + if (allocated(this%w3)) then deallocate(this%w3) end if @@ -438,7 +438,7 @@ subroutine fusedcg_cpld_device_free(this) if (c_associated(this%gs_event2)) then call device_event_destroy(this%gs_event2) end if - + if (c_associated(this%gs_event3)) then call device_event_destroy(this%gs_event3) end if @@ -495,8 +495,8 @@ function fusedcg_cpld_device_solve_coupled(this, Ax, x, y, z, fx, fy, fz, & rtz1 = 1.0_rp p_prev = DEVICE_FUSEDCG_CPLD_P_SPACE p_cur = 1 - - + + call device_rzero(x%x_d, n) call device_rzero(y%x_d, n) call device_rzero(z%x_d, n) @@ -511,7 +511,7 @@ function fusedcg_cpld_device_solve_coupled(this, Ax, x, y, z, fx, fy, fz, & r2_d, r3_d, tmp_d, n) rtr = device_glsc3(tmp_d, coef%mult_d, coef%binv_d, n) - + rnorm = sqrt(rtr)*norm_fac ksp_results%res_start = rnorm ksp_results%res_final = rnorm @@ -531,7 +531,7 @@ function fusedcg_cpld_device_solve_coupled(this, Ax, x, y, z, fx, fy, fz, & beta = rtz1 / rtz2 if (iter .eq. 1) beta = 0.0_rp - + call device_fusedcg_cpld_update_p(p1_d(p_cur), p2_d(p_cur), p3_d(p_cur), & z1_d, z2_d, z3_d, p1_d(p_prev), p2_d(p_prev), p3_d(p_prev), beta, n) @@ -543,15 +543,15 @@ function fusedcg_cpld_device_solve_coupled(this, Ax, x, y, z, fx, fy, fz, & call device_event_sync(this%gs_event1) call device_event_sync(this%gs_event2) call device_event_sync(this%gs_event3) - call bc_list_apply(blstx, w1, n) - call bc_list_apply(blsty, w2, n) - call bc_list_apply(blstz, w3, n) + call blstx%apply(w1, n) + call blsty%apply(w2, n) + call blstz%apply(w3, n) call device_fusedcg_cpld_part1(w1_d, w2_d, w3_d, p1_d(p_cur), & p2_d(p_cur), p3_d(p_cur), tmp_d, n) pap = device_glsc2(tmp_d, coef%mult_d, n) - + alpha(p_cur) = rtz1 / pap rtr = device_fusedcg_cpld_part2(r1_d, r2_d, r3_d, coef%mult_d, & w1_d, w2_d, w3_d, alpha_d, alpha(p_cur), p_cur, n) @@ -595,9 +595,9 @@ function fusedcg_cpld_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) r ksp_results%res_final = 0.0 ksp_results%iter = 0 - + end function fusedcg_cpld_device_solve - + end module fusedcg_cpld_device diff --git a/src/krylov/bcknd/device/fusedcg_device.F90 b/src/krylov/bcknd/device/fusedcg_device.F90 index 2da425d409d..fe1222047e7 100644 --- a/src/krylov/bcknd/device/fusedcg_device.F90 +++ b/src/krylov/bcknd/device/fusedcg_device.F90 @@ -39,7 +39,7 @@ module fusedcg_device use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, abscmp use device_math, only : device_rzero, device_copy, device_glsc3 use device @@ -360,7 +360,7 @@ function fusedcg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result call Ax%compute(w, p(1, p_cur), coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) pap = device_glsc3(w_d, coef%mult_d, this%p_d(p_cur), n) diff --git a/src/krylov/bcknd/device/gmres_device.F90 b/src/krylov/bcknd/device/gmres_device.F90 index 6a5d585e05b..8ceaaf4b935 100644 --- a/src/krylov/bcknd/device/gmres_device.F90 +++ b/src/krylov/bcknd/device/gmres_device.F90 @@ -39,7 +39,7 @@ module gmres_device use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use device_identity, only : device_ident_t use math, only : rone, rzero, abscmp use device_math, only : device_rzero, device_copy, device_glsc3, & @@ -360,7 +360,7 @@ function gmres_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(k call Ax%compute(w, x%x, coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) call device_sub2(r_d, w_d, n) end if @@ -382,7 +382,7 @@ function gmres_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(k call Ax%compute(w, z(1,j), coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) if (NEKO_BCKND_OPENCL .eq. 1) then do i = 1, j diff --git a/src/krylov/bcknd/device/pipecg_device.F90 b/src/krylov/bcknd/device/pipecg_device.F90 index 942b538b6c0..e8186198777 100644 --- a/src/krylov/bcknd/device/pipecg_device.F90 +++ b/src/krylov/bcknd/device/pipecg_device.F90 @@ -39,7 +39,7 @@ module pipecg_device use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, copy, abscmp use device_math, only : device_rzero, device_copy, & device_glsc3, device_vlsc3 @@ -375,7 +375,7 @@ function pipecg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result( call Ax%compute(w, u(1,u_prev), coef, x%msh, x%Xh) call gs_h%op(w, n, GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply(blst, w, n) + call blst%apply(w, n) rtr = device_glsc3(r_d, coef%mult_d, r_d, n) rnorm = sqrt(rtr)*norm_fac @@ -403,7 +403,7 @@ function pipecg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result( call Ax%compute(ni, mi, coef, x%msh, x%Xh) call gs_h%op(ni, n, GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply(blst, ni, n) + call blst%apply(ni, n) call MPI_Wait(request, status, ierr) gamma2 = gamma1 diff --git a/src/krylov/bcknd/sx/cg_sx.f90 b/src/krylov/bcknd/sx/cg_sx.f90 index 0a1cf73b9ac..8c2c52bb3cc 100644 --- a/src/krylov/bcknd/sx/cg_sx.f90 +++ b/src/krylov/bcknd/sx/cg_sx.f90 @@ -39,7 +39,7 @@ module cg_sx use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, add2s1, abscmp implicit none private @@ -167,7 +167,7 @@ function sx_cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resu call Ax%compute(this%w, this%p, coef, x%msh, x%Xh) call gs_h%op(this%w, n, GS_OP_ADD) - call bc_list_apply(blst, this%w, n) + call blst%apply(this%w, n) pap = glsc3(this%w, coef%mult, this%p, n) diff --git a/src/krylov/bcknd/sx/gmres_sx.f90 b/src/krylov/bcknd/sx/gmres_sx.f90 index 16290c78c66..70835919428 100644 --- a/src/krylov/bcknd/sx/gmres_sx.f90 +++ b/src/krylov/bcknd/sx/gmres_sx.f90 @@ -39,7 +39,7 @@ module gmres_sx use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, rzero, rone, copy, cmult2, col2, col3, add2s2, abscmp use comm implicit none @@ -224,7 +224,7 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r call copy (this%r,f,n) call Ax%compute(this%w, x%x, coef, x%msh, x%Xh) call gs_h%op(this%w, n, GS_OP_ADD) - call bc_list_apply(blst, this%w, n) + call blst%apply(this%w, n) call add2s2(this%r,this%w,-one,n) call col2(this%r,this%ml,n) endif @@ -248,7 +248,7 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r call Ax%compute(this%w, this%z(1,j), coef, x%msh, x%Xh) call gs_h%op(this%w, n, GS_OP_ADD) - call bc_list_apply(blst, this%w, n) + call blst%apply(this%w, n) call col2(this%w, this%ml, n) do i = 1, j diff --git a/src/krylov/bcknd/sx/pipecg_sx.f90 b/src/krylov/bcknd/sx/pipecg_sx.f90 index 7bbe191e046..5dee2901f5e 100644 --- a/src/krylov/bcknd/sx/pipecg_sx.f90 +++ b/src/krylov/bcknd/sx/pipecg_sx.f90 @@ -39,7 +39,7 @@ module pipecg_sx use field, only : field_t use coefs, only : coef_t use gather_scatter, only : gs_t, GS_OP_ADD - use bc, only : bc_list_t, bc_list_apply + use bc_list, only : bc_list_t use math, only : glsc3, abscmp use comm implicit none @@ -178,7 +178,7 @@ function sx_pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_ call this%M%solve(this%u, this%r, n) call Ax%compute(this%w, this%u, coef, x%msh, x%Xh) call gs_h%op(this%w, n, GS_OP_ADD) - call bc_list_apply(blst, this%w, n) + call blst%apply(this%w, n) rtr = glsc3(this%r, coef%mult, this%r, n) rnorm = sqrt(rtr)*norm_fac @@ -209,7 +209,7 @@ function sx_pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_ call this%M%solve(this%mi, this%w, n) call Ax%compute(this%ni, this%mi, coef, x%msh, x%Xh) call gs_h%op(this%ni, n, GS_OP_ADD) - call bc_list_apply(blst, this%ni, n) + call blst%apply(this%ni, n) call MPI_Wait(request, status, ierr) gamma2 = gamma1 diff --git a/src/krylov/krylov.f90 b/src/krylov/krylov.f90 index b8c5e537835..375abcdf90a 100644 --- a/src/krylov/krylov.f90 +++ b/src/krylov/krylov.f90 @@ -40,7 +40,7 @@ module krylov use mesh, only : mesh_t use field, only : field_t use utils, only : neko_error, neko_warning - use bc, only : bc_list_t + use bc_list, only : bc_list_t use identity, only : ident_t use device_identity, only : device_ident_t use neko_config diff --git a/src/krylov/pc_hsmg.f90 b/src/krylov/pc_hsmg.f90 index 159fa67ead3..ee8fc4fdb09 100644 --- a/src/krylov/pc_hsmg.f90 +++ b/src/krylov/pc_hsmg.f90 @@ -84,6 +84,8 @@ module hsmg use mesh, only : mesh_t use krylov, only : ksp_t, ksp_monitor_t, KSP_MAX_ITER use krylov_fctry, only : krylov_solver_factory, krylov_solver_destroy + use zero_dirichlet, only : zero_dirichlet_t + use bc_list, only : bc_list_t !$ use omp_lib implicit none private @@ -107,7 +109,7 @@ module hsmg type(space_t) :: Xh_crs, Xh_mg !< spaces for lower levels type(dofmap_t) :: dm_crs, dm_mg type(coef_t) :: c_crs, c_mg - type(dirichlet_t) :: bc_crs, bc_mg, bc_reg + type(zero_dirichlet_t) :: bc_crs, bc_mg, bc_reg type(bc_list_t) :: bclst_crs, bclst_mg, bclst_reg type(schwarz_t) :: schwarz, schwarz_mg, schwarz_crs !< Schwarz decompostions type(field_t) :: e, e_mg, e_crs !< Solve fields @@ -210,28 +212,23 @@ subroutine hsmg_init(this, msh, Xh, coef, dof, gs_h, bclst, crs_pctype) call this%bc_crs%init_base(this%c_crs) call this%bc_mg%init_base(this%c_mg) call this%bc_reg%init_base(coef) - if (bclst%n .gt. 0) then - do i = 1, bclst%n - call this%bc_reg%mark_facets(bclst%bc(i)%bcp%marked_facet) - call this%bc_crs%mark_facets(bclst%bc(i)%bcp%marked_facet) - call this%bc_mg%mark_facets(bclst%bc(i)%bcp%marked_facet) + if (bclst%size() .gt. 0) then + do i = 1, bclst%size() + call this%bc_reg%mark_facets(bclst%items(i)%ptr%marked_facet) + call this%bc_crs%mark_facets(bclst%items(i)%ptr%marked_facet) + call this%bc_mg%mark_facets(bclst%items(i)%ptr%marked_facet) end do end if call this%bc_reg%finalize() - call this%bc_reg%set_g(real(0d0,rp)) - call bc_list_init(this%bclst_reg) - call bc_list_add(this%bclst_reg, this%bc_reg) - call this%bc_crs%finalize() - call this%bc_crs%set_g(real(0d0,rp)) - call bc_list_init(this%bclst_crs) - call bc_list_add(this%bclst_crs, this%bc_crs) - - call this%bc_mg%finalize() - call this%bc_mg%set_g(0.0_rp) - call bc_list_init(this%bclst_mg) - call bc_list_add(this%bclst_mg, this%bc_mg) + + call this%bclst_reg%init() + call this%bclst_crs%init() + call this%bclst_mg%init() + call this%bclst_reg%append(this%bc_reg) + call this%bclst_crs%append(this%bc_crs) + call this%bclst_mg%append(this%bc_mg) call this%schwarz%init(Xh, dof, gs_h, this%bclst_reg, msh) call this%schwarz_mg%init(this%Xh_mg, this%dm_mg, this%gs_mg,& @@ -368,7 +365,7 @@ subroutine hsmg_solve(this, z, r, n) r_d = device_get_ptr(r) !We should not work with the input call device_copy(this%r_d, r_d, n) - call bc_list_apply_scalar(this%bclst_reg, r, n) + call this%bclst_reg%apply_scalar(r, n) !OVERLAPPING Schwarz exchange and solve !! DOWNWARD Leg of V-cycle, we are pretty hardcoded here but w/e @@ -381,10 +378,9 @@ subroutine hsmg_solve(this, z, r, n) this%grids(2)%dof%size(), GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) call device_copy(this%r_d, r_d, n) - call bc_list_apply_scalar(this%bclst_reg, r, n) + call this%bclst_reg%apply_scalar(r, n) call device_copy(this%w_d, this%e%x_d, this%grids(2)%dof%size()) - call bc_list_apply_scalar(this%bclst_mg, this%w, & - this%grids(2)%dof%size()) + call this%bclst_mg%apply_scalar(this%w, this%grids(2)%dof%size()) !OVERLAPPING Schwarz exchange and solve call device_col2(this%w_d, this%grids(2)%coef%mult_d, & this%grids(2)%dof%size()) @@ -393,8 +389,7 @@ subroutine hsmg_solve(this, z, r, n) this%grids(1)%Xh) !Crs solve call device_copy(this%w_d, this%e%x_d, this%grids(2)%dof%size()) - call bc_list_apply_scalar(this%bclst_mg, this%w, & - this%grids(2)%dof%size()) + call this%bclst_mg%apply_scalar(this%w, this%grids(2)%dof%size()) !$omp parallel private(thrdid, nthrds) @@ -414,7 +409,7 @@ subroutine hsmg_solve(this, z, r, n) call this%grids(1)%gs_h%op(this%wf%x, & this%grids(1)%dof%size(), GS_OP_ADD, this%gs_event) call device_event_sync(this%gs_event) - call bc_list_apply_scalar(this%grids(1)%bclst, this%wf%x, & + call this%grids(1)%bclst%apply_scalar(this%wf%x, & this%grids(1)%dof%size()) call profiler_start_region('HSMG coarse-solve', 11) crs_info = this%crs_solver%solve(this%Ax, this%grids(1)%e, this%wf%x, & @@ -423,7 +418,7 @@ subroutine hsmg_solve(this, z, r, n) this%grids(1)%bclst, & this%grids(1)%gs_h, this%niter) call profiler_end_region - call bc_list_apply_scalar(this%grids(1)%bclst, this%grids(1)%e%x,& + call this%grids(1)%bclst%apply_scalar(this%grids(1)%e%x, & this%grids(1)%dof%size()) call profiler_end_region end if @@ -461,8 +456,7 @@ subroutine hsmg_solve(this, z, r, n) !Crs solve call this%grids(1)%gs_h%op(this%r, this%grids(1)%dof%size(), GS_OP_ADD) - call bc_list_apply_scalar(this%grids(1)%bclst, this%r, & - this%grids(1)%dof%size()) + call this%grids(1)%bclst%apply_scalar(this%r, this%grids(1)%dof%size()) call profiler_start_region('HSMG coarse-solve', 11) crs_info = this%crs_solver%solve(this%Ax, this%grids(1)%e, this%r, & this%grids(1)%dof%size(), & @@ -470,7 +464,7 @@ subroutine hsmg_solve(this, z, r, n) this%grids(1)%bclst, & this%grids(1)%gs_h, this%niter) call profiler_end_region - call bc_list_apply_scalar(this%grids(1)%bclst, this%grids(1)%e%x,& + call this%grids(1)%bclst%apply_scalar(this%grids(1)%e%x,& this%grids(1)%dof%size()) diff --git a/src/math/schwarz.f90 b/src/math/schwarz.f90 index a946d0ba04f..50b4656372a 100644 --- a/src/math/schwarz.f90 +++ b/src/math/schwarz.f90 @@ -71,6 +71,7 @@ module schwarz use fdm, only : fdm_t use device use neko_config + use bc_list, only : bc_list_t use, intrinsic :: iso_c_binding implicit none private @@ -425,11 +426,11 @@ subroutine schwarz_compute(this, e, r) call device_event_sync(this%event) call this%gs_h%op(e, n, GS_OP_ADD, this%event) - call bc_list_apply_scalar(this%bclst, e, n) + call this%bclst%apply_scalar(e, n) call device_col2(e_d,this%wt_d, n) call device_stream_wait_event(aux_cmd_queue, this%event, 0) else - call bc_list_apply_scalar(this%bclst, r, n) + call this%bclst%apply_scalar(r, n) call schwarz_toext3d(work1, r, this%Xh%lx, this%msh%nelv) ! exchange interior nodes @@ -454,7 +455,7 @@ subroutine schwarz_compute(this, e, r) ! sum border nodes call this%gs_h%op(e, n, GS_OP_ADD) - call bc_list_apply_scalar(this%bclst, e, n) + call this%bclst%apply_scalar(e, n) call schwarz_wt3d(e, this%wt, this%Xh%lx, this%msh%nelv) end if diff --git a/src/neko.f90 b/src/neko.f90 index 134edf1933b..c31406f52ff 100644 --- a/src/neko.f90 +++ b/src/neko.f90 @@ -56,7 +56,7 @@ module neko use gather_scatter use coefs use bc - use wall + use zero_dirichlet use dirichlet use krylov_fctry use precon_fctry diff --git a/src/scalar/scalar_pnpn.f90 b/src/scalar/scalar_pnpn.f90 index f4aea2ea8e7..a38fcba3af9 100644 --- a/src/scalar/scalar_pnpn.f90 +++ b/src/scalar/scalar_pnpn.f90 @@ -42,8 +42,7 @@ module scalar_pnpn use dirichlet, only : dirichlet_t use neumann, only : neumann_t use field, only : field_t - use bc, only : bc_list_t, bc_list_init, bc_list_free, bc_list_apply_scalar, & - bc_list_add + use bc_list, only : bc_list_t use mesh, only : mesh_t use checkpoint, only : chkp_t use coefs, only : coef_t @@ -68,6 +67,7 @@ module scalar_pnpn use user_intf, only : user_t use material_properties, only : material_properties_t use neko_config, only : NEKO_BCKND_DEVICE + use zero_dirichlet, only : zero_dirichlet_t use time_step_controller, only : time_step_controller_t use scratch_registry, only: neko_scratch_registry implicit none @@ -91,10 +91,12 @@ module scalar_pnpn !> Dirichlet conditions for the residual !! Collects all the Dirichlet condition facets into one bc and applies 0, !! Since the values never change there during the solve. - type(dirichlet_t) :: bc_res + type(zero_dirichlet_t) :: bc_res !> A bc list for the bc_res. Contains only that, essentially just to wrap !! the if statement determining whether to apply on the device or CPU. + !! Also needed since a bc_list is the type that is sent to, e.g. solvers, + !! cannot just send `bc_res` on its own. type(bc_list_t) :: bclst_ds !> Advection operator. @@ -177,23 +179,16 @@ subroutine scalar_pnpn_init(this, msh, coef, gs, params, user, & ! Initialize dirichlet bcs for scalar residual ! todo: look that this works - call this%bc_res%init_base(this%c_Xh) - do i = 1, this%n_dir_bcs - call this%bc_res%mark_facets(this%dir_bcs(i)%marked_facet) + call this%bc_res%init(this%c_Xh, params) + do i = 1, this%n_strong + call this%bc_res%mark_facets(this%bcs%items(i)%ptr%marked_facet) end do - ! Check for user bcs - if (this%user_bc%msk(0) .gt. 0) then - call this%bc_res%mark_facets(this%user_bc%marked_facet) - end if - - call this%bc_res%mark_zones_from_list(msh%labeled_zones, 'd_s', & - this%bc_labels) +! call this%bc_res%mark_zones_from_list('d_s', this%bc_labels) call this%bc_res%finalize() - call this%bc_res%set_g(0.0_rp) - call bc_list_init(this%bclst_ds) - call bc_list_add(this%bclst_ds, this%bc_res) + call this%bclst_ds%init() + call this%bclst_ds%append(this%bc_res) ! Intialize projection space @@ -244,7 +239,7 @@ subroutine scalar_pnpn_free(this) !Deallocate scalar field call this%scheme_free() - call bc_list_free(this%bclst_ds) + call this%bclst_ds%free() call this%proj_s%free() call this%s_res%free() @@ -302,17 +297,8 @@ subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) if_variable_dt => dt_controller%if_variable_dt, & dt_last_change => dt_controller%dt_last_change) - if (neko_log%level_ .ge. NEKO_LOG_DEBUG) then - write(log_buf,'(A,A,E15.7,A,E15.7,A,E15.7)') 'Scalar debug',& - ' l2norm s', glsc2(this%s%x,this%s%x,n),& - ' slag1', glsc2(this%slag%lf(1)%x,this%slag%lf(1)%x,n),& - ' slag2', glsc2(this%slag%lf(2)%x,this%slag%lf(2)%x,n) - call neko_log%message(log_buf) - write(log_buf,'(A,A,E15.7,A,E15.7)') 'Scalar debug2',& - ' l2norm abx1', glsc2(this%abx1%x,this%abx1%x,n),& - ' abx2', glsc2(this%abx2%x,this%abx2%x,n) - call neko_log%message(log_buf) - end if + ! Logs extra information the log level is NEKO_LOG_DEBUG or above. + call print_debug(this) ! Compute the source terms call this%source_term%compute(t, tstep) @@ -323,9 +309,9 @@ subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) else call col2(f_Xh%x, c_Xh%B, n) end if - - ! Apply Neumann boundary conditions - call bc_list_apply_scalar(this%bclst_neumann, this%f_Xh%x, dm_Xh%size()) + ! Apply weak boundary conditions, that contribute to the source terms. + call this%bcs%apply_scalar(this%f_Xh%x, dm_Xh%size(), t, & + tstep, strong=.false.) ! Add the advection operators to the right-hans-side. call this%adv%compute_scalar(u, v, w, s, f_Xh, & @@ -344,12 +330,13 @@ subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) call slag%update() - !> Apply Dirichlet boundary conditions + !> Apply strong boundary conditions. !! We assume that no change of boundary conditions !! occurs between elements. i.e. we do not apply gsop here like in Nek5000 call this%field_dir_bc%update(this%field_dir_bc%field_list, & this%field_dirichlet_bcs, this%c_Xh, t, tstep, "scalar") - call bc_list_apply_scalar(this%bclst_dirichlet, this%s%x, this%dm_Xh%size()) + call this%bcs%apply_scalar(this%s%x, this%dm_Xh%size(), t, tstep, & + strong=.true.) ! Update material properties if necessary @@ -362,8 +349,9 @@ subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) call gs_Xh%op(s_res, GS_OP_ADD) + ! Apply a 0-valued Dirichlet boundary conditions on the ds. - call bc_list_apply_scalar(this%bclst_ds, s_res%x, dm_Xh%size()) + call this%bclst_ds%apply_scalar(s_res%x, dm_Xh%size()) call profiler_end_region @@ -375,8 +363,8 @@ subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) c_Xh, this%bclst_ds, gs_Xh) call profiler_end_region - call this%proj_s%post_solving(ds%x, Ax, c_Xh, & - this%bclst_ds, gs_Xh, n, tstep, dt_controller) + call this%proj_s%post_solving(ds%x, Ax, c_Xh, this%bclst_ds, gs_Xh, & + n, tstep, dt_controller) ! Update the solution if (NEKO_BCKND_DEVICE .eq. 1) then @@ -391,5 +379,25 @@ subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller) call profiler_end_region end subroutine scalar_pnpn_step + subroutine print_debug(this) + class(scalar_pnpn_t), intent(inout) :: this + character(len=LOG_SIZE) :: log_buf + integer :: n + + n = this%dm_Xh%size() + + if (neko_log%level_ .ge. NEKO_LOG_DEBUG) then + write(log_buf,'(A,A,E15.7,A,E15.7,A,E15.7)') 'Scalar debug',& + ' l2norm s', glsc2(this%s%x,this%s%x,n),& + ' slag1', glsc2(this%slag%lf(1)%x,this%slag%lf(1)%x,n),& + ' slag2', glsc2(this%slag%lf(2)%x,this%slag%lf(2)%x,n) + call neko_log%message(log_buf) + write(log_buf,'(A,A,E15.7,A,E15.7)') 'Scalar debug2',& + ' l2norm abx1', glsc2(this%abx1%x,this%abx1%x,n),& + ' abx2', glsc2(this%abx2%x,this%abx2%x,n) + call neko_log%message(log_buf) + end if + end subroutine print_debug + end module scalar_pnpn diff --git a/src/scalar/scalar_scheme.f90 b/src/scalar/scalar_scheme.f90 index b3ba30e6c9f..7e8d93154f4 100644 --- a/src/scalar/scalar_scheme.f90 +++ b/src/scalar/scalar_scheme.f90 @@ -50,10 +50,10 @@ module scalar_scheme use device_jacobi, only : device_jacobi_t use sx_jacobi, only : sx_jacobi_t use hsmg, only : hsmg_t + use bc, only : bc_t + use bc_list, only : bc_list_t use precon_fctry, only : precon_factory, precon_destroy use precon, only : pc_t - use bc, only : bc_t, bc_list_t, bc_list_free, bc_list_init, & - bc_list_apply_scalar, bc_list_add use field_dirichlet, only: field_dirichlet_t, field_dirichlet_update use mesh, only : mesh_t, NEKO_MSH_MAX_ZLBLS, NEKO_MSH_MAX_ZLBL_LEN use facet_zone, only : facet_zone_t @@ -61,17 +61,18 @@ module scalar_scheme use logger, only : neko_log, LOG_SIZE use field_registry, only : neko_field_registry use usr_scalar, only : usr_scalar_t, usr_scalar_bc_eval - use json_utils, only : json_get, json_get_or_default - use json_module, only : json_file + use json_utils, only : json_get, json_get_or_default, json_extract_item + use json_module, only : json_file, json_core, json_value use user_intf, only : user_t use material_properties, only : material_properties_t use utils, only : neko_error use comm, only: NEKO_COMM, MPI_INTEGER, MPI_SUM use scalar_source_term, only : scalar_source_term_t + use field_series, only : field_series_t + use bc_fctry, only : bc_factory use math, only : cfill, add2s2 use device_math, only : device_cfill, device_add2s2 use neko_config, only : NEKO_BCKND_DEVICE - use field_series use time_step_controller, only : time_step_controller_t implicit none @@ -109,24 +110,16 @@ module scalar_scheme integer :: projection_activ_step !> Preconditioner. class(pc_t), allocatable :: pc - !> Dirichlet conditions. - type(dirichlet_t) :: dir_bcs(NEKO_MSH_MAX_ZLBLS) !> Field Dirichlet conditions. type(field_dirichlet_t) :: field_dir_bc !> List of BC objects to pass to user_dirichlet_update type(bc_list_t) :: field_dirichlet_bcs - !> Neumann conditions. - type(neumann_t) :: neumann_bcs(NEKO_MSH_MAX_ZLBLS) + !> Number of strong bcs. + integer :: n_strong = 0 + !> List of boundary conditions, including the user one. + type(bc_list_t) :: bcs !> User Dirichlet conditions. type(usr_scalar_t) :: user_bc - !> Number of Dirichlet bcs. - integer :: n_dir_bcs = 0 - !> Number of Neumann bcs. - integer :: n_neumann_bcs = 0 - !> List of Dirichlet boundary conditions, including the user one. - type(bc_list_t) :: bclst_dirichlet - !> List of Neumann conditions list - type(bc_list_t) :: bclst_neumann !> Case paramters. type(json_file), pointer :: params !> Mesh. @@ -157,7 +150,7 @@ module scalar_scheme !> Validate successful initialization. procedure, pass(this) :: validate => scalar_scheme_validate !> Assigns the evaluation function for `user_bc`. - procedure, pass(this) :: set_user_bc => scalar_scheme_set_user_bc + !procedure, pass(this) :: set_user_bc => scalar_scheme_set_user_bc !> Update variable material properties procedure, pass(this) :: update_material_properties => & scalar_scheme_update_material_properties @@ -230,70 +223,44 @@ end subroutine scalar_scheme_step_intrf contains !> Initialize boundary conditions - !! @param zones List of zones - !! @param bc_labels List of user specified bcs from the parameter file - !! currently dirichlet 'd=X' and 'user' supported - subroutine scalar_scheme_add_bcs(this, zones, bc_labels) + !! @param user The user object binding the user-defined routines. + subroutine scalar_scheme_setup_bcs(this, user) class(scalar_scheme_t), intent(inout) :: this - type(facet_zone_t), intent(inout) :: zones(NEKO_MSH_MAX_ZLBLS) - character(len=NEKO_MSH_MAX_ZLBL_LEN), intent(in) :: bc_labels(:) - character(len=NEKO_MSH_MAX_ZLBL_LEN) :: bc_label - integer :: i - real(kind=rp) :: dir_value, flux + type(user_t), target, intent(in) :: user + integer :: i, j, n_bcs, ierr + real(kind=rp) :: dir_value, flux_value logical :: bc_exists + type(json_core) :: core + type(json_value), pointer :: bc_object + type(json_file) :: bc_subdict + logical :: found - do i = 1, size(bc_labels) - bc_label = trim(bc_labels(i)) - if (bc_label(1:2) .eq. 'd=') then -! The idea of this commented piece of code is to merge bcs with the same -! Dirichlet value into 1 so that one has less kernel launches. Currently -! segfaults, needs investigation. -! bc_exists = .false. -! bc_idx = 0 -! do j = 1, i-1 -! if (bc_label .eq. bc_labels(j)) then -! bc_exists = .true. -! bc_idx = j -! end if -! end do - -! if (bc_exists) then -! call this%dir_bcs(j)%mark_zone(zones(i)) -! else - this%n_dir_bcs = this%n_dir_bcs + 1 - call this%dir_bcs(this%n_dir_bcs)%init_base(this%c_Xh) - call this%dir_bcs(this%n_dir_bcs)%mark_zone(zones(i)) - read(bc_label(3:), *) dir_value - call this%dir_bcs(this%n_dir_bcs)%set_g(dir_value) - call this%dir_bcs(this%n_dir_bcs)%finalize() - end if + call this%bcs%init() - if (bc_label(1:2) .eq. 'n=') then - this%n_neumann_bcs = this%n_neumann_bcs + 1 - call this%neumann_bcs(this%n_neumann_bcs)%init_base(this%c_Xh) - call this%neumann_bcs(this%n_neumann_bcs)%mark_zone(zones(i)) - read(bc_label(3:), *) flux - call this%neumann_bcs(this%n_neumann_bcs)%finalize_neumann(flux) - end if + if (this%params%valid_path('case.scalar.boundary_conditions')) then + call this%params%info('case.scalar.boundary_conditions', n_children=n_bcs) + call this%params%get_core(core) + call this%params%get('case.scalar.boundary_conditions', bc_object, found) + call this%bcs%init(n_bcs) - !> Check if user bc on this zone - if (bc_label(1:4) .eq. 'user') then - call this%user_bc%mark_zone(zones(i)) - end if + ! Gotta set this manually, becase we don't append to the list but + ! rather directly allocate in the factory. + this%bcs%size_ = n_bcs - end do - do i = 1, this%n_dir_bcs - call bc_list_add(this%bclst_dirichlet, this%dir_bcs(i)) - end do + do i=1, n_bcs + ! Create a new json containing just the subdict for this bc + call json_extract_item(core, bc_object, i, bc_subdict) - ! Create list with just Neumann bcs - call bc_list_init(this%bclst_neumann, this%n_neumann_bcs) - do i=1, this%n_neumann_bcs - call bc_list_add(this%bclst_neumann, this%neumann_bcs(i)) - end do + call bc_factory(this%bcs%items(i)%ptr, bc_subdict, & + this%c_Xh, user) - end subroutine scalar_scheme_add_bcs + if (this%bcs%strong(i)) then + this%n_strong = this%n_strong + 1 + end if + end do + end if + end subroutine scalar_scheme_setup_bcs !> Initialize all related components of the current scheme !! @param msh The mesh. @@ -397,19 +364,8 @@ subroutine scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, scheme, user, & ! ! Setup scalar boundary conditions ! - call bc_list_init(this%bclst_dirichlet) - call this%user_bc%init_base(this%c_Xh) - - ! Read boundary types from the case file - allocate(this%bc_labels(NEKO_MSH_MAX_ZLBLS)) - - ! A filler value - this%bc_labels = "not" - - if (params%valid_path('case.scalar.boundary_types')) then - call json_get(params, 'case.scalar.boundary_types', this%bc_labels,& - 'not') - end if + !call bc_list_init(this%bclst_dirichlet) + !call this%user_bc%init_base(this%c_Xh) ! ! Setup right-hand side field. @@ -420,37 +376,34 @@ subroutine scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, scheme, user, & ! Initialize the source term call this%source_term%init(params, this%f_Xh, this%c_Xh, user) - call scalar_scheme_add_bcs(this, msh%labeled_zones, this%bc_labels) + ! Set up boundary conditions + call scalar_scheme_setup_bcs(this, user) + +!! COMMENTING USER STUFF ! Mark BC zones - call this%user_bc%mark_zone(msh%wall) - call this%user_bc%mark_zone(msh%inlet) - call this%user_bc%mark_zone(msh%outlet) - call this%user_bc%mark_zone(msh%outlet_normal) - call this%user_bc%mark_zone(msh%sympln) - call this%user_bc%finalize() - if (this%user_bc%msk(0) .gt. 0) call bc_list_add(this%bclst_dirichlet,& - this%user_bc) +! call this%user_bc%finalize() +! if (this%user_bc%msk(0) .gt. 0) call bc_list_add(this%bclst_dirichlet,& +! this%user_bc) ! Add field dirichlet BCs - call this%field_dir_bc%init_base(this%c_Xh) - call this%field_dir_bc%mark_zones_from_list(msh%labeled_zones, & - 'd_s', this%bc_labels) - call this%field_dir_bc%finalize() - call MPI_Allreduce(this%field_dir_bc%msk(0), integer_val, 1, & - MPI_INTEGER, MPI_SUM, NEKO_COMM, ierr) - if (integer_val .gt. 0) call this%field_dir_bc%init_field('d_s') +! call this%field_dir_bc%init_base(this%c_Xh) +! call this%field_dir_bc%mark_zones_from_list(msh%labeled_zones, & +! 'd_s', this%bc_labels) +! call this%field_dir_bc%finalize() +! call MPI_Allreduce(this%field_dir_bc%msk(0), integer_val, 1, & +! MPI_INTEGER, MPI_SUM, NEKO_COMM, ierr) +! if (integer_val .gt. 0) call this%field_dir_bc%init_field('d_s') - call bc_list_add(this%bclst_dirichlet, this%field_dir_bc) +! call bc_list_add(this%bclst_dirichlet, this%field_dir_bc) ! ! Associate our field dirichlet update to the user one. ! - this%field_dir_bc%update => user%user_dirichlet_update - - call bc_list_init(this%field_dirichlet_bcs, size=1) - call bc_list_add(this%field_dirichlet_bcs, this%field_dir_bc) +! this%field_dir_bc%update => user%user_dirichlet_update +! call bc_list_init(this%field_dirichlet_bcs, size=1) +! call bc_list_add(this%field_dirichlet_bcs, this%field_dir_bc) ! todo parameter file ksp tol should be added call json_get_or_default(params, 'case.fluid.velocity_solver.max_iterations',& @@ -458,7 +411,7 @@ subroutine scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, scheme, user, & call scalar_scheme_solver_factory(this%ksp, this%dm_Xh%size(), & solver_type, integer_val, solver_abstol) call scalar_scheme_precon_factory(this%pc, this%ksp, & - this%c_Xh, this%dm_Xh, this%gs_Xh, this%bclst_dirichlet, solver_precon) + this%c_Xh, this%dm_Xh, this%gs_Xh, this%bcs, solver_precon) call neko_log%end_section() @@ -485,20 +438,15 @@ subroutine scalar_scheme_free(this) deallocate(this%pc) end if - if (allocated(this%bc_labels)) then - deallocate(this%bc_labels) - end if - call this%source_term%free() - call bc_list_free(this%bclst_dirichlet) - call bc_list_free(this%bclst_neumann) + call this%bcs%free() call this%lambda_field%free() call this%slag%free() ! Free everything related to field dirichlet BCs - call bc_list_free(this%field_dirichlet_bcs) + call this%field_dirichlet_bcs%free() call this%field_dir_bc%free() end subroutine scalar_scheme_free @@ -631,5 +579,4 @@ subroutine scalar_scheme_update_material_properties(this) end subroutine scalar_scheme_update_material_properties - end module scalar_scheme