EnvironmentSTUDIO¶

EnvironmentSTUDIO flow (es-flow for short) is a library for analysing and modeling atmospheric and marine boundary layers.

While es-flow can be used for any turbulent boundary layer analysis, its main focus is for wind and tidal site characterisation - generating the ‘best fit’ of analytical models to measured velocity data.

A key strength of es-flow is the adem. This extremely robust method allows users to:

determine detailed turbulence information from instruments like LiDAR

characterise turbulence and shear beyond tip height of even the biggest offshore wind turbines

characterise coherent structure in turbulence, crucially important for fatigue loading in wind turbines.

Aims¶

The es-flow library provides: [1]

Parameterised profiles
- Mean Velocity (using power law, logarithmic law, MOST or Lewkowicz relations)
- Mean Veer (Using Monin-Obukhov approach)
- Reynolds Stress u'w' (using Lewkowicz relations)
- Reynolds Stress u'u', u'v', u'w', v'v', v'w', w'w' (using the adem)
- Spectra Sij (using Kaimal, von Karman)
- Spectra Sij (using adem)
- Integral turbulent intensity and lengthscale I, l
Best fit parameter sets
- To describe the above profiles analytically (given measured velocity data from an instrument).

In future, generation of artificial flow fields for simulation purposes might be considered. This would overlap with - or replace - utilities like TurbSim to produce, for example:

.wnd fields input to BEM or FVM models like Bladed, FAST and TurbineGRID
Inlet boundary conditions for DES or LES codes

Uses¶

At Octue, es-flow is used to:

Provide a basis for developing and validating new processes, like the adem, for characterising Atmospheric Boundary Layers.

Process LiDAR and ADCP datasets from raw instrument files.

Apply windowed analyses to time series data, for flow characterisation.

Generate load cases for FAST, Bladed and our own TurbineGRID aerodynamic wind turbine analysis tools.

We’d like to hear about your use case. Please get in touch!

We use the GitHub Issue Tracker to manage bug reports and feature requests.

Analytical Models¶

Velocity Relations¶

es-flow uses a range of different velocity relations, so you can choose which is preferable. Fitting different profiles returns a metric of accuracy; but this isn’t the be-all and end-all. Always choose which formulation best represents the physics in play!

Power Law¶

Power law

Logarithmic¶

Log law

MOST¶

Most law

Lewkowicz¶

Lewkowicz 1982, modified by P&M. Include the Reynolds Stress profile

Veer Relations¶

Monin obukhov¶

Spectra¶

Kaimal¶

Attached-Detached Eddy Method¶

Describe ADEM here.

Examples¶

Add example usages here.

File Formats¶

Add file format details here.

Installation¶

Third party library installation¶

Intel MKL¶

Attention

If you don’t wish to use Intel MKL, or need to build for a non-intel architecture, please contact Octue.

Download the Intel MKL library packages. Click on the icon and follow installation instructions. You’ll need the administrator password.

The tools are installed in /opt/intel/, the include directory is /opt/intel/include.

matio¶

Whatever you do, don’t try to fork and build from source - the autoconf is complex and not suitable for OSX. Luckily there’s a brew formula:

Please contact Octue for Linux installation instructions.

Please contact Octue for Windows installation instructions.

ceres-solver, eigen and glog¶

Google’s ceres-solver also depends on glog and eigen, so we get three for the price of one:

Please contact Octue for Linux installation instructions.

Please contact Octue for Windows installation instructions.

Third party build requirements¶

Attention

These dependencies are only required if you’re building es-flow from source.

cxxopts¶

To build es-flow, cxxopts must be placed alongside es-flow. From the es-flow root directory:

Please contact Octue for Linux installation instructions.

Please contact Octue for Windows installation instructions.

NumericalIntegration¶

To build es-flow, NumericalIntegration must be placed alongside es-flow. From the es-flow root directory:

Please contact Octue for Linux installation instructions.

Please contact Octue for Windows installation instructions.

License¶

Octue offers free academic licensing for es-flow as well as a commercial subscription. Please contact tom@octue.com for details.

Third Party Libraries¶

es-flow includes or is linked against the following third party libraries:

matio¶

A library for reading and writing MAT files.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY CHRISTOPHER C. HULBERT ‘‘AS IS’’ AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CHRISTOPHER C. HULBERT OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

ceres-solver¶

Ceres Solver is licensed under the New BSD license, whose terms are as follows.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
Neither the name of Google Inc., nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

This software is provided by the copyright holders and contributors “AS IS” and any express or implied warranties, including, but not limited to, the implied warranties of merchantability and fitness for a particular purpose are disclaimed. In no event shall Google Inc. be liable for any direct, indirect, incidental, special, exemplary, or consequential damages (including, but not limited to, procurement of substitute goods or services; loss of use, data, or profits; or business interruption) however caused and on any theory of liability, whether in contract, strict liability, or tort (including negligence or otherwise) arising in any way out of the use of this software, even if advised of the possibility of such damage.

Eigen¶

Eigen is licensed under the Mozilla Public License v2.

NumericalIntegration¶

An extension library for Eigen that allows numerical integration, used under the Mozilla Public License v2.

cumtrapz¶

cumtrapz.h header file utility for numerical integration is adapted from The Biomechanical Toolkit.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

Neither the name(s) of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

cxxopts¶

An argument parser for C++11 under the MIT license.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Version History¶

Origins¶

EnvironmentSTUDIO flow began as an internal tool at Octue, for characterising turbulence at tidal power sites. It was originated in MATLAB (then called the ‘Flow Characterisation Suite’), and grew incrementally.

Unfortunately, MATLAB is extremely unwieldy if deployed in production, requiring either (complicated and badly limited) cross compilation, expensive cloud sever licenses or build and deployment with MATLAB’s 3Gb+ runtime library (prohibitively big for practical use in cloud services).

As the move was made toward object orientation, and usage increased toward requiring a production version that could be deployed without requiring expensive MATLAB licenses, EnvironmentSTUDIO flow project was started in C++.

Gradually, capability is being ported from MATLAB (the original repo having been subsumed within this project and now being progressively deprecated).

0.1.0¶

This version bump was funded by AURA via the University of Hull. The objective of the work was to validate capabilities using measured data, and the library was refactored to allow this work to happen in collaboration with Univ. Hull and the Offshore Renewable Energy Catapult.

New Features¶

Scope of library reduced to preclude specific instrument data readers, and focus on environmental characterisation.
Library API consolidated.
ADEM functionality ported from legacy MATLAB and tested to work with C++.
Added plotting routines (using cpplot).
Documentation and build systems implemented and settled.
Google Test used to create test harness for the majority of the library.

Backward Incompatible API Changes¶

Entire API altered to reflect change in scope; will be much more stable going forward.

Bug Fixes & Minor Changes¶

n/a

0.0.1¶

New Features¶

Development version. Highly unstable.

Backward Incompatible API Changes¶

Everything. API changes daily.

Bug Fixes & Minor Changes¶

Minor???!

Bibliography¶

[Agarwal]

S. Agarwal, N. Snavely, S. M. Seitz and R. Szeliski, Bundle Adjustment in the Large, Proceedings of the European Conference on Computer Vision, pp. 29–42, 2010.

[Bjorck]

A. Bjorck, Numerical Methods for Least Squares Problems, SIAM, 1996

[Brown]

D. C. Brown, A solution to the general problem of multiple station analytical stereo triangulation, Technical Report 43, Patrick Airforce Base, Florida, 1958.

[ByrdNocedal]

R. H. Byrd, J. Nocedal, R. B. Schanbel, Representations of Quasi-Newton Matrices and their use in Limited Memory Methods, Mathematical Programming 63(4):129–-156, 1994.

[ByrdSchnabel]

R.H. Byrd, R.B. Schnabel, and G.A. Shultz, Approximate solution of the trust region problem by minimization over two dimensional subspaces, Mathematical programming, 40(1):247–263, 1988.

[Chen]

Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam, Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate, TOMS, 35(3), 2008.

[Conn]

A.R. Conn, N.I.M. Gould, and P.L. Toint, Trust region methods, Society for Industrial Mathematics, 2000.

[GolubPereyra]

G.H. Golub and V. Pereyra, The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate, SIAM Journal on numerical analysis, 10(2):413–432, 1973.

[HartleyZisserman]

R.I. Hartley & A. Zisserman, Multiview Geometry in Computer Vision, Cambridge University Press, 2004.

[KanataniMorris]

K. Kanatani and D. D. Morris, Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy, IEEE Transactions on Information Theory 47(5):2017-2028, 2001.

[Keys]

R. G. Keys, Cubic convolution interpolation for digital image processing, IEEE Trans. on Acoustics, Speech, and Signal Processing, 29(6), 1981.

[KushalAgarwal]

A. Kushal and S. Agarwal, Visibility based preconditioning for bundle adjustment, In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2012.

[Kanzow]

C. Kanzow, N. Yamashita and M. Fukushima, Levenberg–Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints, Journal of Computational and Applied Mathematics, 177(2):375–397, 2005.

[Levenberg]

K. Levenberg, A method for the solution of certain nonlinear problems in least squares, Quart. Appl. Math, 2(2):164–168, 1944.

[LiSaad]

Na Li and Y. Saad, MIQR: A multilevel incomplete qr preconditioner for large sparse least squares problems, SIAM Journal on Matrix Analysis and Applications, 28(2):524–550, 2007.

[Madsen]

K. Madsen, H.B. Nielsen, and O. Tingleff, Methods for nonlinear least squares problems, 2004.

[Mandel]

J. Mandel, On block diagonal and Schur complement preconditioning, Numer. Math., 58(1):79–93, 1990.

[Marquardt]

D.W. Marquardt, An algorithm for least squares estimation of nonlinear parameters, J. SIAM, 11(2):431–441, 1963.

[Mathew]

T.P.A. Mathew, Domain decomposition methods for the numerical solution of partial differential equations, Springer Verlag, 2008.

[NashSofer]

S.G. Nash and A. Sofer, Assessing a search direction within a truncated newton method, Operations Research Letters, 9(4):219–221, 1990.

[Nocedal]

J. Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Mathematics of Computation, 35(151): 773–782, 1980.

[NocedalWright]

J. Nocedal & S. Wright, Numerical Optimization, Springer, 2004.

[Oren]

S. S. Oren, Self-scaling Variable Metric (SSVM) Algorithms Part II: Implementation and Experiments, Management Science, 20(5), 863-874, 1974.

[Ridders]

C. J. F. Ridders, Accurate computation of F’(x) and F’(x) F”(x), Advances in Engineering Software 4(2), 75-76, 1978.

[RuheWedin]

A. Ruhe and P.Å. Wedin, Algorithms for separable nonlinear least squares problems, Siam Review, 22(3):318–337, 1980.

[Saad]

Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.

[Stigler]

S. M. Stigler, Gauss and the invention of least squares, The Annals of Statistics, 9(3):465-474, 1981.

[TenenbaumDirector]

J. Tenenbaum & B. Director, How Gauss Determined the Orbit of Ceres.

[TrefethenBau]

L.N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997.

[Triggs]

B. Triggs, P. F. Mclauchlan, R. I. Hartley & A. W. Fitzgibbon, Bundle Adjustment: A Modern Synthesis, Proceedings of the International Workshop on Vision Algorithms: Theory and Practice, pp. 298-372, 1999.

[Wiberg]

T. Wiberg, Computation of principal components when data are missing, In Proc. Second Symp. Computational Statistics, pages 229–236, 1976.

[WrightHolt]

S. J. Wright and J. N. Holt, An Inexact Levenberg Marquardt Method for Large Sparse Nonlinear Least Squares, Journal of the Australian Mathematical Society Series B, 26(4):387–403, 1985.

Library API¶

Class Hierarchy¶

- Namespace es
  - Struct LewkowiczSpeedResidual
  - Struct PowerLawSpeedResidual
  - Class AdemData
  - Class BasicLidar
  - Class Bins
  - Class EddySignature
  - Template Class Profile
  - Template Class Reader
  - Class VelocityProfile
  - Enum file_type_check_level
  - Enum timeseries_check_level
- Class spectra

File Hierarchy¶

- Directory source
  - Directory adem
    - File adem.h
    - File signature.h
  - Directory relations
    - File spectra.cpp
    - File spectra.h
    - File stress.h
    - File veer.h
    - File velocity.h
  - Directory utilities
    - File conv.h
    - File cumtrapz.h
    - File filter.h
    - File tensors.h
    - File trapz.h
  - File data_types.h
  - File definitions.h
  - File fit.h
  - File how_to.h
  - File main.cpp
  - File profile.cpp
  - File profile.h
  - File readers.h
  - File variable_readers.h

Full API¶

Namespaces¶

Namespace es¶

Contents

Classes
Enums
Functions

Classes¶

Enums¶

Functions¶

Namespace utilities¶

Contents

Functions

Functions¶

Classes and Structs¶

Struct LewkowiczSpeedResidual¶

Defined in File fit.h

Struct Documentation¶

struct LewkowiczSpeedResidual¶

Cost functor for fitting lewkowicz speed profiles.

Implements operator() as required by ceres-solver. Allows some parameters to be fixed.

Internal use only (see fit_lewkowicz_speed).

Public Functions

LewkowiczSpeedResidual(double z, double u, const Array5b &fixed_params, const Array5d &initial_params)¶

Construct the cost functor.

Parameters, where supplied in input initial_params and masked in input fixed_params are in this order:

pi_coles = params[0]
kappa = params[1]
u_inf = params[2]
shear_ratio = params[3]
delta_c = params[4]

Parameters

z: Observation data point, vertical coordinate in m
u: Observation data point, horizontal speed in m/s
fixed_params: Eigen::Array<bool, 5, 1>, logical mask, true where the parameter is fixed, false where it is allowed to vary
initial_params: Eigen::Array<double, 5, 1> containing initial parameter values. These are used where fixed parameters are specified

template <typename T> bool operator()(T const *const *parameters, T *residual) const¶

Struct PowerLawSpeedResidual¶

Defined in File fit.h

Struct Documentation¶

struct PowerLawSpeedResidual¶

Cost functor for fitting power law speed profiles.

Templated for automatic differentiation.

Internal use only (see fit_power_law_speed).

Public Functions

PowerLawSpeedResidual(double z, double u, double z_ref = 1.0, double u_ref = 1.0)¶

template <typename T> bool operator()(const T *const alpha, T *residual) const¶

Class AdemData¶

Defined in File adem.h

Class Documentation¶

class AdemData¶

Data container for ADEM input parameters and results.

Upcoming refactor

The AdemData class is treated as a struct, whose contents are updated by a number of functions… But, the functions operating on this should be refactored into class methods. That way, appropriate validation of the contents can be undertaken prior to application of each function.

This refactor is captured in issue #34.

Public Functions

void load(std::string file_name, bool print_var = true)¶

Load data from a *.mat file containing eddy signature data.

Parameters

file_name: File name (including relative or absolute path)
print_var: Boolean, default true. Print variables as they are read in (not advised except for debugging!)

void save(std::string filename)¶

Save eddy signature data to a *.mat file.

Parameters

filename: File name (including relative or absolute path)

Public Members

std::vector<std::string> eddy_types = {"A", "B1+B2+B3+B4"}¶: Eddy types used to create the results.

double beta¶

double delta_c¶: Atmospheric boundary layer thickness \( \delta_c \) [m].

double kappa¶: von Karman constant \( \kappa \). Typically \( \kappa = 0.41 \).

double pi_coles¶: Coles wake parameter \( \Pi \).

double shear_ratio¶: Ratio between free-stream and skin friction velocities \( S = U_{inf}/U_{\tau} \).

double u_inf¶: Free-stream velocity \( U_{inf}|_{z = \delta_c} \)) [m/s].

double u_tau¶: Skin friction velocity [m/s].

double zeta¶: Scaled streamwise derivative \( \zeta \) of the Coles wake parameter \( \Pi \).

Eigen::VectorXd z¶: Vertical coordinates used in the analysis [m].

Eigen::VectorXd eta¶: Nondimensionalised vertical coordinates used in the analysis \( \eta = z/\delta_{c} \).

Eigen::VectorXd lambda_e¶: Parameterised nondimensional vertical coordinates used in the analysis.

Eigen::VectorXd u_horizontal¶: Horizontal mean velocity varying with vertical coordinate [m/s].

Eigen::ArrayXXd reynolds_stress¶: Reynolds stress profiles from all eddy types.

Eigen::ArrayXXd reynolds_stress_a¶: Reynolds Stress profiles, contributions from Type A eddies only.

Eigen::ArrayXXd reynolds_stress_b¶: Reynolds Stress profiles, contributions from Type B eddies only.

Eigen::ArrayXXd k1z¶: Wavenumbers \( k_{1}z \) for which spectra are defined at each vertical coordinate.

Eigen::Tensor<double, 3> psi¶: Turbulent Spectra \( \Psi \) (corresponding to wavenumber k1z) at each vertical coordinate.

Eigen::Tensor<double, 3> psi_a¶: Turbulent Spectra (corresponding to wavenumber k1z) at each vertical coordinate from Type A eddies only.

Eigen::Tensor<double, 3> psi_b¶: Turbulent Spectra (corresponding to wavenumber k1z) at each vertical coordinate from Type B eddies only.

Eigen::VectorXd t2wa¶: Negated convolution function \( -(T^2)\omega \) encapsulating variation of eddy strength and scale for Type A eddies.

Eigen::VectorXd t2wb¶: Negated convolution function \( -(T^2)\omega \) encapsulating variation of eddy strength and scale for Type B eddies.

Eigen::VectorXd residual_a¶: Fit residuals from the deconvolution of t2wa

Eigen::VectorXd residual_b¶: Fit residuals from the deconvolution of t2wb

Class BasicLidar¶

Defined in File data_types.h

Class Documentation¶

class BasicLidar¶

Public Functions

void read(mat_t *matfp, bool print_var = true)¶

Public Members

const std::string type = "lidar_basic"¶

VectorXd t¶

VectorXd z¶

Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> u¶

Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> v¶

Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> w¶

Eigen::Vector3d position¶

double half_angle¶

std::string t

std::string z

std::string u

std::string v

std::string w

std::string position

std::string half_angle

struct es::BasicLidar::[anonymous] units¶

Class Bins¶

Defined in File profile.h

Class Documentation¶

class Bins¶

Public Functions

Bins(std::vector<double> z)¶

Bins(std::vector<double> z, std::vector<double> dx, std::vector<double> dy, std::vector<double> dz)¶

Bins(std::vector<double> z, double half_angle_degrees)¶

Bins(std::vector<double> z, double half_angle_degrees, std::vector<double> dz)¶

~Bins()¶

Public Members

unsigned long n_bins¶

std::vector<double> z_ctrs¶

std::vector<double> dx¶

std::vector<double> dy¶

std::vector<double> dz¶

Class EddySignature¶

Defined in File signature.h

Class Documentation¶

class EddySignature¶

Data container for Eddy signature tensors.

Public Functions

void load(std::string file_name, bool print_var = false)¶

Load data from a *.mat file containing eddy signature data.

TODO overload with load(std::vector<std::string> file_names, bool print_var = false){} to load and average multiple signature files

Parameters

file_name: File name (including relative or absolute path)
print_var: Boolean, default true. Print variables as they are read in (not advised except for debugging!)

void save(std::string filename)¶

Save eddy signature data to a *.mat file.

Note

NOT IMPLEMENTED YET

Parameters

filename: File name (including relative or absolute path)

EddySignature operator+(const EddySignature &c) const¶

Define overloaded + (plus) operator for eddy signatures.

Return

A new EddySignature() with combined signatures of the two eddy types.

Parameters

c: The EddySignature to add.

EddySignature operator/(double denom) const¶

Define overloaded / (divide) operator for eddy signatures.

Return

A new EddySignature() whose signature (g, j) is element-wise divided by input denom.

Parameters

denom: A number to divide by

Public Members

std::string eddy_type¶: Eddy types used to create the results.

Eigen::VectorXd lambda¶: Mapped vertical coordinates used in the analysis (e.g. 1 x 50)

Eigen::ArrayXXd k1z¶: Wavenumber (wavenumber space for each vertical coord, e.g. 50 x 801)

Eigen::Tensor<double, 3> g¶: g (6 coefficients at each vertical coord and wavenumber, e.g 50 x 801 x 6)

Eigen::ArrayXXd j¶: J (6 coefficients at each vertical coord, e.g 50 x 6)

Template Class Profile¶

Defined in File profile.h

Class Documentation¶

template <class ProfileType> class Profile¶

Public Functions

Profile(const Bins &bins)¶

Profile(const Bins &bins, double x, double y, double z)¶

~Profile()¶

const std::vector<ProfileType> &getValues() const¶

void setValues(const std::vector<ProfileType> &values)¶

Public Members

template<> double position[3]¶

Bins bins¶

Template Class Reader¶

Defined in File readers.h

Class Documentation¶

template <class DataType> class Reader¶

Public Functions

Reader(const std::string &file)¶

~Reader()¶

std::string logString() const¶

DataType *read(bool print = false)¶

void checkFileType(int level)¶

void checkTimeseries(int level)¶

double getWindowDuration() const¶

void setWindowDuration(double windowDuration)¶

int getWindowSize() const¶

void setWindowSize(int windowSize)¶

void readWindow(const int index)¶

Protected Attributes

std::string file¶

mat_t *matfp¶

std::string file_type = std::string("none")¶

int windowSize¶

double windowDuration¶

DataType data¶

Class VelocityProfile¶

Defined in File profile.h

Inheritance Relationships¶

Base Type¶

public es::Profile< double > (Template Class Profile)

Class Documentation¶

class VelocityProfile : public es::Profile<double>¶

Public Functions

VelocityProfile()¶

virtual ~VelocityProfile()¶

Class spectra¶

Defined in File spectra.h

Class Documentation¶

class spectra¶

Enums¶

Enum file_type_check_level¶

Defined in File readers.h

Enum Documentation¶

enum es::file_type_check_level¶

Values:

NONE = 0¶

OAS_STANDARD = 1¶

Enum timeseries_check_level¶

Defined in File readers.h

Enum Documentation¶

enum es::timeseries_check_level¶

Values:

PRESENT = 0¶

INCREASING = 1¶

MONOTONIC = 2¶

STRICTLY_MONOTONIC = 3¶

Functions¶

Function es::adem¶

Defined in File adem.h

Function Documentation¶

AdemData es::adem(const double beta, const double delta_c, const double kappa, const double pi_coles, const double shear_ratio, const double u_inf, const double zeta, const EddySignature &signature_a, const EddySignature &signature_b)¶

Compute full turbulent properties from given Attached-Detached Eddy Model parameters.

Uses the using the Lewkowicz (1982) formulation (Perry and Marusic eq.9) of the Coles wake function to determine u_h(z), spectra and Reynolds Stresses from input parameters.

Parameters

beta: The Clauser parameter, representing acceleration/decelaration of the boundary layer
delta_c: The boundary layer thickness in m
kappa: The von Karman constant, typically 0.41.
pi_coles: Coles wake parameter Pi
shear_ratio: Ratio between free stream and friction velocity S = u_inf/u_tau
u_inf: The free stream speed in m/s
zeta: Represents a scaled streamwise derivative of Pi
signature_a: EddySignature class loaded with data for type A eddies
signature_b: EddySignature class loaded with data for an ensemble of type B eddies

Template Function es::array_size¶

Defined in File variable_readers.h

Function Documentation¶

template <size_t SIZE, class T> size_t es::array_size(T (&arr)[SIZE])¶

Return the size (i.e. number of elements) in an array variable Comes from https://coderwall.com/p/nb9ngq/better-getting-array-size-in-c

Example: int arr[] = {1, 2, 3, 4, 5}; std::cout << array_size(arr) << std::endl;

Function es::checkVariableType¶

Defined in File variable_readers.h

Function Documentation¶

void es::checkVariableType(matvar_t *mat_var, int matvar_type)¶

Template Function es::coles_wake¶

Defined in File velocity.h

Function Documentation¶

template <typename T_z, typename T_param> T_z es::coles_wake(T_z const &eta, T_param const &pi_coles)¶

Compute coles wake parameter.

Used by Perry and Marusic 1995.

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of z values.

\[\begin{split} W_c[\eta, \Pi_j] & = & 2 \eta^2 \left( 3 - 2\eta \right) - \frac{1}{3\Pi_j}\eta^3 \\ \eta & = & \frac{z+z_0}{\delta + z_0} \end{split}\]

Translated from MATLAB:

function wc = colesWake(eta, Pi)
    wc = 2*eta.^2.*(3-2*eta) - (1/Pi).*eta.^2.*(1-eta).*(1-2*eta);
end

Parameters

eta: Nondimensional height values
pi_coles: The coles wake parameter Pi

Function es::fit_lewkowicz_speed(const Eigen::ArrayXd&, const Eigen::ArrayXd&, const Array5b&, const Array5d&, bool)¶

Defined in File fit.h

Function Documentation¶

Array5d es::fit_lewkowicz_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, const Array5b &fixed_params, const Array5d &initial_params, bool print_report = true)

Determine best fit of Lewkowicz relation to input speed profile.

Robustly fits Lewkowicz speed relation \( u = f(z, \Pi, \kappa, U_{inf}, S, \delta_{c}) \) to input data \( u, z \).

An initial guess is required, and parameters can be fixed to help constrain the fitting process. See also fit_lewkowicz_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, bool print_report = true) for providing a set of default values and fixed parameters.

Parameters are as follows:

pi_coles = params[0]
kappa = params[1]
u_inf = params[2]
shear_ratio = params[3]
delta_c = params[4]

The fitting process uses a Levenberg-Marquadt solver (provided by ceres-solver) with Automatic differentiation (for improved numeric stability and convergence over numerical differentiation approaches), and a Cauchy Loss function of parameter 0.5 (for robustness agains outlying data entries).

Return

Array5d containing fitted values [pi_coles, kappa, u_inf, shear_ratio, delta_c]

Parameters

z: Vertical locations in m (or normalised, if z_ref = 1.0)
u: Observed speeds at corresponding vertical locations, in m/s
fixed_params: Eigen::Array<bool, 5, 1>, logical mask, true where the parameter is fixed, false where it is allowed to vary
initial_params: Eigen::Array<double, 5, 1> containing initial parameter values. These are used where fixed parameters are specified
print_report: bool, default true. If true, the solver prints iterations summary and full final convergence and solution report (for debugging and validation purposes).

Function es::fit_lewkowicz_speed(const Eigen::ArrayXd&, const Eigen::ArrayXd&, bool)¶

Defined in File fit.h

Function Documentation¶

Array5d es::fit_lewkowicz_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, bool print_report = true)

Determine best fit of Lewkowicz relation to input speed profile.

Robustly fits Lewkowicz speed relation \( u = f(z, \Pi, \kappa, U_{inf}, S, \delta_{c}) \) to input data \( u, z \), using ‘sensible’ initial guesses (for Atmospheric Boundary Layer in Northern European conditions) and fixed parameters to constrain the fitting process:

- von Karman constant \form#1 is fixed according to the value in definitions.h.
- Boundary layer thickness \form#13 is fixed to a default value of ``1000 m``.

Return

Array5d containing fitted values [pi_coles, kappa, u_inf, shear_ratio, delta_c]

Parameters

z: Vertical locations in m (or normalised, if z_ref = 1.0)
u: Observed speeds at corresponding vertical locations, in m/s
print_report: bool, default true. If true, the solver prints iterations summary and full final convergence and solution report (for debugging and validation purposes).

Function es::fit_power_law_speed¶

Defined in File fit.h

Function Documentation¶

double es::fit_power_law_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, const double z_ref = 1.0, const double u_ref = 1.0)¶

Determine best fit of power law to input speed profile.

Return

Parameters

z: Vertical locations in m (or normalised, if z_ref = 1.0)
u:
z_ref: Optional reference height (default 1.0) by which the input z is normalised
u_ref: Optional reference velocity (default 1.0) by which the input u is normalised

Template Function es::get_coles_wake¶

Defined in File stress.h

Function Documentation¶

template <typename T> T es::get_coles_wake(const T eta, const double pi_coles)¶: Get the Coles wake distribution \( w_{c} \) from the wake parameter \( \Pi \).

Function es::get_f_integrand¶

Defined in File stress.h

Function Documentation¶

Eigen::ArrayXd es::get_f_integrand(const Eigen::ArrayXd eta, const double kappa, const double pi_coles, const double shear_ratio)¶

Get the integrand \( f \) used in computation of \( R_{13} \) in the modified Lewkowicz method.

Uses equations 2 and 7 Perry and Marusic 1995 Part 1.

Function es::get_mean_speed¶

Defined in File adem.h

Function Documentation¶

void es::get_mean_speed(AdemData &data)¶

Get the mean speed profile and update the data structure with it.

Parameters

data:

Function es::get_reynolds_stresses¶

Defined in File adem.h

Function Documentation¶

void es::get_reynolds_stresses(AdemData &data, const EddySignature &signature_a, const EddySignature &signature_b)¶

Get the Reynolds Stress distributions from T2w and J distributions.

The ouptut Reynolds Stress matrix is of size output_dim_size = input_dim_size - kernel_dim_size + 1 (requires: input_dim_size >= kernel_dim_size). Legacy MATLAB equivalent is:

[R, RA, RB] = getReynoldsStresses(T2wA, T2wB, JA, JB);

Parameters

data:

Function es::get_spectra¶

Defined in File adem.h

Function Documentation¶

void es::get_spectra(AdemData &data, const EddySignature &signature_a, const EddySignature &signature_b)¶

Get the premultiplied power spectra.

Legacy MATLAB equivalent is:

[Psi, PsiA, PsiB] = getSpectra(T2wA, T2wB, gA, gB, U1, S);

Parameters

data:

Function es::get_t2w¶

Defined in File adem.h

Function Documentation¶

void es::get_t2w(AdemData &data, const EddySignature &signature_a, const EddySignature &signature_b)¶

Get the T^2w distributions from the eddy signatures by deconvolution.

These distributions are used for calculation of Spectra and Stress terms.

Updates t2wa, t2wb, lambda_e, residual_a and residual_b in the input data structure.

Legacy MATLAB equivalent is:

[t2wa, t2wb, lambda_e, residual_a, residual_b] = get_tw(pi_coles, s, beta, zeta, JA(:,3), JB(:,3), lambda);

Parameters

data:
signature_a:
signature_b:

Function es::getVariable¶

Defined in File variable_readers.h

Function Documentation¶

matvar_t *es::getVariable(mat_t *matfp, const std::string var_name, bool print_var = true, int max_rank = 2)¶

Get pointer to a variable and check validity. Optionally print variable and check rank.

Return

Parameters

matfp:
var_name:
print_var:
max_rank: The maximum rank of the variable. Default 2 as everything is an array in MATLAB.

Template Function es::lewkowicz_speed¶

Defined in File velocity.h

Function Documentation¶

template <typename T_z, typename T_param> T_z es::lewkowicz_speed(T_z const &z, T_param const &pi_coles, T_param const &kappa, T_param const &u_inf, T_param const &shear_ratio, T_param const &delta_c)¶

Compute Lewkowicz (1982) velocity profile.

Used by Perry and Marusic 1995 (from eqs 2 and 7).

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of z values.

Translated from MATLAB:

function [f] = getf(eta, Pi, kappa, S)
    f = (-1/kappa)*log(eta) + (Pi/kappa)*colesWake(1,Pi)*ones(size(eta)) - (Pi/kappa)*colesWake(eta,Pi);
    f(isinf(f)) = S;
end

Parameters

z: height in m (or Nondimensional heights (eta = z/delta_c) where delta_c = 1.0)
pi_coles: The coles wake parameter Pi
kappa: von Karman constant
u_inf: Speed of flow at z = delta (m/s)
shear_ratio: Shear / skin friction velocity ratio (shear_ratio = u_inf / u_tau)
delta_c: Boundary layer thickness in m, used to normalise z. Defaults to 1.0.

Template Function es::marusic_jones_speed¶

Defined in File velocity.h

Function Documentation¶

template <typename T_z, typename T_pi_j> T_z es::marusic_jones_speed(T_z const &z, T_pi_j const pi_j, const double kappa, const double z_0, const double delta, const double u_inf, const double u_tau)¶

Compute speed profile according to Marusic and Jones’ relations.

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of z values.

Speed is computed as:

\[\begin{split} \frac{\overline{U}}{U_{\tau}} & = & \frac{1}{\kappa} \ln \left( \frac{z+z_0}{k_s} \right) + Br + \frac{\Pi_j}{\kappa} W_c[\eta, \Pi_j] \\ W_c[\eta, \Pi_j] & = & 2 \eta^2 \left( 3 - 2\eta \right) - \frac{1}{3\Pi_j}\eta^3 \\ \eta & = & \frac{z+z_0}{\delta + z_0} \end{split}\]

which reduces to the defecit relation:

\[ U_{D}^{*} = \frac{U_\infty-\overline{U}}{U_{\tau}} = -\frac{1}{\kappa} \ln \left( \eta \right) + \frac{1}{3\kappa} \left(\eta^3 - 1 \right) + 2 \frac{\Pi_j}{\kappa} \left(1 - 3\eta^2 + 2\eta^3 \right) \]

Parameters

z: Height(s) in m at which you want to get speed.
pi_j: Jones’ modification of the Coles wake factor. Double or AutoDiffScalar type acccepted.
kappa: von Karman constant
z_0: Roughness length - represents distance of hypothetical smooth wall from actual rough wall z0 = 0.25k_s (m)
delta: Boundary layer thickness (m)
u_inf: Speed of flow at z = delta (m/s)
u_tau: Shear / skin friction velocity (governed by ratio parameter S = u_inf / u_tau)

Template Function es::most_law_speed¶

Defined in File velocity.h

Function Documentation¶

template <typename T> T es::most_law_speed(T const &z, const double kappa, const double d, const double z0, const double L)¶

Compute speed profile according to the MOST law.

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of z values.

MOST law speed is computed as:

TODO

Parameters

z: Height value(s) at which you want to get speed (m)
kappa: von Karman constant
d: Zero plane offset distance (e.g. for forest canopies) (m)
z0: Roughness length (m)
L: Monin-Obukhov length (m)

Function es::operator<<(std::ostream&, AdemData const&)¶

Defined in File adem.h

Function Documentation¶

std::ostream &es::operator<<(std::ostream &os, AdemData const &data)

Print information about the AdemData attributes to ostream using the << operator.

Return

os The same ostream, with the data representation added to the stream

Parameters

os: The ostream to print to
data: The AdemDataclass instance

Function es::operator<<(::std::ostream&, const Bins&)¶

Defined in File profile.cpp

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “es::operator<<” with arguments (::std::ostream&, const Bins&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- std::ostream &es::operator<<(std::ostream&, AdemData const&)
- std::ostream &es::operator<<(std::ostream&, const Bins&)
- template <class DataType>
  std::ostream &es::operator<<(std::ostream&, const Reader<DataType>&)
- template <class ProfileType>
  std::ostream &es::operator<<(std::ostream&, const Profile<ProfileType>&)

Template Function es::operator<<(::std::ostream&, const Profile<ProfileType>&)¶

Defined in File profile.h

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “es::operator<<” with arguments (::std::ostream&, const Profile<ProfileType>&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- std::ostream &es::operator<<(std::ostream&, AdemData const&)
- std::ostream &es::operator<<(std::ostream&, const Bins&)
- template <class DataType>
  std::ostream &es::operator<<(std::ostream&, const Reader<DataType>&)
- template <class ProfileType>
  std::ostream &es::operator<<(std::ostream&, const Profile<ProfileType>&)

Template Function es::operator<<(::std::ostream&, const Reader<DataType>&)¶

Defined in File readers.h

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “es::operator<<” with arguments (::std::ostream&, const Reader<DataType>&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- std::ostream &es::operator<<(std::ostream&, AdemData const&)
- std::ostream &es::operator<<(std::ostream&, const Bins&)
- template <class DataType>
  std::ostream &es::operator<<(std::ostream&, const Reader<DataType>&)
- template <class ProfileType>
  std::ostream &es::operator<<(std::ostream&, const Profile<ProfileType>&)

Template Function es::power_law_speed¶

Defined in File velocity.h

Function Documentation¶

template <typename T, typename Talf> T es::power_law_speed(T const &z, const double u_ref, const double z_ref, Talf const &alpha)¶

Compute speed profile according to the power law.

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of z values (to obtain shear profile u/dz) and/or alpha values (to obtain variational jacobean for fitting parameter alpha).

Power law speed is computed as:

\[ \frac{\overline{U}}{\overline{U}_{ref}} = \left(\frac{z}{z_{ref}}\right)^{\alpha} \]

Parameters

z: Height(s) in m at which you want to get speed.
u_ref: Reference speed in m/s
z_ref: Reference height in m
alpha: Power law exponent. Must be of same type as input z (allows autodifferentiation).

Function es::readArrayXXd¶

Defined in File variable_readers.h

Function Documentation¶

ArrayXXd es::readArrayXXd(mat_t *matfp, const std::string var_name, bool print_var)¶

Function es::readDouble¶

Defined in File variable_readers.h

Function Documentation¶

double es::readDouble(mat_t *matfp, const std::string var_name, bool print_var)¶

Function es::readString¶

Defined in File variable_readers.h

Function Documentation¶

std::string es::readString(mat_t *matfp, const std::string var_name, bool print_var)¶

Function es::readTensor3d¶

Defined in File variable_readers.h

Function Documentation¶

Tensor<double, 3> es::readTensor3d(mat_t *matfp, const std::string var_name, bool print_var)¶

Function es::readVector3d¶

Defined in File variable_readers.h

Function Documentation¶

Vector3d es::readVector3d(mat_t *matfp, const std::string var_name, bool print_var)¶

Function es::readVectorXd¶

Defined in File variable_readers.h

Function Documentation¶

VectorXd es::readVectorXd(mat_t *matfp, const std::string var_name, bool print_var)¶

Function es::reynolds_stress_13¶

Defined in File stress.h

Function Documentation¶

void es::reynolds_stress_13(Eigen::ArrayXd &r13_a, Eigen::ArrayXd &r13_b, const double beta, const Eigen::ArrayXd &eta, const double kappa, const double pi_coles, const double shear_ratio, const double zeta)¶

Compute Horizontal-vertical Reynolds Stress R13 profile using modified Lewkowicz formulation.

Gets Reynolds Stress profiles due to Type A and B eddies. Adopts the approach of Perry and Marusic 1995, using the modified Lewkowicz formulation (Perry and Marusic eq.51).

TODO - presently, the integrations here are sensitive to the chosen eta distribution. Refactor according to issue #38

TODO - Determine whether it will be an advantage to template this function so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of eta values.

References

[1] Perry AE and Marusic I (1995) A wall-wake model for turbulent boundary layers. Part 1. Extension of the attached eddy hypothesis J Fluid Mech vol 298 pp 361-388

Future Improvements

[1] Optional different mean profile formulations including account for the free surface

[2] Support directional variation with height; i.e. compatible with mean profile formulations using U(y)

[3] Added formulation for contribution of smaller Type C eddies

[4] For the given wall formulation, can f be calculated more efficiently? See schlichting and gersten p. 593

Parameters

beta: Clauser parameter \( \Beta \), representing acceleration/decelaration of the boundary layer
eta: Nondimensional vertical coordinates at which you want to get stress \( R_{13} \). Values must ascend but not necessarily be monotonic.
kappa: von Karman constant.
pi_coles: Coles wake parameter \( \Pi \)
shear_ratio: Ratio between free-stream and skin friction velocities \( S = U_{inf}/U_{\tau} \)
zeta: Scaled streamwise derivative \( \zeta \) of the Coles wake parameter \( \Pi \)

Template Function es::veer_lhs¶

Defined in File veer.h

Function Documentation¶

template <typename T> T es::veer_lhs(T const &ui, const double ui_g, const double phi)¶

Computes the left hand side of the veer relations given u(z) or v(z).

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of u values.

The veer relations are:

\[\begin{split} 2 |\mathbf{\Omega}| sin(\phi) (\overline{v}_g - \overline{v}) & = & \frac{\partial}{\partial z} \bigg( \nu\frac{\partial\overline{u}}{\partial z} - \overline{u'w'}\bigg) \\ 2 |\mathbf{\Omega}| sin(\phi) (\overline{u}_g - \overline{u}) & = & \frac{\partial}{\partial z} \bigg(\nu\frac{\partial\overline{v}}{\partial z} - \overline{v'w'}\bigg) \end{split}\]

Parameters

ui: Mean velocity component at a given height (m/s)
ui_g: Mean geostrophic velocity component outside the atmospheric boundary layer (m/s)
phi: Latitude (degrees)

Template Function es::veer_rhs¶

Defined in File veer.h

Function Documentation¶

template <typename T> T es::veer_rhs(T &ui, T &uiu3_bar, const double nu)¶

Computes the right hand side of the veer relations given u(z) or v(z).

See veer_lhs() for the full relation.

Templated so that it can be called with active scalars (allows use of autodiff), doubles/floats, Eigen::Arrays (directly) or Eigen::VectorXds (via template specialisation) of u values.

Parameters

ui: Mean velocity component at a given height (m/s){
uiu3_bar: Mean cross term of unsteady velocity components \( \overline{u_i'u_3'}} \) (m^2/s^2)
nu: Kinematic viscosity of the fluid (m^2/s)

Function main¶

Defined in File main.cpp

Function Documentation¶

Warning

doxygenfunction: Cannot find function “main” in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml

Function utilities::conv¶

Defined in File conv.h

Function Documentation¶

Eigen::VectorXd utilities::conv(const Eigen::VectorXd &input, const Eigen::VectorXd &kernel)¶

Convolve an input vector with a kernel, returning an output of the same length as the input.

Uses zero-padded fft based output

Parameters

out:
input:
kernel:

Template Function utilities::deconv¶

Defined in File filter.h

Function Documentation¶

template <typename DerivedIn, typename DerivedOut> void utilities::deconv(Eigen::ArrayBase<DerivedOut> &z, const Eigen::ArrayBase<DerivedIn> &b, const Eigen::ArrayBase<DerivedIn> &a)¶

One-dimensional deconvolution.

z = deconv(b, a) deconvolves vector a out of column b, returning a column vector

Template Function utilities::fft_next_good_size¶

Defined in File conv.h

Function Documentation¶

template <typename T> T utilities::fft_next_good_size(const T n)¶

Find the next good size for an fft, to pad with minimal number of zeros.

Return

M Length to pad the signal to, for optimal FFT performance

Parameters

N: Length of a signal

Template Function utilities::filter(Eigen::ArrayBase<DerivedOut>&, const Eigen::ArrayBase<DerivedIn>&, const Eigen::ArrayBase<DerivedIn>&, const Eigen::ArrayBase<DerivedIn>&)¶

Defined in File filter.h

Function Documentation¶

template <typename DerivedOut, typename DerivedIn> void utilities::filter(Eigen::ArrayBase<DerivedOut> &y, const Eigen::ArrayBase<DerivedIn> &b, const Eigen::ArrayBase<DerivedIn> &a, const Eigen::ArrayBase<DerivedIn> &x)

One-dimensional digital filter.

Filters the input 1d array (or std::vector<>) x with the filter described by coefficient arrays b and a.

See https://en.wikipedia.org/wiki/Digital_filter

If a(0) is not equal to 1, the filter coefficients are normalised by a(0).

a(0)*y(n-1) = b(0)*x(n-1) + b(1)*x(n-2) + … + b(nb)*x(n-nb-1) - a(1)*y(n-2) - … - a(na)*y(n-na-1)

Template Function utilities::filter(std::vector<T>&, const std::vector<T>&, const std::vector<T>&, const std::vector<T>&)¶

Defined in File filter.h

Function Documentation¶

template <typename T> void utilities::filter(std::vector<T> &y, const std::vector<T> &b, const std::vector<T> &a, const std::vector<T> &x)

Template Function utilities::tensor_dims¶

Defined in File tensors.h

Function Documentation¶

template <typename T> std::string utilities::tensor_dims(T &tensor)¶

Return a string representation of tensor dimensions.

Return

string The output string like “[2 x 3 x 4]” for a rank 3 tensor

Parameters

tensor: An Eigen::Tensor

Template Function utilities::trapz(Eigen::ArrayBase<OtherDerived> const&, const Eigen::ArrayBase<Derived>&)¶

Defined in File trapz.h

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “utilities::trapz” with arguments (Eigen::ArrayBase<OtherDerived> const&, const Eigen::ArrayBase<Derived>&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- template <typename Derived, typename OtherDerived>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< OtherDerived > const &, const Eigen::ArrayBase< Derived > &)
- template <typename Derived>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< Derived > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< DerivedOut > const &, const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)

Template Function utilities::trapz(const Eigen::ArrayBase<Derived>&)¶

Defined in File trapz.h

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “utilities::trapz” with arguments (const Eigen::ArrayBase<Derived>&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- template <typename Derived, typename OtherDerived>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< OtherDerived > const &, const Eigen::ArrayBase< Derived > &)
- template <typename Derived>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< Derived > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< DerivedOut > const &, const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)

Template Function utilities::trapz(Eigen::ArrayBase<DerivedOut> const&, const Eigen::ArrayBase<DerivedX>&, const Eigen::ArrayBase<DerivedY>&)¶

Defined in File trapz.h

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “utilities::trapz” with arguments (Eigen::ArrayBase<DerivedOut> const&, const Eigen::ArrayBase<DerivedX>&, const Eigen::ArrayBase<DerivedY>&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- template <typename Derived, typename OtherDerived>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< OtherDerived > const &, const Eigen::ArrayBase< Derived > &)
- template <typename Derived>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< Derived > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< DerivedOut > const &, const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)

Template Function utilities::trapz(const Eigen::ArrayBase<DerivedX>&, const Eigen::ArrayBase<DerivedY>&)¶

Defined in File trapz.h

Function Documentation¶

Warning

doxygenfunction: Unable to resolve multiple matches for function “utilities::trapz” with arguments (const Eigen::ArrayBase<DerivedX>&, const Eigen::ArrayBase<DerivedY>&) in doxygen xml output for project “My Project” from directory: ./doxyoutput/xml. Potential matches:

- template <typename Derived, typename OtherDerived>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< OtherDerived > const &, const Eigen::ArrayBase< Derived > &)
- template <typename Derived>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< Derived > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> utilities::trapz(const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)
- template <typename DerivedX, typename DerivedY, typename DerivedOut>
  EIGEN_STRONG_INLINE void utilities::trapz(Eigen::ArrayBase< DerivedOut > const &, const Eigen::ArrayBase< DerivedX > &, const Eigen::ArrayBase< DerivedY > &)

Defines¶

Define acosd¶

Defined in File definitions.h

Define Documentation¶

acosd(x)¶

Define asind¶

Defined in File definitions.h

Define Documentation¶

asind(x)¶

Define cosd¶

Defined in File definitions.h

Define Documentation¶

cosd(x)¶

Define KAPPA_VON_KARMAN¶

Defined in File definitions.h

Define Documentation¶

KAPPA_VON_KARMAN¶

Define OMEGA_WORLD¶

Defined in File definitions.h

Define Documentation¶

OMEGA_WORLD¶

Define sind¶

Defined in File definitions.h

Define Documentation¶

sind(x)¶

Define tand¶

Defined in File definitions.h

Define Documentation¶

tand(x)¶

Typedefs¶

Typedef Array5b¶

Defined in File fit.h

Typedef Documentation¶

typedef Eigen::Array<bool, 5, 1> Array5b¶

Typedef Array5d¶

Defined in File fit.h

Typedef Documentation¶

typedef Eigen::Array<double, 5, 1> Array5d¶

Directories¶

Directory source¶

Directory path: source

Subdirectories¶

Files¶

Directory adem¶

↰ Parent directory (source)

Directory path: source/adem

Files¶

Directory relations¶

↰ Parent directory (source)

Directory path: source/relations

Files¶

Directory utilities¶

↰ Parent directory (source)

Directory path: source/utilities

Files¶

File adem.h¶

↰ Parent directory (source/adem)

Contents

Definition (source/adem/adem.h)
Includes
Namespaces
Classes
Functions

Definition (source/adem/adem.h)¶

Program Listing for File adem.h¶

↰ Return to documentation for file (source/adem/adem.h)

/*
 * adem.h Implementation of the Attached-Detached Eddy Method
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef SOURCE_ADEM_H_
#define SOURCE_ADEM_H_

#include <boost/algorithm/string/classification.hpp> // Include boost::for is_any_of
#include <boost/algorithm/string/split.hpp> // Include for boost::split
#include <Eigen/Dense>
#include <Eigen/Core>
#include <stdexcept>
#include <unsupported/Eigen/CXX11/Tensor>
#include <unsupported/Eigen/FFT>

#include "cpplot.h"

#include "adem/signature.h"
#include "profile.h"
#include "variable_readers.h"
#include "relations/stress.h"
#include "relations/velocity.h"
#include "utilities/conv.h"
#include "utilities/filter.h"
#include "utilities/tensors.h"

using namespace utilities;


namespace es {


class AdemData {
public:

    std::vector<std::string> eddy_types = {"A", "B1+B2+B3+B4"};

    double beta;

    double delta_c;

    double kappa;

    double pi_coles;

    double shear_ratio;

    double u_inf;

    double u_tau;

    double zeta;

    Eigen::VectorXd z;

    Eigen::VectorXd eta;

    Eigen::VectorXd lambda_e;

    Eigen::VectorXd u_horizontal;

    Eigen::ArrayXXd reynolds_stress;

    Eigen::ArrayXXd reynolds_stress_a;

    Eigen::ArrayXXd reynolds_stress_b;

    Eigen::ArrayXXd k1z;

    Eigen::Tensor<double, 3> psi;

    Eigen::Tensor<double, 3> psi_a;

    Eigen::Tensor<double, 3> psi_b;

    Eigen::VectorXd t2wa;

    Eigen::VectorXd t2wb;

    Eigen::VectorXd residual_a;

    Eigen::VectorXd residual_b;

    void load(std::string file_name, bool print_var = true) {
        std::cout << "Reading adem data from file " << file_name << std::endl;

        // Open the MAT file for reading
        mat_t *matfp = Mat_Open(file_name.c_str(), MAT_ACC_RDONLY);
        if (matfp == NULL) {
            std::string msg = "Error reading MAT file: ";
            throw std::invalid_argument(msg + file_name);
        }

        // Use the variable readers to assist
        k1z                 = readArrayXXd(matfp, "k1z", print_var);
        beta                = readDouble(matfp, "beta", print_var);
        delta_c             = readDouble(matfp, "delta_c", print_var);
        kappa               = readDouble(matfp, "kappa", print_var);
        pi_coles            = readDouble(matfp, "pi_coles", print_var);
        shear_ratio         = readDouble(matfp, "shear_ratio", print_var);
        u_inf               = readDouble(matfp, "u_inf", print_var);
        u_tau               = readDouble(matfp, "u_tau", print_var);
        zeta                = readDouble(matfp, "zeta", print_var);
        z                   = readVectorXd(matfp, "z", print_var);
        eta                 = readVectorXd(matfp, "eta", print_var);
        lambda_e            = readVectorXd(matfp, "lambda_e", print_var);
        u_horizontal        = readVectorXd(matfp, "u_horizontal", print_var);
        reynolds_stress     = readArrayXXd(matfp, "reynolds_stress", print_var);
        reynolds_stress_a   = readArrayXXd(matfp, "reynolds_stress_a", print_var);
        reynolds_stress_b   = readArrayXXd(matfp, "reynolds_stress_b", print_var);
        k1z                 = readArrayXXd(matfp, "k1z", print_var);
        psi                 = readTensor3d(matfp, "psi", print_var);
        psi_a               = readTensor3d(matfp, "psi_a", print_var);
        psi_b               = readTensor3d(matfp, "psi_b", print_var);
        t2wa                = readVectorXd(matfp, "t2wa", print_var);
        t2wb                = readVectorXd(matfp, "t2wb", print_var);
        residual_a          = readVectorXd(matfp, "residual_a", print_var);
        residual_b          = readVectorXd(matfp, "residual_b", print_var);

        // Special handling to split the string of types into a vector of strings
        std::string typestr = readString(matfp, "eddy_types", print_var);
        boost::split(eddy_types, typestr, boost::is_any_of(", "), boost::token_compress_on);

        // Close the file
        Mat_Close(matfp);
        std::cout << "Finished reading adem data" << std::endl;
    }

    void save(std::string filename) {
        std::cout << "Writing signature data..." << std::endl;
        throw std::invalid_argument("Error writing mat file - function not implemented");
    }

};


std::ostream &operator<<(std::ostream &os, AdemData const &data) {
    os << "AdemData() with fields:" << std::endl;
    os << "    eddy_types:        ";
    for (std::vector<std::string>::const_iterator i = data.eddy_types.begin(); i != data.eddy_types.end(); ++i) {
        os << *i << ", ";
    }
    os << std::endl;
    os << "    beta:              " << data.beta << std::endl
       << "    delta_c:           " << data.delta_c << std::endl
       << "    kappa:             " << data.kappa << std::endl
       << "    pi_coles:          " << data.pi_coles << std::endl
       << "    shear_ratio:       " << data.shear_ratio << std::endl
       << "    u_inf:             " << data.u_inf << std::endl
       << "    u_tau:             " << data.u_tau << std::endl
       << "    zeta:              " << data.zeta << std::endl
       << "    z:                 [" << data.z.size() << " x 1]" << std::endl
       << "    eta:               [" << data.eta.size() << " x 1]" << std::endl
       << "    lambda_e:          [" << data.lambda_e.size() << " x 1]" << std::endl
       << "    u_horizontal:      [" << data.u_horizontal.size() << " x 1]" << std::endl
       << "    reynolds_stress:   [" << data.reynolds_stress.rows() << " x " << data.reynolds_stress.cols() << "]" << std::endl
       << "    reynolds_stress_a: [" << data.reynolds_stress_a.rows() << " x " << data.reynolds_stress_a.cols() << "]" << std::endl
       << "    reynolds_stress_b: [" << data.reynolds_stress_b.rows() << " x " << data.reynolds_stress_b.cols() << "]" << std::endl
       << "    k1z:               [" << data.k1z.rows() << " x " << data.k1z.cols() << "]" << std::endl
       << "    psi:               [" << tensor_dims(data.psi) << "]" << std::endl
       << "    psi_a:             [" << tensor_dims(data.psi_a) << "]" << std::endl
       << "    psi_b:             [" << tensor_dims(data.psi_b) << "]" << std::endl
       << "    t2wa:              [" << data.t2wa.size() << " x 1]" << std::endl
       << "    t2wb:              [" << data.t2wb.size() << " x 1]" << std::endl
       << "    residual_a:        [" << data.residual_a.size() << " x 1]" << std::endl
       << "    residual_b:        [" << data.residual_b.size() << " x 1]" << std::endl;
    return os;
}


void get_t2w(AdemData& data, const EddySignature& signature_a, const EddySignature& signature_b) {

    // Define a range for lambda_e (lambda_e = ln(delta_c/z) = ln(1/eta))
    Eigen::ArrayXd lambda_e = Eigen::ArrayXd::LinSpaced(10001, 0, 100);

    // Re-express as eta and flip so that eta ascends (required for the integration in reynolds_stress_13)
    Eigen::ArrayXd eta;
    eta = -1.0 * lambda_e; // *-1 inverts 1/eta in the subsequent exp() operator
    eta = eta.exp();
    eta.reverseInPlace();

    // Get the Reynolds Stresses and flip back
    Eigen::ArrayXd r13a;
    Eigen::ArrayXd r13b;
    reynolds_stress_13(r13a, r13b, data.beta, eta, data.kappa, data.pi_coles, data.shear_ratio, data.zeta);
    r13a.reverseInPlace();
    r13b.reverseInPlace();

    // Extract J13 signature terms and trim so that array sizes match after the deconv
    auto len = signature_a.j.rows() - 2;
    Eigen::ArrayXd j13a = signature_a.j.col(2).tail(len);
    Eigen::ArrayXd j13b = signature_b.j.col(2).tail(len);

    // Deconvolve out the A and B structure contributions to the Reynolds Stresses
    // NOTE: it's actually -1*T^2w in this variable
    Eigen::ArrayXd minus_t2wa;
    Eigen::ArrayXd minus_t2wb;
    deconv(minus_t2wa, r13a, j13a);
    deconv(minus_t2wb, r13b, j13b);

    // TODO extend by padding out to the same length as lambda_e. The T^2w distributions converge to a constant at
    //  high lambda (close to the wall) so padding with the last value in the vector is physically valid. This will
    //  result in the Reynolds Stresses and Spectra, which are obtained by convolution, having the same number of points
    //  in the z direction as the axes variables (eta, lambda_e, etc)... very useful!

    // Store in the data object
    data.eta = eta;
    data.lambda_e = lambda_e;
    data.t2wa = minus_t2wa.matrix();
    data.t2wb = minus_t2wb.matrix();

}


void get_mean_speed(AdemData& data) {
    data.u_horizontal = lewkowicz_speed(data.eta, data.pi_coles, data.kappa, data.u_inf, data.u_tau);
}


void get_reynolds_stresses(AdemData& data, const EddySignature& signature_a, const EddySignature& signature_b){

    // Map the input signature data to tensors (shared memory)
    auto input_rows = data.t2wa.rows();
    auto kernel_rows = signature_a.j.rows() - 2; // removes the first two rows of j for stability

    // Compute cumulative length of input and kernel;
    auto N = input_rows + kernel_rows - 1 ;
    auto M = fft_next_good_size(N);

    // Allocate the output and set zero
    data.reynolds_stress_a = Eigen::ArrayXXd(M, 6);
    data.reynolds_stress_b = Eigen::ArrayXXd(M, 6);
    data.reynolds_stress_a.setZero();
    data.reynolds_stress_b.setZero();

    Eigen::FFT<double> fft;

    // Column-wise convolution of T2w with J
    // TODO refactor to use the conv() function
    for (int k = 0; k < 6; k++) {
        Eigen::VectorXd in_a1(M);
        Eigen::VectorXd in_a2(M);
        Eigen::VectorXd in_b1(M);
        Eigen::VectorXd in_b2(M);
        in_a1.setZero();
        in_a2.setZero();
        in_b1.setZero();
        in_b2.setZero();
        in_a1.topRows(input_rows) = data.t2wa.matrix();
        in_b1.topRows(input_rows) = data.t2wb.matrix();
        in_a2.topRows(kernel_rows) = signature_a.j.col(k).bottomRows(kernel_rows).matrix();
        in_b2.topRows(kernel_rows) = signature_b.j.col(k).bottomRows(kernel_rows).matrix();

        // Take the forward ffts
        Eigen::VectorXcd out_a1(M);
        Eigen::VectorXcd out_a2(M);
        Eigen::VectorXcd out_b1(M);
        Eigen::VectorXcd out_b2(M);
        fft.fwd(out_a1, in_a1);
        fft.fwd(out_a2, in_a2);
        fft.fwd(out_b1, in_b1);
        fft.fwd(out_b2, in_b2);

        // Convolve by element-wise multiplication
        Eigen::VectorXcd inter_a = out_a1.array() * out_a2.array();
        Eigen::VectorXcd inter_b = out_b1.array() * out_b2.array();

        // Inverse FFT to deconvolve
        Eigen::VectorXd out_a(M);
        Eigen::VectorXd out_b(M);
        fft.inv(out_a, inter_a);
        fft.inv(out_b, inter_b);
        data.reynolds_stress_a.col(k) = out_a;
        data.reynolds_stress_b.col(k) = out_b;
    }

    // Trim the zero-padded ends
    data.reynolds_stress_a = data.reynolds_stress_a.topRows(data.lambda_e.rows());
    data.reynolds_stress_b = data.reynolds_stress_b.topRows(data.lambda_e.rows());
    data.reynolds_stress = data.reynolds_stress_a + data.reynolds_stress_b;

}


void get_spectra(AdemData& data, const EddySignature& signature_a, const EddySignature& signature_b){

    // Initialise the output spectrum tensors
    Eigen::array<Eigen::Index, 3> dims = signature_a.g.dimensions();
    Eigen::array<Eigen::Index, 3> psi_dims = {data.t2wa.rows(), dims[1], dims[2]};
    Eigen::Tensor<double, 3> psi_a(psi_dims);
    Eigen::Tensor<double, 3> psi_b(psi_dims);

    // Map the t2w arrays as vectors, for use with the conv function
    Eigen::Map<Eigen::VectorXd> t2wa_vec(data.t2wa.data(), data.t2wa.rows());
    Eigen::Map<Eigen::VectorXd> t2wb_vec(data.t2wb.data(), data.t2wb.rows());

    // For each of the 6 auto / cross spectra terms
    for (Eigen::Index j = 0; j < dims[2]; j++) {

        auto page_offset_sig = j * dims[0] * dims[1];
        auto page_offset_psi = j * psi_dims[0] * psi_dims[1];

        // For each of the different heights
        for (Eigen::Index i = 0; i < dims[1]; i++) {

            auto elements_offset_sig = (page_offset_sig + i * dims[0]);
            auto elements_offset_psi = (page_offset_psi + i * psi_dims[0]);

            Eigen::Map<Eigen::VectorXd> g_a_vec((double *)signature_a.g.data() + elements_offset_sig, dims[0]);
            Eigen::Map<Eigen::VectorXd> g_b_vec((double *)signature_b.g.data() + elements_offset_sig, dims[0]);

            Eigen::Map<Eigen::VectorXd> psi_a_vec((double *)psi_a.data() + elements_offset_psi, psi_dims[0]);
            Eigen::Map<Eigen::VectorXd> psi_b_vec((double *)psi_b.data() + elements_offset_psi, psi_dims[0]);

            psi_a_vec = conv(t2wa_vec, g_a_vec);
            psi_b_vec = conv(t2wb_vec, g_b_vec);

        }
    }

    // Premultiply by the u_tau^2 term (see eq. 43)
    auto u_tau_sqd = pow(data.u_tau, 2.0);

    // Don't store the sum of the spectra - they're memory hungry and we can always add the two together
    data.psi_a = u_tau_sqd * psi_a;
    data.psi_b = u_tau_sqd * psi_b;

}


AdemData adem(const double beta,
              const double delta_c,
              const double kappa,
              const double pi_coles,
              const double shear_ratio,
              const double u_inf,
              const double zeta,
              const EddySignature& signature_a,
              const EddySignature& signature_b){

    // Data will contain computed outputs and useful small variables from the large signature files
    AdemData data = AdemData();

    // Add k1z to the data structure
    data.k1z = signature_a.k1z;

    // Add parameter inputs to data structure
    data.beta = beta;
    data.delta_c = delta_c;
    data.kappa = kappa;
    data.pi_coles = pi_coles;
    data.shear_ratio = shear_ratio;
    data.u_inf = u_inf;
    data.u_tau = data.u_inf / data.shear_ratio;
    data.zeta = zeta;

    // Deconvolve for T^2w and update the data structure with the outputs
    get_t2w(data, signature_a, signature_b);

    // Get vertical points through the boundary layer (for which the convolution functions were defined)
    Eigen::VectorXd eta = data.lambda_e;
    eta.array().exp();
    eta.array().inverse();
    data.eta = eta;
    data.z = eta.array() * data.delta_c;

    // Get the mean speed profile at the same vertical points and update the data structure
    get_mean_speed(data);

    // Determine Reynolds Stresses by convolution and add them to the data structure
    get_reynolds_stresses(data, signature_a, signature_b);

    // Determine Spectra by convolution and add them to the data structure
    get_spectra(data, signature_a, signature_b);

    return data;

}

} /* namespace es */

#endif /* SOURCE_ADEM_H_ */

Includes¶

Eigen/Core
Eigen/Dense
adem/signature.h (File signature.h)
boost/algorithm/string/classification.hpp
boost/algorithm/string/split.hpp
cpplot.h
profile.h (File profile.h)
relations/stress.h (File stress.h)
relations/velocity.h (File velocity.h)
stdexcept
unsupported/Eigen/CXX11/Tensor
unsupported/Eigen/FFT
utilities/conv.h (File conv.h)
utilities/filter.h (File filter.h)
utilities/tensors.h (File tensors.h)
variable_readers.h (File variable_readers.h)

Namespaces¶

Namespace es

Classes¶

Class AdemData

Functions¶

File conv.h¶

↰ Parent directory (source/utilities)

Contents

Definition (source/utilities/conv.h)
Includes
Included By
Namespaces
Functions

Definition (source/utilities/conv.h)¶

Program Listing for File conv.h¶

↰ Return to documentation for file (source/utilities/conv.h)

/*
 * conv.h Matlab-like convolution
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_CONV_H
#define ES_FLOW_CONV_H

#include <Eigen/Dense>
#include <Eigen/Core>
#include <unsupported/Eigen/FFT>


namespace utilities {


template<typename T>
T fft_next_good_size(const T n) {
    T result = n;
    if (n <= 2) {
        result = 2;
        return result;
    }
    while (TRUE) {
        T m = result;
        while ((m % 2) == 0) m = m / 2;
        while ((m % 3) == 0) m = m / 3;
        while ((m % 5) == 0) m = m / 5;
        if (m <= 1) {
            return (result);
        }
        result = result + 1;
    }
}


Eigen::VectorXd conv(const Eigen::VectorXd &input, const Eigen::VectorXd &kernel) {

    // Map the input signature data to tensors (shared memory)
    auto input_len = input.rows();
    auto kernel_len = kernel.rows();

    // Compute cumulative length of input and kernel;
    auto N = input_len + kernel_len - 1;
    auto M = fft_next_good_size(N);

    // Create padded FFT inputs
    Eigen::FFT<double> fft;
    Eigen::VectorXd input_padded(M);
    Eigen::VectorXd kernel_padded(M);
    input_padded.setZero();
    kernel_padded.setZero();
    input_padded.topRows(input_len) = input;
    kernel_padded.topRows(kernel_len) = kernel;

    // Take the forward ffts
    Eigen::VectorXcd input_transformed(M);
    Eigen::VectorXcd kernel_transformed(M);
    fft.fwd(input_transformed, input_padded);
    fft.fwd(kernel_transformed, kernel_padded);

    // Convolve by element-wise multiplication
    Eigen::VectorXcd inter = input_transformed.array() * kernel_transformed.array();

    // Inverse FFT to deconvolve
    Eigen::VectorXd out(M);
    out.setZero();
    fft.inv(out, inter);

    // Crop to the size of the input vector
    out = out.topRows(input_len);
    return out;
}

}  /* namespace utilities */

#endif //ES_FLOW_CONV_H

Includes¶

Eigen/Core
Eigen/Dense
unsupported/Eigen/FFT

Included By¶

Namespaces¶

Namespace utilities

Functions¶

File cumtrapz.h¶

↰ Parent directory (source/utilities)

Contents

Definition (source/utilities/cumtrapz.h)
Includes
Included By

Definition (source/utilities/cumtrapz.h)¶

Program Listing for File cumtrapz.h¶

↰ Return to documentation for file (source/utilities/cumtrapz.h)

/*
 * cumtrapz.h Cumulative Trapezoidal Numerical Integration for Eigen
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_CUMTRAPZ_H
#define ES_FLOW_CUMTRAPZ_H

#include <Eigen/Dense>


namespace es {


template<typename Derived, typename OtherDerived>
EIGEN_STRONG_INLINE void cumsum(Eigen::ArrayBase<OtherDerived> const & out, const Eigen::ArrayBase<Derived>& in)
{
    Eigen::ArrayBase<OtherDerived>& out_ = const_cast< Eigen::ArrayBase<OtherDerived>& >(out);
    out_ = in;
    for (int i = 0 ; i < out_.cols() ; ++i)
    {
        for (int j = 1 ; j < out_.rows() ; ++j)
        {
            out_.coeffRef(j,i) += out_.coeff(j-1,i);
        }
    }
}
// Remove template specialisation from doc (causes duplicate) @cond
template<typename Derived>
EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> cumsum(const Eigen::ArrayBase<Derived>& y)
{
    Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> z;
    cumsum(z,y);
    return z;
}
// @endcond


template<typename Derived, typename OtherDerived>
EIGEN_STRONG_INLINE void cumtrapz(Eigen::ArrayBase<OtherDerived> const & out, const Eigen::ArrayBase<Derived>& in)
{
    Eigen::ArrayBase<OtherDerived>& out_ = const_cast< Eigen::ArrayBase<OtherDerived>& >(out);
    out_.derived().setZero(in.rows(), in.cols());
    cumsum(out_.bottomRows(out_.rows()-1), (in.topRows(in.rows()-1) + in.bottomRows(in.rows()-1)) * 0.5);
}
// Remove template specialisation from doc (causes duplicate) @cond
template<typename Derived>
EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> cumtrapz(const Eigen::ArrayBase<Derived>& y)
{
    Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> z;
    cumtrapz(z,y);
    return z;
}
// @endcond


template<typename DerivedX, typename DerivedY, typename DerivedOut>
EIGEN_STRONG_INLINE void cumtrapz(Eigen::ArrayBase<DerivedOut> const & out, const Eigen::ArrayBase<DerivedX>& in_x, const Eigen::ArrayBase<DerivedY>& in_y)
{
    // Input size check
    eigen_assert(in_x.rows() == in_y.rows());
    eigen_assert(in_x.cols() == 1);

    // Get dx for each piece of the integration
    Eigen::Array<typename Eigen::ArrayBase<DerivedX>::Scalar, Eigen::ArrayBase<DerivedX>::RowsAtCompileTime, Eigen::ArrayBase<DerivedX>::ColsAtCompileTime> dx;
    dx = (in_x.bottomRows(in_x.rows()-1) - in_x.topRows(in_x.rows()-1));

    // Get the average heights of the trapezoids
    Eigen::Array<typename Eigen::ArrayBase<DerivedY>::Scalar, Eigen::ArrayBase<DerivedY>::RowsAtCompileTime, Eigen::ArrayBase<DerivedY>::ColsAtCompileTime> inter;
    inter = (in_y.topRows(in_y.rows()-1) + in_y.bottomRows(in_y.rows()-1)) * 0.5;

    // Multiply by trapezoid widths to give areas of the trapezoids.
    // NB Broadcasting with *= only works for arrayX types, not arrayXX types like dx
    //  inter *= dx;
    for (int i = 0 ; i < inter.cols() ; ++i)
    {
        for (int j = 0 ; j < inter.rows() ; ++j)
        {
            inter(j,i) *= dx(j,0);
        }
    }

    // Initialise output
    Eigen::ArrayBase<DerivedOut>& out_ = const_cast< Eigen::ArrayBase<DerivedOut>& >(out);
    out_.derived().setZero(in_y.rows(), in_y.cols());

    // Perform the cumulative sum down the columns
    cumsum(out_.bottomRows(out_.rows()-1), inter);
}
// Remove template specialisation from doc (causes duplicate) @cond
template<typename DerivedX, typename DerivedY, typename DerivedOut>
EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> cumtrapz(const Eigen::ArrayBase<DerivedX>& x, const Eigen::ArrayBase<DerivedY>& y)
{
    Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> z;
    cumtrapz(z, x, y);
    return z;
}

}  /* namespace utilities */

#endif //ES_FLOW_CUMTRAPZ_H

Includes¶

Eigen/Dense

Included By¶

File stress.h

File data_types.h¶

↰ Parent directory (source)

Contents

Definition (source/data_types.h)
Includes
Included By
Namespaces
Classes

Definition (source/data_types.h)¶

Program Listing for File data_types.h¶

↰ Return to documentation for file (source/data_types.h)

/*
 * data_types.h Definitions for different input data types
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_DATA_TYPES_H
#define ES_FLOW_DATA_TYPES_H

#include <string>
#include "matio.h"
#include <eigen3/Eigen/Core>
#include "variable_readers.h"


namespace es {


class BasicLidar {
public:
    const std::string type = "lidar_basic";
    VectorXd t;
    VectorXd z;
    Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> u;
    Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> v;
    Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> w;
    Eigen::Vector3d position;
    double half_angle;
    struct {
        std::string t;
        std::string z;
        std::string u;
        std::string v;
        std::string w;
        std::string position;
        std::string half_angle;
    }units;

    void read(mat_t *matfp, bool print_var = true) {

        // Straightforward reads
        t = readVectorXd(matfp, "t", print_var);
        z = readVectorXd(matfp, "z", print_var);
        u = readArrayXXd(matfp, "u", print_var);
        v = readArrayXXd(matfp, "v", print_var);
        w = readArrayXXd(matfp, "w", print_var);
        half_angle = readDouble(matfp, "half_angle", print_var);

        // Handle initialisation of position as a two element vector, zero padded (elevation = 0) and as a three
        // element vector.
        Eigen::VectorXd pos = readVectorXd(matfp, "position", print_var);
        if (pos.size() == 2) {
            position = Vector3d(pos(0), pos(1), 0.0);
        }else {
            position = Vector3d(pos(0), pos(1), pos(2));
        }

        // TODO read in and tests on units structure

    }

};

} /* namespace es */

#endif // ES_FLOW_DATA_TYPES_H

Includes¶

eigen3/Eigen/Core
matio.h
string
variable_readers.h (File variable_readers.h)

Included By¶

File readers.h

Namespaces¶

Namespace es

Classes¶

Class BasicLidar

File definitions.h¶

↰ Parent directory (source)

Contents

Definition (source/definitions.h)
Included By
Defines

Definition (source/definitions.h)¶

Program Listing for File definitions.h¶

↰ Return to documentation for file (source/definitions.h)

/*
 * definitions.h Constants, definitions and identities
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_DEFINITIONS_H
#define ES_FLOW_DEFINITIONS_H

// Rotational speed of the world in rad/s
#define OMEGA_WORLD 7.2921159e-05

// von Karman constant
#define KAPPA_VON_KARMAN 0.41

// Degrees based trigonometry
#define sind(x) (sin(fmod((x), 360.0) * M_PI / 180.0))
#define cosd(x) (cos(fmod((x), 360.0) * M_PI / 180.0))
#define asind(x) (asin(x) * 180.0 / M_PI)
#define acosd(x) (acos(x) * 180.0 / M_PI)
#define tand(x) (tan(x * M_PI / 180.0) * M_PI / 180.0)

#endif // ES_FLOW_DEFINITIONS_H

Included By¶

Defines¶

File filter.h¶

↰ Parent directory (source/utilities)

Contents

Definition (source/utilities/filter.h)
Includes
Included By
Namespaces
Functions

Definition (source/utilities/filter.h)¶

Program Listing for File filter.h¶

↰ Return to documentation for file (source/utilities/filter.h)

/*
 * filter.h One-dimensional polynomial based digital filter for eigen arrayXd
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_FILTER_H
#define ES_FLOW_FILTER_H

#include <vector>
#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Core>
#include <stdexcept>


namespace utilities {


template<typename DerivedOut, typename DerivedIn>
void filter(Eigen::ArrayBase<DerivedOut> &y,
            const Eigen::ArrayBase<DerivedIn> &b,
            const Eigen::ArrayBase<DerivedIn> &a,
            const Eigen::ArrayBase<DerivedIn> &x) {

    // Check for 1D inputs
    eigen_assert((a.cols() == 1) || (a.rows() == 1));
    eigen_assert((b.cols() == 1) || (a.rows() == 1));
    eigen_assert((x.cols() == 1) || (a.rows() == 1));

    // Initialise output
    Eigen::ArrayBase<DerivedOut> &y_ = const_cast< Eigen::ArrayBase<DerivedOut> & >(y);
    y_.derived().setZero(x.rows(), x.cols());

    // Direct Form II transposed standard difference equation
    typename DerivedOut::Scalar tmp;
    Eigen::Index i;
    Eigen::Index j;
    for (i = 0; i < x.size(); i++) {
        tmp = 0.0;
        j = 0;
        y_(i) = 0.0;

        for (j = 0; j < b.size(); j++) {
            if (i - j < 0) continue;
            tmp += b(j) * x(i - j);
        }

        for (j = 1; j < a.size(); j++) {
            if (i - j < 0) continue;
            tmp -= a(j) * y_(i - j);
        }

        tmp /= a(0);
        y_(i) = tmp;
    }
}

template<typename T>
void filter(std::vector<T> &y, const std::vector<T> &b, const std::vector<T> &a, const std::vector<T> &x) {

    y.resize(0);
    y.resize(x.size());

    for (int i = 0; i < x.size(); i++) {
        T tmp = 0.0;
        int j = 0;
        y[i] = 0.0;
        for (j = 0; j < b.size(); j++) {
            if (i - j < 0) continue;
            tmp += b[j] * x[i - j];
        }

        for (j = 1; j < a.size(); j++) {
            if (i - j < 0) continue;
            tmp -= a[j] * y[i - j];
        }

        tmp /= a[0];
        y[i] = tmp;
    }
}

template<typename DerivedIn, typename DerivedOut>
void deconv(Eigen::ArrayBase<DerivedOut> &z,
            const Eigen::ArrayBase<DerivedIn> &b,
            const Eigen::ArrayBase<DerivedIn> &a) {

    eigen_assert(a(0) != 0);
    eigen_assert((a.cols() == 1) || (a.rows() == 1));
    eigen_assert((b.cols() == 1) || (b.rows() == 1));

    auto na = a.size();
    auto nb = b.size();

    // Initialise output
    Eigen::ArrayBase<DerivedOut> &z_ = const_cast< Eigen::ArrayBase<DerivedOut> & >(z);

    // Cannot deconvolve a longer signal out of a smaller one
    if (na > nb) {
        z_.derived().setZero(1, 1);
        return;
    }
    z_.derived().setZero(nb - na + 1, 1);

    // Deconvolution and polynomial division are the same operations as a digital filter's impulse response B(z)/A(z):
    Eigen::Array<typename DerivedOut::Scalar, Eigen::Dynamic, 1> impulse;
    impulse.setZero(nb - na + 1, 1);
    impulse(0) = 1;
    filter(z_, b, a, impulse);
}

}  /* namespace utilities */

#endif //ES_FLOW_FILTER_H

Includes¶

Eigen/Core
Eigen/Dense
iostream
stdexcept
vector

Included By¶

Namespaces¶

Namespace utilities

Functions¶

File fit.h¶

↰ Parent directory (source)

Contents

Definition (source/fit.h)
Includes
Namespaces
Classes
Functions
Typedefs

Definition (source/fit.h)¶

Program Listing for File fit.h¶

↰ Return to documentation for file (source/fit.h)

/*
 * Fit.h Use ceres-solver to fit mean speed and spectra to measured data
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_FIT_H
#define ES_FLOW_FIT_H

#include <stdio.h>
#include "ceres/ceres.h"
#include <Eigen/Core>

#include "relations/velocity.h"
#include "definitions.h"

using ceres::AutoDiffCostFunction;
using ceres::DynamicAutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;
using ceres::CauchyLoss;


typedef Eigen::Array<bool, 5, 1> Array5b;
typedef Eigen::Array<double, 5, 1> Array5d;


namespace es {


struct PowerLawSpeedResidual {
    PowerLawSpeedResidual(double z, double u, double z_ref=1.0, double u_ref=1.0) : z_(z), u_(u), u_ref_(u_ref), z_ref_(z_ref) {}
    template <typename T> bool operator()(const T* const alpha, T* residual) const {

        residual[0] = u_ - power_law_speed(T(z_), u_ref_, z_ref_, alpha[0]);
        return true;
    }
private:
    // Observations for a sample
    const double u_;
    const double z_;
    const double u_ref_;
    const double z_ref_;
};


double fit_power_law_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, const double z_ref=1.0, const double u_ref=1.0) {

    // Define the variable to solve for with its initial value. It will be mutated in place by the solver.
    double alpha = 0.3;

    // TODO Assert that z.size() == u.size()

    // Build the problem
    Problem problem;
    Solver::Summary summary;
    Solver::Options options;

    options.max_num_iterations = 400;
    options.linear_solver_type = ceres::DENSE_QR;
    options.minimizer_progress_to_stdout = true;
//std::cout << options << std::endl;
    // Set up the only cost function (also known as residual), using cauchy loss function for
    // robust fitting and auto-differentiation to obtain the derivative (jacobian)
    for (auto i = 0; i < z.size(); i++) {
        CostFunction *cost_function = new AutoDiffCostFunction<PowerLawSpeedResidual, 1, 1>(
            new PowerLawSpeedResidual(z[i], u[i], z_ref, u_ref));
        problem.AddResidualBlock(cost_function, new CauchyLoss(0.5), &alpha);
    }

    // Run the solver
    Solve(options, &problem, &summary);
    std::cout << summary.BriefReport() << "\n";
    std::cout << "Initial alpha: " << 0.3 << "\n";
    std::cout << "Final   alpha: " << alpha << "\n";

    return alpha;
}


struct LewkowiczSpeedResidual {

    LewkowiczSpeedResidual(double z, double u, const Array5b &fixed_params, const Array5d &initial_params) : z_(z), u_(u), fixed_params_(fixed_params), initial_params_(initial_params) { }

    template <typename T> bool operator()(T const* const* parameters, T* residual) const {

        // Mask out fixed and variable parameters
        // TODO Once Eigen 3.4.x comes out, reduce this code to use the parameter mask as a logical index.
        //  See this issue which will allow such indexing in eigen: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=329#c27
        int param_ctr = 0;
        std::vector<T> params(5);
        for (auto i = 0; i < 5; i++) {
            if (fixed_params_(i)) {
                params[i] = T(initial_params_(i));
            } else {
                params[i] = parameters[param_ctr][0];
                param_ctr += 1;
            }
        }
        T spd = lewkowicz_speed(T(z_), params[0], params[1], params[2], params[3], params[4]);
        residual[0] = u_ - spd;
        return true;
    }
private:
    const double u_;
    const double z_;
    Array5b fixed_params_;
    Array5d initial_params_;
};


Array5d fit_lewkowicz_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, const Array5b &fixed_params, const Array5d &initial_params, bool print_report = true) {

    // Initialise a results array, and get a vector of pointers to the 'unfixed' parameters (the ones to optimise)
    Array5d final_params(initial_params);
    std::vector<double *> unfixed_param_ptrs;
    for (int p_ctr = 0; p_ctr <5; p_ctr++) {
        if (!fixed_params(p_ctr)) { unfixed_param_ptrs.push_back(final_params.data()+p_ctr); };
    }

    // Build the problem
    Problem problem;
    Solver::Summary summary;
    Solver::Options options;
    options.max_num_iterations = 400;
    options.minimizer_progress_to_stdout = print_report;

    for (auto i = 0; i < z.size(); i++) {

        // Set the stride such that the entire set of derivatives is computed at once (since we have maximum 5)
        auto cost_function = new DynamicAutoDiffCostFunction<LewkowiczSpeedResidual, 5>(
            new LewkowiczSpeedResidual(z[i], u[i], fixed_params, initial_params)
        );

        // Add N parameters, where N is the number of free parameters in the problem
        auto pc = 0;
        for (auto p_ctr = 0; p_ctr < 5; p_ctr++) {
            if (!fixed_params(p_ctr)) {
                cost_function->AddParameterBlock(1);
                pc += 1;
            }
        }
        cost_function->SetNumResiduals(1);
        problem.AddResidualBlock(cost_function, new CauchyLoss(0.5), unfixed_param_ptrs);
    }

    Solve(options, &problem, &summary);

    if (print_report) {
        std::cout << summary.FullReport() << std::endl;
        std::cout << "Initial values: " << initial_params.transpose() << std::endl;
        std::cout << "Final   values: " << final_params.transpose() << std::endl;
    }

    return final_params;
}


Array5d fit_lewkowicz_speed(const Eigen::ArrayXd &z, const Eigen::ArrayXd &u, bool print_report = true) {

    // Set default initial parameters
    double pi_coles = 0.5;
    double kappa = KAPPA_VON_KARMAN;
    double shear_ratio = 20;
    double u_inf = u.maxCoeff();
    double delta_c = 1000;
    Array5d initial_params;
    initial_params << pi_coles, kappa, u_inf, shear_ratio, delta_c;

    // Set default fixed parameters (von karman constant and boundary layer thickness)
    Array5b fixed_params;
    fixed_params << false, true, false, false, true;

    return fit_lewkowicz_speed(z, u, fixed_params, initial_params, print_report);

}

} /* namespace es */


#endif // ES_FLOW_FIT_H

Includes¶

Eigen/Core
ceres/ceres.h
definitions.h (File definitions.h)
relations/velocity.h (File velocity.h)
stdio.h

Namespaces¶

Namespace es

Classes¶

Functions¶

Typedefs¶

File how_to.h¶

↰ Parent directory (source)

Contents

Definition (source/how_to.h)

Definition (source/how_to.h)¶

Program Listing for File how_to.h¶

↰ Return to documentation for file (source/how_to.h)

/*
 * how_to.h  Chunks of code that show us how to do things
 *
 *  May not be complete, may not compile. Not included in any build targets. Just copy, paste, adapt to our own uses.
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */


// Remove template specialisation from doc (causes duplicate) @cond

/***********************************************************************************************************************
 * HOW TO INTEGRATE USING THE NUMERICALINTEGRATION LIBRARY
 **********************************************************************************************************************/
//
#include <Eigen/Dense>
#include <Eigen/Core>
#include <iostream>
#include "NumericalIntegration.h"

/* Integrator for the velocity deficit
 *  We consider the example from:
 *
 *       http://www.gnu.org/software/gsl/manual/html_node/Numerical-integration-examples.html
 *
 *       int_0^1 x^{-1/2} log(x) dx = -4
 *
 *  The integrator expects the user to provide a functor as shown below.
 */

template<typename Scalar>
class IntegrandExampleFunctor
{
public:
    IntegrandExampleFunctor(const Scalar alpha):m_alpha(alpha)
    {
        assert(alpha>0);
    }

    Scalar operator()(const Scalar x) const
    {
        assert(x>0);
        return log(m_alpha*x) / sqrt(x);
    }

    void setAlpha(const Scalar alpha)
    {
        m_alpha = alpha;
    }
private:
    Scalar m_alpha;
};

double do_something(const double arg)
{
    // Define the scalar
    typedef double Scalar;

    // Define the functor
    Scalar alpha=1.;
    IntegrandExampleFunctor<Scalar> inFctr(alpha);

    //define the integrator
    Eigen::Integrator<Scalar> eigIntgtor(200);

    //define a quadrature rule
    Eigen::Integrator<Scalar>::QuadratureRule quadratureRule = Eigen::Integrator<Scalar>::GaussKronrod61;

    //define the desired absolute and relative errors
    Scalar desAbsErr = Scalar(0.);
    Scalar desRelErr = Eigen::NumTraits<Scalar>::epsilon() * 50.;

    //integrate
    Scalar result = eigIntgtor.quadratureAdaptive(inFctr, Scalar(0.),Scalar(1.), desAbsErr, desRelErr, quadratureRule);

    //expected result
    Scalar expected = Scalar(-4.);

    //print output
    size_t outputPrecision  = 18;
    std::cout<<std::fixed;
    std::cout<<"result          = "<<std::setprecision(outputPrecision)<<result<<std::endl;
    std::cout<<"exact result    = "<<std::setprecision(outputPrecision)<<expected<<std::endl;
    std::cout<<"actual error    = "<<std::setprecision(outputPrecision)<<(expected-result)<<std::endl;

    return 0.0;
}

TEST_F(MyTest, test_integrator_example) {

// Test a basic integrator out
double arg = 0.0;
double res = do_something(arg);

}



/***********************************************************************************************************************
 * HOW TO FIT USING CERES-SOLVER
 **********************************************************************************************************************/


//#include <stdlib.h>
//#include <stdio.h>
//#include <unistd.h>
//#include <stdbool.h>
//#include "matio.h"
#include "ceres/ceres.h"

#include <Eigen/Core>


using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;

// A templated cost functor that implements the residual r = 10 -
// x. The method operator() is templated so that we can then use an
// automatic differentiation wrapper around it to generate its
// derivatives.
struct CostFunctor {
    template <typename T> bool operator()(const T* const x, T* residual) const {
        residual[0] = T(10.0) - x[0];
        return true;
    }
};

void my_fitting_func() {

    // Run the ceres solver

    // Define the variable to solve for with its initial value. It will be mutated in place by the solver.
    double x = 0.5;
    const double initial_x = x;

    // Build the problem.
    Problem problem;

    // Set up the only cost function (also known as residual). This uses
    // auto-differentiation to obtain the derivative (jacobian).
    CostFunction *cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
    problem.AddResidualBlock(cost_function, NULL, &x);

    // Run the solver!
    Solver::Options ceroptions;
    ceroptions.minimizer_progress_to_stdout = true;
    Solver::Summary summary;
    Solve(ceroptions, &problem, &summary);
    std::cout << summary.BriefReport() << "\n";
    std::cout << "x : " << initial_x << " -> " << x << "\n";

}


/***********************************************************************************************************************
 * HOW TO DO AN FFT USING EIGEN
 **********************************************************************************************************************/


#include <Eigen/Core>
#include <unsupported/Eigen/FFT>


void my_fft_func() {
    size_t dim_x = 28, dim_y = 126;
    Eigen::FFT<float> fft;
    Eigen::MatrixXf in = Eigen::MatrixXf::Random(dim_x, dim_y);
    Eigen::MatrixXcf out;
    out.setZero(dim_x, dim_y);

    for (int k = 0; k < in.rows(); k++) {
        Eigen::VectorXcf tmpOut(dim_x);
        fft.fwd(tmpOut, in.row(k));
        out.row(k) = tmpOut;
    }

    for (int k = 0; k < in.cols(); k++) {
        Eigen::VectorXcf tmpOut(dim_y);
        fft.fwd(tmpOut, out.col(k));
        out.col(k) = tmpOut;
    }
}


/***********************************************************************************************************************
 * HOW TO DO AN FFT USING MKL DIRECTLY
 **********************************************************************************************************************/

#include "mkl_dfti.h"

void my_mkl_fft() {
    // Run an example FFT

    // Make meaningless data and a size vector
    float xf[200][100];
    MKL_LONG len[2] = {200, 100};

    // Create a decriptor, which is a pattern for what operation an FFT will undertake
    DFTI_DESCRIPTOR_HANDLE fft;
    DftiCreateDescriptor (&fft, DFTI_SINGLE, DFTI_REAL, 2, len);
    DftiCommitDescriptor(fft);

    // Compute a forward transform, in-place on the data
    DftiComputeForward(fft, xf);

    // Free the descriptor
    DftiFreeDescriptor(&fft);

    std::cout << "FFT COMPLETE\n";
}


/***********************************************************************************************************************
 * HOW TO PLOT REYNOLDS STRESS PROFILES, EDDY SIGNATURES AND STRENGTH / SCALE DISTRIBUTIONS WITH CPPLOT
 **********************************************************************************************************************/

#include "cpplot.h"

void my_rs_plotting_func() {

    // Plot the Reynolds Stress profiles (comes from get_spectra)
    cpplot::Figure figa = cpplot::Figure();
    for (auto i = 0; i < 6; i++) {
        cpplot::ScatterPlot p = cpplot::ScatterPlot();
        p.x = Eigen::VectorXd::LinSpaced(data.reynolds_stress_a.rows(), 1, data.reynolds_stress_a.rows());
        p.y = data.reynolds_stress_a.col(i).matrix();
        figa.add(p);
    }
    figa.write("test_t2w_rij_a.json");
    //    legend({'R13A'; 'R13B'})
    //    xlabel('\lambda_E')

    cpplot::Figure figb = cpplot::Figure();
    for (auto i = 0; i < 6; i++) {
        cpplot::ScatterPlot p = cpplot::ScatterPlot();
        p.x = Eigen::VectorXd::LinSpaced(data.reynolds_stress_b.rows(), 1, data.reynolds_stress_b.rows());
        p.y = data.reynolds_stress_b.col(i).matrix();
        figb.add(p);
    }
    figb.write("test_t2w_rij_b.json");
    //    legend({'R13A'; 'R13B'})
    //    xlabel('\lambda_E')


    // Plot the eddy signatures
    cpplot::Figure fig2 = cpplot::Figure();
    cpplot::ScatterPlot ja = cpplot::ScatterPlot();
    ja.x = Eigen::VectorXd::LinSpaced(j13a.rows(), 1, j13a.rows());
    ja.y = j13a.matrix();
    fig2.add(ja);
    cpplot::ScatterPlot jb = cpplot::ScatterPlot();
    jb.x = ja.x;
    jb.y = j13b.matrix();
    fig2.add(jb);
    fig2.write("test_t2w_j13ab.json");
    //    legend({'J13A'; 'J13B'})

    // Plot strength and scale distributions
    cpplot::Figure fig3 = cpplot::Figure();
    cpplot::ScatterPlot twa = cpplot::ScatterPlot();
    twa.x = Eigen::VectorXd::LinSpaced(minus_t2wa.rows(), 1, minus_t2wa.rows());
    twa.y = minus_t2wa.matrix();
    fig3.add(twa);
    cpplot::ScatterPlot twb = cpplot::ScatterPlot();
    twb.x = twa.x;
    twb.y = minus_t2wb.matrix();
    fig3.add(twb);
    fig3.write("test_t2w_t2wab.json");
    //    legend({'T^2\omegaA'; 'T^2\omegaB'})

}

/***********************************************************************************************************************
 * HOW TO DO MAT FILE READ WRITE
 **********************************************************************************************************************/

#include "matio.h"

// Create the output file
mat_t *matfp;
matfp = Mat_CreateVer(options["o"].as<std::string>().c_str(), NULL, MAT_FT_MAT73);
if ( NULL == matfp ) {
std::string msg = "Error creating MAT file: ";
throw  msg + options["o"].as<std::string>();
}

// We haven't written any data - just close the file
Mat_Close(matfp);

std::cout << "MATIO TEST COMPLETE" << std::endl;




// @endcond

File main.cpp¶

↰ Parent directory (source)

Contents

Definition (source/main.cpp)
Includes
Functions

Definition (source/main.cpp)¶

Program Listing for File main.cpp¶

↰ Return to documentation for file (source/main.cpp)

/*
 * main.cpp Command line executable for es-flow
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2016-9 Octue Ltd. All Rights Reserved.
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <stdbool.h>
#include "glog/logging.h"
#include "cxxopts.hpp"


int main(int argc, char* argv[]) {

    try {

        // Handle the input parsing and create the program help page

        // Set the program name (for --help option display) and default option behaviour
        cxxopts::Options options("es-flow", "EnvironmentSTUDIO flow library wrapper");
        bool logging = false;

        // Define the command line options, arranged visually (in the --help output) in two groups:
        options.add_options()
                ("l,log-file",
                 "Switch on logging, optionally specify the directory to save logfiles.",
                 cxxopts::value<std::string>()->implicit_value("logs"), "FILE")
                ("h,help",
                 "Display program help.",
                 cxxopts::value<bool>(), "BOOL");

        options.add_options("Input / output file")
                ("i,input-file", "Name of the input file (read only).", cxxopts::value<std::string>(), "FILE")
                ("o,output-file",
                 "Name of the output results file. Warning - this file will be overwritten if it already exists.",
                 cxxopts::value<std::string>(), "FILE");

        // Parse the input options
        options.parse(argc, argv);

        if (options.count("help")) {
            std::cout << options.help({"", "Input / output file"}) << std::endl;
            exit(0);
        }

        if (options.count("l")) {
            logging = true;
        }

        if (logging) {
            FLAGS_logtostderr = false;
            FLAGS_minloglevel = 0;
            FLAGS_log_dir = options["l"].as<std::string>();
            std::cout << "Logging to: " << FLAGS_log_dir << std::endl;
            google::InitGoogleLogging(argv[0]);
        }


    } catch (const cxxopts::OptionException &e) {
        std::cout << "Error parsing options: " << e.what() << std::endl;
        exit(1);

    } catch (...) {
        // Handle logging of general exceptions
        auto eptr = std::current_exception();
        try {
            // This deliberately rethrows the caught exception in the current scope, so we can log it
            if (eptr) {
                std::rethrow_exception(eptr);
            }
        } catch(const std::exception& e) {
            std::cout << "Caught exception: " << e.what() << std::endl;
        }
        exit(1);

    }

    exit(0);

}

Includes¶

cxxopts.hpp
glog/logging.h
stdbool.h
stdio.h
stdlib.h
unistd.h

Functions¶

Function main

File profile.cpp¶

↰ Parent directory (source)

Contents

Definition (source/profile.cpp)
Includes
Namespaces
Functions

Definition (source/profile.cpp)¶

Program Listing for File profile.cpp¶

↰ Return to documentation for file (source/profile.cpp)

/*
 * profile.cpp Profile handling for (e.g.) Velocity, Reynolds Stress and Spectral Tensor Profile management
 *
 * Function definitions are included in the .h file, because of the templating constraint. See http://stackoverflow.com/questions/495021/why-can-templates-only-be-implemented-in-the-header-file
 *
 * Author:                   Tom Clark (thclark @ github)
 *
 * Copyright (c) 2016-9 Octue Ltd. All Rights Reserved.
 *
 */
#include "profile.h"
#include "math.h"
#include "definitions.h"


namespace es {

Bins::Bins(std::vector<double> z): z_ctrs(z) {
    // Construct from a spacing with default bin size from central difference of the spacing
    int i;
    n_bins = z_ctrs.size();
    dx = std::vector<double>(n_bins);
    dy = std::vector<double>(n_bins);
    dz = std::vector<double>(n_bins);
    dz[0] = fabs(z[1]-z[0]);
    for (i = 1; i < n_bins-1; i++) {
        dz[i] = 0.5*fabs(z[i+1]-z[i-1]);
    }
    dz[n_bins-1] = fabs(z[n_bins-1]-z[n_bins-2]);
    for (i = 0; i < n_bins-1; i++) {
        dx[i] = dz[i];
        dy[i] = dz[i];
    }
}

Bins::Bins(std::vector<double> z, std::vector<double> d_x, std::vector<double> d_y, std::vector<double> d_z): Bins(z) {
    // Construct from a spacing and bin size in x, y, z
    dx = d_x;
    dy = d_y;
    dz = d_z;
    if ((dx.size() != n_bins) ||  (dy.size() != n_bins) || (dz.size() != n_bins)){
        throw std::range_error("Length of dx, dy or dz does not match the number of bins");
    }
}

Bins::Bins(std::vector<double> z, double half_angle_degrees): Bins(z) {
    // Construct from spacing, with bin height from central differencing of the spacing and bin widths using a half angle
    int i;
    double tha;
    tha = tand(half_angle_degrees);
    dx = std::vector<double>(z_ctrs.size());
    dy = std::vector<double>(z_ctrs.size());
    for (i = 0; i < n_bins - 1; i++) {
        dx[i] = 2.0 * z[i] * tha;
        dy[i] = dx[i];
    }
}

Bins::Bins(std::vector<double> z, double half_angle_degrees, std::vector<double> d_z): Bins(z, half_angle_degrees) {
    // Construct from spacing, with bin height given and bin widths using a half angle
    int i;
    dz = d_z;
    if (dz.size() != n_bins){
        throw std::range_error("Length of dx, dy or dz does not match the number of bins");
    }
}

Bins::~Bins() {
    //Destructor
}

::std::ostream& operator<<(::std::ostream& os, const Bins& bins) {
    // Represent in logs or ostream
    return os << "debug statement for bins class";
}


/* Use profile base class constructor and destructor
VelocityProfile::VelocityProfile() {
// TODO Auto-generated constructor stub
}

VelocityProfile::~VelocityProfile() {
// TODO Auto-generated destructor stub
}

VectorXd VelocityProfile::AutoDiff() {
// Differentiate the velocity profile wrt z

    typedef Eigen::AutoDiffScalar<Eigen::VectorXd> AScalar;
    // AScalar stores a scalar and a derivative vector.

    // Instantiate an AutoDiffScalar variable with a normal Scalar
    double s = 0.3;
    AScalar As(s);

    // Get the value from the Instance
    std::cout << "value: " << As.value() << std::endl;

    // The derivative vector is As.derivatives();

    // Resize the derivative vector to the number of dependent variables
    As.derivatives().resize(2);

    // Set the initial derivative vector
    As.derivatives() = Eigen::VectorXd::Unit(2,0);
    std::cout << "Derivative vector : " << As.derivatives().transpose() << std::endl;

    // Instantiate another AScalar
    AScalar Ab(4);
    Ab.derivatives() = Eigen::VectorXd::Unit(2,1);

    // Do the most simple calculation
    AScalar Ac = As * Ab;

    std::cout << "Result/Ac.value()" << Ac.value() << std::endl;
    std::cout << "Gradient: " << Ac.derivatives().transpose() << std::endl;

    Eigen::VectorXd a(1);
    return a;
}
*/

} /* namespace es */

Includes¶

definitions.h (File definitions.h)
math.h
profile.h (File profile.h)

Namespaces¶

Namespace es

Functions¶

Function es::operator<<(::std::ostream&, const Bins&)

File profile.h¶

↰ Parent directory (source)

Contents

Definition (source/profile.h)
Includes
Included By
Namespaces
Classes
Functions

Definition (source/profile.h)¶

Program Listing for File profile.h¶

↰ Return to documentation for file (source/profile.h)

/*
 * profile.h Profile handling for (e.g.) Velocity, Reynolds Stress and Spectral Tensor Profile management
 *
 * Function definitions are included here, rather than the cpp file, because of the templating constraint. See http://stackoverflow.com/questions/495021/why-can-templates-only-be-implemented-in-the-header-file
 *
 * Author:                   Tom Clark (thclark @ github)
 *
 * Copyright (c) 2016-9 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_PROFILE_H
#define ES_FLOW_PROFILE_H

#include <iostream>
#include <vector>
#include <cstdlib>
#include <string>
#include <stdexcept>


namespace es {


class Bins {
public:

    // Construct from a spacing with bin sizes set by central differencing of spacing
    Bins(std::vector<double> z);

    // Construct from spacing with bin sizes specified
    Bins(std::vector<double> z, std::vector<double> dx, std::vector<double> dy, std::vector<double> dz);

    // Construct from spacing and a half angle in degrees, dz set by central differencing of spacing
    Bins(std::vector<double> z, double half_angle_degrees);

    // Construct from spacing and a half angle in degrees, with dz specified
    Bins(std::vector<double> z, double half_angle_degrees, std::vector<double> dz);

    // Destroy
    ~Bins();

    // Vectors of z locations and bin dimensions
    unsigned long n_bins;
    std::vector<double> z_ctrs;
    std::vector<double> dx;
    std::vector<double> dy;
    std::vector<double> dz;

};

// Represent Bins class in logs or ostream
::std::ostream& operator<<(::std::ostream& os, const Bins& bins);

template <class ProfileType>
class Profile {
private:
    std::vector<ProfileType> values;

public:

    // Construct using bins only (zero values and position)
    Profile(const Bins &bins);

    // Construct using bins and a global origin location
    Profile(const Bins &bins, double x, double y, double z);

    // Destroy
    virtual ~Profile();

    // Get the values, whatever type they might be
    const std::vector<ProfileType> &getValues() const;

    void setValues(const std::vector<ProfileType> &values);

    // Properties
    double position[3];
    Bins bins;

};

// Represent Profiles class in logs or ostream
template <class ProfileType>
::std::ostream& operator<<(::std::ostream& os, const Profile<ProfileType>& profile);

template <class ProfileType>
Profile<ProfileType>::Profile(const Bins &bins) : bins(bins) {
    // Construct from just bins, with default position {0.0, 0.0, 0.0}
    position[0] = 0.0;
    position[1] = 0.0;
    position[2] = 0.0;
}

template <class ProfileType>
Profile<ProfileType>::~Profile() {
    // Destructor
}

template <class ProfileType>
Profile<ProfileType>::Profile(const Bins &bins, double x, double y, double z): bins(bins) {
    // Construct from bins and a position
    position[0] = x;
    position[1] = y;
    position[2] = z;
}

template <class ProfileType>
void Profile<ProfileType>::setValues(const std::vector<ProfileType> &values) {
    if (values.size() != bins.n_bins) {
        throw std::out_of_range("size of vector 'values' does not equal the number of bins for this profile");
    }
    Profile::values = values;
}

template <class ProfileType>
const std::vector<ProfileType> &Profile<ProfileType>::getValues() const {
    return values;
}

//template <class ProfileType>
//const std::vector<ProfileType> &Profile<ProfileType>::getZDerivative() const {
//    // Get the first derivative with respect to Z location using Eigen's Autodiff
//    // Some examples can be found in eigen/unsupported/doc/examples/AutoDiff.cpp
//
//
//    return derivatives;
//}

template <class ProfileType>
::std::ostream& operator<<(::std::ostream& os, const Profile<ProfileType>& profile) {
    // Represent in logs or ostream
    return os << "debug statement for profile class";
}


class VelocityProfile: public Profile<double> {
public:
    VelocityProfile();
    virtual ~VelocityProfile();

//    VectorXd AutoDiff();
};


} /* namespace es */


#endif // ES_FLOW_PROFILE_H

Includes¶

cstdlib
iostream
stdexcept
string
vector

Included By¶

Namespaces¶

Namespace es

Classes¶

Functions¶

Template Function es::operator<<(::std::ostream&, const Profile<ProfileType>&)

File readers.h¶

↰ Parent directory (source)

Contents

Definition (source/readers.h)
Includes
Namespaces
Classes
Enums
Functions

Definition (source/readers.h)¶

Program Listing for File readers.h¶

↰ Return to documentation for file (source/readers.h)

/*
 * readers.h File readers for timeseries instrument data
 *
 * Author:                   Tom Clark (thclark @ github)
 *
 * Copyright (c) 2016-9 Octue Ltd. All Rights Reserved.
 *
 */
#ifndef ES_FLOW_READERS_H
#define ES_FLOW_READERS_H

#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <stdbool.h>
#include <iostream>
#include <vector>
#include <cstdlib>
#include <string>
#include <stdexcept>
#include "matio.h"
#include "glog/logging.h"
#include <iostream>
#include <vector>
#include <cstdlib>
#include <string>
#include <stdexcept>
#include <eigen3/Eigen/Core>
#include "data_types.h"


namespace es {


// Level of data validation applied to timeseries. Exception is thrown if checks fail
enum timeseries_check_level {
    PRESENT             = 0,    // Checks that data is present and of correct type. Extracts start and end timestamps.
    INCREASING          = 1,    // As PRESENT + checks that the timestamp always increases.
    MONOTONIC           = 2,    // As INCREASING + checks that timestamp uniformly increases, allows data skip (e.g. instrument turned off then back on later). Extracts sampling frequency.
    STRICTLY_MONOTONIC  = 3,    // AS INCREASING + checks that timestamp uniformly increases, data skipping causes errors. Extracts sampling frequency.
};

// Level of data validation applied to file type. Exception is thrown if checks fail
enum file_type_check_level {
    NONE                = 0,    // No checks
    OAS_STANDARD        = 1,    // Checks that the file contains the 'type' string variable
};

template <class DataType>
class Reader {
public:

    Reader(const std::string &file);

    ~Reader();

    std::string logString() const;

    DataType *read(bool print=false);

    void checkFileType(int level);

    void checkTimeseries(int level);

    double getWindowDuration() const;

    void setWindowDuration(double windowDuration);

    int getWindowSize() const;

    void setWindowSize(int windowSize);

    void readWindow(const int index);

protected:

    std::string file;
    mat_t *matfp;
    std::string file_type = std::string("none");
    int windowSize;
    double windowDuration;
    DataType data;

};

template <class DataType>
Reader<DataType>::Reader(const std::string &file) : file(file) {

    // Open the MAT file for reading and save the pointer
    matfp = Mat_Open(file.c_str(), MAT_ACC_RDONLY);
    if (matfp == NULL) {
        std::string msg = "Error reading MAT file: ";
        throw std::invalid_argument(msg + file);
    }

}

template <class DataType>
Reader<DataType>::~Reader() {

    // Close the file on destruction of the Reader
    Mat_Close(matfp);

}

template <class DataType>
void Reader<DataType>::checkFileType(int level){

    // Get the file type from the reserved 'type' variable in the mat file
    matvar_t *type_var = Mat_VarRead(matfp, "type");
    if (type_var == NULL) {
        if (level == OAS_STANDARD) {
            throw std::invalid_argument("Error reading mat file (most likely not an OAS standard file format - 'type' variable is missing)");
        }
    } else {
        if (type_var->class_type != MAT_C_CHAR) {
            throw std::invalid_argument("Error reading mat file ('type' variable must be a character array)");
        }
        file_type = std::string((const char*)type_var->data, type_var->dims[1]);
    }

}

template <class DataType>
void Reader<DataType>::checkTimeseries(int level) {

    // Check the integrity of the timeseries, automatically determine dt and frequency where possible

}

template <class DataType>
double Reader<DataType>::getWindowDuration() const {
    return windowDuration;
}

template <class DataType>
void Reader<DataType>::setWindowDuration(double windowDuration) {
    Reader<DataType>::windowDuration = windowDuration;
}

template <class DataType>
int Reader<DataType>::getWindowSize() const {
    return windowSize;
}

template <class DataType>
void Reader<DataType>::setWindowSize(int windowSize) {
    Reader<DataType>::windowSize = windowSize;
}

template <class DataType>
std::string Reader<DataType>::logString() const {
    return std::string("Object Reader for type ") + file_type + std::string(", attached to file ") + file;
}

template <class DataType>
DataType *Reader<DataType>::read(bool print) {
    // Simply invoke the read method of the DataType class.
    if (std::strcmp(data.type.c_str(), file_type.c_str())) {
        std::string msg = "Mat file type '" + file_type + "' incompatible with the specified data type '" + data.type + "'";
        throw std::invalid_argument(msg);
    }
    data.read(matfp, print);
    return &data;
}

template <class DataType>
void Reader<DataType>::readWindow(const int index) {
    //<DataType>.read(file_type, matfp);
}

// Represent Reader classes in logs or ostream
template <class DataType>
::std::ostream& operator<<(::std::ostream& os, const Reader<DataType>& reader) {
    // Represent in logs or ostream
    return os << reader.logString();
}

} /* namespace es */



#endif // ES_FLOW_READERS_H

Includes¶

cstdlib
data_types.h (File data_types.h)
eigen3/Eigen/Core
glog/logging.h
iostream
matio.h
stdbool.h
stdexcept
stdio.h
stdlib.h
string
unistd.h
vector

Namespaces¶

Namespace es

Classes¶

Template Class Reader

Enums¶

Functions¶

Template Function es::operator<<(::std::ostream&, const Reader<DataType>&)

File signature.h¶

↰ Parent directory (source/adem)

Contents

Definition (source/adem/signature.h)
Includes
Included By
Namespaces
Classes

Definition (source/adem/signature.h)¶

Program Listing for File signature.h¶

↰ Return to documentation for file (source/adem/signature.h)

/*
 * signature.h Eddy Signatures for use with the Attached-Detached Eddy Method
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef SOURCE_ADEM_SIGNATURE_H_
#define SOURCE_ADEM_SIGNATURE_H_

#include <boost/algorithm/string/classification.hpp> // Include boost::for is_any_of
#include <boost/algorithm/string/split.hpp> // Include for boost::split
#include <Eigen/Dense>
#include <Eigen/Core>
#include <stdexcept>
#include <unsupported/Eigen/CXX11/Tensor>

#include "variable_readers.h"
#include "profile.h"
#include "relations/stress.h"
#include "relations/velocity.h"
#include "utilities/filter.h"

#include "cpplot.h"
#include <unsupported/Eigen/FFT>
#include "utilities/conv.h"
#include "utilities/tensors.h"

using namespace utilities;


namespace es {


class EddySignature {
public:

    std::string eddy_type;

    Eigen::VectorXd lambda;

    Eigen::ArrayXXd k1z;

    Eigen::Tensor<double, 3> g;

    Eigen::ArrayXXd j;

    void load(std::string file_name, bool print_var = false) {
        std::cout << "Reading eddy signature data from file " << file_name << std::endl;

        // Open the MAT file for reading
        mat_t *matfp = Mat_Open(file_name.c_str(), MAT_ACC_RDONLY);
        if (matfp == NULL) {
            std::string msg = "Error reading MAT file: ";
            throw std::invalid_argument(msg + file_name);
        }

        // Use the variable readers to assist
        eddy_type = readString(matfp, "type", print_var);
        lambda = readVectorXd(matfp, "lambda", print_var);
        k1z = readArrayXXd(matfp, "k1z", print_var);
        g = readTensor3d(matfp, "g", print_var);
        j = readArrayXXd(matfp, "J", print_var);

        // Close the file
        Mat_Close(matfp);
        std::cout << "Finished reading eddy signature (Type " + eddy_type + ")" << std::endl;
    }

    void save(std::string filename) {
        std::cout << "Writing signature data..." << std::endl;
        throw std::invalid_argument("Error writing mat file - function not implemented");
    }

    EddySignature operator+(const EddySignature& c) const
    {
        EddySignature result;
        result.eddy_type = this->eddy_type + "+" + c.eddy_type;
        // TODO assert equality of lambda and k1z
        result.lambda = this->lambda;
        result.k1z = this->k1z;
        result.g = (this->g + c.g);
        result.j = (this->j + c.j);
        return result;
    }

    EddySignature operator/(double denom) const
    {
        EddySignature result;
        result.eddy_type = "(" + this->eddy_type + ")/" + std::to_string(denom);
        // TODO assert equality of lambda and k1z
        result.lambda = this->lambda;
        result.k1z = this->k1z;
        result.g = this->g;
        result.j = this->j;
        result.g = result.g / denom;
        result.j = result.j / denom;
        return result;
    }

};

} /* namespace es */

#endif /* SOURCE_ADEM_SIGNATURE_H_ */

Includes¶

Eigen/Core
Eigen/Dense
boost/algorithm/string/classification.hpp
boost/algorithm/string/split.hpp
cpplot.h
profile.h (File profile.h)
relations/stress.h (File stress.h)
relations/velocity.h (File velocity.h)
stdexcept
unsupported/Eigen/CXX11/Tensor
unsupported/Eigen/FFT
utilities/conv.h (File conv.h)
utilities/filter.h (File filter.h)
utilities/tensors.h (File tensors.h)
variable_readers.h (File variable_readers.h)

Included By¶

File adem.h

Namespaces¶

Namespace es

Classes¶

Class EddySignature

File spectra.cpp¶

↰ Parent directory (source/relations)

Contents

Definition (source/relations/spectra.cpp)
Includes

Definition (source/relations/spectra.cpp)¶

Program Listing for File spectra.cpp¶

↰ Return to documentation for file (source/relations/spectra.cpp)

/*
 * spectra.cpp Brief overview sentence
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#include "spectra.h"

Includes¶

spectra.h (File spectra.h)

File spectra.h¶

↰ Parent directory (source/relations)

Contents

Definition (source/relations/spectra.h)
Included By
Classes

Definition (source/relations/spectra.h)¶

Program Listing for File spectra.h¶

↰ Return to documentation for file (source/relations/spectra.h)

/*
 * spectra.h Brief overview sentence
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_SPECTRA_H
#define ES_FLOW_SPECTRA_H

class spectra {

};

#endif //ES_FLOW_SPECTRA_H

Included By¶

File spectra.cpp

Classes¶

Class spectra

File stress.h¶

↰ Parent directory (source/relations)

Contents

Definition (source/relations/stress.h)
Includes
Included By
Namespaces
Functions

Definition (source/relations/stress.h)¶

Program Listing for File stress.h¶

↰ Return to documentation for file (source/relations/stress.h)

/*
 * stress.h Atmospheric Boundary Layer Reynolds Stress relations
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_STRESS_H_
#define ES_FLOW_STRESS_H_

#include <Eigen/Dense>
#include <Eigen/Core>
#include <iostream>
#include <iomanip>
#include <math.h>
#include "NumericalIntegration.h"
#include "profile.h"
#include "utilities/trapz.h"
#include "utilities/cumtrapz.h"

using namespace utilities;


namespace es {


template<typename T>
T get_coles_wake(const T eta, const double pi_coles) {
    T wc;
    T eta_pow_2 = pow(eta, 2.0);
    wc = (2.0 * eta_pow_2 * (3.0 - (2.0 * eta)));
    wc -= (1.0 / pi_coles) * eta_pow_2 * (1.0 - eta) * (1.0 - (2.0 * eta));
    return wc;
}

// Remove template specialisation from doc (causes duplicate) @cond
Eigen::ArrayXd get_coles_wake(const Eigen::ArrayXd &eta, const double pi_coles) {
    Eigen::ArrayXd wc;
    Eigen::ArrayXd eta_pow_2;
    eta_pow_2 = eta.pow(2.0);
    wc = (2.0 * eta_pow_2 * (3.0 - (2.0 * eta)));
    wc -= (1.0 / pi_coles) * eta_pow_2 * (1.0 - eta) * (1.0 - (2.0 * eta));
    return wc;
}
// @endcond


Eigen::ArrayXd get_f_integrand(const Eigen::ArrayXd eta, const double kappa, const double pi_coles, const double shear_ratio) {

    Eigen::ArrayXd f;
    Eigen::ArrayXd ones;
    ones.setOnes(eta.size());
    f = (-1.0 / kappa) * eta.log()
        + (pi_coles/kappa) * get_coles_wake(ones, pi_coles) * ones
        - (pi_coles/kappa) * get_coles_wake(eta, pi_coles);
    for (int k = 0; k < f.size(); k++) {
        if (std::isinf(f[k])) {
            f(k) = shear_ratio; // from eqs 2 and 7 P&M part 1 1995
        }
    }

    return f;
}


void reynolds_stress_13(Eigen::ArrayXd &r13_a, Eigen::ArrayXd &r13_b, const double beta, const Eigen::ArrayXd &eta, const double kappa, const double pi_coles, const double shear_ratio, const double zeta){

    // Get f between eta = 0 and eta = 1 (bounds for C1 integration)
    Eigen::ArrayXd f = get_f_integrand(eta, kappa, pi_coles, shear_ratio);
    const double d_pi = 0.01 * pi_coles;
    Eigen::ArrayXd f_plus  = get_f_integrand(eta, kappa, (pi_coles + d_pi), shear_ratio);
    Eigen::ArrayXd f_minus  = get_f_integrand(eta, kappa, (pi_coles - d_pi), shear_ratio);

    // TODO can we cast this returned value directly?
    Eigen::ArrayXd c1_tmp;
    trapz(c1_tmp, eta, f);
    const double c1 = c1_tmp(0);

    // Do the integrations for ei with central differencing for numerical differentiation of e coefficients 5-7
    Eigen::ArrayXd e1;
    Eigen::ArrayXd e1_plus;
    Eigen::ArrayXd e1_minus;
    Eigen::ArrayXd e2;
    Eigen::ArrayXd e2_plus;
    Eigen::ArrayXd e2_minus;
    Eigen::ArrayXd e3;
    Eigen::ArrayXd e4;
    Eigen::ArrayXd e5;
    Eigen::ArrayXd e6;
    Eigen::ArrayXd e7;
    cumtrapz(e1, eta, f);
    cumtrapz(e1_plus, eta, f_plus);
    cumtrapz(e1_minus, eta, f_minus);
    cumtrapz(e2, eta, f.pow(2.0));
    cumtrapz(e2_plus, eta, f_plus.pow(2.0));
    cumtrapz(e2_minus, eta, f_minus.pow(2.0));
    e3 = f * e1;
    e4 = eta * f;
    e5 = (e1_plus - e1_minus) / (2.0 * d_pi);
    e6 = (e2_plus - e2_minus) / (2.0 * d_pi);
    e7 = f * e5;

    // Get A coefficients from equations A2 a-d
    Eigen::ArrayXd a1 = e2 - e3 + shear_ratio * (e4 - e1);
    Eigen::ArrayXd a2 = (-2.0 * e2) + e3 + (shear_ratio * e1);
    Eigen::ArrayXd a3 = e6 - e7 - (shear_ratio * e5);
    Eigen::ArrayXd a4 = (2.0 * e2) - e3 + (shear_ratio * e4) - (shear_ratio * 3.0 * e1);

    // B coefficients are simply evaluated at eta = 1 (eqn A6)
    // TODO eta must be defined up to 1! check this
    double b1 = a1(eta.rows()-1);
    double b2 = a2(eta.rows()-1);
    double b3 = a3(eta.rows()-1);
    double b4 = a4(eta.rows()-1);

    // E1 from (eqn A4). Can't call it E1 due to name conflict with above.
    double e1_coeff = 1.0 / (kappa * shear_ratio + 1.0);

    // N from (eqn A5) using central differencing as before
    double wc_minus = get_coles_wake(1.0, pi_coles - d_pi);
    double wc_plus =  get_coles_wake(1.0, pi_coles + d_pi);
    double n = get_coles_wake(1.0, pi_coles) + pi_coles * (wc_plus - wc_minus) / (2.0 * d_pi);

    // Compile f_i terms
    Eigen::ArrayXd f1;
    Eigen::ArrayXd f2;
    Eigen::ArrayXd f3;
    f1 = 1
        - a1 / (b1 + e1_coeff * b2)
        - e1_coeff * a2 / (b1 + e1_coeff * b2);

    f2 = (e1_coeff * n * a2 * b1
        + a3 * b1
        - e1_coeff * n * a1 * b2
        + e1_coeff * a3 * b2
        - a1 * b3
        - e1_coeff * a2 * b3) / (b1 + e1_coeff * b2);

    f3 = (e1_coeff * a2 * b1 + a4 * b1 - e1_coeff * a1 * b2 + e1_coeff * a4 * b2 - a1 * b4 - e1_coeff * a2 * b4) / (b1 + e1_coeff * b2);

    // Convert f2 and f3 into g1 and g2 ready for application of eq. 15
    Eigen::ArrayXd g1 = f2 / shear_ratio;
    Eigen::ArrayXd g2 = -f3 / (c1 * shear_ratio);

    // Top 3 boxes of figure 18
    Eigen::ArrayXd minus_r13 = f1 + g1 * zeta + g2 * beta;

    // Get the component due to equilibrium sink flow (OLD version - see P&M eqns 51,53)
    // minusReynStressA = ones(size(eta)) - eta + eta.*log(eta);

    // Lewkowicz 1982 shear stress for equilibrium sink flow, Perry and Marusic eqn. 51
    Eigen::ArrayXd minus_r13_a = 1.0 - (60.0/59.0)*eta - (20.0/59.0)*eta.pow(3.0) + (45.0/59.0)*eta.pow(4.0) - (24.0/59.0)*eta.pow(5.0) + (60.0/59.0)*eta*eta.log();

    // Handle the log(0) singularity
    for (int i = 0; i < minus_r13_a.size(); i++) {
        if (std::isinf(minus_r13_a(i))) {
            minus_r13_a(i) = 1.0;
        }
    }

    // Output correctly signed reynolds stress components
    r13_a = -1.0 * minus_r13_a;
    r13_b = -1.0 * (minus_r13 - minus_r13_a);

}

} /* namespace es */

#endif /* ES_FLOW_STRESS_H_ */

Includes¶

Eigen/Core
Eigen/Dense
NumericalIntegration.h
iomanip
iostream
math.h
profile.h (File profile.h)
utilities/cumtrapz.h (File cumtrapz.h)
utilities/trapz.h (File trapz.h)

Included By¶

Namespaces¶

Namespace es

Functions¶

File tensors.h¶

↰ Parent directory (source/utilities)

Contents

Definition (source/utilities/tensors.h)
Includes
Included By
Namespaces
Functions

Definition (source/utilities/tensors.h)¶

Program Listing for File tensors.h¶

↰ Return to documentation for file (source/utilities/tensors.h)

/*
 * tensors.h Utilities to help with Eigen::Tensors
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_TENSORS_H
#define ES_FLOW_TENSORS_H

#include <string>


namespace utilities {

template<typename T>
std::string tensor_dims(T &tensor) {
    std::stringstream dims;
    for (auto i = tensor.dimensions().begin(); i != tensor.dimensions().end(); ++i) {
        dims << *i << " x ";
    }
    std::string dim_str = dims.str();
    dim_str.pop_back();
    dim_str.pop_back();
    dim_str.pop_back();
    return dim_str;
}

} /* namespace utilities */

#endif //ES_FLOW_TENSORS_H

Includes¶

string

Included By¶

Namespaces¶

Namespace utilities

Functions¶

Template Function utilities::tensor_dims

File trapz.h¶

↰ Parent directory (source/utilities)

Contents

Definition (source/utilities/trapz.h)
Includes
Included By
Namespaces
Functions

Definition (source/utilities/trapz.h)¶

Program Listing for File trapz.h¶

↰ Return to documentation for file (source/utilities/trapz.h)

/*
 * trapz.h Trapezoidal Numerical Integration for Eigen
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_TRAPZ_H
#define ES_FLOW_TRAPZ_H

#include <Eigen/Dense>
#include <stdexcept>


namespace utilities {


template<typename Derived, typename OtherDerived>
EIGEN_STRONG_INLINE void trapz(Eigen::ArrayBase<OtherDerived> const & out, const Eigen::ArrayBase<Derived>& in)
{
    Eigen::ArrayBase<OtherDerived>& out_ = const_cast< Eigen::ArrayBase<OtherDerived>& >(out);
    Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> inter;

    if (in.rows() == 1) {
        out_.derived().setZero(1, in.cols());
    } else {
        inter = (in.topRows(in.rows()-1) + in.bottomRows(in.rows()-1)) * 0.5;
        out_.derived() = inter.colwise().sum();
    }
}

/*
 * Overload method to return result by value.
 */
template<typename Derived>
EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> trapz(const Eigen::ArrayBase<Derived>& y)
{
    Eigen::Array<typename Eigen::ArrayBase<Derived>::Scalar, Eigen::ArrayBase<Derived>::RowsAtCompileTime, Eigen::ArrayBase<Derived>::ColsAtCompileTime> z;
    trapz(z,y);
    return z;
}

template<typename DerivedX, typename DerivedY, typename DerivedOut>
EIGEN_STRONG_INLINE void trapz(Eigen::ArrayBase<DerivedOut> const & out, const Eigen::ArrayBase<DerivedX>& in_x, const Eigen::ArrayBase<DerivedY>& in_y)
{
    // Input size check
    eigen_assert(in_x.rows() == in_y.rows());
    eigen_assert(in_x.cols() == 1);

    // Get dx for each piece of the integration
    Eigen::Array<typename Eigen::ArrayBase<DerivedX>::Scalar, Eigen::ArrayBase<DerivedX>::RowsAtCompileTime, Eigen::ArrayBase<DerivedX>::ColsAtCompileTime> dx;
    dx = (in_x.bottomRows(in_x.rows()-1) - in_x.topRows(in_x.rows()-1));


    // Get the average heights of the trapezoids
    Eigen::Array<typename Eigen::ArrayBase<DerivedY>::Scalar, Eigen::ArrayBase<DerivedY>::RowsAtCompileTime, Eigen::ArrayBase<DerivedY>::ColsAtCompileTime> inter;
    inter = (in_y.topRows(in_y.rows()-1) + in_y.bottomRows(in_y.rows()-1)) * 0.5;

    // Multiply by trapezoid widths. NB Broadcasting with *= only works for arrayX types, not arrayXX types like dx
    //  inter *= dx;
    for (int i = 0 ; i < inter.cols() ; ++i)
    {
        for (int j = 0 ; j < inter.rows() ; ++j)
        {
            inter(j,i) *= dx(j,0);
        }
    }

    // Initialise output
    Eigen::ArrayBase<DerivedOut>& out_ = const_cast< Eigen::ArrayBase<DerivedOut>& >(out);
    out_.derived().setZero(1, in_y.cols());

    // Output the column-wise sum
    out_.derived() = inter.colwise().sum();
}

/*
 * Overload method to return result by value.
 */
template<typename DerivedX, typename DerivedY, typename DerivedOut>
EIGEN_STRONG_INLINE Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> trapz(const Eigen::ArrayBase<DerivedX>& x, const Eigen::ArrayBase<DerivedY>& y)
{
    Eigen::Array<typename Eigen::ArrayBase<DerivedOut>::Scalar, Eigen::ArrayBase<DerivedOut>::RowsAtCompileTime, Eigen::ArrayBase<DerivedOut>::ColsAtCompileTime> z;
    trapz(z, x, y);
    return z;
}

}  /* namespace utilities */

#endif //ES_FLOW_TRAPZ_H

Includes¶

Eigen/Dense
stdexcept

Included By¶

File stress.h

Namespaces¶

Namespace utilities

Functions¶

File variable_readers.h¶

↰ Parent directory (source)

Contents

Definition (source/variable_readers.h)
Includes
Included By
Namespaces
Functions

Definition (source/variable_readers.h)¶

Program Listing for File variable_readers.h¶

↰ Return to documentation for file (source/variable_readers.h)

/*
 * variable_readers.h Helper functions for reading data from mat files to Eigen arrays and vectors.
 *
 * References:
 *
 *   [1] Eigen: http://eigen.tuxfamily.org/dox/index.html
 *
 *   [2] matio: https://sourceforge.net/projects/matio/
 *
 *   [3] eigen-matio: https://github.com/tesch1/eigen-matio
 *
 * Future Improvements:
 *
 *   [1] Extension to structs and a range of different types
 *
 *   [2] More elegant implementation based on type, or possibly use of eigen-matio (ref [3])
 *
 * Author:                   Tom Clark (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_VARIABLE_READERS_H
#define ES_FLOW_VARIABLE_READERS_H

#include <iostream>
#include <string>
#include "matio.h"
#include <eigen3/Eigen/Core>
#include <unsupported/Eigen/CXX11/Tensor>

using Eigen::Array;
using Eigen::ArrayXXd;
using Eigen::Vector3d;
using Eigen::VectorXd;
using Eigen::Tensor;
using Eigen::Dynamic;


namespace es {


template<size_t SIZE, class T> inline size_t array_size(T (&arr)[SIZE]) {
    return SIZE;
}

void checkVariableType(matvar_t *mat_var, int matvar_type) {

    // Throw error if the requested data type does not match the saved data type
    if (mat_var->class_type != matvar_type) {
        std::string msg = "Error reading mat file (" + std::string(mat_var->name) + " incorrect variable type)";
        throw std::invalid_argument(msg);
    }

}

matvar_t * getVariable(mat_t *matfp, const std::string var_name, bool print_var = true, int max_rank = 2) {

    // Get the variable's structure pointer and check validity
    matvar_t *mat_var = Mat_VarRead(matfp, var_name.c_str());
    if (mat_var == NULL) {
        std::string msg = "Error reading mat file (most likely incorrect file type: " + var_name + " variable is missing)";
        throw std::invalid_argument(msg);
    }

    // Optionally print the variable information to terminal
    if (print_var) {
        std::cout << "Reading variable: " << std::endl;
        Mat_VarPrint(mat_var, true);
    }

    // Check the rank
    if (mat_var->rank > max_rank) {
        std::string msg = "Rank of variable " + var_name + " exceeds required rank of " + std::to_string(max_rank) ;
        throw std::invalid_argument(msg);
    }

    return mat_var;
}

std::string readString(mat_t *matfp, const std::string var_name, bool print_var) {

    // Get the variable's structure pointer and check validity
    matvar_t *mat_var = getVariable(matfp, var_name, print_var);

    // Read the const char from the file, instantiate as std::string
    std::string var = std::string((const char *) mat_var->data, mat_var->dims[1]);

    // Free the data pointer and return the new variable
    Mat_VarFree(mat_var);
    return var;
}

Vector3d readVector3d(mat_t *matfp, const std::string var_name, bool print_var) {

    // Get the variable's structure pointer and check validity
    matvar_t *mat_var = getVariable(matfp, var_name, print_var);

    // Check for three elements always
    if (mat_var->dims[1]*mat_var->dims[0] != 3) {
        std::string msg = "Number of elements in variable '" + var_name + "' not equal to 3";
        throw std::invalid_argument(msg);
    }

    // Read a vector, with automatic transpose to column vector always.
    Vector3d var;
    if (mat_var->dims[0] == 1) {
        if (print_var) std::cout << "Converting row vector '" << mat_var->name << "' to column vector" << std::endl;
    }
    var = Vector3d();
    double *var_d = (double *) mat_var->data;
    long int i = 0;
    for (i = 0; i < 3; i++) {
        var[i] = var_d[i];
    }

    // Free the data pointer and return the new variable
    Mat_VarFree(mat_var);
    return var;

}

VectorXd readVectorXd(mat_t *matfp, const std::string var_name, bool print_var) {

    // Get the variable's structure pointer and check validity
    matvar_t *mat_var = getVariable(matfp, var_name, print_var);

    // Read a vector, with automatic transpose to column vector always.
    VectorXd var;
    if (mat_var->dims[0] == 1) {
        if (print_var) std::cout << "Converting row vector '" << mat_var->name << "' to column vector" << std::endl;
        var = VectorXd(mat_var->dims[1], mat_var->dims[0]);
    } else {
        var = VectorXd(mat_var->dims[0], mat_var->dims[1]);
    }

    // Copy the data into the native Eigen types. However, we can also map to data already in
    // memory which avoids a copy. Read about mapping here:
    // http://eigen.tuxfamily.org/dox/group__TutorialMapClass.html
    // TODO consider mapping to reduce peak memory overhead
    double *var_d = (double *) mat_var->data;
    long int i = 0;
    for (i = 0; i < var.rows()*var.cols(); i++) {
        var[i] = var_d[i];
    }

    // Free the data pointer and return the new variable
    Mat_VarFree(mat_var);
    return var;

}

ArrayXXd readArrayXXd(mat_t *matfp, const std::string var_name, bool print_var) {

    // Get the variable's structure pointer and check validity
    matvar_t *mat_var = getVariable(matfp, var_name, print_var);

    // Read an array. We always read as eigen arrays, not matrices, as these are far more flexible.
    // Can easily cast to matrix type later if linear algebra functionality is needed.

    ArrayXXd var;
    var = ArrayXXd(mat_var->dims[0],  mat_var->dims[1]);

    // Copy the data into the native Eigen types. However, we can also map to data already in
    // memory which avoids a copy. Read about mapping here:
    // http://eigen.tuxfamily.org/dox/group__TutorialMapClass.html
    // TODO consider mapping to reduce peak memory overhead
    double *var_d = (double *) mat_var->data;
    long int i = 0;
    long int j = 0;
    long int ind = 0;
    for (j = 0; j < var.cols(); j++) {
        for (i = 0; i < var.rows(); i++) {
            var(i, j) = var_d[ind];
            ind = ind+1;
        }
    }

    // Free the data pointer and return the new variable
    Mat_VarFree(mat_var);
    return var;

}

Tensor<double, 3> readTensor3d(mat_t *matfp, const std::string var_name, bool print_var) {

    // Get the variable's structure pointer and check validity for rank 3
    matvar_t *mat_var = getVariable(matfp, var_name, print_var, 3);

    // Read the tensor
    // TODO this code typecasts from size_t to long int, meaning possible data loss for very large arrays. Add a check.
    Eigen::Index dim0 = mat_var->dims[0];
    Eigen::Index dim1 = mat_var->dims[1];
    Eigen::Index dim2 = mat_var->dims[2];
    // std::cout  << "Initialising Tensor of size: " << dim0 << ", " << dim1 << ", "  << dim2 << std::endl;
    Eigen::Tensor<double, 3> var = Eigen::Tensor<double, 3>(dim0, dim1, dim2);
    var.setZero();

    // Copy the data into the native Eigen types. However, we can also map to data already in
    // memory which avoids a copy. Read about mapping here:
    // http://eigen.tuxfamily.org/dox/group__TutorialMapClass.html
    // TODO consider mapping to reduce peak memory overhead
    double *var_d = (double *) mat_var->data;
    long int i = 0;
    long int j = 0;
    long int k = 0;
    long int ind = 0;
    for (k = 0; k < var.dimensions()[2]; k++) {
        for (j = 0; j < var.dimensions()[1]; j++) {
            for (i = 0; i < var.dimensions()[0]; i++) {
                var(i, j, k) = var_d[ind];
                ind = ind + 1;
            }
        }
    }

    // Free the data pointer and return the new variable
    Mat_VarFree(mat_var);
    return var;

}

double readDouble(mat_t *matfp, const std::string var_name, bool print_var){

    // Get the variable's structure pointer and check validity
    matvar_t *mat_var = getVariable(matfp, var_name, print_var);

    // Read a single element double
    double *var_d = (double *) mat_var->data;
    double var = var_d[0];

    // Free the data pointer and return the new variable
    Mat_VarFree(mat_var);
    return var;
}


} /* end namespace */


#endif // ES_FLOW_VARIABLE_READERS_H

Includes¶

eigen3/Eigen/Core
iostream
matio.h
string
unsupported/Eigen/CXX11/Tensor

Included By¶

Namespaces¶

Namespace es

Functions¶

File veer.h¶

↰ Parent directory (source/relations)

Contents

Definition (source/relations/veer.h)
Includes
Namespaces
Functions

Definition (source/relations/veer.h)¶

Program Listing for File veer.h¶

↰ Return to documentation for file (source/relations/veer.h)

/*
 * veer.h Atmospheric Boundary Layer Veer relations
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2015-9 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_VEER_H_
#define ES_FLOW_VEER_H_

#include <Eigen/Dense>
#include <Eigen/Core>

#include "definitions.h"
#include "profile.h"


namespace es {


template <typename T>
T veer_lhs(T const & ui, const double ui_g, const double phi){
    T lhs = 2.0*OMEGA_WORLD*sind(phi)*(ui_g - ui);
    return lhs;
}


template <typename T>
T veer_rhs(T & ui, T & uiu3_bar, const double nu){
    T rhs, mixing_term;
    mixing_term = nu*ui.getZDerivative() + uiu3_bar;
    rhs = mixing_term.getZDerivative();
    return rhs;
}

} /* namespace es */

#endif /* ES_FLOW_VEER_H_ */

Includes¶

Eigen/Core
Eigen/Dense
definitions.h (File definitions.h)
profile.h (File profile.h)

Namespaces¶

Namespace es

Functions¶

File velocity.h¶

↰ Parent directory (source/relations)

Contents

Definition (source/relations/velocity.h)
Includes
Included By
Namespaces
Functions

Definition (source/relations/velocity.h)¶

Program Listing for File velocity.h¶

↰ Return to documentation for file (source/relations/velocity.h)

/*
 * velocity.h Atmospheric Boundary Layer velocity profile relations
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2013-9 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_VELOCITY_H_
#define ES_FLOW_VELOCITY_H_

#include <Eigen/Dense>
#include <Eigen/Core>
#include "profile.h"
#include <iostream>


namespace es {

template <typename T, typename Talf>
T power_law_speed(T const & z, const double u_ref, const double z_ref, Talf const & alpha){
    T z_norm = z / z_ref;
    T speed = pow(z_norm, alpha) * u_ref;
    return speed;
};

// Remove template specialisations from doc (causes duplicate) @cond
Eigen::VectorXd power_law_speed(Eigen::VectorXd const & z, const double u_ref, const double z_ref, const double alpha) {
    Eigen::VectorXd z_norm = z.array() / z_ref;
    Eigen::VectorXd speed = pow(z_norm.array(), alpha) * u_ref;
    return speed;
};
// @endcond


template <typename T>
T most_law_speed(T const & z, const double kappa, const double d, const double z0, const double L){
    std::cout << "MOST Law not implemented yet" << std::endl;
    T speed;
    return speed;
}


// Remove template specialisation from doc (causes duplicate) @cond
template <>
Eigen::VectorXd most_law_speed(Eigen::VectorXd const & z, const double kappa, const double d, const double z0, const double L){
    std::cout << "MOST Law not implemented yet" << std::endl;
    Eigen::VectorXd speed;
    return speed;
};
// @endcond


template <typename T_z, typename T_pi_j>
T_z marusic_jones_speed(T_z const & z, T_pi_j const pi_j, const double kappa, const double z_0,
                      const double delta, const double u_inf, const double u_tau){
    T_z eta = (z + z_0) / (delta + z_0);
    T_z eta_cubed = pow(eta, 3.0);
    T_z term1 = log(eta) / kappa;
    T_z term2 = (eta_cubed - 1.0) / (3.0 * kappa);
    T_z term3 = 2.0 * pi_j * (1.0 - pow(eta, 2.0) * 3.0 + eta_cubed * 2.0) / kappa;
    T_z u_deficit = term2 - term1 + term3;
    T_z speed = u_inf - u_deficit * u_tau;
    return speed;
};
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_pi_j>
Eigen::VectorXd marusic_jones_speed(Eigen::VectorXd const & z, T_pi_j const pi_j, const double kappa, const double z_0,
                             const double delta, const double u_inf, const double u_tau){
    Eigen::VectorXd eta = (z.array() + z_0) / (delta + z_0);
    Eigen::VectorXd eta_cubed = eta.array().cube();
    Eigen::VectorXd term1 = eta.array().log() / kappa;
    Eigen::VectorXd term2 = (eta_cubed.array() - 1.0) / (3.0 * kappa);
    Eigen::VectorXd term3 = 2.0 * pi_j * (1.0 - eta.array().square() * 3.0 + eta_cubed.array() * 2.0) / kappa;
    Eigen::VectorXd u_deficit = term2 - term1 + term3;
    Eigen::VectorXd speed = u_inf - u_deficit.array() * u_tau;
    return speed;
};
//@endcond


template <typename T_z, typename T_param>
T_z coles_wake(T_z const &eta, T_param const &pi_coles){
    T_z wc, eta_sqd;
    eta_sqd = pow(eta, 2.0);
    wc = 2.0 * eta_sqd * (3.0 - 2.0 * eta)
        - eta_sqd * (1.0 - eta) * (1.0 - 2.0*eta) / pi_coles;
    return wc;
};
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_param>
Eigen::VectorXd coles_wake(Eigen::VectorXd const &eta, T_param const &pi_coles){
    Eigen::VectorXd wc, eta_sqd;
    eta_sqd = eta.array().pow(2.0);
    wc = 2.0 * eta_sqd.array() * (3.0 - 2.0 * eta.array())
        - eta_sqd.array() * (1.0 - eta.array()) * (1.0 - 2.0*eta.array()) / pi_coles;
    return wc;
};
//@endcond


// Do not document @cond
/* Template for isinf to kill off problems where ceres::Jet is used (autodifferentiation) instead of standard types
 */
bool isinf(double in) {
    return std::isinf(in);
};
template <typename T>
bool isinf(T in) {
    return false;
};
// @endcond


template <typename T_z, typename T_param>
T_z lewkowicz_speed(T_z const & z, T_param const & pi_coles, T_param const & kappa, T_param const & u_inf, T_param const & shear_ratio, T_param const &delta_c) {
    T_z f, speed, eta;
    eta = z / delta_c;
    T_param u_tau = u_inf / shear_ratio;
    T_z term1 = log(eta) / (-1.0*kappa);
    T_z term2 = pi_coles * coles_wake(T_z(1.0), pi_coles) / kappa;
    T_z term3 = pi_coles * coles_wake(eta, pi_coles) / kappa;
    f = term1 + term2 - term3;
    // TODO FIX This only works for doubles, floats - need to template for autodiff
    if (isinf(f)) {
        f = u_inf / u_tau;
    };
    speed = u_inf - f * u_tau;
    return speed;
};
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_param>
Eigen::VectorXd lewkowicz_speed(Eigen::VectorXd const & z, T_param const &pi_coles, T_param const &kappa, T_param const &u_inf, T_param const &shear_ratio, T_param const &delta_c=1.0){
    Eigen::VectorXd f, speed, eta;
    eta = z.array() / delta_c;
    T_param u_tau = u_inf/shear_ratio;
    Eigen::VectorXd term1 = eta.array().log() / (-1.0*kappa);
    double term2 = pi_coles * coles_wake(1.0, pi_coles) / kappa;
    Eigen::VectorXd term3 = pi_coles * coles_wake(eta, pi_coles).array() / kappa;
    f = term1.array() + term2 - term3.array();
    for (int k = 0; k < f.size(); k++) {
        if (std::isinf(f[k])) {
            f(k) = u_inf / u_tau;
        }
    }
    speed = u_inf - f.array() * u_tau;
    return speed;
};
template <typename T_param>
Eigen::ArrayXd lewkowicz_speed(Eigen::ArrayXd const & z, T_param const &pi_coles, T_param const &kappa, T_param const &u_inf, T_param const &shear_ratio, T_param const &delta_c=1.0){
    Eigen::ArrayXd f, speed, eta;
    eta = z / delta_c;
    T_param u_tau = u_inf/shear_ratio;
    Eigen::ArrayXd term1 = eta.log() / (-1.0*kappa);
    double term2 = pi_coles * coles_wake(1.0, pi_coles) / kappa;
    Eigen::ArrayXd term3 = pi_coles * coles_wake(eta, pi_coles) / kappa;
    f = term1 + term2 - term3;
    for (int k = 0; k < f.size(); k++) {
        if (std::isinf(f[k])) {
            f(k) = u_inf / u_tau;
        }
    }
    speed = u_inf - f * u_tau;
    return speed;
};
// @endcond


} /* namespace es */


#endif /* ES_FLOW_VELOCITY_H_ */

Includes¶

Eigen/Core
Eigen/Dense
iostream
profile.h (File profile.h)

Included By¶

Namespaces¶

Namespace es

Functions¶

Footnotes

[1]	Not all of this functionality is implemented yet!