Skip to content

Commit

Permalink
clang
Browse files Browse the repository at this point in the history
  • Loading branch information
@llibert94
llibert94 committed Jul 29, 2023
1 parent da1774f commit bbc08dc
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 74 deletions.
44 changes: 22 additions & 22 deletions Examples/PerfectFluid/Fluxes.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,50 +20,50 @@ vars_t<data_t> compute_flux(const vars_t<data_t> &vars, const int idir)
vars_t<data_t> out;
out.D = (vars.lapse * vars.vi[idir] - vars.shift[idir]) * vars.D;
Tensor<1, data_t> Sj_U;
FOR(i) {
Sj_U[i] = 0;
FOR(j) Sj_U[i] += vars.chi * h_UU[i][j] * vars.Sj[j];
FOR(i)
{
Sj_U[i] = 0;
FOR(j) Sj_U[i] += vars.chi * h_UU[i][j] * vars.Sj[j];
}

data_t chi_regularised = simd_max(1e-6, vars.chi);
Tensor<1, data_t> vi_D;
FOR(i)
{
vi_D[i] = 0.;
FOR(j)
{
vi_D[i] += vars.h[i][j] * vars.vi[j] / chi_regularised;
}
FOR(j) { vi_D[i] += vars.h[i][j] * vars.vi[j] / chi_regularised; }
}

data_t v2 = 0.;
FOR(i) v2 += vars.vi[i] * vi_D[i];

data_t P = vars.rho * (1. + vars.eps) / 3.; //for now we assume a conformal fluid

data_t P =
vars.rho * (1. + vars.eps) / 3.; // for now we assume a conformal fluid
data_t WW = 1. / (1. - v2);
data_t hh = 1. + vars.eps + P / vars.rho;
FOR(j)

FOR(j)
{
out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] -
vars.shift[idir] * vars.Sj[j];
FOR(k) out.Sj[j] += vars.lapse * P * h_UU[idir][k] * vars.h[j][k];
out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] -
vars.shift[idir] * vars.Sj[j];
FOR(k) out.Sj[j] += vars.lapse * P * h_UU[idir][k] * vars.h[j][k];
}

out.tau = vars.lapse * (Sj_U[idir] - vars.D * vars.vi[idir]) -
vars.shift[idir] * vars.tau;
vars.shift[idir] * vars.tau;

return out;
}
template <class data_t, template <typename> class vars_t>
vars_t<data_t> compute_num_flux(const vars_t<data_t> &vars, const int idir, const double lambda) {
vars_t<data_t> compute_num_flux(const vars_t<data_t> &vars, const int idir,
const double lambda)
{
vars_t<data_t> out = compute_flux(vars, idir);
out.D -= lambda * vars.D;
FOR(j) out.Sj[j] -= lambda * vars.Sj[j];
out.tau -= lambda * vars.tau;
return out;

}

}
} // namespace Fluxes
#endif /* FLUXES_HPP_ */
90 changes: 46 additions & 44 deletions Examples/PerfectFluid/InitialFluidData.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,39 +11,38 @@
#include "Coordinates.hpp"
#include "FluidCCZ4RHS.hpp"
#include "PerfectFluid.hpp"
#include "PrimitiveRecovery.hpp"
#include "Tensor.hpp"
#include "UserVariables.hpp" //This files needs NUM_VARS - total no. components
#include "VarsTools.hpp"
#include "simd.hpp"
#include "PrimitiveRecovery.hpp"

//! Class which sets the initial fluid matter config
// aka Kevin data
class InitialFluidData
{
// Now the non grid ADM vars
template <class data_t>
using MetricVars = ADMConformalVars::VarsWithGauge<data_t>;
// Now the non grid ADM vars
template <class data_t>
using MetricVars = ADMConformalVars::VarsWithGauge<data_t>;

public:
//! A structure for the input params for fluid properties and initial
//! conditions
struct params_t
{
std::array<double, CH_SPACEDIM>
center; //!< Centre of perturbation in initial SF bubble
center; //!< Centre of perturbation in initial SF bubble
double rho0;
double uflow;
double amplitude;
double awidth;
double sigma;
std::array<double, 2>
ycenter;
std::array<double, 2> ycenter;
};

//! The constructor
InitialFluidData(params_t a_params, double a_dx)
: m_dx(a_dx), m_params(a_params)
InitialFluidData(params_t a_params, double a_dx)
: m_dx(a_dx), m_params(a_params)
{
}

Expand All @@ -53,47 +52,50 @@ class InitialFluidData
// where am i?
Coordinates<data_t> coords(current_cell, m_dx, m_params.center);

MetricVars<data_t> metric_vars;

data_t x = coords.x;
double y = coords.y;
double z = coords.z;

Tensor<1, data_t> vi, Sj;
data_t chi_regularised = simd_max(metric_vars.chi, 1e-6);

vi[0] = m_params.uflow * (tanh((y - m_params.ycenter[0]) / m_params.awidth)
- tanh((y - m_params.ycenter[1]) / m_params.awidth) - 1.);
vi[1] = m_params.amplitude * sin(2. * M_PI * x)
* (exp(-pow((y - m_params.ycenter[0]) / m_params.sigma,2))
+ exp(-pow((y - m_params.ycenter[1]) / m_params.sigma,2)));
vi[2] = 0.;

data_t nn = 1. + 0.5 * (tanh((y - m_params.ycenter[0]) / m_params.awidth)
- tanh((y - m_params.ycenter[1]) / m_params.awidth));
data_t eps = 0.;
MetricVars<data_t> metric_vars;

data_t x = coords.x;
double y = coords.y;
double z = coords.z;

Tensor<1, data_t> vi, Sj;
data_t chi_regularised = simd_max(metric_vars.chi, 1e-6);

vi[0] = m_params.uflow *
(tanh((y - m_params.ycenter[0]) / m_params.awidth) -
tanh((y - m_params.ycenter[1]) / m_params.awidth) - 1.);
vi[1] = m_params.amplitude * sin(2. * M_PI * x) *
(exp(-pow((y - m_params.ycenter[0]) / m_params.sigma, 2)) +
exp(-pow((y - m_params.ycenter[1]) / m_params.sigma, 2)));
vi[2] = 0.;

data_t nn =
1. + 0.5 * (tanh((y - m_params.ycenter[0]) / m_params.awidth) -
tanh((y - m_params.ycenter[1]) / m_params.awidth));
data_t eps = 0.;

data_t rho = m_params.rho0;
data_t v2 = 0.;
FOR(i,j) v2 += metric_vars.h[i][j] * vi[i] * vi[j] / chi_regularised;
data_t P = rho * (1. + eps) / 3.;
data_t WW = 1. / (1. - v2);
data_t v2 = 0.;
FOR(i, j) v2 += metric_vars.h[i][j] * vi[i] * vi[j] / chi_regularised;
data_t P = rho * (1. + eps) / 3.;
data_t WW = 1. / (1. - v2);
data_t hh = 1. + eps + P / rho;
data_t D = rho * sqrt(WW);
data_t tau = rho * hh * WW - P - D;
FOR(i) Sj[i] = rho * hh * WW * vi[i];

data_t D = rho * sqrt(WW);
data_t tau = rho * hh * WW - P - D;
FOR(i) Sj[i] = rho * hh * WW * vi[i];

// store the vars
current_cell.store_vars(rho, c_rho);
current_cell.store_vars(vi, GRInterval<c_vi1, c_vi3>());
current_cell.store_vars(rho, c_rho);
current_cell.store_vars(vi, GRInterval<c_vi1, c_vi3>());
current_cell.store_vars(D, c_D);
current_cell.store_vars(Sj, GRInterval<c_Sj1, c_Sj3>());
current_cell.store_vars(tau, c_tau);
current_cell.store_vars(Sj, GRInterval<c_Sj1, c_Sj3>());
current_cell.store_vars(tau, c_tau);
}
protected:
double m_dx;
const params_t m_params;

protected:
double m_dx;
const params_t m_params;
};

#endif /* INITIALFLUIDDATA_HPP_ */
20 changes: 12 additions & 8 deletions Examples/PerfectFluid/PerfectFluid.impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,18 +83,22 @@ void PerfectFluid<eos_t>::add_matter_rhs(
data_t advec_chi = 0.;
FOR(i) advec_chi += vars.shift[i] * d1.chi[i] / chi_regularised;
// source - vars * \partial_t\sqrt{\gamma}/\sqrt{\gamma}
rhs.D = source.D - vars.D * (vars.lapse * vars.K
- divshift + GR_SPACEDIM / 2. * advec_chi);
FOR(i) rhs.Sj[i] = source.Sj[i] - vars.Sj[i] * (vars.lapse *
vars.K - divshift + GR_SPACEDIM / 2. * advec_chi);
rhs.tau = source.tau - vars.tau * (vars.lapse * vars.K
- divshift + GR_SPACEDIM / 2. * advec_chi);
rhs.D = source.D - vars.D * (vars.lapse * vars.K - divshift +
GR_SPACEDIM / 2. * advec_chi);
FOR(i)
rhs.Sj[i] = source.Sj[i] - vars.Sj[i] * (vars.lapse * vars.K - divshift +
GR_SPACEDIM / 2. * advec_chi);
rhs.tau = source.tau - vars.tau * (vars.lapse * vars.K - divshift +
GR_SPACEDIM / 2. * advec_chi);

// - F^i\partial_i\sqrt{\gamma}/\sqrt{\gamma}
FOR(i) {
FOR(i)
{
vars_t<data_t> flux = Fluxes::compute_flux(vars, i);
rhs.D -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.D;
FOR(j) rhs.Sj[j] -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.Sj[j];
FOR(j)
rhs.Sj[j] -=
GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.Sj[j];
rhs.tau -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.tau;
}

Expand Down

0 comments on commit bbc08dc

Please sign in to comment.