From 914d6ba0e7eb27803572b0e8b6665ff6b4bdd500 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Sat, 17 Jun 2017 16:17:49 +0200
Subject: [PATCH] Testing solvers CG
---
example/Numerics/Vortex_in_cell/CG.hpp | 6 ++-
.../Vortex_in_cell/main_vic_petsc.cpp | 49 ++++++++++++++-----
2 files changed, 40 insertions(+), 15 deletions(-)
diff --git a/example/Numerics/Vortex_in_cell/CG.hpp b/example/Numerics/Vortex_in_cell/CG.hpp
index 98eedd529..3868718b6 100644
--- a/example/Numerics/Vortex_in_cell/CG.hpp
+++ b/example/Numerics/Vortex_in_cell/CG.hpp
@@ -145,6 +145,8 @@ void CG(void (* A)(grid_type_scal & in, grid_type_scal & out), grid_type_scal &
A(x,Ap);
x_plus_alpha_p(r,Ap,-1);
+ std::cout << "It init: " << " " << sqrt(norm(r)) << std::endl;
+
// r0 = p0
copy(r,p);
@@ -157,9 +159,9 @@ void CG(void (* A)(grid_type_scal & in, grid_type_scal & out), grid_type_scal &
float r_norm = norm(rn);
- std::cout << "It: " << i << " " << r_norm << std::endl;
+ std::cout << "It: " << i << " " << sqrt(r_norm) << std::endl;
- if (r_norm < 0.1)
+ if (sqrt(r_norm) < 0.1)
return;
float beta_c = beta(r,rn);
diff --git a/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp b/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
index 0179103a0..9978ef3e2 100644
--- a/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
+++ b/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
@@ -133,6 +133,10 @@ float nu = 0.0001535;
float dt = 0.006;
+#ifdef SE_CLASS1
+DIOCANE
+#endif
+
// All the properties by index
constexpr unsigned int vorticity = 0;
constexpr unsigned int velocity = 0;
@@ -516,14 +520,21 @@ void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
solver.setSolver(KSPCG);
// Set the absolute tolerance to determine that the method is converged
- solver.setAbsTol(0.001);
+ solver.setAbsTol(0.1);
// Set the maximum number of iterations
solver.setMaxIter(500);
+ solver.log_monitor();
+
+ 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());
+ tm_solve.stop();
+
//////////////// CG Call //////////////////////////
// Here we create a distributed grid to store the result of the helmotz projection
@@ -540,8 +551,16 @@ void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
++zit;
}
+ timer tm_solve2;
+
+ tm_solve2.start();
+
CG(OpA,sol,psi);
+ tm_solve2.stop();
+
+ std::cout << "Time to solve: " << tm_solve2.getwct() << " " << tm_solve.getwct() << std::endl;
+
///////////////////////////////////////////////////
// copy the solution x to the grid psi
@@ -753,11 +772,16 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
// Get the vector that represent the right-hand-side
Vector<double,PETSC_BASE> & b = fd.getB();
+ timer tm_solve;
+ tm_solve.start();
+
// 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);
+ tm_solve.stop();
+
//////////////// CG Call //////////////////////////
// Here we create a distributed grid to store the result of the helmotz projection
@@ -774,8 +798,16 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
++zit;
}
+ timer tm_solve2;
+ tm_solve2.start();
+
CG(OpA,sol,gr_ps);
+ tm_solve2.stop();
+
+ std::cout << "Time to solve: " << tm_solve2.getwct() << " " << tm_solve.getwct() << std::endl;
+
+
///////////////////////////////////////////////////
// Calculate the residual
@@ -790,13 +822,6 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
// copy the solution to grid
fd.template copy<phi>(phi_s[i],gr_ps);
- /////////////// DEBUG ///////////////
-
- gr_ps.write("Solution_petsc_" + std::to_string(i));
- sol.write("Solution_my_" + std::to_string(i));
-
- /////////////////////////////////////
-
//! \cond [solve_poisson_comp] \endcond
//! \cond [copy_to_phi_v] \endcond
@@ -816,8 +841,6 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
//! \cond [copy_to_phi_v] \endcond
}
- exit(0);
-
//! \cond [curl_phi_v] \endcond
phi_v.ghost_get<phi>();
@@ -1303,7 +1326,7 @@ int main(int argc, char* argv[])
}
// Time Integration
- for ( ; i < 10001 ; i++)
+ for ( ; i < 1 ; i++)
{
// do step 4-5-6-7
do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s);
@@ -1312,7 +1335,7 @@ int main(int argc, char* argv[])
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);
// do step 13
rk_step2(particles);
@@ -1330,7 +1353,7 @@ int main(int argc, char* argv[])
// every 100 steps write the output
if (i % 100 == 0) {check_point_and_save(particles,g_vort,g_vel,g_dvort,i);}
-
+*/
}
openfpm_finalize();
--
GitLab