Commit 914d6ba0 authored by incardon's avatar incardon
Browse files

Testing solvers CG

parent 3b4db8ae
......@@ -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);
......
......@@ -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();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment