Skip to content

Commit

Permalink
Updated Geliosphere input table, updated equation changes in Geliosph… (
Browse files Browse the repository at this point in the history 
#9)

* Updated implementations of Geliosphere models according to updated equations

* Updated evaluation of 2D models

* Added evaluation for 2D models matching AMS bins

* Added logging for tilt angle and polarity for 2D models
  • Loading branch information
@msolanik
msolanik authored Jun 11, 2023
1 parent 5c25191 commit 8a6820f
Show file tree
Hide file tree
Showing 8 changed files with 595 additions and 543 deletions.
16 changes: 16 additions & 0 deletions Algorithm/include/ResultConstants.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,4 +100,20 @@ const double SPbins[32] = { 0, 0.01, 0.015, 0.0225, 0.03375, 0.050625,
*/
const double UlyssesBins[4] = { 0, 0.125, 0.250, 2.0 };

/**
* @brief Array of bins for 2D models.
*
*/
const double AmsBins[72] = { 4.924000e-01, 6.207000e-01, 7.637000e-01, 9.252000e-01, 1.105000e+00,
1.303000e+00, 1.523000e+00, 1.765000e+00, 2.034000e+00, 2.329000e+00, 2.652000e+00, 3.005000e+00,
3.390000e+00, 3.810000e+00, 4.272000e+00, 4.774000e+00, 5.317000e+00, 5.906000e+00, 6.546000e+00,
7.236000e+00, 7.981000e+00, 8.787000e+00, 9.653000e+00, 1.060000e+01, 1.160000e+01, 1.264000e+01,
1.379000e+01, 1.504000e+01, 1.639000e+01, 1.784000e+01, 1.938000e+01, 2.103000e+01, 2.283000e+01,
2.478000e+01, 2.683000e+01, 2.903000e+01, 3.138000e+01, 3.387000e+01, 3.657000e+01, 3.947000e+01,
4.257000e+01, 4.587000e+01, 4.942000e+01, 5.322000e+01, 5.727000e+01, 6.162000e+01, 6.632000e+01,
7.137000e+01, 7.677000e+01, 8.257000e+01, 8.882000e+01, 9.557000e+01, 1.031000e+02, 1.111000e+02,
1.196000e+02, 1.291000e+02, 1.401000e+02, 1.526000e+02, 1.666000e+02, 1.826000e+02, 2.006000e+02,
2.211000e+02, 2.451000e+02, 2.741000e+02, 3.096000e+02, 3.536000e+02, 4.091000e+02, 4.821000e+02,
5.831000e+02, 7.316000e+02, 9.751000e+02, 1.464000e+03 };

#endif
15 changes: 14 additions & 1 deletion Algorithm/src/TwoDimensionBpResults.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ void TwoDimensionBpResults::runAlgorithm(ParamsCarrier *singleTone)
ResultsUtils *resultsUtils = new ResultsUtils();
double w, Rig, p1AU, Tkin, r, p, Tkinw, Rig1AU, Tkininj, theta, thetainj, tt, t2, beta;
double tem6, tem5, jlis, tem, wJGR, jlisJGR, jlisBurger, wBurger;
double speSP[31] = {0}, speJGR[31] = {0}, speN[31] = {0}, speBurger[31] = {0}, ulyssesBin[4] = {0}, ulyssesBinN[4] = {0};
double speSP[31] = {0}, speJGR[31] = {0}, speN[31] = {0}, speBurger[31] = {0}, ulyssesBin[4] = {0}, ulyssesBinN[4] = {0}, amsBin[72] = {0}, amsBinN[72] = {0};
FILE *inputFile = fopen("log.dat", "r");
int numberOfIterations = resultsUtils->countLines(inputFile) - 1;
int targetArray[] = {numberOfIterations};
Expand Down Expand Up @@ -94,6 +94,14 @@ void TwoDimensionBpResults::runAlgorithm(ParamsCarrier *singleTone)
ulyssesBinN[ii]++;
}
}
for (int ii = 1; ii < 71; ii++)
{
if ((Tkininj > AmsBins[ii]) && (Tkininj < AmsBins[ii+1]))
{
amsBin[ii] = amsBin[ii] + wBurger;
amsBinN[ii]++;
}
}
}
fclose(inputFile);

Expand All @@ -117,4 +125,9 @@ void TwoDimensionBpResults::runAlgorithm(ParamsCarrier *singleTone)
spectrumOutput.isCsv = singleTone->getInt("csv", 0);
resultsUtils->writeSpectrum(&spectrumOutput, ulyssesBinN, ulyssesBin, SPECTRUM_ULYSSES);
spdlog::info("Spectrum with Ulysses bins based on Burger has been written to file.");

spectrumOutput.fileName = "Ams";
spectrumOutput.isCsv = singleTone->getInt("csv", 0);
resultsUtils->writeSpectrum(&spectrumOutput, amsBinN, amsBin, SPECTRUM_AMS);
spdlog::info("Spectrum with AMS bins based on Burger has been written to file.");
}
19 changes: 9 additions & 10 deletions CUDAKernel/src/GeliosphereModel.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,19 +159,18 @@ __global__ void trajectorySimulationGeliosphere(trajectoryHistoryGeliosphere *hi
COmega = COmega + (2.0f * r * deltarh2);

// Equation 36
dKper = ratio * K0 * beta * Rig * ((2.0f * r * sqrtf(Cb)) - (r * r * COmega / (2.0f * sqrtf(Cb)))) / Cb;
dKper = ratio * K0 * beta * Rig * ((2.0f * r * sqrtf(Cb)) - (r * r * COmega / (2.0f * sqrtf(Cb)))) / (3.0f * (5.0f / 3.4f) * Cb);

// Equation 35
dKrr = dKper + ((1.0f - ratio) * K0 * beta * Rig * ((2.0f * r * powf(Cb, 1.5f)) - (r * r * COmega * 3.0f * sqrtf(Cb) / 2.0f)) / powf(Cb, 3.0f));
dKrr = dKrr * 5.0f / (3.0f * 3.4f);
dKrr = dKper + ((1.0f - ratio) * K0 * beta * Rig * ((2.0f * r * powf(Cb, 1.5f)) - (r * r * COmega * 3.0f * sqrtf(Cb) / 2.0f)) / ( 3.0f * (5.0f / 3.4f) * powf(Cb, 3.0f)));

if ((theta > (1.7f * Pi / 180.0f)) && (theta < (178.3f * Pi / 180.0f)))
{
// Equation 37
CKtt = sinf(theta) * cosf(theta) * (omega * omega * r * r / (V * V));

// Equation 38
dKtt1 = (-1.0f * ratio * K0 * beta * Rig * r * r * CKtt) / powf(Cb, 1.5f);
dKtt1 = (-1.0f * ratio * K0 * beta * Rig * r * r * CKtt) / (3.0f * (5.0f / 3.4f) * powf(Cb, 1.5f));

// Equation 39
dKtt2 = (1.0f - ratio) * K0 * beta * Rig * r * r * r * r * deltarh2;
Expand All @@ -192,34 +191,34 @@ __global__ void trajectorySimulationGeliosphere(trajectoryHistoryGeliosphere *hi
CKtt = CKtt - (r * r * delta0 * delta0 * cos(theta) / (rh * rh * sin3));

// Equation 38
dKtt1 = (-1.0f * ratio * K0 * beta * Rig * r * r * CKtt) / powf(Cb, 1.5f);
dKtt1 = (-1.0f * ratio * K0 * beta * Rig * r * r * CKtt) / (3.0f * (5.0f / 3.4f) * powf(Cb, 1.5f));

// Equation 39
dKtt2 = (1.0f - ratio) * K0 * beta * Rig * r * r * r * r * delta0 * delta0 / (rh * rh);

// Equation 40
dKtt3 = -2.0f * (cosf(theta) / sin3) / powf(Cb, 1.5f);
dKtt3 = -2.0f * (cosf(theta) / sin3) / (3.0f * (5.0f / 3.4f) * powf(Cb, 1.5f));

// Equation 41
dKtt4 = 3.0f * CKtt / (sin2 * powf(Cb, 2.5f));
dKtt4 = 3.0f * (CKtt / sin2) / (3.0f * (5.0f / 3.4f) * powf(Cb, 2.5f));

// Equation 42
dKtt = dKtt1 + (dKtt2 * (dKtt3 - dKtt4));
}

// Equation 43
dKrtr = (1.0f - ratio) * K0 * beta * Rig * deltarh * 3.0 * r * r / powf(Cb, 2.5f);
dKrtr = (1.0f - ratio) * K0 * beta * Rig * deltarh * r * r / (3.0f * (5.0f / 3.4f) * powf(Cb, 2.5f));

// Equation 44
if ((theta > (1.7f * Pi / 180.0f)) && (theta < (178.3f * Pi / 180.0f)))
{
dKrtt = (1.0f - ratio) * K0 * beta * Rig * r * r * r / (rh * powf(Cb, 2.5f));
dKrtt = (1.0f - ratio) * K0 * beta * Rig * r * r * r / (3.0f * (5.0f / 3.4f) * rh * powf(Cb, 2.5f));
dKrtt = -1.0f * dKrtt * delta;
dKrtt = dKrtt * 3.0f * gamma2 * cosf(theta) / sinf(theta);
}
else
{
dKrtt = (1.0f - ratio) * K0 * beta * Rig * r * r * r / (rh * powf(Cb, 2.5f));
dKrtt = (1.0f - ratio) * K0 * beta * Rig * r * r * r / (3.0f * (5.0f / 3.4f) * rh * powf(Cb, 2.5f));
dKrtt = -1.0f * dKrtt * delta0 * cosf(theta) / (sinf(theta) * sinf(theta));
dKrtt = dKrtt * (1.0f - (2.0f * r * r * deltarh2) + (4.0f * gamma2));
}
Expand Down
19 changes: 9 additions & 10 deletions CpuImplementations/src/GeliosphereCpuModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,19 +175,18 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr
COmega = COmega + (2.0 * r * deltarh2);

// Equation 36
dKper = ratio * K0 * beta * Rig * ((2.0 * r * sqrt(Cb)) - (r2 * COmega / (2.0 * sqrt(Cb)))) / Cb;
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)) / pow(Cb, 3.0));
dKrr = dKrr * 5. / (3. * 3.4);
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) / pow(Cb, 1.5);
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;
Expand All @@ -208,34 +207,34 @@ void GeliosphereCpuModel::simulation(int threadNumber, unsigned int availableThr
CKtt = CKtt - (r2 * delta0 * delta0 * cos(theta) / (rh * rh * sin3));

// Equation 38
dKtt1 = (-1.0 * ratio * K0 * beta * Rig * r2 * CKtt) / pow(Cb, 1.5);
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) / pow(Cb, 1.5);
dKtt3 = -2.0 * (cos(theta) / sin3) / (3.0 * (5.0 / 3.4) * pow(Cb, 1.5));

// Equation 41
dKtt4 = 3.0 * CKtt / (sin2 * pow(Cb, 2.5));
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 * 3.0 * r2 / pow(Cb, 2.5);
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 / (rh * pow(Cb, 2.5));
dKrtt = (1.0 - ratio) * K0 * beta * Rig * r2 * r / ((3.0 * (5.0 / 3.4) * rh * pow(Cb, 2.5)));
dKrtt = -1.0 * dKrtt * delta;
dKrtt = dKrtt * 3.0 * gamma2 * cos(theta) / sin(theta);
}
else
{
dKrtt = (1.0 - ratio) * K0 * beta * Rig * r2 * r / (rh * pow(Cb, 2.5));
dKrtt = (1.0 - ratio) * K0 * beta * Rig * r2 * r / ((3.0 * (5.0 / 3.4) * rh * pow(Cb, 2.5)));
dKrtt = -1.0 * dKrtt * delta0 * cos(theta) / (sin(theta) * sin(theta));
dKrtt = dKrtt * (1.0 - (2.0 * r2 * deltarh2) + (4.0 * gamma2));
}
Expand Down
Loading

0 comments on commit 8a6820f

Please sign in to comment.