From 6f2388aab59ec96b949b5418ec5b2c2f05e44092 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Mon, 20 Apr 2020 17:32:33 +0200
Subject: [PATCH] Writer for SparseMatrix + Fixing DCPSE.hpp for parallel
---
mpi_include | 2 +-
mpi_libs | 2 +-
src/CMakeLists.txt | 2 +-
src/DCPSE_op/DCPSE_Solver.hpp | 60 +++++++++-
src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp | 5 +-
src/DCPSE_op/DCPSE_op_test.cpp | 124 ++++++++++++++++++++-
src/DCPSE_op/EqnsStructEigen.hpp | 12 ++
src/DCPSE_op/EqnsStructPetsc.hpp | 24 +++-
src/Matrix/SparseMatrix_petsc.hpp | 41 +++++++
9 files changed, 253 insertions(+), 19 deletions(-)
diff --git a/mpi_include b/mpi_include
index 16981421..4439b1ca 100644
--- a/mpi_include
+++ b/mpi_include
@@ -1 +1 @@
--I
\ No newline at end of file
+-I/home/i-bird/MPI/include
\ No newline at end of file
diff --git a/mpi_libs b/mpi_libs
index 0519ecba..3a37a5cd 100644
--- a/mpi_libs
+++ b/mpi_libs
@@ -1 +1 @@
-
\ No newline at end of file
+-Wl,-rpath -Wl,/home/i-bird/MPI/lib -Wl,--enable-new-dtags -pthread /home/i-bird/MPI/lib/libmpi.so
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 124f6db8..44efcfaa 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
########################### Executables
-add_executable(numerics DCPSE_op/DCPSE_op_test3d.cpp DCPSE_op/DCPSE_op_test2.cpp DCPSE_op/DCPSE_op_test.cpp main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp Vector/Vector_unit_tests.cpp Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp DCPSE_op/DCPSE_Solver.hpp MatrixAssembler/MatrixAssembler.hpp DCPSE_op/EqnsStruct.hpp)
+add_executable(numerics DCPSE_op/DCPSE_op_test3d.cpp DCPSE_op/DCPSE_op_test2.cpp DCPSE_op/DCPSE_op_test.cpp main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp Vector/Vector_unit_tests.cpp Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp )
###########################
diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index ce939d08..1d93b4c4 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -295,7 +295,7 @@ class DCPSE_scheme: public MatrixAssembler
while (it.isNext())
{
auto p = it.get();
- exp.value(p) = x(p.getKey()*Sys_eqs::nvar + comp);
+ exp.value(p) = x(p.getKey()*Sys_eqs::nvar + comp + s_pnt*Sys_eqs::nvar);
++it;
}
}
@@ -353,19 +353,61 @@ public:
*/
template<typename ... expr_type>
void solve(expr_type ... exps)
+ {
+ if (sizeof...(exps) != Sys_eqs::nvar)
+ {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
+ dimensionality, I am expecting " << Sys_eqs::nvar <<
+ " properties " << std::endl;};
+ typename Sys_eqs::solver_type solver;
+// umfpack_solver<double> solver;
+ auto x = solver.solve(getA(opt),getB(opt));
+
+ unsigned int comp = 0;
+ copy_nested(x,comp,exps ...);
+ }
+
+ /*! \brief Solve an equation
+ *
+ * \warning exp must be a scalar type
+ *
+ * \param exp where to store the result
+ *
+ */
+ template<typename SolverType, typename ... expr_type>
+ void solve_with_solver(SolverType & solver, expr_type ... exps)
{
if (sizeof...(exps) != Sys_eqs::nvar)
{std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
dimensionality, I am expecting " << Sys_eqs::nvar <<
" properties " << std::endl;};
- umfpack_solver<double> solver;
auto x = solver.solve(getA(opt),getB(opt));
unsigned int comp = 0;
copy_nested(x,comp,exps ...);
}
+ /*! \brief Solve an equation
+ *
+ * \warning exp must be a scalar type
+ *
+ * \param exp where to store the result
+ *
+ */
+ template<typename SolverType, typename ... expr_type>
+ void try_solve_with_solver(SolverType & solver, expr_type ... exps)
+ {
+ if (sizeof...(exps) != Sys_eqs::nvar)
+ {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
+ dimensionality, I am expecting " << Sys_eqs::nvar <<
+ " properties " << std::endl;};
+
+ auto x = solver.try_solve(getA(opt),getB(opt));
+
+ unsigned int comp = 0;
+ copy_nested(x,comp,exps ...);
+ }
+
/*! \brief Constructor
*
* The periodicity is given by the grid b_g
@@ -506,6 +548,14 @@ public:
trpl.add(t2);
}
+ triplet t3;
+
+ t3.col() = tot * Sys_eqs::nvar;
+ t3.row() = tot * Sys_eqs::nvar;
+ t3.value() = 0;
+
+ trpl.add(t3);
+
row_b++;
row++;
}
@@ -577,6 +627,12 @@ public:
// get the particle
auto key = it.get();
+ if (key == 298 && create_vcluster().rank() == 1)
+ {
+ int debug = 0;
+ debug++;
+ }
+
// Calculate the non-zero colums
typename Sys_eqs::stype coeff = 1.0;
op.template value_nz<Sys_eqs>(p_map,key,cols,coeff,0);
diff --git a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
index 8dfd39a4..ad34d8e0 100644
--- a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
+++ b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
@@ -11,6 +11,7 @@
#include <Space/Shape/Point.hpp>
#include <Vector/vector_dist.hpp>
#include "Support.hpp"
+#include <utility>
enum support_options
{
@@ -23,7 +24,7 @@ class SupportBuilder
{
private:
vector_type &domain;
- CellList<vector_type::dims, typename vector_type::stype, Mem_fast<HeapMemory, local_index_>> cellList;
+ decltype(std::declval<vector_type>().getCellList(0.0)) cellList;
const Point<vector_type::dims, unsigned int> differentialSignature;
typename vector_type::stype rCut;
@@ -57,7 +58,7 @@ SupportBuilder<vector_type>::SupportBuilder(vector_type &domain, const Point<vec
differentialSignature(differentialSignature),
rCut(rCut)
{
- cellList = domain.template getCellList<CellList<vector_type::dims, typename vector_type::stype, Mem_fast<HeapMemory, local_index_>>>(rCut);
+ cellList = domain.getCellList(rCut);
}
template<typename vector_type>
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 59e23e1d..0c4c19ac 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -974,7 +974,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin_anal) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- const size_t sz[2] = {161,161};
+ const size_t sz[2] = {50,50};
Box<2, double> box({0, 0}, {0.5, 0.5});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -1000,15 +1000,96 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
++it;
}
+
+ // Add multi res patch 1
+
+ {
+ const size_t sz2[2] = {40,40};
+ Box<2,double> bx({0.25 + it.getSpacing(0)/4.0,0.25 + it.getSpacing(0)/4.0},{sz2[0]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0, sz2[1]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0});
+ openfpm::vector<size_t> rem;
+
+ auto it = domain.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto k = it.get();
+
+ Point<2,double> xp = domain.getPos(k);
+
+ if (bx.isInside(xp) == true)
+ {
+ rem.add(k.getKey());
+ }
+
+ ++it;
+ }
+
+ domain.remove(rem);
+
+ auto it2 = domain.getGridIterator(sz2);
+ while (it2.isNext()) {
+ domain.add();
+
+ auto key = it2.get();
+ double x = key.get(0) * spacing/2.0 + 0.25 + spacing/4.0;
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * spacing/2.0 + 0.25 + spacing/4.0;
+ domain.getLastPos()[1] = y;
+
+ ++it2;
+ }
+ }
+
+ // Add multi res patch 2
+
+ {
+ const size_t sz2[2] = {40,40};
+ Box<2,double> bx({0.25 + 21.0*spacing/8.0,0.25 + 21.0*spacing/8.0},{sz2[0]*spacing/4.0 + 0.25 + 21.0*spacing/8.0, sz2[1]*spacing/4.0 + 0.25 + 21*spacing/8.0});
+ openfpm::vector<size_t> rem;
+
+ auto it = domain.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto k = it.get();
+
+ Point<2,double> xp = domain.getPos(k);
+
+ if (bx.isInside(xp) == true)
+ {
+ rem.add(k.getKey());
+ }
+
+ ++it;
+ }
+
+ domain.remove(rem);
+
+ auto it2 = domain.getGridIterator(sz2);
+ while (it2.isNext()) {
+ domain.add();
+
+ auto key = it2.get();
+ double x = key.get(0) * spacing/4.0 + 0.25 + 21*spacing/8.0;
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * spacing/4.0 + 0.25 + 21*spacing/8.0;
+ domain.getLastPos()[1] = y;
+
+ ++it2;
+ }
+ }
+
+ ///////////////////////
+
BOOST_TEST_MESSAGE("Sync domain across processors...");
domain.map();
domain.ghost_get<0>();
- Derivative_x Dx(domain, 2, rCut,1.9,support_options::RADIUS);
- Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
+ Derivative_x Dx(domain, 2, rCut / 3.0 ,1.9/*,support_options::RADIUS*/);
+ Derivative_y Dy(domain, 2, rCut / 3.0 ,1.9/*,support_options::RADIUS*/);
//Gradient Grad(domain, 2, rCut);
- Laplacian Lap(domain, 2, rCut,1.9,support_options::RADIUS);
+ Laplacian Lap(domain, 2, rCut / 3.0 ,1.9/*,support_options::RADIUS*/);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
@@ -1090,6 +1171,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
++it2;
}
+
+ domain.ghost_get<1,3>();
+
DCPSE_scheme<equations2d1,decltype(domain)> Solver( domain);
auto Poisson = Lap(v);
auto D_x = Dx(v);
@@ -1099,8 +1183,35 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.impose(D_x, r_p, 0);
Solver.impose(v, dw_p, 0);
Solver.impose(v, l_p, 0);
- Solver.solve(sol);
+
+ petsc_solver<double> solver;
+
+ solver.setRestart(500);
+
+// auto A = Solver.getA();
+// A.write("Matrix_anal_sol");
+
+// solver.setSolver(KSPGMRES);
+// solver.setPreconditioner(PCKSP);
+// solver.setRestart(300);
+
+/* solver.setPreconditioner(PCHYPRE_BOOMERAMG);
+ solver.setPreconditionerAMG_nl(6);
+ solver.setPreconditionerAMG_maxit(10);
+ solver.setPreconditionerAMG_relax("SOR/Jacobi");
+ solver.setPreconditionerAMG_cycleType("V",0,4);
+ solver.setPreconditionerAMG_coarsenNodalType(0);
+ solver.setPreconditionerAMG_coarsen("HMIS");*/
+
+// solver.log_monitor();
+
+ Solver.solve_with_solver(solver,sol);
+
+ domain.ghost_get<2>();
+
DCPSE_sol=Lap(sol);
+ domain.ghost_get<5>();
+
double worst1 = 0.0;
v=abs(DCPSE_sol-RHS);
@@ -1115,6 +1226,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
+ domain.ghost_get<4>();
domain.write("Robin_anasol");
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1123,7 +1235,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(dcpse_poisson_Dirichlet_anal) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- const size_t sz[2] = {81,81};
+ const size_t sz[2] = {512,512};
Box<2, double> box({0, 0}, {1, 1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
diff --git a/src/DCPSE_op/EqnsStructEigen.hpp b/src/DCPSE_op/EqnsStructEigen.hpp
index c5c954e4..36d2e3bc 100644
--- a/src/DCPSE_op/EqnsStructEigen.hpp
+++ b/src/DCPSE_op/EqnsStructEigen.hpp
@@ -27,6 +27,8 @@ struct equations2d1 {
//! type of Vector for the linear solver
typedef Vector<double> Vector_type;
+
+ typedef umfpack_solver<double> solver_type;
};
struct equations2d2 {
@@ -49,6 +51,8 @@ struct equations2d2 {
//! type of Vector for the linear solver
typedef Vector<double> Vector_type;
+
+ typedef umfpack_solver<double> solver_type;
};
@@ -72,6 +76,8 @@ struct equations2d1p {
//! type of Vector for the linear solver
typedef Vector<double> Vector_type;
+
+ typedef umfpack_solver<double> solver_type;
};
struct equations2d2p {
@@ -94,6 +100,8 @@ struct equations2d2p {
//! type of Vector for the linear solver
typedef Vector<double> Vector_type;
+
+ typedef umfpack_solver<double> solver_type;
};
struct equations3d3 {
@@ -116,6 +124,8 @@ struct equations3d3 {
//! type of Vector for the linear solver
typedef Vector<double> Vector_type;
+
+ typedef umfpack_solver<double> solver_type;
};
struct equations3d1 {
@@ -138,6 +148,8 @@ struct equations3d1 {
//! type of Vector for the linear solver
typedef Vector<double> Vector_type;
+
+ typedef umfpack_solver<double> solver_type;
};
diff --git a/src/DCPSE_op/EqnsStructPetsc.hpp b/src/DCPSE_op/EqnsStructPetsc.hpp
index 51edf411..5ea1e106 100644
--- a/src/DCPSE_op/EqnsStructPetsc.hpp
+++ b/src/DCPSE_op/EqnsStructPetsc.hpp
@@ -26,7 +26,9 @@ struct equations2d1 {
typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
//! type of Vector for the linear solver
- typedef Vector<double> Vector_type;
+ typedef Vector<double, PETSC_BASE> Vector_type;
+
+ typedef petsc_solver<double> solver_type;
};
struct equations2d2 {
@@ -48,7 +50,9 @@ struct equations2d2 {
typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
//! type of Vector for the linear solver
- typedef Vector<double> Vector_type;
+ typedef Vector<double, PETSC_BASE> Vector_type;
+
+ typedef petsc_solver<double> solver_type;
};
@@ -71,7 +75,9 @@ struct equations2d1p {
typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
//! type of Vector for the linear solver
- typedef Vector<double> Vector_type;
+ typedef Vector<double, PETSC_BASE> Vector_type;
+
+ typedef petsc_solver<double> solver_type;
};
struct equations2d2p {
@@ -93,7 +99,9 @@ struct equations2d2p {
typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
//! type of Vector for the linear solver
- typedef Vector<double> Vector_type;
+ typedef Vector<double, PETSC_BASE> Vector_type;
+
+ typedef petsc_solver<double> solver_type;
};
struct equations3d3 {
@@ -115,7 +123,9 @@ struct equations3d3 {
typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
//! type of Vector for the linear solver
- typedef Vector<double> Vector_type;
+ typedef Vector<double, PETSC_BASE> Vector_type;
+
+ typedef petsc_solver<double> solver_type;
};
struct equations3d1 {
@@ -137,7 +147,9 @@ struct equations3d1 {
typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
//! type of Vector for the linear solver
- typedef Vector<double> Vector_type;
+ typedef Vector<double, PETSC_BASE> Vector_type;
+
+ typedef petsc_solver<double> solver_type;
};
diff --git a/src/Matrix/SparseMatrix_petsc.hpp b/src/Matrix/SparseMatrix_petsc.hpp
index 3534924b..a03b0671 100644
--- a/src/Matrix/SparseMatrix_petsc.hpp
+++ b/src/Matrix/SparseMatrix_petsc.hpp
@@ -12,6 +12,7 @@
#include "Vector/map_vector.hpp"
#include <boost/mpl/int.hpp>
#include <petscmat.h>
+#include "VTKWriter/VTKWriter.hpp"
#define PETSC_BASE 2
@@ -352,6 +353,46 @@ public:
return 0;
}
+
+ /* Write matrix on vtk
+ *
+ *
+ */
+ bool write(std::string out)
+ {
+ Vcluster<> & v_cl = create_vcluster();
+
+ openfpm::vector<Point<2, double>> row_col;
+ openfpm::vector<aggregate<double>> values;
+
+ row_col.resize(trpl.size());
+ values.resize(trpl.size());
+
+ for (int i = 0 ; i < trpl.size() ; i++)
+ {
+ row_col.template get<0>(i)[1] = trpl.get(i).row();
+ row_col.template get<0>(i)[0] = trpl.get(i).col();
+
+ values.template get<0>(i) = trpl.get(i).value();
+ }
+
+ auto ft = file_type::BINARY;
+
+ // VTKWriter for a set of points
+ VTKWriter<boost::mpl::pair<openfpm::vector<Point<2, double>>,
+ openfpm::vector<aggregate<double>>>,
+ VECTOR_POINTS> vtk_writer;
+
+ vtk_writer.add(row_col,values,row_col.size());
+
+ std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + std::to_string(".vtk"));
+
+ openfpm::vector<std::string> prp_names;
+ prp_names.add("value");
+
+ // Write the VTK file
+ return vtk_writer.write(output,prp_names,"particles","",ft);
+ }
};
--
GitLab