NREL5MW Rigid Benchmark Updates

exawind/NREL_5MW_Turbine/rigid/README.md

@@ -1,74 +1,110 @@
 <!-- This file is automatically compiled into the website. Please copy linked files into .website_src/ paths to enable website rendering -->
 
-<style>
-    img,
-    .whitebackground img {
-        color: #ffffff !important;
-    }
-
-    h1,h2,h3 {
-        margin-top: 30px;
-    }
-</style>
-
 # NREL 5MW Rotor in Atmospheric Boundary Layer With Rigid Bodies and Constant Rotation
 This benchmark contains a geometry and boundary-layer resolved model of the NREL 5MW reference turbine [1] in a developed near-neutral convective atmospheric boundary layer with a mean hub-height streamwise velocity of 11.4m/s.
 
 - Hub-height: 90m
 - Blade Radius: 63m
 - Rated Power: 5 MW
 - Type: Upwind 3 Blade
+- Simulation dt: 0.003443526170799 s
+- OpenFAST dt: 0.0008608815426997245 s
 
 <i>Note on run-time: Current results are for 125s simulation time, but by the time of benchmark release, a minimum of 360s will be posted.</i>
 
 ## Simulation Setup
 - ExaWind driver SHA: [cba5259fc43ddeca67329630d3c84faad90e91bb](https://github.com/Exawind/exawind-driver/commit/cba5259fc43ddeca67329630d3c84faad90e91bb)
 - Nalu-Wind SHA: [b9e4ae654b646ecd0501dd6391dc7537239c82db](https://github.com/Exawind/nalu-wind/commit/b9e4ae654b646ecd0501dd6391dc7537239c82db)
-- AMR-Wind SHA: [091b07fa3840af98925792290fb8788696290a5e](https://github.com/Exawind/amr-wind/commit/091b07fa3840af98925792290fb8788696290a5e) 
+- AMR-Wind SHA: [091b07fa3840af98925792290fb8788696290a5e](https://github.com/Exawind/amr-wind/commit/091b07fa3840af98925792290fb8788696290a5e)
 
 ## Freestream Conditions
 A near-neutral convective boundary layer precursor was run to feed the domain for this case. Full details are posted here: [Convective ABL for NREL5MW](../../../amr-wind/atmospheric_boundary_layer/convective_abl_nrel5mw/README.md)
 
 10 minute flow statistics:
-- Hub-height streamwise velocity: 11.63 m/s
+- Hub-height streamwise velocity: 11.4 m/s
 - Hub-height flow direction: ~240.0 degrees
 
-<img src="../../../amr-wind/atmospheric_boundary_layer/convective_abl_nrel5mw/results/avgmetmast_0600.png" alt="10-Minute Average Virtual Met-mast Precursor Results" style="width:800px; background-color: #ffffff;"/>
+![avgmetmest60](../../../amr-wind/atmospheric_boundary_layer/convective_abl_nrel5mw/results/avgmetmast_0600.png)
 
 ## CFD Mesh
 
-<img src="figures_and_scripts/mesh-nrel5mw-abl.png" alt="Mesh layout" style="width:800px; background-color: #ffffff;"/>
+![meshpicture](figures_and_scripts/split_mesh_layout.png)
+**Total number of cells: 155,363,046**
 
-**Total number of cells: 177,729,426**
-
-The near-body (Nalu-Wind) mesh was created using PGL, pyHyp, and Pointwise. Three blades are connected at the hub, and a tower is included.
-- Structured hex mesh
-- Cell count: 15,593,362
+The near-body (Nalu-Wind) mesh was created using a proprietary surface mesher, pyHyp, and Pointwise. Three blades are split at the hub, and an unconnected tower is included.
+- Hexahedral cells
+- Cell count: 13,436,646
 
 Off-body (AMR-Wind) mesh was generated using the built-in capability of AMR-Wind. Off-body mesh information is summarized below 
-- Mesh topology: Cartesian with AMR
+- Mesh topology: structured hex with nested refinements
 - Domain in x= 0 to 5000m, y=0m to 5000m, z=0m to 1000m
 - Initial grid size: 10m
 - Finest cell size: 0.625m with 4 AMR levels
-- Cell count: 162,136,064
+- Cell count: 141,926,400
 
 ## Results
 
-<img src="figures_and_scripts/rigid_output.png" alt="Force and Moment Output" style="width:800px; background-color: #ffffff;"/>
+#### Postprocessing Procedure
+
+Pressure and viscous force vectors are read at all blade surface faces. The integrated values, as well as moments around the hub, are recorded every four timesteps. 
+
+Yaw angle: $\phi=30.0 \degree\\$
+Tilt angle: $\psi=5.0 \degree\\$
+Generator Efficiency: $E_g = 0.944\\$
+Rotor Speed: $\Omega = 12.1 rpm$
+
+$ f_x = f_{p_x}+f_{v_x}\\$
+$ f_y = f_{p_y}+f_{v_y}\\$
+$ f_z = f_{p_z}+f_{v_z}\\$
+
+$ f_{\phi x} = f_x*cos(-\phi)-f_y*sin(-\phi)\\$
+$ f_{\phi y} = f_x*sin(-\phi)+f_y*cos(-\phi)\\$
+$ f_{\phi z} = f_z\\$
+
+$ m_{\phi x} = m_x*cos(-\phi)-m_y*sin(-\phi)\\$
+$ m_{\phi y} = m_x*sin(-\phi)+m_y*cos(-\phi)\\$
+$ m_{\phi z} = m_z\\$
+
+$ Torque = m_{\phi x}*cos(-\psi)+m_{\phi z}*sin(-\psi)\\$
+$ Thrust =f_{\phi x}*cos(-\psi)+f_{\phi z}*sin(-\psi)\\$
+$ Power = Torque*\Omega*E_g/1000\\$
+
+![rigidoutput](figures_and_scripts/rigid_output.png)
+
+| variable | value |
+| -------- | ------- |
+|Generator Power|4756.45 kW|
+|Rotor Torque|3976.46 kN-m|
+|Rotor Thrust|650.25 kN|
+|Rotor Speed|12.1 rpm|
 
 
 ## Simulation Timings
 
-This benchmark was run on Sandia National Laboratories' machine "Flight" which is comprised of 112 Intel(R) Xeon(R) Platinum 8480+ CPU cores per node. 
+This benchmark was run at Sandia National Laboratories on a machine comprised of 112 Intel(R) Xeon(R) Platinum 8480+ CPU cores per node. 
 
 ### Processor Decomposition: 
 - Nalu-Wind: 672 ranks
-- AMR-Wind: 6496 ranks
+- AMR-Wind: 4928 ranks
 
 ### Timings:
 - Mean wall-clock time per timestep for entire simulation: 7.35s 
 - Mean wall-clock time per timestep per cell 4.135e-08s
 
+## Quick Exawind Simulation Guide
+### Step 0: Run the ABL precursor in AMR-Wind
+This simulation is driven by an ABL precursor, run in AMR-Wind, with recorded boundary planes, initial conditions, and average temperature. (See case linked above)
+### Step 1: Develop and decompose overset mesh
+It is necessary to develop a near-body mesh surrounding the blade surface in Nalu-Wind using external tools. This
+mesh should be in the Exodus II format and decomposed to approximately 20K cells per core available to Nalu-Wind. (Current mesh provided via dvc/github assets) 
+### Step 2: Run the OpenFAST precursor for 5760 OpenFAST timesteps (approximately 1 rotation)
+``srun -n 1 openfastcpp iyaml``
+### Step 3: Run the full Exawind suite
+``srun -N 50 -n 5600 exawind --nwind 672 --awind 4928 nrel5mw.yaml &> log``
+### Step 4: Setup and Run the full Exawind suite restart(s)
+``srun -N 50 -n 5600 exawind --nwind 672 --awind 4928 nrel5mw_r1.yaml &> log_1``
+
+
 ## References
 
 [1]: Jonkman, J. et al, "Definition of a 5-MW Reference Wind Turbine for Offshore System Development" NREL/TP-500-38060, 2009. https://www.nrel.gov/docs/fy09osti/38060.pdf

exawind/NREL_5MW_Turbine/rigid/deploy.yaml

@@ -1,5 +1,5 @@
 # Location where files should exist
-destdir: ../../.website_src/exawind/nrel5mw/rigid
+destdir: ../../../.website_src/exawind/nrel5mw/rigid
 # Define the valid file types to copy over (optional)
 validtypes: ['.html', '.md', '.ipynb', '.png', '.jpg', '.rst', ] 
 

exawind/NREL_5MW_Turbine/rigid/figures_and_scripts/get_timings.py

@@ -117,13 +117,18 @@ def processline(inputline):
 def main():
 
     casedir = '/pscratch/ndeveld/hfm-2025-q1'
-    casename = 'nrel5mw_rigid_abl_noopenfast2'
+    casename = 'rigid_abl_splitmsh_nate2'
     casepath = os.path.join(casedir,casename)
 
+    plotdir = os.path.join(casepath,'plots')
+
     exlogfile = casepath+"/log"
     amrlogfile = casepath+"/nrel5mw_amr.log"
     nalulogfile = casepath+"/nrel5mw_nalu.log"
+    prefix = 'rign2'
 
+    #cmd="grep '^Exawind::Total' "+exlogfile+" | awk '{print $3}' > "+casepath+"/exatimestep.dat"
+    #result = sp.run(cmd, shell=True, capture_output=True, text=True)
 
     cmd="grep '^Exawind::Total' "+exlogfile+" | awk '{print $3}'"
     result = sp.run(cmd, shell=True, capture_output=True, text=True)
@@ -137,9 +142,30 @@ def main():
     result = sp.run(cmd, shell=True, capture_output=True, text=True)
     amr_timesteps = [float(x) for x in result.stdout.replace('\n',' ').split()]
 
+    cmd="grep '^  MAC_projection ' "+amrlogfile+" | awk '{print $2}'"
+    result = sp.run(cmd, shell=True, capture_output=True, text=True)
+    amr_mac= [float(x) for x in result.stdout.replace('\n',' ').split()]
+
+    cmd="grep '^  Nodal_projection ' "+amrlogfile+" | awk '{print $2}'"
+    result = sp.run(cmd, shell=True, capture_output=True, text=True)
+    amr_nodal= [float(x) for x in result.stdout.replace('\n',' ').split()]
+
+    cmd="grep '^        MomentumEQS ' "+nalulogfile+" | awk '{print $2}'"
+    result = sp.run(cmd, shell=True, capture_output=True, text=True)
+    nalu_mom= [float(x) for x in result.stdout.replace('\n',' ').split()]
+
+    cmd="grep '^        ContinuityEQS ' "+nalulogfile+" | awk '{print $2}'"
+    result = sp.run(cmd, shell=True, capture_output=True, text=True)
+    nalu_cont= [float(x) for x in result.stdout.replace('\n',' ').split()]
+
 
     #tsdata = pd.read_csv(casepath+'/avgtimestep.dat',header=None)
     ts = range(len(total_timesteps))
+    nts = range(len(nalu_timesteps))
+    ats = range(len(amr_timesteps))
+    amrts = range(len(amr_mac)) 
+    naluts = range(len(nalu_mom))
+    conts =  range(len(nalu_cont))
 
     print('Mean timestep',np.mean(total_timesteps))
 
@@ -168,14 +194,23 @@ def main():
 
     print('Mean Timestep per cell',np.mean(total_timesteps)/(amr_cells+nalu_cells))
 
+    cs = 5  # Plot the 5th occurence of continuity iters
+    ms = 4  # Plot the 4th occurence of momentum iters
+
+    nmd_filt = np.array(nalu_mom)[0::ms].copy()
+    ncd_filt = np.array(nalu_cont)[0::cs].copy()
+    nmts_filt = np.arange(len(nmd_filt))
+    ncts_filt = np.arange(len(ncd_filt))
+
+    # Plot main exawind
     plt.rcParams.update({'font.size': 18})
 
     fig, ax = plt.subplots(1,3,figsize=(13,4))
     ax[0].scatter(ts,total_timesteps,s=0.3)
     ax[0].set_title('Exawind')
-    ax[1].scatter(ts,nalu_timesteps,s=0.3)
+    ax[1].scatter(nts,nalu_timesteps,s=0.3)
     ax[1].set_title('Nalu-Wind')
-    ax[2].scatter(ts,amr_timesteps,s=0.3)
+    ax[2].scatter(ats,amr_timesteps,s=0.3)
     ax[2].set_title('AMR-Wind')
 
     for i in range(3):
@@ -186,6 +221,30 @@ def main():
     fig.tight_layout()
     fig.savefig('timepertimestep.png')
 
+    # Plot AMR
+    plt.rcParams.update({'font.size': 18})
+
+    in_min = min(len(amrts),len(amr_nodal))-1
+
+    fig, ax = plt.subplots(1,2,figsize=(13,4))
+    ax[0].scatter(amrts[0:in_min],amr_mac[0:in_min],s=0.3,label="AMR MAC Projection")
+    ax[0].scatter(amrts[0:in_min],amr_nodal[0:in_min],s=0.3,label="AMR Nodal Projection")
+    ax[0].set_title('AMR-Wind')
+
+    ax[1].scatter(nmts_filt,nmd_filt,s=0.3,label="Nalu Momentum")
+    ax[1].scatter(ncts_filt,ncd_filt,s=0.3,label="Nalu Continuity")
+    ax[1].set_title('Nalu-Wind')
+
+    for i in range(2):
+        ax[i].set_ylabel('N')
+        ax[i].set_xlabel('Timestep')
+        ax[i].legend()
+        ax[i].set_ylim([0,50])
+
+    fig.tight_layout()
+    fig.savefig('iters.png')
+
+
 
 
 

exawind/NREL_5MW_Turbine/rigid/figures_and_scripts/timeseries_process.py

@@ -58,45 +58,74 @@ def main():
 
     args = parser.parse_args()
 
-    case_list = ['nrel5mw_rigid_abl_noopenfast2']
-    force_file_names = ['forces01.dat']
+    case_list = ['rigid_abl_splitmsh_nate2']
+    force_file_names = ['forcesBlades.dat']
     case_lab = ['NREL 5MW Rigid']
 
-    omega = 1.25663706
+    omega = 1.26710903694788
     diffn = 60000
     rotaxis = [0.862729916,0.498097349,-0.087155742]
 
     matplotlib.rcParams['font.size'] = 16
 
-    fig = plt.figure(constrained_layout=True,figsize=(13,4))
+    fig = plt.figure(constrained_layout=True,figsize=(12,4))
     subfigs = fig.subfigures(nrows=1, ncols=1)
-    ax = subfigs.subplots(nrows=1, ncols=2)
+    ax = subfigs.subplots(nrows=1, ncols=3)
 
     for i,c in enumerate(case_list):
 
         print('Processing: ',c)
 
-        fullpath = os.path.join(args.directory,c,force_file_names[i])
-        print(fullpath)
+        this_data = []
 
-        this_data = pd.read_csv(fullpath,sep='\s+',skipinitialspace=True)
+        for j in range(len(force_file_names)):
+            fullpath = os.path.join(args.directory,c,force_file_names[j])
+            print(fullpath)
 
-        this_data['Thrust'] = (this_data['Fpx']*rotaxis[0] + this_data['Fpx']*rotaxis[1])/1000.0
+            this_data.append(pd.read_csv(fullpath,sep='\s+',skipinitialspace=True))
 
-        ax[0].plot(this_data['Time'], this_data['Thrust'], label=case_lab[i])
-        #plt.axhline(y = ofpower, color = 'k', linestyle = ':')
+        yawangle = 30.0*np.pi/180.0
+        tiltangle = 5.0*np.pi/180.0
+
+
+        all_data = pd.concat(this_data, ignore_index=True)
+
+
+        all_data['fx'] = all_data['Fpx']+all_data['Fvx']
+        all_data['fy'] = all_data['Fpy']+all_data['Fvy']
+        all_data['fz'] = all_data['Fpz']+all_data['Fvz']
+
+        all_data['frot30x'] = all_data['fx']*np.cos(-yawangle)-all_data['fy']*np.sin(-yawangle)
+        all_data['frot30y'] = all_data['fx']*np.sin(-yawangle)+all_data['fy']*np.cos(-yawangle)
+        all_data['frot30z'] = all_data['fz']       
+
+        all_data['mrot30x'] = all_data['Mtx']*np.cos(-yawangle)-all_data['Mty']*np.sin(-yawangle)
+        all_data['mrot30y'] = all_data['Mtx']*np.sin(-yawangle)+all_data['Mty']*np.cos(-yawangle)
+        all_data['mrot30z'] = all_data['Mtz']
+
+        all_data['Torque'] = (all_data['mrot30x']*np.cos(-tiltangle)+all_data['mrot30z']*np.sin(-tiltangle))/1000
+        all_data['Thrust'] = (all_data['frot30x']*np.cos(-tiltangle)+all_data['frot30z']*np.sin(-tiltangle))/1000
+
+        all_data['Power'] = all_data['Torque']*omega*0.944
+
+        ax[0].plot(all_data['Time'], all_data['Thrust'], label=case_lab[i])
         ax[0].set_xlabel("Time [s]")
         ax[0].set_ylabel("Thrust [kN]")
-        ax[0].set_ylim([150,900])
-        #ax.set_xlim([0,30])
-        #ax[0].legend()
 
-        ax[1].plot(this_data['Time'], this_data['Mtx']/1000.0, label="x-dir Moment")
-        ax[1].plot(this_data['Time'], this_data['Mty']/1000.0, label="y-dir Moment")
-        ax[1].plot(this_data['Time'], this_data['Mtz']/1000.0, label="z-dir Moment")
+        ax[1].plot(all_data['Time'], all_data['Torque'], label=case_lab[i])
         ax[1].set_xlabel("Time [s]")
-        ax[1].set_ylabel("Moment [kN-m]")
-        ax[1].legend()
+        ax[1].set_ylabel("Torque [kN-m]")
+
+        ax[2].plot(all_data['Time'], all_data['Power'], label=case_lab[i])
+        ax[2].set_xlabel("Time [s]")
+        ax[2].set_ylabel("Power [kW]")
+
+    mean_data = all_data[all_data['Time']>30.0].mean()
+
+    print('Mean Power: ',mean_data.Power,'kW')
+    print('Mean Thrust: ',mean_data.Thrust,'kN')
+    print('Mean Torque: ',mean_data.Torque,'kN-m')
+
 
     plt.savefig('rigid_output.png')
     plt.close()