Memory Overflow? std::badalloc

Question

I have a program that solves generally for 1D brownian motion using an Euler's Method. Being a stochastic process, I want to average it over many particles. But I find that as I ramp up the number of particles, it overloads and i get the std::badalloc error, which I understand is a memory error.

Here is my full code

#include <iostream>
#include <vector>
#include <fstream>
#include <cmath>
#include <cstdlib>
#include <limits>
#include <ctime>

using namespace std;

// Box-Muller Method to generate gaussian numbers
double generateGaussianNoise(double mu, double sigma) {
    const double epsilon = std::numeric_limits<double>::min();
    const double tau = 2.0 * 3.14159265358979323846;

    static double z0, z1;
    static bool generate;
    generate = !generate;

    if (!generate) return z1 * sigma + mu;

    double u1, u2;
    do {
        u1 = rand() * (1.0 / RAND_MAX);
        u2 = rand() * (1.0 / RAND_MAX);
    } while (u1 <= epsilon);

    z0 = sqrt(-2.0 * log(u1)) * cos(tau * u2);
    z1 = sqrt(-2.0 * log(u1)) * sin(tau * u2);
    return z0 * sigma + mu;
}

int main() {
    // Initialize Variables
    double gg;  // Gaussian Number Picked  from distribution

    // Integrator
    double t0 = 0;  // Setting the Time Window
    double tf = 10;
    double n = 5000;           // Number of Steps
    double h = (tf - t0) / n;  // Time Step Size

    // Set Constants
    const double pii = atan(1) * 4;  // pi
    const double eta = 1;            // viscous constant
    const double m = 1;              // mass
    const double aa = 1;             // radius
    const double Temp = 30;          // Temperature in Kelvins
    const double KB = 1;             // Boltzmann Constant
    const double alpha = (6 * pii * eta * aa);

    // More Constants
    const double mu = 0;  // Gaussian Mean
    const double sigma = 1;  // Gaussian Std Deviation
    const double ng = n;  // No. of pts to generate for Gauss distribution
    const double npart = 1000;  // No. of Particles

    // Initial Conditions
    double x0 = 0;
    double y0 = 0;
    double t = t0;
    // Vectors
    vector<double> storX;  // Vector that keeps displacement values
    vector<double> storY;  // Vector that keeps velocity values

    vector<double> storT;  // Vector to store time
    vector<double> storeGaussian;  // Vector to store Gaussian numbers generated

    vector<double> holder;  // Placeholder Vector for calculation   operations
    vector<double> mainstore;  // Vector that holds the final value desired

    storT.push_back(t0);

    // Prepares mainstore
    for (int z = 0; z < (n+1); z++) {
        mainstore.push_back(0);
    }

    for (int NN = 0; NN < npart; NN++) {
        holder.clear();
        storX.clear();
        storY.clear();
        storT.clear();
        storT.push_back(0);

        // Prepares holder
        for (int z = 0; z < (n+1); z++) {
            holder.push_back(0);

            storX.push_back(0);

            storY.push_back(0);
        }

        // Gaussian Generator
        srand(time(NULL));
        for (double iiii = 0; iiii < ng; iiii++) {
            gg = generateGaussianNoise(0, 1);  // generateGaussianNoise(mu,sigma)
            storeGaussian.push_back(gg);
        }

        // Solver
        for (int ii = 0; ii < n; ii++) {
            storY[ii + 1] =
                storY[ii] - (alpha / m) * storY[ii] * h +
                (sqrt(2 * alpha * KB * Temp) / m) * sqrt(h) * storeGaussian[ii];
            storX[ii + 1] = storX[ii] + storY[ii] * h;
            holder[ii + 1] =
                pow(storX[ii + 1], 2);  // Finds the displacement squared

            t = t + h;
            storT.push_back(t);
        }

        // Updates the Main Storage
        for (int z = 0; z < storX.size(); z++) {
            mainstore[z] = mainstore[z] + holder[z];
        }
    }

    // Average over the number of particles
    for (int z = 0; z < storX.size(); z++) {
        mainstore[z] = mainstore[z] / (npart);
    }

    // Outputs the data
    ofstream fout("LangevinEulerTest.txt");
    for (int jj = 0; jj < storX.size(); jj++) {
        fout << storT[jj] << '\t' << mainstore[jj] << '\t' << storX[jj] << endl;
    }

    return 0;
}

As you can see, npart is the variable that I change to vary the number of particles. But after each iteration, I do clear my storage vectors like storX,storY... So on paper, the number of particles should not affect memory? I am only just calling the compiler to repeat many more times, and add onto the main storage vector mainstore. I am running my code on a computer with 4GB ram.

Would greatly appreciate it if anyone could point out my errors in logic or suggest improvements.

Edit: Currently the number of particles is set to npart = 1000. So when I try to ramp it up to like npart = 20000 or npart = 50000, it gives me memory errors.

Edit2 I've edited the code to allocate an extra index to each of the storage vectors. But it does not seem to fix the memory overflow

Answer 1

There is an out of bounds exception in the solver part. storY has size n and you access ii+1 where i goes up to n-1. So for your code provided. storY has size 5000. It is allowed to access with indices between 0 and 4999 (including) but you try to access with index 5000. The same for storX, holder and mainstore.

Also, storeGaussian does not get cleared before adding new variables. It grows by n for each npart loop. You access only the first n values of it in the solver part anyway.

Answer 2

Please note, that vector::clear removes all elements from the vector, but does not necessarily change the vector's capacity (i.e. it's storage array), see the documentation.

This won't cause the problem here, because you'll reuse the same array in the next runs, but it's something to be aware when using vectors.