Victor/enable xp

.gitignore

@@ -1,4 +1,6 @@
 *.mod
+*.txt
+*.case
 *.o
 *.smod
 *.f0*
@@ -10,3 +12,8 @@ fst_spectrum.dat
 *.chkp
 sphere.dat
 sphere.csv
+FST/example/fort.0
+FST/example/neko
+FST/example/usr_driver.f90
+tripping/neko
+FST/src/generate_FST

FST/README.md

@@ -15,28 +15,79 @@ The initialization, generation, application of the FST is driven by `07_fst_bc_d
 - `fst_bc_driver_apply()`
 - `fst_bc_driver_finalize()`
 
-An example of usage is given in the user file `example.f90`. Note that to apply the boundary condition we use the `field_dirichlet_update` function which
-requires the use of the `user_velocity` boundary condition on the desired boundary (see `example.case`).
+An example of usage is given in the user file `example/user.f90`. Note that to apply the boundary condition we use the `field_dirichlet_update` function which
+requires the use of the `user_velocity` boundary condition on the desired boundary (see `example/run.case`).
 
 ## Case file
 
-The driver module uses some parameters that should be given in the case file. Below is the JSON object taken from `example.case` that shows which parameters to use:
+The driver module uses some parameters that should be given in the case file. Below is the JSON object taken from `example/run.case` that shows which parameters to use:
 
+With file regeneration:
 ```.json
 "FST": {
       "enabled": true,   // default is true
       "t_start": 0.0001, // Time at which to start applying FST
       "t_ramp": 0.001,   // Length of the linear ramp in time
       "alpha": 0.2,      // see below for full explanation of what this is
-      "ystart": -0.01,  // Lower bound for the fringe function
-      "yend": 0.01,     // High bound for the fringe function
-      "periodic_z": true // Self-explanatory. If periodic in y add "periodic_y": true
+      "ystart": -0.01,   // Lower bound for the fringe function (Also exists for z, if y is periodic)
+      "yend": 0.01,      // High bound for the fringe function  (Also exists for z, if y is periodic)
+      "periodic_z": true, // Self-explanatory. If periodic in y add "periodic_y": true
+      "regen_files": true, // Set to true to generate wavenumbers etc. See below for further explanation
+      "Uinf": 1.0,       // Free-stream velocity. Only read if "regen_files" is false.      "fst_path": "src" // Path where the fst files should be written.
 }
 ```
 
+Without file regeneration (reuse previously written files)
+```.json
+"FST": {
+      "enabled": true,   // default is true
+      "t_start": 0.0001, // Time at which to start applying FST
+      "t_ramp": 0.001,   // Length of the linear ramp in time
+      "alpha": 0.2,      // see below for full explanation of what this is
+      "ystart": -0.01,   // Lower bound for the fringe function (Also exists for z, if y is periodic)
+      "yend": 0.01,      // High bound for the fringe function  (Also exists for z, if y is periodic)
+      "periodic_z": true, // Self-explanatory. If periodic in y add "periodic_y": true
+      "regen_files": false, // Set to true to generate wavenumbers etc. See below for further explanation
+      "Uinf": 1.0,       // Free-stream velocity. Only read if "regen_files" is false.
+      "fst_path": "src" // Path to the fst files from which to read.
+}
+```
+
+@note That `fst_path` is interpreted differently based on the value of 
+`regen_files`.
+
+### FST generation
+
+In the original implementation, FST wavenumbers and amplitudes are generated
+on-the-fly at the beginning of the simulation. This is the default behavior in
+the present implementation. Generating wavenumbers/amplitudes on-the-fly will
+create three separate files:
+- `bb.txt`, which contains the random phases,
+- `fst_spectrum.csv`, which contains wavenumbers, amplitudes and the random,
+  unitary, divergence-free vectors, and
+- `sphere.dat`, which contains information about # shells, points per shell,
+  and wavenumber discretization parameters.
+
+You also have the possibility to reuse previously generated files to keep
+the same FST parameters across two simulations. To do that, use the parameter
+`FST.regen_files` and set it to `false`. Be careful that the 3 files mentioned
+above must be present in the same folder as your executable. 
+
+The format of the fst_spectrum.csv file is a bit special and differs
+from previous implementations. This version of the code requires it to have 8
+columns, which are as follows:
+
+```.csv
+ShellNo,kx,ky,kz,amp,u_hat_pn1,u_hat_pn2,u_hat_pn3
+```
+
+If `regen_files` is `false`, you must also specify the free-stream velocity
+since nothing from `01_global_params.f90` will be used. This can be done via
+the parameters `FST.Uinf`.
+
 ### Spatial fringe parameters
 
-A smooth fringe function is applied on the 2D inlet plane, which at the moment is assumed to be `(y,z)`.
+A smoothing function in space is applied on the 2D inlet boundary.
 The shape of this fringe is the one used in SIMSON and by lots of other people:
 
 $$
@@ -51,9 +102,20 @@ $$
 
 Note that $\lambda_u = 1$ if the direction `u` is set to be periodic.
 
-`_start` and `_end` parameters need to be set by the user, which represent geometrical coordinates. 
-By default, and only if the direction is not periodic, `_start` will be set to the minimum coordinate on the boundary (in that direction). 
-The same goes for `_end`, it will be by default set to the maximum value.
+
+Below is an example of fringe function in 1 dimension.
+![An example of fringe function with different alphas](fringe.png "Example of fringe function")
+
+And here is an example of how it looks like on an inlet boundary. In this case, we are showing
+the z-component of velocity where the baseflow has 0 velocity. The z-direction is set to be 
+periodic and therefore there is no smoothing along the z-direction. The y-direction is not 
+periodic. We have exaggerated the extents of the fringe function to make it obvious, but
+this will very much be case dependent.
+![A practical example of an inlet boundary](fst.png "Example of z-component of FST on an inlet boundary")
+
+The `_start` and `_end` parameters can be set by the user. By default, and only if the direction is not 
+periodic, `_start`/`_end` will be set to the minimum/maximum coordinate on the boundary (in that direction). 
+
 The quantities $\delta_{u,*}$ are computed as a percentage $\alpha$ of the total boundary length 
 in the direction `u`: $\delta_{u,rise} = \delta_{u,fall} = \alpha * L_u$, where 
 $L_u$ is the total domain length at the inlet in the direction u. 

FST/example/makeneko_cpu

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+usage() {
+
+  echo "Usage"
+  echo "./makeneko_cpu <user file> <path to neko install> <path to src folder>"
+  echo ""
+  echo "Example"
+  echo "./makeneko_cpu swept-wing.f90 /cfs/klemming/scratch/y/ycl/Neko/neko ./src"
+
+}
+
+# Check if exactly 4 arguments are given
+if [ "$#" -ne 3 ]; then
+  echo "Error: Exactly 4 arguments are required."
+  usage
+  exit
+fi
+
+NEKO_INSTALL_PATH=${2}
+src_path=${3}
+
+backend="${src_path}/bcknd/cpu/opr_fst_cpu.f90 ${src_path}/bcknd/device/opr_fst_device.F90"
+interface="${src_path}/fst_operator.f90"
+fst_core="${src_path}/FST-core/0*"
+driver="${src_path}/drivers/0*"
+
+echo ${NEKO_INSTALL_PATH}/bin/makeneko $backend $interface $fst_core $driver $1
+${NEKO_INSTALL_PATH}/bin/makeneko $backend $interface $fst_core $driver $1

FST/example/makeneko_gpu

@@ -0,0 +1,79 @@
+#!/bin/bash
+
+usage() {
+
+  echo "Usage"
+  echo "./makeneko_gpu <user file> <path to neko install> <hip or cuda> <path to src folder>"
+  echo ""
+  echo "Example"
+  echo "./makeneko_gpu swept-wing.f90 /cfs/klemming/scratch/y/ycl/Neko/neko hip ./src"
+
+}
+
+NEKO_INSTALL_PATH=${2}
+which_gpu=${3}
+src_path=${4}
+
+# Check if exactly 4 arguments are given
+if [ "$#" -ne 4 ]; then
+  echo "Error: Exactly 4 arguments are required."
+  usage
+  exit
+fi
+
+# Check if $3 is 'cuda' or 'hip'
+if [ "$3" != "cuda" ] && [ "$3" != "hip" ]; then
+  echo "Error: Argument 3 must be 'cuda' or 'hip'."
+  usage
+fi
+
+echo "[INFO] Neko install path: $NEKO_INSTALL_PATH"
+echo "[INFO] GPU backend      : $which_gpu"
+echo "[INFO] User source path : $src_path"
+echo
+
+# First, copy makeneko to a dummy one
+cp "${NEKO_INSTALL_PATH}/bin/makeneko" dummy_makeneko
+
+# Add relevant flag depending on CUDA or HIP
+if [ "${which_gpu}" == "cuda" ]; then
+   HAVE_CUDA=$(grep "HAVE_CUDA_BCKND=" dummy_makeneko | cut -d'=' -f2)
+   ADD_CUSTOM_DEF="-DHAVE_CUDA=${HAVE_CUDA}"
+else
+   HAVE_HIP=$(grep "HAVE_HIP_BCKND=" dummy_makeneko | cut -d'=' -f2)
+   ADD_CUSTOM_DEF="-DHAVE_HIP=${HAVE_HIP}"
+fi
+
+# Modify the dummy makeneko and add the relevant variables HAVE_CUDA
+# and HAVE_HIP to the FCFLAGS
+FCFLAGS_CURRENT=$(grep "FCFLAGS=" dummy_makeneko | cut -d"'" -f2)
+FCFLAGS_NEW="${FCFLAGS_CURRENT} ${ADD_CUSTOM_DEF}"
+
+sed -i "s|^FCFLAGS=.*|FCFLAGS='${FCFLAGS_NEW}'|g" dummy_makeneko
+
+# Set the correct files to compile
+backend_path_gpu="${src_path}/bcknd/device"
+backend_path_cpu="${src_path}/bcknd/cpu"
+
+backend_cpu="$backend_path_cpu/opr_fst_cpu.f90"
+
+if [ "${which_gpu}" == "cuda" ]; then
+   backend_device="$backend_path_gpu/cuda/opr_fst.cu"
+else
+   backend_device="$backend_path_gpu/hip/opr_fst.hip"
+fi
+
+backend_interface="$backend_path_gpu/opr_fst_device.F90"
+
+backend="$backend_device $backend_cpu $backend_interface"
+interface="${src_path}/fst_operator.f90"
+fst_core="${src_path}/FST-core/0*"
+driver="${src_path}/drivers/0*"
+
+# Put here any extra files to compile
+extras="${src_path}/extras/0*"
+
+echo $NEKO_INSTALL_PATH/bin/makeneko $backend $interface $fst_core $driver $extras $1
+./dummy_makeneko $backend $interface $fst_core $driver $extras $1
+
+rm dummy_makeneko

FST/example/run.case

@@ -13,9 +13,8 @@
       "enabled": false
     },
     "time": {
-      "end_time": 1.5,
-      "variable_timestep": true,
-      "target_cfl": 0.5,
+      "end_time": 0.3,
+      "variable_timestep": false,
       "timestep": 1e-3,
       "max_timestep": 0.25e-2
     },
@@ -27,10 +26,14 @@
 
     "FST": {
       "enabled": true, // !!! FST IS DISABLED !!
-      "t_start": 0.3,
+      "t_start": 0.05,
       "t_ramp": 0.001,
-      "alpha": 0.3, // 15% distance to the y boundaries
-      "periodic_z": true
+      "alpha": 0.15, // 15% distance to the y boundaries
+      "periodic_z": true,
+      "y_start": -0.1,
+      "y_end": 0.1,
+      "regen_files": false,
+      "Uinf": 1.0
     },
 
     "fluid": {

FST/example/user.f90

@@ -5,7 +5,7 @@ module user
 
 contains
 
-  ! Register user defined functions (see user_intf.f90)
+! Register user defined functions (see user_intf.f90)
   subroutine user_setup(u)
     type(user_t), intent(inout) :: u
     u%initialize => initialize ! Initialize parameters
@@ -28,17 +28,16 @@ subroutine initialize(time)
     p => neko_registry%get_field("p")
 
     coef => neko_user_access%case%fluid%c_Xh
-
-    call fst_bc_driver_initialize(time%t,u,v,w,p,coef,neko_user_access%case%params) 
 
-  end subroutine initialize
+    call fst_bc_driver_initialize(time%t,u,v,w,p,coef,neko_user_access%case%params)
 
-  !> Set boundary conditions
+  end subroutine initialize
+
+!> Set boundary conditions
   subroutine user_bc(fields, bc, time)
     type(field_list_t), intent(inout) :: fields
     type(field_dirichlet_t), intent(in) :: bc
     type(time_state_t), intent(in) :: time
-
     type(field_t), pointer :: ui,vi,wi,pi
     type(coef_t)   , pointer :: coef
 
@@ -61,8 +60,8 @@ subroutine user_bc(fields, bc, time)
          ! At the first time step, apply velocity BC by copying
          ! what has been applied in the initial condition
          !
-         if (time%tstep .eq. 1) then    
-          
+         if (time%tstep .eq. 1) then
+
             ! Get solution fields (u,v,w) which contain the initial condition
             ! if we are at the first timestep 
             ui => neko_registry%get_field("u")
@@ -81,12 +80,14 @@ subroutine user_bc(fields, bc, time)
 
          end if
 
-
          !
          ! Apply FST
          !
          coef => neko_user_access%case%fluid%c_Xh
-         call fst_bc_driver_apply(u, v, w, bc, coef, time%t, time%tstep, 0.0_rp, .false.)
+         call fst_bc_driver_apply(u, v, w, bc, coef, time%t, time%tstep, &
+                 0.0_xp, & ! Angle in the x-y direction
+                 .false.)  ! Wether to do FST on CPU
+
        end associate
 
     !
@@ -95,10 +96,10 @@ subroutine user_bc(fields, bc, time)
     else if (trim(fields%items(1)%ptr%name) .eq. "p") then
 
        associate(p => fields%items(1)%ptr)
-         
-         if (time%tstep .eq. 1) then    
+
+         if (time%tstep .eq. 1) then
             pi => neko_registry%get_field("p")
-         
+
             ! Pressure can be handled directly on GPU
             if (neko_bcknd_device .eq. 1) then
                call device_masked_copy_0(p%x_d, pi%x_d, bc%msk_d, p%dof%size(), bc%msk(0))

FST/src/FST-core/01_global_params.f90

@@ -25,7 +25,7 @@ module global_params
 
   integer :: k_length
   integer :: shell(fst_modes)
-  integer :: seed
+  !integer :: seed
 
   real(kind=rp) :: k_num(fst_modes, 3)
   real(kind=rp) :: k_num_all(fst_modes, 3)
@@ -44,7 +44,7 @@ subroutine print_param(name, value)
     real(kind=rp), intent(in) :: value
     character(len=LOG_SIZE) :: log_buf
 
-    write(log_buf, *) name, ": ", value
+    write(log_buf, *) "[FST] ", name, ": ", value
     call neko_log%message(log_buf)
 
   end subroutine print_param
@@ -57,7 +57,7 @@ subroutine print_int(name, value)
     integer, intent(in) :: value
     character(len=LOG_SIZE) :: log_buf
 
-    write(log_buf, *) name, ": ", value
+    write(log_buf, *) "[FST] ", name, ": ", value
     call neko_log%message(log_buf)
 
   end subroutine print_int

FST/src/FST-core/02_sphere.f90

@@ -1,7 +1,6 @@
 module sphere
   use neko
   use global_params
-  ! use global_params
   implicit none
 
 contains