lp_solver.h

Go to the documentation of this file.

#ifndef __LP_SOLVER_H__
#define __LP_SOLVER_H__

#include <cstddef>
#include <vector>
#include <string>
#include <fstream>

class Initializer;
class ParticleData;
class NeighbourSearcher;
class EOS;


class LPSolver {

public:
    virtual ~LPSolver() {}

    virtual int solve(double dt) = 0;
    
    virtual double getMinParticleSpacing() const {return m_fMinParticleSpacing;}
    
    virtual double getMaxSoundSpeed() const {return m_fMaxSoundSpeed;}

    virtual double getMaxFluidVelocity() const {return m_fMaxFluidVelocity;}

protected:

    double m_fMinParticleSpacing; 
    double m_fMaxSoundSpeed; 
    double m_fMaxFluidVelocity; 
    bool m_iIfDebug;
    std::ofstream debug;
};


class HyperbolicLPSolver : public LPSolver {

public:
    HyperbolicLPSolver(const Initializer& init, ParticleData* pData, NeighbourSearcher* ns);
    
    virtual int solve(double dt);   

private:

    //-----------------------------------------Data----------------------------------------

    //--------------------------Info got from input argument list---------------------------

    ParticleData* m_pParticleData; 
    NeighbourSearcher* m_pNeighbourSearcher; 
    
    //--------------------------------------------------------------------------------------    

    //--------------------------Info get from Initializer class------------------------------
    
    EOS* m_pEOS; 
    int m_iNumThreads; 
    bool m_iIfMultiThreads; 
    int m_iDimension; 
    bool m_iRandomDirSplitOrder; 
    int m_iLPFOrder; 
    std::size_t m_iNumRow2ndOrder; 
    std::size_t m_iNumRow1stOrder; 
    std::size_t m_iNumCol2ndOrder; 
    std::size_t m_iNumCol1stOrder; 
    bool m_iMovingBoxForGhostParticle; 
    double m_fInitParticleSpacing; 
    double m_fGravity; 
    double m_fInvalidPressure; 
    double m_fInvalidVolume; 
    double m_fNeiSearchRadius; 
    size_t m_iNumParticleWithinSearchRadius; 
    double m_fContactLength; 
    
    //--------------------------------------------------------------------------------------
    
    //---------------------------------Other parameters-------------------------------------

    int m_iNumPhase; 

    std::vector<std::vector<int> > m_vDirSplitTable; 

    int m_iDirSplitOrder;
    
    double m_fDt; 
    
    int m_iContactAlert; 
    
    //std::ofstream debug;///< output information for debugging
    //bool m_iIfDebug;///< if true then print debug info
    //-------------------------------------------------------------------------------------




    //-------------------------------------Methods-----------------------------------------
    
    void computeSetupsForNextIteration(); 
    
        
    void updateFluidBoundingBox();
    
    void checkForContactAlert();
    
    void generateGhostParticle();
    
    bool isValidGhostParticle(double x, double y, double z, int* neiList, size_t numNei, int objectTag); 
    
    void searchNeighbourForAllParticle();   
    
    void setUpwindNeighbourList();
    
    void resetLPFOrder();
    
    void computeMinParticleSpacing();
    
    void computeMaxSoundSpeed();
    
    void computeMaxFluidVelocity();
    
    
    void changeFluidObjectTagIfInContact(int index, size_t numNeiFound, const double* neiListDist);
    
    void changeNeighbourhoodToIncludeOnlyFluidNeiFromSameObject(int index, size_t numNeiFound);
    
    void changeNeighbourhoodToIncludeOnlyFluidNei(int index, size_t numNeiFound);

    
    bool directionalSplitting(int phase);
    
    void setNeighbourListPointers(int dir, // input
        const int **neighbourList0, const int **neighbourList1, // output
        const int **neighbourListSize0, const int **neighbourListSize1);
    
    void setInAndOutDataPointers(int phase, int dir,
        const double** inVelocity, const double** inPressure, const double** inVolume, const double** inSoundSpeed, 
        double** outVelocity, double** outPressure, double** outVolume, double** outSoundSpeed);
        
    void setLPFOrderPointers(int dir, // input
        int** LPFOrder0, int** LPFOrder1, std::vector<int*>& LPFOrderOther); // output

    
    void computeSpatialDer(int dir, size_t index, // input
        int offset, void (HyperbolicLPSolver::*computeA) (size_t, const int *, const int*, size_t, size_t,double*),
        const double* inPressure, const double* inVelocity,
        const int *neighbourList, const int *neighbourListSize, 
        int* LPFOrder, double* vel_d, double* vel_dd, double* p_d, double* p_dd); // output
    
    
    void timeIntegration(double real_dt, double multiplier1st, double multiplier2nd, 
        double gravity, double inVolume, double inVelocity, double inPressure,
        double inSoundSpeed, 
        double vel_d_0, double vel_dd_0, double p_d_0, double p_dd_0,
        double vel_d_1, double vel_dd_1, double p_d_1, double p_dd_1,
        double* outVolume, double* outVelocity, double* outPressure); // output
    
    
    void printInvalidState(int phase, int dir,  int index, double positionX, double positionY, double positionZ,
        double vel_d_0, double vel_dd_0, double p_d_0, double p_dd_0, 
        double vel_d_1, double vel_dd_1, double p_d_1, double p_dd_1);
    
    
    
    bool lowerLPFOrder(int index, const std::vector<int*>& LPFOrderOther, // input
        int* LPFOrder0, int* LPFOrder1); // output  
    
    void computeNumRowAndNumColAndLPFOrder(size_t index, // input
        const int *neighbourList, const int *neighbourListSize, size_t numRow2nd, size_t numRow1st,
        int* LPFOrder, size_t *numRow, size_t *numCol); // output


    void computeA2D(size_t index, const int *neighbourList, const int* LPFOrder, size_t numRow, size_t numCol, // input
                    double *A); // output 
    
    void computeA3D(size_t index, const int *neighbourList, const int* LPFOrder, size_t numRow, size_t numCol, // input
                    double *A); // output
    
    void computeB(size_t index, const int *neighbourList, size_t numRow, const double* inData, 
                  double *b); // output
    
    
    void setBoundaryPressureAndVelocity(int phase);
    
    void setGhostPressureAndVelocity(int phase); 
    
    void computeOthOrderWeightedLPF(std::vector<const double*>& position, 
                                    std::size_t startIndex, std::size_t numParticle,
                                    std::vector<double*>& data);
    
    void updateFluidState();
    
    void moveFluidParticle();
    
    void updateFluidVelocity();


    void testNeighbourSearch();
    
    void searchNeighbourBruteForce(int, int* neighbourList, int* neighbourListSize); 
    
    void searchNeighbourBruteForce(double x, double y, double z, int* neighbourList, size_t& numNeiFound); 
    
    void generateGhostParticleByBruteForceNeighbourSearch();
    
    void searchNeighbourForAllParticleByBruteForceNeighbourSearch();
    
    bool checkUpwindNeighbourList();
    
    std::string rightFlush(size_t writeStep, size_t numDigits); 
    
    int writeResult(double time, size_t writeStep, size_t startIndex, size_t numParticle);
};





class HyperbolicLPSolver1D : public LPSolver {

public:
    HyperbolicLPSolver1D(const Initializer& init, ParticleData* pData);
    
    virtual int solve(double dt);   

private:

    //-----------------------------------------Data----------------------------------------

    //--------------------------Info got from input argument list---------------------------

    ParticleData* m_pParticleData; 
    
    //--------------------------------------------------------------------------------------    

    //--------------------------Info get from Initializer class------------------------------
    
    EOS* m_pEOS; 
    int m_iDimension; 
    int m_iLPFOrder; 
    std::size_t m_iNumRow2ndOrder; 
    std::size_t m_iNumRow1stOrder; 
    std::size_t m_iNumCol2ndOrder; 
    std::size_t m_iNumCol1stOrder; 
    double m_fInitParticleSpacing; 
    double m_fInvalidPressure; 
    double m_fInvalidVolume; 
    
    std::string m_sBoundaryType; 
    bool m_iUseLimiter;
    double m_fThresholdP;
    //--------------------------------------------------------------------------------------
    
    //---------------------------------Other parameters-------------------------------------
    
    double m_fPeriodicLeftBoundary; 
    double m_fPeriodicRightBoundary; 
    double m_fDisBetweenPeriodicBoundaries; 
    
    double m_fDt; 
    //bool m_iIfDebug;///< if true then print debug info
    //std::ofstream debug;///< output information for debugging     
    //-------------------------------------------------------------------------------------



    //-------------------------------------Methods-----------------------------------------
    
    void computeSetupsForNextIteration(); 


    void updateFreeBoundaryLocation();


    void updateFreeBoundaryPressureAndVelocity();


    void updateSolidBoundaryPressureAndVelocity();
    
    
    void computeMinParticleSpacing();


    void computeMaxSoundSpeed();


    void computeMaxFluidVelocity();

    
    void updateLimiter();

    void computeLPFOrder(std::size_t index, int& left_order, int& right_order); 

    
    void computeSpatialDer(int order, int direction, int i, 
                           const double *local_u, const double *local_x,
                           double& dudx, double& dudxx); //output 
    
    void solveByQR(int order, int num_nei, const double *local_u, const double *local_x, 
                  double &dudx, double &dudxx); //output

    
    void solveByCramer(int order, int num_nei, const double *local_u, const double *local_x, 
                       double &dudx, double &dudxx);    //output
    

    /*
     * \brief returns a constant weight independent of the input
     *
     * \note A helper function for solveByCramer(), which always returns 1 
     */
    double constantWeight(double h);

    void timeIntegration(int i, const double* V_old, const double* up_old,
                         const double* p_old, const double* cs_old,
                         double ux_left,  double uxx_left,
                         double px_left,  double pxx_left,
                         double ux_right, double uxx_right,
                         double px_right, double pxx_right,
                         double* V, double* up, double* p); // output
    
    
    void printInvalidState(int i, 
                           double ux_left, double uxx_left, double px_left, double pxx_left, 
                           double ux_right, double uxx_right, double px_right, double pxx_right);
    
            
    void updateFluidState(); 
    
    
    void moveFluidParticle(); 
    
    
    //void updateFluidVelocity(double*& up, double*& up_old);
    

//  bool lowerLPFOrder(int index, const std::vector<int*>& LPFOrderOther, // input
//      int* LPFOrder0, int* LPFOrder1); // output
    
//  std::string rightFlush(size_t writeStep, size_t numDigits); 
    
//  int writeResult(double time, size_t writeStep, size_t startIndex, size_t numParticle);


};





#endif // __LP_SOLVER_H__