velocityProfileSemilog.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from .uniqueFinder import uniqueFinder
from .velocityProfile import velocityProfile
from .Rxx import Rxx
from .Rxy import Rxy
from .calcUstar import calcUstar
from .outlierRemove import outlierRemove
from .GuoProfileFunction import GuoProfileFunction
from matplotlib.ticker import LogLocator, LogFormatter
import matplotlib.ticker as ticker
import pickle

from cycler import cycler ## https://ranocha.de/blog/colors/
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))  ## https://ranocha.de/blog/colors/

markers = ['o', 's', 'D', '^', 'v', 'p', '*', 'h', '+', 'x']

def velocityProfileSemilog(dataPickle, positionList, plotTitle, subplot=False, roughWall=False, lenRemove=None):
    '''
    This function aims at producing the velocity profiles in the y+ vs u+ format.
    INPUT : dataPickle, positionList, plotTitle, subplot=False, roughWall=False
            lenRemove = [int] number of data in the dataList will be skipped in case of calculating the shear stress at the carpet. To use in combination with roughWall=True.
    '''

    dataPickle = dataPickle
    positionList = positionList
    plotTitle = plotTitle
    fz = 15 # fontsize

    with open(dataPickle, 'rb') as f:
        dataList = pickle.load(f)

    colors = plt.cm.tab10(np.linspace(0, 1, len(dataList))) # create colors for scatter plot

    #plt.rc("axes", prop_cycle=line_cycler)

    ## Theoretical profile
    yplus_theoretical1 =np.linspace(0,14);
    u_theoretical1 = yplus_theoretical1
    eq1 = (r"\begin{eqnarray*}"
           r"U^+ = y^+"
           r"\end{eqnarray*}")
    #axes.text(1,3, eq1, color="k", fontsize=16)

    kappa = 0.412
    #axes.text(100, 15, eq2, color="k", fontsize=16)
    #axes.plot(yplus_theoretical1, u_theoretical1, '-.k', linewidth=2) #label='$2.5\ln(y^+)+5.45$'
    #axes.plot(yplus_theoretical2, u_theoretical2, '-.k', linewidth=2)## for theoreitcal profile

    if subplot == False and roughWall == False:  #Option 1
        fig, axes = plt.subplots(1,1, layout='constrained', sharey=True, sharex=True)
        axes.set_xscale('log')
        for i, data in reversed(list(enumerate(dataList))):
            print('----Program Executing--',i,'---')
            Ustar = calcUstar(data) # function calling to get the shear velocity
            data = velocityProfile(data) # function calling to get the averaged velocity

            Uplus = data[0].reset_index()['U']/Ustar
            yplus = data[0].reset_index()['depth[cm]']/100*Ustar/1e-6   # 1e-6 is the kinematic viscosity of the water

            axes.scatter(yplus, Uplus, s=4,facecolors='k', marker=markers[i])

            axes.set_xlabel('$y^+$')
            axes.set_ylabel('$U^+$')
            axes.set_ylim(bottom=0)
            axes.grid()

        plt.savefig(plotTitle+'.pdf')

    if subplot == True and roughWall == True:   #Option 2 : length scaled by the roughness height ks
        kappa = 0.412
        ks = 0.01 # Characteristic height of the roughness elements. !!! if you change ks here, please change it in calcUstar.py as well!!!
        UList = [] # a list to store the shear velocity values (for plotting purpose) and finally to pick up the max; depth
        plt.cla() ## clear the current plot
        rel_x = 0.1  # relative x-position within the axes (0 to 1)
        rel_y = 0.3  # relative y-position within the axes (0 to 1)
        # parameter for Subplots #
        num_cols = 4
        # num_cols = min(len(dataList),5)
        lengthDataList = len(dataList)-lenRemove
        num_rows = -(-(lengthDataList)//num_cols) #
        fig, axes = plt.subplots(nrows=num_rows, ncols=num_cols, layout='constrained', sharey=True, sharex=True, figsize=(8, 3.25*num_rows))
        #fig.subplots_adjust(wspace=0.4) # bottom=0.1, right=0.95, top=0.9, left=0.075
        axes_flat = axes.flatten()

        file = open('U*'+plotTitle+'topWall.txt', 'w')  ## Create a file to store results
        file.write("Position\tU*\ty0/ks\tCf\n")

        #dataList = dataList[:8]
        for i, data in reversed(list(enumerate(dataList))):
            print('----Program Executing--',i,'---')
            row = i // num_cols ## commented for subplots array
            col = i % num_cols

            try:
                Ustar, B = calcUstar(data, vonKarman=True) # function calling to get the shear velocity
                print('U*=', Ustar) # U* is obtained by fitting the line, keeping kappa=0.412, for the flow near the carpet
                y0byks = 1/np.exp(kappa*B) # Calculate y0/ks

                data = velocityProfile(data) # function calling to get the averaged velocity
                Uplus = data[0].reset_index()['U'][::-1]/Ustar
                #yplus = (data[1]-data[0].reset_index()['depth[cm]'][::-1])/100*Ustar/1e-6
                yplus = (data[1]-1-data[0].reset_index()['depth[cm]'][::-1])/100/ks  # This is not y+ value actually. It in fact means y/ks. data[-1] is used for changing refernce frame
                                                                                # at the free surface instead of the flume bottom. 1 is subtracted because we measured from the half line of the carpet.
                UList.append(np.max(Uplus)) # storing Uplus values for all positions
                ## Theoretical profile
                yplus_theoretical = np.linspace(yplus.to_list()[0],yplus.to_list()[-1]);
                u_theoretical = 1/kappa*np.log(yplus_theoretical)+B

                Cf = 2*(Ustar/data[-1])**2 # Skin Friction coefficient; Ref; Wim's lecture notes - data[-1] is depth-averaged velocity

                file.write(f"{positionList[i]}\t{Ustar}\t{y0byks}\t{Cf}\n") ## Writing U* in the file

                if num_rows > 1:
                    axes[row, col].set_yscale('log')
                    axes[row, col].scatter(Uplus, yplus, facecolors='None', s=4, color='k') ## scatter plot measured data
                    axes[row, col].tick_params(axis='x', rotation=90)  # Rotate x-ticks vertically
                    axes[row, col].tick_params(axis='y', rotation=90)  # Rotate y-ticks horizontally
                    axes[row, col].set_title('at '+str(positionList[i])+' m', fontsize=fz)
                    axes[row, col].set_box_aspect(2) # aspect ratio for subplot
                    axes[row, col].set_xlim([np.max(UList), 0]) # limiting the x-lim using max; depth

                    eq = (r"\begin{eqnarray*}"
                        r"U^+ = \frac{1}{\kappa} {\ln} \frac{y}{k_s} + %.2f"
                        r"\end{eqnarray*}" %B)
                    axes[row, col].text(rel_x, rel_y, eq, fontsize=fz, rotation=90, transform=axes[row, col].transAxes)
                    axes[row, col].plot(u_theoretical, yplus_theoretical, '-.', color="gray" ) ## Profile for theoreitcal curve
                    axes[row, col].xaxis.set_major_locator(ticker.MultipleLocator(10))
                    axes[row, col].xaxis.set_minor_locator(ticker.MultipleLocator(2))
                    #axes[row, col].tick_params(axis='y', labelright=True)
                    axes[row, col].tick_params(which='both', direction='in')
                    axes[row, col].tick_params(which='minor', direction='in')

                if num_rows == 1:
                    axes[col].set_yscale('log')
                    axes[col].scatter(Uplus, yplus, facecolors='None', s=4, color='k') ## scatter plot measured data
                    axes[col].tick_params(axis='x', rotation=90)  # Rotate x-ticks vertically
                    axes[col].tick_params(axis='y', rotation=90)  # Rotate y-ticks horizontally
                    axes[col].set_title('at '+str(positionList[i])+' m', fontsize=fz)
                    axes[col].set_box_aspect(2) # aspect ratio for subplot
                    axes[col].set_xlim([np.max(UList), 0]) # limiting the x-lim using max; depth

                    eq = (r"\begin{eqnarray*}"
                        r"U^+ = \frac{1}{\kappa} {\ln} \frac{y}{k_s} + %.2f"
                        r"\end{eqnarray*}" %B)
                    axes[col].text(rel_x, rel_y, eq, fontsize=fz, rotation=90, transform=axes[col].transAxes)
                    axes[col].plot(u_theoretical, yplus_theoretical, '-.', color="gray" ) ## Profile for theoreitcal curve
                    axes[col].xaxis.set_major_locator(ticker.MultipleLocator(10))
                    axes[col].xaxis.set_minor_locator(ticker.MultipleLocator(2))
                    #axes[row, col].tick_params(axis='y', labelright=True)
                    axes[col].tick_params(which='both', direction='in')
                    axes[col].tick_params(which='minor', direction='in')
            except np.linalg.LinAlgError:
                # Handle the exception if polyfit fails
                print("Polynomial fitting failed. Skipping the polyfit.")
                #axes[row, col].cla()
                #fig.delaxes(axes[row, col])

        # fig.supxlabel('$U^+$', rotation='vertical')  # Rotate y-label horizontally
        #fig.text(0.96, 0.5, '$y^+$', va='center', rotation='vertical') ## this is y-label

        fig.supxlabel(r'$U^+$')  ## for scaling by ks
        #fig.supxlabel(r'$U^+, \frac{1}{%.2f}\ln \frac{y}{k_s}+%.2f $' % (kappa, B))  ## for scaling by ks
        fig.supylabel('$y/k_s$')
        #fig.supylabel('$y/k_s, k_s=%.3f$'%ks)
        #fig.supxlabel(r'$U^+, \frac{1}{%.2f}\ln y^+ +%.2f $' % (kappa, Cont)) # for scaling by kinematic viscosity (nu)
        #fig.supylabel('$y^+$')

        num_empty = num_cols * num_rows - lengthDataList ## Clean the extra plot
        for i in range(num_empty):
            ax = axes_flat[-1]
            ax.axis('off')
            axes_flat = axes_flat[:-1]

        file.close() ## Closing the file
        plt.savefig(plotTitle+'roughWall'+'.pdf')

    if subplot == True and roughWall == False:  #Option 3
        UList = [] # a list to store the shear velocity values (for plotting purpose) and finally to pick up the max; depth
        plt.cla() ## clear the current plot
        # parameter for Subplots #
        num_cols = 5
        rel_x = 0.5  # relative x-position within the axes (0 to 1)
        rel_y = 0.3  # relative y-position within the axes (0 to 1)
        # num_cols = min(len(dataList),5)
        num_rows = -(-len(dataList)//num_cols)
        fig, axes = plt.subplots(nrows=num_rows, ncols=num_cols, layout='constrained', sharey=True, sharex=True, figsize=(8, 3.25*num_rows))
        axes_flat = axes.flatten()

        file = open('U*'+plotTitle+'bottom.txt', 'w')  ## Create a file to store results
        file.write("Position\tU*\tCf\n")

        for i, data in reversed(list(enumerate(dataList))):
            print('----Program Executing--',i,'---')
            row = i // num_cols ## commented for subplots array
            col = i % num_cols
            Ustar, A = calcUstar(data, bottomWall=True, vonKarman=True) # function calling to get the shear velocity
            data = velocityProfile(data) # function calling to get the averaged velocity
            Cf = 2*(Ustar/data[-1])**2 # Skin Friction coefficient; Ref; Wim's lecture notes - data[-1] is depth-averaged velocity
            file.write(f"{positionList[i]}\t{Ustar}\t{Cf}\n") ## Writing U* in the file
            Uplus = data[0].reset_index()['U']/Ustar
            yplus = data[0].reset_index()['depth[cm]']/100*Ustar/1e-6   # 1e-6 is the kinematic viscosity of the water

            axes[row, col].set_yscale('log')
            axes[row, col].tick_params(axis='x', rotation=90)  # Rotate x-ticks vertically
            axes[row, col].tick_params(axis='y', rotation=90)  # Rotate y-ticks horizontally

            axes[row, col].scatter(Uplus, yplus, facecolors='None', s=4, color='k')
            ## Theoretical profile
            kappa = 0.412
            #Cont = 7.5 # Integration constant
            yplus_theoretical = np.linspace(yplus.to_list()[0],yplus.to_list()[-1]);
            u_theoretical = 1/kappa*np.log(yplus_theoretical)+A
            axes[row, col].plot(u_theoretical, yplus_theoretical, '-.', color="gray" )

            eq = (r"\begin{eqnarray*}"
                r"U^+ = \frac{1}{\kappa} {\ln} y^+ + %.2f"
                r"\end{eqnarray*}" %A)
            axes[row, col].text(rel_x, rel_y, eq, fontsize=fz, rotation=90, transform=axes[row, col].transAxes)

            axes[row, col].set_title('at '+str(positionList[i])+' m', fontsize=fz)
            axes[row, col].set_box_aspect(2) # aspect ratio for subplot

            UList.append(np.max(Uplus)) # storing Uplus values for all positions
            axes[row, col].set_xlim([np.max(UList), 0]) # limiting the x-lim using max; depth
            #axes[row, col].set_ylim([np.min(yplus), np.max(yplus)]) # limiting the y-lim using max; depth

            axes[row, col].xaxis.set_major_locator(ticker.MultipleLocator(10))
            axes[row, col].xaxis.set_minor_locator(ticker.MultipleLocator(2))

            #axes[row, col].tick_params(axis='y', labelright=True)
            axes[row, col].tick_params(which='both', direction='in')
            axes[row, col].tick_params(which='minor', direction='in')

            # axes[row, col].tick_params(which='both', direction='in', length=6, top=True, bottom=True, left=True, right=True)
            # axes[row, col].tick_params(which='minor', direction='in', length=4, top=True, bottom=True, left=True, right=True)

        # fig.supxlabel('$U^+$', rotation='vertical')  # Rotate y-label horizontally
        #fig.text(0.96, 0.5, '$y^+$', va='center', rotation='vertical') ## this is y-label
        #fig.supxlabel(r'$U^+, \frac{1}{%.2f}\ln y^+ +%.2f $' % (kappa, Cont))
        fig.supxlabel(r'$U^+$')
        fig.supylabel('$y^+$')

        num_empty = num_cols * num_rows - len(dataList) ## Clean the extra plot
        for i in range(num_empty):
            ax = axes_flat[-1]
            ax.axis('off')
            axes_flat = axes_flat[:-1]
        file.close() ## Closing the file
        plt.savefig(plotTitle+'bottomWall'+'.pdf')