Added example and compilation instructions

examples/lid_insitu/README.md

@@ -0,0 +1,6 @@
+## Simulation of a lid-driven cavity
+In this case we simulate a lid-driven cavity with smoothened belt velocity to fulfil the continuity equaiton in the corners. There is a 3D example (with 3 elements in the spanwise $z$ direction) and a pseudo-2D example (with a 2D mesh that then generates a case with one element in the spanwise direction).
+
+The `lid.box` file needs to be converted to a .re2 file using the Nek5000 tool genbox, and then using rea2nbin into the `lid.nmsh` Neko mesh file.
+
+The cavity flow is stable (steady) up to a Reynolds number of about 8000. The mesh size may be changed in the box file.

examples/lid_insitu/insitu_task.py

@@ -0,0 +1,53 @@
+# Import general modules
+import sys
+import os
+import json
+import numpy as np
+
+# Import MPI
+from mpi4py import MPI #equivalent to the use of MPI_init() in C
+
+# Split communicator for MPI - MPMD
+worldcomm = MPI.COMM_WORLD
+worldrank = worldcomm.Get_rank()
+worldsize = worldcomm.Get_size()
+col = 1
+comm = worldcomm.Split(col,worldrank)
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+import adios2
+
+# Import data streamer
+from pynektools.io.adios2.stream import DataStreamer
+from pynektools.io.utils import get_fld_from_ndarray
+
+# Instance the streamer (start the stream)
+ds = DataStreamer(comm)
+
+# Recieve the data from fortran
+x = get_fld_from_ndarray(ds.recieve(), ds.lx, ds.ly, ds.lz, ds.nelv) # Recieve and reshape x
+y = get_fld_from_ndarray(ds.recieve(), ds.lx, ds.ly, ds.lz, ds.nelv) # Recieve and reshape y
+z = get_fld_from_ndarray(ds.recieve(), ds.lx, ds.ly, ds.lz, ds.nelv) # Recieve and reshape z
+
+# Finalize the stream 
+ds.finalize()
+
+# Now that the data is here. Create a mesh object, an empty field and write it to disk
+# For this, the pyNek tools are needed
+from pynektools.io.ppymech.neksuite import pwritenek
+from pynektools.datatypes.msh import Mesh
+from pynektools.datatypes.field import Field
+from pynektools.datatypes.utils import create_hexadata_from_msh_fld
+
+## Instance the mesh
+msh = Mesh(comm, x = x, y = y, z = z)
+## Create an empty field and update its metadata
+out_fld = Field(comm)
+out_fld.fields["scal"].append(np.ones_like(msh.x))
+out_fld.update_vars()
+## Create the hexadata to write out
+out_data = create_hexadata_from_msh_fld(msh = msh, fld = out_fld)
+## Write out a file
+fname = "python_test0.f"+str(1).zfill(5)
+pwritenek("./"+fname,out_data, comm)

examples/lid_insitu/lid.box

@@ -0,0 +1,9 @@
+orig.rea
+ -3                      spatial dimension
+ 1                     number of fields
+Box
+6     6     3       nelx,nely,nelz for Box)
+0. 0.06 0.25 0.5 0.75 0.94 1.   	  x0 x1 ratio 
+0. 0.06 0.25 0.5 0.75 0.94 1.     	  y0 y1 ratio 
+0. 1. 2. 3.
+W  ,W  ,W  ,v  ,P  ,P  ,    V bc's  ! NB:  3 characters each !

examples/lid_insitu/lid.case

@@ -0,0 +1,48 @@
+{
+  "version": 1.0,
+  "case":
+  {
+  "mesh_file": "lid.nmsh",
+  "end_time": 0.01,
+  "timestep": 1e-3,
+  "numerics": {
+    "time_order": 3,
+    "polynomial_order": 7,
+    "dealias": true
+  },
+  "fluid": {
+    "scheme": "pnpn",
+    "Re": 5000,
+    "inflow_condition": {
+      "type": "user"
+    },
+    "initial_condition": {
+      "type": "user"
+    },
+    "velocity_solver": {
+      "type": "cg",
+      "preconditioner": "jacobi",
+      "absolute_tolerance": 1e-7,
+      "max_iterations": 200
+    },
+    "pressure_solver": {
+      "type": "gmres",
+      "preconditioner": "hsmg",
+      "projection_space_size": 8,
+      "absolute_tolerance": 1e-7,
+      "max_iterations": 200
+    },
+    "output_control": "simulationtime",
+    "output_value": 5,
+    "boundary_types": ["w", "v"]
+  },
+  "simulation_components":
+  [
+    {
+      "type": "vorticity",
+      "compute_control": "tsteps",
+      "compute_value": 50
+    }
+  ]
+  }
+}

examples/lid_insitu/lid.f90

@@ -0,0 +1,179 @@
+! Lid-driven cavity
+!
+! Time-integration of the lid-driven cavity with smoothened
+! belt velocity to fulfil continuity equation.
+!
+module user
+  use neko
+  implicit none
+
+  type(data_streamer_t) :: dstream
+
+  ! Global user variables
+  type(field_t) :: w1
+
+  type(file_t) output_file ! output file
+  type(vector_t) :: vec_out    ! will store our output data
+
+ contains
+
+  ! Register user-defined functions (see user_intf.f90)
+  subroutine user_setup(user)
+    type(user_t), intent(inout) :: user
+    user%fluid_user_ic => user_ic
+    user%fluid_user_if => user_bc
+    user%user_check => user_calc_quantities
+    user%user_init_modules => user_initialize
+    user%user_finalize_modules => user_finalize
+  end subroutine user_setup
+
+  ! user-defined boundary condition
+  subroutine user_bc(u, v, w, x, y, z, nx, ny, nz, ix, iy, iz, ie, t, tstep)
+    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), intent(in) :: t
+    integer, intent(in) :: tstep
+
+    real(kind=rp) lsmoothing
+    lsmoothing = 0.05_rp    ! length scale of smoothing at the edges
+
+    u = step( x/lsmoothing ) * step( (1._rp-x)/lsmoothing )
+    v = 0._rp
+    w = 0._rp
+
+  end subroutine user_bc
+
+  ! 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(json_file), intent(inout) :: params
+    u = 0._rp
+    v = 0._rp
+    w = 0._rp
+    p = 0._rp
+  end subroutine user_ic
+
+  ! User-defined initialization called just before time loop starts
+  subroutine user_initialize(t, u, v, w, p, coef, params)
+    real(kind=rp) :: t
+    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(coef_t), intent(inout) :: coef
+    type(json_file), intent(inout) :: params
+
+    integer tstep
+    ! Initialize the file object and create the output.csv file
+    ! in the working directory
+    output_file = file_init("output.csv")
+    call output_file%set_header("t,Ekin,enst")
+    call vec_out%init(2) ! Initialize our vector with 2 elements (Ekin, enst)
+
+    ! initialize work arrays for postprocessing
+    call w1%init(u%dof, 'work1')
+
+    ! call usercheck also for tstep=0
+    tstep = 0
+    call user_calc_quantities(t, tstep, u, v, w, p, coef, params)
+
+    ! Initialize the streamer
+    call dstream%init(coef)
+
+    ! Stream the mesh
+    call dstream%stream(coef%dof%x)
+    call dstream%stream(coef%dof%y)
+    call dstream%stream(coef%dof%z)
+
+  end subroutine user_initialize
+
+  ! User-defined routine called at the end of every time step
+  subroutine user_calc_quantities(t, tstep, u, v, w, p, coef, params)
+    real(kind=rp), intent(in) :: t
+    integer, intent(in) :: tstep
+    type(coef_t), intent(inout) :: coef
+    type(json_file), intent(inout) :: 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(field_t), pointer :: om1, om2, om3
+    integer :: ntot, i
+    real(kind=rp) :: e1, e2
+
+    if (mod(tstep,50).ne.0) return
+
+    om1 => neko_field_registry%get_field("omega_x")
+    om2 => neko_field_registry%get_field("omega_y")
+    om3 => neko_field_registry%get_field("omega_z")
+
+    ntot = u%dof%size()
+
+    call col3(w1%x,u%x,u%x,ntot)
+    call addcol3(w1%x,v%x,v%x,ntot)
+    call addcol3(w1%x,w%x,w%x,ntot)
+    e1 = 0.5 * glsc2(w1%x,coef%B,ntot) / coef%volume
+
+    call col3(w1%x,om1%x,om1%x,ntot)
+    call addcol3(w1%x,om2%x,om2%x,ntot)
+    call addcol3(w1%x,om3%x,om3%x,ntot)
+    e2 = 0.5 * glsc2(w1%x,coef%B,ntot) / coef%volume
+
+    ! Output all this to file
+    call neko_log%message("Writing csv file")
+    vec_out%x = (/e1, e2/)
+    call output_file%write(vec_out, t)
+
+  end subroutine user_calc_quantities
+
+  ! Smooth step function
+  function step(x)
+  ! Taken from Simson code
+  ! x<=0 : step(x) = 0
+  ! x>=1 : step(x) = 1
+    real(kind=rp) :: step, x
+
+    if (x.le.0.02_rp) then
+       step = 0.0_rp
+    else
+       if (x.le.0.98_rp) then
+           step = 1._rp/( 1._rp + exp(1._rp/(x-1._rp) + 1._rp/x) )
+       else
+           step = 1._rp
+       end if
+    end if
+
+  end function step
+
+  ! User-defined finalization routine called at the end of the simulation
+  subroutine user_finalize(t, params)
+    real(kind=rp) :: t
+    type(json_file), intent(inout) :: params
+
+    ! Finalize the stream
+    call dstream%free()
+
+    ! Deallocate the fields
+    call w1%free()
+
+    ! Deallocate output file and vector
+    call file_free(output_file)
+    call vec_out%free
+
+  end subroutine user_finalize
+
+end module user

examples/lid_insitu/lid2d.box

@@ -0,0 +1,8 @@
+orig.rea
+ -2                      spatial dimension
+ 1                     number of fields
+Box
+6     6     3       nelx,nely,nelz for Box)
+0. 0.06 0.25 0.5 0.75 0.94 1.   	  x0 x1 ratio 
+0. 0.06 0.25 0.5 0.75 0.94 1.     	  y0 y1 ratio 
+W  ,W  ,W  ,v  ,    V bc's  ! NB:  3 characters each !

examples/lid_insitu/lid2d.case

@@ -0,0 +1,48 @@
+{
+  "version": 1.0,
+  "case":
+  {
+  "mesh_file": "lid2d.nmsh",
+  "end_time": 50,
+  "timestep": 1e-3,
+  "numerics": {
+    "time_order": 3,
+    "polynomial_order": 7,
+    "dealias": true
+  },
+  "fluid": {
+    "scheme": "pnpn",
+    "Re": 5000,
+    "inflow_condition": {
+      "type": "user"
+    },
+    "initial_condition": {
+      "type": "user"
+    },
+    "velocity_solver": {
+      "type": "cg",
+      "preconditioner": "jacobi",
+      "absolute_tolerance": 1e-7,
+      "max_iterations": 200
+    },
+    "pressure_solver": {
+      "type": "gmres",
+      "preconditioner": "hsmg",
+      "projection_space_size": 8,
+      "absolute_tolerance": 1e-7,
+      "max_iterations": 200,
+      "boundary_types": ["w", "v"]
+    },
+    "output_control": "simulationtime",
+    "output_value": 5
+  },
+  "simulation_components":
+  [
+    {
+      "type": "vorticity",
+      "compute_control": "tsteps",
+      "compute_value": 50
+    }
+  ]
+  }
+}