Wake center curl model FAST.Farm

[MYMahfouz]

I am trying to extract the wake center from the vtk file of FAST.Farm output. To find the wake center of the

I am getting a wake center lower than the hub height of the turbine. The hub height is at 135m and the height of the wake center of the figure below is at 112m.

I believe I followed this example but maybe I overlooked something.

I have also attached the code i am using

from samwich.waketrackers import track
from samwich.dataloaders import PlanarData
import os
import numpy as np
import pyFAST.input_output.vtk_file as vtk_file

# Temporary inputs - Used only to make the specification of file names easier
# Main inputs
numPlans = 1 # Number of planes to be processed. Assumed to be equally spaced
dt = 3.6       # Time step between vtk files
folderPath = r'G:\FAST_Farm_Template\farm_9T_Mann_baseline_avg_V_curl_v2\Cond00_v10.0_PL0.2_TI6\Case00_wdir0.0_fixed_kcurl6\Seed_0' # Path to the folder with the FAST.Farm simulation
rootNameVTK = 'FFarm_mod.Low.DisYZ1.'     # Root name of the test case
outFolder = os.path.join(folderPath, 'wakescenter') # Output folder. Output files will have the same name as the vtks, but with .csv extensions.
rotorDiam = 240 # Rotor diameter
rootPath = os.path.join(folderPath, 'vtk_ff') # Subfolder with vtk's generated by FAST.Farm



# Create the output folder if it does not exist
if not os.path.exists(outFolder):
    os.makedirs(outFolder)





# List all vtk files

planes_xPos = np.zeros(1)+2519.22

listNames = []
listXPos  = []
for i in range(1):
    # Some versions of FAST.Farm write VTk names as 01 and others 1
    nameVTK = f'{rootNameVTK}{i+1000:05d}'
    nameVTK_wo0 = f'{rootNameVTK}{i+1000}'

    # List all files with nameVTK in the name
    thisListOfNames = [f for f in os.listdir(rootPath) if os.path.isfile(os.path.join(rootPath, f)) and nameVTK in f or nameVTK_wo0 in f]
    listNames.extend(thisListOfNames)
    listXPos.extend([planes_xPos[i]]*len(thisListOfNames))    

list_vtk = []
for vtkName in listNames[0:1]:
    filevtk=vtk_file.VTKFile(rootPath + '/' + vtkName)
    filevtk.yp_grid = filevtk.yp_grid[3:14]
    filevtk.zp_grid = filevtk.zp_grid[0:9]
    filevtk.point_data_grid['Velocity']=filevtk.point_data_grid['Velocity'][:,3:14,0:9,:]
    list_vtk.append(filevtk)


# Check how many files
print(len(list_vtk))




# Sort listNames alphabetically. That is enough to make sure that things are sorted by time
listNames = sorted(listNames)

for i in range(numPlans):
    # Get indices of files for this plan
    ind4plan = [j for j, nm in enumerate(listNames) if f'{rootNameVTK}{i+1000:05d}' in nm or f'{rootNameVTK}{i+1000+1}' in nm]
    # Write to a csv file
    outFlPath = os.path.join(outFolder, f'{rootNameVTK}{i+1000:05d}.csv')    

    with open(outFlPath, 'w') as f:
        f.write('t,x,y,z')        

        dict4wake = {'y':[], 'z':[], 'u':[], 'v':[], 'w':[]}   # Used by SAMWICH
        for iTime, iP in enumerate(ind4plan):            
            time = dt*int(listNames[iP].split('.')[-2]) # TIme
            vtk = list_vtk[iP]            

            # 2D arrays with shape len(vtk.yp_grid): rows, len(vtk.zp_grid): columns
            dict4wake['u'] = vtk.point_data_grid['Velocity'][0, :, :, 0]
            dict4wake['v'] = vtk.point_data_grid['Velocity'][0, :, :, 1]
            dict4wake['w'] = vtk.point_data_grid['Velocity'][0, :, :, 2]
            xc = listXPos[i] # X position of the wake center (constant and equal to plane position)

            # Need 2D arrays with same shape as u, v, w
            dict4wake['z'], dict4wake['y'] = np.meshgrid(vtk.zp_grid, vtk.yp_grid)
            inp=dict4wake.copy()
            wakedata = PlanarData(inp)
            wake = track(wakedata.sliceI(), method='ConstantArea', verbose=True)
            
            wake.clear_plot() 
            wake.plot_contour(vmin=1,vmax=13)
            
            # Compute wake centers
            refArea = np.pi*240**2/4
            
            wake.remove_shear(alpha=0.2, zref=135, Uref=10)
            wake.clear_plot()  # force replot for python notebook
            wake.plot_contour(vmin=-5,vmax=5)             
            yc, zc = wake.find_centers(refArea,min_contour_points=5,verbosity=2) # y and z coordinates of the wake center

            
            # Write to file
            f.write(f'{time:.3f},{xc:.3f},{yc[0]:.3f},{zc[0]:.3f}')
            f.flush()
            wake.clear_plot()  
            wake.plot_contour(vmin=-5,vmax=5,outline=True)

I would really appreciate your help.

Best regards,

[ewquon]

Hello, it looks like you're using the "ConstantArea" wake tracker for a 240-m diameter rotor. Your velocity deficit region does appear to have a comparable size, the background wind appears to have been removed, and overall your usage appears reasonable. Could you share one of your wake data snapshots?

Have you tried the Gaussian2D tracker? That's probably the fastest and most robust.

[salicto]

FAST.Farm.Low.DisYZ01.0546.csv

Hello everyone,
I’m encountering difficulties using the various functions in the samwich library, specifically with contour wake tracking. I started with the same approach as MYMahfouz, converting FAST.Farm .vtk files to .csv format, matching the structure of the example file provided.
I’ve tried using the _find_contour_center function, but none of the methods seem to work as expected. Here’s the code I’m using:

`from samwich.dataloaders import RawData
from samwich.waketrackers import ContourWakeTracker
import samwich.waketrackers as waketrackers
from samwich.waketrackers import Plotter, track
import inspect
import os

for name, obj in inspect.getmembers(waketrackers):
if inspect.isfunction(obj):
print(f"Fonction : {name}")
elif inspect.isclass(obj):
print(f"Classe : {name}")

D=284
exampleData=os.path.join(outFolder,'FAST.Farm.Low.DisYZ01.0546.csv')
with open(exampleData,'r') as f:
for _ in range(5):
print(f.readline().strip())

data=np.loadtxt(exampleData,skiprows=1,delimiter= ',')
NX = len(np.unique(data[:, 0]))
NY = len(np.unique(data[:, 1]))
NZ = len(np.unique(data[:, 2]))
row,col=data.shape

x = data[:, 0]
y = data[:, 1]
z = data[:, 2]
u = data[:, 3]
v = data[:, 4]
w = data[:, 5]

y_unique = np.unique(y)
z_unique = np.unique(z)
Y, Z = np.meshgrid(y_unique, z_unique, indexing='ij')

U = np.zeros((len(y_unique), len(z_unique)))
V = np.zeros((len(y_unique), len(z_unique)))
W = np.zeros((len(y_unique), len(z_unique)))

for i, y_val in enumerate(y_unique):
for j, z_val in enumerate(z_unique):
mask = (y == y_val) & (z == z_val)
if np.any(mask):
U[i, j] = u[mask][0]
V[i, j] = v[mask][0]
W[i, j] = w[mask][0]

datadict = {
'x': np.full((1, len(y_unique), len(z_unique)), x[0]), # x est constant plan y0z
'y': Y.reshape((1, len(y_unique), len(z_unique))),
'z': Z.reshape((1, len(y_unique), len(z_unique))),
'u': U.reshape((1, len(y_unique), len(z_unique))),
'v': V.reshape((1, len(y_unique), len(z_unique))),
'w': W.reshape((1, len(y_unique), len(z_unique)))
}
c1,c2 =(np.min(data[:,3]),np.max(data[:,3]))
planar_data = PlanarData(datadict, center_x=False, center_y=True)

'''tracker = '1-D Gaussian'
wake = track(planar_data.sliceI(), method='Gaussian2D', verbose=True)
wake.remove_shear(alpha=0.13, zref=172, Uref=8)
wake.plot_contour(vmin=c1,vmax=c2)
plotter = Plotter(y=Y, z=Z, u=U, figsize=(10, 6), dpi=100, cmap='viridis', vmin=c1, vmax=c2)'''

wake = ContourWakeTracker(planar_data.sliceI())

'''wakecenter=track(wake,method='Gaussian2D',verbose=True)
#wake.apply_filter(filtertype='gaussian', N=100, size=1, minsize=2)
wake.clear_plot()
wake.plot_contour(vmin=c1,vmax=c2)
wake.remove_shear(alpha=0.13, zref=172, Uref=8)
wake.clear_plot()
wake.plot_contour(vmin=0,vmax=6)'''

refArea=np.pi*D**2/4
itime = 0
target_value =refArea
weighted_center = True
contour_closure = False
min_contour_points = 10
Ntest0 = 20
Ntest = 4
umax = np.max(planar_data.data)
tol = 0.01
func = None
fields = ('u_tot',)
vdcheck = True
debug = True

wake._find_contour_center(itime=itime,
target_value=target_value,
weighted_center=weighted_center,
contour_closure=contour_closure,
min_contour_points=min_contour_points,
Ntest0=Ntest0,
Ntest=Ntest,
umax=umax,
tol=tol,
func=func,
fields=fields,
vdcheck=vdcheck,
debug=debug)

print(f"Centre du sillage : yc = {yc}, zc = {zc}")`

I try a another method from the classe contours but this time the algorihm find four different contours and non of them is correct.

`exampleData=os.path.join(outFolder,'FAST.Farm.Low.DisYZ01.0546.csv')
data=np.loadtxt(exampleData,skiprows=1,delimiter= ',')

x = data[:, 0]
y = data[:, 1]
z = data[:, 2]
u = data[:, 3]

umin=np.min(u)
umax=np.max(u)
print(umax,umin)

x_target = 1284.0
mask = x == x_target
y_filtered = y[mask]
z_filtered = z[mask]
u_filtered = u[mask]

y_unique = np.unique(y_filtered)
z_unique = np.unique(z_filtered)
Y, Z = np.meshgrid(y_unique, z_unique, indexing='ij')

U = np.zeros((len(y_unique), len(z_unique)))
for i, y_val in enumerate(y_unique):
for j, z_val in enumerate(z_unique):
idx = np.where((y_filtered == y_val) & (z_filtered == z_val))[0]
if len(idx) > 0:
U[i, j] = u_filtered[idx[0]]

contours = Contours(Y, Z, u=U-np.mean(U), umin=umin, umax=umax)
Clevel = 0.4* np.max(U)

print(f"Nombre de contours fermés trouvés : {len(closed_paths)}")

centers = []
for contour in closed_paths:
xc, yc = contours.calc_weighted_center(contour)
centers.append((xc, yc))

print(centers)

plt.figure(figsize=(10, 6))
plt.contourf(Y,Z,U-np.mean(U), levels=20, cmap='viridis')
plt.colorbar(label='Déficit de vitesse')
for contour, (xc, yc) in zip(closed_paths, centers):
x_contour, y_contour = contours.to_coords(contour, closed=True)
plt.plot(x_contour, y_contour, 'r-', linewidth=1)
plt.scatter(xc, yc, color='red', marker='x', s=50)

plt.title('Contours du sillage et leurs centres')
plt.xlabel('y')
plt.ylabel('z')
plt.grid(True)
plt.show()`

Any guidance or suggestions on how to resolve these issues would be greatly appreciated!
Kind regards,