diff --git a/OpenMP/Matrix.h b/OpenMP/Matrix.h
index 39c666f..455f4a2 100644
--- a/OpenMP/Matrix.h
+++ b/OpenMP/Matrix.h
@@ -1,141 +1,134 @@
 //
 // Created by Joachim Jacques Koerfer on 06.05.2019.
 //
 
 #ifndef MATRIX_H
 #define MATRIX_H
 
 #include <string>
 #include <fstream>
 #include <iostream>
 #include <vector>
 #include <chrono>
 
 #include "exception.h"
 
 /**
  * This matrix class is just a simple wrapper that tries to be as simple as possible without performances sacrifices.
  * It stores the element in a row-first format
  */
 class Matrix {
 private:
     using iterator_ = double*;
 
     int size_x_;
     int size_y_;
     double *data_;
     bool init_ = false;
 public:
     /**
      * TODO : avoid empty declaration
      */
     Matrix() {};
     /**
      * Constructor
      * @param size_x
      * @param size_y
      * @param init if init=false, no value in the matrix is initialized
      */
     Matrix(int size_x, int size_y, bool init = true) : size_x_(size_x), size_y_(size_y), init_(true) {
         data_ = new double[size_x*size_y];
         if(init) {
             for(int i = 0;i < size_x*size_y;++i) {
                 data_[i] = 0;
             }
         }
     };
     /**
      * Construct Matrix from a binary file
      * @param file_name path to the file
      * @param size_x
      * @param size_y
      * @param init if init=false, no value in the matrix is initialized
      */
     Matrix(const std::string &file_name, int size_x, int size_y) : size_x_(size_x), size_y_(size_y), init_(true) {
         //Test reading speed
         auto startTime = std::chrono::high_resolution_clock::now();
         //Check for file existence
         std::ifstream file(file_name.c_str(), std::ios::in | std::ios::binary | std::ios::ate);
         if(!file.good()) {
             throw FileOpenError(file_name);
         }
 
         std::streampos size;
         size = file.tellg();
 
         char* memblock = new char [size];
         file.seekg(0, std::ios::beg);
         file.read(memblock, size);
         file.close();
 
         auto endTime = std::chrono::high_resolution_clock::now();
         auto speed = std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count();
         //std::cout << (double)size / speed / 1024.  << std::endl;
 
         //Create the matrix, reinterpret the data as doubles
         data_ = (double*)memblock;
     }
     Matrix(const Matrix &copy) {
         size_x_ = copy.size_x_;
         size_y_ = copy.size_y_;
         init_ = true;
         data_ = new double[size_x_*size_y_];
         for(int i = 0;i <size_x_*size_y_;i++) {
             data_[i] = copy.data_[i];
         }
     }
     ~Matrix() {
         if(init_)
             delete[] data_ ;
     };
 
-    void update_row(int y, const double *row, int start_x, int end_x) {
-        std::cout << start_x << std::endl;
-        for(int i = start_x;i < end_x;++i) {
-            //data_[size_x_*y + x] = row[x];
-        }
-    }
-
     void save(const std::string &file_name) {
        std::ios_base::sync_with_stdio(false);
 
         auto file = std::fstream(file_name.c_str(), std::ios::out | std::ios::binary);
         if(!file.good()) {
             throw FileWriteError(file_name);
         }
         double size = size_x_*size_y_*8;
         file.write((char*)&data_[0], size);
         file.close();
     }
 
     /**
      * Gets an element from the matrix by reference
      * This function doesn't verify if the indices are out of bound.
      * @param x index x
      * @param y index y
      * @return reference to the element
      */
     double& operator()(int x, int y) { return data_[size_x_*y + x];};
 
     /**
      * Gets an element in the array of the matrix by reference
      * This function doesn't verify if the indices are out of bound.
      * @param index index which range from 0 to size_x*size_y
      * @return reference to the element
      */
     double& operator[](int index) { return data_[index]; };
 
     /**
      * Iterator function that is there for usage of for(auto d : Matrix) {}
      * @return reference to the beginning of the array of doubles
      */
     iterator_ begin() { return &data_[0]; }
     /**
      * Iterator function that is there for usage of for(auto px : Matrix) {}
      * @return reference to the end of the array of doubles
      */
     iterator_ end() { return &data_[size_x_*size_y_]; };
 };
 
 
 #endif //MATRIX_H
diff --git a/OpenMP/Simulation.cc b/OpenMP/Simulation.cc
index bec5ee6..e1d413b 100644
--- a/OpenMP/Simulation.cc
+++ b/OpenMP/Simulation.cc
@@ -1,227 +1,239 @@
 //
 // Created by Joachim Koerfer on 10.05.2019.
 //
 
 #include <cmath>
 #include <omp.h>
 #include <iostream>
 #include <algorithm>
+#include <random>
 #include <iomanip>
 
 #include "Simulation.h"
 
 
-Simulation::Simulation(std::string path, int nx, int size, double Tend, bool verbose) : nx_(nx), size_(size), Tend_(Tend), delta_nx_((double)size/nx), verbose_(verbose) {
+Simulation::Simulation(std::string path, int nx, int size, double Tend, bool verbose, bool random_values) : nx_(nx), size_(size), Tend_(Tend), delta_nx_((double)size/nx), verbose_(verbose) {
+    int num_threads = 1;
+    #pragma omp parallel
+    {
+        num_threads = omp_get_num_threads();
+    }
+
     std::string file_name = path + "Data_nx" + std::to_string(nx) + "_" + std::to_string(size) + "km_T0.2";
+    if(random_values)
+        file_name = path + "Data_nx8001_" + std::to_string(size) + "km_T0.2";
 
     Clock::time_point startTime = Clock::now();
     H_ = std::make_unique<Matrix>(file_name + "_h.bin", nx, nx);
 
     Hu_ = std::make_unique<Matrix>(file_name + "_hu.bin", nx, nx);
     Hv_ = std::make_unique<Matrix>(file_name + "_hv.bin", nx, nx);
 
     Zdx_ = std::make_unique<Matrix>(file_name + "_Zdx.bin", nx, nx);
     Zdy_ = std::make_unique<Matrix>(file_name + "_Zdy.bin", nx, nx);
     double reading_time = std::chrono::duration_cast<std::chrono::milliseconds>(Clock::now() - startTime).count();
 
-    int num_threads = 1;
-    #pragma omp parallel
-    {
-        num_threads = omp_get_num_threads();
+    if(random_values) {
+        if (verbose_) {
+            std::cout << "Initilization of simulation with random values\n"
+                         "Parameters nx=" << nx_ << ", Tend=" << Tend_ << ", size=" << size_ << "km, #threads.="
+                      << num_threads << std::endl;
+        }
     }
-
-    if(verbose_) {
-        std::cout << "Initilization of simulation\n"
-                     "Parameters nx=" << nx_ << ", Tend=" << Tend_ << ", size=" << size_ << "km, #threads.=" << num_threads << std::endl
-                  << "Reading speed=" << (nx_*(double)nx_*8.)/reading_time/1e6*5. << "GB/s" << std::endl;
+    else {
+        if (verbose_) {
+            std::cout << "Initilization of simulation\n"
+                         "Parameters nx=" << nx_ << ", Tend=" << Tend_ << ", size=" << size_ << "km, #threads.="
+                      << num_threads << std::endl
+                      << "Reading speed=" << (nx_ * (double) nx_ * 8.) / reading_time / 1e6 * 5. << "GB/s" << std::endl;
+        }
     }
     //Temp storage
     Ht_ = std::make_unique<Matrix>(*H_);
     Hut_ = std::make_unique<Matrix>(*Hu_);
     Hvt_ = std::make_unique<Matrix>(*Hv_);
 
     //Begin to measure time
     start_time_ = Clock::now();
 
     //Initialisation of the variables necessary to the simulation
     T_ = 0.;
     niter_ = 0;
 }
 
 void Simulation::compute() {
     while(compute_step());
 }
 bool Simulation::compute_step() {
     if (T_ > Tend_) {
         //End the measurement of time
         end_time_ = Clock::now();
         measured_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time_ - start_time_).count();
         nflops_ = (15 + 2 + 11 + 30 + 30 + 1. )*(double)nx_*(double)nx_*niter_/measured_time*1e3;
         return false;
     }
 
     Clock::time_point startTime = Clock::now();
 
     delta_t_ = compute_delta_t();
     cst_ = delta_t_/(2*delta_nx_);
 
     enforce_boundary_conditions();
 
     update_matrices();
 
     copy_old_borders();
 
     H_.swap(Ht_);
     Hu_.swap(Hut_);
     Hv_.swap(Hvt_);
 
     impose_tolerances();
     niter_++;
 
     if (T_ + delta_t_ > Tend_ && Tend_ != T_)
         delta_t_ = Tend_ - T_;
     T_ += delta_t_;
 
     double measure_time = std::chrono::duration_cast<std::chrono::milliseconds>(Clock::now() - startTime).count();
     double ops = (15 + 2 + 11 + 30 + 30 + 1. )*nx_*nx_;
     if (verbose_)
         std::cout << "Current T = " << T_ << "("  << T_/Tend_*100 << "%) perf. = " << ops/measure_time/1e6 << "Gflops" <<
                   " , time = " << measure_time << "ms" << std::endl;
 
     return true;
 }
 
 
 void Simulation::update_matrices() {
     //Avoid dereferencing the pointers for each x,y
     Matrix &Zdx = *Zdx_;
     Matrix &Zdy = *Zdy_;
 
     Matrix &H = *H_;
     Matrix &Hu = *Hu_;
     Matrix &Hv = *Hv_;
 
     Matrix &Ht = *Ht_;
     Matrix &Hut = *Hut_;
     Matrix &Hvt = *Hvt_;
 
     //Use Lax-Friedrichs method to discretize the problem
     #pragma omp parallel for schedule(static, 256)
     for(int y = nx_ - 2;y >= 1;y--) {
         for (int x = 1; x < nx_ - 1; x++) {
             Ht(x,y) = 0.25 * (H(x, y - 1) + H(x, y + 1) + H(x - 1, y) + H(x + 1, y))
                        + cst_*(Hu(x, y - 1) - Hu(x, y + 1) + Hv(x - 1, y) - Hv(x + 1, y));
             Hut(x, y) = 0.25 * (Hu(x, y - 1) + Hu(x, y + 1) + Hu(x - 1, y) + Hu(x + 1, y))
                         - delta_t_ * g_ * Ht(x, y) * Zdx(x, y)
                         + cst_ * (Hu(x, y - 1) * Hu(x, y - 1) / H(x, y - 1) + .5 * g_ * H(x, y - 1) * H(x, y - 1)
                                   - Hu(x, y + 1) * Hu(x, y + 1) / H(x, y + 1) - .5 * g_ * H(x, y + 1) * H(x, y + 1)
                                   + Hu(x - 1, y) * Hv(x - 1, y) / H(x - 1, y) -
                                   Hu(x + 1, y) * Hv(x + 1, y) / H(x + 1, y));
             Hvt(x,y) = 0.25*(Hv(x,y-1) + Hv(x,y+1) + Hv(x-1,y) + Hv(x+1,y))
                        - delta_t_*g_*Ht(x,y)*Zdy(x,y)
                        + cst_*(Hv(x-1,y)*Hv(x-1,y)/H(x-1,y)+ .5*g_*H(x-1,y)*H(x-1,y)
                                - Hv(x+1,y)*Hv(x+1,y)/H(x+1,y) - .5*g_*H(x+1,y)*H(x+1,y)
                                + Hu(x,y-1)*Hv(x,y-1)/H(x,y-1) - Hu(x,y+1)*Hv(x,y+1)/H(x,y+1));
         }
     }
 }
 
 double Simulation::compute_delta_t() {
     double nu;
     double max_nu = 0;
 
     Matrix &H = *H_;
     Matrix &Hu = *Hu_;
     Matrix &Hv = *Hv_;
 
-
-#pragma omp parallel for schedule(static,256) private(nu) reduction(max:max_nu)
+    #pragma omp parallel for schedule(static,256) private(nu) reduction(max:max_nu)
     for(int y = 0;y < nx_;y++) {
         for (int x = 0; x < nx_; x++) {
             double frac_hu = Hu(x, y) / H(x, y);
             double frac_hv = Hv(x, y) / H(x, y);
             double gh = std::sqrt(H(x, y) * g_);
             double max1 = std::max(std::abs(frac_hu - gh), std::abs(frac_hu + gh));
             double max2 = std::max(std::abs(frac_hv - gh), std::abs(frac_hv + gh));
             nu = std::sqrt(max1 * max1 + max2 * max2);
             if (nu > max_nu) {
                 max_nu = nu;
             }
         }
     }
     return delta_nx_/(std::sqrt(2)*max_nu);
 }
 
 void Simulation::enforce_boundary_conditions() {
     for(int y = 0; y <  nx_; y++){
         //Boundary conditions
         (*H_)(0,y) = (*H_)(1,y);
         (*H_)(nx_ - 1, y) = (*H_)(nx_ - 2, y);
 
         (*Hu_)(0,y) = (*Hu_)(1,y);
         (*Hu_)(nx_ - 1, y) = (*Hu_)(nx_ - 2, y);
 
         (*Hv_)(0,y) = (*Hv_)(1,y);
         (*Hv_)(nx_ - 1, y) = (*Hv_)(nx_ - 2, y);
 
     }
     for(int x = 0; x < nx_; x++){
         //Boundary conditions
         (*H_)(x,0) = (*H_)(x,1);
         (*H_)(x, nx_ - 1) = (*H_)(x, nx_ - 2);
 
         (*Hu_)(x,0) = (*Hu_)(x,1);
         (*Hu_)(x, nx_ - 1) = (*Hu_)(x, nx_ - 2);
 
         (*Hv_)(x,0) = (*Hv_)(x,1);
         (*Hv_)(x, nx_ - 1) = (*Hv_)(x, nx_ - 2);
     }
 }
 void Simulation::copy_old_borders() {
     for(int y = 0; y <  nx_; y++){
         //Old borders have to be copied from temp. storage before swapping pointers
         (*H_)(0,y) = (*Ht_)(0, y);
         (*Hu_)(0,y) = (*Hut_)(0, y);
         (*Hv_)(0,y) = (*Hvt_)(0, y);
         (*H_)(nx_ - 1,y) = (*Ht_)(nx_ - 1, y);
         (*Hu_)(nx_ - 1,y) = (*Hut_)(nx_ - 1, y);
         (*Hv_)(nx_ - 1,y) = (*Hvt_)(nx_ - 1, y);
     }
     for(int x = 0; x < nx_; x++){
         //Old borders have to be copied from temp. storage before swapping pointers
         (*H_)(x,0) = (*Ht_)(x, 0);
         (*Hu_)(x,0) = (*Hut_)(x, 0);
         (*Hv_)(x,0) = (*Hvt_)(x, 0);
         (*H_)(x, nx_ -  1) = (*Ht_)(x, nx_ -  1);
         (*Hu_)(x, nx_ -  1) = (*Hut_)(x, nx_ -  1);
         (*Hv_)(x, nx_ -  1) = (*Hvt_)(x, nx_ -  1);
     }
 }
 
 void Simulation::save(const std::string &file_name) {
     Clock::time_point startTime = Clock::now();
     (*H_).save(file_name);
     double writing_time = std::chrono::duration_cast<std::chrono::milliseconds>(Clock::now() - startTime).count();
 
     if(verbose_) {
         std::cout << "\nSolution saved at " << file_name << "\n"
             "Writing speed=" <<(nx_*(double)nx_*8.)/writing_time/1e6 << "GB/s" << std::endl;
     }
 }
 
 void Simulation::impose_tolerances() {
     Matrix &H = *H_;
     Matrix &Hu = *Hu_;
     Matrix &Hv = *Hv_;
 
     for(int y = 0;y < nx_;y++) {
         for (int x = 0; x < nx_; x++) {
             if(H(x,y) <= 5*1e-4) {
                 Hu(x,y) = 0;
                 Hv(x,y) = 0;
                 if(H(x,y) < 0)
                     H(x,y) = 1e-5;
             }
         }
     }
 }
diff --git a/OpenMP/Simulation.h b/OpenMP/Simulation.h
index 872e800..742530b 100644
--- a/OpenMP/Simulation.h
+++ b/OpenMP/Simulation.h
@@ -1,74 +1,74 @@
 //
 // Created by Joachim Koerfer on 10.05.2019.
 //
 
 #ifndef SIMULATION_H
 #define SIMULATION_H
 
 #include <memory>
 #include <string>
 #include "Matrix.h"
 
 #define SQRT_2 1.414213562373095
 
 struct State {
     int niter = 0;
     double T = 0;
     double max_nu = 0;
     double delta_t = 0;
 
 };
 
 typedef std::chrono::high_resolution_clock Clock;
 
 class Simulation {
 private:
     std::unique_ptr<Matrix> H_;
     std::unique_ptr<Matrix> Hu_;
     std::unique_ptr<Matrix> Hv_;
 
     std::unique_ptr<Matrix> Zdx_;
     std::unique_ptr<Matrix> Zdy_;
 
     //Temp storage
     std::unique_ptr<Matrix> Ht_;
     std::unique_ptr<Matrix> Hut_;
     std::unique_ptr<Matrix> Hvt_;
 
     int size_;
     int nx_;
     double delta_nx_;
     double Tend_;
     const double g_ = 127267.2;
     bool verbose_;
 
     //Variables relative to the state of the simulation
     double T_;
     double delta_t_;
     double cst_;
     int niter_;
 
     //Variables relative to the measure of performance
     Clock::time_point start_time_;
     Clock::time_point end_time_;
     double measured_time = 0;
     double nflops_ = 0;
 
     double compute_delta_t();
     void update_matrices();
     void impose_tolerances();
     void enforce_boundary_conditions();
     void copy_old_borders();
 
 
 public:
-    Simulation(std::string path, int nx, int size, double Tend, bool verbose=true);
+    Simulation(std::string path, int nx, int size, double Tend, bool verbose=true, bool random_values=false);
     void compute();
     bool compute_step();
     double get_flops_performance() { return nflops_; }
     double get_TTS() { return measured_time; }
     void save(const std::string &file_name);
 };
 
 
 #endif //SIMULATION_H
diff --git a/OpenMP/exception.h b/OpenMP/exception.h
index 2483943..8ba2539 100644
--- a/OpenMP/exception.h
+++ b/OpenMP/exception.h
@@ -1,24 +1,35 @@
 //
 // Created by Joachim Koerfer on 08.05.2019.
 //
 #ifndef SEQUENTIAL_EXCEPTION_H
 #define SEQUENTIAL_EXCEPTION_H
 
+#include <string>
 #include <exception>
 
 class FileOpenError: public std::exception
 {
-    virtual const char* what() const throw()
-    {
-        return "Failed to open file (not existing or corrupted)";
+private:
+    std::string msg_;
+public:
+    FileOpenError(std::string file_name) {
+        msg_ = "Failed to open file '" + file_name + "' (not existing or corrupted)";
+    }
+    virtual const char* what() const throw() {
+        return msg_.c_str();
     }
 };
 class FileWriteError: public std::exception
 {
-    virtual const char* what() const throw()
-    {
-        return "Failed to write to file";
+private:
+    std::string msg_;
+public:
+    FileWriteError(std::string file_name) {
+        msg_ = "Failed to write file at '" + file_name + "'";
+    }
+    virtual const char* what() const throw() {
+        return msg_.c_str();
     }
 };
 
 #endif //SEQUENTIAL_EXCEPTION_H
diff --git a/OpenMP/main.cpp b/OpenMP/main.cpp
index b7f769d..343b20a 100644
--- a/OpenMP/main.cpp
+++ b/OpenMP/main.cpp
@@ -1,60 +1,65 @@
 #include <omp.h>
 #include <iostream>
 
 #include "Simulation.h"
 
 #include <algorithm>
 
 int main(int argc, char *argv[]) {
     if (argc < 9) {
         std::cerr << "Usage: " << argv[0] << "-i <data folder> -nx <size of problem> -t <end time of simulation> \n"
-                                             "-o <output folder> \n\n";
+                                             "-o <output folder> [-r : init. with random values] \n\n";
         return 1;
     }
 
     std::string input_folder;
     std::string output_folder;
+    bool random_values = false;
     int nx;
     double tend;
     std::string flag;
 
     //Assigning all the correct arguments
     for(int i = 1;i < argc;++i) {
         if(std::string(argv[i]) == "-o") {
             flag = "o";
             continue;
         }
         else if(std::string(argv[i]) == "-i") {
             flag = "i";
             continue;
         }
         else if(std::string(argv[i]) == "-nx") {
             flag = "nx";
             continue;
         }
         else if(std::string(argv[i]) == "-t") {
             flag = "t";
             continue;
         }
+        else if(std::string(argv[i]) == "-r") {
+            random_values = true;
+            continue;
+        }
 
         if(flag == "o")
             output_folder = std::string(argv[i]);
         else if(flag == "i")
             input_folder = std::string(argv[i]);
         else if(flag == "nx")
             nx = std::stoi(argv[i]);
         else if(flag == "t")
             tend = std::stod(argv[i]);
     }
 
-    Simulation sim(input_folder, nx, 500, tend);
+    Simulation sim(input_folder, nx, 500, tend, true, random_values);
 
     sim.compute();
 
     std::cout << "Perf. = " << sim.get_flops_performance()/1e9 << "Gflops / TTS : " << sim.get_TTS()/1e3 << "s" << std::endl;
 
     std::string file_name = output_folder + "Data_nx" + std::to_string(nx) + "_500km_T0.2_h_Sol.bin";
-    sim.save(file_name);
+    //sim.save(file_name);
 
     return 0;
 }
\ No newline at end of file