1 #include "Grid/grid_dist_id.hpp"
2 #include "data_type/aggregate.hpp"
54 extern "C" void update_new(
const int* lo,
const int* hi,
55 double* u,
const int* ulo,
const int* uhi,
56 double* v,
const int* vlo,
const int* vhi,
57 double* flu,
const int* fulo,
const int* fuhi,
58 double* flv,
const int* fvlo,
const int* fvhi,
59 const double * dt,
const double * uFactor,
const double * vFactor,
const double *
F,
71 auto it = OldU.getDomainIterator();
89 long int x_start = OldU.size(0)*1.55f/domain.
getHigh(0);
90 long int y_start = OldU.size(1)*1.55f/domain.
getHigh(1);
91 long int z_start = OldU.size(1)*1.55f/domain.
getHigh(2);
93 long int x_stop = OldU.size(0)*1.85f/domain.
getHigh(0);
94 long int y_stop = OldU.size(1)*1.85f/domain.
getHigh(1);
95 long int z_stop = OldU.size(1)*1.85f/domain.
getHigh(2);
99 auto it_init = OldU.getSubDomainIterator(start,stop);
101 while (it_init.isNext())
105 OldU.get(
key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0;
106 OldV.get(
key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0;
120 double uFactor_s,
double vFactor_s,
double deltaT,
double F,
double K)
122 #ifndef FORTRAN_UPDATE
126 Vc::double uFactor = uFactor_s;
127 Vc::double vFactor = vFactor_s;
129 WHILE_M(OldU,star_stencil_3D)
130 auto & U_old = GET_GRID_M(OldU);
131 auto & V_old = GET_GRID_M(OldV);
133 auto & U_new = GET_GRID_M(NewU);
134 auto & V_new = GET_GRID_M(NewV);
135 ITERATE_3D_M(Vc::double_v::Size)
138 auto Cp = it.getStencil<0>();
141 auto mx = it.getStencil<1>();
142 auto px = it.getStencil<2>();
143 auto my = it.getStencil<3>();
144 auto py = it.getStencil<4>();
145 auto mz = it.getStencil<5>();
146 auto pz = it.getStencil<6>();
149 Vc::double_v u_c(&U_old.get<0>(Cp),Vc::Unaligned);
150 Vc::double_v u_mz(&U_old.get<0>(mz),Vc::Unaligned);
151 Vc::double_v u_pz(&U_old.get<0>(pz),Vc::Unaligned);
152 Vc::double_v u_my(&U_old.get<0>(my),Vc::Unaligned);
153 Vc::double_v u_py(&U_old.get<0>(py),Vc::Unaligned);
154 Vc::double_v u_mx(&U_old.get<0>(mx),Vc::Unaligned);
155 Vc::double_v u_px(&U_old.get<0>(px),Vc::Unaligned);
158 Vc::double_v v_c(&V_old.get<0>(Cp),Vc::Unaligned);
159 Vc::double_v v_mz(&V_old.get<0>(mz),Vc::Unaligned);
160 Vc::double_v v_pz(&V_old.get<0>(pz),Vc::Unaligned);
161 Vc::double_v v_my(&V_old.get<0>(my),Vc::Unaligned);
162 Vc::double_v v_py(&V_old.get<0>(py),Vc::Unaligned);
163 Vc::double_v v_mx(&V_old.get<0>(mx),Vc::Unaligned);
164 Vc::double_v v_px(&V_old.get<0>(px),Vc::Unaligned);
166 Vc::double_v out1 = u_c + uFactor * (u_mz + u_pz +
170 - deltaT * u_c * v_c * v_c
171 - deltaT * F * (u_c - 1.0);
173 Vc::double_v out2 = v_c + vFactor * (v_mz + v_pz +
177 deltaT * u_c * v_c * v_c +
178 - deltaT * (F+K) * v_c;
180 out1.store(&U_new.get<0>(Cp),Vc::Unaligned);
181 out2.store(&V_new.get<0>(Cp),Vc::Unaligned);
182 END_LOOP_M(Vc::double_v::Size)
190 auto & ginfo = OldU.getLocalGridsInfo();
192 for (
size_t i = 0 ; i < OldU.getN_loc_grid() ; i++)
194 auto & U_old = OldU.get_loc_grid(i);
195 auto & V_old = OldV.get_loc_grid(i);
197 auto & U_new = NewU.get_loc_grid(i);
198 auto & V_new = NewV.get_loc_grid(i);
200 int lo[3] = {(int)ginfo.get(i).Dbox.getLow(0),(int)ginfo.get(i).Dbox.getLow(1),(int)ginfo.get(i).Dbox.getLow(2)};
201 int hi[3] = {(int)ginfo.get(i).Dbox.getHigh(0),(int)ginfo.get(i).Dbox.getHigh(1),(int)ginfo.get(i).Dbox.getHigh(2)};
203 int ulo[3] = {0,0,0};
204 int uhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
205 int nulo[3] = {0,0,0};
206 int nuhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
208 int vlo[3] = {0,0,0};
209 int vhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
210 int nvlo[3] = {0,0,0};
211 int nvhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
214 (
double *)U_old.getPointer(),ulo,uhi,
215 (
double *)V_old.getPointer(),vlo,vhi,
216 (
double *)U_new.getPointer(),nulo,nuhi,
217 (
double *)V_new.getPointer(),nulo,nvhi,
218 &deltaT, &uFactor_s, &vFactor_s,&F,&K);
228 int main(
int argc,
char* argv[])
230 openfpm_init(&argc,&argv);
236 size_t sz[3] = {256,256,256};
254 size_t timeSteps = 5000;
284 double spacing[3] = {OldU.
spacing(0),OldU.spacing(1),OldU.spacing(2)};
286 init(OldU,OldV,NewU,NewV,domain);
292 OldU.template ghost_get<0>();
293 OldV.template ghost_get<0>();
297 double uFactor = deltaT * du/(spacing[x]*spacing[x]);
298 double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
311 for (
size_t i = 0; i < timeSteps; ++i)
314 std::cout <<
"STEP: " << i << std::endl;
351 step(OldU,OldV,NewU,NewV,star_stencil_3D,uFactor,vFactor,deltaT,F,K);
358 step(NewU,NewV,OldU,OldV,star_stencil_3D,uFactor,vFactor,deltaT,F,K);
380 OldU.save(
"output_u_" + std::to_string(count));
381 OldV.save(
"output_v_" + std::to_string(count));
390 if (create_vcluster().rank() == 0)
391 {std::cout <<
"Total simulation: " << tot_sim.
getwct() << std::endl;}
grid_key_dx is the key to access any element in the grid
T getHigh(int i) const
get the high interval of the box
double getwct()
Return the elapsed real time.
mem_id get(size_t i) const
Get the i index.
This is a distributed grid.
void start()
Start the timer.
This class represent an N-dimensional box.
This class is a trick to indicate the compiler a specific specialization pattern. ...
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Class for cpu time benchmarking.
St spacing(size_t i) const
Get the spacing of the grid in direction i.
void stop()
Stop the timer.
[v_transform metafunction]