Why does my code fail when calling curand too many times?

Question

The Monte-Carlo code below fails ("process exited with code -2147483645", the output of the price function does not get printed) when I increase numPaths to 1000000, unless I decrease numSteps at the same time. If I remove all references to curand and use fixed numbers instead the code always works so I assume the issue has to do with generating too many random numbers. I have tried calling curand_init with a new seed per thread rather than a new subsequence, but it did not resolve the issue.

Also, I have encountered the same issue when using Thrust rather than curand to generate random numbers.

Can someone help?

#include <math.h>
#include <curand_kernel.h>

#include "cuda_runtime.h"
#include <thrust/device_vector.h>
#include <thrust/reduce.h>
#include <thrust/transform_reduce.h>
#include <thrust/transform.h>
#include <thrust/functional.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/zip_iterator.h>
#include <device_launch_parameters.h>
#include <thrust/execution_policy.h>

using State3D = thrust::tuple<double, double, double>;

template<typename State>
struct AbstractModel : public thrust::binary_function<State, int, State>
{
private:
    double seed = 0;
    int numSteps;
public:
    double dt, sqdt;

    AbstractModel(double T, int numSteps)
        : numSteps(numSteps), dt(T / numSteps), sqdt(sqrt(T / numSteps)) {}

    __device__ virtual void nextStep(State& state, curandState &rand_state) const = 0;

    __device__ State operator() (const State &state, const int &i) const
    {
        curandState rand_state;
        curand_init(seed, i, 0, &rand_state);

        State currentState = state;
        for (unsigned i = 0; i < numSteps; i++)
            nextStep(currentState, rand_state);
        return currentState;
    }
};

struct Model : public AbstractModel<State3D>
{
    double lambda, m, w, theta, eta;

    Model(double T, int numSteps, double lambda, double m, double w, double theta, double eta)
        : AbstractModel(T, numSteps), lambda(lambda), m(m), w(w), theta(theta), eta(eta) {}

    __device__ void nextStep(State3D &state, curandState &rand_state) const
    {
        double level = fmax(state.get<2>(), 0.0),
            sqlevel = sqrt(level), vol = lambda * sqlevel;
        state.get<0>() += vol * sqdt * curand_normal_double(&rand_state);
        state.get<1>() += vol * vol * dt;
        state.get<2>() += theta * (1.0 - level) * dt + eta * sqlevel * sqdt * curand_normal_double(&rand_state);
    }
};

struct ModelInitializer : public thrust::unary_function<int, State3D>
{
    __device__ State3D operator() (const int idx) const
    {
        return thrust::make_tuple<double, double, double>(0.0, 0.0, 0.0);
    }
};

struct ModelPayoffTest : public thrust::unary_function<thrust::tuple<double, double, double>, double>
{
    double operator()(const thrust::tuple<double, double, double> &state) const
    {
        return state.get<2>();
    }
};

template<typename Model, typename Payoff>
double price(unsigned numPaths, Model &model, Payoff &payoff)
{
    thrust::device_vector<State3D> markov_states(numPaths);
    //initialize paths with start value
    thrust::transform(thrust::make_counting_iterator(0), thrust::make_counting_iterator(0) + numPaths, markov_states.begin(), ModelInitializer()/*model*/);

    //diffuse paths
    thrust::transform(markov_states.begin(), markov_states.end(), thrust::make_counting_iterator(0)/*rand_states.begin()*/, markov_states.begin(), model);

    return thrust::transform_reduce(markov_states.begin(), markov_states.end(), payoff, 0.0, thrust::plus<double>()) / numPaths;
}

int main()
{
    double lambda = 0.05, m = 20, w = 0.5, theta = 0.05, eta = 0.1;
    double T = 10;
    int numSteps = 10 * 365;

    Model cheyette(T, numSteps, lambda, m, w, theta, eta);
    ModelPayoffTest payoff;

    int numPaths = 100000;
    std::cout << price<Model, ModelPayoffTest>(numPaths, cheyette, payoff) << "\n";
}