OpenLB 1.7
Loading...
Searching...
No Matches
random.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2021 Nicolas Hafen, Mathias J. Krause
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24#ifndef RANDOM_HH
25#define RANDOM_HH
26
27
28#include <type_traits>
29
30namespace olb {
31
32namespace util {
33
34template<typename T,bool useStored>
36{
37 static_assert(!useStored, "ERROR: no sequence provided");
38}
39
40template<typename T,bool useStored>
42 : _storedSequence(sequence)
43{}
44
45
46template<typename T,bool useStored>
47Randomizer<T,useStored>::Randomizer(std::string filePathSequence,
48 bool enforceStored)
49{
50 [[maybe_unused]] std::size_t count = 0;
51 int rank = singleton::mpi().getRank();
52 if (rank == 0){
53 //Check whether file exists and create, if not
54 std::ifstream filePre(filePathSequence.c_str(), std::ios::in);
55 if (!filePre.good()){
56 if (!enforceStored){
57 writeSequence(1000,filePathSequence.c_str());
58 } else {
59 throw std::runtime_error("ERROR: Sequence not provided for Randomizer!");
60 }
61 }
62 //Read actual file
63 std::ifstream file(filePathSequence.c_str(), std::ios::in);
64 std::string strVal;
65 while (std::getline(file, strVal))
66 {
67 _storedSequence.push_back(std::stod(strVal));
68 ++count;
69 }
70 }
71
72#ifdef PARALLEL_MODE_MPI
73 //Align count on all cores
74 singleton::mpi().bCast<std::size_t>( &count, 1 );
75
76 //Align vector size on all cores
77 if (rank != 0){
78 _storedSequence.resize(count);
79 }
80
81 //Align vector content on all cores
82 singleton::mpi().bCast( _storedSequence.data(), count );
83#endif
84}
85
86
87
88template<typename T,bool useStored>
89template<typename O>
91{
92 if constexpr (std::is_arithmetic<O>::value) {
93 O output = generateScalar();
94 return output;
95 } else {
96 O output;
97 for(unsigned iDim=0; iDim<O::d; ++iDim){
98 output[iDim] = generateScalar();
99 }
100 return output;
101 }
102}
103
104template<typename T,bool useStored>
106{
107 T A = generateScalar();
108 T B = generateScalar();
109 T x = util::cos(2 * M_PI * A) * util::sqrt(-2 * util::log(B));
110 T output = avg + x * stdDev;
111 return output;
112}
113
114template<typename T,bool useStored>
116{
117 T scalar = generateScalarNormal( avg, stdDev );
118 if (scalar<(avg-cutoff) || scalar>(avg+cutoff)){
119 scalar = generateScalarNormal(avg, stdDev, cutoff );
120 }
121 return scalar;
122}
123
124template<typename T,bool useStored>
125void Randomizer<T,useStored>::writeSequence( std::size_t numOfValues,
126 std::string filePathSequence,
127 int precision )
128{
129 int rank = singleton::mpi().getRank();
130 if (rank == 0){
131 std::ofstream writer;
132 writer.open( filePathSequence );
133 writer << std::setprecision(precision);
134 for (std::size_t i=0; i<numOfValues; ++i){
135 writer << generateScalarRandom() << " " << std::endl;
136 }
137 writer.close();
138 }
139}
140
141
142template<typename T,bool useStored>
144{
145 // Non-deterministic random number generator
146 std::random_device rd;
147 // Pseudo-random number engine: Mersenne Twister 19937 generator
148 std::mt19937 engine(rd());
149
150 // Distribution to generate numbers in the range [0.0, 1.0)
151 std::uniform_real_distribution<T> distribution(0.0, 1.0);
152
153 // Generate a random number
154 return distribution(engine);
155}
156
157template<typename T,bool useStored>
158T Randomizer<T,useStored>::generateScalarStored()
159{
160 static_assert(useStored, "ERROR: no sequence provided");
161 auto idx = _count%_storedSequence.size();
162 T scalar = _storedSequence[idx];
163 ++_count;
164 return scalar;
165}
166
167template<typename T,bool useStored>
168T Randomizer<T,useStored>::generateScalar()
169{
170 if constexpr (useStored){
171 return generateScalarStored();
172 } else {
173 T scalar = generateScalarRandom();
174#ifdef PARALLEL_MODE_MPI
175 singleton::mpi().bCast<T>( &scalar, 1 );
176#endif
177 return scalar;
178 }
179}
180
181} //namespace util
182
183} //namespace olb
184
185#endif
#define M_PI
void bCast(T *sendBuf, int sendCount, int root=0, MPI_Comm comm=MPI_COMM_WORLD)
Broadcast data from one processor to multiple processors.
int getRank() const
Returns the process ID.
void writeSequence(std::size_t numOfValues, std::string filePathSequence="./randomSequence.dat", int precision=5)
Write sequence to file for later retrieval.
Definition random.hh:125
T generateScalarNormal(T avg, T stdDev)
Generate scalar leading to normal distribution based on the Box-Muller approach.
Definition random.hh:105
O generate()
Generate scalar or vector filled with scalars.
Definition random.hh:90
Randomizer()
Constructor for (useStored=false)
Definition random.hh:35
MpiManager & mpi()
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
ADf< T, DIM > log(const ADf< T, DIM > &a)
Definition aDiff.h:475
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Definition aDiff.h:578
Top level namespace for all of OpenLB.