Final fixes to AM and Coriolis body forces

body_forces/am_f2/am_f2.inc

@@ -0,0 +1,25 @@
+! Suppresses superposition of large scales at the wall
+! Parameters are hardcoded for Poiseuille MSU, retau 1000
+
+subroutine config_body_force()
+    if (has_terminal) print *, "Using bodyforce, am_f2 (suppression of modulation)"
+    if (has_terminal) print *, ""
+    call set_body_force()
+end subroutine
+
+subroutine set_body_force()
+integer(C_INT) :: ix,iz,iy,iV
+real(C_DOUBLE) :: yp
+    do iV=1,3
+        do ix=nx0,nxN
+            do iz=-nz,nz
+                do iy=ny0-2,nyn+2
+                    yp = MERGE(ymax-y(iy), y(iy), y(iy)>1)*1000
+                    if ( (yp <= 40) .AND. (ABS(iz) >= 7) .AND. (MOD(ABS(iz),7) /= 0) ) then 
+                        F(iy,iz,ix,iV)=-amp*V(iy,iz,ix,iV)
+                    end if
+                end do
+            end do
+        end do
+    end do
+end subroutine

body_forces/am_f2/am_pardec.inc

@@ -0,0 +1 @@
+real(C_DOUBLE) :: amp=1000

body_forces/coriolis/coriolis.in

@@ -1,3 +1,3 @@
-0.02    ! 2*Ro(tation number)
+0.02    ! -Ro(tation number)
 9999999 ! kz_cutoff (for large/small scale discrimination)
 1.0     ! y_threshold (defines maximum y, wrt bottom plate, at which force is active)

body_forces/coriolis/coriolis.inc

@@ -14,7 +14,7 @@ subroutine config_body_force()
     y_threshold_top = ymax - y_threshold_bot
     iz_thr = min(nz,floor(kz_cutoff/beta0))
     ! print out parameters
-    if (has_terminal) print *, "Ro(tation number)", omega2/2
+    if (has_terminal) print *, "Ro(tation number)", -omega2
     if (has_terminal) print *, "kz_cutoff, iz_cutoff", kz_cutoff, iz_thr
     if (has_terminal) print *, "y_threshold_bottom", y_threshold_bot
     if (has_terminal) print *, ""

body_forces/coriolis/coriolis_pardec.inc

@@ -1,4 +1,4 @@
-real(C_DOUBLE) :: omega2 ! 2 * rotation number
+real(C_DOUBLE) :: omega2 ! minus one times rotation number (that's right: it's -Ro)
 real(C_DOUBLE) :: kz_cutoff ! for large/small scale discrimination
-real(C_DOUBLE) :: y_threshold_top, y_threshold_bot ! for large/small scale discrimination
+real(C_DOUBLE) :: y_threshold_top, y_threshold_bot ! for confining Coriolis force to a wall-normal region
 integer :: iz_thr

postpro/tke/uiuj_spectra.f90

@@ -336,8 +336,6 @@ program uiuj_largesmall
         close(15)
     end if
 
-    ! FIN QUI
-
     ! write to disk
     currfname = "uiuj_spectra.bin"
     call MPI_File_open(MPI_COMM_WORLD, trim(currfname), IOR(MPI_MODE_WRONLY, MPI_MODE_CREATE), MPI_INFO_NULL, fh)
@@ -357,10 +355,10 @@ program uiuj_largesmall
         CALL MPI_File_set_view(fh, offset, MPI_DOUBLE_PRECISION, profiles_write_type, 'native', MPI_INFO_NULL)
         CALL MPI_File_write_all(fh, uiujprofiles, 1, profiles_inmem_type, MPI_STATUS_IGNORE)
 
-        ! write convs data
-        offset = offset + ( int(ny+3,MPI_OFFSET_KIND) * 60_MPI_OFFSET_KIND * sizeof(uiujprofiles(1,1,1)) ) ! DON'T DO SIZEOF(uiujprofiles)! Program is PARALLEL!!!
-        CALL MPI_File_set_view(fh, offset, MPI_DOUBLE_PRECISION, convs_write_type, 'native', MPI_INFO_NULL)
-        CALL MPI_File_write_all(fh, convs, 1, convs_inmem_type, MPI_STATUS_IGNORE)
+        ! write convs data - NOT IMPLEMENTED YET!
+        !offset = offset + ( int(ny+3,MPI_OFFSET_KIND) * 60_MPI_OFFSET_KIND * sizeof(uiujprofiles(1,1,1)) ) ! DON'T DO SIZEOF(uiujprofiles)! Program is PARALLEL!!!
+        !CALL MPI_File_set_view(fh, offset, MPI_DOUBLE_PRECISION, convs_write_type, 'native', MPI_INFO_NULL)
+        !CALL MPI_File_write_all(fh, convs, 1, convs_inmem_type, MPI_STATUS_IGNORE)
 
     call MPI_File_close(fh)
 
@@ -545,7 +543,7 @@ subroutine print_help()
             print *, ""
             print *, "This program is meant to be used on plane channels."
             print *, ""
-            print *, "Results are output to uiuj_spectra/uiuj_spectra.bin."
+            print *, "Results are output to uiuj_spectra.bin. Input (dns.in) is read from parent folder."
             print *, ""
             print *, "Mean TKE budget terms are calculated as:"
             print *, "INST    --> dK/dt"

utilities/compress_uiuj_spectra.py

@@ -0,0 +1,28 @@
+from argparse import ArgumentParser
+import channel as ch
+import numpy as np
+import os
+import gc
+
+parser = ArgumentParser(prog='compress_uiuj_spectra', description='Compresses uiuj_spectra for online distribution.')
+parser.add_argument('path_to_uiuj_spectra', type=str, nargs=1)
+stngs = parser.parse_args()
+
+bin_file = stngs.path_to_uiuj_spectra[0]
+dnsin_file = bin_file.replace('uiuj_spectra.bin','../dns.in')
+
+m = ch.mesh(ch.read_dnsin(dnsin_file))
+
+new_size =  ((m.ny+3)*7) + ((m.ny+3)*((2*m.nz)+1)*10*6) + ((m.ny+3)*10*6) 
+old_size =  ((m.ny+3)*7) + ((m.ny+3)*((2*m.nz)+1)*10*6) + ((m.ny+3)*10*6) + (4*((2*m.nz)+1)*(m.nz+1)*(m.ny+3))
+
+old_size_bytes = os.stat(bin_file).st_size
+print(f'Old size in bytes: {old_size_bytes}')
+print(f'Expected size in bytes: {old_size*8}')
+if not (old_size_bytes == (old_size*8)):
+    raise Exception('input file is unexpectedly large.')
+
+old_file = np.memmap(bin_file,  dtype=np.float64, mode='r', offset=0, shape=(old_size))
+new_file = np.memmap(bin_file.replace('uiuj_spectra.bin','compressed_uiuj_spectra.bin'), dtype=np.float64, mode='w+', shape=(new_size))
+
+new_file[:] = old_file[:new_size]