110 #include "interpolation/interpolation.hpp"
111 #include "Grid/grid_dist_id.hpp"
112 #include "Vector/vector_dist.hpp"
113 #include "Matrix/SparseMatrix.hpp"
114 #include "Vector/Vector.hpp"
115 #include "FiniteDifference/FDScheme.hpp"
116 #include "Solvers/petsc_solver.hpp"
117 #include "interpolation/mp4_kernel.hpp"
118 #include "Solvers/petsc_solver_AMG_report.hpp"
133 float sigma = 1.0/3.523;
136 float tgtre = 3000.0;
139 float nu = 1.0/tgtre;
146 const unsigned int nsteps = 10;
148 const unsigned int nsteps = 10001;
152 constexpr
unsigned int vorticity = 0;
153 constexpr
unsigned int velocity = 0;
154 constexpr
unsigned int p_vel = 1;
155 constexpr
int rhs_part = 2;
156 constexpr
unsigned int old_vort = 3;
157 constexpr
unsigned int old_pos = 4;
161 template<
typename gr
id>
void calc_and_print_max_div_and_int(
grid & g_vort)
165 g_vort.template ghost_get<vorticity>();
166 auto it5 = g_vort.getDomainIterator();
168 double max_vort = 0.0;
170 double int_vort[3] = {0.0,0.0,0.0};
174 auto key = it5.get();
178 tmp = 1.0f/g_vort.spacing(x)/2.0f*(g_vort.template get<vorticity>(
key.move(x,1))[x] - g_vort.template get<vorticity>(
key.move(x,-1))[x] ) +
179 1.0f/g_vort.spacing(y)/2.0f*(g_vort.template get<vorticity>(
key.move(y,1))[y] - g_vort.template get<vorticity>(
key.move(y,-1))[y] ) +
180 1.0f/g_vort.spacing(z)/2.0f*(g_vort.template get<vorticity>(
key.move(z,1))[z] - g_vort.template get<vorticity>(
key.move(z,-1))[z] );
182 int_vort[x] += g_vort.template get<vorticity>(
key)[x];
183 int_vort[y] += g_vort.template get<vorticity>(
key)[y];
184 int_vort[z] += g_vort.template get<vorticity>(
key)[z];
192 Vcluster & v_cl = create_vcluster();
194 v_cl.
sum(int_vort[0]);
195 v_cl.
sum(int_vort[1]);
196 v_cl.
sum(int_vort[2]);
200 {std::cout <<
"Max div for vorticity " << max_vort <<
" Integral: " << int_vort[0] <<
" " << int_vort[1] <<
" " << int_vort[2] << std::endl;}
239 constexpr
int nk = 32;
244 for (
size_t i = 0 ; i < nk ; i++)
246 ak[i] = rand()/RAND_MAX;
247 bk[i] = rand()/RAND_MAX;
251 float gamma = nu * tgtre;
252 float rinv2 = 1.0f/(sigma*sigma);
253 float max_vorticity = gamma*rinv2/M_PI;
260 auto key_d = it.get();
261 auto key = it.getGKey(key_d);
266 float theta1 = atan2((ty-domain.
getHigh(1)/2.0f),(tz-domain.
getHigh(2)/2.0f));
269 float rad1r = sqrt((ty-domain.
getHigh(1)/2.0f)*(ty-domain.
getHigh(1)/2.0f) + (tz-domain.
getHigh(2)/2.0f)*(tz-domain.
getHigh(2)/2.0f)) - ringr1;
270 float rad1t = tx - 1.0f;
271 float rad1sq = rad1r*rad1r + rad1t*rad1t;
272 float radstr = -exp(-rad1sq*rinv2)*rinv2*gamma/M_PI;
273 gr.template get<vorticity>(key_d)[x] = 0.0f;
274 gr.template get<vorticity>(key_d)[y] = -radstr * cos(theta1);
275 gr.template get<vorticity>(key_d)[z] = radstr * sin(theta1);
279 float rad1r_ = sqrt((ty-domain.
getHigh(1)/2.0f)*(ty-domain.
getHigh(1)/2.0f) + (tz-domain.
getHigh(2)/2.0f)*(tz-domain.
getHigh(2)/2.0f)) + ringr1;
280 float rad1sqTILDA = rad1sq*rinv2;
281 radstr = exp(-rad1sq*rinv2)*rinv2*gamma/M_PI;
282 gr.template get<vorticity>(key_d)[x] = 0.0f;
283 gr.template get<vorticity>(key_d)[y] = -radstr * cos(theta1);
284 gr.template get<vorticity>(key_d)[z] = radstr * sin(theta1);
303 static const unsigned int dims = 3;
305 static const unsigned int nvar = 1;
317 static const int grid_type = NORMAL_GRID;
446 constexpr
unsigned int phi = 0;
462 gr.template ghost_get<vorticity>();
477 psi.template get<phi>(
key) = 1.0f/gr.
spacing(x)/2.0f*(gr.template get<vorticity>(
key.move(x,1))[x] - gr.template get<vorticity>(
key.move(x,-1))[x] ) +
478 1.0f/gr.
spacing(y)/2.0f*(gr.template get<vorticity>(
key.move(y,1))[y] - gr.template get<vorticity>(
key.move(y,-1))[y] ) +
479 1.0f/gr.
spacing(z)/2.0f*(gr.template get<vorticity>(
key.move(z,1))[z] - gr.template get<vorticity>(
key.move(z,-1))[z] );
487 calc_and_print_max_div_and_int(gr);
503 fd.template impose_dit<0>(poisson(),psi,psi.getDomainIterator());
537 solver.
solve(fd.getA(),x_,fd.getB());
544 solver.
solve(x_,fd.getB());
549 fd.template copy<phi>(x_,psi);
555 psi.template ghost_get<phi>();
563 auto key = it2.get();
565 gr.template get<vorticity>(
key)[x] -= 1.0f/2.0f/psi.
spacing(x)*(psi.template get<phi>(
key.move(x,1))-psi.template get<phi>(
key.move(x,-1)));
566 gr.template get<vorticity>(
key)[y] -= 1.0f/2.0f/psi.
spacing(y)*(psi.template get<phi>(
key.move(y,1))-psi.template get<phi>(
key.move(y,-1)));
567 gr.template get<vorticity>(
key)[z] -= 1.0f/2.0f/psi.
spacing(z)*(psi.template get<phi>(
key.move(z,1))-psi.template get<phi>(
key.move(z,-1)));
574 calc_and_print_max_div_and_int(gr);
608 auto key = git.get_dist();
619 vd.template getLastProp<vorticity>()[x] = gr.template get<vorticity>(
key)[x];
620 vd.template getLastProp<vorticity>()[y] = gr.template get<vorticity>(
key)[y];
621 vd.template getLastProp<vorticity>()[z] = gr.template get<vorticity>(
key)[z];
691 constexpr
unsigned int phi = 0;
713 calc_and_print_max_div_and_int(g_vort);
716 for (
size_t i = 0 ; i < 3 ; i++)
721 auto it2 = gr_ps.getDomainIterator();
726 auto key = it2.get();
729 gr_ps.get<vorticity>(
key) = g_vort.template get<vorticity>(
key)[i];
741 fd.template impose_dit<phi>(poisson(),gr_ps,gr_ps.getDomainIterator());
754 solver.
solve(phi_s[i],b);
763 Vcluster & v_cl = create_vcluster();
765 {std::cout <<
"Solved component " << i <<
" Error: " << serr.
err_inf << std::endl;}
768 fd.template copy<phi>(phi_s[i],gr_ps);
774 auto it3 = gr_ps.getDomainIterator();
779 auto key = it3.get();
781 phi_v.get<velocity>(
key)[i] = gr_ps.get<phi>(
key);
791 phi_v.ghost_get<phi>();
793 float inv_dx = 0.5f / g_vort.
spacing(0);
794 float inv_dy = 0.5f / g_vort.
spacing(1);
795 float inv_dz = 0.5f / g_vort.
spacing(2);
797 auto it3 = phi_v.getDomainIterator();
802 auto key = it3.get();
804 float phi_xy = (phi_v.get<phi>(
key.move(y,1))[x] - phi_v.get<phi>(
key.move(y,-1))[x])*inv_dy;
805 float phi_xz = (phi_v.get<phi>(
key.move(z,1))[x] - phi_v.get<phi>(
key.move(z,-1))[x])*inv_dz;
806 float phi_yx = (phi_v.get<phi>(
key.move(x,1))[y] - phi_v.get<phi>(
key.move(x,-1))[y])*inv_dx;
807 float phi_yz = (phi_v.get<phi>(
key.move(z,1))[y] - phi_v.get<phi>(
key.move(z,-1))[y])*inv_dz;
808 float phi_zx = (phi_v.get<phi>(
key.move(x,1))[z] - phi_v.get<phi>(
key.move(x,-1))[z])*inv_dx;
809 float phi_zy = (phi_v.get<phi>(
key.move(y,1))[z] - phi_v.get<phi>(
key.move(y,-1))[z])*inv_dy;
811 g_vel.template get<velocity>(
key)[x] = phi_zy - phi_yz;
812 g_vel.template get<velocity>(
key)[y] = phi_xz - phi_zx;
813 g_vel.template get<velocity>(
key)[z] = phi_yx - phi_xy;
831 gr.template get<prp>(
key)[0] = 0.0;
832 gr.template get<prp>(
key)[1] = 0.0;
833 gr.template get<prp>(
key)[2] = 0.0;
849 vd.template getProp<prp>(
key)[0] = 0.0;
850 vd.template getProp<prp>(
key)[1] = 0.0;
851 vd.template getProp<prp>(
key)[2] = 0.0;
873 template<
typename gr
id>
void calc_rhs(
grid & g_vort,
grid & g_vel,
grid & g_dwp)
876 constexpr
int rhs = 0;
880 float fac1 = 2.0f*nu/(g_vort.spacing(0)*g_vort.spacing(0));
881 float fac2 = 2.0f*nu/(g_vort.spacing(1)*g_vort.spacing(1));
882 float fac3 = 2.0f*nu/(g_vort.spacing(2)*g_vort.spacing(2));
884 float fac4 = 0.5f/(g_vort.spacing(0));
885 float fac5 = 0.5f/(g_vort.spacing(1));
886 float fac6 = 0.5f/(g_vort.spacing(2));
888 auto it = g_dwp.getDomainIterator();
890 g_vort.template ghost_get<vorticity>();
896 g_dwp.template get<rhs>(
key)[x] = fac1*(g_vort.template get<vorticity>(
key.move(x,1))[x]+g_vort.template get<vorticity>(
key.move(x,-1))[x])+
897 fac2*(g_vort.template get<vorticity>(
key.move(y,1))[x]+g_vort.template get<vorticity>(
key.move(y,-1))[x])+
898 fac3*(g_vort.template get<vorticity>(
key.move(z,1))[x]+g_vort.template get<vorticity>(
key.move(z,-1))[x])-
899 2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(
key)[x]+
900 fac4*g_vort.template get<vorticity>(
key)[x]*
901 (g_vel.template get<velocity>(
key.move(x,1))[x] - g_vel.template get<velocity>(
key.move(x,-1))[x])+
902 fac5*g_vort.template get<vorticity>(
key)[y]*
903 (g_vel.template get<velocity>(
key.move(y,1))[x] - g_vel.template get<velocity>(
key.move(y,-1))[x])+
904 fac6*g_vort.template get<vorticity>(
key)[z]*
905 (g_vel.template get<velocity>(
key.move(z,1))[x] - g_vel.template get<velocity>(
key.move(z,-1))[x]);
907 g_dwp.template get<rhs>(
key)[y] = fac1*(g_vort.template get<vorticity>(
key.move(x,1))[y]+g_vort.template get<vorticity>(
key.move(x,-1))[y])+
908 fac2*(g_vort.template get<vorticity>(
key.move(y,1))[y]+g_vort.template get<vorticity>(
key.move(y,-1))[y])+
909 fac3*(g_vort.template get<vorticity>(
key.move(z,1))[y]+g_vort.template get<vorticity>(
key.move(z,-1))[y])-
910 2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(
key)[y]+
911 fac4*g_vort.template get<vorticity>(
key)[x]*
912 (g_vel.template get<velocity>(
key.move(x,1))[y] - g_vel.template get<velocity>(
key.move(x,-1))[y])+
913 fac5*g_vort.template get<vorticity>(
key)[y]*
914 (g_vel.template get<velocity>(
key.move(y,1))[y] - g_vel.template get<velocity>(
key.move(y,-1))[y])+
915 fac6*g_vort.template get<vorticity>(
key)[z]*
916 (g_vel.template get<velocity>(
key.move(z,1))[y] - g_vel.template get<velocity>(
key.move(z,-1))[y]);
919 g_dwp.template get<rhs>(
key)[z] = fac1*(g_vort.template get<vorticity>(
key.move(x,1))[z]+g_vort.template get<vorticity>(
key.move(x,-1))[z])+
920 fac2*(g_vort.template get<vorticity>(
key.move(y,1))[z]+g_vort.template get<vorticity>(
key.move(y,-1))[z])+
921 fac3*(g_vort.template get<vorticity>(
key.move(z,1))[z]+g_vort.template get<vorticity>(
key.move(z,-1))[z])-
922 2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(
key)[z]+
923 fac4*g_vort.template get<vorticity>(
key)[x]*
924 (g_vel.template get<velocity>(
key.move(x,1))[z] - g_vel.template get<velocity>(
key.move(x,-1))[z])+
925 fac5*g_vort.template get<vorticity>(
key)[y]*
926 (g_vel.template get<velocity>(
key.move(y,1))[z] - g_vel.template get<velocity>(
key.move(y,-1))[z])+
927 fac6*g_vort.template get<vorticity>(
key)[z]*
928 (g_vel.template get<velocity>(
key.move(z,1))[z] - g_vel.template get<velocity>(
key.move(z,-1))[z]);
980 auto key = it2.get();
1036 while (it2.isNext())
1038 auto key = it2.get();
1074 template<
typename gr
id,
typename vector>
void do_step(vector & particles,
1083 constexpr
int rhs = 0;
1088 inte.template p2m<vorticity,vorticity>(particles,g_vort);
1090 g_vort.template ghost_put<add_,vorticity>();
1097 comp_vel(domain,g_vort,g_vel,phi_s,solver);
1098 calc_rhs(g_vort,g_vel,g_dvort);
1104 g_dvort.template ghost_get<rhs>();
1105 g_vel.template ghost_get<velocity>();
1108 inte.template m2p<rhs,rhs_part>(g_dvort,particles);
1109 inte.template m2p<velocity,p_vel>(g_vel,particles);
1114 template<
typename vector,
typename gr
id>
void check_point_and_save(vector & particles,
1120 Vcluster & v_cl = create_vcluster();
1124 particles.write(
"part_out_" + std::to_string(i),VTK_WRITER | FORMAT_BINARY);
1125 g_vort.template ghost_get<vorticity>();
1126 g_vel.template ghost_get<velocity>();
1127 g_vel.write_frame(
"grid_velocity",i, VTK_WRITER | FORMAT_BINARY);
1128 g_vort.write_frame(
"grid_vorticity",i, VTK_WRITER | FORMAT_BINARY);
1137 while (it_s.isNext())
1139 auto p = it_s.get();
1141 float vort_magn = sqrt(particles.template getProp<vorticity>(p)[0] * particles.template getProp<vorticity>(p)[0] +
1142 particles.template getProp<vorticity>(p)[1] * particles.template getProp<vorticity>(p)[1] +
1143 particles.template getProp<vorticity>(p)[2] * particles.template getProp<vorticity>(p)[2]);
1145 if (vort_magn < 0.1)
1153 part_save.getLastPos()[0] = particles.getPos(p)[0];
1154 part_save.getLastPos()[1] = particles.getPos(p)[1];
1155 part_save.getLastPos()[2] = particles.getPos(p)[2];
1157 part_save.template getLastProp<0>() = vort_magn;
1165 particles.save(
"check_point");
1186 int main(
int argc,
char* argv[])
1189 openfpm_init(&argc,&argv);
1202 long int sz[] = {96,96,96};
1203 size_t szu[] = {(size_t)sz[0],(
size_t)sz[1],(size_t)sz[2]};
1208 grid_type g_vel(g_vort.getDecomposition(),szu,g);
1209 grid_type g_dvort(g_vort.getDecomposition(),szu,g);
1218 Vcluster & v_cl = create_vcluster();
1222 {std::cout <<
"SPACING " << g_vort.spacing(0) <<
" " << g_vort.spacing(1) <<
" " << g_vort.spacing(2) << std::endl;}
1225 init_ring(g_vort,domain);
1227 x_.
resize(g_vort.size(),g_vort.getLocalDomainSize());
1234 helmotz_hodge_projection(g_vort,domain,solver,x_,
true);
1237 remesh(particles,g_vort,domain);
1240 size_t tot_part = particles.size_local();
1245 phi_s[0].resize(g_vort.size(),g_vort.getLocalDomainSize());
1246 phi_s[1].resize(g_vort.size(),g_vort.getLocalDomainSize());
1247 phi_s[2].resize(g_vort.size(),g_vort.getLocalDomainSize());
1263 {g_vort.write(
"grid_vorticity_init", VTK_WRITER | FORMAT_BINARY);}
1274 std::string restart(argv[1]);
1277 i = std::stoi(restart);
1280 particles.load(
"check_point_" + std::to_string(i));
1285 {std::cout <<
"Restarting from " << i << std::endl;}
1289 for ( ; i < nsteps ; i++)
1292 do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s,solver);
1295 rk_step1(particles);
1298 do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s,solver);
1301 rk_step2(particles);
1305 inte.template p2m<vorticity,vorticity>(particles,g_vort);
1306 g_vort.template ghost_put<add_,vorticity>();
1309 helmotz_hodge_projection(g_vort,domain,solver,x_,
false);
1311 remesh(particles,g_vort,domain);
1315 {std::cout <<
"Step " << i << std::endl;}
1318 if (i % 100 == 0) {check_point_and_save(particles,g_vort,g_vel,g_dvort,i);}
void sum(T &num)
Sum the numbers across all processors and get the result.
auto get(const grid_dist_key_dx< dim > &v1) const -> typename std::add_lvalue_reference< decltype(loc_grid.get(v1.getSub()).template get< p >(v1.getKey()))>::type
Get the reference of the selected element.
grid_dist_iterator< dim, device_grid, FREE > getDomainIterator() const
It return an iterator that span the full grid domain (each processor span its local domain) ...
void setPreconditioner(PCType type)
Set the preconditioner of the linear solver.
auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get< id >(vec_key.getKey()))
Get the property of an element.
T getLow(int i) const
get the i-coordinate of the low bound interval of the box
size_t getProcessUnitID()
Get the process unit id.
void execute()
Execute all the requests.
const grid_sm< dim, T > & getGridInfo() const
Get an object containing the grid informations.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support...
static const unsigned int dims
Number of dimansions of the problem.
auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get< 0 >(vec_key.getKey()))
Get the position of an element.
void setMaxIter(PetscInt n)
Set the maximum number of iteration for Krylov solvers.
In case T does not match the PETSC precision compilation create a stub structure. ...
SparseMatrix< double, int, PETSC_BASE > SparseMatrix_type
Sparse matrix used to sove the linear system (we use PETSC)
T getHigh(int i) const
get the high interval of the box
Meta-function to return a compatible zero-element.
It contain statistic of the error of the calculated solution.
void setPreconditionerAMG_cycleType(const std::string &cycle_type, int sweep_up=-1, int sweep_dw=-1, int sweep_crs=-1)
It set the type of cycle and optionally the number of sweep.
float stype
type of the spatial coordinates
void setAbsTol(PetscReal abstol_)
Set the absolute tolerance as stop criteria.
Vector< double, PETSC_BASE > solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, bool initial_guess=false)
Here we invert the matrix and solve the system.
void max(T &num)
Get the maximum number across all processors (or reduction with infinity norm)
grid_dist_id_iterator_dec< Decomposition > getGridIterator(const size_t(&sz)[dim])
Vector< double, PETSC_BASE > Vector_type
Vector to solve the system (PETSC)
Implementation of VCluster class.
void setPreconditionerAMG_maxit(int nit)
Set the maximum number of V or W cycle for algebraic-multi-grid.
Sparse Matrix implementation.
This is a distributed grid.
void resize(size_t row, size_t row_n)
stub resize
static const unsigned int nvar
We have only one scalar unknown.
const size_t(& getSize() const)[N]
Return the size of the grid as an array.
void setPreconditionerAMG_nl(int nl)
Set the number of levels for the algebraic-multigrid preconditioner.
This class is able to do Matrix inversion in parallel with PETSC solvers.
void clear()
remove all the elements
void setSolver(KSPType type)
Set the Petsc solver.
void start()
Start the timer.
Laplacian second order on h (spacing)
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor...
void add()
Add local particle.
This class is a trick to indicate the compiler a specific specialization pattern. ...
grid_dist_id< 3, float, aggregate< float > > b_grid
grid that store
grid_dist_iterator< dim, device_grid, FIXED > getDomainGhostIterator() const
It return an iterator that span the grid domain + ghost part.
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.
void setPreconditionerAMG_coarsen(const std::string &type)
Set the coarsening method for the algebraic-multi-grid preconditioner.
void setPreconditionerAMG_relax(const std::string &type, int k=REL_ALL)
Set the relaxation method for the algebraic-multi-grid preconditioner.
auto getLastPos() -> decltype(v_pos.template get< 0 >(0))
Get the position of the last element.
static const bool boundary[]
specify the boundary conditions
solError get_residual_error(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &x, const Vector< double, PETSC_BASE > &b)
It return the resiual error.
size_t getProcessingUnits()
Get the total number of processors.
Decomposition & getDecomposition()
Get the object that store the information about the decomposition.
Class for cpu time benchmarking.
void setPreconditionerAMG_coarsenNodalType(int norm)
Set the block coarsening norm type.
Main class for interpolation Particle to mest p2m and Mesh to particle m2p.
St spacing(size_t i) const
Get the spacing of the grid in direction i.
void stop()
Stop the timer.
PetscReal err_inf
infinity norm of the error