OpenFPM  5.2.0
Project that contain the implementation of distributed structures
Lid driven cavity with Pressure Correction and FD

Lid Driven Cavity Problem with Pressure Correction and FD

In this example, we solve the incompressible Navier-Stokes equation in a 2D Square Box:

\[ \mathbf{v}\cdot(\nabla \mathbf{v})-\frac{1}{\text{Re}}\mathrm{\Delta} \mathbf{v}=-\nabla \Pi \label{eq:NS} \\ \nabla\cdot \mathbf{v}=0 \\ \mathbf{v}(x_b,y_b)=(0,0), \text{ except } \mathbf{v}(x_b,1)=(1,0)\, , \]

We do that by solving the implicit stokes equation and and employing an iterative pressure correction scheme:

Output: Steady State Solution to the Lid Driven Cavity Problem.

Including the headers

These are the header files that we need to include:

// Include Vector Expression,Vector Expressions for Subset,FD , and Solver header files
#include "config.h"
#include <iostream>
#include "FiniteDifference/FD_Solver.hpp"
#include "FiniteDifference/FD_op.hpp"
#include "FiniteDifference/FD_expressions.hpp"
#include "FiniteDifference/util/EqnsStructFD.hpp"

Initialization

  • Initialize the library
  • Define some useful constants
  • define Ghost size
  • Non-periodic boundary conditions
openfpm_init(&argc,&argv);
using namespace FD;
timer tt2;
tt2.start();
size_t gd_sz = 81;
constexpr int x = 0;
constexpr int y = 1;
const size_t szu[2] = {gd_sz,gd_sz};
int sz[2]={int(gd_sz),int(gd_sz)};
Box<2, double> box({0, 0}, {1,1});
periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
double spacing;
spacing = 1.0 / (sz[0] - 1);
Ghost<2,long int> ghost(1);
auto &v_cl = create_vcluster();
typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double> LidCavity;
grid_dist_id<2, double, LidCavity> domain(szu, box,ghost,bc);
double x0, y0, x1, y1;
x0 = box.getLow(0);
y0 = box.getLow(1);
x1 = box.getHigh(0);
y1 = box.getHigh(1);
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
Ghost
Definition: Ghost.hpp:40
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:28
grid_dist_id
This is a distributed grid.
Definition: grid_dist_id.hpp:279
timer
Class for cpu time benchmarking.
Definition: timer.hpp:28
timer::start
void start()
Start the timer.
Definition: timer.hpp:90
cub::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition: dispatch_radix_sort.cuh:336
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:221
periodicity
Boundary conditions.
Definition: common.hpp:22

Creating particles in the 2D domain

We set the appropriate subset number 0 for bulk and other for boundary. Note that for different walls we need different subsets as for the pressure, we need normal derivative zero condition.

auto it = domain.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
auto gkey = it.getGKey(key);
double x = gkey.get(0) * domain.spacing(0);
double y = gkey.get(1) * domain.spacing(1);
if (y==1)
{domain.get<1>(key) = 1.0;}
else
{domain.get<1>(key) = 0.0;}
++it;
}
domain.ghost_get<0>();

Creating Subsets and Vector Expressions for Fields and Differential Operators for easier encoding.

Derivative_x Dx;
Derivative_y Dy;
Derivative_xx Dxx;
Derivative_yy Dyy;
auto P = getV<0>(domain);
auto V = getV<1>(domain);
auto RHS = getV<2>(domain);
auto dV = getV<3>(domain);
auto div = getV<4>(domain);
auto V_star = getV<5>(domain);
auto H = getV<6>(domain);
auto dP = getV<7>(domain);
Point_test
Test structure used for several test.
Definition: Point_test.hpp:106

Creating a 3D implicit solver for the given set of particles and iteratively solving wit pressure correction.

eq_id x_comp,y_comp;
x_comp.setId(0);
y_comp.setId(1);
double sum, sum1, sum_k,V_err_old;
auto StokesX=(V[x]*Dx(V_star[x])+V[y]*Dy(V_star[x]))-(1.0/100.0)*(Dxx(V_star[x])+Dyy(V_star[x]));
auto StokesY=(V[x]*Dx(V_star[y])+V[y]*Dy(V_star[y]))-(1.0/100.0)*(Dxx(V_star[y])+Dyy(V_star[y]));
RHS=0;
dV=0;
double V_err=1;
int n;
while (V_err >= 0.05 && n <= 50) {
if (n%5==0){
domain.ghost_get<0,1>(SKIP_LABELLING);
domain.write_frame("LID",n,BINARY);
domain.ghost_get<0>();
}
domain.ghost_get<0>(SKIP_LABELLING);
RHS[x] = -Dx(P);
RHS[y] = -Dy(P);
FD_scheme<equations2d2,decltype(domain)> Solver(ghost,domain);
Solver.impose(StokesX, {1,1},{sz[0]-2,sz[1]-2}, RHS[x], x_comp);
Solver.impose(StokesY, {1,1},{sz[0]-2,sz[1]-2}, RHS[y], y_comp);
Solver.impose(V_star[x], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 1.0, x_comp);
Solver.impose(V_star[y], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 0.0, y_comp);
Solver.impose(V_star[x], {0,1},{0,sz[1]-2}, 0.0, x_comp);
Solver.impose(V_star[y], {0,1},{0,sz[1]-2}, 0.0, y_comp);
Solver.impose(V_star[x], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, x_comp);
Solver.impose(V_star[y], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, y_comp);
Solver.impose(V_star[x], {0,0},{sz[0]-1,0}, 0.0, x_comp);
Solver.impose(V_star[y], {0,0},{sz[0]-1,0}, 0.0, y_comp);
Solver.solve(V_star[x], V_star[y]);
domain.ghost_get<5>(SKIP_LABELLING);
div = (Dx(V_star[x]) + Dy(V_star[y]));
FD_scheme<equations2d1E,decltype(domain)> SolverH(ghost,domain);
auto Helmholtz = Dxx(H)+Dyy(H);
SolverH.impose(Helmholtz,{1,1},{sz[0]-2,sz[1]-2},div);
SolverH.impose(Dy(H), {0,sz[0]-1},{sz[0]-1,sz[1]-1},0);
SolverH.impose(Dx(H), {sz[0]-1,1},{sz[0]-1,sz[1]-2},0);
SolverH.impose(Dx(H), {0,0},{sz[0]-1,0},0,x_comp,true);
SolverH.impose(H, {0,0},{0,0},0);
SolverH.impose(Dy(H), {0,1},{0,sz[1]-2},0);
SolverH.solve(H);
domain.ghost_get<1,4,6>(SKIP_LABELLING);
P = P - 0.01*(div-0.5*(V[x]*Dx(H)+V[y]*Dy(H)));
V_star[0] = V_star[0] - Dx(H);
V_star[1] = V_star[1] - Dy(H);
sum = 0;
sum1 = 0;
auto it2 = domain.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
sum += (domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) *
(domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) +
(domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]) *
(domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]);
sum1 += domain.getProp<5>(p)[0] * domain.getProp<5>(p)[0] +
domain.getProp<5>(p)[1] * domain.getProp<5>(p)[1];
++it2;
}
V[x] = V_star[x];
V[y] = V_star[y];
v_cl.sum(sum);
v_cl.sum(sum1);
v_cl.execute();
sum = sqrt(sum);
sum1 = sqrt(sum1);
V_err_old = V_err;
V_err = sum / sum1;
n++;
if (v_cl.rank() == 0) {
std::cout << "Rel l2 cgs err in V = " << V_err << "."<< std::endl;
}
}
domain.write("LID_final");
tt2.stop();
if (v_cl.rank() == 0) {
std::cout << "The entire solver took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
<< "(Wall) Seconds.";
}
}
openfpm_finalize();
FD_scheme
Finite Differences.
Definition: FD_Solver.hpp:128
eq_id
Definition: eq_solve_common.hpp:15
timer::stop
void stop()
Stop the timer.
Definition: timer.hpp:119
timer::getcputime
double getcputime()
Return the cpu time.
Definition: timer.hpp:142
timer::getwct
double getwct()
Return the elapsed real time.
Definition: timer.hpp:130
equations2d1E
Specify the general characteristic of system to solve.
Definition: EqnsStruct.hpp:648
sum
It model an expression expr1 + ... exprn.
Definition: sum.hpp:93

Full code

//
// Created by Abhinav Singh on 15.11.2021.
//
// Include Vector Expression,Vector Expressions for Subset,FD , and Solver header files
#include "config.h"
#include <iostream>
#include "FiniteDifference/FD_Solver.hpp"
#include "FiniteDifference/FD_op.hpp"
#include "FiniteDifference/FD_expressions.hpp"
#include "FiniteDifference/util/EqnsStructFD.hpp"
int main(int argc, char* argv[])
{
{
openfpm_init(&argc,&argv);
using namespace FD;
timer tt2;
tt2.start();
size_t gd_sz = 81;
constexpr int x = 0;
constexpr int y = 1;
const size_t szu[2] = {gd_sz,gd_sz};
int sz[2]={int(gd_sz),int(gd_sz)};
Box<2, double> box({0, 0}, {1,1});
periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
double spacing;
spacing = 1.0 / (sz[0] - 1);
Ghost<2,long int> ghost(1);
auto &v_cl = create_vcluster();
typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double> LidCavity;
grid_dist_id<2, double, LidCavity> domain(szu, box,ghost,bc);
double x0, y0, x1, y1;
x0 = box.getLow(0);
y0 = box.getLow(1);
x1 = box.getHigh(0);
y1 = box.getHigh(1);
auto it = domain.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
auto gkey = it.getGKey(key);
double x = gkey.get(0) * domain.spacing(0);
double y = gkey.get(1) * domain.spacing(1);
if (y==1)
{domain.get<1>(key) = 1.0;}
else
{domain.get<1>(key) = 0.0;}
++it;
}
domain.ghost_get<0>();
Derivative_x Dx;
Derivative_y Dy;
Derivative_xx Dxx;
Derivative_yy Dyy;
auto P = getV<0>(domain);
auto V = getV<1>(domain);
auto RHS = getV<2>(domain);
auto dV = getV<3>(domain);
auto div = getV<4>(domain);
auto V_star = getV<5>(domain);
auto H = getV<6>(domain);
auto dP = getV<7>(domain);
eq_id x_comp,y_comp;
x_comp.setId(0);
y_comp.setId(1);
double sum, sum1, sum_k,V_err_old;
auto StokesX=(V[x]*Dx(V_star[x])+V[y]*Dy(V_star[x]))-(1.0/100.0)*(Dxx(V_star[x])+Dyy(V_star[x]));
auto StokesY=(V[x]*Dx(V_star[y])+V[y]*Dy(V_star[y]))-(1.0/100.0)*(Dxx(V_star[y])+Dyy(V_star[y]));
RHS=0;
dV=0;
double V_err=1;
int n;
while (V_err >= 0.05 && n <= 50) {
if (n%5==0){
domain.ghost_get<0,1>(SKIP_LABELLING);
domain.write_frame("LID",n,BINARY);
domain.ghost_get<0>();
}
domain.ghost_get<0>(SKIP_LABELLING);
RHS[x] = -Dx(P);
RHS[y] = -Dy(P);
FD_scheme<equations2d2,decltype(domain)> Solver(ghost,domain);
Solver.impose(StokesX, {1,1},{sz[0]-2,sz[1]-2}, RHS[x], x_comp);
Solver.impose(StokesY, {1,1},{sz[0]-2,sz[1]-2}, RHS[y], y_comp);
Solver.impose(V_star[x], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 1.0, x_comp);
Solver.impose(V_star[y], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 0.0, y_comp);
Solver.impose(V_star[x], {0,1},{0,sz[1]-2}, 0.0, x_comp);
Solver.impose(V_star[y], {0,1},{0,sz[1]-2}, 0.0, y_comp);
Solver.impose(V_star[x], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, x_comp);
Solver.impose(V_star[y], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, y_comp);
Solver.impose(V_star[x], {0,0},{sz[0]-1,0}, 0.0, x_comp);
Solver.impose(V_star[y], {0,0},{sz[0]-1,0}, 0.0, y_comp);
Solver.solve(V_star[x], V_star[y]);
domain.ghost_get<5>(SKIP_LABELLING);
div = (Dx(V_star[x]) + Dy(V_star[y]));
FD_scheme<equations2d1E,decltype(domain)> SolverH(ghost,domain);
auto Helmholtz = Dxx(H)+Dyy(H);
SolverH.impose(Helmholtz,{1,1},{sz[0]-2,sz[1]-2},div);
SolverH.impose(Dy(H), {0,sz[0]-1},{sz[0]-1,sz[1]-1},0);
SolverH.impose(Dx(H), {sz[0]-1,1},{sz[0]-1,sz[1]-2},0);
SolverH.impose(Dx(H), {0,0},{sz[0]-1,0},0,x_comp,true);
SolverH.impose(H, {0,0},{0,0},0);
SolverH.impose(Dy(H), {0,1},{0,sz[1]-2},0);
SolverH.solve(H);
domain.ghost_get<1,4,6>(SKIP_LABELLING);
P = P - 0.01*(div-0.5*(V[x]*Dx(H)+V[y]*Dy(H)));
V_star[0] = V_star[0] - Dx(H);
V_star[1] = V_star[1] - Dy(H);
sum = 0;
sum1 = 0;
auto it2 = domain.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
sum += (domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) *
(domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) +
(domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]) *
(domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]);
sum1 += domain.getProp<5>(p)[0] * domain.getProp<5>(p)[0] +
domain.getProp<5>(p)[1] * domain.getProp<5>(p)[1];
++it2;
}
V[x] = V_star[x];
V[y] = V_star[y];
v_cl.sum(sum);
v_cl.sum(sum1);
v_cl.execute();
sum = sqrt(sum);
sum1 = sqrt(sum1);
V_err_old = V_err;
V_err = sum / sum1;
n++;
if (v_cl.rank() == 0) {
std::cout << "Rel l2 cgs err in V = " << V_err << "."<< std::endl;
}
}
domain.write("LID_final");
tt2.stop();
if (v_cl.rank() == 0) {
std::cout << "The entire solver took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
<< "(Wall) Seconds.";
}
}
openfpm_finalize();
}