bpls segfaults

[NAThompson]

With

$ bpls -v --version
blps: ADIOS file introspection utility

Build configuration:
ADIOS version: 2.5.0.677 (e4522cbf)
C++ Compiler:  AppleClang 11.0.0.11000033
Target OS:     Darwin-18.7.0
Target Arch:   x86_64

I run this code:

#include <iostream>
#include <string>
#include <cmath>
#include <chrono>
#include <adios2.h>

void display_progress(double progress)
{
    int barWidth = 70;

    std::cout << "[";
    int pos = barWidth * progress;
    for (int i = 0; i < barWidth; ++i) {
        if (i < pos) std::cout << "=";
        else if (i == pos) std::cout << ">";
        else std::cout << " ";
    }
    std::cout << "] "
              << int(progress * 100.0)
              << "%\r";
    std::cout.flush();
}

// This momentum is a conserved quantity of the numerical method.
// Use it to sanity check the solution:
template<typename Real>
Real momentum(std::vector<Real> const & u)
{
    Real p = 0;
    for (int64_t i = 0; i < u.size(); ++i)
    {
        p += u[i];
    }
    return p;
}

template<typename Real>
void KdV(int64_t N, Real dt, Real t_max)
{
    using std::ceil;
    using std::cos;
    using std::abs;
    using std::isnan;
    // Sad! Looking forward to std::numeric::pi;
    const Real pi = 4*std::atan(Real(1));

    if (N <= 0)
    {
        throw std::domain_error("N > 0 is required");
    }
    if (dt > 1)
    {
        throw std::domain_error("time step is too big");
    }
    if (dt <= 0)
    {
        throw std::domain_error("dt > 0 is required");
    }

    Real dx = Real(2)/(N);

    int64_t M = ceil(t_max/dt);
    std::cout << "Solving the KdV equation for dx = " << dx << ", dt = " << dt << " and t_max = " << t_max << "\n";
    adios2::ADIOS adios;
    adios2::IO io = adios.DeclareIO("myio");
    auto u_variable = io.DefineVariable<Real>("u(x, t)", {size_t(N)}, {0}, {size_t(N)}, adios2::ConstantDims);
    io.DefineAttribute<Real>("x0", 0);
    io.DefineAttribute<Real>("dx", dx);
    io.DefineAttribute<std::string>("interpretation", "Equispaced");
    adios2::Engine adios_engine = io.Open("kdv.bp", adios2::Mode::Write);

    std::vector<Real> u0(N);
    for (int64_t i = 0; i < N; ++i)
    {
        u0[i] = cos(pi*i*dx);
    }

    adios_engine.BeginStep();
    adios_engine.Put(u_variable, u0.data(), adios2::Mode::Sync);
    adios_engine.EndStep();
    Real original_momentum = momentum(u0);
    std::vector<Real> u1(N);
    for (int64_t i = 0; i < N; ++i)
    {
        Real cdt = cos(pi*i*dx)*dt;
        u1[i] = cos(pi*(i*dx - cdt));
    }

    // Now the update:
    Real k1 = dt/(3*dx);
    Real delta = 0.022;
    Real k2 = delta*delta*dt/(dx*dx*dx);
    Real t1, t2;
    std::vector<Real> u2(N);
    for (int64_t j = 1; j < M - 1; ++j)
    {
        adios_engine.BeginStep();
        t1 = (u1[1] + u1[0] + u1[N-1])*(u1[1] - u1[N-1]);
        t2 = u1[2]- 2*u1[1] + 2*u1[N-1] - u1[N-2];
        u2[0] = u0[0] - k1*t1 - k2*t2;
        t1 = (u1[2] + u1[1] + u1[0])*(u1[2] - u1[0]);
        t2 = u1[3] - 2*u1[2] + 2*u1[0] - u1[N-1];
        u2[1] = u0[1] - k1*t1 - k2*t2;
        for (int64_t i = 2; i < N-2; ++i)
        {
        t1 = (u1[i+1] + u1[i] + u1[i-1])*(u1[i+1] - u1[i-1]);
        t2 = u1[i+2] - 2*u1[i+1] + 2*u1[i-1] - u1[i-2];
        u2[i] = u0[i] - k1*t1 - k2*t2;
        }
        u2[N-2] = u0[N-2] - k1*(u1[N-1] + u1[N-2] + u1[N-3])*(u1[N-1] - u1[N-3]) - k2*(u1[0] - 2*u1[N-1] + 2*u1[N-3] - u1[N-4]);
        u2[N-1] = u0[N-1] - k1*(u1[0  ] + u1[N-1] + u1[N-2])*(u1[0  ] - u1[N-2]) - k2*(u1[1] - 2*u1[0  ] + 2*u1[N-2] - u1[N-3]);

        Real p = momentum(u2);
        if (abs(p) > sqrt(std::numeric_limits<Real>::epsilon()) || isnan(p))
        {
            std::cout << "Solution diverged at t = " << (j+1)*dt << "\n";
            std::cout << "Momentum = " << p << "\n";
            return;
        }
        if ( (j+1) % 10000*4000 == 0)
        {
            display_progress(double(j+1)/(M-1));
            adios_engine.Put(u_variable, u2.data(), adios2::Mode::Sync);
        }
        u0 = u1;
        u1 = u2;
        adios_engine.EndStep();
    }
    Real p = momentum(u2);
    std::cout << "Final momentum = " << p << "\n";
    adios_engine.Close();
}

int main(int argc, char** argv)
{
    using Real = double;
    int64_t N = 256;
    Real t_max = 5;
    if (argc > 1)
    {
        std::string dx_str = argv[1];
        if (dx_str == "-h")
        {
            std::cout << "Usage: ./KdV.x N t_max, where N is number of spacial gridpoints (dt chosen from dx via stability conditions) and t_max is max simulation time; e.g., ./KdV.x 512 10\n";
            return 0;
        }
        N = std::stoi(dx_str);
    }
    if (argc > 2)
    {
        t_max = Real(std::stod(argv[2]));
    }

    Real dx = Real(1)/N;
    Real dt = 27*dx*dx*dx/4;
    KdV<Real>(N, dt, t_max);
}

Then I run

$ bpls kdv.bp
[1]    88208 segmentation fault  bpls kdv.bp

[pnorbert]

I never expected anyone to call adios BeginStep and EndStep without doing I/O. Apparently BP4 cannot handle empty steps. Is this a bug or unintended use of the API? It should be at line 123:

            display_progress(double(j+1)/(M-1));
            adios_engine.BeginStep();
            adios_engine.Put(u_variable, u2.data(), adios2::Mode::Sync);
            adios_engine.EndStep();

[eisenhauer]

Hmm. That's a case we should have a test for... I'll self-assign that.

[NAThompson]

@pnorbert : Your change fixes the issue, and you are correct that it is not an intended use of the API.

So, depending on your priorities, might it be sensible to throw an exception on an empty step? In any case, I've added the low priority label.

The actual reason I had the step outside the if check was to give the Put more time to operate (I didn't have it in sync mode before the bug report), and this gave me an opportunity to screw up!