Skip to content

Commit 157dde7

Browse files
MSABuschmannMSABuschmann
committed
Merge remote-tracking branch 'origin/development'
Merge developmental changes into the main branch.
2 parents e9c268b + 4180ec3 commit 157dde7

36 files changed

+2028
-1123
lines changed

‎Make.sledgehamr‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,14 @@ CEXE_headers += $(shell ls $(SLEDGEHAMR_PROJECT_PATH)/*/*.h | xargs -n 1 basenam
3434
include $(SLEDGEHAMR_HOME_ABS)/Make.projects
3535

3636
# SWFFT
37-
include $(AMREX_HOME)/Src/Extern/SWFFT/Make.package
38-
INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
39-
VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
37+
#include $(AMREX_HOME)/Src/Extern/SWFFT/Make.package
38+
#INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
39+
#VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
40+
41+
include $(AMREX_HOME)/Src/FFT/Make.package
42+
INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/FFT
43+
VPATH_LOCATIONS += $(AMREX_HOME)/Src/FFT
44+
4045
LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3_omp -lfftw3
4146

4247
# Boost
@@ -49,6 +54,7 @@ LIBRARIES += -lhdf5
4954
# AMReX
5055
USE_MPI = TRUE
5156
USE_OMP = TRUE
57+
USE_FFT = TRUE
5258
DIM = 3
5359

5460
include $(AMREX_HOME)/Tools/GNUMake/Make.defs

‎projects/AxionOnly/AxionOnly.cpp‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
#include "AxionOnly.h"
2+
3+
namespace AxionOnly {
4+
5+
void AxionOnly::Init() {
6+
ParseVariables();
7+
SetProjections();
8+
SetSpectrum();
9+
}
10+
11+
void AxionOnly::ParseVariables() {
12+
amrex::ParmParse pp_prj("project");
13+
pp_prj.get("n", n);
14+
pp_prj.get("eta_c", eta_c);
15+
pp_prj.get("eta_star", eta_star);
16+
pp_prj.get("N_QCD", N_QCD);
17+
}
18+
19+
void AxionOnly::SetProjections() {
20+
io_module->projections.emplace_back(dtheta_prime2, "dtheta_prime2");
21+
}
22+
23+
void AxionOnly::SetSpectrum() {
24+
io_module->spectra.emplace_back(theta_spectrum, "theta");
25+
io_module->spectra.emplace_back(dtheta_spectrum, "dtheta");
26+
}
27+
28+
void AxionOnly::SetParamsRhs(std::vector<double> &params, const double time,
29+
const int lev) {
30+
const double eta = time;
31+
const double ma_sq = std::pow(std::min(eta, eta_c) / eta_star, n) / N_QCD;
32+
params.push_back(ma_sq);
33+
}
34+
35+
}; // namespace AxionOnly

‎projects/AxionOnly/AxionOnly.h‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
#ifndef PROJECTS_AXION_ONLY_H_
2+
#define PROJECTS_AXION_ONLY_H_
3+
4+
#include <sledgehamr.h>
5+
6+
#include "kernels_energy_densities.h"
7+
#include "kernels_rhs.h"
8+
#include "kernels_tagging.h"
9+
10+
namespace AxionOnly {
11+
12+
SLEDGEHAMR_FINISH_SETUP
13+
14+
/** @brief Class to simulate axion strings.
15+
*/
16+
class AxionOnly : public sledgehamr::Sledgehamr {
17+
public:
18+
SLEDGEHAMR_INITIALIZE_PROJECT(AxionOnly)
19+
20+
/** @brief Override Init function and let cosmology module handle the setup.
21+
*/
22+
void Init() override;
23+
void SetParamsRhs(std::vector<double> &params, const double time,
24+
const int lev) override;
25+
26+
private:
27+
void ParseVariables();
28+
void SetProjections();
29+
void SetSpectrum();
30+
31+
double eta_c;
32+
double eta_star;
33+
double n;
34+
double N_QCD;
35+
};
36+
37+
}; // namespace AxionOnly
38+
39+
#endif // PROJECTS_AXION_ONLY_H_

‎projects/AxionOnly/kernels_energy_densities.h‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
#ifndef PROJECTS_AXION_ONLY_KERNELS_ENERGY_DENSITIES_H_
2+
#define PROJECTS_AXION_ONLY_KERNELS_ENERGY_DENSITIES_H_
3+
4+
#include "setup.h"
5+
6+
namespace AxionOnly {
7+
8+
/** @brief Kernel function computing axion energy density a'^2.
9+
* @param state Data from which to calculate RHS (current state).
10+
* @param i i-th cell index.
11+
* @param j j-th cell index.
12+
* @param k k-th cell index.
13+
* @param lev Current level.
14+
* @param time Current time.
15+
* @param dt Time step size.
16+
* @param dx Grid spacing.
17+
* @param params Optional parameters.
18+
* @return a'^2.
19+
*/
20+
AMREX_FORCE_INLINE
21+
double dtheta_prime2(amrex::Array4<amrex::Real const> const &state, const int i,
22+
const int j, const int k, const int lev, const double time,
23+
const double dt, const double dx,
24+
const std::vector<double> &params) {
25+
double dtheta = state(i, j, k, Scalar::dtheta);
26+
return dtheta * dtheta;
27+
}
28+
29+
/** @brief Kernel function computing axion energy density a'^2.
30+
* @param state Data from which to calculate RHS (current state).
31+
* @param i i-th cell index.
32+
* @param j j-th cell index.
33+
* @param k k-th cell index.
34+
* @param lev Current level.
35+
* @param time Current time.
36+
* @param dt Time step size.
37+
* @param dx Grid spacing.
38+
* @param params Optional parameters.
39+
* @return a'^2.
40+
*/
41+
AMREX_FORCE_INLINE
42+
double theta_spectrum(amrex::Array4<amrex::Real const> const &state,
43+
const int i, const int j, const int k, const int lev,
44+
const double time, const double dt, const double dx,
45+
const std::vector<double> &params) {
46+
return state(i, j, k, Scalar::theta);
47+
}
48+
49+
AMREX_FORCE_INLINE
50+
double dtheta_spectrum(amrex::Array4<amrex::Real const> const &state,
51+
const int i, const int j, const int k, const int lev,
52+
const double time, const double dt, const double dx,
53+
const std::vector<double> &params) {
54+
return state(i, j, k, Scalar::dtheta);
55+
}
56+
57+
}; // namespace AxionOnly
58+
59+
#endif // PROJECTS_AXION_ONLY_KERNELS_ENERGY_DENSITIES_H_

‎projects/AxionOnly/kernels_rhs.h‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
#ifndef PROJECTS_AXION_ONLY_KERNELS_RHS_H_
2+
#define PROJECTS_AXION_ONLY_KERNELS_RHS_H_
3+
4+
#include "setup.h"
5+
#include <sledgehamr_utils.h>
6+
7+
namespace AxionOnly {
8+
9+
/** @brief Function that calculates the RHS of the EOM at a single cell.
10+
* @param rhs Container to be filled with RHS.
11+
* @param state Data from which to calculate RHS (current state).
12+
* @param i i-th cell index.
13+
* @param j j-th cell index.
14+
* @param k k-th cell index.
15+
* @param lev Current level.
16+
* @param time Current time.
17+
* @param dt Time step size.
18+
* @param dx Grid spacing.
19+
* @param params Optional parameters.
20+
*/
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) {
25+
// Fetch field values.
26+
double theta = state(i, j, k, Scalar::theta);
27+
double dtheta = state(i, j, k, Scalar::dtheta);
28+
double eta = time;
29+
30+
// Compute Laplacians.
31+
constexpr int order = 2;
32+
double laplacian_theta = sledgehamr::utils::Laplacian<order>(
33+
state, i, j, k, Scalar::theta, dx * dx);
34+
35+
// Compute EOM.
36+
double ma_sq = params[0];
37+
double potential = ma_sq * eta * eta * std::sin(theta);
38+
39+
rhs(i, j, k, Scalar::theta) = dtheta;
40+
rhs(i, j, k, Scalar::dtheta) =
41+
-dtheta * 2. / eta + laplacian_theta - potential;
42+
}
43+
44+
}; // namespace AxionOnly
45+
46+
#endif // PROJECTS_AXION_ONLY_KERNELS_RHS_H_

‎projects/AxionOnly/kernels_tagging.h‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
#ifndef PROJECTS_AXION_ONLY_KERNELS_TAGGING_H_
2+
#define PROJECTS_AXION_ONLY_KERNELS_TAGGING_H_
3+
4+
#include "setup.h"
5+
6+
namespace AxionOnly {
7+
8+
/** @brief Modifies the truncation error criteria for Pi1 and Pi2 from its
9+
* default \tau > \tau_{crit} to \tau * \Delta t_{\ell} > \tau_{crit}.
10+
* @param truncation_error \tau
11+
* @return f(\tau) for criteria f(\tau) > \tau_{crit}.
12+
*/
13+
#define AXION_ONLY_TRUNCATION_MODIFIER(x) \
14+
template <> \
15+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE double TruncationModifier<x>( \
16+
const amrex::Array4<const double> &state, const int i, const int j, \
17+
const int k, const int lev, const double time, const double dt, \
18+
const double dx, const double truncation_error, \
19+
const double *params) { \
20+
return truncation_error * dt / time; \
21+
}
22+
23+
#define AXION_ONLY_TRUNCATION_MODIFIER2(x) \
24+
template <> \
25+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE double TruncationModifier<x>( \
26+
const amrex::Array4<const double> &state, const int i, const int j, \
27+
const int k, const int lev, const double time, const double dt, \
28+
const double dx, const double truncation_error, \
29+
const double *params) { \
30+
return truncation_error / time; \
31+
}
32+
33+
AXION_ONLY_TRUNCATION_MODIFIER(Scalar::dtheta)
34+
AXION_ONLY_TRUNCATION_MODIFIER2(Scalar::theta)
35+
36+
}; // namespace AxionOnly
37+
38+
#endif // PROJECTS_AXION_ONLY_KERNELS_TAGGING_H_

‎projects/AxionOnly/setup.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
#ifndef PROJECTS_AXION_ONLY_SETUP_H_
2+
#define PROJECTS_AXION_ONLY_SETUP_H_
3+
4+
namespace AxionOnly {
5+
6+
SLEDGEHAMR_ADD_SCALARS(theta)
7+
SLEDGEHAMR_ADD_CONJUGATE_MOMENTA(dtheta)
8+
9+
}; // namespace AxionOnly
10+
11+
#endif // PROJECTS_AXION_ONLY_SETUP_H_

‎projects/AxionStrings/cosmology.h‎

Lines changed: 14 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,9 @@
33

44
#include <sledgehamr.h>
55

6-
#include "setup.h"
7-
#include "kernels_tagging.h"
86
#include "kernels_energy_densities.h"
7+
#include "kernels_tagging.h"
8+
#include "setup.h"
99

1010
namespace AxionStrings {
1111

@@ -14,7 +14,7 @@ namespace AxionStrings {
1414
*/
1515
class Cosmology {
1616
public:
17-
void Init(sledgehamr::Sledgehamr* owner);
17+
void Init(sledgehamr::Sledgehamr *owner);
1818

1919
/** @brief We only want to create a new level if string width is below
2020
* threshold.
@@ -23,41 +23,37 @@ class Cosmology {
2323
* @return Whether we want to create level lev.
2424
*/
2525
bool CreateLevelIf(const int lev, const double time) {
26-
return StringWidth(lev-1, time) <= string_width_threshold;
26+
return StringWidth(lev - 1, time) <= string_width_threshold;
2727
}
2828

2929
/** @brief Current radial mode mass.
3030
* @param eta Current time eta.
3131
* @return m_r(eta)
3232
*/
33-
double Mr(const double eta) {
34-
return std::sqrt(2. * lambda) * eta;
35-
}
33+
double Mr(const double eta) { return std::sqrt(2. * lambda) * eta; }
3634

3735
/** @brief Current hubble time.
3836
* @param eta Current time eta.
3937
* @return H(eta).
4038
*/
41-
double H(const double eta) {
42-
return 1./eta;
43-
}
39+
static double H(const double eta) { return 1. / eta; }
4440

4541
/** @brief Current string width in units of the grid spacing.
4642
* @param lev Current level.
4743
* @param eta Current time eta.
4844
* @return String width.
4945
*/
5046
double StringWidth(const int lev, const double eta) {
51-
return 1./(Mr(eta) * sim->GetDx(lev));
47+
return 1. / (Mr(eta) * sim->GetDx(lev));
5248
}
5349

5450
/** @brief Returns time at which a level will be introduced.
5551
* @param lev Level.
5652
* @return Refinement time for level lev.
5753
*/
5854
double RefinementTime(const int lev) {
59-
return sim->GetDimN(lev) / (sqrt(2.*lambda)
60-
* string_width_threshold * sim->GetL());
55+
return sim->GetDimN(lev) /
56+
(sqrt(2. * lambda) * string_width_threshold * sim->GetL());
6157
}
6258

6359
/** @brief Computes scale separation.
@@ -67,7 +63,7 @@ class Cosmology {
6763
double Log(const double eta) {
6864
if (eta <= 0)
6965
return -DBL_MAX;
70-
return std::log( Mr(eta) / H(eta) );
66+
return std::log(Mr(eta) / H(eta));
7167
}
7268

7369
/** @brief Computes the physical box length.
@@ -90,17 +86,15 @@ class Cosmology {
9086
* @return Physical Hubble parameter.
9187
*/
9288
double Hubble(const double T, const double mpl, const double gStar) {
93-
return std::sqrt(4.*std::pow(M_PI, 3) / 45. * gStar * std::pow(T, 4)
94-
/ std::pow(mpl, 2));
89+
return std::sqrt(4. * std::pow(M_PI, 3) / 45. * gStar * std::pow(T, 4) /
90+
std::pow(mpl, 2));
9591
}
9692

9793
double XiTime(double T, double mpl, double gStar) {
9894
return 0.3012 / std::sqrt(gStar) * mpl / std::pow(T, 2);
9995
}
10096

101-
double XiTemp(double eta, double T1) {
102-
return T1 / eta;
103-
}
97+
double XiTemp(double eta, double T1) { return T1 / eta; }
10498

10599
double Xi(const int lev, const double eta);
106100

@@ -132,7 +126,7 @@ class Cosmology {
132126

133127
/** @brief Pointer to the simulation.
134128
*/
135-
sledgehamr::Sledgehamr* sim;
129+
sledgehamr::Sledgehamr *sim;
136130
};
137131

138132
}; // namespace AxionStrings

0 commit comments

Comments
 (0)