Line data Source code
1 : // Copyright (C) 2015 Technische Universitaet Muenchen
2 : // This file is part of the Mamico project. For conditions of distribution
3 : // and use, please see the copyright notice in Mamico's main folder, or at
4 : // www5.in.tum.de/mamico
5 : #ifndef _MOLECULARDYNAMICS_COUPLING_SOLVERS_NUMERICALSOLVER_H_
6 : #define _MOLECULARDYNAMICS_COUPLING_SOLVERS_NUMERICALSOLVER_H_
7 :
8 : #include "coupling/CouplingMDDefinitions.h"
9 : #include "tarch/la/Vector.h"
10 : #include <cmath>
11 : #include <fstream>
12 : #include <iomanip>
13 : #include <sstream>
14 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
15 : #include <mpi.h>
16 : #endif
17 : #include "coupling/datastructures/CouplingCell.h"
18 : #include "coupling/datastructures/FlexibleCellContainer.h"
19 : #include "coupling/indexing/IndexingService.h"
20 : #include "coupling/services/ParallelTimeIntegrationService.h"
21 : #include "coupling/solvers/CouetteSolver.h"
22 : namespace coupling {
23 : namespace solvers {
24 : class NumericalSolver;
25 : }
26 : } // namespace coupling
27 :
28 : /** The setup is, a 3d solver. A channel flow for the Couette scenario, where
29 : * the moving wall is located at the lower wall (z=0)
30 : * @brief is a virtual base class for the interface for a numerical fluid
31 : * solver for the Couette scenario
32 : * @author Philipp Neumann & Helene Wittenberg */
33 : class coupling::solvers::NumericalSolver : public coupling::solvers::AbstractCouetteSolver<3> {
34 : public:
35 : /** @brief a simple constructor
36 : * @param channelheight the width and height of the channel in y and z
37 : * direction
38 : * @param dx the spacial step size, and equidistant grid is applied
39 : * @param dt the time step
40 : * @param kinVisc the kinematic viscosity of the fluid
41 : * @param plotEveryTimestep the time step interval for plotting data;
42 : * 4 means, every 4th time step is plotted
43 : * @param filestem the name of the plotted file
44 : * @param processes defines on how many processes the solver will run;
45 : * 1,1,1 - sequential run - 1,2,2 = 1*2*2 = 4 processes */
46 116 : NumericalSolver(const double channelheight, const double dx, const double dt, const double kinVisc, const int plotEveryTimestep, const std::string filestem,
47 : const tarch::la::Vector<3, unsigned int> processes, const Scenario* scen = nullptr)
48 232 : : coupling::solvers::AbstractCouetteSolver<3>(), _channelheight(channelheight), _dx(dx), _dt(dt), _kinVisc(kinVisc), _processes(processes),
49 1044 : _plotEveryTimestep(plotEveryTimestep), _filestem(filestem), _scen(scen) {
50 116 : _vel = new double[3 * (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2)];
51 116 : _density = new double[(_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2)];
52 116 : _flag = new Flag[(_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2)];
53 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
54 116 : _sendBufferX = new double[5 * (_domainSizeY + 2) * (_domainSizeZ + 2)];
55 116 : _recvBufferX = new double[5 * (_domainSizeY + 2) * (_domainSizeZ + 2)];
56 116 : _sendBufferY = new double[5 * (_domainSizeX + 2) * (_domainSizeZ + 2)];
57 116 : _recvBufferY = new double[5 * (_domainSizeX + 2) * (_domainSizeZ + 2)];
58 116 : _sendBufferZ = new double[5 * (_domainSizeX + 2) * (_domainSizeY + 2)];
59 116 : _recvBufferZ = new double[5 * (_domainSizeX + 2) * (_domainSizeY + 2)];
60 : #endif
61 : // zero velocity, unit density; flags are set to FLUID
62 309604 : for (int i = 0; i < (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2); i++) {
63 1237952 : for (int d = 0; d < 3; d++) {
64 928464 : _vel[i * 3 + d] = (double)0.0;
65 : }
66 309488 : _density[i] = 1.0;
67 309488 : _flag[i] = FLUID;
68 : }
69 : // determine parallel neighbours
70 116 : determineParallelNeighbours();
71 : // correct boundary flags based on physical description (Couette scenario)
72 : // bottom - moving wall
73 14500 : for (int i = 0; i < (_domainSizeX + 2) * (_domainSizeY + 2); i++) {
74 14384 : _flag[i] = MOVING_WALL;
75 : }
76 : // top - noslip
77 14500 : for (int i = (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 1); i < (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2); i++) {
78 14384 : _flag[i] = NO_SLIP;
79 : }
80 : // front - periodic
81 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
82 13340 : for (int x = 0; x < _domainSizeX + 2; x++) {
83 12760 : _flag[get(x, 0, z)] = PERIODIC;
84 : }
85 : }
86 : // back - periodic
87 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
88 13340 : for (int x = 0; x < _domainSizeX + 2; x++) {
89 12760 : _flag[get(x, _domainSizeY + 1, z)] = PERIODIC;
90 : }
91 : }
92 : // left - periodic
93 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
94 12180 : for (int y = 1; y < _domainSizeY + 1; y++) {
95 11600 : _flag[get(0, y, z)] = PERIODIC;
96 : }
97 : }
98 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
99 12180 : for (int y = 1; y < _domainSizeY + 1; y++) {
100 11600 : _flag[get(_domainSizeX + 1, y, z)] = PERIODIC;
101 : }
102 : }
103 : // correct boundary flags in case of MPI-parallel simulations (Couette
104 : // scenario)
105 116 : setParallelBoundaryFlags();
106 116 : }
107 :
108 : /** @brief a simple destructor */
109 116 : virtual ~NumericalSolver() {
110 116 : if (_vel != NULL) {
111 0 : delete[] _vel;
112 0 : _vel = NULL;
113 : }
114 116 : if (_density != NULL) {
115 0 : delete[] _density;
116 0 : _density = NULL;
117 : }
118 116 : if (_flag != NULL) {
119 0 : delete[] _flag;
120 0 : _flag = NULL;
121 : }
122 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
123 116 : if (_sendBufferX != NULL) {
124 0 : delete[] _sendBufferX;
125 0 : _sendBufferX = NULL;
126 : }
127 116 : if (_sendBufferY != NULL) {
128 0 : delete[] _sendBufferY;
129 0 : _sendBufferY = NULL;
130 : }
131 116 : if (_sendBufferZ != NULL) {
132 0 : delete[] _sendBufferZ;
133 0 : _sendBufferZ = NULL;
134 : }
135 116 : if (_recvBufferX != NULL) {
136 0 : delete[] _recvBufferX;
137 0 : _recvBufferX = NULL;
138 : }
139 116 : if (_recvBufferY != NULL) {
140 0 : delete[] _recvBufferY;
141 0 : _recvBufferY = NULL;
142 : }
143 116 : if (_recvBufferZ != NULL) {
144 0 : delete[] _recvBufferZ;
145 0 : _recvBufferZ = NULL;
146 : }
147 : #endif
148 116 : }
149 :
150 : /** @brief flags the domain boundary cells.
151 : * @param mdDomainOffset lower/left/front corner of the MD domain
152 : * @param mdDomainSize total 3d size of the md domain
153 : * @param overlapStrip the number of cells in the overlap layer;
154 : * The overlap of md and macro cells
155 : * @param recvIndice the coupling cell indices that will be received
156 : * @param size the number of cells that will be received */
157 : void setMDBoundary(tarch::la::Vector<3, double> mdDomainOffset, tarch::la::Vector<3, double> mdDomainSize, unsigned int overlapStrip) {
158 : if (skipRank()) {
159 : return;
160 : }
161 :
162 : for (int d = 0; d < 3; d++) {
163 : _offset[d] = (floor(mdDomainOffset[d] / _dx + 0.5));
164 : if (fabs(_offset[d] * _dx - mdDomainOffset[d]) / _dx > 1.0e-8) {
165 : std::cout << "ERROR NumericalSolver::setMDBoundary(): offset does not match!" << std::endl;
166 : exit(EXIT_FAILURE);
167 : }
168 : _globalNumberCouplingCells[d] = (floor(mdDomainSize[d] / _dx + 0.5));
169 : if (fabs(_globalNumberCouplingCells[d] * _dx - mdDomainSize[d]) / _dx > 1.0e-8) {
170 : std::cout << "ERROR NumericalSolver::setMDBoundary(): globalNumber "
171 : "does not match!"
172 : << std::endl;
173 : exit(EXIT_FAILURE);
174 : }
175 : }
176 : // flag local domain
177 : for (int z = 0; z < _domainSizeZ + 2; z++) {
178 : for (int y = 0; y < _domainSizeY + 2; y++) {
179 : for (int x = 0; x < _domainSizeX + 2; x++) {
180 : // determine global cell coordinates of the process-local (sub)domain
181 : tarch::la::Vector<3, int> globalCoords(x - 1 + _coords[0] * _avgDomainSizeX, y - 1 + _coords[1] * _avgDomainSizeY,
182 : z - 1 + _coords[2] * _avgDomainSizeZ);
183 : bool isMDCell = true;
184 : for (int d = 0; d < 3; d++) {
185 : isMDCell = isMDCell && (globalCoords[d] > _offset[d] + (int)overlapStrip - 1) &&
186 : (globalCoords[d] < _offset[d] + _globalNumberCouplingCells[d] - (int)overlapStrip);
187 : }
188 : if (isMDCell) {
189 : _flag[get(x, y, z)] = MD_BOUNDARY;
190 : }
191 : }
192 : }
193 : }
194 : }
195 :
196 : /** @brief applies the values received from the MD-solver within the
197 : * conntinuum solver
198 : * @param md2macroBuffer holds the data from the md solver
199 : * coupling cells*/
200 : virtual void setMDBoundaryValues(coupling::datastructures::FlexibleCellContainer<3>& md2macroBuffer) = 0;
201 :
202 : /** @brief returns the number of process, regards parallel runs
203 : * @returns the number of processes */
204 24 : tarch::la::Vector<3, unsigned int> getNumberProcesses() const { return _processes; }
205 :
206 : /** @brief returns the average number of cells on each process
207 : * @returns the average number of cells */
208 8 : tarch::la::Vector<3, unsigned int> getAvgNumberLBCells() const {
209 8 : tarch::la::Vector<3, unsigned int> avgCells(_avgDomainSizeX, _avgDomainSizeY, _avgDomainSizeZ);
210 8 : return avgCells;
211 : }
212 :
213 : /** @brief returns the density for a given position
214 : * @param pos position for which the density will be returned
215 : * @returns the density */
216 : virtual double getDensity(tarch::la::Vector<3, double> pos) const = 0;
217 :
218 : /** @brief changes the velocity at the moving wall (z=0)
219 : * @param wallVelocity new wall velocity to apply */
220 : virtual void setWallVelocity(const tarch::la::Vector<3, double> wallVelocity) = 0;
221 :
222 : /** determines the "avg" domain size which is the domain size on each MPI
223 : * process, except for potentially the last one (the last one may include
224 : * additional cells) */
225 348 : static int getAvgDomainSize(double channelheight, double dx, tarch::la::Vector<3, unsigned int> processes, int d) {
226 348 : int globalDomainSize = floor((channelheight + 0.5) / dx);
227 348 : return globalDomainSize / processes[d];
228 : }
229 :
230 : private:
231 : /** @brief determines the process coordinates
232 : * @returns the coordinates of the current process */
233 116 : tarch::la::Vector<3, unsigned int> getProcessCoordinates() const {
234 116 : tarch::la::Vector<3, unsigned int> coords(0);
235 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
236 116 : constexpr int dim = 3;
237 116 : unsigned int rank = IDXS.getRank();
238 : // determine rank coordinates
239 116 : coords[2] = rank / (_processes[0] * _processes[1]);
240 116 : coords[1] = (rank - coords[2] * _processes[0] * _processes[1]) / _processes[0];
241 116 : coords[0] = rank - coords[2] * _processes[0] * _processes[1] - coords[1] * _processes[0];
242 : #endif
243 116 : return coords;
244 : }
245 :
246 : /** @brief determines the neighbour relation between the processes*/
247 116 : void determineParallelNeighbours() {
248 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
249 116 : int rank;
250 116 : int size;
251 116 : MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
252 116 : MPI_Comm_size(coupling::indexing::IndexingService<3>::getInstance().getComm(), &size);
253 : // check if enough ranks are available
254 116 : if (_processes[0] * _processes[1] * _processes[2] > (unsigned int)size) {
255 0 : std::cout << "ERROR NumericalSolver::determineParallelNeighbours(): Not "
256 0 : "enough ranks available!"
257 0 : << std::endl;
258 0 : exit(EXIT_FAILURE);
259 : }
260 : // neighbour dependencies based on Couette problem
261 : // left,right: periodic
262 116 : _parallelNeighbours[LEFT] = ((_coords[0] + _processes[0] - 1) % _processes[0]) + _processes[0] * (_coords[1] + _processes[1] * _coords[2]);
263 116 : _parallelNeighbours[RIGHT] = ((_coords[0] + _processes[0] + 1) % _processes[0]) + _processes[0] * (_coords[1] + _processes[1] * _coords[2]);
264 : // back,front: periodic
265 116 : _parallelNeighbours[FRONT] = _coords[0] + _processes[0] * (((_coords[1] + _processes[1] - 1) % _processes[1]) + _processes[1] * _coords[2]);
266 116 : _parallelNeighbours[BACK] = _coords[0] + _processes[0] * (((_coords[1] + _processes[1] + 1) % _processes[1]) + _processes[1] * _coords[2]);
267 : // top: either neighbour or MPI_PROC_NULL
268 116 : if (_coords[2] == _processes[2] - 1) {
269 29 : _parallelNeighbours[TOP] = MPI_PROC_NULL;
270 : } else {
271 87 : _parallelNeighbours[TOP] = _coords[0] + _processes[0] * (_coords[1] + _processes[1] * (_coords[2] + 1));
272 : }
273 : // bottom either neighbour or MPI_PROC_NULL
274 116 : if (_coords[2] == 0) {
275 29 : _parallelNeighbours[BOTTOM] = MPI_PROC_NULL;
276 : } else {
277 87 : _parallelNeighbours[BOTTOM] = _coords[0] + _processes[0] * (_coords[1] + _processes[1] * (_coords[2] - 1));
278 : }
279 : // std::cout << "Parallel neighbours for rank " << rank << ": " <<
280 : // _parallelNeighbours << std::endl;
281 : #endif
282 116 : }
283 :
284 : /** @brief sets parallel boundary flags according to Couette scenario */
285 116 : void setParallelBoundaryFlags() {
286 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
287 : // bottom - moving wall
288 116 : if (_coords[2] != 0) {
289 435 : for (int i = 0; i < (_domainSizeX + 2) * (_domainSizeY + 2); i++) {
290 348 : _flag[i] = PARALLEL_BOUNDARY;
291 : }
292 : }
293 : // top - noslip
294 116 : if (_coords[2] != _processes[2] - 1) {
295 435 : for (int i = (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 1); i < (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2); i++) {
296 348 : _flag[get(i)] = PARALLEL_BOUNDARY;
297 : }
298 : }
299 : // for all other boundaries front,back,left,right, we use parallel
300 : // boundaries. If we send from one processes to itself, this is still the
301 : // same as periodic conditions front - periodic
302 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
303 13340 : for (int x = 0; x < _domainSizeX + 2; x++) {
304 12760 : _flag[get(x, 0, z)] = PARALLEL_BOUNDARY;
305 : }
306 : }
307 : // back - periodic
308 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
309 13340 : for (int x = 0; x < _domainSizeX + 2; x++) {
310 12760 : _flag[get(x, _domainSizeY + 1, z)] = PARALLEL_BOUNDARY;
311 : }
312 : }
313 : // left - periodic
314 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
315 12180 : for (int y = 1; y < _domainSizeY + 1; y++) {
316 11600 : _flag[get(0, y, z)] = PARALLEL_BOUNDARY;
317 : }
318 : }
319 696 : for (int z = 1; z < _domainSizeZ + 1; z++) {
320 12180 : for (int y = 1; y < _domainSizeY + 1; y++) {
321 11600 : _flag[get(_domainSizeX + 1, y, z)] = PARALLEL_BOUNDARY;
322 : }
323 : }
324 : #endif // COUPLING_MD_PARALLEL
325 116 : }
326 :
327 : /** @brief determines the local domain size on this rank where channelheight
328 : * is the domain length in direction d.
329 : * @returns the size of the domain */
330 348 : int getDomainSize(double channelheight, double dx, tarch::la::Vector<3, unsigned int> processes, int d) const {
331 348 : int globalDomainSize = floor((channelheight + 0.5) / dx);
332 348 : tarch::la::Vector<3, unsigned int> coords(0);
333 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
334 348 : int rank;
335 348 : MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
336 : // if this rank is outside the range given by processes: return 0
337 : // -> cannot use method skipRank() here, since _processes may not be
338 : // initialized yet!
339 348 : if ((unsigned int)rank > processes[0] * processes[1] * processes[2] - 1) {
340 : return 0;
341 : }
342 : // determine rank coordinates
343 87 : coords[2] = ((unsigned int)rank) / (processes[0] * processes[1]);
344 87 : coords[1] = (((unsigned int)rank) - coords[2] * processes[0] * processes[1]) / processes[0];
345 87 : coords[0] = ((unsigned int)rank) - coords[2] * processes[0] * processes[1] - coords[1] * processes[0];
346 : #endif
347 : // if this is not the last process along this direction: just return avg.
348 : // number of cells
349 87 : if (coords[d] < processes[d] - 1) {
350 0 : return globalDomainSize / processes[d];
351 : // otherwise: add the cells that have not been distributed yet
352 : } else {
353 87 : return globalDomainSize / processes[d] + globalDomainSize % processes[d];
354 : }
355 : }
356 :
357 : protected:
358 : /** @brief returns i and performs checks in debug mode
359 : * @returns i */
360 : int get(int i) const {
361 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
362 : if (i < 0) {
363 : std::cout << "ERROR NumericalSolver::get(i): i<0!" << std::endl;
364 : exit(EXIT_FAILURE);
365 : }
366 : if (i > (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2)) {
367 : std::cout << "ERROR NumericalSolver::get(i): i>max. Value!" << std::endl;
368 : exit(EXIT_FAILURE);
369 : }
370 : #endif
371 : return i;
372 : }
373 :
374 : /** @brief returns linearized index and performs checks in debug mode
375 : * @returns the linearized index */
376 1705680 : int get(int x, int y, int z) const {
377 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
378 : if (x < 0) {
379 : std::cout << "ERROR NumericalSolver::get(x,y,z): x<0!" << std::endl;
380 : exit(EXIT_FAILURE);
381 : }
382 : if (x > _domainSizeX + 1) {
383 : std::cout << "ERROR NumericalSolver::get(x,y,z): x>max. Value!" << std::endl;
384 : exit(EXIT_FAILURE);
385 : }
386 : if (y < 0) {
387 : std::cout << "ERROR NumericalSolver::get(x,y,z): y<0!" << std::endl;
388 : exit(EXIT_FAILURE);
389 : }
390 : if (y > _domainSizeY + 1) {
391 : std::cout << "ERROR NumericalSolver::get(x,y,z): y>max. Value!" << std::endl;
392 : exit(EXIT_FAILURE);
393 : }
394 : if (z < 0) {
395 : std::cout << "ERROR NumericalSolver::get(x,y,z): z<0!" << std::endl;
396 : exit(EXIT_FAILURE);
397 : }
398 : if (z > _domainSizeZ + 1) {
399 : std::cout << "ERROR NumericalSolver::get(x,y,z): z>max. Value!" << std::endl;
400 : exit(EXIT_FAILURE);
401 : }
402 : #endif
403 1705380 : return x + (_domainSizeX + 2) * (y + (_domainSizeY + 2) * z);
404 : }
405 :
406 : /** @brief returns index in 2D parallel buffer with buffer dimensions
407 : * lengthx+2,lengthy+2. Performs checks in debug mode
408 : * @returns the index in the buffer */
409 1103640 : int getParBuf(int x, int y, int lengthx, int lengthy) const {
410 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
411 : if (x < 0 || x > lengthx + 1) {
412 : std::cout << "ERROR NumericalSolver::getParBuf(...): x out of range!" << std::endl;
413 : exit(EXIT_FAILURE);
414 : }
415 : if (y < 0 || y > lengthy + 1) {
416 : std::cout << "ERROR NumericalSolver::getParBuf(...): y out of range!" << std::endl;
417 : exit(EXIT_FAILURE);
418 : }
419 : #endif
420 1103640 : return x + (lengthx + 2) * y;
421 : }
422 :
423 51 : void plot() const {
424 51 : int rank = 0; // rank in MPI-parallel simulations
425 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
426 51 : MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
427 : #endif
428 51 : std::stringstream ss;
429 102 : ss << _filestem << "_r" << rank << "_c" << _counter;
430 51 : if (_scen != nullptr) {
431 0 : auto ts = _scen->getTimeIntegrationService();
432 0 : if (ts != nullptr) {
433 0 : if (ts->isPintEnabled())
434 0 : ss << "_i" << ts->getIteration();
435 : }
436 : }
437 51 : ss << ".vtk";
438 51 : plot(ss.str());
439 51 : }
440 :
441 : /** @brief create vtk plot if required */
442 51 : void plot(std::string filename) const {
443 : // only plot output if this is the correct timestep
444 51 : if (_plotEveryTimestep < 1) {
445 51 : return;
446 : }
447 0 : if (_counter % _plotEveryTimestep != 0) {
448 : return;
449 : }
450 0 : std::ofstream file(filename.c_str());
451 0 : if (!file.is_open()) {
452 0 : std::cout << "ERROR NumericalSolver::plot(): Could not open file " << filename << "!" << std::endl;
453 0 : exit(EXIT_FAILURE);
454 : }
455 0 : std::stringstream flag, density, velocity;
456 :
457 0 : file << "# vtk DataFile Version 2.0" << std::endl;
458 0 : file << "MaMiCo NumericalSolver" << std::endl;
459 0 : file << "ASCII" << std::endl << std::endl;
460 0 : file << "DATASET STRUCTURED_GRID" << std::endl;
461 0 : file << "DIMENSIONS " << _domainSizeX + 3 << " " << _domainSizeY + 3 << " " << _domainSizeZ + 3 << std::endl; // everything +1 cause of change in index
462 0 : file << "POINTS " << (_domainSizeX + 3) * (_domainSizeY + 3) * (_domainSizeZ + 3) << " float" << std::endl;
463 :
464 0 : flag << "CELL_DATA " << (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2) << std::endl;
465 0 : flag << "SCALARS flag float 1" << std::endl;
466 0 : flag << "LOOKUP_TABLE default" << std::endl;
467 :
468 0 : density << std::setprecision(12);
469 0 : density << "SCALARS density float 1 " << std::endl;
470 0 : density << "LOOKUP_TABLE default" << std::endl;
471 :
472 0 : velocity << std::setprecision(12);
473 0 : velocity << "VECTORS velocity float" << std::endl;
474 : // loop over domain (incl. boundary) and write point coordinates
475 0 : for (int z = -1; z < _domainSizeZ + 2; z++) {
476 0 : for (int y = -1; y < _domainSizeY + 2; y++) {
477 0 : for (int x = -1; x < _domainSizeX + 2; x++) {
478 0 : file << ((int)(_coords[0] * _avgDomainSizeX) + x) * _dx << " " << ((int)(_coords[1] * _avgDomainSizeY) + y) * _dx << " "
479 0 : << ((int)(_coords[2] * _avgDomainSizeZ) + z) * _dx << std::endl;
480 : }
481 : }
482 : }
483 : // loop over domain (incl. boundary)
484 0 : for (int z = 0; z < _domainSizeZ + 1 + 1; z++) {
485 0 : for (int y = 0; y < _domainSizeY + 1 + 1; y++) {
486 0 : for (int x = 0; x < _domainSizeX + 1 + 1; x++) { // CHANGE: start index used to be one
487 0 : const int index = get(x, y, z);
488 : // write information to streams
489 0 : flag << _flag[index] << std::endl;
490 0 : density << _density[index] << std::endl;
491 0 : velocity << _vel[3 * index] << " " << _vel[3 * index + 1] << " " << _vel[3 * index + 2] << std::endl;
492 : }
493 : }
494 : }
495 0 : file << std::endl;
496 0 : file << flag.str() << std::endl << std::endl;
497 0 : flag.str("");
498 0 : file << density.str() << std::endl;
499 0 : density.str("");
500 0 : file << velocity.str() << std::endl;
501 0 : velocity.str("");
502 0 : file.close();
503 0 : }
504 :
505 : /** @brief returns true, if this rank is not of relevance for the LB
506 : * simulation
507 : * @returns a bool, which indicates if the rank shall not do anything (true)
508 : * or not (false) */
509 401 : bool skipRank() const {
510 401 : int rank = 0;
511 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
512 401 : MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
513 : #endif
514 401 : return ((unsigned int)rank > _processes[0] * _processes[1] * _processes[2] - 1);
515 : }
516 :
517 : /** @brief for every cell exists a flag entry, upon this is defined how the
518 : * cell is handled */
519 : enum Flag {
520 : FLUID = 0, ///< @brief a normal fluid cell
521 : NO_SLIP = 1, ///< @brief a cell on the no slip (non-moving) wall
522 : MOVING_WALL = 2, ///< @brief a cell on the moving wall
523 : PERIODIC = 3, ///< @brief a cell on a periodic boundary
524 : MD_BOUNDARY = 4, ///< @brief a cell on the boundary to md
525 : PARALLEL_BOUNDARY = 5 ///< @brief a cell on a inner boundary of a splitted
526 : ///< domain in a parallel run
527 : };
528 : /** @brief The flags are used on parallel boundaries to define in which
529 : * direction the boundary goes */
530 : enum NbFlag {
531 : LEFT = 0, ///< @brief a parallel boundary to the left
532 : RIGHT = 1, ///< @brief a parallel boundary to the right
533 : BACK = 2, ///< @brief a parallel boundary to the back
534 : FRONT = 3, ///< @brief a parallel boundary to the front
535 : BOTTOM = 4, ///< @brief a parallel boundary to the bottom
536 : TOP = 5 ///< @brief a parallel boundary to the top
537 : };
538 : /** @brief the height and width of the channel in z and y direction */
539 : const double _channelheight; //
540 : /** @brief mesh size, dx=dy=dz */
541 : const double _dx;
542 : /** @brief time step*/
543 : const double _dt;
544 : /** @brief kinematic viscosity of the fluid */
545 : const double _kinVisc;
546 : /** @brief domain decomposition on MPI rank basis; total number is given by
547 : * multipling all entries*/
548 : tarch::la::Vector<3, unsigned int> _processes;
549 : /** @brief number of time steps between vtk plots */
550 : const int _plotEveryTimestep;
551 : /** @brief file stem for vtk plot */
552 : const std::string _filestem;
553 : /** @brief domain size in x-direction */
554 : const int _domainSizeX{getDomainSize(_channelheight, _dx, _processes, 0)};
555 : /** @brief domain size in y-direction */
556 : const int _domainSizeY{getDomainSize(_channelheight, _dx, _processes, 1)};
557 : /** @brief domain size in z-direction */
558 : const int _domainSizeZ{getDomainSize(_channelheight, _dx, _processes, 2)};
559 : /** @brief avg. domain size in MPI-parallel simulation in x-direction */
560 : const int _avgDomainSizeX{getAvgDomainSize(_channelheight, _dx, _processes, 0)};
561 : /** @brief avg. domain size in MPI-parallel simulation in y-direction */
562 : const int _avgDomainSizeY{getAvgDomainSize(_channelheight, _dx, _processes, 1)}; //
563 : /** @brief avg. domain size in MPI-parallel simulation in z-direction */
564 : const int _avgDomainSizeZ{getAvgDomainSize(_channelheight, _dx, _processes, 2)}; //
565 : /** @brief coordinates of this process (=1,1,1, unless parallel run of the
566 : * solver )*/
567 : const tarch::la::Vector<3, unsigned int> _coords{getProcessCoordinates()};
568 : /** @brief time step counter */
569 : int _counter{0};
570 : /** @brief velocity field */
571 : double* _vel{NULL};
572 : /** @brief density field */
573 : double* _density{NULL};
574 : /** @brief flag field */
575 : Flag* _flag{NULL};
576 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
577 : /** @brief buffer to send data from left/right to right/left neighbour */
578 : double* _sendBufferX{NULL};
579 : /** @brief buffer to receive data from from left/right neighbour */
580 : double* _recvBufferX{NULL};
581 : /** @brief buffer to send data from front/back to front/back neighbour */
582 : double* _sendBufferY{NULL};
583 : /** @brief buffer to receive data from from front/back neighbour */
584 : double* _recvBufferY{NULL};
585 : /** @brief buffer to send data from top/buttom to top/buttom neighbour */
586 : double* _sendBufferZ{NULL};
587 : /** @brief buffer to receive data from from top/buttom neighbour */
588 : double* _recvBufferZ{NULL};
589 : #endif
590 : /** @brief offset for y-direction (lexicographic grid ordering) */
591 : const int _xO{_domainSizeX + 2};
592 : /** @brief offset for z-direction */
593 : const int _yO{(_domainSizeX + 2) * (_domainSizeY + 2)};
594 : /** @brief neighbour ranks */
595 : tarch::la::Vector<6, int> _parallelNeighbours{(-1)};
596 : /** @brief offset of the md domain */
597 : tarch::la::Vector<3, int> _offset{(-1)};
598 : /** @brief the total number of coupling cells of the coupled simulation */
599 : tarch::la::Vector<3, int> _globalNumberCouplingCells{(-1)};
600 : const Scenario* _scen;
601 : };
602 :
603 : #endif // _MOLECULARDYNAMICS_COUPLING_SOLVERS_NUMERICALSOLVER_H_
|
