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_ */