Pressure possion equation problem of 2D flow around a cylinder. There exists pressure oscillation near the cylindrical wall surface. #288
Description
Hi @ddemidov ,
I am using AMGCL to solve a pressure possion equation problem of 2D flow around a cylinder. The processing of cylindrical boundaries adopts the virtual point embedded boundary method.
But when I use the AMGCL+BICGSTAB to solve the problem, the result is different from the exact solution obtained by the direct method.
When I use nullspace, it looks as if it doesn't work.
The data file contains coordinates, matrices, and right-side items. The first picture is the result of the exact solution. The second image is the result using amgcl but not nullspace, and the third image is the result using nullspace.
I would appreciate if you can check my code or try to solve it.
Thanks in advantage. The following is the amgcl part of my program.
#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/ilu0.hpp>
#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/solver/cg.hpp>
#include <vector>
#include <omp.h>
#include "bubble.h"
#include <amgcl/io/mm.hpp>
#include <amgcl/profiler.hpp>
#include <amgcl/adapter/eigen.hpp>
#include <amgcl/coarsening/rigid_body_modes.hpp>
using namespace Eigen;
Eigen::MatrixXd Possion(
Eigen::MatrixXd& ustar, Eigen::MatrixXd& vstar,
Eigen::SparseMatrix<double>& L,
double dt,
Eigen::MatrixXd& p,
double dx, double dy,
std::vector<double>& xm, std::vector<double>& ym,
double xc, double yc, double r,
int nx, int ny, Eigen::MatrixXd& num, int numb, std::vector<double>& b)
{
amgcl::profiler<> prof("Serena");
typedef amgcl::backend::builtin<double> Backend;
typedef amgcl::make_solver<
amgcl::amg<
Backend,
amgcl::coarsening::smoothed_aggregation,
amgcl::relaxation::ilu0
>,
amgcl::solver::bicgstab<Backend>
> Solver;
// OMP
omp_set_num_threads(omp_get_max_threads());
std::cout << "Using " << omp_get_max_threads() << " OpenMP threads.\n";
const size_t n = L.rows();
//std::vector<double> coo;
//for (int j = 1; j <= ny; ++j) {
// for (int i = 1; i <= nx; ++i) {
// double dist = std::hypot(xm[i] - xc, ym[j] - yc);
// if (dist >= r ) {
// coo.push_back(xm[i]); // x
// coo.push_back(ym[j]); // y
// }
// }
//}
//const size_t points = coo.size() / 2; //
// params
Solver::params prm;
prm.solver.tol = 1e-10;
prm.solver.maxiter = 2000;
prm.precond.coarsening.aggr.eps_strong = 0;
//prm.precond.coarsening.nullspace.cols = amgcl::coarsening::rigid_body_modes(
// 2, coo, prm.precond.coarsening.nullspace.B);
// 5. transfer to AMGCL format
std::vector<ptrdiff_t> ptr(L.outerIndexPtr(), L.outerIndexPtr() + L.outerSize() + 1);
std::vector<ptrdiff_t> col(L.innerIndexPtr(), L.innerIndexPtr() + L.nonZeros());
std::vector<double> val(L.valuePtr(), L.valuePtr() + L.nonZeros());
auto A = std::tie(n, ptr, col, val);
Solver solve(A, prm);
// x_initial and RHSb
std::vector<double> x_initial = flattenTensor(p, nx, ny, num, numb, xm, ym, xc, yc, r);
computeRHS(ustar, vstar, xm, ym, dt, dx, dy, xc, yc, r, nx, ny, num, numb, b);
int iterations;
double residual;
prof.tic("solve");
std::tie(iterations, residual) = solve(A,b, x_initial);
prof.toc("solve");
std::cout << "Iterations: " << iterations << ", Residual: " << residual << "\n";
return reshapeToTensor(x_initial, nx, ny, num, xm, ym, xc, yc, r);
}