Open
Description
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
Activity