quadgrid 0.1
simple cartesian quad grid with particles for c++/octave
Loading...
Searching...
No Matches
taylor_dispersion.cpp File Reference
#include "taylor_dispersion.h"
#include <thrust/extrema.h>
#include <thrust/iterator/transform_iterator.h>

Go to the source code of this file.

Macros

#define gpuGetDeviceCount   cudaGetDeviceCount
#define gpuSetDevice   cudaSetDevice
#define gpuGetDevice   cudaGetDevice

Functions

int main ()

Macro Definition Documentation

◆ gpuGetDevice

#define gpuGetDevice   cudaGetDevice

Definition at line 14 of file taylor_dispersion.cpp.

◆ gpuGetDeviceCount

#define gpuGetDeviceCount   cudaGetDeviceCount

Definition at line 12 of file taylor_dispersion.cpp.

◆ gpuSetDevice

#define gpuSetDevice   cudaSetDevice

Definition at line 13 of file taylor_dispersion.cpp.

Function Documentation

◆ main()

int main ( )

Definition at line 18 of file taylor_dispersion.cpp.

18 {
19
20#ifndef THRUST_CPU
21 int num_gpus;
22 auto err = gpuGetDeviceCount (&num_gpus); if (err) return err;
23 std::cerr << "num_gpus=" << num_gpus <<std::endl;
24 int device;
25 err = gpuSetDevice (num_gpus - 1); if (err) return err;
26 err = gpuGetDevice (&device); if (err) return err;
27 std::cerr << "running on gpu n. " << device << std::endl;
28#endif
29
30 using idx_t = particles_t::idx_t;
31
32 cdf::timer::timer_t timer;
33
34 // read data from file
35 constexpr auto filename = "td.json";
36
37 nlohmann::json j;
38 std::ifstream inbuf (filename);
39 inbuf >> j;
40
41 quadgrid_t<vector_t<real_t>> qg(j["grid_properties"]);
42 particles_t p(j,qg);
43
44 p.build_mass();
45 p.init_particle_mesh();
46
47 std::map<std::string, vector_t<real_t>> vars=
48 j["grid_vars"].get<std::map<std::string, vector_t<real_t>>> ();
49
50 p.memcpy_host_to_device();
51
52 for (const auto & g : vars)
53 p.device_grid_vars[g.first] = g.second;
54
55 inbuf.close ();
56
57
58 p2g_step1 Jdrift_x(p.device_x.cbegin(), p.device_y.cbegin(), thrust::raw_pointer_cast(p.device_grid_M.data()),
59 thrust::raw_pointer_cast(p.device_grid_vars["Jdrift_x"].data()), p.device_ptcl_to_grd.cbegin(), qg.num_rows(), qg.hx(),
60 qg.hy(), p.device_dprops["BETAx"].cbegin(), p.device_dprops["M"].cbegin(), false);
61
62 p2g_step1 Jdrift_y(p.device_x.cbegin(), p.device_y.cbegin(), thrust::raw_pointer_cast(p.device_grid_M.data()),
63 thrust::raw_pointer_cast(p.device_grid_vars["Jdrift_y"].data()), p.device_ptcl_to_grd.cbegin(), qg.num_rows(), qg.hx(),
64 qg.hy(), p.device_dprops["BETAy"].cbegin(), p.device_dprops["M"].cbegin(), false);
65
66 p2gd_step2 Jdiff_x(p.device_x.cbegin(), p.device_y.cbegin(), thrust::raw_pointer_cast(p.device_grid_M.data()),
67 thrust::raw_pointer_cast(p.device_grid_vars["Jdiff_x"].data()), p.device_ptcl_to_grd.cbegin(), qg.num_rows(), qg.hx(),
68 qg.hy(), 1., p.device_dprops["M"].cbegin(), p.device_dprops["zero"].cbegin(), false);
69
70 p2gd_step2 Jdiff_y(p.device_x.cbegin(), p.device_y.cbegin(), thrust::raw_pointer_cast(p.device_grid_M.data()),
71 thrust::raw_pointer_cast(p.device_grid_vars["Jdiff_y"].data()), p.device_ptcl_to_grd.cbegin(), qg.num_rows(), qg.hx(),
72 qg.hy(), 1., p.device_dprops["zero"].cbegin(), p.device_dprops["M"].cbegin(), false);
73
74 g2p_step3 VX(thrust::raw_pointer_cast(p.device_x.data()), thrust::raw_pointer_cast(p.device_y.data()), p.device_grid_M.cbegin(),
75 p.device_grid_vars["Jdrift_x"].cbegin(), p.device_grid_vars["Jdiff_x"].cbegin(), p.device_ptcl_to_grd.cbegin(), qg.num_rows (), qg.hx (), qg.hy (),
76 thrust::raw_pointer_cast(p.device_dprops["VX"].data()), false);
77
78 g2p_step3 VY(thrust::raw_pointer_cast(p.device_x.data()), thrust::raw_pointer_cast(p.device_y.data()), p.device_grid_M.cbegin(),
79 p.device_grid_vars["Jdrift_y"].cbegin(), p.device_grid_vars["Jdiff_y"].cbegin(), p.device_ptcl_to_grd.cbegin(), qg.num_rows (), qg.hx (), qg.hy (),
80 thrust::raw_pointer_cast(p.device_dprops["VY"].data()), false);
81
82 boundary bc(thrust::raw_pointer_cast(p.device_grid_vars["Jdiff_y"].data()), qg.num_cols(), qg.num_rows());
83
84 stepper step(p.device_x.begin(), p.device_y.begin(), p.device_dprops["VX"].begin(), p.device_dprops["VY"].begin(), 1e-8);
85
86 thrust::counting_iterator<idx_t> first_p(0), last_p(p.num_particles), first_n(0), last_n(qg.num_cols());
87
88 constexpr int nsave = 10;
89 constexpr double tmax = 0.02;
90 constexpr double dtsave = tmax / nsave;
91 double t = 0.;
92
93 // Saving iteration stepping
94 for(int isave=0; isave < nsave; ++isave){
95
96 //INTERMIDIATE COPY
97 p.p2g(p.device_grid_vars, {"M"}, {"rho"}, true);
98 p.memcpy_device_to_host();
99 //The following copy cannot be included in the memcpy function since vars is not defined in particles.h
100 thrust::copy (p.device_grid_vars["rho"].cbegin(), p.device_grid_vars["rho"].cend(), vars.at("rho").begin());
101 thrust::copy (p.device_grid_vars["Jdrift_x"].cbegin(), p.device_grid_vars["Jdrift_x"].cend(), vars.at("Jdrift_x").begin());
102 thrust::copy (p.device_grid_vars["Jdrift_y"].cbegin(), p.device_grid_vars["Jdrift_y"].cend(), vars.at("Jdrift_y").begin());
103 thrust::copy (p.device_grid_vars["Jdiff_x"].cbegin(), p.device_grid_vars["Jdiff_x"].cend(), vars.at("Jdiff_x").begin());
104 thrust::copy (p.device_grid_vars["Jdiff_y"].cbegin(), p.device_grid_vars["Jdiff_y"].cend(), vars.at("Jdiff_y").begin());
105
106 // write particle data to file
107 const std::string ofilename = "particle";
108 const std::string ofileext = ".csv";
109 const std::string numfile = std::string(".") + std::to_string(isave);
110
111 std::ofstream outbuf (ofilename + numfile + ofileext);
112 p.print<particles_t::output_format::csv> (outbuf);
113
114 outbuf.close ();
115
116 // write grid data to file
117 const std::string gfilename = std::string("grid.") + std::to_string(isave) + std::string(".vts");
118 qg.vtk_export (gfilename.c_str(), vars);
119
120 //Time Stepping
121
122 while(t < dtsave*(isave +1)){
123 thrust::fill(p.device_grid_vars["Jdrift_x"].begin(), p.device_grid_vars["Jdrift_x"].end(), 0.0);
124 thrust::fill(p.device_grid_vars["Jdrift_y"].begin(), p.device_grid_vars["Jdrift_y"].end(), 0.0);
125 thrust::fill(p.device_grid_vars["Jdiff_x"].begin(), p.device_grid_vars["Jdiff_x"].end(), 0.0);
126 thrust::fill(p.device_grid_vars["Jdiff_y"].begin(), p.device_grid_vars["Jdiff_y"].end(), 0.0);
127 thrust::fill(p.device_grid_vars["rho"].begin(), p.device_grid_vars["rho"].end(), 0.0);
128 thrust::fill(p.device_dprops["BETAx"].begin(), p.device_dprops["BETAx"].end(), 0.0);
129 thrust::fill(p.device_dprops["BETAy"].begin(), p.device_dprops["BETAy"].end(), 0.0);
130 thrust::fill(p.device_dprops["VX"].begin(), p.device_dprops["VX"].end(), 0.0);
131 thrust::fill(p.device_dprops["VY"].begin(), p.device_dprops["VY"].end(), 0.0);
132
133 //G2P of drift velocity
134 timer.tic("g2p_1");
135 p.g2p(p.device_grid_vars, {"betax", "betay"}, {"BETAx", "BETAy"}, false);
136 timer.toc("g2p_1");
137
138 //P2G of Jdrift
139 timer.tic("p2g");
140 p.p2g(p.device_grid_vars, {"M"}, {"rho"}, false);
141 p.p2g(p.device_grid_vars, Jdrift_x);
142 p.p2g(p.device_grid_vars, Jdrift_y);
143
144 thrust::transform(p.device_grid_vars["Jdrift_x"].begin(), p.device_grid_vars["Jdrift_x"].end(),
145 p.device_grid_vars["rho"].begin(), p.device_grid_vars["Jdrift_x"].begin(), safe_divide());
146
147 thrust::transform(p.device_grid_vars["Jdrift_y"].begin(), p.device_grid_vars["Jdrift_y"].end(),
148 p.device_grid_vars["rho"].begin(), p.device_grid_vars["Jdrift_y"].begin(), safe_divide());
149 timer.toc("p2g");
150
151 //P2GD of Jdiff
152 timer.tic("p2gd");
153 p.p2gd(p.device_grid_vars, Jdiff_x);
154 p.p2gd(p.device_grid_vars, Jdiff_y);
155
156 thrust::transform(p.device_grid_vars["Jdiff_x"].begin(), p.device_grid_vars["Jdiff_x"].end(),
157 p.device_grid_vars["rho"].begin(), p.device_grid_vars["Jdiff_x"].begin(), safe_divide());
158
159 thrust::transform(p.device_grid_vars["Jdiff_y"].begin(), p.device_grid_vars["Jdiff_y"].end(),
160 p.device_grid_vars["rho"].begin(), p.device_grid_vars["Jdiff_y"].begin(), safe_divide());
161 //set BC
162 thrust::for_each(thrust::device, first_n, last_n, bc);
163 timer.toc("p2gd");
164
165 // G2P of VX and VY
166 timer.tic("g2p_2");
167 p.g2p(p.device_grid_vars, VX);
168 p.g2p(p.device_grid_vars, VY);
169 timer.toc("g2p_2");
170
171 //Compute dt
172 timer.tic("dt");
173 auto abs_vx = thrust::make_transform_iterator(p.device_dprops["VX"].begin(), abs_no_nan());
174 auto abs_vy = thrust::make_transform_iterator(p.device_dprops["VY"].begin(), abs_no_nan());
175 real_t maxvx = *thrust::max_element(abs_vx, abs_vx + p.num_particles);
176 real_t maxvy = *thrust::max_element(abs_vy, abs_vy + p.num_particles);
177 real_t dtx = maxvx > 0.0 ? .5 * qg.hx() / maxvx : dtsave * (isave + 1) - t;
178 real_t dty = maxvy > 0.0 ? .5 * qg.hy() / maxvy : dtsave * (isave + 1) - t;
179 step.dt = fmin (dtx, dty);
180 if (t + step.dt > dtsave * (isave + 1)) step.dt = dtsave * (isave + 1) - t;
181 timer.toc("dt");
182
183 //Moving Particles
184 timer.tic("move partcles");
185 thrust::for_each(thrust::device, first_p, last_p, step);
186 p.update_ptcl_to_grd<particles_t::update_ptcl_to_grd_device>();
187 timer.toc("move partcles");
188
189 t+=step.dt;
190
191 }
192 }
193
194 //FINAL COPY
195 p.p2g(p.device_grid_vars, {"M"}, {"rho"}, true);
196 p.memcpy_device_to_host();
197 //The following copies cannot be included in the memcpy function since vars is not defined in particles.h
198 thrust::copy (p.device_grid_vars["rho"].cbegin(), p.device_grid_vars["rho"].cend(), vars.at("rho").begin());
199 thrust::copy (p.device_grid_vars["Jdrift_x"].cbegin(), p.device_grid_vars["Jdrift_x"].cend(), vars.at("Jdrift_x").begin());
200 thrust::copy (p.device_grid_vars["Jdrift_y"].cbegin(), p.device_grid_vars["Jdrift_y"].cend(), vars.at("Jdrift_y").begin());
201 thrust::copy (p.device_grid_vars["Jdiff_x"].cbegin(), p.device_grid_vars["Jdiff_x"].cend(), vars.at("Jdiff_x").begin());
202 thrust::copy (p.device_grid_vars["Jdiff_y"].cbegin(), p.device_grid_vars["Jdiff_y"].cend(), vars.at("Jdiff_y").begin());
203
204 // write particle data to file
205 const std::string ofilename = "particle";
206 const std::string ofileext = ".csv";
207 const std::string numfile = std::string(".") + std::to_string(nsave);
208
209 std::ofstream outbuf (ofilename + numfile + ofileext);
210 p.print<particles_t::output_format::csv> (outbuf);
211
212 outbuf.close ();
213
214 // write grid data to file
215 const std::string gfilename = std::string("grid.") + std::to_string(nsave) + std::string(".vts");
216 qg.vtk_export (gfilename.c_str(), vars);
217
218
219
220
221 timer.print_report();
222 return 0;
223}
g2p_step3
Definition drift_diffusion.h:156
p2g_step1
Definition drift_diffusion.h:59
p2gd_step2
Definition drift_diffusion.h:106
idx_t
stepper
Functor class for moving particles.
Definition grid_velocity_field.cpp:22
real_t
double real_t
Definition quadgrid_config.h:48
abs_no_nan
Definition taylor_dispersion.h:25
particles_t
Class to represent particles embedded in a grid.
Definition particles.h:29
particles_t::idx_t
quadgrid_t< vector_t< real_t > >::idx_t idx_t
datatype for indexing into vectors of properties
Definition particles.h:32
particles_t::output_format::csv
@ csv
Definition particles.h:65
safe_divide
Definition taylor_dispersion.h:18
gpuGetDevice
#define gpuGetDevice
Definition taylor_dispersion.cpp:14
gpuSetDevice
#define gpuSetDevice
Definition taylor_dispersion.cpp:13
gpuGetDeviceCount
#define gpuGetDeviceCount
Definition taylor_dispersion.cpp:12