main_petsc.cpp

 
 
#include "config.h"
 
#ifdef HAVE_PETSC
 
#include "Grid/grid_dist_id.hpp"
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include "FiniteDifference/FDScheme.hpp"
#include "FiniteDifference/util/common.hpp"
#include "FiniteDifference/eq.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "Solvers/petsc_solver.hpp"
 
struct lid_nn
{
    // dimensionaly of the equation ( 3D problem ...)
    static const unsigned int dims = 3;
 
    // number of fields in the system v_x, v_y, v_z, P so a total of 4
    static const unsigned int nvar = 4;
 
    // boundary conditions PERIODIC OR NON_PERIODIC
    static const bool boundary[];
 
    // type of space float, double, ...
    typedef float stype;
 
    // type of base grid, it is the distributed grid that will store the result
    // Note the first property is a 3D vector (velocity), the second is a scalar (Pressure)
    typedef grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> b_grid;
 
    // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
    // that you are using, in case of umfpack using <double,int> it is the only possible choice
    // We select PETSC SparseMatrix implementation
    typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
 
    // type of Vector for the linear system, this parameter is bounded by the solver
    // that you are using, in case of umfpack using <double> it is the only possible choice
    // We select PETSC implementation
    typedef Vector<double,PETSC_BASE> Vector_type;
 
    // Define that the underline grid where we discretize the system of equation is staggered
    static const int grid_type = STAGGERED_GRID;
};
 
const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
 
 
 
 
// Constant Field
struct eta
{
    typedef void const_field;
 
    static float val()  {return 1.0;}
};
 
// Model the equations
 
constexpr unsigned int v[] = {0,1,2};
constexpr unsigned int P = 3;
constexpr unsigned int ic = 3;
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
 
typedef Field<v[x],lid_nn> v_x;
typedef Field<v[y],lid_nn> v_y;
typedef Field<v[z],lid_nn> v_z;
typedef Field<P,lid_nn> Prs;
 
// Eq1 V_x
 
typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx; // Step1
typedef D<x,Prs,lid_nn> p_x; // Step 2
typedef minus<p_x,lid_nn> m_p_x; // Step3
typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq; // Step4
 
// Eq2 V_y
 
typedef mul<eta,Lap<v_y,lid_nn>,lid_nn> eta_lap_vy;
typedef D<y,Prs,lid_nn> p_y;
typedef minus<p_y,lid_nn> m_p_y;
typedef sum<eta_lap_vy,m_p_y,lid_nn> vy_eq;
 
// Eq3 V_z
 
typedef mul<eta,Lap<v_z,lid_nn>,lid_nn> eta_lap_vz;
typedef D<z,Prs,lid_nn> p_z;
typedef minus<p_z,lid_nn> m_p_z;
typedef sum<eta_lap_vz,m_p_z,lid_nn> vz_eq;
 
// Eq4 Incompressibility
 
typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step5
typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step6
typedef D<z,v_z,lid_nn,FORWARD> dz_vz; // Step 7
typedef sum<dx_vx,dy_vy,dz_vz,lid_nn> ic_eq; // Step8
 
 
 
// Directional Avg
typedef Avg<x,v_y,lid_nn> avg_x_vy;
typedef Avg<z,v_y,lid_nn> avg_z_vy;
 
typedef Avg<y,v_x,lid_nn> avg_y_vx;
typedef Avg<z,v_x,lid_nn> avg_z_vx;
 
typedef Avg<y,v_z,lid_nn> avg_y_vz;
typedef Avg<x,v_z,lid_nn> avg_x_vz;
 
typedef Avg<x,v_y,lid_nn,FORWARD> avg_x_vy_f;
typedef Avg<z,v_y,lid_nn,FORWARD> avg_z_vy_f;
 
typedef Avg<y,v_x,lid_nn,FORWARD> avg_y_vx_f;
typedef Avg<z,v_x,lid_nn,FORWARD> avg_z_vx_f;
 
typedef Avg<y,v_z,lid_nn,FORWARD> avg_y_vz_f;
typedef Avg<x,v_z,lid_nn,FORWARD> avg_x_vz_f;
 
 
#define EQ_1 0
#define EQ_2 1
#define EQ_3 2
#define EQ_4 3
 
 
#include "Vector/vector_dist.hpp"
#include "data_type/aggregate.hpp"
 
int main(int argc, char* argv[])
{
 
    // Initialize
    openfpm_init(&argc,&argv);
    {
    // velocity in the grid is the property 0, pressure is the property 1
    constexpr int velocity = 0;
    constexpr int pressure = 1;
 
    // Domain
    Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
 
    // Ghost (Not important in this case but required)
    Ghost<3,float> g(0.01);
 
    // Grid points on x=36 and y=12 z=12
    long int sz[] = {36,12,12};
    size_t szu[3];
    szu[0] = (size_t)sz[0];
    szu[1] = (size_t)sz[1];
    szu[2] = (size_t)sz[2];
 
 
    // We need one more point on the left and down part of the domain
    // This is given by the boundary conditions that we impose.
    //
    Padding<3> pd({1,1,1},{0,0,0});
 
 
 
    grid_dist_id<3,float,aggregate<float[3],float>> g_dist(szu,domain,g);
 
 
 
    // It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
    Ghost<3,long int> stencil_max(1);
 
    // Finite difference scheme
    FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
 
 
 
    // start and end of the bulk
 
    fd.impose(ic_eq(),0.0, EQ_4, {0,0,0},{sz[0]-2,sz[1]-2,sz[2]-2},true);
    fd.impose(Prs(),  0.0, EQ_4, {0,0,0},{0,0,0});
    fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2,sz[2]-2});
    fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
    fd.impose(vz_eq(),0.0, EQ_3, {0,0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
 
    // v_x
    // R L (Right,Left)
    fd.impose(v_x(),0.0, EQ_1, {0,0,0},      {0,sz[1]-2,sz[2]-2});
    fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-2});
 
    // T B (Top,Bottom)
    fd.impose(avg_y_vx_f(),0.0, EQ_1, {0,-1,0},     {sz[0]-1,-1,sz[2]-2});
    fd.impose(avg_y_vx(),0.0, EQ_1,   {0,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-2});
 
    // A F (Forward,Backward)
    fd.impose(avg_z_vx_f(),0.0, EQ_1, {0,-1,-1},     {sz[0]-1,sz[1]-1,-1});
    fd.impose(avg_z_vx(),0.0, EQ_1, {0,-1,sz[2]-1},{sz[0]-1,sz[1]-1,sz[2]-1});
 
    // v_y
    // R L
    fd.impose(avg_x_vy_f(),0.0, EQ_2,  {-1,0,0},     {-1,sz[1]-1,sz[2]-2});
    fd.impose(avg_x_vy(),1.0, EQ_2,    {sz[0]-1,0,0},{sz[0]-1,sz[1]-1,sz[2]-2});
 
    // T B
    fd.impose(v_y(), 0.0, EQ_2, {0,0,0},      {sz[0]-2,0,sz[2]-2});
    fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1,0},{sz[0]-2,sz[1]-1,sz[2]-2});
 
    // F A
    fd.impose(avg_z_vy(),0.0, EQ_2,   {-1,0,sz[2]-1}, {sz[0]-1,sz[1]-1,sz[2]-1});
    fd.impose(avg_z_vy_f(),0.0, EQ_2, {-1,0,-1},      {sz[0]-1,sz[1]-1,-1});
 
    // v_z
    // R L
    fd.impose(avg_x_vz_f(),0.0, EQ_3, {-1,0,0},     {-1,sz[1]-2,sz[2]-1});
    fd.impose(avg_x_vz(),1.0, EQ_3,   {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-1});
 
    // T B
    fd.impose(avg_y_vz(),0.0, EQ_3, {-1,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-1});
    fd.impose(avg_y_vz_f(),0.0, EQ_3, {-1,-1,0},   {sz[0]-1,-1,sz[2]-1});
 
    // F A
    fd.impose(v_z(),0.0, EQ_3, {0,0,0},      {sz[0]-2,sz[1]-2,0});
    fd.impose(v_z(),0.0, EQ_3, {0,0,sz[2]-1},{sz[0]-2,sz[1]-2,sz[2]-1});
 
    // When we pad the grid, there are points of the grid that are not
    // touched by the previous condition. Mathematically this lead
    // to have too many variables for the conditions that we are imposing.
    // Here we are imposing variables that we do not touch to zero
    //
 
    // L R
    fd.impose(Prs(), 0.0, EQ_4, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
    fd.impose(Prs(), 0.0, EQ_4, {sz[0]-1,-1,-1},{sz[0]-1,sz[1]-1,sz[2]-1});
 
    // T B
    fd.impose(Prs(), 0.0, EQ_4, {0,sz[1]-1,-1}, {sz[0]-2,sz[1]-1,sz[2]-1});
    fd.impose(Prs(), 0.0, EQ_4, {0,-1     ,-1}, {sz[0]-2,-1,     sz[2]-1});
 
    // F A
    fd.impose(Prs(), 0.0, EQ_4, {0,0,sz[2]-1}, {sz[0]-2,sz[1]-2,sz[2]-1});
    fd.impose(Prs(), 0.0, EQ_4, {0,0,-1},      {sz[0]-2,sz[1]-2,-1});
 
    // Impose v_x  v_y v_z padding
    fd.impose(v_x(), 0.0, EQ_1, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
    fd.impose(v_y(), 0.0, EQ_2, {-1,-1,-1},{sz[0]-1,-1,sz[2]-1});
    fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});
 
 
 
    // Create an PETSC solver
    petsc_solver<double> solver;
 
    // Warning try many solver and collect statistics require a lot of time
    // To just solve you can comment this line
//  solver.best_solve();
 
    // Give to the solver A and b, return x, the solution
    auto x = solver.solve(fd.getA(),fd.getB());
 
 
 
    // copy the solution to grid
    fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
 
    g_dist.write("lid_driven_cavity_p_petsc");
 
 
 
    }
    openfpm_finalize();
 
 
 
}
 
#else
 
int main(int argc, char* argv[])
{
        return 0;
}
 
#endif