diff --git a/classes/sph0/louis/README.md b/classes/sph0/louis/README.md index b9c11af0..7f6b6439 100644 --- a/classes/sph0/louis/README.md +++ b/classes/sph0/louis/README.md @@ -78,3 +78,10 @@ cmake . -B build && cmake --build build --config Release && run.py -k 10 --cpp t ``` cmake . -DCMAKE_BUILD_TYPE=Release -B build && make -C build -j 10 && ./run.py tests/small.py -k 10 --cpp ``` + + +## Notes résultat code C++ + +* les particules fixes ont une masse de 0 dans paraview +* la densité des particules fixes est constante et vaut rho0 dans paraview +* max(mu_ab) est nul dans paraview (c'est un bug dans les 2 codes où mu_ab est toujours negatif!) \ No newline at end of file diff --git a/classes/sph0/louis/src_cpp/MobileParticle.cpp b/classes/sph0/louis/src_cpp/MobileParticle.cpp index 386d8c85..ca4bd56f 100644 --- a/classes/sph0/louis/src_cpp/MobileParticle.cpp +++ b/classes/sph0/louis/src_cpp/MobileParticle.cpp @@ -103,12 +103,19 @@ MobileParticle::compute_viscosity(Particle *neigh, double alpha, double beta) if (u_ab.dot(x_ab) < 0.0) { mu_ab = this->h * u_ab.dot(x_ab) / (x_ab.dot(x_ab) + 0.01 * this->h * this->h); + double c_ab = 0.5 * (this->c[RKstep] + neigh->c[RKstep]); double rho_ab = 0.5 * (this->rho[RKstep] + neigh->rho[RKstep]); viscosity = (-alpha * c_ab * mu_ab + beta * mu_ab * mu_ab) / rho_ab; + // TODO: remove negative sign and make mu_ab positive! } // update of max_mu_ab for the calculation of the timestep + // TODO: mu_ab is negative; thus max_mu_ab always 0!! + // This is a bug in the code (in Fortran too). + // we should use the absolute value of mu_ab or change the sign of mu_ab + // Anyway, the viscosity is correctly calculated. + // This only impacts the time step calculation. if ((RKstep == 0) && (mu_ab > this->max_mu_ab)) this->max_mu_ab = mu_ab;