From 7cb5820d48dd0ef7adbe666bbcdf8fd45cad4374 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Mon, 3 Jul 2017 00:31:16 +0200
Subject: [PATCH] Latest modules
---
example/Numerics/Vortex_in_cell/main_vic.cpp | 2 +-
.../Vortex_in_cell/main_vic_petsc.cpp | 75 ++++++++++---------
.../Vortex_in_cell/main_vic_petsc2.cpp | 2 +-
.../Vortex_in_cell/main_vic_petsc3.cpp | 2 +-
openfpm_numerics | 2 +-
5 files changed, 45 insertions(+), 38 deletions(-)
diff --git a/example/Numerics/Vortex_in_cell/main_vic.cpp b/example/Numerics/Vortex_in_cell/main_vic.cpp
index 2605cea41..7481f1cad 100644
--- a/example/Numerics/Vortex_in_cell/main_vic.cpp
+++ b/example/Numerics/Vortex_in_cell/main_vic.cpp
@@ -18,7 +18,7 @@
#include "Solvers/petsc_solver.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "interpolation/mp4_kernel.hpp"
-#include "interpolation/interpolation.hpp"
+#include "../../../openfpm_numerics/src/interpolation/interpolation.hpp"
constexpr int x = 0;
constexpr int y = 1;
diff --git a/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp b/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
index 3d8071b36..ed2a4c73b 100644
--- a/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
+++ b/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
@@ -96,6 +96,7 @@
//! \cond [inclusion] \endcond
+#include "../../../openfpm_numerics/src/interpolation/interpolation.hpp"
#include "Grid/grid_dist_id.hpp"
#include "Vector/vector_dist.hpp"
#include "Matrix/SparseMatrix.hpp"
@@ -103,7 +104,6 @@
#include "FiniteDifference/FDScheme.hpp"
#include "Solvers/petsc_solver.hpp"
#include "interpolation/mp4_kernel.hpp"
-#include "interpolation/interpolation.hpp"
#include "Solvers/petsc_solver_AMG_report.hpp"
constexpr int x = 0;
@@ -293,6 +293,7 @@ const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};
//! \cond [poisson_syseq] \endcond
+
/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
*
* # Step 2: Helmotz-hodge projection # {#vic_hlm_proj}
@@ -406,7 +407,7 @@ const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};
* domain simulation domain
*
*/
-void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
+void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain, petsc_solver<double> & solver, petsc_solver<double>::return_type & x_ , bool init)
{
///////////////////////////////////////////////////////////////
@@ -476,25 +477,29 @@ void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
//! \cond [solve_petsc] \endcond
- // Create a PETSC solver to get the solution x
- petsc_solver<double> solver;
-
- // Set the conjugate-gradient as method to solve the linear solver
- solver.setSolver(KSPBCGS);
+ if (init == true)
+ {
+ // Set the conjugate-gradient as method to solve the linear solver
+ solver.setSolver(KSPBCGS);
- // Set the absolute tolerance to determine that the method is converged
- solver.setAbsTol(0.1);
+ // Set the absolute tolerance to determine that the method is converged
+ solver.setAbsTol(0.1);
- // Set the maximum number of iterations
- solver.setMaxIter(500);
+ // Set the maximum number of iterations
+ solver.setMaxIter(500);
- timer tm_solve;
- tm_solve.start();
+ timer tm_solve;
+ tm_solve.start();
- // Give to the solver A and b, return x, the solution
- auto x_ = solver.solve(fd.getA(),fd.getB());
+ // Give to the solver A and b, return x, the solution
+ solver.solve(fd.getA(),x_,fd.getB());
- tm_solve.stop();
+ tm_solve.stop();
+ }
+ else
+ {
+ solver.solve(x_,fd.getB());
+ }
// copy the solution x to the grid psi
fd.template copy<phi>(x_,psi);
@@ -633,7 +638,7 @@ void remesh(particles_type & vd, grid_type & gr,Box<3,float> & domain)
*
*/
-void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc_solver<double>::return_type (& phi_s)[3])
+void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc_solver<double>::return_type (& phi_s)[3],petsc_solver<double> & solver)
{
//! \cond [poisson_obj_eq] \endcond
@@ -690,16 +695,7 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
// impose the poisson equation using gr_ps = vorticity for the right-hand-side (on the full domain)
fd.template impose_dit<phi>(poisson(),gr_ps,gr_ps.getDomainIterator());
- // Create an PETSC solver
- petsc_solver<double> solver;
-
- solver.setSolver(KSPBCGS);
solver.setAbsTol(0.01);
- solver.setMaxIter(500);
-
- // Get the sparse matrix that represent the left-hand-side
- // of the equation
- SparseMatrix<double,int,PETSC_BASE> & A = fd.getA();
// Get the vector that represent the right-hand-side
Vector<double,PETSC_BASE> & b = fd.getB();
@@ -710,14 +706,14 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
// Give to the solver A and b in this case we are giving
// an intitial guess phi_s calculated in the previous
// time step
- solver.solve(A,phi_s[i],b);
+ solver.solve(phi_s[i],b);
tm_solve.stop();
// Calculate the residual
solError serr;
- serr = solver.get_residual_error(A,phi_s[i],b);
+ serr = solver.get_residual_error(phi_s[i],b);
Vcluster & v_cl = create_vcluster();
if (v_cl.getProcessUnitID() == 0)
@@ -1034,7 +1030,8 @@ template<typename grid, typename vector> void do_step(vector & particles,
grid & g_dvort,
Box<3,float> & domain,
interpolate<particles_type,grid_type,mp4_kernel<float>> & inte,
- petsc_solver<double>::return_type (& phi_s)[3])
+ petsc_solver<double>::return_type (& phi_s)[3],
+ petsc_solver<double> & solver)
{
constexpr int rhs = 0;
@@ -1049,7 +1046,7 @@ template<typename grid, typename vector> void do_step(vector & particles,
//! \cond [step_56] \endcond
// Calculate velocity from vorticity
- comp_vel(domain,g_vort,g_vel,phi_s);
+ comp_vel(domain,g_vort,g_vel,phi_s,solver);
calc_rhs(g_vort,g_vel,g_dvort);
//! \cond [step_56] \endcond
@@ -1163,6 +1160,7 @@ int main(int argc, char* argv[])
// It store the solution to compute velocity
// It is used as initial guess every time we call the solver
Vector<double,PETSC_BASE> phi_s[3];
+ Vector<double,PETSC_BASE> x_;
// Parallel object
Vcluster & v_cl = create_vcluster();
@@ -1174,8 +1172,14 @@ int main(int argc, char* argv[])
// initialize the ring step 1
init_ring(g_vort,domain);
+ x_.resize(g_vort.size(),g_vort.getLocalDomainSize());
+ x_.setZero();
+
+ // Create a PETSC solver to get the solution x
+ petsc_solver<double> solver;
+
// Do the helmotz projection step 2
- helmotz_hodge_projection(g_vort,domain);
+ helmotz_hodge_projection(g_vort,domain,solver,x_,true);
// initialize the particle on a mesh step 3
remesh(particles,g_vort,domain);
@@ -1233,13 +1237,16 @@ int main(int argc, char* argv[])
for ( ; i < 10001 ; i++)
{
// do step 4-5-6-7
- do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s);
+ do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s,solver);
+
+ openfpm_finalize();
+ return 0;
// do step 8
rk_step1(particles);
// do step 9-10-11-12
- do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s);
+ do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s,solver);
// do step 13
rk_step2(particles);
@@ -1250,7 +1257,7 @@ int main(int argc, char* argv[])
g_vort.template ghost_put<add_,vorticity>();
// helmotz-hodge projection
- helmotz_hodge_projection(g_vort,domain);
+ helmotz_hodge_projection(g_vort,domain,solver,x_,false);
remesh(particles,g_vort,domain);
diff --git a/example/Numerics/Vortex_in_cell/main_vic_petsc2.cpp b/example/Numerics/Vortex_in_cell/main_vic_petsc2.cpp
index 22d153abc..9ab52d9a2 100644
--- a/example/Numerics/Vortex_in_cell/main_vic_petsc2.cpp
+++ b/example/Numerics/Vortex_in_cell/main_vic_petsc2.cpp
@@ -18,7 +18,7 @@
#include "Solvers/petsc_solver.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "interpolation/mp4_kernel.hpp"
-#include "interpolation/interpolation.hpp"
+#include "../../../openfpm_numerics/src/interpolation/interpolation.hpp"
constexpr int x = 0;
constexpr int y = 1;
diff --git a/example/Numerics/Vortex_in_cell/main_vic_petsc3.cpp b/example/Numerics/Vortex_in_cell/main_vic_petsc3.cpp
index 56452eab9..13972fb05 100644
--- a/example/Numerics/Vortex_in_cell/main_vic_petsc3.cpp
+++ b/example/Numerics/Vortex_in_cell/main_vic_petsc3.cpp
@@ -18,7 +18,7 @@
#include "Solvers/petsc_solver.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "interpolation/mp4_kernel.hpp"
-#include "interpolation/interpolation.hpp"
+#include "../../../openfpm_numerics/src/interpolation/interpolation.hpp"
constexpr int x = 0;
constexpr int y = 1;
diff --git a/openfpm_numerics b/openfpm_numerics
index 0f84683a5..e8cb410eb 160000
--- a/openfpm_numerics
+++ b/openfpm_numerics
@@ -1 +1 @@
-Subproject commit 0f84683a560d059a1de5fe83554794c21f71a710
+Subproject commit e8cb410eba653c121428df63a440da58920d5f27
--
GitLab