Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
clvisc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file clvisc.cc
1 /*******************************************************************************
2  * Copyright (c) 2018-2019 LongGang Pang, lgpang@qq.com
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a
5  * copy of this software and/or associated documentation files (the
6  * "Materials"), to deal in the Materials without restriction, including
7  * without limitation the rights to use, copy, modify, merge, publish,
8  * distribute, sublicense, and/or sell copies of the Materials, and to
9  * permit persons to whom the Materials are furnished to do so, subject to
10  * the following conditions:
11  *
12  * The above copyright notice and this permission notice shall be included
13  * in all copies or substantial portions of the Materials.
14  *
15  * THE MATERIALS ARE PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18  * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19  * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20  * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21  * MATERIALS OR THE USE OR OTHER DEALINGS IN THE MATERIALS.
22  ******************************************************************************/
23 
24 #include <cmath>
25 #include <cstring>
26 #include <ctime>
27 #include<cstdlib>
28 #include <algorithm>
29 #include "include/clvisc.h"
30 #include "include/error_msgs.h"
31 
32 namespace clvisc {
33 CLVisc::CLVisc(const Config & cfg, std::string device_type,
34  int device_id):cfg_(cfg),
35  ideal_(CLIdeal(cfg, device_type, device_id)),
36  backend_(ideal_.get_backend()),
37  bulkinfo_(cfg.nx, cfg.ny, cfg.nz, cfg.nxskip,
38  cfg.nyskip, cfg.nzskip, backend_,
39  ideal_.get_compile_option())
40 {
41  // set current time to initial time
42  tau0_ = cfg.tau0;
43  tau_ = tau0_;
44 
45  size_ = cfg.nx * cfg.ny * cfg.nz;
46 
47  compile_option_ = ideal_.get_compile_option();
48 
49  std::cout << compile_option_ << std::endl;
50 
51  try {
52  // build kernels for hydrodynamic evolution
53  auto prg = backend_.BuildProgram("clvisc_kernel/kernel_IS.cl",
54  compile_option_);
55  // is == Israel Stewart here
56  kernel_is_initialize_ = cl::Kernel(prg, "visc_initialize");
57  kernel_is_src_christoffel_ = cl::Kernel(prg, "visc_src_christoffel");
58  kernel_is_src_alongx_ = cl::Kernel(prg, "visc_src_alongx");
59  kernel_is_src_alongy_ = cl::Kernel(prg, "visc_src_alongy");
60  kernel_is_src_alongz_ = cl::Kernel(prg, "visc_src_alongz");
61  kernel_is_update_pimn_ = cl::Kernel(prg, "update_pimn");
62  kernel_is_get_udiff_ = cl::Kernel(prg, "get_udiff");
63 
64  auto prg2 = backend_.BuildProgram("clvisc_kernel/kernel_visc.cl",
65  compile_option_);
66  kernel_visc_src_christoffel_ = cl::Kernel(prg2, "kt_src_christoffel");
67  kernel_visc_kt_src_alongx_ = cl::Kernel(prg2, "kt_src_alongx");
68  kernel_visc_kt_src_alongy_ = cl::Kernel(prg2, "kt_src_alongy");
69  kernel_visc_kt_src_alongz_ = cl::Kernel(prg2, "kt_src_alongz");
70  kernel_visc_update_ev_ = cl::Kernel(prg2, "update_ev");
71  } catch (cl::Error & err ){
72  std::cerr<<"Error:"<<err.what()<<"("<<err.err()<<")\n";
73  std::cerr<<"@" << __FILE__ << ":line " << __LINE__ << std::endl;
74  std::cerr<<ErrorMessage(err.err())<<std::endl;
75  throw(err);
76  }
77 }
78 
79 
80 void CLVisc::initialize_gpu_buffer_() {
81  bool read_only_option = false;
82  // ed_bytes is the size of energy density in bytes
83  // d_ev_[i] has size (4*ed_bytes)
84  // d_shear_pi_[i] has size (10*ed_bytes)
85  size_t ed_bytes = sizeof(cl_real) * size_;
86 
87  for (int i = 0; i < 3; i++) {
88  d_shear_pi_[i] = backend_.CreateBufferByCopyVector(h_shear_pi_, read_only_option);
89  d_bulk_pi_[i] = backend_.CreateBufferByCopyVector(h_bulk_pi_, read_only_option);
90  d_net_charge_[i] = backend_.CreateBufferByCopyVector(h_net_charge_, read_only_option);
91  }
92 
93  // initialize the d_udz_ to vector of (real4)
94  // in case the nz = 1 and one wants to skip z direction calculation
95  std::vector<cl_real4> h_udz_(size_, (cl_real4){{0.0f, 0.0f, 0.0f, 0.0f}});
96 
97  d_is_shear_pi_src_ = backend_.CreateBuffer(10*ed_bytes);
98  d_is_bulk_pi_src_ = backend_.CreateBuffer(ed_bytes);
99  d_udx_ = backend_.CreateBuffer(4*ed_bytes); // cl_real4 type
100  d_udy_ = backend_.CreateBuffer(4*ed_bytes);
101  d_udz_ = backend_.CreateBufferByCopyVector(h_udz_, read_only_option);
102  d_udiff_ = backend_.CreateBuffer(4*ed_bytes);
103  d_goodcell_ = backend_.CreateBuffer(ed_bytes);
104 }
105 
106 
107 void CLVisc::israel_stewart_initialize_() {
108  kernel_is_initialize_.setArg(0, d_shear_pi_[1]);
109  kernel_is_initialize_.setArg(1, d_goodcell_);
110  kernel_is_initialize_.setArg(2, d_udiff_);
111  kernel_is_initialize_.setArg(3, ideal_.d_ev_[1]);
112  kernel_is_initialize_.setArg(4, static_cast<cl_real> (tau_));
113  kernel_is_initialize_.setArg(5, ideal_.eos_table_);
114  backend_.enqueue_run(kernel_is_initialize_, cl::NDRange(size_), cl::NullRange);
115 }
116 
117 
118 void CLVisc::update_udiff_() {
119  kernel_is_get_udiff_.setArg(0, d_udiff_);
120  kernel_is_get_udiff_.setArg(1, ideal_.d_ev_[0]);
121  kernel_is_get_udiff_.setArg(2, ideal_.d_ev_[1]);
122  backend_.enqueue_run(kernel_is_get_udiff_,
123  cl::NDRange(size_), cl::NullRange);
124 }
125 
126 void CLVisc::half_step_israel_stewart_(int step) {
127  kernel_is_src_christoffel_.setArg(0, d_is_shear_pi_src_);
128  kernel_is_src_christoffel_.setArg(1, d_shear_pi_[step]);
129  kernel_is_src_christoffel_.setArg(2, ideal_.d_ev_[step]);
130  kernel_is_src_christoffel_.setArg(3, static_cast<cl_real> (tau_));
131  kernel_is_src_christoffel_.setArg(4, step);
132  backend_.enqueue_run(kernel_is_src_christoffel_,
133  cl::NDRange(size_), cl::NullRange);
134 
135  kernel_is_src_alongx_.setArg(0, d_is_shear_pi_src_);
136  kernel_is_src_alongx_.setArg(1, d_udx_);
137  kernel_is_src_alongx_.setArg(2, d_shear_pi_[step]);
138  kernel_is_src_alongx_.setArg(3, ideal_.d_ev_[step]);
139  kernel_is_src_alongx_.setArg(4, ideal_.eos_table_);
140  kernel_is_src_alongx_.setArg(5, static_cast<cl_real> (tau_));
141  backend_.enqueue_run(kernel_is_src_alongx_,
142  cl::NDRange(cfg_.block_size, cfg_.ny, cfg_.nz),
143  cl::NDRange(cfg_.block_size, 1, 1));
144 
145  kernel_is_src_alongy_.setArg(0, d_is_shear_pi_src_);
146  kernel_is_src_alongy_.setArg(1, d_udy_);
147  kernel_is_src_alongy_.setArg(2, d_shear_pi_[step]);
148  kernel_is_src_alongy_.setArg(3, ideal_.d_ev_[step]);
149  kernel_is_src_alongy_.setArg(4, ideal_.eos_table_);
150  kernel_is_src_alongy_.setArg(5, static_cast<cl_real> (tau_));
151  backend_.enqueue_run(kernel_is_src_alongy_,
152  cl::NDRange(cfg_.nx, cfg_.block_size, cfg_.nz),
153  cl::NDRange(1, cfg_.block_size, 1));
154 
155  if (cfg_.nz != 1) {
156  kernel_is_src_alongz_.setArg(0, d_is_shear_pi_src_);
157  kernel_is_src_alongz_.setArg(1, d_udz_);
158  kernel_is_src_alongz_.setArg(2, d_shear_pi_[step]);
159  kernel_is_src_alongz_.setArg(3, ideal_.d_ev_[step]);
160  kernel_is_src_alongz_.setArg(4, ideal_.eos_table_);
161  kernel_is_src_alongz_.setArg(5, static_cast<cl_real> (tau_));
162  backend_.enqueue_run(kernel_is_src_alongz_,
163  cl::NDRange(cfg_.nx, cfg_.ny, cfg_.block_size),
164  cl::NDRange(1, 1, cfg_.block_size));
165  }
166 
167  kernel_is_update_pimn_.setArg(0, d_shear_pi_[3-step]);
168  kernel_is_update_pimn_.setArg(1, d_goodcell_);
169  kernel_is_update_pimn_.setArg(2, d_shear_pi_[1]);
170  kernel_is_update_pimn_.setArg(3, d_shear_pi_[step]);
171  kernel_is_update_pimn_.setArg(4, ideal_.d_ev_[1]);
172  kernel_is_update_pimn_.setArg(5, ideal_.d_ev_[2]);
173  kernel_is_update_pimn_.setArg(6, d_udiff_);
174  kernel_is_update_pimn_.setArg(7, d_udx_);
175  kernel_is_update_pimn_.setArg(8, d_udy_);
176  kernel_is_update_pimn_.setArg(9, d_udz_);
177  kernel_is_update_pimn_.setArg(10, d_is_shear_pi_src_);
178  kernel_is_update_pimn_.setArg(11, ideal_.eos_table_);
179  kernel_is_update_pimn_.setArg(12, static_cast<cl_real> (tau_));
180  kernel_is_update_pimn_.setArg(13, step);
181  backend_.enqueue_run(kernel_is_update_pimn_,
182  cl::NDRange(size_), cl::NullRange);
183 
184 }
185 
186 void CLVisc::half_step_visc_(int step) {
187  kernel_visc_src_christoffel_.setArg(0, ideal_.d_src_);
188  kernel_visc_src_christoffel_.setArg(1, ideal_.d_ev_[step]);
189  kernel_visc_src_christoffel_.setArg(2, d_shear_pi_[step]);
190  kernel_visc_src_christoffel_.setArg(3, ideal_.eos_table_);
191  kernel_visc_src_christoffel_.setArg(4, static_cast<cl_real> (tau_));
192  kernel_visc_src_christoffel_.setArg(5, step);
193  backend_.enqueue_run(kernel_visc_src_christoffel_,
194  cl::NDRange(size_), cl::NullRange);
195 
196  kernel_visc_kt_src_alongx_.setArg(0, ideal_.d_src_);
197  kernel_visc_kt_src_alongx_.setArg(1, ideal_.d_ev_[step]);
198  kernel_visc_kt_src_alongx_.setArg(2, d_shear_pi_[step]);
199  kernel_visc_kt_src_alongx_.setArg(3, ideal_.eos_table_);
200  kernel_visc_kt_src_alongx_.setArg(4, static_cast<cl_real> (tau_));
201  backend_.enqueue_run(kernel_visc_kt_src_alongx_,
202  cl::NDRange(cfg_.block_size, cfg_.ny, cfg_.nz),
203  cl::NDRange(cfg_.block_size, 1, 1));
204 
205  kernel_visc_kt_src_alongy_.setArg(0, ideal_.d_src_);
206  kernel_visc_kt_src_alongy_.setArg(1, ideal_.d_ev_[step]);
207  kernel_visc_kt_src_alongy_.setArg(2, d_shear_pi_[step]);
208  kernel_visc_kt_src_alongy_.setArg(3, ideal_.eos_table_);
209  kernel_visc_kt_src_alongy_.setArg(4, static_cast<cl_real> (tau_));
210  backend_.enqueue_run(kernel_visc_kt_src_alongy_,
211  cl::NDRange(cfg_.nx, cfg_.block_size, cfg_.nz),
212  cl::NDRange(1, cfg_.block_size, 1));
213 
214 
215  if (cfg_.nz != 1) {
216  kernel_visc_kt_src_alongz_.setArg(0, ideal_.d_src_);
217  kernel_visc_kt_src_alongz_.setArg(1, ideal_.d_ev_[step]);
218  kernel_visc_kt_src_alongz_.setArg(2, d_shear_pi_[step]);
219  kernel_visc_kt_src_alongz_.setArg(3, ideal_.eos_table_);
220  kernel_visc_kt_src_alongz_.setArg(4, static_cast<cl_real> (tau_));
221  backend_.enqueue_run(kernel_visc_kt_src_alongz_,
222  cl::NDRange(cfg_.nx, cfg_.ny, cfg_.block_size),
223  cl::NDRange(1, 1, cfg_.block_size));
224  }
225 
226  kernel_visc_update_ev_.setArg(0, ideal_.d_ev_[3-step]);
227  kernel_visc_update_ev_.setArg(1, ideal_.d_ev_[1]);
228  kernel_visc_update_ev_.setArg(2, d_shear_pi_[0]);
229  kernel_visc_update_ev_.setArg(3, d_shear_pi_[3-step]);
230  kernel_visc_update_ev_.setArg(4, ideal_.d_src_);
231  kernel_visc_update_ev_.setArg(5, ideal_.eos_table_);
232  kernel_visc_update_ev_.setArg(6, static_cast<cl_real> (tau_));
233  kernel_visc_update_ev_.setArg(7, step);
234  backend_.enqueue_run(kernel_visc_update_ev_,
235  cl::NDRange(size_), cl::NullRange);
236 } // end half_step_visc_
237 
238 
239 void CLVisc::half_step_(int step) {
240  half_step_israel_stewart_(step);
241  half_step_visc_(step);
242 }
243 
244 void CLVisc::one_step() {
245  half_step_(1);
246  tau_ += cfg_.dt;
247  half_step_(2);
248 }
249 
250 
251 template <typename ValueType>
252 void CLVisc::read_ini(const std::vector<ValueType> & ed) {
253  ideal_.read_ini(ed);
254  h_shear_pi_ = std::vector<cl_real> (10*size_, 0.0);
255  h_bulk_pi_ = std::vector<cl_real> (size_, 0.0);
256  cl_real4 zero4 = (cl_real4){{0.0, 0.0, 0.0, 0.0}};
257  h_net_charge_ = std::vector<cl_real4> (size_, zero4);
258  initialize_gpu_buffer_();
259 }
260 
261 template <typename ValueType>
262 void CLVisc::read_ini(const std::vector<ValueType> & ed,
263  const std::vector<ValueType> & vx,
264  const std::vector<ValueType> & vy,
265  const std::vector<ValueType> & vz)
266 {
267  ideal_.read_ini(ed, vx, vy, vz);
268  h_shear_pi_ = std::vector<cl_real> (10*size_, 0.0);
269  h_bulk_pi_ = std::vector<cl_real> (size_, 0.0);
270  cl_real4 zero4 = (cl_real4){{0.0, 0.0, 0.0, 0.0}};
271  h_net_charge_ = std::vector<cl_real4> (size_, zero4);
272  initialize_gpu_buffer_();
273 }
274 
275 template <typename ValueType>
276 void CLVisc::read_ini(const std::vector<ValueType> & ed,
277  const std::vector<ValueType> & vx,
278  const std::vector<ValueType> & vy,
279  const std::vector<ValueType> & vz,
280  const std::vector<ValueType> & pi00,
281  const std::vector<ValueType> & pi01,
282  const std::vector<ValueType> & pi02,
283  const std::vector<ValueType> & pi03,
284  const std::vector<ValueType> & pi11,
285  const std::vector<ValueType> & pi12,
286  const std::vector<ValueType> & pi13,
287  const std::vector<ValueType> & pi22,
288  const std::vector<ValueType> & pi23,
289  const std::vector<ValueType> & pi33) {
290  ideal_.read_ini(ed, vx, vy, vz);
291  h_shear_pi_.clear();
292  for (size_t i = 0; i < size_; i++) {
293  h_shear_pi_.push_back(static_cast<cl_real>(pi00.at(i)));
294  h_shear_pi_.push_back(static_cast<cl_real>(pi01.at(i)));
295  h_shear_pi_.push_back(static_cast<cl_real>(pi02.at(i)));
296  h_shear_pi_.push_back(static_cast<cl_real>(pi03.at(i)));
297  h_shear_pi_.push_back(static_cast<cl_real>(pi11.at(i)));
298  h_shear_pi_.push_back(static_cast<cl_real>(pi12.at(i)));
299  h_shear_pi_.push_back(static_cast<cl_real>(pi13.at(i)));
300  h_shear_pi_.push_back(static_cast<cl_real>(pi22.at(i)));
301  h_shear_pi_.push_back(static_cast<cl_real>(pi23.at(i)));
302  h_shear_pi_.push_back(static_cast<cl_real>(pi33.at(i)));
303  }
304  h_bulk_pi_ = std::vector<cl_real> (size_, 0.0);
305  cl_real4 zero4 = (cl_real4){{0.0, 0.0, 0.0, 0.0}};
306  h_net_charge_ = std::vector<cl_real4> (size_, zero4);
307  initialize_gpu_buffer_();
308 }
309 
310 
311 
312 
313 void CLVisc::evolve() {
314  int max_loops = 5000;
315  float total_exec_time = 0.0;
316  std::time_t timer1, timer2;
317  std::time(&timer1);
318  tau_ = tau0_;
319  try {
320  israel_stewart_initialize_();
321  ideal_.predict_first_step();
322  for (int loop = 0; loop < max_loops; loop++) {
323  backend_.enqueue_copy(d_shear_pi_[1], d_shear_pi_[0], 10*size_*sizeof(cl_real));
324  backend_.enqueue_copy(ideal_.d_ev_[1], ideal_.d_ev_[0], size_*sizeof(cl_real4));
325  one_step();
326  update_udiff_();
327  if (loop % cfg_.ntskip == 0) {
328  float max_ed = ideal_.max_energy_density();
329  std::cout << "tau = " << tau_ << " fm; ";
330  std::cout << "max_ed = " << max_ed << " ";
331  std::time(&timer2);
332  total_exec_time = std::difftime(timer2, timer1);
333  //std::cout << "Total computing time: " << total_exec_time << " s; ";
334  //std::cout << std::endl;
335  bulkinfo_.add_data(ideal_.d_ev_[0], d_shear_pi_[0],
336  d_bulk_pi_[0], ideal_.eos_table_);
337  if ( max_ed < 0.05 ) break;
338  }
339  }
340  std::cout << std::endl;
341 
342  std::cout << "Total computing time: " << total_exec_time << " s; ";
343  } catch (cl::Error & err) {
344  std::cout << err.what() << " " << err.err() << std::endl;
345  std::cerr<<"@" << __FILE__ << ":line " << __LINE__ << std::endl;
346  std::cout << ErrorMessage(err.err()) << std::endl;
347  throw(err);
348  }
349 }
350 
351 CLVisc::~CLVisc() {
352 }
353 
354 // in case the input array is a vector of double instead of float
355 template void CLVisc::read_ini(const std::vector<float> & ed);
356 template void CLVisc::read_ini(const std::vector<double> & ed);
357 
358 template void CLVisc::read_ini(const std::vector<float> & ed,
359  const std::vector<float> & vx,
360  const std::vector<float> & vy,
361  const std::vector<float> & vz);
362 
363 template void CLVisc::read_ini(const std::vector<double> & ed,
364  const std::vector<double> & vx,
365  const std::vector<double> & vy,
366  const std::vector<double> & vz);
367 
368 
369 template void CLVisc::read_ini(const std::vector<float> & ed,
370  const std::vector<float> & vx,
371  const std::vector<float> & vy,
372  const std::vector<float> & vz,
373  const std::vector<float> & pi00,
374  const std::vector<float> & pi01,
375  const std::vector<float> & pi02,
376  const std::vector<float> & pi03,
377  const std::vector<float> & pi11,
378  const std::vector<float> & pi12,
379  const std::vector<float> & pi13,
380  const std::vector<float> & pi22,
381  const std::vector<float> & pi23,
382  const std::vector<float> & pi33);
383 
384 template void CLVisc::read_ini(const std::vector<double> & ed,
385  const std::vector<double> & vx,
386  const std::vector<double> & vy,
387  const std::vector<double> & vz,
388  const std::vector<double> & pi00,
389  const std::vector<double> & pi01,
390  const std::vector<double> & pi02,
391  const std::vector<double> & pi03,
392  const std::vector<double> & pi11,
393  const std::vector<double> & pi12,
394  const std::vector<double> & pi13,
395  const std::vector<double> & pi22,
396  const std::vector<double> & pi23,
397  const std::vector<double> & pi33);
398 
399 
400 } // end namespace clvisc