inamuronewtonraphsondynamics.hh

来自「open lattice boltzmann project www.open」· HH 代码 · 共 644 行 · 第 1/2 页

        {
            approxMomentum[iDim] += L::c[knownIndexes[iPop]][iDim] *
                cell[knownIndexes[iPop]];
        }
        for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
        {
            approxMomentum[iDim] += L::c[missingIndexes[iPop]][iDim] *
                computeEquilibrium(missingIndexes[iPop],xi[0],csVel,csVelSqr);
        }
    }
}

template<typename T, template<typename U> class Lattice, typename Dynamics, int direction, int orientation>
T InamuroNewtonRaphsonDynamics<T,Lattice,Dynamics,direction,orientation>::
    computeError(const T &rho, const T u[Lattice<T>::d], const T approxMomentum[Lattice<T>::d])
{
    typedef Lattice<T> L;
    
    T err = T();
    for (int iDim = 0; iDim < L::d; ++iDim)
    {
        err += (rho * u[iDim]-approxMomentum[iDim]) * (rho * u[iDim]-approxMomentum[iDim]);
    }
    return sqrt(err);
}

template<typename T, template<typename U> class Lattice, typename Dynamics, int direction, int orientation>
void InamuroNewtonRaphsonDynamics<T,Lattice,Dynamics,direction,orientation>::computeGradGradError(
        T gradGradError[Lattice<T>::d][Lattice<T>::d],T gradError[Lattice<T>::d],
        const T &rho, const T u[Lattice<T>::d],const T xi[Lattice<T>::d],
        const T approxMomentum[Lattice<T>::d],
        const std::vector<int> missingIndexes)
{
    typedef Lattice<T> L;
    
    T csVel[L::d];
    int counter = 0;
    for (int iDim = 0; iDim < L::d; ++iDim)
    {
        if (direction == iDim)
        {
            ++counter;
            csVel[iDim] = u[iDim];
        }
        else
        {
            csVel[iDim] = u[iDim] + xi[iDim+1-counter];
        }
    }
    T csVelSqr = util::normSqr<T,L::d>(csVel);

    std::vector<T> df[L::d];

    for (int iXi = 0; iXi < L::d; ++iXi)
    {
        for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
        {
            df[iXi].push_back(T());
        }
    }
    
    // all the null terms are allready contained in df (no need for else in the ifs)
    
    for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
    {
        T cu = T();
        for (int iDim = 0; iDim < L::d; ++iDim)
        {
            cu += L::c[missingIndexes[iPop]][iDim] * csVel[iDim];
        }
        df[0][iPop] = L::t[missingIndexes[iPop]]
                      * ((T)1+L::invCs2*cu
                         + 0.5 * L::invCs2 * L::invCs2 * cu * cu
                         - 0.5 * L::invCs2 * csVelSqr);
    }

    counter = 0;
    for (int iDim = 0; iDim < L::d; ++iDim)
    {
        if (direction == iDim)
        {
            ++counter;
        }
        else
        {
            for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
            {
                T temp = T();
                for (int jDim = 0; jDim < L::d; ++jDim)
                {
                    temp += L::c[missingIndexes[iPop]][jDim] * csVel[jDim];
                }

                df[iDim+1-counter][iPop] = xi[0]*L::t[missingIndexes[iPop]] *
                    (L::invCs2*L::c[missingIndexes[iPop]][iDim]
                     + L::invCs2*L::invCs2*L::c[missingIndexes[iPop]][iDim]*temp
                     - L::invCs2*csVel[iDim]);
            }
        }
    }

    std::vector<T> ddf[L::d][L::d];

    for (int iAlpha = 0; iAlpha < L::d; ++iAlpha)
    {
        for (int iBeta = 0; iBeta < L::d; ++iBeta)
        {
            for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
            {
                ddf[iAlpha][iBeta].push_back(T());
            }
        }
    }

    // ddf contains allready all the zero terms

    counter = 0;
    for (int iDim = 0; iDim < L::d; ++iDim)
    {
        if (direction == iDim)
        {
            ++counter;
        }
        else
        {
            for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
            {
                T temp = T();
                for (int jDim = 0; jDim < L::d; ++jDim)
                {
                    temp += L::c[missingIndexes[iPop]][jDim] * csVel[jDim];
                }

                T d_rho_sa = L::t[missingIndexes[iPop]] *
                            (L::invCs2*L::c[missingIndexes[iPop]][iDim]
                             + L::invCs2*L::invCs2*L::c[missingIndexes[iPop]][iDim]*temp
                             - L::invCs2*csVel[iDim]);

                ddf[iDim+1-counter][0][iPop] = d_rho_sa;
                ddf[0][iDim+1-counter][iPop] = d_rho_sa;
            }
        }
    }

    for (int iAlpha = 1; iAlpha < L::d; ++iAlpha)
    {
        for (int iBeta = 1; iBeta < L::d; ++iBeta)
        {
            for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
            {
                ddf[iAlpha][iBeta][iPop] = L::t[missingIndexes[iPop]]*xi[0] *
                    L::invCs2*L::invCs2*L::c[missingIndexes[iPop]][iAlpha]*L::c[missingIndexes[iPop]][iBeta];
                if (iAlpha == iBeta)
                    ddf[iAlpha][iBeta][iPop] -= L::t[missingIndexes[iPop]]*xi[0] * L::invCs2;
            }
        }
    }

    // computation of the vector gradError
    counter = 0;
    for (int iDim = 0; iDim < L::d; ++iDim)
    {
        T du[L::d];
        gradError[iDim] = T();
        for (int jDim = 0; jDim < L::d; ++jDim)
        {
            du[jDim] = T();
            for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
            {
                du[jDim] += L::c[missingIndexes[iPop]][jDim] * df[iDim][iPop];
            }
            gradError[iDim] += (approxMomentum[jDim]- rho * u[jDim]) * du[jDim];
        }
        gradError[iDim] *= (T)2;
    }

    // computation of the matrix gradGradError

    for (int iAlpha = 0; iAlpha < L::d; ++iAlpha)
    {
        for (int iBeta = 0; iBeta < L::d; ++iBeta)
        {
            gradGradError[iAlpha][iBeta] = T();

            T duAlpha[L::d], duBeta[L::d], ddu[L::d];
            for (int iDim = 0; iDim < L::d; ++iDim)
            {
                duAlpha[iDim] = T();
                duBeta[iDim] = T();
                ddu[iDim] = T();
                for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop)
                {
                    duAlpha[iDim] += L::c[missingIndexes[iPop]][iDim]
                        * df[iAlpha][iPop];

                    duBeta[iDim] += L::c[missingIndexes[iPop]][iDim]
                        * df[iBeta][iPop];

                    ddu[iDim] += L::c[missingIndexes[iPop]][iDim]
                        * ddf[iAlpha][iBeta][iPop];
                }
                gradGradError[iAlpha][iBeta] += duAlpha[iDim] * duBeta[iDim]
                    + (approxMomentum[iDim]-rho * u[iDim]) * ddu[iDim];
            }
            gradGradError[iAlpha][iBeta] *= (T)2;
// gradGradError = 2*sum_alpha(du_alpha/dxi_beta*du_alpha/dxi_gamma +
//                  approxMomentum_alpha-u_alpha*d^2u_alpha/(dxi_gamma*dxi_beta))
        }
    }
}

template<typename T, template<typename U> class Lattice, typename Dynamics, int direction, int orientation>
bool InamuroNewtonRaphsonDynamics<T,Lattice,Dynamics,direction,orientation>::
    newtonRaphson(T xi[Lattice<T>::d],
        const T gradError[Lattice<T>::d],
        const T gradGradError[Lattice<T>::d][Lattice<T>::d])
{
    typedef Lattice<T> L;

    T invGradGradError[L::d][L::d];
    bool inversion = invert(gradGradError,invGradGradError);

    for (int iXi = 0; iXi < L::d; ++iXi)
    {
        T correction = T();
        for (int iAlpha = 0; iAlpha < L::d; ++iAlpha)
            correction += invGradGradError[iXi][iAlpha] * gradError[iAlpha];

        xi[iXi] -= correction;
    }

    return inversion;
}

template<typename T, template<typename U> class Lattice, typename Dynamics, int direction, int orientation>
bool InamuroNewtonRaphsonDynamics<T,Lattice,Dynamics,direction,orientation>::
    invert(const T a[2][2],T invA[2][2])
{
    T detA = a[0][0]*a[1][1] - a[0][1]*a[1][0];

    if (fabs(detA) < 1.0e-13)
    {
        std::cout << "error detA too small! = " << detA << "\n";
        for (int iAlpha = 0; iAlpha < 2; ++iAlpha)
        {
            for (int iBeta = 0; iBeta < 2; ++iBeta)
            {
                std::cout << a[iAlpha][iBeta] << "\t";
            }
            std::cout << "\n";
        }
        return false;
    }
    else
    {
        invA[0][0] = a[1][1];
        invA[1][1] = a[0][0];

        invA[0][1] = -a[1][0];
        invA[1][0] = -a[0][1];

        for (int iPop = 0; iPop < 2; ++iPop)
        {
            for (int jPop = 0; jPop < 2; ++jPop)
            {
                invA[iPop][jPop] /= detA;
            }
        }
        return true;
    }
}

template<typename T, template<typename U> class Lattice, typename Dynamics, int direction, int orientation>
bool InamuroNewtonRaphsonDynamics<T,Lattice,Dynamics,direction,orientation>::invert(const T a[3][3],T invA[3][3])
{
    T detA = a[0][0]*a[1][1]*a[2][2] + a[1][0]*a[2][1]*a[0][2] + a[2][0]*a[0][1]*a[1][2]
           - a[0][0]*a[2][1]*a[1][2] - a[2][0]*a[1][1]*a[0][2] - a[1][0]*a[0][1]*a[2][2];


    if (fabs(detA) < 1.0e-13)
    {
        std::cout << "error detA too small! = " << detA << "\n";
        for (int iAlpha = 0; iAlpha < 3; ++iAlpha)
        {
            for (int iBeta = 0; iBeta < 3; ++iBeta)
            {
                std::cout << a[iAlpha][iBeta] << "\t";
            }
            std::cout << "\n";
        }
        return false;
    }
    else
    {
        invA[0][0] = a[1][1]*a[2][2]-a[1][2]*a[2][1];
        invA[0][1] = a[0][2]*a[2][1]-a[0][1]*a[2][2];
        invA[0][2] = a[0][1]*a[1][2]-a[0][2]*a[1][1];

        invA[1][0] = a[1][2]*a[2][0]-a[1][0]*a[2][2];
        invA[1][1] = a[0][0]*a[2][2]-a[0][2]*a[2][0];
        invA[1][2] = a[0][2]*a[1][0]-a[0][0]*a[1][2];

        invA[2][0] = a[1][0]*a[2][1]-a[1][1]*a[2][0];
        invA[2][1] = a[0][1]*a[2][0]-a[0][0]*a[2][1];
        invA[2][2] = a[0][0]*a[1][1]-a[0][1]*a[1][0];

        for (int iPop = 0; iPop < 3; ++iPop)
        {
            for (int jPop = 0; jPop < 3; ++jPop)
            {
                invA[iPop][jPop] /= detA;
            }
        }
        return true;
    }
}

}  // namespace olb


#endif