1D shock hydrodynamics test

euler1d.cpp

/* euler1d.cpp - A 1D shock hydrodynamics example
 * Copyright 2017 Jason Sewall
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
// compile as $(CXX) -std=c++11 -Wall euler1d.cpp -o euler1d

#include <cmath>
#include <cassert>
#include <cstring>
#include <cstdio>
#include <algorithm>

template <int N>
struct q_base
{
    float * asfloat() { return data_;}
    enum { order = N };

    float & operator[](int i) { return data_[i]; }
    const float & operator[](int i) const { return data_[i]; }

    void correct() {}

    float data_[N];

    void print(FILE *o) const
    {
        for(int i = 0; i < N; ++i)
        {
            fprintf(o, "%f ", data_[i]);
        }
    }
};

template <typename Q>
struct riemann_solution
{
    Q right_fluctuation;
    Q left_fluctuation;
    float speed[Q::order];
    float abs_speed[Q::order];
    Q wave[Q::order];
};

struct euler1d_q : q_base<3>
{
    euler1d_q() {}
    euler1d_q(float rho, float rhou, float E) { data_[0] = rho; data_[1] = rhou; data_[2] = E;}
    typedef q_base<3> base;
    enum { order = base::order };
    float & rho()             { return base::data_[0]; }
    const float & rho() const { return base::data_[0]; }
    float & rhou()             { return base::data_[1]; }
    const float & rhou() const { return base::data_[1]; }
    float & E()             { return base::data_[2]; }
    const float & E() const { return base::data_[2]; }

    float u()             { return rhou()/rho(); }
    float u() const { return rhou()/rho(); }

    static const char * const fieldnames[];
};

struct euler_phi
{
    euler_phi(const float & ul, const float & rhol, const float & pl,
              const float & ur, const float & rhor, const float & pr,
              const float & gamma)
        : ul_(ul), rhol_(rhol), pl_(pl), ur_(ur), rhor_(rhor), pr_(pr), gamma_(gamma), beta_((gamma_ + 1.0f)/(gamma_ - 1.0f))
    {}

    float operator()(float p) const
    {
        if(p <= 0.0f)
        {
            printf("Non-positive pressure! (%f)\n", p);
        }

        float phil = phi_left(p);
        float phir = phi_right(p);

        return phil - phir;
    }

    float phi_left(float p) const
    {
        float cl = std::sqrt(gamma_*pl_/rhol_);

        if(p <  pl_)
            return ul_ + 2.0f*cl/(gamma_-1.0f)*(1.0f - std::pow(p/pl_, (gamma_-1.0f)/(2.0f*gamma_)));
        else
            return ul_ + 2.0f*cl/std::sqrt(2.0f*gamma_*(gamma_-1.0f))*((1.0f - p/pl_)/std::sqrt(1.0f + beta_*p/pl_));
    }

    float phi_right(float p) const
    {
        float cr = std::sqrt(gamma_*pr_/rhor_);

        if(p <  pr_)
            return ur_ - 2.0f*cr/(gamma_ - 1.0f)*(1.0f - std::pow(p/pr_, (gamma_ - 1.0f)/(2.0f*gamma_)));
        else
            return ur_ - 2.0f*cr/std::sqrt(2.0f*gamma_*(gamma_-1.0f))*((1.0f - p/pr_)/std::sqrt(1.0f + beta_*p/pr_));
    }

    float rho_star_left(float p) const
    {
        return (1.0f + beta_*p/pl_)/(p/pl_ + beta_)*rhol_;
    }

    float rho_star_right(float p) const
    {
        return (1.0f + beta_*p/pr_)/(p/pr_ + beta_)*rhor_;
    }

    const float & ul_;
    const float & rhol_;
    const float & pl_;
    const float & ur_;
    const float & rhor_;
    const float & pr_;
    const float & gamma_;
    const float beta_;
};

template <class f>
static float secant(const f & func, float x0, float x1, float bottom, float top, float tol, int maxiter)
{
    float xn_1 = x0;
    float xn   = x1;
    float fn_1 = func(xn_1);
    float fn   = func(xn);
    float denom = fn-fn_1;
    for(int i = 0; std::abs(fn) > tol && std::abs(denom) > tol && i < maxiter; ++i)
    {
        float newxn = xn - (xn - xn_1)/denom * fn;
        while(newxn <= bottom)
            newxn = (newxn + xn)*0.5f;
        while(newxn >= top)
            newxn = (newxn + xn)*0.5f;
        xn_1 = xn;
        xn = newxn;
        fn_1 = fn;
        fn = func(xn);
        denom = (fn-fn_1);
    }
    return xn;
}

struct euler1d_equation
{
    static constexpr float gamma = 1.6;
    static constexpr float gammaminus1 = 0.4;

    typedef euler1d_q q;
    typedef riemann_solution<q> rs;

    static float p(const q * state) { return (state->E() - 0.5f*state->rhou()*state->rhou()/state->rho())*gammaminus1; }

    static float c(const q * state)
    {
        return std::sqrt(gamma*p(state)/state->rho());
    }

    static q conservative(float rho, float u, float p)
    {
        return q(rho, rho*u, p/gammaminus1 + 0.5f*rho*u*u);
    }

    static q flux(const q * state)
    {
        float pressure = p(state);
        return q(state->rhou(), state->rhou()*state->rhou()/state->rho() + pressure, (state->E() + pressure)*state->rhou()/state->rho());
    }

    static inline void riemann(rs *__restrict__ soln,
                               const q *__restrict__ qi,
                               const q *__restrict__ qiminus1,
                               int dim)
    {
        assert(qiminus1->rho() > 0.0f);
        assert(qi->rho() > 0.0f);

        float ul = qiminus1->u();
        float ur = qi->u();

        float pl = p(qiminus1);
        float pr = p(qi);

        euler_phi phi_m(ul, qiminus1->rho(), pl,
                        ur, qi->rho(), pr,
                        gamma);

        float pm = secant(phi_m, pl, pr, 0.0f, 100.0f, 1e-6, 100);
        float um = phi_m.phi_left(pm);

        float rhostar_l = phi_m.rho_star_left(pm);
        float rhostar_r = phi_m.rho_star_right(pm);

        q star_left (conservative(rhostar_l, um, pm));
        q star_right(conservative(rhostar_r, um, pm));

        float cl = c(&star_left);
        float cr = c(&star_right);

        soln->speed[0] = um - cl;
        soln->speed[1] = um;
        soln->speed[2] = um + cr;

        for(int i =0; i < q::order; ++i)
            soln->wave[0][i] = star_left[i] - (*qi)[i];

        soln->wave[1][0] = rhostar_r-rhostar_l;
        soln->wave[1][1] = 0.0f;
        soln->wave[1][2] = 0.0f;

        for(int i =0; i < q::order; ++i)
            soln->wave[2][i] = (*qiminus1)[i] - star_right[i];

        q fluxi(flux(qi));
        q fluximinus1(flux(qiminus1));

        q flux_star_l(flux(&star_left));
        q flux_star_r(flux(&star_right));

        if(soln->speed[0] > 0.0f)
        {
            for(int i = 1; i < q::order; ++i)
            {
                soln->left_fluctuation [i] = 0.0f;
                soln->right_fluctuation[i] = fluxi[i] - fluximinus1[i];
            }
        }
        else if(soln->speed[2] < 0.0f)
        {
            for(int i = 1; i < q::order; ++i)
            {
                soln->left_fluctuation [i] = fluximinus1[i] - fluxi[i];
                soln->right_fluctuation[i] = 0.0f;
            }
        }
        else
        {
            for(int i = 1; i < q::order; ++i)
            {
                soln->left_fluctuation [i] = flux_star_l[i] - fluximinus1[i];
                soln->right_fluctuation[i] = fluxi[i]       - flux_star_r[i];
            }
        }

        soln->left_fluctuation [0] = 0.0f;
        soln->right_fluctuation[0] = soln->speed[1]*soln->wave[1][0];

        if(soln->speed[1] < 0.0f)
            std::swap(soln->left_fluctuation[0], soln->right_fluctuation[0]);

        for(int i = 0; i < q::order; ++i)
            soln->abs_speed[i] = std::abs(soln->speed[i]);
    }
};

struct euler1d_approx_equation
{
    static constexpr float gamma = 1.6;
    static constexpr float gammaminus1 = 0.4;

    typedef euler1d_q q;
    typedef riemann_solution<q> rs;

    static float p(const q * state) { return (state->E() - 0.5f*state->rhou()*state->rhou()/state->rho())*gammaminus1; }

    static float c(const q * state)
    {
        return std::sqrt(gamma*p(state)/state->rho());
    }

    static q conservative(float rho, float u, float p)
    {
        return q(rho, rho*u, p/gammaminus1 + 0.5f*rho*u*u);
    }

    static q flux(const q * state)
    {
        float pressure = p(state);
        return q(state->rhou(), state->rhou()*state->rhou()/state->rho() + pressure, (state->E() + pressure)*state->rhou()/state->rho());
    }

    static inline void riemann(rs *__restrict__ soln,
                               const q *__restrict__ qi,
                               const q *__restrict__ qiminus1,
                               int dim)
    {
        assert(qiminus1->rho() > 0.0f);
        assert(qi->rho() > 0.0f);

        float ul = qiminus1->u();
        float ur = qi->u();

        float pl = p(qiminus1);
        float pr = p(qi);

        float rhol = qiminus1->rho();
        float rhor = qi->rho();

        float El = qiminus1->E();
        float Er = qi->E();

        float ubar = (std::sqrt(rhol)*ul + std::sqrt(rhor)*ur)/(std::sqrt(rhol)+std::sqrt(rhor));
        float Hbar = ((El + pl)/std::sqrt(rhol) + (Er + pr)/std::sqrt(rhor))/(std::sqrt(rhol)+std::sqrt(rhor));
        float cbar = std::sqrt(gammaminus1*(Hbar - 0.5f*ubar*ubar));

        soln->speed[0] = ubar - cbar;
        soln->speed[1] = ubar;
        soln->speed[2] = ubar + cbar;

        q r[3];
        r[0] = q(1.0f, soln->speed[0], Hbar - ubar*cbar);
        r[1] = q(1.0f, soln->speed[1], 0.5f*ubar*ubar);
        r[2] = q(1.0f, soln->speed[2], Hbar + ubar*cbar);

        q delta;
        for(int i = 0; i < q::order; ++i)
            delta[i] = (*qi)[i] - (*qiminus1)[i];

        float alpha[3];

        alpha[1] = gammaminus1*((Hbar - ubar*ubar)*delta[0] + ubar*delta[1] - delta[2])/(cbar*cbar);
        alpha[2] = (delta[1] + (cbar - ubar)*delta[0] - cbar*alpha[1])/(2*cbar);
        alpha[0] = delta[0] - alpha[1] - alpha[2];

        for(int i = 0; i < q::order; ++i)
            for(int j = 0; j < q::order; ++j)
                soln->wave[i][j] = alpha[i]*r[i][j];

        memset(&(soln->left_fluctuation[0]), 0, sizeof (float) * q::order);
        memset(&(soln->right_fluctuation[0]), 0, sizeof (float) * q::order);

        for(int i = 0; i < q::order; ++i)
            if(soln->speed[i] < 0.0f)
                for(int j = 0; j < q::order; ++j)
                    soln->left_fluctuation[j] += soln->speed[i]*soln->wave[i][j];
            else
                for(int j = 0; j < q::order; ++j)
                    soln->right_fluctuation[j] += soln->speed[i]*soln->wave[i][j];

        for(int i = 0; i < q::order; ++i)
            soln->abs_speed[i] = std::abs(soln->speed[i]);
    }
};

using roe = euler1d_approx_equation;

using full = euler1d_equation;

template <typename SOLVER>
void test(const typename SOLVER::q *left, const typename SOLVER::q *right)
{
    typename SOLVER::rs res;
    SOLVER::riemann(&res, right, left, 0);

    for(int j = 0; j < 3; ++j)
    {
        printf("wave[%d] = ", j);
        res.wave[j].print(stdout);
        printf("\n");
    }
}

int main()
{
    {
        printf("Given parameters\n");

        const roe::q left(1.0, 0.0, 2.0);
        const roe::q right(1.0, 200.0, 6.0);

        printf("left = "); left.print(stdout); printf("\n");
        printf("right = "); right.print(stdout); printf("\n");

        printf("\nFull solution:\n");
        test<full>(&left, &right);
        printf("\nRoe solution:\n");
        test<roe>(&left, &right);
    }
    printf("\n");
    {
        printf("Sod shock tube test\n");

        const roe::q left(roe::conservative(1.0, 0.0, 1.0));
        const roe::q right(roe::conservative(0.125, 0.0, 0.1));

        printf("left = "); left.print(stdout); printf("\n");
        printf("right = "); right.print(stdout); printf("\n");

        printf("\nFull solution:\n");
        test<full>(&left, &right);
        printf("\nRoe solution:\n");
        test<roe>(&left, &right);
    }
};