fix wrong formulae in Louis' code

rboman · Mar 26, 2024 · cf2f8a6 · cf2f8a6
1 parent ff3a545
commit cf2f8a6
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 17 deletions.
diff --git 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
@@ -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<double>::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
@@ -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<dTf) then
                 dTf = dTftemp
             end if