diff --git a/Examples/PerfectFluid/.#PerfectFluid.impl.hpp b/Examples/PerfectFluid/.#PerfectFluid.impl.hpp new file mode 120000 index 000000000..3a402fc9b --- /dev/null +++ b/Examples/PerfectFluid/.#PerfectFluid.impl.hpp @@ -0,0 +1 @@ +dtraykova@sakura01.11661:1659939318 \ No newline at end of file diff --git a/Examples/PerfectFluid/DefaultEoS.hpp b/Examples/PerfectFluid/DefaultEoS.hpp new file mode 100644 index 000000000..9c71429de --- /dev/null +++ b/Examples/PerfectFluid/DefaultEoS.hpp @@ -0,0 +1,31 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef DEFAULTEOS_HPP_ +#define DEFAULTEOS_HPP_ + +#include "Tensor.hpp" +#include "simd.hpp" + +class DefaultEoS +{ + public: + //! The constructor + DefaultEoS() {} + + //! Set the potential function for the scalar field here to zero + template class vars_t> + void compute_eos(data_t &P_of_rho0, data_t &dPdrho0, + const vars_t &vars) const + { + // The pressure value at rho0 + P_of_rho0 = 0.0; + + // The pressure gradient at rho0 + dPdrho0 = 0.0; + } +}; + +#endif /* DEFAULTEOS_HPP_ */ diff --git a/Examples/PerfectFluid/DiagnosticVariables.hpp b/Examples/PerfectFluid/DiagnosticVariables.hpp new file mode 100644 index 000000000..cafda91b2 --- /dev/null +++ b/Examples/PerfectFluid/DiagnosticVariables.hpp @@ -0,0 +1,28 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef DIAGNOSTICVARIABLES_HPP +#define DIAGNOSTICVARIABLES_HPP + +// assign an enum to each variable +// todo: here add pressure, density, etc. +enum +{ + c_Ham, + + c_Mom, + + NUM_DIAGNOSTIC_VARS +}; + +namespace DiagnosticVariables +{ +static const std::array variable_names = { + "Ham", + + "Mom"}; +} + +#endif /* DIAGNOSTICVARIABLES_HPP */ diff --git a/Examples/PerfectFluid/EoS.hpp b/Examples/PerfectFluid/EoS.hpp new file mode 100644 index 000000000..eebf6f7a0 --- /dev/null +++ b/Examples/PerfectFluid/EoS.hpp @@ -0,0 +1,40 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef EOS_HPP_ +#define EOS_HPP_ + +#include "simd.hpp" + +class EoS +{ + public: + struct params_t + { + double eos_w; + }; + + private: + params_t m_params; + + public: + //! The constructor + EoS(params_t a_params) : m_params(a_params) {} + + //! Set the potential function for the scalar field here + template class vars_t> + void compute_eos(data_t &P_of_rho0, data_t &dPdrho0, + const vars_t &vars) const + { + // The pressure value at rho0 + // w rho0 + P_of_rho0 = m_params.eos_w * vars.rho0; + + // The pressure gradient at rho0 + dPdrho0 = m_params.eos_w; + } +}; + +#endif /* EOS_HPP_ */ diff --git a/Examples/PerfectFluid/GNUmakefile b/Examples/PerfectFluid/GNUmakefile new file mode 100644 index 000000000..111bb2375 --- /dev/null +++ b/Examples/PerfectFluid/GNUmakefile @@ -0,0 +1,25 @@ +# -*- Mode: Makefile -*- + +### This makefile produces an executable for each name in the `ebase' +### variable using the libraries named in the `LibNames' variable. + +# Included makefiles need an absolute path to the Chombo installation +# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) + +GRCHOMBO_SOURCE = ../../Source + +ebase := Main_PerfectFluid + +LibNames := AMRTimeDependent AMRTools BoxTools + +src_dirs := $(GRCHOMBO_SOURCE)/utils \ + $(GRCHOMBO_SOURCE)/simd \ + $(GRCHOMBO_SOURCE)/CCZ4 \ + $(GRCHOMBO_SOURCE)/Matter \ + $(GRCHOMBO_SOURCE)/BoxUtils \ + $(GRCHOMBO_SOURCE)/GRChomboCore \ + $(GRCHOMBO_SOURCE)/TaggingCriteria \ + $(GRCHOMBO_SOURCE)/AMRInterpolator \ + $(GRCHOMBO_SOURCE)/InitialConditions/BlackHoles + +include $(CHOMBO_HOME)/mk/Make.test diff --git a/Examples/PerfectFluid/InitialFluidData.hpp b/Examples/PerfectFluid/InitialFluidData.hpp new file mode 100644 index 000000000..25392abe7 --- /dev/null +++ b/Examples/PerfectFluid/InitialFluidData.hpp @@ -0,0 +1,69 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef INITIALFLUIDDATA_HPP_ +#define INITIALFLUIDDATA_HPP_ + +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "MatterCCZ4RHS.hpp" +#include "PerfectFluid.hpp" +#include "Tensor.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total no. components +#include "VarsTools.hpp" +#include "simd.hpp" +//#include "PrimitiveRecovety.hpp" + +//! Class which sets the initial scalar field matter config +class InitialFluidData +{ + public: + //! A structure for the input params for fluid properties and initial + //! conditions + struct params_t + { + double amplitude; // + std::array + center; //!< Centre of perturbation in initial SF bubble + double delta; + double width; //!< Width of bump in initial SF bubble + }; + + //! The constructor + InitialFluidData(params_t a_params, double a_dx) + : m_dx(a_dx), m_params(a_params) + { + } + + //! Function to compute the value of all the initial vars on the grid + template void compute(Cell current_cell) const + { + // where am i? + Coordinates coords(current_cell, m_dx, m_params.center); + data_t rr = coords.get_radius(); + data_t rr2 = rr * rr; + + // calculate the field value + data_t rho0 = m_params.amplitude * + (exp(-pow(rr / m_params.width, 2.0))) + m_params.delta; + data_t v2 = 0.; + data_t eps = 0.; + data_t P = rho0*(1.+eps)/3.; + data_t WW = 1./(1. - v2); + data_t hh = 1. + eps + P/rho0; + + data_t D0 = rho0*sqrt(WW); + data_t Ec0 = rho0 * hh * WW - P - D0; + // store the vars + current_cell.store_vars(D0, c_D); + current_cell.store_vars(Ec0, c_Ec); + } + + protected: + double m_dx; + const params_t m_params; //!< The matter initial condition params +}; + +#endif /* INITIALFLUIDDATA_HPP_ */ diff --git a/Examples/PerfectFluid/Main_PerfectFluid.cpp b/Examples/PerfectFluid/Main_PerfectFluid.cpp new file mode 100644 index 000000000..59a90f220 --- /dev/null +++ b/Examples/PerfectFluid/Main_PerfectFluid.cpp @@ -0,0 +1,64 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +// Chombo includes +#include "parstream.H" //Gives us pout() + +// System includes +#include + +// Our general includes +#include "DefaultLevelFactory.hpp" +#include "GRAMR.hpp" +#include "GRParmParse.hpp" +#include "SetupFunctions.hpp" +#include "SimulationParameters.hpp" + +// Problem specific includes: +#include "PerfectFluidLevel.hpp" + +// Chombo namespace +#include "UsingNamespace.H" + +int runGRChombo(int argc, char *argv[]) +{ + // Load the parameter file and construct the SimulationParameter class + // To add more parameters edit the SimulationParameters file. + char *in_file = argv[1]; + GRParmParse pp(argc - 2, argv + 2, NULL, in_file); + SimulationParameters sim_params(pp); + + if (sim_params.just_check_params) + return 0; + + // The line below selects the problem that is simulated + // (To simulate a different problem, define a new child of AMRLevel + // and an associated LevelFactory) + GRAMR gr_amr; + DefaultLevelFactory perfect_fluid_level_fact(gr_amr, + sim_params); + setupAMRObject(gr_amr, perfect_fluid_level_fact); + + // Engage! Run the evolution + gr_amr.run(sim_params.stop_time, sim_params.max_steps); + gr_amr.conclude(); + + return 0; +} + +int main(int argc, char *argv[]) +{ + mainSetup(argc, argv); + + int status = runGRChombo(argc, argv); + + if (status == 0) + pout() << "GRChombo finished." << std::endl; + else + pout() << "GRChombo failed with return code " << status << std::endl; + + mainFinalize(); + return status; +} diff --git a/Examples/PerfectFluid/Minkowski.hpp b/Examples/PerfectFluid/Minkowski.hpp new file mode 100644 index 000000000..e1d9a664c --- /dev/null +++ b/Examples/PerfectFluid/Minkowski.hpp @@ -0,0 +1,69 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef MINKOWSKI_HPP_ +#define MINKOWSKI_HPP_ + +#include "ADMConformalVars.hpp" +#include "Cell.hpp" +#include "CoordinateTransformations.hpp" +#include "Coordinates.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components +#include "VarsTools.hpp" +#include "simd.hpp" + +//! Class which computes the Kerr initial conditions per arXiv 1401.1548 +class Minkowski +{ + // Use the variable definition in CCZ4 + template + using Vars = ADMConformalVars::VarsWithGauge; + + public: + //! Stuct for the params of the Kerr BH + struct params_t + { + double mass; //!<< The mass of the Kerr BH + std::array center; //!< The center of the Kerr BH + double spin; //!< The spin param a = J/M, so 0 <= |a| <= M + }; + + protected: + double m_dx; + params_t m_params; + + public: + Minkowski(params_t a_params, double a_dx) : m_dx(a_dx), m_params(a_params) + + { + // check this spin param is sensible + if (std::abs(m_params.spin) > m_params.mass) + { + MayDay::Error("The spin parameter must satisfy |a| <= M"); + } + } + + template void compute(Cell current_cell) const; + + protected: + //! Function which computes the components of the metric in spherical coords + template + void compute_kerr( + Tensor<2, data_t> + &spherical_g, //!<< The spatial metric in spherical coords + Tensor<2, data_t> + &spherical_K, //!<< The extrinsic curvature in spherical coords + Tensor<1, data_t> + &spherical_shift, //!<< The spherical components of the shift + data_t &kerr_lapse, //!<< The lapse for the kerr solution + const Coordinates coords //!<< Coords of current cell + ) const; +}; + +#include "Minkowski.impl.hpp" + +#endif /* MINKOWSKI_HPP_ */ diff --git a/Examples/PerfectFluid/Minkowski.impl.hpp b/Examples/PerfectFluid/Minkowski.impl.hpp new file mode 100644 index 000000000..27ccf4934 --- /dev/null +++ b/Examples/PerfectFluid/Minkowski.impl.hpp @@ -0,0 +1,151 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(MINKOWSKI_HPP_) +#error "This file should only be included through Minkowski.hpp" +#endif + +#ifndef MINKOWSKI_IMPL_HPP_ +#define MINKOWSKI_IMPL_HPP_ + +#include "DimensionDefinitions.hpp" + +// Computes semi-isotropic Kerr solution as detailed in Liu, Etienne and Shapiro +// 2010, arxiv gr-qc/1001.4077 +template void Minkowski::compute(Cell current_cell) const +{ + // set up vars for the metric and extrinsic curvature, shift and lapse in + // spherical coords + Tensor<2, data_t> spherical_g; + Tensor<2, data_t> spherical_K; + Tensor<1, data_t> spherical_shift; + data_t kerr_lapse; + + // The cartesian variables and coords + Vars vars; + Coordinates coords(current_cell, m_dx, m_params.center); + + // Compute the components in spherical coords as per 1401.1548 + compute_kerr(spherical_g, spherical_K, spherical_shift, kerr_lapse, coords); + + // work out where we are on the grid + data_t x = coords.x; + double y = coords.y; + double z = coords.z; + + using namespace CoordinateTransformations; + // Convert spherical components to cartesian components using coordinate + // transforms + vars.A = spherical_to_cartesian_LL(spherical_K, x, y, z); + vars.shift = spherical_to_cartesian_U(spherical_shift, x, y, z); + + using namespace TensorAlgebra; + // Convert to BSSN vars + FOR(i, j) { vars.h[i][j] = delta(i, j); } + data_t deth = compute_determinant(vars.h); + auto h_UU = compute_inverse_sym(vars.h); + vars.chi = pow(deth, -1. / 3.); + + // transform extrinsic curvature into A and TrK - note h is still non + // conformal version which is what we need here + vars.K = 0.0; //compute_trace(vars.A, h_UU); + make_trace_free(vars.A, vars.h, h_UU); + + // Make conformal + FOR(i, j) + { + vars.h[i][j] *= vars.chi; + vars.A[i][j] *= vars.chi; + } + + // use a pre collapsed lapse, could also use analytic one + // vars.lapse = kerr_lapse; + vars.lapse = pow(vars.chi, 0.5); + + // Populate the variables on the grid + // NB We stil need to set Gamma^i which is NON ZERO + // but we do this via a separate class/compute function + // as we need the gradients of the metric which are not yet available + current_cell.store_vars(vars); +} + +template +void Minkowski::compute_kerr(Tensor<2, data_t> &spherical_g, + Tensor<2, data_t> &spherical_K, + Tensor<1, data_t> &spherical_shift, + data_t &kerr_lapse, + const Coordinates coords) const +{ + // Kerr black hole params - mass M and spin a + double M = m_params.mass; + double a = m_params.spin; + + // work out where we are on the grid + data_t x = coords.x; + double y = coords.y; + double z = coords.z; + + // the radius, subject to a floor + data_t r = coords.get_radius(); + data_t r2 = r * r; + + // the radius in xy plane, subject to a floor + data_t rho2 = simd_max(x * x + y * y, 1e-12); + data_t rho = sqrt(rho2); + + // calculate useful position quantities + data_t cos_theta = z / r; + data_t sin_theta = rho / r; + data_t cos_theta2 = cos_theta * cos_theta; + data_t sin_theta2 = sin_theta * sin_theta; + + // calculate useful metric quantities + double r_plus = M + sqrt(M * M - a * a); + double r_minus = M - sqrt(M * M - a * a); + + // The Boyer-Lindquist coordinate + data_t r_BL = r * pow(1.0 + 0.25 * r_plus / r, 2.0); + + // Other useful quantities per 1001.4077 + data_t Sigma = r_BL * r_BL + a * a * cos_theta2; + data_t Delta = r_BL * r_BL - 2.0 * M * r_BL + a * a; + // In the paper this is just 'A', but not to be confused with A_ij + data_t AA = pow(r_BL * r_BL + a * a, 2.0) - Delta * a * a * sin_theta2; + // The rr component of the conformal spatial matric + data_t gamma_rr = + Sigma * pow(r + 0.25 * r_plus, 2.0) / (r * r2 * (r_BL - r_minus)); + + // Metric in semi isotropic Kerr-Schild coordinates, r, theta (t or th), phi + // (p) + FOR(i, j) { spherical_g[i][j] = 0.0; } + spherical_g[0][0] = gamma_rr; // gamma_rr + spherical_g[1][1] = Sigma; // gamma_tt + spherical_g[2][2] = AA / Sigma * sin_theta2; // gamma_pp + + // Extrinsic curvature + FOR(i, j) { spherical_K[i][j] = 0.0; } + + // set non zero elements of Krtp - K_rp, K_tp + //spherical_K[0][2] = + // a * M * sin_theta2 / (Sigma * sqrt(AA * Sigma)) * + // (3.0 * pow(r_BL, 4.0) + 2 * a * a * r_BL * r_BL - pow(a, 4.0) - + // a * a * (r_BL * r_BL - a * a) * sin_theta2) * + // (1.0 + 0.25 * r_plus / r) / sqrt(r * r_BL - r * r_minus); + //spherical_K[2][0] = spherical_K[0][2]; + //spherical_K[2][1] = -2.0 * pow(a, 3.0) * M * r_BL * cos_theta * sin_theta * + // sin_theta2 / (Sigma * sqrt(AA * Sigma)) * + // (r - 0.25 * r_plus) * sqrt(r_BL / r - r_minus / r); + //spherical_K[1][2] = spherical_K[2][1]; + + // set the analytic lapse + kerr_lapse = sqrt(Delta * Sigma / AA); + + // set the shift (only the phi component is non zero) + spherical_shift[0] = 0.0; + spherical_shift[1] = 0.0; + spherical_shift[2] = -2.0 * M * a * r_BL / AA; +} + +#endif /* MINKOWSKI_IMPL_HPP_ */ diff --git a/Examples/PerfectFluid/PerfectFluid.hpp b/Examples/PerfectFluid/PerfectFluid.hpp new file mode 100644 index 000000000..6d66de871 --- /dev/null +++ b/Examples/PerfectFluid/PerfectFluid.hpp @@ -0,0 +1,114 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef PERFECTFLUID_HPP_ +#define PERFECTFLUID_HPP_ + +#include "CCZ4Geometry.hpp" +#include "DefaultEoS.hpp" +#include "DimensionDefinitions.hpp" +#include "FourthOrderDerivatives.hpp" +#include "WENODerivatives.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS, total num of components +#include "VarsTools.hpp" + +//! Calculates the matter type specific elements such as the EMTensor and +// matter evolution +/*! + This class is an example of a matter_t object which calculates the + matter type specific elements for the RHS update and the evaluation + of the constraints. This includes the Energy Momentum Tensor, and + the matter evolution terms. In this case, a scalar field, + the matter elements are phi and (minus) its conjugate momentum, Pi. + It is templated over a potential function potential_t which the + user must specify in a class, although a default is provided which + sets dVdphi and V_of_phi to zero. + It assumes minimal coupling of the field to gravity. + \sa MatterCCZ4(), ConstraintsMatter() +*/ +template class PerfectFluid +{ + protected: + //! The local copy of the potential + eos_t my_eos; + + public: + //! Constructor of class ScalarField, inputs are the matter parameters. + PerfectFluid(const eos_t a_eos) : my_eos(a_eos) {} + + //! Structure containing the rhs variables for the matter fields + template struct Vars + { + data_t D; + Tensor<1, data_t> Sj; + data_t Ec; + Tensor<1, data_t> vi; + data_t rho0; + data_t eps; + + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + using namespace VarsTools; // define_enum_mapping is part of + // VarsTools + define_enum_mapping(mapping_function, c_D, D); + define_enum_mapping(mapping_function, + GRInterval(), Sj); + define_enum_mapping(mapping_function, c_Ec, Ec); + define_enum_mapping(mapping_function, + GRInterval(), vi); + define_enum_mapping(mapping_function, c_rho0, rho0); + define_enum_mapping(mapping_function, c_eps, eps); + } + }; + + //! Structure containing the rhs variables for the matter fields requiring + //! 2nd derivs + template struct Diff2Vars + { + data_t rho0; + + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + VarsTools::define_enum_mapping(mapping_function, c_rho0, rho0); + // VarsTools::define_enum_mapping( + // mapping_function, GRInterval(), Avec); + } + }; + + //! The function which calculates the EM Tensor, given the vars and + //! derivatives + template class vars_t> + emtensor_t compute_emtensor( + const vars_t &vars, //!< the value of the variables + const vars_t> &d1, //!< the value of the 1st derivs + const Tensor<2, data_t> &h_UU, //!< the inverse metric (raised indices) + const Tensor<3, data_t> &chris_ULL) + const; //!< the conformal christoffel symbol + + //! The function which adds in the RHS for the matter field vars + template class vars_t, + template class diff2_vars_t, + template class rhs_vars_t> + void add_matter_rhs( + rhs_vars_t &rhs, //!< value of the RHS for all vars + const vars_t &vars, //!< value of the variables + const vars_t> &d1, //!< value of the 1st derivs + const diff2_vars_t> &d2, //!< value of the 2nd derivs + const vars_t &advec) + const; //!< the value of the advection terms + +}; + +#include "PerfectFluid.impl.hpp" + +#endif /* PERFECTFLUID_HPP_ */ diff --git a/Examples/PerfectFluid/PerfectFluid.impl.hpp b/Examples/PerfectFluid/PerfectFluid.impl.hpp new file mode 100644 index 000000000..1cc5b19f7 --- /dev/null +++ b/Examples/PerfectFluid/PerfectFluid.impl.hpp @@ -0,0 +1,96 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(PERFECTFLUID_HPP_) +#error "This file should only be included through PerfectFluid.hpp" +#endif + +#ifndef PERFECTFLUID_IMPL_HPP_ +#define PERFECTFLUID_IMPL_HPP_ + +// Calculate the stress energy tensor elements +template +template class vars_t> +emtensor_t PerfectFluid::compute_emtensor( + const vars_t &vars, const vars_t> &d1, + const Tensor<2, data_t> &h_UU, const Tensor<3, data_t> &chris_ULL) const +{ + emtensor_t out; + + // Set up EoS + data_t P_of_rho0 = 0.0; + data_t dPdrho0 = 0.0; + my_eos.compute_eos(P_of_rho0, dPdrho0, vars); + + // Useful quantities + data_t v2 = 0.0; + FOR(i, j) + { + v2 += vars.h[i][j] * vars.vi[i] * vars.vi[j] / vars.chi; + } + data_t WW = 1./(1. - v2); + data_t hh = 1. + vars.eps + P_of_rho0/vars.rho0; + + Tensor<1, data_t> vi_D; + FOR(i) + { + vi_D[i] = 0; + FOR(j) + { + vi_D[i] += vars.h[i][j] * vars.vi[j] / vars.chi; + } + } + + // Calculate components of EM Tensor + // S_ij = T_ij + FOR(i, j) + { + out.Sij[i][j] = vars.rho0 * hh * WW * vi_D[i] * vi_D[j] + + vars.h[i][j] * P_of_rho0 / vars.chi; + } + + // S = Tr_S_ij + out.S = vars.chi * TensorAlgebra::compute_trace(out.Sij, h_UU); + + // S_i (note lower index) = - n^a T_ai + FOR(i) { out.Si[i] = vars.rho0 * hh * WW * vi_D[i]; } + + // rho = n^a n^b T_ab + out.rho = vars.rho0 * hh * WW - P_of_rho0; +} + +template +template class vars_t, + template class diff2_vars_t, + template class rhs_vars_t> +void PerfectFluid::add_matter_rhs( + rhs_vars_t &rhs, const vars_t &vars, + const vars_t> &d1, + const diff2_vars_t> &d2, + const vars_t &advec) const +{ + using namespace TensorAlgebra; + + const auto h_UU = compute_inverse_sym(vars.h); + const auto chris = compute_christoffel(d1.h, h_UU); + + data_t P_of_rho0 = 0.0; + data_t dPdrho0 = 0.0; + my_eos.compute_eos(P_of_rho0, dPdrho0, vars); + + // evolution equations for the fluid conservative variables + rhs.D = 0.0; + rhs.Ec = 0.0; + rhs.rho0 = 0.0; + rhs.eps = 0.0; + + FOR(i) + { + rhs.Sj[i] = 0.0; + rhs.vi[i] = 0.0; + } +} + +#endif /* SCALARFIELD_IMPL_HPP_ */ diff --git a/Examples/PerfectFluid/PerfectFluidLevel.cpp b/Examples/PerfectFluid/PerfectFluidLevel.cpp new file mode 100644 index 000000000..571f08061 --- /dev/null +++ b/Examples/PerfectFluid/PerfectFluidLevel.cpp @@ -0,0 +1,132 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +// General includes common to most GR problems +#include "PerfectFluidLevel.hpp" +#include "BoxLoops.hpp" +#include "NanCheck.hpp" +#include "PositiveChiAndAlpha.hpp" +#include "SixthOrderDerivatives.hpp" +#include "WENODerivatives.hpp" +#include "TraceARemoval.hpp" + +// For RHS update +#include "MatterCCZ4RHS.hpp" + +// For constraints calculation +#include "NewMatterConstraints.hpp" + +// For tag cells +#include "FixedGridsTaggingCriterion.hpp" + +// Problem specific includes +#include "ComputePack.hpp" +#include "GammaCalculator.hpp" +#include "InitialFluidData.hpp" +#include "Minkowski.hpp" +#include "EoS.hpp" +#include "PerfectFluid.hpp" +#include "SetValue.hpp" + +// Things to do at each advance step, after the RK4 is calculated +void PerfectFluidLevel::specificAdvance() +{ + // Enforce trace free A_ij and positive chi and alpha + BoxLoops::loop( + make_compute_pack(TraceARemoval(), + PositiveChiAndAlpha(m_p.min_chi, m_p.min_lapse)), + m_state_new, m_state_new, INCLUDE_GHOST_CELLS); + + // Check for nan's + if (m_p.nan_check) + BoxLoops::loop(NanCheck(), m_state_new, m_state_new, + EXCLUDE_GHOST_CELLS, disable_simd()); +} + +// Initial data for field and metric variables +void PerfectFluidLevel::initialData() +{ + CH_TIME("PerfectFluidLevel::initialData"); + if (m_verbosity) + pout() << "PerfectFluidLevel::initialData " << m_level << endl; + + // First set everything to zero then initial conditions for scalar field - + // here a Kerr BH and a scalar field profile + BoxLoops::loop( + make_compute_pack(SetValue(0.), Minkowski(m_p.kerr_params, m_dx), + InitialFluidData(m_p.initial_params, m_dx)), + m_state_new, m_state_new, INCLUDE_GHOST_CELLS); + + fillAllGhosts(); + BoxLoops::loop(GammaCalculator(m_dx), m_state_new, m_state_new, + EXCLUDE_GHOST_CELLS); +} + +#ifdef CH_USE_HDF5 +// Things to do before outputting a checkpoint file +void PerfectFluidLevel::prePlotLevel() +{ + fillAllGhosts(); + EoS eos(m_p.eos_params); + PerfectFluidEoS perfect_fluid(eos); + BoxLoops::loop( + MatterConstraints( + perfect_fluid, m_dx, m_p.G_Newton, c_Ham, Interval(c_Mom, c_Mom)), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS); +} +#endif + +// Things to do in RHS update, at each RK4 step +void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, + const double a_time) +{ + // Enforce trace free A_ij and positive chi and alpha + BoxLoops::loop( + make_compute_pack(TraceARemoval(), + PositiveChiAndAlpha(m_p.min_chi, m_p.min_lapse)), + a_soln, a_soln, INCLUDE_GHOST_CELLS); + + // Calculate MatterCCZ4 right hand side with matter_t = ScalarField + EoS eos(m_p.eos_params); + PerfectFluidEoS perfect_fluid(eos); + if (m_p.max_spatial_derivative_order == 4) + { + MatterCCZ4RHS + my_ccz4_matter(perfect_fluid, m_p.ccz4_params, m_dx, m_p.sigma, + m_p.formulation, m_p.G_Newton); + BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS); + } + else if (m_p.max_spatial_derivative_order == 6) + { + MatterCCZ4RHS + my_ccz4_matter(perfect_fluid, m_p.ccz4_params, m_dx, m_p.sigma, + m_p.formulation, m_p.G_Newton); + BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS); + } +} + +// Things to do at ODE update, after soln + rhs +void PerfectFluidLevel::specificUpdateODE(GRLevelData &a_soln, + const GRLevelData &a_rhs, Real a_dt) +{ + // Enforce trace free A_ij + BoxLoops::loop(TraceARemoval(), a_soln, a_soln, INCLUDE_GHOST_CELLS); +} + +void PerfectFluidLevel::preTagCells() +{ + // we don't need any ghosts filled for the fixed grids tagging criterion + // used here so don't fill any +} + +void PerfectFluidLevel::computeTaggingCriterion(FArrayBox &tagging_criterion, + const FArrayBox ¤t_state) +{ + BoxLoops::loop( + FixedGridsTaggingCriterion(m_dx, m_level, m_p.L, m_p.center), + current_state, tagging_criterion); +} diff --git a/Examples/PerfectFluid/PerfectFluidLevel.hpp b/Examples/PerfectFluid/PerfectFluidLevel.hpp new file mode 100644 index 000000000..a501030be --- /dev/null +++ b/Examples/PerfectFluid/PerfectFluidLevel.hpp @@ -0,0 +1,60 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef PERFECTFLUIDLEVEL_HPP_ +#define PERFECTFLUIDLEVEL_HPP_ + +#include "DefaultLevelFactory.hpp" +#include "GRAMRLevel.hpp" +// Problem specific includes +#include "EoS.hpp" +#include "PerfectFluid.hpp" + +//! A class for the evolution of a scalar field, minimally coupled to gravity +/*! + The class takes some initial data for a scalar field (variables phi and Pi) + and evolves it using the CCZ4 equations. It is possible to specify an + initial period of relaxation for the conformal factor chi, for non analytic + initial conditions (for example, a general field configuration at a moment of + time symmetry assuming conformal flatness). \sa MatterCCZ4(), + ConstraintsMatter(), ScalarField(), RelaxationChi() +*/ +class PerfectFluidLevel : public GRAMRLevel +{ + friend class DefaultLevelFactory; + // Inherit the contructors from GRAMRLevel + using GRAMRLevel::GRAMRLevel; + + // Typedef for fluid + typedef PerfectFluid PerfectFluidEoS; + + //! Things to do at the end of the advance step, after RK4 calculation + virtual void specificAdvance(); + + //! Initialize data for the field and metric variables + virtual void initialData(); + +#ifdef CH_USE_HDF5 + //! routines to do before outputting plot file + virtual void prePlotLevel(); +#endif + + //! RHS routines used at each RK4 step + virtual void specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, + const double a_time); + + //! Things to do in UpdateODE step, after soln + rhs update + virtual void specificUpdateODE(GRLevelData &a_soln, + const GRLevelData &a_rhs, Real a_dt); + + /// Things to do before tagging cells (i.e. filling ghosts) + virtual void preTagCells() override; + + //! Tell Chombo how to tag cells for regridding + virtual void computeTaggingCriterion(FArrayBox &tagging_criterion, + const FArrayBox ¤t_state); +}; + +#endif /* PERFECTFLUIDLEVEL_HPP_ */ diff --git a/Examples/PerfectFluid/PrimitiveRecovey.hpp b/Examples/PerfectFluid/PrimitiveRecovey.hpp new file mode 100644 index 000000000..a13c3ff45 --- /dev/null +++ b/Examples/PerfectFluid/PrimitiveRecovey.hpp @@ -0,0 +1,69 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef PRIMITIVERECOVERY_HPP_ +#define PRIMITIVERECOVERY_HPP_ + +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "simd.hpp" + +namespace PrimitiveRecovery +{ +// Primitive to conservative variables +template +Tensor<1, data_t> PtoC(const double p) +{ + Tensor<1, data_t> q; + const double eps = p[0]; + const double ux = p[1]; + const double uy = p[2]; + const double nn = p[3]; + double ux2 = ux*ux; + double uy2 = uy*uy; + double ut = sqrt(1.+ux2+uy2); + + // Ttt + q[0] = (eps*(3.+4.*(ux2+uy2)))/3.; + // Ttx + q[1] = (4.*ux*ut*eps)/3.; + // Tty + q[2] = (4.*uy*ut*eps)/3.; + // Jt + q[3] = nn*ut; + + return q; +} + +// Conservative to primitive variables + template +Tensor<1, data_t> CtoP(const double q) +{ + Tensor<1, data_t> p; + const double Ttt = q[0]; + const double Ttx = q[1]; + const double Tty = q[2]; + const double Jt = q[3]; + double Ttt2 = Ttt*Ttt; + double Ttx2 = Ttx*Ttx; + double Tty2 = Tty*Tty; + double sqrt1 = sqrt(4.*Ttt2-3.*(Ttx2+Tty2)); + double sqrt2 = sqrt(2.*Ttt2-3.*(Ttx2+Tty2)+Ttt*sqrt1); + double ut; + + // eps + p[0] = -Ttt + sqrt1; + // ux + p[1] = (3.*Ttx)/(2.*sqrt2); + // uy + p[2] = (3.*Tty)/(2.*sqrt2); + // ut + ut = sqrt(1.+p[1]*p[1]+p[2]*p[2]); + p[3] = Jt/ut; + return p; +} +} +#endif /* PRIMITIVERECOVERY_HPP_ */ diff --git a/Examples/PerfectFluid/SimulationParameters.hpp b/Examples/PerfectFluid/SimulationParameters.hpp new file mode 100644 index 000000000..018e2c608 --- /dev/null +++ b/Examples/PerfectFluid/SimulationParameters.hpp @@ -0,0 +1,80 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMULATIONPARAMETERS_HPP_ +#define SIMULATIONPARAMETERS_HPP_ + +// General includes +#include "GRParmParse.hpp" +#include "SimulationParametersBase.hpp" + +// Problem specific includes: +#include "InitialFluidData.hpp" +#include "Minkowski.hpp" +#include "EoS.hpp" + +class SimulationParameters : public SimulationParametersBase +{ + public: + SimulationParameters(GRParmParse &pp) : SimulationParametersBase(pp) + { + // read the problem specific params + read_params(pp); + check_params(); + } + + void read_params(GRParmParse &pp) + { + // Initial scalar field data + initial_params.center = + center; // already read in SimulationParametersBase + pp.load("G_Newton", G_Newton, + 0.0); // for now the example neglects backreaction + pp.load("fluid_amplitude", initial_params.amplitude, 0.1); + pp.load("fluid_delta", initial_params.delta, 0.0); + pp.load("fluid_width", initial_params.width, 1.0); + pp.load("eos_w", eos_params.eos_w, 1./3.); + + // Initial Kerr data + pp.load("kerr_mass", kerr_params.mass, 1.0); + pp.load("kerr_spin", kerr_params.spin, 0.0); + pp.load("kerr_center", kerr_params.center, center); + } + + void check_params() + { + // warn_parameter("scalar_mass", potential_params.scalar_mass, + // potential_params.scalar_mass < + // 0.2 / coarsest_dx / dt_multiplier, + // "oscillations of scalar field do not appear to be " + // "resolved on coarsest level"); + //warn_parameter("scalar_width", initial_params.width, + // initial_params.width < 0.5 * L, + // "is greater than half the domain size"); + warn_parameter("kerr_mass", kerr_params.mass, kerr_params.mass >= 0.0, + "should be >= 0.0"); + check_parameter("kerr_spin", kerr_params.spin, + std::abs(kerr_params.spin) <= kerr_params.mass, + "must satisfy |a| <= M = " + + std::to_string(kerr_params.mass)); + FOR(idir) + { + std::string name = "kerr_center[" + std::to_string(idir) + "]"; + warn_parameter( + name, kerr_params.center[idir], + (kerr_params.center[idir] >= 0) && + (kerr_params.center[idir] <= (ivN[idir] + 1) * coarsest_dx), + "should be within the computational domain"); + } + } + + // Initial data for matter and potential and BH + double G_Newton; + InitialFluidData::params_t initial_params; + EoS::params_t eos_params; + Minkowski::params_t kerr_params; +}; + +#endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Examples/PerfectFluid/UserVariables.hpp b/Examples/PerfectFluid/UserVariables.hpp new file mode 100644 index 000000000..82b49d4ff --- /dev/null +++ b/Examples/PerfectFluid/UserVariables.hpp @@ -0,0 +1,45 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef USERVARIABLES_HPP +#define USERVARIABLES_HPP + +#include "ArrayTools.hpp" +#include "CCZ4UserVariables.hpp" +#include "DiagnosticVariables.hpp" + +// assign an enum to each variable +enum +{ + // todo: Here add conservative variables + // Note that it is important that the first enum value is set to 1 more than + // the last CCZ4 var enum + c_D = NUM_CCZ4_VARS, // matter field added + c_Sj1, + c_Sj2, + c_Sj3, + c_Ec, + c_vi1, + c_vi2, + c_vi3, + c_rho0, + c_eps, + + NUM_VARS +}; + +namespace UserVariables +{ +static const std::array +user_variable_names = {"D", "Sj1", "Sj2", "Sj3", "Ec", + "vi1", "vi2", "vi3", "rho0", "eps"}; + +static const std::array variable_names = + ArrayTools::concatenate(ccz4_variable_names, user_variable_names); +} // namespace UserVariables + +#include "UserVariables.inc.hpp" + +#endif /* USERVARIABLES_HPP */ diff --git a/Examples/PerfectFluid/WENODerivatives.hpp b/Examples/PerfectFluid/WENODerivatives.hpp new file mode 100644 index 000000000..c7804f752 --- /dev/null +++ b/Examples/PerfectFluid/WENODerivatives.hpp @@ -0,0 +1,415 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef WENODERIVATIVES_HPP_ +#define WENODERIVATIVES_HPP_ + +#include "Cell.hpp" +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "UserVariables.hpp" +#include + +class WENODerivatives +{ + public: + const double m_dx; + + private: + const double m_one_over_dx; + const double m_one_over_dx2; + + public: + WENODerivatives(double dx) + : m_dx(dx), m_one_over_dx(1 / dx), m_one_over_dx2(1 / (dx * dx)) + { + } + + template + ALWAYS_INLINE data_t diff1(const double *in_ptr, const int idx, + const int stride) const + { + auto in = SIMDIFY(in_ptr); + + data_t weight_far = 8.33333333333333333333e-2; + data_t weight_near = 6.66666666666666666667e-1; + + // NOTE: if you have been sent here by the debugger because of + // EXC_BAD_ACCESS or something similar you might be trying to take + // derivatives without ghost points. + return (weight_far * in[idx - 2 * stride] - + weight_near * in[idx - stride] + + weight_near * in[idx + stride] - + weight_far * in[idx + 2 * stride]) * + m_one_over_dx; + } + + // Writes directly into the vars object - use this wherever possible + template class vars_t> + void diff1(vars_t> &d1, const Cell ¤t_cell, + int direction) const + { + const int stride = + current_cell.get_box_pointers().m_in_stride[direction]; + const int in_index = current_cell.get_in_index(); + d1.enum_mapping([&](const int &ivar, Tensor<1, data_t> &var) { + var[direction] = + diff1(current_cell.get_box_pointers().m_in_ptr[ivar], + in_index, stride); + }); + } + + /// Calculates all first derivatives and returns as variable type specified + /// by the template parameter + template