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 <vector>
#include <iostream>
#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 <unsupported/Eigen/Splines>

#include "cpplot.h"

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


using namespace utilities;
using namespace cpplot;


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 r13a_analytic;

    Eigen::ArrayXXd r13b_analytic;

    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;

    Eigen::Index start_idx;

    // TODO consider whether we can tidy these variables up

    Eigen::ArrayXXd ja_fine;

    Eigen::ArrayXXd jb_fine;

    Eigen::ArrayXd r13a_analytic_fine;

    Eigen::ArrayXd r13b_analytic_fine;

    Eigen::ArrayXd minus_t2wa_fine;

    Eigen::ArrayXd minus_t2wb_fine;

    Eigen::ArrayXd lambda_fine;

    Eigen::ArrayXd eta_fine;

    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) {

    /* On the largest eddy scales:
     * lambda_e is the scale of a given eddy, which contributes to the stresses and spectra in the boundary layer.
     *
     * As far as the signatures are concerned, lambda is mapped for a domain spanning the actual eddy.
     * For a signature, values lambda < 0 are valid; since an eddy exerts an influence a short distance above itself
     * (eddies accelerate flow above them, and decelerate flow below).
     *
     * As far as the deconvolution is concerned, lambda is mapped over the whole boundary layer.
     * For a boundary layer profile, values < 0 are not meaningful: R_ij and S_ij must definitively be 0 outside
     * the boundary layer.
     *
     * Since values of lambda < 0 are theoretically valid for an eddy signature, the signatures as produced by
     * EddySignature().computeSignatures() extend out beyond lambda_e (to z/delta = 1.5, in our implementation).
     * Having an eddy whose scale is of the same size as the boundary layer would therefore create an influence outside
     * the boundary layer, breaking all assumptions about boundary layers ever.
     *
     * For us to maintain the fundamental boundary layer top condition that U = U_inf and u' = 0 for all z/delta >= 1,
     * we assume that eddies can have no influence above their own scale (following Perry & Marusic 1995 pp. 371-2).
     *
     * Clipping the signatures like this is the equivalent of setting the eddy signatures to zero where lambda < 0.
     *
     */
    data.start_idx = 0;
    for (Eigen::Index i=0; i<signature_a.lambda.rows(); i++) {
        if (signature_a.lambda(i) >= 0) {
            data.start_idx = i;
            break;
        }
    }
    Eigen::Index n_lambda_signature = signature_a.lambda.rows() - data.start_idx;
    Eigen::ArrayXd lambda_signature = signature_a.lambda.bottomRows(n_lambda_signature);

    /* On the smallest eddy scales:
     *
     * So the smallest scale we want to go to is lambda_1, which can be defined by the physics of the problem (i.e. the
     * cutoff beyond which viscosity takes over). This scale is defined by the Karman number K_tau.
     *  P&M 1995 again give us an approximate relation: lambda_1 ~= 100 nu/u_tau.
     *
     * As far as this function goes, we take lambda_1 as the smallest scale in the eddy signature array, which assumes
     * that signatures were computed on an appropriate grid.
     *
     */
    double d_lambda_fine = signature_a.domain_spacing(2);
    double lambda_min = lambda_signature(0);
    double lambda_max = lambda_signature(n_lambda_signature-1);
    Eigen::Index n_lambda_fine = round((lambda_max - lambda_min)/d_lambda_fine);

    Eigen::ArrayXd lambda_fine = Eigen::ArrayXd::LinSpaced(n_lambda_fine, lambda_min, lambda_max);

    // Hard limit on lambda max
    double mymax = log(1.0/0.01);
    Eigen::Index up_to = 0;
    for (Eigen::Index i=0; i<lambda_fine.rows(); i++) {
        if (lambda_fine(i) >= mymax) {
            up_to = i;
            break;
        }
    }
    lambda_fine = lambda_fine.topRows(up_to);

    Eigen::ArrayXd eta_fine = lambda_fine.exp().inverse();

    /* On the signatures:
     *
     * To avoid warping the reconstruction, we need to do the convolution and the deconvolution with
     * both the input and the signature on the same, linearly spaced, basis.
     *
     * To do this, we interpolate to a fine distribution and store the fine values.
     *
     */
    Eigen::ArrayXXd ja_fine = signature_a.getJ(lambda_fine);
    Eigen::ArrayXXd jb_fine = signature_b.getJ(lambda_fine);
    Eigen::ArrayXd j13a_fine = ja_fine.col(2);
    Eigen::ArrayXd j13b_fine = jb_fine.col(2);

    // Get an ascending version of eta, containing the point eta=0, for use in determining analytic R13 distributions
    Eigen::ArrayXd bounded_eta = Eigen::ArrayXd(eta_fine.rows()+2);
    bounded_eta.setZero();
    bounded_eta.middleRows(1, eta_fine.rows()) = eta_fine.reverse();
    bounded_eta.bottomRows(1) = 1.0;

    // Get Reynolds Stresses, trim the zero point, reverse back so the ordering is consistent with lambda coordinates
    Eigen::ArrayXd r13a_fine;
    Eigen::ArrayXd r13b_fine;
    reynolds_stress_13(r13a_fine, r13b_fine, data.beta, bounded_eta, data.kappa, data.pi_coles, data.shear_ratio, data.zeta);
    r13a_fine = r13a_fine.middleRows(1, eta_fine.rows());
    r13b_fine = r13b_fine.middleRows(1, eta_fine.rows());
    r13a_fine.reverseInPlace();
    r13b_fine.reverseInPlace();

    // Produce visual check that eta is ascending exponentially
    ScatterPlot pc = ScatterPlot();
    pc.x = Eigen::ArrayXd::LinSpaced(bounded_eta.rows(), 1, bounded_eta.rows());
    pc.y = bounded_eta;
    Layout layc = Layout("Check that eta ascends exponentially");
    layc.xTitle("Row number");
    layc.yTitle("$\\eta$");
    Figure figc = Figure();
    figc.add(pc);
    figc.setLayout(layc);
    figc.write("check_that_eta_ascends_exponentially.json");

    // Produce a visual check that lambda is ascending linearly
    ScatterPlot pc1 = ScatterPlot();
    pc1.x = Eigen::ArrayXd::LinSpaced(lambda_fine.rows(), 1, lambda_fine.rows());
    pc1.y = lambda_fine;
    Layout layc1 = Layout("Check that lambda ascends linearly");
    layc1.xTitle("Row number");
    layc1.yTitle("$\\lambda$");
    Figure figc1 = Figure();
    figc1.add(pc1);
    figc1.setLayout(layc1);
    figc1.write("check_that_lambda_ascends.json");

    // Double check J13 and interpolations
    ScatterPlot pja = ScatterPlot();
    pja.x = lambda_signature;
    pja.y = signature_a.j.bottomRows(n_lambda_signature).col(2);
    pja.name = "Type A (signature)";
    ScatterPlot pjaf = ScatterPlot();
    pjaf.x = lambda_fine;
    pjaf.y = j13a_fine;
    pjaf.name = "Type A (fine)";
    ScatterPlot pjb = ScatterPlot();
    pjb.x = lambda_signature;
    pjb.y = signature_b.j.bottomRows(n_lambda_signature).col(2);
    pjb.name = "Type B (signature)";
    ScatterPlot pjbf = ScatterPlot();
    pjbf.x = lambda_fine;
    pjbf.y = j13b_fine;
    pjbf.name = "Type B (fine)";
    Layout layj = Layout();
    layj.xTitle("$\\lambda$");
    layj.yTitle("$J_{13}$");
    Figure figj = Figure();
    figj.add(pja);
    figj.add(pjaf);
    figj.add(pjb);
    figj.add(pjbf);
    figj.setLayout(layj);
    figj.write("check_j13.json");

    // 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_fine;
    Eigen::ArrayXd minus_t2wb_fine;
//    double stability = 0.005;
//    minus_t2wa_fine = utilities::lowpass_fft_deconv(r13a_fine, j13a_fine, "Type_A", stability);
//    minus_t2wb_fine = utilities::lowpass_fft_deconv(r13b_fine, j13b_fine, "Type_B", stability);

    double alpha = 0.1;
    minus_t2wa_fine = utilities::diagonal_loading_deconv(r13a_fine, j13a_fine, alpha);
    minus_t2wb_fine = utilities::diagonal_loading_deconv(r13b_fine, j13b_fine, alpha);

    // Test the roundtripping convolution without any of the remapping
    Eigen::ArrayXd r13a_fine_recon = conv(minus_t2wa_fine.matrix(), j13a_fine.matrix()).matrix();
    Eigen::ArrayXd r13b_fine_recon = conv(minus_t2wb_fine.matrix(), j13b_fine.matrix()).matrix();

    // Show Reynolds Stresses behaving as part of validation
    ScatterPlot pa = ScatterPlot();
    pa.x = lambda_fine;
    pa.y = -1.0*r13a_fine;
    pa.name = "Type A (input)";

    ScatterPlot pb = ScatterPlot();
    pb.x = lambda_fine;
    pb.y = -1.0*r13b_fine;
    pb.name = "Type B (input)";

    ScatterPlot par = ScatterPlot();
    par.x = lambda_fine;
    par.y = -1.0*r13a_fine_recon;
    par.name = "Type A (reconstructed)";

    ScatterPlot pbr = ScatterPlot();
    pbr.x = lambda_fine;
    pbr.y = -1.0*r13b_fine_recon;
    pbr.name = "Type B (reconstructed)";

    Layout lay = Layout();
    lay.xTitle("$\\lambda$");
    lay.yTitle("$-\\overline{u_1u_3}/U_\\tau^2$");

    Figure fig = Figure();
    fig.add(pa);
    fig.add(pb);
    fig.add(par);
    fig.add(pbr);
    fig.setLayout(lay);
    fig.write("check_r13_fine_reconstruction.json");

    // Create a plot to show the t2w terms
    Figure figt = Figure();
    ScatterPlot pt = ScatterPlot();
    pt.x = lambda_fine;
    pt.y = minus_t2wa_fine;
    pt.name = "$-T_{A}^2(\\lambda - \\lambda_E)\\omega_{A}(\\lambda-\\lambda_E)$";
    figt.add(pt);
    ScatterPlot pt2 = ScatterPlot();
    pt2.x = lambda_fine;
    pt2.y = minus_t2wb_fine;
    pt2.name = "$-T_{B}^2(\\lambda - \\lambda_E)\\omega_{B}(\\lambda-\\lambda_E)$";
    figt.add(pt2);
    Layout layt = Layout();
    layt.xTitle("$\\lambda - \\lambda_E$");
    figt.setLayout(layt);
    figt.write("check_t2w.json");
    std::cout << "Wrote figure check_t2w" << std::endl;


    // 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)
//    double pad_value_a = minus_t2wa(minus_t2wa.rows()-1);
//    double pad_value_b = minus_t2wb(minus_t2wb.rows()-1);
//
//    auto last_n = lambda_e.rows() - minus_t2wa.rows();
//    minus_t2wa.conservativeResize(lambda_e.rows(),1);
//    minus_t2wb.conservativeResize(lambda_e.rows(),1);
//    minus_t2wa.tail(last_n) = pad_value_a;
//    minus_t2wb.tail(last_n) = pad_value_b;

    // Store in the data object
    data.lambda_fine = lambda_fine;
    data.r13a_analytic_fine = r13a_fine;
    data.r13b_analytic_fine = r13b_fine;
    data.minus_t2wa_fine = minus_t2wa_fine;
    data.minus_t2wb_fine = minus_t2wb_fine;
    data.ja_fine = ja_fine;
    data.jb_fine = jb_fine;
    data.eta_fine = eta_fine;
    data.eta = eta_fine;

}


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){

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

    // Column-wise convolution of T2w with J
    for (int k = 0; k < 6; k++) {
        data.reynolds_stress_a.col(k) = conv(data.minus_t2wa_fine.matrix(), data.ja_fine.col(k).matrix());
        data.reynolds_stress_b.col(k) = conv(data.minus_t2wb_fine.matrix(), data.jb_fine.col(k).matrix());
    }

    // Trim the zero-padded ends
    data.reynolds_stress = data.reynolds_stress_a + data.reynolds_stress_b;

    // Check that R13 stacks up against the analytical solution
    Figure fig = Figure();
    Layout lay = Layout();
    ScatterPlot paa = ScatterPlot();
    paa.x = data.eta_fine;
    paa.y = -1.0*data.r13a_analytic_fine;
    paa.name = "Type A - analytic";
    ScatterPlot pba = ScatterPlot();
    pba.x = data.eta_fine;
    pba.y = -1.0*data.r13b_analytic_fine;
    pba.name = "Type B - analytic";
    ScatterPlot par = ScatterPlot();
    par.x = data.eta_fine;
    par.y = -1.0*(data.reynolds_stress_a.col(2));
    par.name = "Type A - reconstructed";
    ScatterPlot pbr = ScatterPlot();
    pbr.x = data.eta_fine;
    pbr.y = -1.0*(data.reynolds_stress_b.col(2));
    pbr.name = "Type B - reconstructed";
    fig.add(paa);
    fig.add(pba);
    fig.add(par);
    fig.add(pbr);
    lay.xTitle("$z/\\delta_{c}$");
    lay.yTitle("$-\\overline{u_1u_3}/U_\\tau^2$");
    fig.setLayout(lay);
    fig.write("check_r13_analytic_and_reconstructed.json");
}


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());

    std::cout << "DIMS " << dims[0] << " " << dims[1] << " " << dims[2] << std::endl;
    std::cout << "PSI_DIMS " << psi_dims[0] << " " << psi_dims[1] << " " << psi_dims[2] << std::endl;

    // 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 wavenumbers
        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).reverse();
            psi_b_vec = conv(t2wb_vec, g_b_vec).reverse();

        }
    }

    // Don't store the sum of the spectra - they're memory hungry and we can always add the two together

    // Premultiply by the u_tau^2 term (see eq. 43)
    // auto u_tau_sqd = pow(data.u_tau, 2.0);
    // data.psi_a = u_tau_sqd * psi_a;
    // data.psi_b = u_tau_sqd * psi_b;
    data.psi_a = psi_a;
    data.psi_b = 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,
              bool compute_spectra=true){

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

    // 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)
    data.z = data.eta.array() * data.delta_c;

    // Add k1z to the data structure
    Eigen::ArrayXd eta = data.eta.array();
    data.k1z = signature_a.k1z(eta);

    // 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
    if (compute_spectra) {
        get_spectra(data, signature_a, signature_b);
    }

    return data;

}

} /* namespace es */

#endif /* SOURCE_ADEM_H_ */