The second derivative of poly6 kernel function.

[ZeusYang]

Hello, douyubkim. I just discoverd the following problems. It is about the second derivative of poly6 kernel function. In your code:

inline double SphStdKernel3::secondDerivative(double distance) const {
    if (distance * distance >= h2) {
        return 0.0;
    } else {
        double x = distance * distance / h2;
        return 945.0 / (32.0 * kPiD * h5) * (1 - x) * (3 * x - 1);
    }
}

But I calculated it by hands for serveral times and wrote a different version:

inline double SphStdKernel3::secondDerivative(double distance) const {
    if (distance * distance >= h2) {
        return 0.0;
    } else {
        double x = distance * distance / h2;
        return 945.0 / (32.0 * kPiD * h5) * (1 - x) * (5 * x - 1);
    }
}

It's just 3 should be replaced with 5 in that equation. I know it was missing from any test because we just use the second derivative of spiky kernel function instead of poly6 kernel function.

Meanwhile, there is also a small mistake on Page 133 of your book. It's about formula (2.9), the equation of state. I check out the paper (Becker M, Teschner M. Weakly compressible SPH for free surface flows[C]//Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation. Eurographics Association, 2007: 209-217.) and find that you just put the gamma exponent in the wrong place. The code shown in the book is right.

Thanks.

[doyubkim]

Yes, you are right! Thanks for finding the bug and the typo in the book. For the typos, I have a web page that collects errors from the book, and I will update that page with your findings.

[doyubkim]

PR is up #255

[doyubkim]

The Errata page has been updated: https://fluidenginedevelopment.org/errata/

@yanwc, please let me know if you want to use different name/homepage link for the web page.

[ZeusYang]

@doyubkim, I'm so glad to help you. For the first problem I mentioned above, I think you also need to publish an update in your error collection web page since the code is shown on Page 127.

I got another two typos of your book. On Page 137, the formula (2.14), you forgot the square symbol of m. But the code in the page is right. On Page 134, the formula (2.10), you also forgot the square symbol of c and the code is right.

And I feel a little confused about the EOS scale coefficient of the equation of state. In the paper: http://cg.informatik.uni-freiburg.de/publications/2007_SCA_SPH.pdf The coefficient is $$ \rho_0\frac{c_s^2}{\gamma} $$ But in the book, according the formula (2.9) and the formula (2.10), the EOS scale becomes $$ \rho_0\frac{c_s^2}{\gamma^2} $$ That is, the latter formula is divided by square of gamma . So is this a mistake or innovation?

[doyubkim]

For the EOS scale, the equation is definitely mixed up. My intension was to reference the 2007 SCA paper. Fixing this will affect other parameters in the code to reproduce the similar results, so let me open another issue to track that problem specially.

Thanks for finding all these bugs and typos!!

[ZeusYang]

@doyubkim Yeah, I also think it might be some sort of mixed up. Finally, I made a experiment about that coefficient. I just wanted to compare the difference between them. And so I removed the divisor eosExpoent from the original code in the function below:

inline double computePressureFromEos(
    double density,
    double targetDensity,
    double eosScale,
    double eosExponent,
    double negativePressureScale) {
    // Equation of state
    // (http://www.ifi.uzh.ch/vmml/publications/pcisph/pcisph.pdf)

    // Before
    //double p = eosScale / eosExponent
    //  * (std::pow((density / targetDensity), eosExponent) - 1.0);

    // After
    double p = eosScale
        * (std::pow((density / targetDensity), eosExponent) - 1.0);

    // Negative pressure scaling
    if (p < 0) {
        p *= negativePressureScale;
    }

    return p;
}

Then I used the WCSPH solver to simulate the dam breaking scene, and I got the following two results. I don't think there exists obvious difference. It looks like the solver isn't sensitive to this modification at all. The first GIF shows no modification version while the second one shows the modification version I made above.

[doyubkim]

Thanks for testing the fix with full sim! Yes, it doesn't make any noticeable difference from this particular tests which is a good sign. I think factor of 7 (the default gamma value) didn't make Tait equation stiff enough for given time-steps, which makes sense because the time-step is controlled dynamically.

[doyubkim]

Errata page has been updated to reflect the latest changes.

[ZeusYang]

Thanks for testing the fix with full sim! Yes, it doesn't make any noticeable difference from this particular tests which is a good sign. I think factor of 7 (the default gamma value) didn't make Tait equation stiff enough for given time-steps, which makes sense because the time-step is controlled dynamically.

Yes, it does. I can't agree more. Adaptive time-step controlling plays a significant role. But removing that divisor from the code will make simulation process more time-consuming because of the stiffer EOS scale, won't it? It's my great honour to help this out. Fluid simulation is so funny and I am really grateful to you for your book and code.

[doyubkim]

Thanks! Please feel free to report any other errata or bugs.