2-D transition case: added the convergence history (#22)

nalu-wind/2D_airfoil_Transition/NLF1-0416/README.md

@@ -53,6 +53,13 @@ Option 2 is activated only if fsti is explicitly specified in the Nalu-Wind inpu
 Based on the grid sensitivity results, a full sweep of angles of attack was performed using the Fine mesh level. The two figures above compare the lift and drag polar with the experimental measurements[^5]. For the lift, the transition simulation slightly over-predicts the lift coefficient in the linear range of the lift curve, a similar behavior also observed in transition predictions using other transition models and other flow solvers. For the drag polar, the transition simulation predicts lower drag across the range of angles of attack than the fully turbulent simulation and better compares with the experimental data. 
 Specifically, the errors in the predicted drag coefficient at AoA=5° are 2.87% for the transition simulation and 57.56% for the fully turbulent simualtion.
 
+## Results: Convergence
+
+### Time history of the residuals and aerodynamic coefficients
+<img src="aoa_5/figures_and_scripts/time_history.png" alt="Cf" width="1000">
+
+Convergence behaviors of the transition simulation are presented in the above figure at AoA=5°. Left figure shows the time history of the non-linear residuals from Nalu-Wind for the mean flow (momentum), turbulence model (𝑘), and transition model (𝛾) at the last Picard iteration of each main iteration. In the figure, all residuals decline smoothly by 3.5 to 4 orders of magnitude until they finally stall. Similary, both the lift and drag converge well without oscillations, with less than 0.1% difference from the converged value around iteration 4,000.
+
 Each case with the "Fine" mesh took approximately 40 minutes to 10,000 iterations, using 4 Picard iterations per time step, on 26 cores of NREL's Kestrel HPC cluster. The number of cores per case was not determined by Nalu-Wind’s scalability on Kestrel, but simply to accommodate 4 cases on a single node of Kestrel.
 
 ## References

nalu-wind/2D_airfoil_Transition/NLF1-0416/aoa_5/figures_and_scripts/plot_time_hist_clcd.py

@@ -0,0 +1,84 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.backends.backend_pdf import PdfPages
+import glob, pathlib
+import pandas as pd
+import math
+from mpl_toolkits.mplot3d import Axes3D
+
+
+# Figure size format
+size=15
+params = {'legend.fontsize': size-2,
+          'axes.labelsize': size,
+          'axes.titlesize': size,
+          'xtick.labelsize': size,
+          'ytick.labelsize': size}
+plt.rcParams.update(params)
+
+# Forc outpute file from Nalu-Wind
+force_dir='../performance/forces.dat'
+
+# Flow conditions 
+u_infty = 341*0.1
+rho     = 1.225
+dt      = 1/(u_infty*2.0)
+alpha   = 5.0
+
+
+# Read the force file
+data = pd.read_csv(force_dir,sep="\s+",skiprows=1,
+    names=[
+        "Time",
+        "Fpx",
+        "Fpy",
+        "Fpz",
+        "Fvx",
+        "Fvy",
+        "Fvz",
+        "Mtx",
+        "Mty",
+        "Mtz",
+        "Y+min",
+        "Y+max",
+    ],
+)
+
+# Non-dimensionalization
+dyn_pres = 0.5 * rho * (u_infty ** 2)
+cfz = (data["Fpy"] + data["Fvy"]) / dyn_pres
+cfx = (data["Fpx"] + data["Fvx"]) / dyn_pres
+
+cos_aoa=math.cos(math.radians(alpha))
+sin_aoa=math.sin(math.radians(alpha))
+
+cl = cfz*cos_aoa - cfx*sin_aoa
+cd = cfz*sin_aoa + cfx*cos_aoa
+
+# Plot time time history
+fig, ax1 = plt.subplots(1)
+
+# Adjust plot legend
+cl_l=cl.iloc[-1]*0.985
+cl_u=cl.iloc[-1]*1.015
+
+cd_l=cd.iloc[-1]*0.99
+cd_u=cd.iloc[-1]*1.01
+
+x= data['Time']/dt
+ax2 = ax1.twinx()
+ax1.plot(x, cl, 'g')
+ax2.plot(x, cd, 'b')
+ax1.set_xlabel('Iteration')
+ax1.set_ylabel('Lift coefficient, $C_{l}$', color='g')
+ax2.set_ylabel('Drag coefficient, $C_{d}$', color='b')
+ax1.set_xlim(0, 10000)
+ax1.set_ylim(cl_l, cl_u)
+ax2.set_ylim(cd_l, cd_u)
+ax2.set_ylim(0.0075,0.0082)
+ax1.tick_params(axis='y', colors='g')
+ax2.tick_params(axis='y', colors='b')
+plt.tight_layout()
+plt.savefig("time_clcd_nlf0416_F_aoa_5.png",dpi=300)
+plt.show()
+

nalu-wind/2D_airfoil_Transition/NLF1-0416/aoa_5/figures_and_scripts/plot_time_hist_resid.py

@@ -0,0 +1,114 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.backends.backend_pdf import PdfPages
+import glob, pathlib
+import pandas as pd
+from mpl_toolkits.mplot3d import Axes3D
+
+
+# Plot format
+size=15
+params = {'legend.fontsize': size-2,
+          'axes.labelsize': size,
+          'axes.titlesize': size,
+          'xtick.labelsize': size,
+          'ytick.labelsize': size}
+plt.rcParams.update(params)
+
+# Nalu-Wind input
+trans   = 1 # transition model on
+log_file = "nlf0416_F_aoa_5.0.log" # Nalu-Wind log file
+
+
+def read_resid(log_file,niter,line_num):
+
+    # Find residual values
+    dum=[]
+    file=open(log_file,'r')
+    for lines in file:
+       dum += [lines.split()]
+    file.close()
+    n_line=int(len(dum))
+
+    # Save residuals in arrays
+    resid_mom_eqs_x = np.zeros((niter))
+    resid_mom_eqs_y = np.zeros((niter))
+    resid_mom_eqs_z = np.zeros((niter))
+    resid_mom_eqs = np.zeros((niter))
+    resid_rho_eqs = np.zeros((niter))
+    resid_k_eqs = np.zeros((niter))
+    resid_w_eqs = np.zeros((niter))
+    it = np.zeros((niter))
+    resid_gamma_eqs = np.zeros((niter))
+
+    nlines=9
+    if(trans==1):
+        nlines=10
+
+    for j in range(niter):
+        i0 = line_num[j+2]
+
+        # Read time step count number
+        ii = i0+1
+        it[j] = dum[ii][3]
+
+        dum_lines = 8
+        n_solver = 4
+        i  =i0+ dum_lines + (n_solver-1)*nlines
+        resid_mom_eqs[j]   = dum[i+5][3]
+        resid_rho_eqs[j]   = dum[i+6][3]
+        resid_k_eqs[j]     = dum[i+8][3]
+        resid_w_eqs[j]     = dum[i+9][3]
+        if(trans==1):
+           resid_gamma_eqs[j]     = dum[i+10][3]
+
+    return it,resid_mom_eqs,resid_rho_eqs,resid_k_eqs,resid_w_eqs,resid_gamma_eqs
+
+def find_index_location(log_file):
+    # Find line numbers
+    with open(log_file, 'r') as file:
+        lines = file.readlines()
+    file.close()
+
+    n_line=len(lines)
+
+    line_num=[]
+    for i in range(n_line):
+        if(lines[i][0]=="*"):
+            line_num.append(i)
+    del lines
+
+    niter=len(line_num)-4
+    print(niter)
+
+    return niter, line_num
+
+
+def create_nalu_hist(niter,line_num,log_file):
+    print(log_file)
+
+    [it,resid_mom_eqs,resid_rho_eqs,resid_k_eqs,resid_w_eqs,resid_gamma]=read_resid(log_file,niter,line_num)
+
+    plt.figure
+
+    plt.semilogy(it, resid_mom_eqs)
+    plt.semilogy(it, resid_k_eqs,'-.')
+    plt.semilogy(it, resid_gamma,'--')
+    plt.xlabel('Iteration')
+    plt.ylabel('Residuals')
+    plt.xlim([0, 8000])
+    plt.xticks(np.arange(0, 10000.1, 2000))
+    plt.ylim([1e-7, 1e-1])
+    #plt.legend(["Mean flow","Turbulence model:k","Turbuelnce model: omega","Transition model"])
+    plt.legend(["Mean flow","Turbulence model","Transition model"])
+    plt.tight_layout()
+    plt.savefig("resid_nlf0416_F_aoa_5.png",dpi=300)
+    plt.show()
+
+    del it,resid_mom_eqs,resid_rho_eqs,resid_k_eqs,resid_w_eqs
+
+if __name__=="__main__":
+
+    # Nalu-Wind log file
+    [niter, line_num] = find_index_location(log_file)
+    create_nalu_hist(niter,line_num,log_file)

nalu-wind/2D_airfoil_Transition/README.md

@@ -8,4 +8,5 @@ They include the following airfoils and flow parameters:
 - NASA NLF(1)-0416 airfoil at M=0.1, Re=4M, Tu=0.15%
 - TU Delft DU00-W-212 airfoil at M=0.1, Re=3M, Tu=0.0864%
 
-A paper is currently in draft containing these results.
+Read the NLF(1)-0416 case first, then the DU00-W-212 case.
+A paper containing detailed validation is currently under review.