lorenz.cpp

// Copyright (C) 2012, 2013 Rhys Ulerich
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

/** @file
 * Generate (t, x, y, z) data from the Lorenz attractor on standard output.
 * Intended to be piped into another program, e.g. arsel.cpp, for analysis.
 */

#include <sys/time.h>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <vector>

#include "optionparser.h"
#include "real.hpp"

#define STRINGIFY_HELPER(x) #x
#define STRINGIFY(x) STRINGIFY_HELPER(x)

// Declarations for argument checking logic based upon optionparser.h examples
struct Arg : public option::Arg
{
    static option::ArgStatus Int      (const option::Option& opt, bool msg);
    static option::ArgStatus IntNonNeg(const option::Option& opt, bool msg);
    static option::ArgStatus Double   (const option::Option& opt, bool msg);
    static option::ArgStatus DoublePos(const option::Option& opt, bool msg);
};

// Command line argument declarations for optionparser.h usage
enum OptionIndex {
    UNKNOWN, BETA, BURN, DT, EVERY, MORE, SCHEME, RHO, SEED, SIGMA,
    TFINAL, INITX, INITY, INITZ, HELP
};
enum SchemeType {
    RK3 = 0, RK4, RK2, RK1
};
const option::Descriptor usage[] = {
    {UNKNOWN, 0, "", "", Arg::None,
     "Usage: lorenz [OPTION]...\n"
     "Output tab-separated (t, x, y, z) data from the Lorenz equations ("
         /* Working precision */ STRINGIFY(REAL) "):\n"
     "  d/dt {x,y,z} = {sigma*(y - x), x*(rho - z) - y, x*y - beta*z}\n"
     "Advance using Runge--Kutta, either classical (Euler & RK4) or\n"
     "Gottlieb and Shu 1998 (doi:10.1090/S0025-5718-98-00913-2) (RK2 & RK3).\n"
     "Individual columns may be extracted by piping to cut(1) utility.\n"
     "\n"
     "Options:" },
    {0,0,"","",Arg::None,0}, // table break
    {BETA,     0, "B", "beta",     Arg::Double,
     "  -B \t--beta=BETA  \t Beta coefficient defaulting to 8/3"              },
    {BURN,     0, "b", "burn",     Arg::Double,
     "  -b \t--burn=BURN  \t \"Burn-in\" for 0 <= t < BURN defaulting to 500" },
    {DT,       0, "d", "dt",       Arg::DoublePos,
     "  -d \t--dt=DT      \t Fixed time step size defaulting to 0.01"         },
    {EVERY,    0, "e", "every",    Arg::IntNonNeg,
     "  -e \t--every=N    \t Output every Nth step defaulting to 1"           },
    {MORE,      0, "m", "more",    Arg::None,
     "  -m \t--more       \t Output more columns (xx, xy, xz, yy, yz, zz)"    },
    {RHO,      0, "R", "rho",      Arg::Double,
     "  -R \t--rho=RHO    \t Rho coefficient defaulting to 28"                },
    {SEED,     0, "s", "seed",     Arg::Double,
     "  -s \t--seed=SEED  \t Random seed defaulting to gettimeofday tv_usec"  },
    {SIGMA,    0, "S", "sigma",    Arg::Double,
     "  -S \t--sigma=SIGMA\t Sigma coefficient defaulting to 10"              },
    {TFINAL,   0, "t", "tfinal",   Arg::DoublePos,
     "  -t \t--tfinal=T   \t Advance time until t >= T defaulting to 3000"    },
    {INITX,    0, "x", "initx",    Arg::Double,
     "  -x \t--initx=POS  \t Initial x at t = 0 defaulting to rand()/RAND_MAX"},
    {INITY,    0, "y", "inity",    Arg::Double,
     "  -y \t--inity=POS  \t Initial y at t = 0 defaulting to rand()/RAND_MAX"},
    {INITZ,    0, "z", "initz",    Arg::Double,
     "  -z \t--initz=POS  \t Initial z at t = 0 defaulting to rand()/RAND_MAX"},
    {SCHEME, RK1, "1", "euler",    Arg::None,
     "  -1 \t--euler      \t Advance with 1st order Forward Euler scheme"     },
    {SCHEME, RK2, "2", "rk2",      Arg::None,
     "  -2 \t--rk2        \t Advance with 2nd order TVD Runge--Kutta scheme"  },
    {SCHEME, RK3, "3", "rk3",      Arg::None,
     "  -3 \t--rk3        \t Advance with 3rd order TVD Runge--Kutta scheme"  },
    {SCHEME, RK4, "4", "rk4",      Arg::None,
     "  -4 \t--rk4        \t Advance with classical 4th order Runge--Kutta"   },
    {HELP,     0, "h", "help",     Arg::None,
     "  -h \t--help       \t Display this help message and immediately exit"  },
    {0,0,0,0,0,0}
};

// Enable strict IEEE behavior for reproducibility
// Inlining should help offset at least some of the performance loss
#if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) > 40305
# pragma GCC push_options
#elif _MSC_VER > 1400
# pragma float_control(push)
# pragma float_control(precise, on)
#else
# warning Pushing precise floating point behavior unimplemented for compiler
#endif

/** Computation of the Lorenz equation right hand side. */
static inline void lorenz(
    const real  beta, const real rho,   const real sigma,
    const real  x,    const real y,     const real z,
          real &dxdt,       real &dydt,       real &dzdt)
{
     dxdt = sigma*(y - x);
     dydt = x*(rho - z) - y;
     dzdt = x*y - beta*z;
}

/** Advance by \c dt using one step of 1st order Forward Euler. */
static void euler(
    const real dt, const real beta, const real rho, const real sigma,
    long double &t, real &x, real &y, real &z)
{
    real dxdt, dydt, dzdt;
    lorenz(beta, rho, sigma, x, y, z, dxdt, dydt, dzdt);
    x += dt*dxdt;
    y += dt*dydt;
    z += dt*dzdt;
    t += dt;
}

/**
 * Advance by \c dt using one step of 2nd order TVD Runge--Kutta.
 * This optimal scheme appears in Proposition 3.1 of Gottlieb and Shu 1998.
 */
static void tvd_rk2(
    const real dt, const real beta, const real rho, const real sigma,
    long double &t, real &x, real &y, real &z)
{
    real dxdt, dydt, dzdt;

    lorenz(beta, rho, sigma, x, y, z, dxdt, dydt, dzdt);
    real u1x = x + dt*dxdt;
    real u1y = y + dt*dydt;
    real u1z = z + dt*dzdt;

    lorenz(beta, rho, sigma, u1x, u1y, u1z, dxdt, dydt, dzdt);
    x = (x + u1x + dt*dxdt)/2;
    y = (y + u1y + dt*dydt)/2;
    z = (z + u1z + dt*dzdt)/2;
    t += dt;
}

/**
 * Advance by \c dt using one step of 3rd order TVD Runge--Kutta.
 * The optimal scheme appears in Proposition 3.2 of Gottlieb and Shu 1998.
 */
static void tvd_rk3(
    const real dt, const real beta, const real rho, const real sigma,
    long double &t, real &x, real &y, real &z)
{
    real dxdt, dydt, dzdt;

    lorenz(beta, rho, sigma, x, y, z, dxdt, dydt, dzdt);
    real u1x = x + dt*dxdt;
    real u1y = y + dt*dydt;
    real u1z = z + dt*dzdt;

    lorenz(beta, rho, sigma, u1x, u1y, u1z, dxdt, dydt, dzdt);
    real u2x = (3*x + u1x + dt*dxdt)/4;
    real u2y = (3*y + u1y + dt*dydt)/4;
    real u2z = (3*z + u1z + dt*dzdt)/4;

    lorenz(beta, rho, sigma, u2x, u2y, u2z, dxdt, dydt, dzdt);
    x = (x + 2*u2x + 2*dt*dxdt)/3;
    y = (y + 2*u2y + 2*dt*dydt)/3;
    z = (z + 2*u2z + 2*dt*dzdt)/3;
    t += dt;
}

/**
 * Advance by \c dt using one step of classical 4th order Runge--Kutta.
 */
static void std_rk4(
    const real dt, const real beta, const real rho, const real sigma,
    long double &t, real &x, real &y, real &z)
{
    real dxdt, dydt, dzdt;

    lorenz(beta, rho, sigma, x, y, z, dxdt, dydt, dzdt);
    real ux = x + dt*dxdt/2;
    real uy = y + dt*dydt/2;
    real uz = z + dt*dzdt/2;

    real dx = dt*dxdt/6;
    real dy = dt*dydt/6;
    real dz = dt*dzdt/6;

    lorenz(beta, rho, sigma, ux, uy, uz, dxdt, dydt, dzdt);
    ux = x + dt*dxdt/2;
    uy = y + dt*dydt/2;
    uz = z + dt*dzdt/2;

    dx += dt*dxdt/3;
    dy += dt*dydt/3;
    dz += dt*dzdt/3;

    lorenz(beta, rho, sigma, ux, uy, uz, dxdt, dydt, dzdt);
    ux = x + dt*dxdt;
    uy = y + dt*dydt;
    uz = z + dt*dzdt;

    dx += dt*dxdt/3;
    dy += dt*dydt/3;
    dz += dt*dzdt/3;

    lorenz(beta, rho, sigma, ux, uy, uz, dxdt, dydt, dzdt);
    x += (dx + dt*dxdt/6);
    y += (dy + dt*dydt/6);
    z += (dz + dt*dzdt/6);
    t += dt;
}

// Disable strict IEEE behavior for remainder of file
#if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) > 40305
# pragma GCC pop_options
#elif _MSC_VER > 1400
# pragma float_control(pop)
#else
# warning Popping precise floating point behavior unimplemented for compiler
#endif

// TODO Simulation time t drifts slowly due to accumulation error

int main(int argc, char *argv[])
{
    using namespace std;

    real        beta   = 8./3;
    real        burn   = 500;
    real        dt     = 0.01;
    long        every  = 1;
    real        rho    = 28;
    real        sigma  = 10;
    long double t      = 0;
    real        tfinal = 3000;
    bool        more   = false;
    SchemeType  scheme;
    real        x;
    real        y;
    real        z;
    {
        option::Stats stats(usage, argc-(argc>0), argv+(argc>0));

        vector<option::Option> opts(stats.options_max + stats.buffer_max);

        option::Parser parse(usage, argc-(argc>0), argv+(argc>0),
                             &opts[0], &opts[stats.options_max]);

        if (parse.error() || opts[UNKNOWN]) {
            for (option::Option* o = opts[UNKNOWN]; o; o = o->next()) {
                cerr << "Unknown option: " << o->name << "\n";
            }
            return EXIT_FAILURE;
        }

        if (opts[HELP]) {
            printUsage(cout, usage);
            return EXIT_SUCCESS;
        }

        // SEED must come before any other randomly-generated option, e.g. INITX
        if (opts[SEED]) {
            srandom((unsigned) strtol(opts[SEED].last()->arg, NULL, 10));
        } else {
            struct timeval tv;
            struct timezone tz;
            gettimeofday(&tv, &tz);
            srand((unsigned) tv.tv_usec);
        }

        // INITX, INITY, INITZ either specified or randomized
        x = opts[INITX] ? strtod(opts[INITX].last()->arg, NULL)
                        : (real)rand()/(real)RAND_MAX;
        y = opts[INITY] ? strtod(opts[INITY].last()->arg, NULL)
                        : (real)rand()/(real)RAND_MAX;
        z = opts[INITZ] ? strtod(opts[INITZ].last()->arg, NULL)
                        : (real)rand()/(real)RAND_MAX;

        // Default timestepping scheme is whichever SchemeType has value zero
        scheme = static_cast<SchemeType>(opts[SCHEME].last()->type());

        // Parse remaining options
        if (opts[BETA  ]) beta   = strtod(opts[BETA  ].last()->arg, NULL);
        if (opts[BURN  ]) burn   = strtod(opts[BURN  ].last()->arg, NULL);
        if (opts[DT    ]) dt     = strtod(opts[DT    ].last()->arg, NULL);
        if (opts[EVERY ]) every  = strtol(opts[EVERY ].last()->arg, NULL, 10);
        if (opts[RHO   ]) rho    = strtod(opts[RHO   ].last()->arg, NULL);
        if (opts[SIGMA ]) sigma  = strtod(opts[SIGMA ].last()->arg, NULL);
        if (opts[TFINAL]) tfinal = strtod(opts[TFINAL].last()->arg, NULL);
        if (opts[MORE  ]) more   = true;

        // Warn whenever burn >= tfinal
        if (burn >= tfinal) {
            cerr << "Warning: burn >= tfinal suppresses output ("
                 << burn << " >= " << tfinal << ") \n";
        }
    }

    // Discard 0 <= t < burn
    switch (scheme) {  // Switch designed to avoid spurious jumps
        default:  cerr << "Sanity failure: unknown scheme at "
                       << __FILE__ << ":" << __LINE__ << '\n';
                  return EXIT_FAILURE;
        case RK1: while (t < burn) euler  (dt, beta, rho, sigma, t, x, y, z);
                  break;
        case RK2: while (t < burn) tvd_rk2(dt, beta, rho, sigma, t, x, y, z);
                  break;
        case RK3: while (t < burn) tvd_rk3(dt, beta, rho, sigma, t, x, y, z);
                  break;
        case RK4: while (t < burn) std_rk4(dt, beta, rho, sigma, t, x, y, z);
                  break;
    }

    // Output (t, x, y, z) during burn < t <= tfinal at periodic intervals
    cout.precision(numeric_limits<real>::digits10 + 2);
    cout << showpos;
    do {

        cout << t << '\t' << x << '\t' << y << '\t' << z;
        if (more) {
            cout << '\t' << x*x << '\t' << x*y << '\t' << x*z
                                << '\t' << y*y << '\t' << y*z
                                               << '\t' << z*z;
        }
        cout << '\n';

        switch (scheme) {
        default:
            cerr << "Sanity failure: unknown scheme at "
                 << __FILE__ << ":" << __LINE__ << '\n';
            return EXIT_FAILURE;
        case RK1:
            for (long i = 0; i < every; ++i) {
                euler  (dt, beta, rho, sigma, t, x, y, z);
            }
            break;
        case RK2:
            for (long i = 0; i < every; ++i) {
                tvd_rk2(dt, beta, rho, sigma, t, x, y, z);
            }
            break;
        case RK3:
            for (long i = 0; i < every; ++i) {
                tvd_rk3(dt, beta, rho, sigma, t, x, y, z);
            }
            break;
        case RK4:
            for (long i = 0; i < every; ++i) {
                std_rk4(dt, beta, rho, sigma, t, x, y, z);
            }
            break;
        }

    } while (t < tfinal);

    return EXIT_SUCCESS;
}

option::ArgStatus Arg::Int(const option::Option& opt, bool msg)
{
    char *p = 0;
    if (opt.arg) {
        long v = strtol(opt.arg, &p, 10); (void) v; // Discard
        if (p != opt.arg && !*p) return option::ARG_OK;
    }
    if (msg) {
        (std::cerr << "Option ").write(opt.name, opt.namelen)
                    << " requires an integer argument\n";
    }
    return option::ARG_ILLEGAL;
}

option::ArgStatus Arg::IntNonNeg(const option::Option& opt, bool msg)
{
    char *p = 0;
    if (opt.arg) {
        double v = strtol(opt.arg, &p, 10);
        if (p != opt.arg && !*p && v >= 0) return option::ARG_OK;
    }
    if (msg) {
        (std::cerr << "Option ").write(opt.name, opt.namelen)
                    << " requires a nonnegative integer argument\n";
    }
    return option::ARG_ILLEGAL;
}

option::ArgStatus Arg::Double(const option::Option& opt, bool msg)
{
    char *p = 0;
    if (opt.arg) {
        double v = strtod(opt.arg, &p); (void) v; // Discard
        if (p != opt.arg && !*p) return option::ARG_OK;
    }
    if (msg) {
        (std::cerr << "Option ").write(opt.name, opt.namelen)
                   << " requires a floating point argument\n";
    }
    return option::ARG_ILLEGAL;
}

option::ArgStatus Arg::DoublePos(const option::Option& opt, bool msg)
{
    char *p = 0;
    if (opt.arg) {
        double v = strtod(opt.arg, &p);
        if (p != opt.arg && !*p && v > 0) return option::ARG_OK;
    }
    if (msg) {
        (std::cerr << "Option ").write(opt.name, opt.namelen)
                   << " requires a positive floating point argument\n";
    }
    return option::ARG_ILLEGAL;
}