Example Python Code

This jupyter notebook is an example for generating a valid netCDF file for uploading to perform the DLC response analysis.

Code similar to this example is intended to be used after tank testing or a dynamic simulation of your WEC device is performed, using the surface elevation time series for your realized sea states downloaded through the DLC Generator tool.

# imports - NOTE: xarray issues a warning for a backend engine that isn't needed for this code
import numpy as np
import xarray as xr
import pandas as pd
from scipy.integrate import solve_ivp

Data Structure

The data structure used is a netCDF file, where the coordinates are time and case, and the variables are the response value arrays that are indexed by case and time.

Data structure requirements

1. Sample rate (time variable) needs to be consistent across all cases and response variables.
2. Data variables (device response arrays) need to exist across all cases. IE Case 1 having 4 response variables and Case 2 having 5 is invalid. Both would need 5 response variables.
3. Time array is in seconds (simulation/ tank time), NOT a datetime or string.
4. No NaN's/ nulls can be present in response data variables.

Data that would come from tank testing or dynamic simulation

Suppose you've calculated contour statistics for a data source and computed 2 sea states from the Hm50 contour line. Having downloaded the surface elevation traces, we can perform a simulation to generate some WEC data.

# the Hm0, Tp values for the  downloaded cases
sea_states = [
    {"Hm0": 11.92, "Tp": 10.20},
    {"Hm0": 21.33, "Tp": 17.70},
]
# Helper function for file loading
def decimal_to_p(f):
    return "{:.2f}".format(f).replace(".", "p")

# Combine all the wave cases into a single DataFrame
elevations = []
for sea_state in sea_states:
    Hm0 = sea_state["Hm0"]
    Tp = sea_state["Tp"]
    filename = "Hm0-{}_Period-{}.csv".format(decimal_to_p(Hm0), decimal_to_p(Tp))
    case = "RM3 WEC Response Hm0: {} Tp: {}".format(Hm0, Tp)
    eta = pd.read_csv(filename, index_col=0, names=["time", case], header=0)

    elevations.append(eta)

elevation=pd.concat(elevations,axis=1)
    

Perform simple simulation

Here we perform a very simple simulation for a WEC that experiences a force proportional the surface elevation.

Set up ODE

m = 500 # kg
d = 1000 # Ns/m
k = 600 # N/m
def ode(t, x, case):
    eta = elevation[case]
    # External force
    f_external = 400 * np.interp(t, eta.index.values, eta.values)
    pos = x[0]
    vel = x[1]

    acc = 1/m * (f_external - d * vel - k * pos)

    return [vel, acc]

Set up quantities we want to output

These have the same shape as the original surface elevation data

pto_displacement = np.zeros(elevation.shape)
pto_force = np.zeros(elevation.shape)

Perform the simulation for each case

for i, data in enumerate(elevations):
    time = data.index.values
    sol = solve_ivp(ode, t_span=[time[0], time[-1]], y0 = [0.0, 0.0], t_eval=time, args=(data))

    pos = sol.y[0, :]
    vel = sol.y[1, :]

    force = d*vel + k*pos
    
    pto_displacement[:, i] = pos
    pto_force[:, i] = force

Option 1 - Create a combined xarray Dataset

data = xr.Dataset(
    # Create the data variables, indicating they're indexed by case then time.
    # Additionally a dictionary of attributes can be provided: 'units' and 'name'
    # will be used to display the data in the DLC UI.
    data_vars=dict(
        pto_displacement=(["case", "time"], pto_displacement.T, {"units":"m", "name":"PTO Displacement"}),
        pto_force=(["case", "time"], pto_force.T, {"units":"N", "name":"PTO Force"}),
    ),
    # The coordinates must be supplied in this form.
    coords=dict(
        time=time,
        case=(["case"], elevation.columns)
    )
)

# show data structure
print(data)

<xarray.Dataset>
Dimensions:           (case: 2, time: 216000)
Coordinates:
  * time              (time) float64 0.0 0.05 0.1 ... 1.08e+04 1.08e+04 1.08e+04
  * case              (case) object 'DLC 6.2 Hm0: 11.92 Tp: 10.2' 'DLC 6.2 Hm...
Data variables:
    pto_displacement  (case, time) float64 0.0 -0.0002085 ... -1.745 -1.757
    pto_force         (case, time) float64 0.0 -8.202 ... -1.304e+03 -1.286e+03

# Save data for upload
data.to_netcdf("response-upload.nc")

Option 2 - Create a CSV for each case

for i, data in enumerate(elevations):
    name = data.columns[0].replace(" ", "_").replace(":", "")
    d = np.array([pto_displacement[:,i], pto_force[:,i]]).T
    df = pd.DataFrame(data=d, index=data.index.values, columns=('pto_displacement', "pto_force"))
    
    df.to_csv(name+".csv", index_label="time")