Skip to content

Commit 6fc126f

Browse files
MSABuschmannMSABuschmann
committed
Ability to preevolve axion string sim from truly random initial state
1 parent fd85a25 commit 6fc126f

File tree

3 files changed

+61
-28
lines changed

3 files changed

+61
-28
lines changed

‎projects/AxionStringsPreevolution/AxionStringsPreevolution.cpp‎

Lines changed: 34 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,17 @@
22

33
namespace AxionStringsPreevolution {
44

5-
void AxionStringsPreevolution::SetParamsRhs(
6-
std::vector<double>& params, const double time, const int lev) {
5+
void AxionStringsPreevolution::SetParamsRhs(std::vector<double> &params,
6+
const double time, const int lev) {
77
params.push_back(eta_0);
88
}
99

1010
void AxionStringsPreevolution::Init() {
1111
cosmo.Init(this);
1212
ParseConstants();
13+
14+
if (random_state > 0)
15+
SetRandomState();
1316
}
1417

1518
bool AxionStringsPreevolution::StopRunning(const double time) {
@@ -19,12 +22,40 @@ bool AxionStringsPreevolution::StopRunning(const double time) {
1922
return (xi <= xi_0 && time >= min_eta);
2023
}
2124

25+
void AxionStringsPreevolution::SetRandomState() {
26+
sledgehamr::LevelData &state = GetLevelData(0);
27+
amrex::InitRandom(random_state + amrex::ParallelDescriptor::MyProc());
28+
amrex::Print() << "Set Random Initial Conditions" << std::endl;
29+
30+
// #pragma omp parallel
31+
for (amrex::MFIter mfi(state, true); mfi.isValid(); ++mfi) {
32+
const amrex::Array4<double> &state_arr = state.array(mfi);
33+
const amrex::Box &bx = mfi.tilebox();
34+
35+
const amrex::Dim3 lo = amrex::lbound(bx);
36+
const amrex::Dim3 hi = amrex::ubound(bx);
37+
38+
for (int k = lo.z; k <= hi.z; ++k) {
39+
for (int j = lo.y; j <= hi.y; ++j) {
40+
AMREX_PRAGMA_SIMD
41+
for (int i = lo.x; i <= hi.x; ++i) {
42+
state_arr(i, j, k, Scalar::Psi1) =
43+
amrex::Random() * 2. - 1.;
44+
state_arr(i, j, k, Scalar::Psi2) =
45+
amrex::Random() * 2. - 1.;
46+
}
47+
}
48+
}
49+
}
50+
}
51+
2252
void AxionStringsPreevolution::ParseConstants() {
2353
amrex::ParmParse pp("project");
54+
pp.get("random", random_state);
2455
pp.get("starting_log", log_0);
2556
pp.get("starting_xi", xi_0);
2657
pp.get("min_eta", min_eta);
27-
eta_0 = std::sqrt(std::exp(log_0)/std::sqrt(2.));
58+
eta_0 = std::sqrt(std::exp(log_0) / std::sqrt(2.));
2859
}
2960

3061
} // namespace AxionStringsPreevolution

‎projects/AxionStringsPreevolution/AxionStringsPreevolution.h‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
#ifndef PROJECTS_AXION_STRINGS_PREEVOLUTION_H_
22
#define PROJECTS_AXION_STRINGS_PREEVOLUTION_H_
33

4-
#include <sledgehamr.h>
54
#include "../AxionStrings/cosmology.h"
65
#include "kernels_rhs.h"
6+
#include <sledgehamr.h>
77

88
namespace AxionStringsPreevolution {
99

@@ -15,20 +15,23 @@ class AxionStringsPreevolution : public sledgehamr::Sledgehamr {
1515
public:
1616
SLEDGEHAMR_INITIALIZE_PROJECT(AxionStringsPreevolution)
1717

18-
void SetParamsRhs(std::vector<double>& params, const double time,
18+
void SetParamsRhs(std::vector<double> &params, const double time,
1919
const int lev) override;
2020

2121
void Init() override;
2222

2323
bool StopRunning(const double time) override;
2424

25+
void SetRandomState();
26+
2527
private:
2628
void ParseConstants();
2729

2830
double log_0 = 2;
2931
double eta_0 = 2.3;
3032
double xi_0 = 0.18;
3133
double min_eta = 2;
34+
int random_state = -1;
3235
Cosmology cosmo;
3336
};
3437

‎projects/AxionStringsPreevolution/kernels_rhs.h‎

Lines changed: 22 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
#ifndef PROJECTS_AXION_STRINGS_PREEVOLUTION_KERNELS_RHS_H_
22
#define PROJECTS_AXION_STRINGS_PREEVOLUTION_KERNELS_RHS_H_
33

4-
#include <sledgehamr_utils.h>
54
#include "../AxionStrings/cosmology.h"
5+
#include <sledgehamr_utils.h>
66

77
namespace AxionStringsPreevolution {
88
using namespace AxionStrings;
@@ -18,41 +18,40 @@ using namespace AxionStrings;
1818
* @param dt Time step size.
1919
* @param dx Grid spacing.
2020
*/
21-
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
22-
void Rhs(const amrex::Array4<double>& rhs,
23-
const amrex::Array4<const double>& state,
24-
const int i, const int j, const int k, const int lev,
25-
const double time, const double dt, const double dx,
26-
const double* params) {
21+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
22+
Rhs(const amrex::Array4<double> &rhs, const amrex::Array4<const double> &state,
23+
const int i, const int j, const int k, const int lev, const double time,
24+
const double dt, const double dx, const double *params) {
2725
// Fetch field values.
2826
double Psi1 = state(i, j, k, Scalar::Psi1);
2927
double Psi2 = state(i, j, k, Scalar::Psi2);
30-
double Pi1 = state(i, j, k, Scalar::Pi1);
31-
double Pi2 = state(i, j, k, Scalar::Pi2);
28+
double Pi1 = state(i, j, k, Scalar::Pi1);
29+
double Pi2 = state(i, j, k, Scalar::Pi2);
3230

3331
double eta = time;
3432

3533
// Compute Laplacians.
3634
constexpr int order = 2;
3735
double laplacian_Psi1 = sledgehamr::utils::Laplacian<order>(
38-
state, i, j, k, Scalar::Psi1, dx*dx);
36+
state, i, j, k, Scalar::Psi1, dx * dx);
3937
double laplacian_Psi2 = sledgehamr::utils::Laplacian<order>(
40-
state, i, j, k, Scalar::Psi2, dx*dx);
38+
state, i, j, k, Scalar::Psi2, dx * dx);
4139

4240
// Compute EOM.
4341
const double eta_0 = params[0];
44-
const double eta_sq = eta*eta;
45-
const double lambda = eta_0*eta_0;
46-
const double drag = eta_0;
47-
48-
double potential = lambda/eta_sq*( Psi1*Psi1 + Psi2*Psi2 - 1. );
49-
50-
rhs(i, j, k, Scalar::Psi1) = Pi1;
51-
rhs(i, j, k, Scalar::Psi2) = Pi2;
52-
rhs(i, j, k, Scalar::Pi1) = - Pi1*3./eta + laplacian_Psi1/eta_sq/drag
53-
- Psi1*potential;
54-
rhs(i, j, k, Scalar::Pi2) = - Pi2*3./eta + laplacian_Psi2/eta_sq/drag
55-
- Psi2*potential;
42+
const double eta_sq = eta * eta;
43+
const double lambda = eta_0 * eta_0;
44+
const double drag = std::sqrt(eta_0);
45+
// const double drag = eta_0;
46+
47+
double potential = lambda / eta_sq * (Psi1 * Psi1 + Psi2 * Psi2 - 1.);
48+
49+
rhs(i, j, k, Scalar::Psi1) = Pi1;
50+
rhs(i, j, k, Scalar::Psi2) = Pi2;
51+
rhs(i, j, k, Scalar::Pi1) =
52+
-Pi1 * 3. / eta + laplacian_Psi1 / eta_sq / drag - Psi1 * potential;
53+
rhs(i, j, k, Scalar::Pi2) =
54+
-Pi2 * 3. / eta + laplacian_Psi2 / eta_sq / drag - Psi2 * potential;
5655
}
5756

5857
}; // namespace AxionStringsPreevolution

0 commit comments

Comments
 (0)