Random number generator does not work in c++

Actually I should get 10 different numbers with 10 calls of the function HH_model - I do, but they are always the same numbers - is something wrong with the generation of the random numbers? why do I always get the same rowfogle?

#include<math.h>
#include<iostream>
#include<random>
#include<vector>
#include<algorithm>
#include<fstream>
#include<omp.h>

// parameters
constexpr double v_Rest = -65.0;
constexpr double gNa = 1200.0;
constexpr double gK = 360.0;
constexpr double gL = 3.0;
constexpr double vNa = 115.0;
constexpr double vK = -12.0;
constexpr double vL = 10.6;
constexpr double c = 1.0;
constexpr double knoise = 0.0005;

// stepsize PFs
constexpr int steps = 10;

int store[steps];
const bool root1 = true;

// time constants
constexpr double t_end = 1.0;
constexpr double delay = 0.1;
constexpr double duration = 0.1;
constexpr double dt = 0.0025;
constexpr int t_steps = t_end/dt;

constexpr int runs = 1000;

double doubleRand() {
    thread_local std::mt19937 generator(std::random_device{}());
    std::normal_distribution<double> distribution(0.0, 1.0);
    return distribution(generator);
}

    
    double HH_model(const double I, const double area_factor){
    //some calculations here - interesting line is: 
// const double Inoise = (doubleRand() * knoise * sqrt(gNa * A));
    
int main(){
    for(int i=0; i < steps; i++){
        store[i] = HH_model(0,1);
        std::cout << store[i] << std::endl;}
}

EDIT: I am providing the full code:

it's strange, because every time i get exactly the same ouput... so something in the generation of the random numbers just doesn't work, but i can't explain why... when i compile the code online, i actually always get different random numbers.

#include<math.h>
#include<iostream>
#include<random>
#include<vector>
#include<algorithm>
#include<fstream>
#include<omp.h>

// parameters
constexpr double v_Rest = -65.0;
constexpr double gNa = 1200.0;
constexpr double gK = 360.0;
constexpr double gL = 3.0;
constexpr double vNa = 115.0;
constexpr double vK = -12.0;
constexpr double vL = 10.6;
constexpr double c = 1.0;
constexpr double knoise = 0.0005;

bool print = false;
bool bisection = false;
bool test = true;

// stepsize PFs
constexpr int steps = 3;

int store[steps];
const bool root1 = true;

// time constants
constexpr double t_end = 1.0;
constexpr double delay = 0.1;
constexpr double duration = 0.1;
constexpr double dt = 0.0025;
constexpr int t_steps = t_end/dt;

constexpr int runs = 5000;

double doubleRand() {
    thread_local std::mt19937 generator(std::random_device{}());
    std::normal_distribution<double> distribution(0.0, 1.0);
    return distribution(generator);
}

double alphaM(const double v){ return 12.0 * ((2.5 - 0.1 * (v)) / (exp(2.5 - 0.1 * (v)) - 1.0)); }

double betaM(const double v){ return 12.0 * (4.0 * exp(-(v) / 18.0)); }

double betaH(const double v){ return 12.0 * (1.0 / (exp(3.0 - 0.1 * (v)) + 1.0)); }

double alphaH(const double v){ return 12.0 * (0.07 * exp(-(v) / 20.0)); }

double alphaN(const double v){ return 12.0 * ((1.0 - 0.1 * (v)) / (10.0 * (exp(1.0 - 0.1 * (v)) - 1.0))); }

double betaN(const double v){ return 12.0 * (0.125 * exp(-(v) / 80.0)); }

double HH_model(const double I, const double area_factor){

    const double A = 1.0e-8 * area_factor;
    const double C = c*A;

    const double v0 = 0.0;
    const double m0 = alphaM(v0)/(alphaM(v0)+betaM(v0));
    const double h0 = alphaH(v0)/(alphaH(v0)+betaH(v0));
    const double n0 = alphaN(v0)/(alphaN(v0)+betaN(v0));

    int count = 0;

    for(int j=0; j<runs; j++){

        double vT = v0;
        double mT = m0;
        double hT = h0;
        double nT = n0;

        for(int i=0; i<t_steps; i++){

            double IStim = 0.0;
            if ((delay / dt <= (double)i) && ((double)i <= (delay + duration) / dt))
               IStim = I;

            mT = (mT + dt * alphaM(vT)) / (1.0 + dt * (alphaM(vT) + betaM(vT)));
            hT = (hT + dt * alphaH(vT)) / (1.0 + dt * (alphaH(vT) + betaH(vT)));
            nT = (nT + dt * alphaN(vT)) / (1.0 + dt * (alphaN(vT) + betaN(vT)));

            const double iNa = gNa * pow(mT, 3.0) * hT * (vT - vNa);
            const double iK = gK * pow(nT, 4.0) * (vT - vK);
            const double iL = gL * (vT-vL);
            const double Inoise = (doubleRand() * knoise * sqrt(gNa * A));
            const double IIon = ((iNa + iK + iL) * A) + Inoise;

            vT += ((-IIon + IStim) / C) * dt;

            if(vT > 60.0) {
                count++;
                break;
            }
        }
    }
    return count;
}

double g(double x){

return HH_model(x,10)-runs/2;
}

int main(){

if (test){ std::cout << doubleRand() << std::endl;

}

if (print){

    for(int i=0; i < steps; i++){
        store[i] = HH_model(0,1);
        std::cout << store[i] << std::endl;}
}
if (bisection){

     /* Declaring required variables */
     float x0, x1, x, f0, f1, f, e;
     int step = 1;

     /* Setting precision and writing floating point values in fixed-point notation. */


     /* Inputs */
     up:
     std::cout<<"Enter first guess: ";
     std::cin>>x0;
     std::cout<<"Enter second guess: ";
     std::cin>>x1;
     std::cout<<"Enter tolerable error: ";
     std::cin>>e;

     /* Calculating Functional Value */
     f0 = g(x0);
     f1 = g(x1);

     /* Checking whether given guesses brackets the root or not. */
     if( f0 * f1 > 0.0)
     {
          std::cout<<"Incorrect Initial Guesses."<< std::endl;
          goto up;
     }
   /* Implementing Bisection Method */
     std::cout<< std::endl<<"****************"<< std::endl;
     std::cout<<"Bisection Method"<< std::endl;
     std::cout<<"****************"<< std::endl;
     do
     {
          /* Bisecting Interval */
          x = (x0 + x1)/2;
          f = g(x);

          std::cout<<"Iteration-"<< step<<":\t x = "<<  x<<" and f(x) = "<< g(x)<< std::endl;

          if( f0 * f < 0)
          {
               x1 = x;
          }
          else
          {
               x0 = x;
          }
          step = step + 1;
     }while(fabs(f)>e);

     std::cout<< std::endl<< "Root is: "<<  x<< std::endl;
     std::cout << (HH_model(x,10)/runs)*100 << std::endl;


     return 0;

}


}