OPALX (Object Oriented Parallel Accelerator Library for Exascal) MINIorX
OPALX
FieldSolver.cpp
Go to the documentation of this file.
1#include "PartBunch/FieldSolver.hpp"
2
3#include <iomanip>
4#include <fstream>
5
6#include <boost/filesystem.hpp>
7#include <boost/format.hpp>
8
9#include "Utilities/Util.h"
10#include "AbstractObjects/OpalData.h"
11
12extern Inform* gmsg;
13
14template <>
15void FieldSolver<double,3>::dumpVectField(std::string what) {
16 /*
17 what == ef
18 */
19
20 Inform m("FS::dumpVectorField() ");
21
22 // std::variant<Field_t<3>*, VField_t<double, 3>* > field;
23
24 if (ippl::Comm->size() > 1 || call_counter_m<2) {
25 return;
26 }
27
28/* Save the files in the output directory of the simulation. The file
29 * name of vector fields is
30 *
31 * 'basename'-'name'_field-'******'.dat
32 *
33 * and of scalar fields
34 *
35 * 'basename'-'name'_scalar-'******'.dat
36 *
37 * with
38 * 'basename': OPAL input file name (*.in)
39 * 'name': field name (input argument of function)
40 * '******': call_counter_m padded with zeros to 6 digits
41 */
42
43 std::string dirname = "data/";
44
45 std::string type;
46 std::string unit;
47
48 if (Util::toUpper(what) == "EF") {
49 type = "vector";
50 unit = "";
51 }
52
53 VField_t<double, 3>* field = this->getE();
54
55 auto localIdx = field->getOwned();
56 auto mesh_mp = &(field->get_mesh());
57 auto spacing = mesh_mp->getMeshSpacing();
58 auto origin = mesh_mp->getOrigin();
59
60 auto fieldV = field->getView();
61 auto field_hostV = field->getHostMirror();
62 Kokkos::deep_copy(field_hostV, fieldV);
63
64 boost::filesystem::path file(dirname);
65 boost::format filename("%1%-%2%-%|3$06|.dat");
66 std::string basename = OpalData::getInstance()->getInputBasename();
67 filename % basename % (what + std::string("_") + type) % call_counter_m;
68 file /= filename.str();
69 std::ofstream fout(file.string(), std::ios::out);
70
71 fout << std::setprecision(9);
72
73 fout << "# " << Util::toUpper(what) << " " << type << " data on grid" << std::endl
74 << "# origin= " << std::fixed << origin << " h= " << std::fixed << spacing << std::endl
75 << std::setw(5) << "i"
76 << std::setw(5) << "j"
77 << std::setw(5) << "k"
78 << std::setw(17) << "x [m]"
79 << std::setw(17) << "y [m]"
80 << std::setw(17) << "z [m]";
81
82 fout << std::setw(10) << what << "x [" << unit << "]"
83 << std::setw(10) << what << "y [" << unit << "]"
84 << std::setw(10) << what << "z [" << unit << "]";
85
86 fout << std::endl;
87
88 for (int i = localIdx[0].first() + 1; i <= localIdx[0].last() +1 ; i++) {
89 for (int j = localIdx[1].first() + 1 ; j <= localIdx[1].last() +1 ; j++) {
90 for (int k = localIdx[2].first() + 1 ; k <= localIdx[2].last() +1 ; k++) {
91
92 // define the physical points (cell-centered)
93 double x = i * spacing[0] + origin[0];
94 double y = j * spacing[1] + origin[1];
95 double z = k * spacing[2] + origin[2];
96
97 fout << std::setw(5) << i-1
98 << std::setw(5) << j-1
99 << std::setw(5) << k-1
100 << std::setw(17) << x
101 << std::setw(17) << y
102 << std::setw(17) << z
103 << std::scientific
104 << "\t" << field_hostV(i,j,k)[0]
105 << "\t" << field_hostV(i,j,k)[1]
106 << "\t" << field_hostV(i,j,k)[2]
107 << std::endl;
108 }
109 }
110 }
111 fout.close();
112 m << "*** FINISHED DUMPING " + Util::toUpper(what) + " FIELD *** to " << file.string() << endl;
113}
114
115template <>
116void FieldSolver<double,3>::dumpScalField(std::string what) {
117
118 /*
119 what == phi | rho
120 */
121
122 Inform m("FS::dumpScalField() ");
123
124 if (ippl::Comm->size() > 1 || call_counter_m<2) {
125 return;
126 }
127
128/* Save the files in the output directory of the simulation. The file
129 * name of vector fields is
130 *
131 * 'basename'-'name'_field-'******'.dat
132 *
133 * and of scalar fields
134 *
135 * 'basename'-'name'_scalar-'******'.dat
136 *
137 * with
138 * 'basename': OPAL input file name (*.in)
139 * 'name': field name (input argument of function)
140 * '******': call_counter_m padded with zeros to 6 digits
141 */
142
143
144 std::string dirname = "data/";
145
146 std::string type;
147 std::string unit;
148 bool isVectorField = false;
149
150 if (Util::toUpper(what) == "RHO") {
151 type = "scalar";
152 unit = "Cb/m^3";
153 } else if (Util::toUpper(what) == "PHI") {
154 type = "scalar";
155 unit = "V";
156 }
157
158
159 Field_t<3>* field = this->getRho(); // both rho and phi are in the same variable (in place computation)
160
161 auto localIdx = field->getOwned();
162 auto mesh_mp = &(field->get_mesh());
163 auto spacing = mesh_mp->getMeshSpacing();
164 auto origin = mesh_mp->getOrigin();
165
166 auto fieldV = field->getView();
167 auto field_hostV = field->getHostMirror();
168 Kokkos::deep_copy(field_hostV, fieldV);
169
170 boost::filesystem::path file(dirname);
171 boost::format filename("%1%-%2%-%|3$06|.dat");
172 std::string basename = OpalData::getInstance()->getInputBasename();
173 filename % basename % (what + std::string("_") + type) % call_counter_m;
174 file /= filename.str();
175 std::ofstream fout(file.string(), std::ios::out);
176
177 fout << std::setprecision(9);
178
179 fout << "# " << Util::toUpper(what) << " " << type << " data on grid" << std::endl
180 << "# origin= " << std::fixed << origin << " h= " << std::fixed << spacing << std::endl
181 << std::setw(5) << "i"
182 << std::setw(5) << "j"
183 << std::setw(5) << "k"
184 << std::setw(17) << "x [m]"
185 << std::setw(17) << "y [m]"
186 << std::setw(17) << "z [m]";
187
188 if (isVectorField) {
189 fout << std::setw(10) << what << "x [" << unit << "]"
190 << std::setw(10) << what << "y [" << unit << "]"
191 << std::setw(10) << what << "z [" << unit << "]";
192 } else {
193 fout << std::setw(13) << what << " [" << unit << "]";
194 }
195
196 fout << std::endl;
197
198 if (Util::toUpper(what) == "RHO") {
199 for (int i = localIdx[0].first(); i <= localIdx[0].last(); i++) {
200 for (int j = localIdx[1].first(); j <= localIdx[1].last(); j++) {
201 for (int k = localIdx[2].first(); k <= localIdx[2].last(); k++) {
202
203 // define the physical points (cell-centered)
204 double x = i * spacing[0] + origin[0];
205 double y = j * spacing[1] + origin[1];
206 double z = k * spacing[2] + origin[2];
207
208 fout << std::setw(5) << i
209 << std::setw(5) << j
210 << std::setw(5) << k
211 << std::setw(17) << x
212 << std::setw(17) << y
213 << std::setw(17) << z
214 << std::scientific << "\t" << field_hostV(i,j,k)
215 << std::endl;
216 }
217 }
218 }
219 }
220 else {
221 for (int i = localIdx[0].first() + 1; i <= localIdx[0].last() +1 ; i++) {
222 for (int j = localIdx[1].first() + 1 ; j <= localIdx[1].last() +1 ; j++) {
223 for (int k = localIdx[2].first() + 1 ; k <= localIdx[2].last() +1 ; k++) {
224
225 // define the physical points (cell-centered)
226 double x = i * spacing[0] + origin[0];
227 double y = j * spacing[1] + origin[1];
228 double z = k * spacing[2] + origin[2];
229
230 fout << std::setw(5) << i-1
231 << std::setw(5) << j-1
232 << std::setw(5) << k-1
233 << std::setw(17) << x
234 << std::setw(17) << y
235 << std::setw(17) << z
236 << std::scientific << "\t" << field_hostV(i,j,k)
237 << std::endl;
238 }
239 }
240 }
241 }
242 fout.close();
243 m << "*** FINISHED DUMPING " + Util::toUpper(what) + " FIELD *** to " << file.string() << endl;
244}
245
246template <>
247void FieldSolver<double,3>::initOpenSolver() {
248 ippl::ParameterList sp;
249 sp.add("output_type", OpenSolver_t<double, 3>::SOL_AND_GRAD);
250 sp.add("use_heffte_defaults", false);
251 sp.add("use_pencils", true);
252 sp.add("use_reorder", false);
253 sp.add("use_gpu_aware", true);
254 sp.add("comm", ippl::p2p_pl);
255 sp.add("r2c_direction", 0);
256 sp.add("algorithm", OpenSolver_t<double, 3>::HOCKNEY);
257 initSolverWithParams<OpenSolver_t<double, 3>>(sp);
258}
259
260template <>
261void FieldSolver<double,3>::initSolver() {
262 Inform m;
263 if (this->getStype() == "FFT") {
264 initOpenSolver();
265 } else if (this->getStype() == "FFTOPEN") {
266 initOpenSolver();
267 } else if (this->getStype() == "NONE") {
268 initNullSolver();
269 }
270 else {
271 m << "No solver matches the argument: " << this->getStype() << endl;
272 throw std::runtime_error("No solver match");
273 }
274}
275
276template <>
277void FieldSolver<double,3>::setPotentialBCs() {
278 // CG requires explicit periodic boundary conditions while the periodic Poisson solver
279 // simply assumes them
280 if (this->getStype() == "CG") {
281 typedef ippl::BConds<Field<double, 3>, 3> bc_type;
282 bc_type allPeriodic;
283 for (unsigned int i = 0; i < 2 * 3; ++i) {
284 allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<double, 3>>>(i);
285 }
286 phi_m->setFieldBC(allPeriodic);
287 }
288 }
289
290template<>
291void FieldSolver<double,3>::runSolver() {
292 constexpr int Dim = 3;
293
294 if (this->getStype() == "CG") {
295 CGSolver_t<double, 3>& solver = std::get<CGSolver_t<double, 3>>(this->getSolver());
296 solver.solve();
297
298 if (ippl::Comm->rank() == 0) {
299 std::stringstream fname;
300 fname << "data/CG_";
301 fname << ippl::Comm->size();
302 fname << ".csv";
303
304 Inform log(NULL, fname.str().c_str(), Inform::APPEND);
305 int iterations = solver.getIterationCount();
306 // Assume the dummy solve is the first call
307 if (iterations == 0) {
308 log << "residue,iterations" << endl;
309 }
310 // Don't print the dummy solve
311 if (iterations > 0) {
312 log << solver.getResidue() << "," << iterations << endl;
313 }
314 }
315 ippl::Comm->barrier();
316 } else if (this->getStype() == "FFT") {
317 if constexpr (Dim == 2 || Dim == 3) {
318#ifdef OPALX_FIELD_DEBUG
319 this->dumpScalField("rho");
320 call_counter_m++;
321#endif
322
323 std::get<OpenSolver_t<double, 3>>(this->getSolver()).solve();
324#ifdef OPALX_FIELD_DEBUG
325 this->dumpScalField("phi");
326 this->dumpVectField("ef");
327 call_counter_m++;
328#endif
329 }
330 } else if (this->getStype() == "P3M") {
331 if constexpr (Dim == 3) {
332 std::get<FFTTruncatedGreenSolver_t<double, 3>>(this->getSolver()).solve();
333 }
334 } else if (this->getStype() == "FFTOPEN") {
335 if constexpr (Dim == 3) {
336#ifdef OPALX_FIELD_DEBUG
337 this->dumpScalField("rho");
338 call_counter_m++;
339#endif
340 std::get<OpenSolver_t<double, 3>>(this->getSolver()).solve();
341#ifdef OPALX_FIELD_DEBUG
342 this->dumpScalField("phi");
343 this->dumpVectField("ef");
344 call_counter_m++;
345#endif
346 }
347 } else if (this->getStype() == "NONE") {
348 std::get<NullSolver_t<T, Dim>>(this->getSolver()).solve();
349 } else {
350 throw std::runtime_error("Unknown solver type");
351 }
352}
353
354template<>
355void FieldSolver<double,3>::initNullSolver() {
356 ippl::ParameterList sp;
357 if constexpr (Dim == 2 || Dim == 3) {
358 initSolverWithParams<NullSolver_t<T, Dim>>(sp);
359 } else {
360 throw std::runtime_error("Unsupported dimensionality for Null solver");
361 }
362}
363
364/*
365template <>
366void FieldSolver<double,3>::initCGSolver() {
367 ippl::ParameterList sp;
368 sp.add("output_type", CGSolver_t<double, 3>::GRAD);
369 // Increase tolerance in the 1D case
370 sp.add("tolerance", 1e-10);
371
372 initSolverWithParams<CGSolver_t<double, 3>>(sp);
373}
374
375template<>
376void FieldSolver<double,3>::initP3MSolver() {
377 // if constexpr (Dim == 3) {
378 ippl::ParameterList sp;
379 sp.add("output_type", P3MSolver_t<double, 3>::GRAD);
380 sp.add("use_heffte_defaults", false);
381 sp.add("use_pencils", true);
382 sp.add("use_reorder", false);
383 sp.add("use_gpu_aware", true);
384 sp.add("comm", ippl::p2p_pl);
385 sp.add("r2c_direction", 0);
386
387 initSolverWithParams<P3MSolver_t<double, 3>>(sp);
388 // } else {
389 // throw std::runtime_error("Unsupported dimensionality for P3M solver");
390 // }
391}
392
393*/
gmsg
Inform * gmsg
Definition changes.cpp:7
OpenSolver_t
ConditionalType< Dim==3, ippl::FFTOpenPoissonSolver< VField_t< T, Dim >, Field_t< Dim > > > OpenSolver_t
Definition #PartBunch.hpp#:54
Field_t
ippl::Field< double, Dim, ViewArgs... > Field_t
Definition PBunchDefs.h:30
VField_t
ippl::Field< ippl::Vector< T, Dim >, Dim, ViewArgs... > VField_t
Definition PBunchDefs.h:33
Dim
constexpr unsigned Dim
Definition datatypes.h:6
Util::toUpper
std::string toUpper(const std::string &str)
Definition Util.cpp:145
OpalData::getInputBasename
std::string getInputBasename()
get input file name without extension
Definition OpalData.cpp:685
OpalData::getInstance
static OpalData * getInstance()
Definition OpalData.cpp:195
FieldSolver::getE
VField_t< T, Dim > * getE() const
Definition FieldSolver.hpp:35
FieldSolver::dumpVectField
void dumpVectField(std::string what)
FieldSolver::initSolver
void initSolver() override
FieldSolver::runSolver
void runSolver() override
FieldSolver::initSolverWithParams
void initSolverWithParams(const ippl::ParameterList &sp)
Definition FieldSolver.hpp:58
FieldSolver::getRho
Field_t< Dim > * getRho()
Definition FieldSolver.hpp:28
FieldSolver::initOpenSolver
void initOpenSolver()
FieldSolver::initNullSolver
void initNullSolver()
FieldSolver::call_counter_m
unsigned int call_counter_m
Definition FieldSolver.hpp:15
FieldSolver::dumpScalField
void dumpScalField(std::string what)
FieldSolver::setPotentialBCs
void setPotentialBCs()
FieldSolver::phi_m
Field_t< Dim > * phi_m
Definition FieldSolver.hpp:14