diff --git a/CpuImplementations/include/GeliosphereCpuModel.hpp b/CpuImplementations/include/GeliosphereCpuModel.hpp index 8050976..0ff1154 100644 --- a/CpuImplementations/include/GeliosphereCpuModel.hpp +++ b/CpuImplementations/include/GeliosphereCpuModel.hpp @@ -1,24 +1,13 @@ -/** - * @file GeliosphereCpuModel.hpp - * @author Michal Solanik - * @brief CPU implementation for Geliosphere 2D model - * @version 0.2 - * @date 2022-07-07 - * - * @copyright Copyright (c) 2022 - * - */ - #ifndef BP_CPU_THREE_DIMENSION_SIMULATION_H #define BP_CPU_THREE_DIMENSION_SIMULATION_H #include "AbstractCpuModel.hpp" +#include "IGeliosphereCpuModel.hpp" #include #include -struct SimulationOutput -{ +struct SimulationOutput { double Tkininj; double Tkin; double r; @@ -27,21 +16,20 @@ struct SimulationOutput double theta; }; -/** - * @brief Class implements @ref AbstractAlgorithm "AbstractAlgorithm" interface - * to define support functions for running implementation of Geliosphere 2D B-p model. - * - */ -class GeliosphereCpuModel : public AbstractCpuModel -{ +class GeliosphereCpuModel : public AbstractCpuModel, public IGeliosphereCpuModel { public: - /** - * @brief Definition of simulation runner. - * - * @param singleTone data structure containing input parameters. - */ void runSimulation(ParamsCarrier *singleTone); +protected: + double Beta(double Tkin) override; + double RigFromTkin(double Tkin) override; + double W(double p) override; + double GetMomentum(double Rig) override; + double LarmorRadius(double Rig, double Bfield) override; + double AlphaH(double Larmor, double r) override; + double ComputeF(double theta, double alphaH) override; + double ComputeFPrime(double theta, double alphaH) override; + private: void simulation(int threadNumber, unsigned int availableThreads, int iteration); std::mutex outputMutex; diff --git a/CpuImplementations/include/IGeliosphereCpuModel.hpp b/CpuImplementations/include/IGeliosphereCpuModel.hpp new file mode 100644 index 0000000..610f83f --- /dev/null +++ b/CpuImplementations/include/IGeliosphereCpuModel.hpp @@ -0,0 +1,19 @@ +#ifndef I_GELIOSPHERE_CPU_MODEL_H +#define I_GELIOSPHERE_CPU_MODEL_H + +class IGeliosphereCpuModel { +protected: + virtual double Beta(double Tkin) = 0; + virtual double RigFromTkin(double Tkin) = 0; + virtual double W(double p) = 0; + virtual double GetMomentum(double Rig) = 0; + virtual double LarmorRadius(double Rig, double Bfield) = 0; + virtual double AlphaH(double Larmor, double r) = 0; + virtual double ComputeF(double theta, double alphaH) = 0; + virtual double ComputeFPrime(double theta, double alphaH) = 0; + +public: + virtual ~IGeliosphereCpuModel() = default; +}; + +#endif // I_GELIOSPHERE_CPU_MODEL_H diff --git a/CpuImplementations/include/IOneDimensionBpCpuModel.hpp b/CpuImplementations/include/IOneDimensionBpCpuModel.hpp new file mode 100644 index 0000000..bb79b8b --- /dev/null +++ b/CpuImplementations/include/IOneDimensionBpCpuModel.hpp @@ -0,0 +1,17 @@ +#ifndef I_ONE_DIMENSION_BP_CPU_MODEL_H +#define I_ONE_DIMENSION_BP_CPU_MODEL_H + +class IOneDimensionBpCpuModel { +protected: + virtual double Beta(double Tkin) = 0; + virtual double RigFromTkin(double Tkin) = 0; + virtual double Kdiffr(double beta, double Rig) = 0; + virtual double Dp(double V, double p, double r) = 0; + virtual double Dr(double V, double Kdiff, double r, double dt, double rand) = 0; + virtual double W(double p) = 0; + +public: + virtual ~IOneDimensionBpCpuModel() = default; +}; + +#endif // I_ONE_DIMENSION_BP_CPU_MODEL_H diff --git a/CpuImplementations/include/IOneDimensionFpCpuModel.hpp b/CpuImplementations/include/IOneDimensionFpCpuModel.hpp new file mode 100644 index 0000000..b0ef0de --- /dev/null +++ b/CpuImplementations/include/IOneDimensionFpCpuModel.hpp @@ -0,0 +1,20 @@ +#ifndef I_ONE_DIMENSION_FP_CPU_MODEL_H +#define I_ONE_DIMENSION_FP_CPU_MODEL_H + +class IOneDimensionFpCpuModel { +protected: + virtual double Beta(double Tkin) = 0; + virtual double RigFromTkin(double Tkin) = 0; + virtual double RigFromTkinJoule(double Tkin) = 0; + virtual double Kdiffr(double beta, double Rig) = 0; + virtual double Dp(double V, double p, double r) = 0; + virtual double Dr(double V, double Kdiff, double r, double dt, double rand) = 0; + virtual double Cfactor(double V, double r) = 0; + virtual double TkinFromRig(double Rig) = 0; + virtual double W(double p) = 0; + +public: + virtual ~IOneDimensionFpCpuModel() = default; +}; + +#endif // I_ONE_DIMENSION_FP_CPU_MODEL_H diff --git a/CpuImplementations/include/OneDimensionBpCpuModel.hpp b/CpuImplementations/include/OneDimensionBpCpuModel.hpp index 000ee4a..eac8aa0 100644 --- a/CpuImplementations/include/OneDimensionBpCpuModel.hpp +++ b/CpuImplementations/include/OneDimensionBpCpuModel.hpp @@ -13,6 +13,7 @@ #define BP_CPU_ONE_DIMENSION_SIMULATION_H #include "AbstractCpuModel.hpp" +#include "IOneDimensionBpCpuModel.hpp" #include #include @@ -22,6 +23,10 @@ struct SimulationOutput { double Tkininj; double r; double w; + + // Constructor + //SimulationOutput(double Tkin_, double Tkininj_, double r_, double w_) + //: Tkin(Tkin_), Tkininj(Tkininj_), r(r_), w(w_) {} }; /** @@ -29,7 +34,7 @@ struct SimulationOutput { * to define support functions for running implementation of 1D B-p model. * */ -class OneDimensionBpCpuModel : public AbstractCpuModel +class OneDimensionBpCpuModel : public AbstractCpuModel, public IOneDimensionBpCpuModel { public: /** @@ -38,6 +43,18 @@ class OneDimensionBpCpuModel : public AbstractCpuModel * @param singleTone data structure containing input parameters. */ void runSimulation(ParamsCarrier *singleTone); + +protected: + // Implemented from interface + double Beta(double Tkin) override; + double RigFromTkin(double Tkin) override; + //double RigFromMomentum(double p) override; + double Kdiffr(double beta, double Rig) override; + double Dp(double V, double p, double r) override; + double Dr(double V, double Kdiff, double r, double dt, double rand) override; + //double TkinFromRig(double Rig) override; + double W(double p) override; + private: void simulation(int threadNumber, unsigned int availableThreads, int iteration); std::mutex outputMutex; diff --git a/CpuImplementations/include/OneDimensionFpCpuModel.hpp b/CpuImplementations/include/OneDimensionFpCpuModel.hpp index 0c7bfcc..b99d9dd 100644 --- a/CpuImplementations/include/OneDimensionFpCpuModel.hpp +++ b/CpuImplementations/include/OneDimensionFpCpuModel.hpp @@ -13,6 +13,7 @@ #define FP_CPU_ONE_DIMENSION_SIMULATION_H #include "AbstractCpuModel.hpp" +#include "IOneDimensionFpCpuModel.hpp" #include #include @@ -30,7 +31,7 @@ struct SimulationOutput { * to define support functions for running implementation of 1D F-p model. * */ -class OneDimensionFpCpuModel : public AbstractCpuModel +class OneDimensionFpCpuModel : public AbstractCpuModel, public IOneDimensionFpCpuModel { public: /** @@ -39,6 +40,19 @@ class OneDimensionFpCpuModel : public AbstractCpuModel * @param singleTone data structure containing input parameters. */ void runSimulation(ParamsCarrier *singleTone); + +protected: + // Implemented from interface + double Beta(double Tkin) override; + double RigFromTkin(double Tkin) override; + double RigFromTkinJoule(double Tkin) override; + double Kdiffr(double beta, double Rig) override; + double Dp(double V, double p, double r) override; + double Dr(double V, double Kdiff, double r, double dt, double rand) override; + double Cfactor(double V, double r) override; + double TkinFromRig(double Rig) override; + double W(double p) override; + private: void simulation(int threadNumber, unsigned int availableThreads, int iteration); std::mutex outputMutex; diff --git a/CpuImplementations/src/GeliosphereCpuModel.cpp b/CpuImplementations/src/GeliosphereCpuModel.cpp index 5d723ec..e3443ac 100644 --- a/CpuImplementations/src/GeliosphereCpuModel.cpp +++ b/CpuImplementations/src/GeliosphereCpuModel.cpp @@ -65,46 +65,38 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr double Kphph, Krt, Krph, Ktph, B11Temp, B11, B12, B13, B22, B23; double sin2, sin3, dKtt1, dKtt2; double dKrtr, dKrtt; + thread_local std::random_device rd{}; thread_local std::mt19937 generator(rd()); thread_local std::normal_distribution distribution(0.0f, 1.0f); + std::vector localOutputs; + localOutputs.reserve(30 * 500); + for (int energy = 0; energy < 30; energy++) { for (int particlePerEnergy = 0; particlePerEnergy < 500; particlePerEnergy++) { - Tkininj = (useUniformInjection) - ? getTkinInjection(((availableThreads * iteration + threadNumber) * 500) + particlePerEnergy, 0.0001, uniformEnergyInjectionMaximum, 10000) - : SPbins[energy]; - Tkin = Tkininj; + Tkininj = (useUniformInjection) + ? getTkinInjection(((availableThreads * iteration + threadNumber) * 500) + particlePerEnergy, 0.0001, uniformEnergyInjectionMaximum, 10000) + : SPbins[energy]; + Tkin = Tkininj; Tkinw = Tkin * 1e9 * q; Rig = sqrt(Tkinw * (Tkinw + (2.0 * T0w))) / q; - p = Rig * q / c; - - // Equation under Equation 6 from - // Yamada et al. "A stochastic view of the solar modulation phenomena of cosmic rays" GEOPHYSICAL RESEARCH LETTERS, VOL. 25, NO.13, PAGES 2353-2356, JULY 1, 1998 - // https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/98GL51869 - w = (m0 * m0 * c * c * c * c) + (p * p * c * c); - w = pow(w, -1.85) / p; - w = w / 1e45; + p = GetMomentum(Rig); + w = W(p); r = rInit; theta = thetainj; while (r < 100.0) { - // Equation 5 - beta = sqrtf(Tkin * (Tkin + T0 + T0)) / (Tkin + T0); - - // Equation 8 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf - Rig = sqrt(Tkin * (Tkin + (2.0 * T0))); + beta = Beta(Tkin); + Rig = RigFromTkin(Tkin); r2 = r * r; - // Equation 25 if (theta < (1.7 * Pi / 180.) || theta > (178.3 * Pi / 180.0)) { delta = 0.003; @@ -117,16 +109,13 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr deltarh = delta / rh; deltarh2 = deltarh * deltarh; - // Equation 24 gamma = (r * omega) * sin(theta) / V; gamma2 = gamma * gamma; - // Equation 26 Cb = 1.0 + gamma2 + (r2 * deltarh2); Cb2 = Cb * Cb; Bfactor = (5. / 3.4) * r2 / sqrt(Cb); - // Equation 32 if (Rig < 0.1) { Kpar = K0 * beta * 0.1 * Bfactor / 3.0; @@ -136,65 +125,43 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr Kpar = K0 * beta * Rig * Bfactor / 3.0; } - // Equation 33 Kper = ratio * Kpar; - // Equation 27 Krr = Kper + ((Kpar - Kper) / Cb); - - // Equation 28 Ktt = Kper + (r2 * deltarh2 * (Kpar - Kper) / Cb); Kphph = 1.0; - // Equation 29 Krt = deltarh * (Kpar - Kper) * r / Cb; Krph = 0.0; Ktph = 0.0; - // Equation 16, where Krph = Ktph = 0, and Kphph = 1 B11Temp = (Kphph * Krt * Krt) - (2.0 * Krph * Krt * Ktph) + (Krr * Ktph * Ktph) + (Ktt * Krph * Krph) - (Krr * Ktt * Kphph); B11 = 2.0 * B11Temp / ((Ktph * Ktph) - (Ktt * Kphph)); B11 = sqrt(B11); - // Equation 17 B12 = ((Krph * Ktph) - (Krt * Kphph)) / ((Ktph * Ktph) - (Ktt * Kphph)); B12 = B12 * sqrt(2.0 * (Ktt - (Ktph * Ktph / Kphph))); - // Equation 20 B13 = sqrt(2.0) * Krph / sqrt(Kphph); - // Equation 18 B22 = Ktt - (Ktph * Ktph / Kphph); B22 = sqrt(2.0 * B22) / r; - // Equation 20 B23 = Ktph * sqrt(2.0 / Kphph) / r; - // Equation 34 COmega = 2.0 * r * omega * omega * sin(theta) * sin(theta) / (V * V); COmega = COmega + (2.0 * r * deltarh2); - - // Equation 36 + dKper = ratio * K0 * beta * Rig * ((2.0 * r * sqrt(Cb)) - (r2 * COmega / (2.0 * sqrt(Cb)))) / (3.0 * (5.0 / 3.4) * Cb); - // Equation 35 - dKrr = dKper + ((1.0 - ratio) * K0 * beta * Rig * ((2.0 * r * pow(Cb, 1.5)) - (r2 * COmega * 3.0 * sqrt(Cb) / 2.0)) / ( 3.0 * (5.0 / 3.4) * pow(Cb, 3.0))); + dKrr = dKper + ((1.0 - ratio) * K0 * beta * Rig * ((2.0 * r * pow(Cb, 1.5)) - (r2 * COmega * 3.0 * sqrt(Cb) / 2.0)) / (3.0 * (5.0 / 3.4) * pow(Cb, 3.0))); if ((theta > (1.7 * Pi / 180.)) && (theta < (178.3 * Pi / 180.0))) { - // Equation 37 CKtt = sin(theta) * cos(theta) * (omega * omega * r2 / (V * V)); - - // Equation 38 dKtt1 = (-1.0 * ratio * K0 * beta * Rig * r2 * CKtt) / (3.0 * (5.0 / 3.4) * pow(Cb, 1.5)); - - // Equation 39 dKtt2 = (1.0 - ratio) * K0 * beta * Rig * r2 * r2 * deltarh2; - - // Equation 41 dKtt4 = 3.0 * CKtt / pow(Cb, 2.5); - - // Equation 42 dKtt = dKtt1 - (dKtt2 * dKtt4); } else @@ -202,30 +169,17 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr sin2 = sin(theta) * sin(theta); sin3 = sin(theta) * sin(theta) * sin(theta); - // Equation 37 CKtt = sin(theta) * cos(theta) * (omega * omega * r2 / (V * V)); CKtt = CKtt - (r2 * delta0 * delta0 * cos(theta) / (rh * rh * sin3)); - - // Equation 38 dKtt1 = (-1.0 * ratio * K0 * beta * Rig * r2 * CKtt) / (3.0 * (5.0 / 3.4) * pow(Cb, 1.5)); - - // Equation 39 dKtt2 = (1.0 - ratio) * K0 * beta * Rig * r2 * r2 * delta0 * delta0 / (rh * rh); - - // Equation 40 dKtt3 = -2.0 * (cos(theta) / sin3) / (3.0 * (5.0 / 3.4) * pow(Cb, 1.5)); - - // Equation 41 dKtt4 = 3.0 * (CKtt / sin2) / (3.0 * (5.0 / 3.4) * pow(Cb, 2.5)); - - // Equation 42 dKtt = dKtt1 + (dKtt2 * (dKtt3 - dKtt4)); } - // Equation 43 dKrtr = (1.0 - ratio) * K0 * beta * Rig * deltarh * r2 / (3.0 * (5.0 / 3.4) * pow(Cb, 2.5)); - // Equation 44 if ((theta > (1.7 * Pi / 180.)) && (theta < (178.3 * Pi / 180.0))) { dKrtt = (1.0 - ratio) * K0 * beta * Rig * r2 * r / ((3.0 * (5.0 / 3.4) * rh * pow(Cb, 2.5))); @@ -239,70 +193,31 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr dKrtt = dKrtt * (1.0 - (2.0 * r2 * deltarh2) + (4.0 * gamma2)); } - // Equation 21 dr = ((-1.0 * V) + (2.0 * Krr / r) + dKrr) * dt; dr = dr + (dKrtt * dt / r) + (Krt * cos(theta) * dt / (r * sin(theta))); dr = dr + (distribution(generator) * B11 * sqrt(dt)); dr = dr + (distribution(generator) * B12 * sqrt(dt)); dr = dr + (distribution(generator) * B13 * sqrt(dt)); - // Equation 22 dtheta = (Ktt * cos(theta)) / (r2 * sin(theta)); dtheta = (dtheta * dt) + (dKtt * dt / r2); dtheta = dtheta + (dKrtr * dt) + (2.0 * Krt * dt / r); dtheta = dtheta + (distribution(generator) * B22 * sqrt(dt)) + (distribution(generator) * B23 * sqrt(dt)); - // Equation 23 - alfa = (Tkin + T0 + T0)/(Tkin + T0); + alfa = (Tkin + T0 + T0)/(Tkin + T0); dTkin = -2.0 * V * alfa * Tkin * dt / (3.0 * r); Bfield = A * sqrt(Cb) / (r * r); - - // Equation 26 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf - Larmor = 0.0225 * Rig / Bfield; - - // Equation 34 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf - alphaH = Pi / sin(alphaM + (2.0 * Larmor * Pi / (r * 180.0))); - alphaH = alphaH - 1.0; - alphaH = 1.0 / alphaH; - alphaH = acos(alphaH); - - // Equation 32 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf - arg = (1. - (2. * theta / Pi)) * tan(alphaH); - f = atan(arg) / alphaH; - - // Equation 35 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf + Larmor = LarmorRadius(Rig, Bfield); + alphaH = AlphaH(Larmor, r); + f = ComputeF(theta, alphaH); + fprime = ComputeFPrime(theta, alphaH); + DriftR = polarity * konvF * (2.0 / (3.0 * A)) * Rig * beta * r * cos(theta) * gamma * f / (Cb2 * sin(theta)); - - // Equation 36 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf DriftTheta = -1.0 * polarity * konvF * (2.0 / (3.0 * A)) * Rig * beta * r * gamma * (2.0 + (gamma * gamma)) * f / Cb2; - - // Equation 33 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf - fprime = 1.0 + (arg * arg); - fprime = tan(alphaH) / fprime; - fprime = -1.0 * fprime * 2.0 / (Pi * alphaH); - - // Equation 37 from - // Kappl, Rolf. "SOLARPROP: Charge-sign dependent solar modulation for everyone." Computer Physics Communications 207 (2016): 386-399. - // https://arxiv.org/pdf/1511.07875.pdf DriftSheetR = polarity * konvF * (1.0 / (3.0 * A)) * Rig * beta * r * gamma * fprime / Cb; - // Equation 21 dr = dr + ((DriftR + DriftSheetR) * dt); - - // Equation 22 dtheta = dtheta + (DriftTheta * dt / r); r = r + dr; @@ -329,12 +244,58 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr if (r > 100.0) { - outputMutex.lock(); - outputQueue.push({Tkininj, Tkin, r, w, thetainj, theta}); - outputMutex.unlock(); + + localOutputs.emplace_back(SimulationOutput{Tkininj, Tkin, r, w, thetainj, theta}); + break; } } } } + std::lock_guard lock(outputMutex); + for (const auto& output : localOutputs) + { + outputQueue.push(output); + } +} + + + + +double GeliosphereCpuModel::Beta(double Tkin) { + return sqrt(Tkin * (Tkin + 2.0 * T0)) / (Tkin + T0); +} + +double GeliosphereCpuModel::RigFromTkin(double Tkin) { + return sqrt(Tkin * (Tkin + 2.0 * T0)); +} + +double GeliosphereCpuModel::W(double p) { + double w = (m0 * m0 * c * c * c * c) + (p * p * c * c); + return pow(w, -1.85) / p / 1e45; +} + +double GeliosphereCpuModel::GetMomentum(double Rig) { + return Rig * q / c; +} + +double GeliosphereCpuModel::LarmorRadius(double Rig, double Bfield) { + return 0.0225 * Rig / Bfield; +} + +double GeliosphereCpuModel::AlphaH(double Larmor, double r) { + double alphaH = Pi / sin(alphaM + (2.0 * Larmor * Pi / (r * 180.0))); + alphaH = 1.0 / (alphaH - 1.0); + return acos(alphaH); +} + +double GeliosphereCpuModel::ComputeF(double theta, double alphaH) { + double arg = (1. - (2. * theta / Pi)) * tan(alphaH); + return atan(arg) / alphaH; +} + +double GeliosphereCpuModel::ComputeFPrime(double theta, double alphaH) { + double arg = (1. - (2. * theta / Pi)) * tan(alphaH); + double fprime = tan(alphaH) / (1.0 + arg * arg); + return -1.0 * fprime * 2.0 / (Pi * alphaH); } \ No newline at end of file diff --git a/CpuImplementations/src/OneDimensionBpCpuModel.cpp b/CpuImplementations/src/OneDimensionBpCpuModel.cpp index 945a27a..a160f00 100644 --- a/CpuImplementations/src/OneDimensionBpCpuModel.cpp +++ b/CpuImplementations/src/OneDimensionBpCpuModel.cpp @@ -7,6 +7,8 @@ #include #include + + void OneDimensionBpCpuModel::runSimulation(ParamsCarrier *singleTone) { srand(time(NULL)); @@ -53,13 +55,28 @@ void OneDimensionBpCpuModel::runSimulation(ParamsCarrier *singleTone) void OneDimensionBpCpuModel::simulation(int threadNumber, unsigned int availableThreads, int iteration) { - double r, dr, arnum; - double Tkin, Tkininj, Rig, beta; - double w; - double p, dp, pp, Kdiff; + double r = 0, dr = 0, arnum = 0; + double Tkin = 0, Tkininj = 0, Rig = 0, beta = 0; + double w = 0; + double p = 0, dp = 0, pp = 0, Kdiff = 0; + + + const double C_Q = c / q; + const double C_Q_1E9 = C_Q / 1e9; + const double t0Term = T0 * T0 * q * q * 1e18; + const double q1e9 = q * 1e9; + thread_local std::random_device rd{}; thread_local std::mt19937 generator(rd()); thread_local std::normal_distribution distribution(0.0f, 1.0f); + + + auto& gen = generator; + auto& dist = distribution; + + std::vector localOutputs; + localOutputs.reserve(101 * 250); + for (int energy = 0; energy < 101; energy++) { for (int particlePerEnergy = 0; particlePerEnergy < 250; particlePerEnergy++) @@ -67,64 +84,92 @@ void OneDimensionBpCpuModel::simulation(int threadNumber, unsigned int available Tkininj = getTkinInjection(((availableThreads * iteration + threadNumber) * 250) + particlePerEnergy, 0.0001, uniformEnergyInjectionMaximum, 10000); Tkin = Tkininj; - Rig = sqrt(Tkin * (Tkin + (2.0 * T0))); - p = Rig * 1e9 * q / c; + Rig = RigFromTkin(Tkin); + + p = Rig / C_Q_1E9; r = rInit; while (r < 100.0002) { - // Equation 5 - beta = sqrtf(Tkin * (Tkin + T0 + T0)) / (Tkin + T0); - - // Equation 6 in GeV - Rig = (p * c / q) / 1e9; - pp = p; - - // Equation 14 - dp = (2.0f * V * pp * dt / (3.0f * r)); + beta = Beta(Tkin); + + Rig = p * C_Q_1E9; + pp = p; + dp = Dp(V, pp, r); p -= dp; - // Equation 7 - Kdiff = K0 * beta * Rig; - arnum = distribution(generator); - - // Equation 13 - dr = (V + (2.0 * Kdiff / r)) * dt + (arnum * sqrt(2.0 * Kdiff * dt)); - r += dr; - - // Equation 6 in J - Rig = p * c / q; - Tkin = (sqrt((T0 * T0 * q * q * 1e9 * 1e9) + (q * q * Rig * Rig)) - (T0 * q * 1e9)) / (q * 1e9); - - // Equation 5 - beta = sqrtf(Tkin * (Tkin + T0 + T0)) / (Tkin + T0); + Kdiff = Kdiffr(beta, Rig); + arnum = dist(gen); + + dr = Dr(V, Kdiff, r, dt, arnum); + r += dr; + + + Rig = p * C_Q; + + Tkin = (sqrt(t0Term + (q * q * Rig * Rig)) - (T0 * q1e9)) / q1e9; + beta = Beta(Tkin); if (beta > 0.01f && Tkin < 200.0f) - { - if ((r > 100.0f) && ((r - dr) < 100.0f)) - { - // Equation under Equation 6 from - // Yamada et al. "A stochastic view of the solar modulation phenomena of cosmic rays" GEOPHYSICAL RESEARCH LETTERS, VOL. 25, NO.13, PAGES 2353-2356, JULY 1, 1998 - // https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/98GL51869 - w = (m0 * m0 * c * c * c * c) + (p * p * c * c); - w = (pow(w, -1.85) / p) / 1e45; - - outputMutex.lock(); - outputQueue.push({Tkin, Tkininj, r, w}); - outputMutex.unlock(); - break; + { + if ((r > 100.0f) && ((r - dr) < 100.0f)) + { + w = W(p); + + localOutputs.emplace_back(SimulationOutput{Tkin, Tkininj, r, w}); + + break; } } - else if (beta < 0.01f) - { - break; - } - if (r < 0.3f) - { - r -= dr; - p = pp; - } - } + else if (beta < 0.01f) + { + break; + } + + if (r < 0.3f) + { + r -= dr; + p = pp; + } + } } - } -} \ No newline at end of file + } + { + std::lock_guard lock(outputMutex); + for (const auto& output : localOutputs) + { + outputQueue.push(output); + } + } +} + +double OneDimensionBpCpuModel::Beta(double Tkin) +{ + return sqrt(Tkin * (Tkin + 2.0 * T0)) / (Tkin + T0); +} + +double OneDimensionBpCpuModel::RigFromTkin(double Tkin) +{ + return sqrt(Tkin * (Tkin + 2.0 * T0)); +} + +double OneDimensionBpCpuModel::Kdiffr(double beta, double Rig) +{ + return K0 * beta * Rig; +} + +double OneDimensionBpCpuModel::Dp(double V, double p, double r) +{ + return (2.0 * V * p * dt / (3.0 * r)); +} + +double OneDimensionBpCpuModel::Dr(double V, double Kdiff, double r, double dt, double rand) +{ + return (V + (2.0 * Kdiff / r)) * dt + (rand * sqrt(2.0 * Kdiff * dt)); +} + +double OneDimensionBpCpuModel::W(double p) +{ + double w = (m0 * m0 * c * c * c * c) + (p * p * c * c); + return (pow(w, -1.85) / p) / 1e45; +} diff --git a/CpuImplementations/src/OneDimensionFpCpuModel.cpp b/CpuImplementations/src/OneDimensionFpCpuModel.cpp index f7deeec..efbe6b5 100644 --- a/CpuImplementations/src/OneDimensionFpCpuModel.cpp +++ b/CpuImplementations/src/OneDimensionFpCpuModel.cpp @@ -27,7 +27,7 @@ void OneDimensionFpCpuModel::runSimulation(ParamsCarrier *singleTone) unsigned int nthreads = std::thread::hardware_concurrency(); int targetIterations = ceil((double)singleTone->getInt("millions", 1) * 1000000.0 / ((double)nthreads * 101.0 * 250.0)); setContants(singleTone); - for (int iteration = 0; targetIterations < targetIterations; iteration++) + for (int iteration = 0; iteration < targetIterations; iteration++) { spdlog::info("Processed: {:03.2f}%", (float) iteration / ((float) targetIterations / 100.0)); std::vector threads; @@ -56,86 +56,116 @@ void OneDimensionFpCpuModel::simulation(int threadNumber, unsigned int available double Tkin, Tkininj, Rig, beta; double w; double Tkinw, p, rp, dp, pp, pinj, cfactor, sumac; + thread_local std::random_device rd{}; thread_local std::mt19937 generator(rd()); thread_local std::normal_distribution distribution(0.0f, 1.0f); - for (int energy = 0; energy < 101; energy++) - { - for (int particlePerEnergy = 0; particlePerEnergy < 250; particlePerEnergy++) - { + + auto& gen = generator; + auto& dist = distribution; + + std::vector localOutputs; + localOutputs.reserve(101); + + for (int energy = 0; energy < 101; energy++) { + for (int particlePerEnergy = 0; particlePerEnergy < 250; particlePerEnergy++) { Tkininj = getTkinInjection(((availableThreads * iteration + threadNumber) * 250) + particlePerEnergy, 0.0001, uniformEnergyInjectionMaximum, 10000); Tkin = Tkininj; Tkinw = Tkin * 1e9 * q; - Rig = sqrt(Tkinw * (Tkinw + (2.0 * T0w))) / q; + Rig = RigFromTkinJoule(Tkin); p = Rig * q / c; pinj = p; - // Equation under Equation 6 from - // Yamada et al. "A stochastic view of the solar modulation phenomena of cosmic rays" GEOPHYSICAL RESEARCH LETTERS, VOL. 25, NO.13, PAGES 2353-2356, JULY 1, 1998 - // https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/98GL51869 - w = (m0 * m0 * c * c * c * c) + (p * p * c * c); - w = pow(w, -1.85) / p; - w = w / 1e45; - + w = W(p); sumac = 0.0; r = 100.0001; - while (r < 100.0002) - { - // Equation 5 - beta = sqrt(Tkin * (Tkin + T0 + T0)) / (Tkin + T0); - - // Equation 6 - Rig = sqrt(Tkin * (Tkin + (2.0 * T0))); - - // Equation 7 - Kdiff = K0 * beta * Rig; - - arnum = distribution(generator); - - // Equation 9 - dr = (V + (2.0 * Kdiff / r)) * dt; - dr = dr + (arnum * sqrt(2.0 * Kdiff * dt)); - - // Equation 10 - dp = 2.0 * V * p * dt / (3.0 * r); + while (r < 100.0002) { + beta = Beta(Tkin); + Rig = RigFromTkin(Tkin); + Kdiff = Kdiffr(beta, Rig); + arnum = dist(gen); - // Equation 11 - cfactor = 4.0 * V / (3.0 * r); - sumac = sumac + (cfactor * dt); + dr = Dr(V, Kdiff, r, dt, arnum); + dp = Dp(V, p, r); + cfactor = Cfactor(V, r); + sumac += cfactor * dt; rp = r; pp = p; - r = r + dr; - p = p - dp; + r += dr; + p -= dp; - // Equation 6 in J Rig = p * c / q; - Tkin = sqrt((T0 * T0 * q * q * 1e9 * 1e9) + (q * q * Rig * Rig)) - (T0 * q * 1e9); - Tkin = Tkin / (q * 1e9); - - if (beta > 0.01 & Tkin < 100.0) - { - if ((r - 1.0) / ((r - dr) - 1.0) < 0.0) - { - outputMutex.lock(); - outputQueue.push({pinj, p, r, w, sumac}); - outputMutex.unlock(); + Tkin = TkinFromRig(Rig); + + if (beta > 0.01 && Tkin < 100.0) { + if ((r - 1.0) / ((r - dr) - 1.0) < 0.0) { + + localOutputs.emplace_back(SimulationOutput{pinj, p, r, w, sumac}); + } } + if (beta < 0.01) break; + if (r < 0.1) { r = rp; p = pp; } + } + } + } + + + std::lock_guard lock(outputMutex); + for (const auto& output : localOutputs) + { + outputQueue.push(output); + } + +} - if (beta < 0.01) - { - break; - } - if (r < 0.1) - { - r = rp; - p = pp; - } - } - } - } +double OneDimensionFpCpuModel::Beta(double Tkin) +{ + return sqrt(Tkin * (Tkin + 2.0 * T0)) / (Tkin + T0); +} + +double OneDimensionFpCpuModel::RigFromTkin(double Tkin) +{ + return sqrt(Tkin * (Tkin + 2.0 * T0)); +} + +double OneDimensionFpCpuModel::RigFromTkinJoule(double Tkin) +{ + double Tkinw = Tkin * 1e9 * q; + return sqrt(Tkinw * (Tkinw + 2.0 * T0w)) / q; +} + +double OneDimensionFpCpuModel::Kdiffr(double beta, double Rig) +{ + return K0 * beta * Rig; +} + +double OneDimensionFpCpuModel::Dp(double V, double p, double r) +{ + return 2.0 * V * p * dt / (3.0 * r); +} + +double OneDimensionFpCpuModel::Dr(double V, double Kdiff, double r, double dt, double rand) +{ + return (V + 2.0 * Kdiff / r) * dt + (rand * sqrt(2.0 * Kdiff * dt)); +} + +double OneDimensionFpCpuModel::Cfactor(double V, double r) +{ + return 4.0 * V / (3.0 * r); +} + +double OneDimensionFpCpuModel::TkinFromRig(double Rig) +{ + return (sqrt((T0 * T0 * q * q * 1e9 * 1e9) + (q * q * Rig * Rig)) - (T0 * q * 1e9)) / (q * 1e9); +} + +double OneDimensionFpCpuModel::W(double p) +{ + double w = (m0 * m0 * c * c * c * c) + (p * p * c * c); + return pow(w, -1.85) / p / 1e45; } \ No newline at end of file diff --git a/Factory/src/CosmicFactory.cpp b/Factory/src/CosmicFactory.cpp index adb5125..f15f709 100644 --- a/Factory/src/CosmicFactory.cpp +++ b/Factory/src/CosmicFactory.cpp @@ -18,6 +18,7 @@ AbstractAlgorithm *CosmicFactory::getAlgorithm(std::string name) else if (name.compare("2D SolarProp-like") == 0) { return new SolarPropLikeAlgorithm(); + } else if (name.compare("2D Geliosphere") == 0) {