Merge commit '439bf24' into bug

classes/sph0/louis/src_cpp/CMakeLists.txt

@@ -22,9 +22,9 @@ SET(CMAKE_CXX_STANDARD_REQUIRED ON)
 # increase warning level
 INCLUDE(CheckCXXCompilerFlag)
 CHECK_CXX_COMPILER_FLAG("-Wall" WALL)
-IF(WALL)
-    ADD_COMPILE_OPTIONS(-Wall)
-ENDIF()
+# IF(WALL)
+#     ADD_COMPILE_OPTIONS(-Wall)
+# ENDIF()
 
 # find OpenMP
 FIND_PACKAGE(OpenMP)
@@ -56,5 +56,5 @@ ENDIF()
 
 TARGET_COMPILE_DEFINITIONS(louis++ PUBLIC EIGEN_DONT_PARALLELIZE)
 # try to increase Eigen performance in Debug mode (inlining and no assertions)
-TARGET_COMPILE_DEFINITIONS(louis++ PUBLIC $<$<CONFIG:Debug>: EIGEN_NO_DEBUG >)
-TARGET_COMPILE_DEFINITIONS(louis++ PUBLIC $<$<CONFIG:Debug>: "EIGEN_STRONG_INLINE=__attribute__((always_inline)) inline" > )
+# TARGET_COMPILE_DEFINITIONS(louis++ PUBLIC $<$<CONFIG:Debug>: EIGEN_NO_DEBUG >)
+# TARGET_COMPILE_DEFINITIONS(louis++ PUBLIC $<$<CONFIG:Debug>: "EIGEN_STRONG_INLINE=__attribute__((always_inline)) inline" > )

classes/sph0/louis/src_cpp/DisplayHook.h

@@ -7,6 +7,7 @@ class DisplayHook
 {
 public:
     DisplayHook();
+    virtual ~DisplayHook() = default;
     virtual void display() = 0;
 };
 

classes/sph0/louis/src_cpp/EqState.h

@@ -0,0 +1,77 @@
+#ifndef EQSTATE_H
+#define EQSTATE_H
+
+#include "sph.h"
+
+class EqState
+{
+public:
+    const double rho0; ///< reference density
+protected:
+    const double c0;   ///< speed of sound
+public:
+    EqState(double _rho0, double _c0) : rho0(_rho0), c0(_c0) {}
+    virtual ~EqState() = default;
+    virtual double pressure(double rho) const = 0;
+    virtual double speed_of_sound(double rho) const = 0;
+    virtual double h_factor() const = 0;
+    virtual double dt_factor() const = 0;
+
+};
+
+class IdealGas : public EqState
+{
+    double M;
+
+public:
+    IdealGas(double _rho0, double _c0, double _M)
+        : EqState(_rho0, _c0), M(_M)
+    {
+    }
+    virtual double pressure(double rho) const override
+    {
+        static const double RT = 8.3144621 * 293.15; // J/(mol K) * K
+        return RT / M * (rho / rho0 - 1.0); // eq (3.24)
+    }
+    virtual double speed_of_sound(double rho) const override
+    {
+        return c0;
+    }
+    virtual double h_factor() const override
+    {
+        return 1.1;
+    }
+    virtual double dt_factor() const override
+    {
+        return 5.0;
+    }
+};
+
+class QincFluid : public EqState
+{
+    double gamma;
+
+public:
+    QincFluid(double _rho0, double _c0, double _gamma)
+        : EqState(_rho0, _c0), gamma(_gamma) {}
+
+    virtual double pressure(double rho) const override
+    {
+        double B = c0 * c0 * rho0 / gamma; // eq (3.27)
+        return B * (pow(rho / rho0, gamma) - 1.0);
+    }
+    virtual double speed_of_sound(double rho) const override
+    {
+        return c0 * pow(rho / rho0, (gamma - 1.0) / 2.0);
+    }
+    virtual double h_factor() const override
+    {
+        return 1.02;
+    }
+    virtual double dt_factor() const override
+    {
+        return 1.0;
+    }
+};
+
+#endif // EQSTATE_H

classes/sph0/louis/src_cpp/FixedParticle.cpp

@@ -1,4 +1,5 @@
 #include "FixedParticle.h"
+#include "EqState.h"
 #include "Model.h"
 #include <fstream>
 #include <iostream>
@@ -31,14 +32,26 @@ FixedParticle::update_vars()
         this->rho[2] = this->rho[0] + drho_dt * dt / 2.0; // prepare second step
         this->speed[1] = this->speed[0];
         this->coord[1] = this->coord[0];
+        double p = this->model.eqState->pressure(this->rho[1]);
+        double c = this->model.eqState->speed_of_sound(this->rho[1]);
         this->p[1] = this->compute_pressure(this->rho[1]);
-        this->c[1] = this->compute_speedofsound(this->rho[1]);
+        this->c[1] = this->compute_speedofsound(this->rho[1]); 
+        if(p != this->p[1] || c != this->c[1])
+        {
+            std::cout << "Pressure or speed of sound changed" << std::endl;
+            std::cout << "p: " << p << " c: " << c << std::endl;
+            std::cout << "p1: " << this->p[1] << " c1: " << this->c[1] << std::endl;
+        }
+
+
     }
     else // 2nd RK step
     {
         this->rho[2] = this->rho[2] + drho_dt * dt / 2.0;
         this->speed[2] = this->speed[1];
         this->coord[2] = this->coord[1];
+        // this->p[2] = this->model.eqState->pressure(this->rho[2]);
+        // this->c[2] = this->model.eqState->speed_of_sound(this->rho[2]);
         this->p[2] = this->compute_pressure(this->rho[2]);
         this->c[2] = this->compute_speedofsound(this->rho[2]);
     }

classes/sph0/louis/src_cpp/MobileParticle.cpp

@@ -1,5 +1,6 @@
 #include "MobileParticle.h"
 #include "Model.h"
+#include "EqState.h"
 #include <iostream>
 
 MobileParticle::MobileParticle(Model &m) : FixedParticle(m)
@@ -72,7 +73,9 @@ MobileParticle::update_vars()
         this->coord[1] = this->coord[0] + this->speed[0] * dt;
         this->coord[2] = this->coord[0] + this->speed[0] * dt / 2.0;
         this->p[1] = this->compute_pressure(this->rho[1]);
-        this->c[1] = this->compute_speedofsound(this->rho[1]);
+        this->c[1] = this->compute_speedofsound(this->rho[1]);        
+        // this->p[1] = this->model.eqState->pressure(this->rho[1]);
+        // this->c[1] = this->model.eqState->speed_of_sound(this->rho[1]);
     }
     else // 2nd RK step
     {
@@ -81,6 +84,8 @@ MobileParticle::update_vars()
         this->coord[2] = this->coord[2] + this->speed[1] * dt / 2.0;
         this->p[2] = this->compute_pressure(this->rho[2]);
         this->c[2] = this->compute_speedofsound(this->rho[2]);
+        // this->p[2] = this->model.eqState->pressure(this->rho[2]);
+        // this->c[2] = this->model.eqState->speed_of_sound(this->rho[2]);
     }
 }
 

classes/sph0/louis/src_cpp/Model.cpp

@@ -8,20 +8,25 @@
 #include "MobileParticle.h"
 #include "Sorter.h"
 #include "Kernels.h"
+#include "EqState.h"
+#include "DisplayHook.h"
 
-Model::Model() : sorter(*this)
+Model::Model()
+    : sorter(*this), kernel(nullptr),
+      eqState(nullptr), displayHook(nullptr)
 {
     this->timeStep = 1.0e-15;
     this->currentTime = 0.0;
     this->RKstep = 0;
-    this->kernel = nullptr;
 }
 
 Model::~Model()
 {
     for (int i = 0; i < this->numPart; i++)
         delete this->particles[i];
     delete this->kernel;
+    delete this->eqState;
+    delete this->displayHook;
 }
 
 void
@@ -172,21 +177,30 @@ Model::load_parameters(std::string const &param_path)
     file >> this->numFP;
     file >> this->numMP;
     file >> this->h_0;
-    file >> this->c_0;
-    file >> this->rho_0;
+    double c0, rho0;
+    file >> c0;
+    file >> rho0;
     file >> this->dom_dim;
     int kernelKind;
     file >> kernelKind;
     file >> this->alpha;
     file >> this->beta;
-    file >> this->eqnState;
-    file >> this->state_gamma;
-    file >> this->molMass;
+    int eqnState; 
+    double gamma, molMass;
+    file >> eqnState;
+    file >> gamma;
+    file >> molMass;
     file >> this->kernelCorrection;
     file >> this->maxTime;
     file >> this->saveInt;
     file.close();
 
+    this->c_0 = c0;
+    this->rho_0 = rho0;
+    this->eqnState = eqnState;
+    this->state_gamma = gamma;
+    this->molMass = molMass;
+
     // create kernel
     switch (kernelKind)
     {
@@ -202,7 +216,20 @@ Model::load_parameters(std::string const &param_path)
     default:
         throw std::runtime_error("Bad value of kernel kind");
     }
-    this->kappa = this->kernel->kappa;    
+    this->kappa = this->kernel->kappa;
+
+    // create equation of state
+    switch (eqnState)
+    {
+    case LAW_IDEAL_GAS:
+        this->eqState = new IdealGas(rho0, c0, molMass);
+        break;
+    case LAW_QINC_FLUID:
+        this->eqState = new QincFluid(rho0, c0, gamma);
+        break;
+    default:
+        throw std::runtime_error("Bad value of equation of state");
+    }
 }
 
 /// Save a particle set onto disk.
@@ -213,7 +240,7 @@ Model::load_parameters(std::string const &param_path)
 
 void
 Model::save_particles(std::string const &name, int ite,
-                                int start, int end) const
+                      int start, int end) const
 {
     timers["save"].start();
     // build filename from name and ite
@@ -262,8 +289,10 @@ Model::update_dt()
     this->timeStep = std::min(0.4 * dTf, 0.25 * dTcv);
 
     // possibility to change the timestep if we use the ideal gas law
-    if (this->eqnState == LAW_IDEAL_GAS)
-        this->timeStep = 5 * this->timeStep;
+    this->timeStep *= this->eqState->dt_factor();
+
+    // if (this->eqnState == LAW_IDEAL_GAS)
+    //     this->timeStep = 5 * this->timeStep;
 
     timers["update_dt"].stop();
 }
@@ -283,7 +312,7 @@ Model::update_h()
     mean_rho = mean_rho / this->numPart;
 
     // calculation of the new smoothing length
-    double new_h = this->h_0 * pow(this->rho_0 / mean_rho, 1.0 / 3.0);
+    double new_h = this->h_0 * pow(this->eqState->rho0 / mean_rho, 1.0 / 3.0);
 
     // if the smoothing length is too large, it is limited
     if (new_h > 0.5 * this->sorter.dx)

classes/sph0/louis/src_cpp/Model.h

@@ -13,8 +13,10 @@
 class Model
 {
 public:
-    Sorter sorter; ///< sorter machine
+    Sorter sorter;
     Kernel *kernel;
+    EqState *eqState;
+    DisplayHook *displayHook;
 
     std::vector<Particle *> particles; ///< array of pointers toward particles
 
@@ -25,20 +27,24 @@ class Model
     int kappa;            ///< kappa linked to the eqn state
     double alpha;         ///< weighting factor in the artificial viscosity formulation
     double beta;          ///< weighting factor in the artificial viscosity formulation
+
     int eqnState;         ///< equation of state
                           ///<   1 = ideal gas law
                           ///<   2 = quasi-incompressible fluid
     int state_gamma;      ///< power in eqn State 2.
                           ///< often taken around 7
     double molMass;       ///< Molar mass of the fluid for the prefect gas law
+
     int kernelCorrection; ///< correction of the kernel
                           ///<   0 = no correction
                           ///<   1 = correction enabled
     double maxTime;       ///< simulation time in seconds
     double saveInt;       ///< saving interval
     double h_0;           ///< initial smoothing length
+
     double rho_0;         ///< density of the fluid at free surface
     double c_0;           ///< speed of sound in normal conditions
+
     double timeStep;      ///< timestep (not constant)
     double currentTime;   ///< current time
     int RKstep;           ///< used to know in which RK iteration we are      [RB] (1 or 2)

classes/sph0/louis/src_cpp/Particle.cpp

@@ -1,6 +1,7 @@
 #include "Particle.h"
 #include "Model.h"
 #include "Kernels.h"
+#include "EqState.h"
 #include <fstream>
 #include <iostream>
 
@@ -22,7 +23,9 @@ Particle::load(std::ifstream &ufile, double h_0)
     this->m = m;
     this->h = h_0;
     this->p[0] = this->compute_pressure(rho);
-    this->c[0] = this->compute_speedofsound(rho);
+    this->c[0] = this->compute_speedofsound(rho);   
+    // this->p[0] = this->model.eqState->pressure(rho);
+    // this->c[0] = this->model.eqState->speed_of_sound(rho);
     this->max_mu_ab = 0.0;
 }
 

classes/sph0/louis/src_cpp/Sorter.cpp

@@ -2,6 +2,7 @@
 #include "Model.h"
 #include "FixedParticle.h"
 #include "Neighbour.h"
+#include "EqState.h"
 #include <fstream>
 #include <iostream>
 
@@ -92,7 +93,7 @@ Sorter::compute_hmax()
 {
     double hmax = 0.0;
     for(auto &p : this->model.particles)
-        if (p->h > h_max)
+        if (p->h > hmax)
             hmax = p->h;
 
     // Increase of h_max in order to have a security if h changes.

classes/sph0/louis/src_cpp/Sorter.h

@@ -10,7 +10,7 @@
 class Sorter
 {
     Model &model;
-    double h_max;    ///< maximum smoothing length
+   //double h_max;    ///< maximum smoothing length
 
 public:
     double dx; ///< length of a side of a cube

classes/sph0/louis/src_cpp/sph.h

@@ -11,6 +11,10 @@ class Kernel;
 class CubicSplineKernel;
 class QuadraticKernel;
 class QuinticSplineKernel;
+class DisplayHook;
+class EqState;
+class IdealGas;
+class QIncFluid;
 
 enum KernelKind
 {