From cf2f8a60f4efe346ed09a8d0b2f967f166a50fd9 Mon Sep 17 00:00:00 2001 From: Romain Boman Date: Tue, 26 Mar 2024 11:05:19 +0100 Subject: [PATCH] fix wrong formulae in Louis' code --- classes/sph0/louis/src_cpp/MobileParticle.cpp | 21 ++++++++----------- classes/sph0/louis/src_cpp/Model.cpp | 4 +++- classes/sph0/louis/src_f/SPH_module.f90 | 11 ++++++---- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/classes/sph0/louis/src_cpp/MobileParticle.cpp b/classes/sph0/louis/src_cpp/MobileParticle.cpp index c919a076..13037023 100644 --- a/classes/sph0/louis/src_cpp/MobileParticle.cpp +++ b/classes/sph0/louis/src_cpp/MobileParticle.cpp @@ -41,7 +41,6 @@ MobileParticle::update_vars() drho_dt += this->m * u_ab.dot(this->vec_gradW[i]); double rho2b = neigh->rho[RKstep] * neigh->rho[RKstep]; du_dt += this->m * (neigh->p[RKstep] / rho2b + this->p[RKstep] / rho2a + pi_ab) * this->vec_gradW_mod[i]; - // TODO: BUG: le 2e terme semble etre faux. Il doit impliquer la densité de la particule et pas celle du voisin } break; } @@ -57,7 +56,6 @@ MobileParticle::update_vars() drho_dt += this->m * u_ab.dot(this->vec_gradW[i]); double rho2b = neigh->rho[RKstep] * neigh->rho[RKstep]; du_dt += this->m * (neigh->p[RKstep] / rho2b + this->p[RKstep] / rho2a + pi_ab) * this->vec_gradW[i]; - // TODO: BUG: le 2e terme semble etre faux. Il doit impliquer la densité de la particule et pas celle du voisin } break; } @@ -104,22 +102,21 @@ MobileParticle::compute_viscosity(Particle *neigh, double alpha, double beta) Eigen::Vector3d x_ab = this->coord[RKstep] - neigh->coord[RKstep]; double mu_ab = 0.0; // this term represents a kind of viscosity - if (u_ab.dot(x_ab) < 0.0) + double uab_xab = u_ab.dot(x_ab); + if (uab_xab < 0.0) { - mu_ab = this->h * u_ab.dot(x_ab) / (x_ab.dot(x_ab) + 0.01 * this->h * this->h); + // mistake in the formula in Louis' thesis: missing negative sign + // (see Monagan-1989 eq 4.11 p8) + // => max(mu_ab) calculated later for the timestep will always be 0! + // note: some papers use a negative mu_ab with an absolute value when finding the max. + mu_ab = - this->h * uab_xab / (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! + viscosity = (alpha * c_ab * mu_ab + beta * mu_ab * mu_ab) / rho_ab; + // negative sign removed here due to sign correction in mu_ab } - // 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; diff --git a/classes/sph0/louis/src_cpp/Model.cpp b/classes/sph0/louis/src_cpp/Model.cpp index ae2cb3aa..87dbe25c 100644 --- a/classes/sph0/louis/src_cpp/Model.cpp +++ b/classes/sph0/louis/src_cpp/Model.cpp @@ -284,12 +284,14 @@ Model::update_dt() timers["update_dt"].start(); // timestep relative to the body forces + // mistake in Louis' thesis: the square root is missing + double sq_g = sqrt(9.81); double dTf = std::numeric_limits::max(); for (int i = this->numFP + 1; i < this->numPart; i++) { Particle const *p = this->particles[i]; - double dt = sqrt(p->h / 9.81); + double dt = p->h / sq_g; if (dt < dTf) dTf = dt; } diff --git a/classes/sph0/louis/src_f/SPH_module.f90 b/classes/sph0/louis/src_f/SPH_module.f90 index 25584b97..daf44c30 100644 --- a/classes/sph0/louis/src_f/SPH_module.f90 +++ b/classes/sph0/louis/src_f/SPH_module.f90 @@ -818,7 +818,7 @@ function artificialViscosity(this, neighObj, alpha, beta) uab_xab = dot_product(u_ab, x_ab) if(uab_xab<0.0d0) then - ! mistake in the formula in Louis'thesis: missing negative sign + ! mistake in the formula in Louis' thesis: missing negative sign ! (see Monagan-1989 eq 4.11 p8) ! => max(mu_ab) calculated later for the timestep will always be 0! ! note: some papers use a negative mu_ab with an absolute value when finding the max. @@ -1040,12 +1040,15 @@ subroutine timeStepUpdate(this) real(DP) :: dTcv, dTcvtemp !< time step relative to the viscous forces and Courrant number integer :: i class(fixed_particle), pointer :: cur_ptr - + real(DP) :: sq_g + + sq_g = sqrt(9.81d0) + ! computes the timestep relative to the body forces - dTf = sqrt( this%part(this%numFP+1)%ptr%h / 9.81d0 ) + dTf = this%part(this%numFP+1)%ptr%h / sq_g do i = this%numFP+2, this%numPart cur_ptr => this%part(i)%ptr - dTftemp = sqrt(cur_ptr%h / 9.81d0) + dTftemp = cur_ptr%h / sq_g if(dTftemp