Merge pull request #544 from ExtremeFLOW/release/0.4

NEKO v0.4.1

Makefile.am

@@ -44,6 +44,12 @@ EXTRA_DIST = \
 	examples/poisson/ax_poisson.f90\
 	examples/poisson/driver.f90\
 	examples/poisson/setup.f90\
+	examples/cyl_boundary_layer/cyl_bl_rot.case\
+	examples/cyl_boundary_layer/cyl_bl.f90\
+	examples/cyl_boundary_layer/README.md\
+	examples/cyl_boundary_layer/cyl_bl_user.case\
+	examples/cyl_boundary_layer/cyl_bl_basic.case\
+	examples/cyl_boundary_layer/cyl.nmsh\
 	reframe/checks.py\
 	reframe/settings.py\
 	reframe/README.md\

configure.ac

@@ -1,4 +1,4 @@
-AC_INIT([neko],[0.4.0])
+AC_INIT([neko],[0.4.1])
 AM_INIT_AUTOMAKE([foreign subdir-objects])
 AM_MAINTAINER_MODE
 AC_CONFIG_MACRO_DIR([m4])

examples/cyl_boundary_layer/README.md

@@ -0,0 +1,8 @@
+## Simulation of a wall-mounted cylinder
+In this case we simulate a cylinder mounted on a no-slip wall in a open channel. We provide 3 different case files that use the same mesh.
+
+* cyl_bl_basic.case can be run with the usual'neko' executable
+* cyl_bl_user.case uses the user file cyl_bl.f90 and the neko executable that is generated with makeneko.
+* cyl_bl_rot.case also needs the the user file and changes the boundary conditions to simulate a rotating cylinder with dong outflow conditions.
+
+The Reynolds number can probably be increased in these simulations to get some more interesting results.

examples/cyl_boundary_layer/cyl_bl.f90

@@ -0,0 +1,150 @@
+!>User file to set up at rotating cylinder. Martin Karp 30/5-2022
+
+module user
+  use neko
+  implicit none
+  ! Case parameters
+  real(kind=rp), parameter :: h = 1.0
+  real(kind=rp), parameter :: gam = 20.0
+  real(kind=rp), parameter :: rad = h/gam
+  real(kind=rp), parameter :: n = 7
+  real(kind=rp), parameter :: pw = 1/n
+  real(kind=rp), parameter :: ucl = 1
+  real(kind=rp), parameter :: alpha = 1
+  real(kind=rp), parameter :: u_th2 = ucl*alpha
+  real(kind=rp), parameter :: u_rho = 0.0
+  real(kind=rp), parameter :: u_axial = 0.0
+  real(kind=rp), parameter :: y0 = 0.0
+  real(kind=rp), parameter :: y1 = 0.0
+  real(kind=rp), parameter :: delta = 0.005*h
+contains
+
+
+  ! Register user defined functions (see user_intf.f90)
+  subroutine user_setup(u)
+    type(user_t), intent(inout) :: u
+    u%fluid_usr_ic => user_ic
+    u%fluid_usr_if => user_inflow_eval
+    !u%usr_chk => user_do_stuff
+    u%usr_msh_setup => cylinder_deform
+  end subroutine user_setup
+
+  subroutine cylinder_deform(msh)
+    type(mesh_t), intent(inout) :: msh
+    msh%apply_deform => cylinder_gen_curve
+  end subroutine cylinder_deform
+
+  !> Make the cylinder more circular!
+  subroutine cylinder_gen_curve(msh, x, y, z, lx, ly, lz)
+    class(mesh_t) :: msh
+    integer, intent(in) :: lx, ly, lz
+    real(kind=rp), intent(inout) :: x(lx, lx, lx, msh%nelv)
+    real(kind=rp), intent(inout) :: y(lx, lx, lx, msh%nelv)
+    real(kind=rp), intent(inout) :: z(lx, lx, lx, msh%nelv)
+    type(tuple_i4_t) :: el_and_facet
+    real(kind=rp) :: th
+    integer :: e, i, j ,k, l,  facet
+
+    !The cylinders zone number is 7
+    do l = 1,msh%labeled_zones(7)%size
+       el_and_facet = msh%labeled_zones(7)%facet_el(l)
+       facet = el_and_facet%x(1)
+       e = el_and_facet%x(2)
+       do k = 1, lz
+          do j = 1, ly
+              do i = 1, lx
+                 if (index_is_on_facet(i,j,k,lx,ly,lz, facet)) then
+                    th = atan2(z(i,j,k,e), x(i,j,k,e))
+                    x(i,j,k,e) = rad * cos(th)
+                    z(i,j,k,e) = rad * sin(th) 
+                 end if
+              end do
+          end do
+       end do
+    end do
+  end subroutine cylinder_gen_curve
+
+  subroutine user_inflow_eval(u, v, w, x, y, z, nx, ny, nz, ix, iy, iz, ie)
+    real(kind=rp), intent(inout) :: u
+    real(kind=rp), intent(inout) :: v
+    real(kind=rp), intent(inout) :: w
+    real(kind=rp), intent(in) :: x
+    real(kind=rp), intent(in) :: y
+    real(kind=rp), intent(in) :: z
+    real(kind=rp), intent(in) :: nx
+    real(kind=rp), intent(in) :: ny
+    real(kind=rp), intent(in) :: nz
+    integer, intent(in) :: ix
+    integer, intent(in) :: iy
+    integer, intent(in) :: iz
+    integer, intent(in) :: ie
+    real(kind=rp) ::  u_th,dist,th, yy
+    real(kind=rp) ::  arg
+
+!   Two different regions (inflow & cyl) have the label 'v  '
+!   Let compute the distance from the (0,0) in the x-y plane
+!   to identify the proper one
+    dist = sqrt(x**2 + z**2)
+
+! --- INFLOW
+    if (dist .gt. 1.1*rad) then
+       u =  ucl*y**pw
+    end if
+! --- 
+
+    w = 0.0
+    v = 0.0
+! --- SPINNING CYLINDER
+
+    if (dist.lt.1.5*rad .and. y.gt. 0.1) then                      
+       th = atan2(z,x)
+       u = cos(th)*u_rho - sin(th)*u_th2
+       w = sin(th)*u_rho + cos(th)*u_th2   
+    end if     
+
+! --- 
+
+
+!     Smoothing function for the velocity u_th on the spinning cylinder
+!     to avoid gap in the at the bottom wall
+
+!     u_th is smoothed if z0 < z < delta
+!     u_th=1 if z >= delta
+
+
+    yy = y + abs(y0) ! coordinate shift 
+
+    if (dist .lt. 1.5*rad) then 
+       if (yy.lt.delta) then
+          arg  = yy/delta
+          u_th = u_th2/(1.0_rp+exp(1.0_rp/(arg-1.0_rp)+1.0_rp/arg))
+       else
+          u_th = u_th2
+       endif
+
+       th = atan2(z,x)
+
+       u = cos(th)*u_rho - sin(th)*u_th
+       w = sin(th)*u_rho + cos(th)*u_th  
+    end if
+  end subroutine user_inflow_eval
+
+  ! User defined initial condition
+  subroutine user_ic(u, v, w, p, params)
+    type(field_t), intent(inout) :: u
+    type(field_t), intent(inout) :: v
+    type(field_t), intent(inout) :: w
+    type(field_t), intent(inout) :: p
+    type(param_t), intent(inout) :: params
+    integer :: i
+    real(kind=rp) :: y
+
+    do i = 1, u%dof%size()
+       y = u%dof%y(i,1,1,1)
+       u%x(i,1,1,1) =  ucl*y**pw
+       v%x(i,1,1,1) = 0.0
+       w%x(i,1,1,1) = 0.0
+    end do
+  end subroutine user_ic
+
+end module user

examples/cyl_boundary_layer/cyl_bl_basic.case

@@ -0,0 +1,25 @@
+&NEKO_CASE
+mesh_file= 'cyl.nmsh'
+fluid_scheme='pnpn'
+lx = 6
+source_term='noforce'
+initial_condition = 'uniform'
+/
+&NEKO_PARAMETERS
+dt = 8d-4
+T_end = 100
+nsamples = 400
+dealias=.true.
+uinf= 1.0,0.0,0.0
+output_bdry = .true.
+rho = 1
+Re = 500
+abstol_prs = 1d-5
+abstol_vel = 1d-8
+pc_vel = 'jacobi'
+pc_prs = 'hsmg'
+proj_prs_dim = 20
+proj_vel_dim = 3
+bc_labels='v','o','sym','w','on','on','w'
+fluid_inflow="default"
+/

examples/cyl_boundary_layer/cyl_bl_rot.case

@@ -0,0 +1,27 @@
+&NEKO_CASE
+mesh_file= 'cyl.nmsh'
+fluid_scheme='pnpn'
+lx = 6
+source_term='noforce'
+initial_condition = 'user'
+/
+&NEKO_PARAMETERS
+dt = 5d-4
+T_end = 100
+nsamples = 400
+dealias=.true.
+uinf= 1.0,0.0,0.0
+output_bdry = .true.
+rho = 1
+Re = 500
+abstol_prs = 1d-5
+abstol_vel = 1d-8
+pc_vel = 'jacobi'
+pc_prs = 'hsmg'
+proj_prs_dim = 20
+proj_vel_dim = 3
+bc_labels='v','o+dong','sym','w','on','on+dong','v'
+fluid_inflow="user"
+dong_uchar = 1.0
+dong_delta = 0.01
+/

examples/cyl_boundary_layer/cyl_bl_user.case

@@ -0,0 +1,25 @@
+&NEKO_CASE
+mesh_file= 'cyl.nmsh'
+fluid_scheme='pnpn'
+lx = 6
+source_term='noforce'
+initial_condition = 'user'
+/
+&NEKO_PARAMETERS
+dt = 8d-4
+T_end = 100
+nsamples = 400
+dealias=.true.
+uinf= 1.0,0.0,0.0
+output_bdry = .true.
+rho = 1
+Re = 500
+abstol_prs = 1d-5
+abstol_vel = 1d-8
+pc_vel = 'jacobi'
+pc_prs = 'hsmg'
+proj_prs_dim = 10
+proj_vel_dim = 3
+bc_labels='v','o','sym','w','on','on','w'
+fluid_inflow="user"
+/

src/io/file.f90

@@ -132,7 +132,12 @@ end subroutine file_read
   subroutine file_set_counter(this, n)
     class(file_t), intent(inout) :: this
     integer, intent(in) :: n
-    call this%file_type%set_counter(n)
+
+    select type(ft => this%file_type)
+    class is (generic_file_t)
+       call ft%set_counter(n)
+    end select
+
   end subroutine file_set_counter
 
 end module file