diff --git a/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp b/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp new file mode 100644 index 000000000..47dfb2574 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp @@ -0,0 +1,42 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef CONSERVEDQUANTITIES_HPP_ +#define CONSERVEDQUANTITIES_HPP_ + +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "simd.hpp" + +namespace ConservedQuantities +{ +// Primitive to conservative variables +template class vars_t> +void PtoC(const data_t P_of_rho, vars_t &vars) +{ + 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; } + } + + data_t v2 = 0.; + FOR(i) v2 += vars.vi[i] * vi_D[i]; + + data_t WW = 1. / (1. - v2); + data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + + vars.D = vars.rho * sqrt(WW); + vars.tau = vars.rho * hh * WW - P_of_rho - vars.D; + + // S_j (note lower index) = - n^a T_ai + FOR(i) { vars.Sj[i] = vars.rho * hh * WW * vi_D[i]; } +} + +} // namespace ConservedQuantities +#endif /* CONSERVEDQUANTITIES_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/DefaultEoS.hpp b/Examples/Fluid_Kerr_merged/DefaultEoS.hpp new file mode 100644 index 000000000..41b53e4f9 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/DefaultEoS.hpp @@ -0,0 +1,26 @@ +/* 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 pressure of the perfect fluid here to zero + template class vars_t> + void compute_eos(data_t &P_of_rho, const vars_t &vars) const + { + P_of_rho = vars.rho * (1. + vars.eps) / 3.; + } +}; + +#endif /* DEFAULTEOS_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/DiagnosticVariables.hpp b/Examples/Fluid_Kerr_merged/DiagnosticVariables.hpp new file mode 100644 index 000000000..d4c0c864a --- /dev/null +++ b/Examples/Fluid_Kerr_merged/DiagnosticVariables.hpp @@ -0,0 +1,29 @@ +/* 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_Mom1, + c_Mom2, + c_Mom3, + + NUM_DIAGNOSTIC_VARS +}; + +namespace DiagnosticVariables +{ +static const std::array variable_names = { + "Ham", + + "Mom1", "Mom2", "Mom3"}; +} + +#endif /* DIAGNOSTICVARIABLES_HPP */ diff --git a/Examples/Fluid_Kerr_merged/EoS.hpp b/Examples/Fluid_Kerr_merged/EoS.hpp new file mode 100644 index 000000000..8aa3a2303 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/EoS.hpp @@ -0,0 +1,29 @@ +/* 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: + //! The constructor + EoS() {} + + //! Set the pressure of Gamma-2 polytrope + template class vars_t> + void compute_eos(data_t &P_of_rho, const vars_t &vars) const + { + // The pressure value as a function of rho + const double K = 1./3.;//100.; + const double n = 1.; + const double Gamma = 1.;// + 1. / n; + P_of_rho = K * pow(vars.rho, Gamma); + } +}; + +#endif /* EOS_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/FluidCCZ4RHS.hpp b/Examples/Fluid_Kerr_merged/FluidCCZ4RHS.hpp new file mode 100644 index 000000000..fd05fc932 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/FluidCCZ4RHS.hpp @@ -0,0 +1,119 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef FLUIDCCZ4RHS_HPP_ +#define FLUIDCCZ4RHS_HPP_ + +#include "CCZ4Geometry.hpp" +#include "CCZ4RHS.hpp" +#include "Cell.hpp" +#include "FourthOrderDerivatives.hpp" +#include "MovingPunctureGauge.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components +#include "VarsTools.hpp" +#include "WENODerivatives.hpp" +#include "simd.hpp" + +//! Calculates RHS using CCZ4 including matter terms, and matter variable +//! evolution +/*! + The class calculates the RHS evolution for all the variables. It inherits + from the CCZ4RHS class, which it uses to do the non matter evolution of + variables. It then adds in the additional matter terms to the CCZ4 evolution + (those including the stress energy tensor), and calculates the evolution of + the matter variables. It does not assume a specific form of matter but is + templated over a matter class matter_t. Please see the class ScalarField as + an example of a matter_t. \sa CCZ4RHS(), ScalarField() +*/ + +template +class FluidCCZ4RHS : public CCZ4RHS +{ + public: + // Use this alias for the same template instantiation as this class + using CCZ4 = CCZ4RHS; + + using params_t = CCZ4_params_t; + + template + using MatterVars = typename matter_t::template Vars; + + // template + // using MatterDiff2Vars = typename matter_t::template Diff2Vars; + + template + using CCZ4Vars = typename CCZ4::template Vars; + + template + using CCZ4Diff2Vars = typename CCZ4::template Diff2Vars; + + // Inherit the variable definitions from CCZ4RHS + matter_t + template + struct Vars : public CCZ4Vars, public MatterVars + { + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + CCZ4Vars::enum_mapping(mapping_function); + MatterVars::enum_mapping(mapping_function); + } + }; + + template struct Diff2Vars : public CCZ4Diff2Vars + // public MatterDiff2Vars + { + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + CCZ4Diff2Vars::enum_mapping(mapping_function); + // MatterDiff2Vars::enum_mapping(mapping_function); + } + }; + + //! Constructor of class MatterCCZ4 + /*! + Inputs are the grid spacing, plus the CCZ4 evolution parameters and a + matter object. It also takes the dissipation parameter sigma, and allows + the formulation to be toggled between CCZ4 and BSSN. The default is CCZ4. + It allows the user to set the value of Newton's constant, which is set to + one by default. + */ + FluidCCZ4RHS(matter_t a_matter, params_t a_params, double a_dx, + double a_sigma, int a_formulation = CCZ4RHS<>::USE_CCZ4, + double a_G_Newton = 1.0); + + //! The compute member which calculates the RHS at each point in the box + //! \sa matter_rhs_equation() + template void compute(Cell current_cell) const; + + protected: + //! The function which adds in the EM Tensor terms to the CCZ4 rhs \sa + //! compute() + template + void add_emtensor_rhs( + Vars + &matter_rhs, //!< the RHS data for each variable at that point. + const Vars &vars, //!< the value of the variables at the point. + const Vars> + &d1 //!< the value of the first derivatives of the variables. + ) const; + + // Class members + weno_t m_weno; + matter_t my_matter; //!< The matter object, e.g. a scalar field. + const double m_G_Newton; //!< Newton's constant, set to one by default. +}; + +#include "FluidCCZ4RHS.impl.hpp" + +#endif /* FLUIDCCZ4RHS_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/FluidCCZ4RHS.impl.hpp b/Examples/Fluid_Kerr_merged/FluidCCZ4RHS.impl.hpp new file mode 100644 index 000000000..3708fac34 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/FluidCCZ4RHS.impl.hpp @@ -0,0 +1,120 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(FLUIDCCZ4RHS_HPP_) +#error "This file should only be included through FluidCCZ4RHS.hpp" +#endif + +#ifndef FLUIDCCZ4RHS_IMPL_HPP_ +#define FLUIDCCZ4RHS_IMPL_HPP_ +#include "DimensionDefinitions.hpp" + +template +FluidCCZ4RHS::FluidCCZ4RHS( + matter_t a_matter, CCZ4_params_t a_params, + double a_dx, double a_sigma, int a_formulation, double a_G_Newton) + : CCZ4RHS(a_params, a_dx, a_sigma, a_formulation, + 0.0 /*No cosmological constant*/), + my_matter(a_matter), m_G_Newton(a_G_Newton) +{ +} + +template +template +void FluidCCZ4RHS::compute( + Cell current_cell) const +{ + // copy data from chombo gridpoint into local variables + const auto matter_vars = current_cell.template load_vars(); + const auto d1 = this->m_deriv.template diff1(current_cell); + const auto d2 = this->m_deriv.template diff2(current_cell); + const auto advec = + this->m_deriv.template advection(current_cell, matter_vars.shift); + const auto lm = this->m_weno.template get_Pface( + current_cell, WENODerivatives::LEFT_MINUS); + // (current_cell,0); + const auto lp = this->m_weno.template get_Pface( + current_cell, WENODerivatives::LEFT_PLUS); + // (current_cell,1); + const auto rm = this->m_weno.template get_Pface( + current_cell, WENODerivatives::RIGHT_MINUS); + //(current_cell,2); + const auto rp = this->m_weno.template get_Pface( + current_cell, WENODerivatives::RIGHT_PLUS); + //(current_cell,3); + + // Call CCZ4 RHS - work out RHS without matter, no dissipation + Vars matter_rhs; + this->rhs_equation(matter_rhs, matter_vars, d1, d2, advec); + + // add RHS matter terms from EM Tensor + add_emtensor_rhs(matter_rhs, matter_vars, d1); + + // add evolution of matter fields themselves + my_matter.add_matter_rhs(matter_rhs, matter_vars, d1, lm, lp, rm, rp); + + // Add dissipation to all terms + this->m_deriv.add_dissipation(matter_rhs, current_cell, this->m_sigma); + + // Write the rhs into the output FArrayBox + current_cell.store_vars(matter_rhs); +} + +// Function to add in EM Tensor matter terms to CCZ4 rhs +template +template +void FluidCCZ4RHS::add_emtensor_rhs( + Vars &matter_rhs, const Vars &matter_vars, + const Vars> &d1) const +{ + using namespace TensorAlgebra; + + const auto h_UU = compute_inverse_sym(matter_vars.h); + const auto chris = compute_christoffel(d1.h, h_UU); + + // Calculate elements of the decomposed stress energy tensor + const auto emtensor = + my_matter.compute_emtensor(matter_vars, d1, h_UU, chris.ULL); + + // Update RHS for K and Theta depending on formulation + if (this->m_formulation == CCZ4RHS<>::USE_BSSN) + { + matter_rhs.K += 4.0 * M_PI * m_G_Newton * matter_vars.lapse * + (emtensor.S + emtensor.rho); + matter_rhs.Theta += 0.0; + } + else + { + matter_rhs.K += 4.0 * M_PI * m_G_Newton * matter_vars.lapse * + (emtensor.S - 3 * emtensor.rho); + matter_rhs.Theta += + -8.0 * M_PI * m_G_Newton * matter_vars.lapse * emtensor.rho; + } + + // Update RHS for other variables + Tensor<2, data_t> Sij_TF = emtensor.Sij; + make_trace_free(Sij_TF, matter_vars.h, h_UU); + + FOR(i, j) + { + matter_rhs.A[i][j] += -8.0 * M_PI * m_G_Newton * matter_vars.chi * + matter_vars.lapse * Sij_TF[i][j]; + } + + FOR(i) + { + data_t matter_term_Gamma = 0.0; + FOR(j) + { + matter_term_Gamma += -16.0 * M_PI * m_G_Newton * matter_vars.lapse * + h_UU[i][j] * emtensor.Si[j]; + } + + matter_rhs.Gamma[i] += matter_term_Gamma; + matter_rhs.B[i] += matter_term_Gamma; + } +} + +#endif /* FLUIDCCZ4RHS_IMPL_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/Fluxes.hpp b/Examples/Fluid_Kerr_merged/Fluxes.hpp new file mode 100644 index 000000000..9114bb7f8 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/Fluxes.hpp @@ -0,0 +1,71 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef FLUXES_HPP_ +#define FLUXES_HPP_ + +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "simd.hpp" + +namespace Fluxes +{ +template class vars_t> +vars_t compute_flux(const data_t P_of_rho, const vars_t &vars, + const int idir) +{ + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + vars_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]; + } + + 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; } + } + + data_t v2 = 0.; + FOR(i) v2 += vars.vi[i] * vi_D[i]; + + data_t WW = 1. / (1. - v2); + data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + + 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_of_rho * 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; + + return out; +} +template class vars_t> +vars_t compute_num_flux(const data_t P_of_rho, + const vars_t &vars, const int idir, + const double lambda, const int sign) +{ + vars_t out = compute_flux(P_of_rho, vars, idir); + out.D += sign * lambda * vars.D; + FOR(j) out.Sj[j] += sign * lambda * vars.Sj[j]; + out.tau += sign * lambda * vars.tau; + // out.Jt += sign * lambda * vars.Jt; + return out; +} + +} // namespace Fluxes +#endif /* FLUXES_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/GNUmakefile b/Examples/Fluid_Kerr_merged/GNUmakefile new file mode 100644 index 000000000..111bb2375 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/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/Fluid_Kerr_merged/InitialFluidData.hpp b/Examples/Fluid_Kerr_merged/InitialFluidData.hpp new file mode 100644 index 000000000..5ef4f4e2f --- /dev/null +++ b/Examples/Fluid_Kerr_merged/InitialFluidData.hpp @@ -0,0 +1,88 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef INITIALFLUIDDATA_HPP_ +#define INITIALFLUIDDATA_HPP_ + +#include "ADMConformalVars.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "FluidCCZ4RHS.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 + +//! Class which sets the initial fluid matter config +class InitialFluidData : public EoS +{ + template + using MetricVars = ADMConformalVars::VarsWithGauge; + + public: + //! A structure for the input params for fluid properties and initial + //! conditions + struct params_t + { + std::array + center; //!< Centre of perturbation in initial SF bubble + double rho0; + double awidth; + double delta; + double spacing; + double *rho_1D; + }; + + //! 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 + //! The constructor + template void compute(Cell current_cell) const; + + //! Structure containing the rhs variables for the matter fields + template struct Vars + { + data_t D; + Tensor<1, data_t> Sj; + data_t tau; + // data_t Jt; + Tensor<1, data_t> vi; + data_t rho; + data_t eps; + // data_t nn; + + /// 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_tau, tau); + // define_enum_mapping(mapping_function, c_Jt, Jt); + define_enum_mapping(mapping_function, GRInterval(), + vi); + define_enum_mapping(mapping_function, c_rho, rho); + define_enum_mapping(mapping_function, c_eps, eps); + // define_enum_mapping(mapping_function, c_nn, nn); + } + }; + + protected: + double m_dx; + const params_t m_params; +}; + +#include "InitialFluidData.impl.hpp" + +#endif /* INITIALFLUIDDATA_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp b/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp new file mode 100644 index 000000000..8f1c575e4 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp @@ -0,0 +1,65 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(INITIALFLUIDDATA_HPP_) +#error "This file should only be included through PrimitiveRecovery.hpp" +#endif + +#ifndef INITIALFLUIDDATA_IMPL_HPP_ +#define INITIALFLUIDDATA_IMPL_HPP_ + +template +void InitialFluidData::compute(Cell current_cell) const +{ + // load vars + auto metric_vars = current_cell.template load_vars(); + auto matter_vars = current_cell.template load_vars(); + VarsTools::assign(matter_vars, 0.); + + // where am i? + Coordinates coords(current_cell, m_dx, m_params.center); + data_t rr = coords.get_radius(); + data_t rr2 = rr * rr; + + 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); + + // calculate the field value + matter_vars.rho = + m_params.rho0;// * (exp(-pow(rr / 2. / m_params.awidth, 2.0))) + m_params.delta; + data_t v2 = 0.; + FOR(i, j) + v2 += metric_vars.h[i][j] * matter_vars.vi[i] * matter_vars.vi[j] / + chi_regularised; + + data_t P_of_rho = 0.; + EoS::compute_eos(P_of_rho, matter_vars); + + data_t WW = 1. / (1. - v2); + data_t hh = 1. + matter_vars.eps + P_of_rho / matter_vars.rho; + + matter_vars.D = matter_vars.rho * sqrt(WW); + matter_vars.tau = matter_vars.rho * hh * WW - P_of_rho - matter_vars.D; + FOR(i) + { + matter_vars.Sj[i] = 0.; + FOR(j) + matter_vars.Sj[i] += matter_vars.rho * hh * WW * metric_vars.h[i][j] * + matter_vars.vi[j] / chi_regularised; + } + + // store the vars + current_cell.store_vars(matter_vars); + using namespace TensorAlgebra; + + metric_vars.K += 24. * M_PI * matter_vars.rho; + current_cell.store_vars(metric_vars); +} + +#endif /* INITIALFLUIDDATA_IMPL_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/Main_PerfectFluid.cpp b/Examples/Fluid_Kerr_merged/Main_PerfectFluid.cpp new file mode 100644 index 000000000..374df002b --- /dev/null +++ b/Examples/Fluid_Kerr_merged/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/Fluid_Kerr_merged/PerfectFluid.hpp b/Examples/Fluid_Kerr_merged/PerfectFluid.hpp new file mode 100644 index 000000000..7157a4b8d --- /dev/null +++ b/Examples/Fluid_Kerr_merged/PerfectFluid.hpp @@ -0,0 +1,107 @@ +/* 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 "ConservedQuantities.hpp" +#include "DefaultEoS.hpp" +#include "DimensionDefinitions.hpp" +#include "Fluxes.hpp" +#include "FourthOrderDerivatives.hpp" +#include "Sources.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS, total num of components +#include "VarsTools.hpp" +#include "WENODerivatives.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. + It is templated over an equation of state function eos_t, which the + user must specify in a class, although a default is provided which + sets P_of_rho to rho * (1 + eps) / 3. +*/ +template class PerfectFluid +{ + protected: + //! The local copy of the potential + double m_dx; + double m_lambda; + eos_t my_eos; + + public: + //! Constructor of class PerfectFluid, inputs are the matter parameters. + PerfectFluid(double a_dx, const double a_lambda, const eos_t a_eos) + : m_dx(a_dx), m_lambda(a_lambda), 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 tau; + // data_t Jt; + Tensor<1, data_t> vi; + data_t rho; + data_t eps; + // data_t nn; + + /// 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_tau, tau); + // define_enum_mapping(mapping_function, c_Jt, Jt); + define_enum_mapping(mapping_function, GRInterval(), + vi); + define_enum_mapping(mapping_function, c_rho, rho); + define_enum_mapping(mapping_function, c_eps, eps); + // define_enum_mapping(mapping_function, c_nn, nn); + } + }; + + //! 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 first derivative + const vars_t> &lm, //!< value of the left_minus weno + const vars_t> &lp, //!< value of the left_plus weno + const vars_t> &rm, //!< value of the right_minus weno + const vars_t> &rp) //!< value of the right_plus weno + const; //!< the value of the advection terms +}; + +#include "PerfectFluid.impl.hpp" + +#endif /* PERFECTFLUID_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp b/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp new file mode 100644 index 000000000..8ef31ec0b --- /dev/null +++ b/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp @@ -0,0 +1,163 @@ +/* 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_rho = 0.; + my_eos.compute_eos(P_of_rho, vars); + + // Useful quantities + data_t chi_regularised = simd_max(1e-6, vars.chi); + data_t v2 = 0.; + FOR(i, j) + { + v2 += vars.h[i][j] * vars.vi[i] * vars.vi[j] / chi_regularised; + } + data_t WW = 1. / (1. - v2); + data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + + 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; } + } + + // Calculate components of EM Tensor + // S_ij = T_ij + FOR(i, j) + { + out.Sij[i][j] = vars.rho * hh * WW * vi_D[i] * vi_D[j] + + vars.h[i][j] * P_of_rho / chi_regularised; + } + + // 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.rho * hh * WW * vi_D[i]; } + + // rho = n^a n^b T_ab + out.rho = vars.rho * hh * WW - P_of_rho; + + return out; +} + +template +template class 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 vars_t> &lm, + const vars_t> &lp, const vars_t> &rm, + const vars_t> &rp) const +{ + using namespace TensorAlgebra; + + data_t P_of_rho = 0.0; + my_eos.compute_eos(P_of_rho, vars); + + vars_t source = Sources::compute_source(P_of_rho, vars, d1); + // evolution equations for the fluid conservative variables + data_t divshift = TensorAlgebra::compute_trace(d1.shift); + data_t chi_regularised = simd_max(vars.chi, 1e-6); + data_t advec_chi = 0.; + FOR(i) advec_chi += vars.shift[i] * d1.chi[i] / chi_regularised; + 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); + FOR(i) + { + vars_t flux = Fluxes::compute_flux(P_of_rho, 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]; + rhs.tau += GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.tau; + } + + FOR(idir) + { + vars_t vars_right_p = vars; + vars_right_p.rho = rp.rho[idir]; + vars_right_p.eps = rp.eps[idir]; + FOR(j) { vars_right_p.vi[j] = rp.vi[j][idir]; } + my_eos.compute_eos(P_of_rho, vars_right_p); + ConservedQuantities::PtoC(P_of_rho, vars_right_p); + vars_t flux_right_p = Fluxes::compute_num_flux( + P_of_rho, vars_right_p, idir, m_lambda, -1); + + vars_t vars_right_m = vars; + vars_right_m.rho = rm.rho[idir]; + vars_right_m.eps = rm.eps[idir]; + FOR(j) vars_right_m.vi[j] = rm.vi[j][idir]; + my_eos.compute_eos(P_of_rho, vars_right_m); + ConservedQuantities::PtoC(P_of_rho, vars_right_m); + vars_t flux_right_m = + Fluxes::compute_num_flux(P_of_rho, vars_right_m, idir, m_lambda, 1); + + rhs.D += -1. / (2. * m_dx) * (flux_right_p.D + flux_right_m.D); + FOR(j) + { + rhs.Sj[j] += + -1. / (2. * m_dx) * (flux_right_p.Sj[j] + flux_right_m.Sj[j]); + } + rhs.tau += -1. / (2. * m_dx) * (flux_right_p.tau + flux_right_m.tau); + } + + FOR(idir) + { + vars_t vars_left_p = vars; + vars_left_p.rho = lp.rho[idir]; + vars_left_p.eps = lp.eps[idir]; + FOR(j) { vars_left_p.vi[j] = lp.vi[j][idir]; } + my_eos.compute_eos(P_of_rho, vars_left_p); + ConservedQuantities::PtoC(P_of_rho, vars_left_p); + vars_t flux_left_p = + Fluxes::compute_num_flux(P_of_rho, vars_left_p, idir, m_lambda, -1); + + vars_t vars_left_m = vars; + vars_left_m.rho = lm.rho[idir]; + vars_left_m.eps = lm.eps[idir]; + FOR(j) { vars_left_m.vi[j] = lm.vi[j][idir]; } + my_eos.compute_eos(P_of_rho, vars_left_m); + ConservedQuantities::PtoC(P_of_rho, vars_left_m); + vars_t flux_left_m = + Fluxes::compute_num_flux(P_of_rho, vars_left_m, idir, m_lambda, 1); + + rhs.D += 1. / (2. * m_dx) * (flux_left_p.D + flux_left_m.D); + FOR(j) + { + rhs.Sj[j] += + 1. / (2. * m_dx) * (flux_left_p.Sj[j] + flux_left_m.Sj[j]); + } + rhs.tau += 1. / (2. * m_dx) * (flux_left_p.tau + flux_left_m.tau); + } + + rhs.rho = 0.; + rhs.eps = 0.; + FOR(i) { rhs.vi[i] = 0.; } +} + +#endif /* PERFECTFLUID_IMPL_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp new file mode 100644 index 000000000..4d2f177a4 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp @@ -0,0 +1,138 @@ +/* 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 "PrimitiveRecovery.hpp" +#include "SixthOrderDerivatives.hpp" +#include "TraceARemoval.hpp" +#include "WENODerivatives.hpp" + +// For RHS update +#include "FluidCCZ4RHS.hpp" + +// For constraints calculation +#include "NewMatterConstraints.hpp" + +// For tag cells +// #include "FixedGridsTaggingCriterion.hpp" +#include "ChiTaggingCriterion.hpp" + +// Problem specific includes +#include "ComputePack.hpp" +#include "EoS.hpp" +#include "GammaCalculator.hpp" +#include "InitialFluidData.hpp" +#include "KerrBH.hpp" +#include "PerfectFluid.hpp" +#include "PositiveDensity.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(), + PositiveDensity()), + m_state_new, m_state_new, INCLUDE_GHOST_CELLS, disable_simd()); + + // 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.), KerrBH(m_p.kerr_params, m_dx), + InitialFluidData(m_p.initial_params, m_dx)), + m_state_new, m_state_new, INCLUDE_GHOST_CELLS, disable_simd()); + + fillAllGhosts(); + BoxLoops::loop(GammaCalculator(m_dx), m_state_new, m_state_new, + EXCLUDE_GHOST_CELLS, disable_simd()); +} + +#ifdef CH_USE_HDF5 +// Things to do before outputting a checkpoint file +void PerfectFluidLevel::prePlotLevel() +{ + fillAllGhosts(); + EoS eos; + PerfectFluidEoS perfect_fluid(m_dx, m_p.lambda, eos); + BoxLoops::loop(MatterConstraints(perfect_fluid, m_dx, + m_p.G_Newton, c_Ham, + Interval(c_Mom1, c_Mom3)), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS, disable_simd()); +} +#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(), PositiveDensity(), + PositiveChiAndAlpha(), + PrimitiveRecovery()), + a_soln, a_soln, INCLUDE_GHOST_CELLS/*, disable_simd()*/); + + // Calculate MatterCCZ4 right hand side with matter_t = ScalarField + EoS eos; + PerfectFluidEoS perfect_fluid(m_dx, m_p.lambda, eos); + if (m_p.max_spatial_derivative_order == 4) + { + FluidCCZ4RHS + // FourthOrderDerivatives> + 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, disable_simd()); + } + else if (m_p.max_spatial_derivative_order == 6) + { + FluidCCZ4RHS + // FourthOrderDerivatives> + 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, disable_simd()); + } +} + +// 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(make_compute_pack(TraceARemoval(), PositiveDensity(), + PositiveChiAndAlpha(), + PrimitiveRecovery()), + a_soln, a_soln, INCLUDE_GHOST_CELLS/*, disable_simd()*/); +} + +void PerfectFluidLevel::preTagCells() +{ + // We only use chi in the tagging criterion so only fill the ghosts for chi + fillAllGhosts(VariableType::evolution, Interval(c_chi, c_chi)); +} + +void PerfectFluidLevel::computeTaggingCriterion( + FArrayBox &tagging_criterion, const FArrayBox ¤t_state, + const FArrayBox ¤t_state_diagnostics) +{ + BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion, disable_simd()); +} diff --git a/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp new file mode 100644 index 000000000..40cffd497 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp @@ -0,0 +1,61 @@ +/* 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() override; + + //! Initialize data for the field and metric variables + virtual void initialData() override; + +#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) override; + + //! Things to do in UpdateODE step, after soln + rhs update + virtual void specificUpdateODE(GRLevelData &a_soln, + const GRLevelData &a_rhs, Real a_dt) override; + + /// 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, + const FArrayBox ¤t_state_diagnostics) override; +}; + +#endif /* PERFECTFLUIDLEVEL_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/PositiveDensity.hpp b/Examples/Fluid_Kerr_merged/PositiveDensity.hpp new file mode 100644 index 000000000..46fc960d4 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/PositiveDensity.hpp @@ -0,0 +1,71 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +// This compute class enforces the dominant energy conditions +#ifndef POSITIVEDENSITY_HPP_ +#define POSITIVEDENSITY_HPP_ + +#include "Cell.hpp" +#include "UserVariables.hpp" +#include "simd.hpp" + +class PositiveDensity +{ + // Only variables needed are chi, D, tau, Sj, h + template struct Vars + { + data_t chi; + data_t D; + data_t tau; + Tensor<1, data_t> Sj; + Tensor<2, data_t> h; + + 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_chi, chi); + define_enum_mapping(mapping_function, c_D, D); + define_enum_mapping(mapping_function, c_tau, tau); + define_enum_mapping( + mapping_function, GRInterval(), Sj); + define_symmetric_enum_mapping( + mapping_function, GRInterval(), h); + } + }; + private: + const double m_min_D; + const double m_min_v; + + public: + //! Constructor for class + PositiveDensity(const double a_min_D = 1e-12, const double a_min_v = 0.) + : m_min_D(a_min_D), m_min_v(a_min_v) + { + } + + template void compute(Cell current_cell) const + { + auto vars = current_cell.template load_vars(); + const data_t chi_regularised = simd_max(1e-6, vars.chi); + + // Take into account that conservative variables are conformally rescaled + + //D >= min_D + vars.D = simd_max(vars.D, m_min_D / pow(chi_regularised, 1.5)); + + data_t S2_over_chi = 0.; + FOR(i,j) S2_over_chi += vars.h[i][j] * vars.Sj[i] * vars.Sj[j]; + + // S^2 <= (D + tau)^2 + data_t factor = simd_min(1., vars.chi * (vars.D + vars.tau) / sqrt(S2_over_chi)); + FOR(i) vars.Sj[i] *= factor; + + current_cell.store_vars(vars.D, c_D); + current_cell.store_vars(vars.Sj, GRInterval()); + } +}; + +#endif /* POSITIVEDENSITY_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp b/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp new file mode 100644 index 000000000..b521e5073 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp @@ -0,0 +1,109 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +// This class converts conservative to primitive variables +#ifndef PRIMITIVERECOVERY_HPP_ +#define PRIMITIVERECOVERY_HPP_ + +#include "CCZ4Geometry.hpp" +#include "Cell.hpp" +#include "DefaultEoS.hpp" +#include "EoS.hpp" +#include "simd.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" +#include "UsingNamespace.H" +#include "VarsTools.hpp" + +// template +class PrimitiveRecovery : public EoS +{ + public: + template struct Vars + { + Tensor<2, data_t> h; + data_t chi, D, tau, rho, eps; + Tensor<1, data_t> Sj, vi; + + template + void enum_mapping(mapping_function_t mapping_function); + }; + + template void compute(Cell current_cell) const + { + auto vars = current_cell.template load_vars(); + + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + + data_t tolerance = 1e-8; + data_t Wa2, Wa, xn, diff; + + data_t P_of_rho = 0.; + + data_t S2 = 0.; + FOR(i, j) S2 += vars.chi * h_UU[i][j] * vars.Sj[i] * vars.Sj[j]; + + data_t q = vars.tau / vars.D; + data_t r = S2 / pow(vars.D, 2.); + data_t xa = 1.5 * (1. + q); + Wa = sqrt(pow(xa, 2.) / (pow(xa, 2.) - r)); + + vars.rho = vars.D / Wa; + vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q); + EoS::compute_eos(P_of_rho, vars); + xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); + diff = abs(xn - xa); + + int i = 0; + + data_t empty; + while (!simd_all_false(simd_compare_lt(tolerance, diff), empty)) + { + i++; + Wa = sqrt(pow(xa, 2.) / (pow(xa, 2.) - r)); + + vars.rho = vars.D / Wa; + vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q); + EoS::compute_eos(P_of_rho, vars); + xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); + diff = abs(xn - xa); + xa = xn; + if (i >= 100) + break; + } + data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + Tensor<1, data_t> vi_D; + FOR(i) vi_D[i] = vars.Sj[i] / vars.rho / hh / Wa / Wa; + + FOR(i) + { + vars.vi[i] = 0.; + FOR(j) vars.vi[i] += vars.chi * h_UU[i][j] * vi_D[j]; + } + + current_cell.store_vars(vars); + } +}; + +template +template +void PrimitiveRecovery::Vars::enum_mapping( + mapping_function_t mapping_function) +{ + VarsTools::define_symmetric_enum_mapping(mapping_function, + GRInterval(), h); + VarsTools::define_enum_mapping(mapping_function, c_chi, chi); + VarsTools::define_enum_mapping(mapping_function, c_D, D); + VarsTools::define_enum_mapping(mapping_function, GRInterval(), + Sj); + VarsTools::define_enum_mapping(mapping_function, c_tau, tau); + VarsTools::define_enum_mapping(mapping_function, GRInterval(), + vi); + VarsTools::define_enum_mapping(mapping_function, c_rho, rho); + VarsTools::define_enum_mapping(mapping_function, c_eps, eps); +} + +#endif /* PRIMITIVERECOVERY_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/SimulationParameters.hpp b/Examples/Fluid_Kerr_merged/SimulationParameters.hpp new file mode 100644 index 000000000..9c591037e --- /dev/null +++ b/Examples/Fluid_Kerr_merged/SimulationParameters.hpp @@ -0,0 +1,88 @@ +/* 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 "EoS.hpp" +#include "InitialFluidData.hpp" +#include "KerrBH.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 fluid 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_rho0", initial_params.rho0, 1.0); + pp.load("fluid_width", initial_params.awidth, 0.1); + pp.load("fluid_delta", initial_params.delta, 0.2); + pp.load("lambda", lambda, 1.); // eigenvalue for numerical flux + + // Reading data + // pp.load("spacing", initial_params.spacing); + //pp.load("lines", lines); + //double rho_1D[lines]; + + //double tmp_data; + //ifstream read_file("rho_vs_r.csv"); + + //for (int i = 0; i < lines; ++i) + //{ + // read_file >> tmp_data; + // rho_1D[i] = tmp_data; + //} + //initial_params.rho_1D = rho_1D; + + // Initial Kerr data + pp.load("kerr_mass", kerr_params.mass); + pp.load("kerr_spin", kerr_params.spin); + pp.load("kerr_center", kerr_params.center, center); + } + + void check_params() + { + 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; + double lambda; + //int lines; + InitialFluidData::params_t initial_params; + KerrBH::params_t kerr_params; +}; + +#endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/Sources.hpp b/Examples/Fluid_Kerr_merged/Sources.hpp new file mode 100644 index 000000000..55c3e8a43 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/Sources.hpp @@ -0,0 +1,63 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SOURCES_HPP_ +#define SOURCES_HPP_ + +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "simd.hpp" + +namespace Sources +{ +// Primitive to conservative variables +template class vars_t> +vars_t compute_source(const data_t P_of_rho, const vars_t &vars, + const vars_t> &d1) +{ + data_t chi_regularised = simd_max(1e-6, vars.chi); + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + vars_t out; + + data_t v2 = 0.; + FOR(i, j) v2 += vars.h[i][j] * vars.vi[j] * vars.vi[i] / chi_regularised; + + data_t WW = 1. / (1. - v2); + data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + + out.D = 0.; + FOR(j) + { + out.Sj[j] = -(vars.tau + vars.D) * d1.lapse[j]; + FOR(i) + { + out.Sj[j] += vars.Sj[i] * d1.shift[i][j]; + FOR(k) + { + out.Sj[j] += vars.lapse / 2. * + (d1.h[i][k][j] - + vars.h[i][k] * d1.chi[j] / chi_regularised) * + (vars.rho * hh * WW * vars.vi[i] * vars.vi[k] / + chi_regularised + + P_of_rho * h_UU[i][k]); + } + } + } + out.tau = 0.; + FOR(i, j) + { + out.tau += vars.lapse * (vars.A[i][j] + vars.h[i][j] / 3. * vars.K) * + (vars.rho * hh * WW * vars.vi[i] * vars.vi[j] / + chi_regularised + + P_of_rho * h_UU[i][j]) - + vars.chi * h_UU[i][j] * vars.Sj[i] * d1.lapse[j]; + } + + return out; +} + +} // namespace Sources +#endif /* SOURCES_HPP_ */ diff --git a/Examples/Fluid_Kerr_merged/UserVariables.hpp b/Examples/Fluid_Kerr_merged/UserVariables.hpp new file mode 100644 index 000000000..ad9e00927 --- /dev/null +++ b/Examples/Fluid_Kerr_merged/UserVariables.hpp @@ -0,0 +1,47 @@ +/* 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_tau, + // c_Jt, + c_vi1, + c_vi2, + c_vi3, + c_rho, + c_eps, + // c_nn, + + NUM_VARS +}; + +namespace UserVariables +{ +static const std::array + user_variable_names = {"D", "Sj1", "Sj2", "Sj3", "tau", + "vi1", "vi2", "vi3", "rho", "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/Fluid_Kerr_merged/WENODerivatives.hpp b/Examples/Fluid_Kerr_merged/WENODerivatives.hpp new file mode 100644 index 000000000..ff958975c --- /dev/null +++ b/Examples/Fluid_Kerr_merged/WENODerivatives.hpp @@ -0,0 +1,165 @@ +/* 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: + enum + { + LEFT_MINUS, + LEFT_PLUS, + RIGHT_MINUS, + RIGHT_PLUS, + }; + + // public: + // WENODerivatives(double a_eW) : m_eW(a_eW){} + WENODerivatives() {} + template + ALWAYS_INLINE data_t get_Pface(const double *in_ptr, const int idx, + const int stride, int dir_switch) const + { + data_t beta[3] = {0., 0., 0.}; + const double dd[3] = {3. / 10., 3. / 5., 1. / 10.}; + const double eW = 1.; + data_t alpha[3] = {0., 0., 0.}; + data_t sum_alpha; + data_t weights[3] = {0., 0., 0.}; + data_t v[3] = {0., 0., 0.}; + // double u[3] = {0.,0.,0.}; + data_t pim2, pim1, pi0, pip1, pip2; + + auto in = SIMDIFY(in_ptr); + + if (dir_switch == 0) // LEFT_MINUS + // Negative fluxes to compute primitive variables + // at the left cell boundary i-1/2: p^{L-}_{i-1/2} + { + pim2 = in[idx - 3 * stride]; + pim1 = in[idx - 2 * stride]; + pi0 = in[idx - stride]; + pip1 = in[idx]; + pip2 = in[idx + stride]; + + // ENO polynomials + v[0] = (2. * pim2 - 7. * pim1 + 11. * pi0) / 6.; + v[1] = (-pim1 + 5. * pi0 + 2. * pip1) / 6.; + v[2] = (2. * pi0 + 5. * pip1 - pip2) / 6.; + } + + if (dir_switch == 1) // LEFT_PLUS + { + // Positive fluxes to compute primitive variables + // at the left cell boundary i-1/2: p^{L+}_{i-1/2} + pim2 = in[idx - 2 * stride]; + pim1 = in[idx - stride]; + pi0 = in[idx]; + pip1 = in[idx + stride]; + pip2 = in[idx + 2 * stride]; + + // ENO polynomials + v[0] = (-pim2 + 5. * pim1 + 2. * pi0) / 6.; + v[1] = (2. * pim1 + 5. * pi0 - pip1) / 6.; + v[2] = (11. * pi0 - 7. * pip1 + 2. * pip2) / 6.; + } + + if (dir_switch == 2) // RIGHT_MINUS + // Negative fluxes to compute primitive variables + // at the right cell boundary i+1/2: p^{R-}_{i+1/2} + { + pim2 = in[idx - 2 * stride]; + pim1 = in[idx - stride]; + pi0 = in[idx]; + pip1 = in[idx + stride]; + pip2 = in[idx + 2 * stride]; + + // ENO polynomials + v[0] = (2. * pim2 - 7. * pim1 + 11. * pi0) / 6.; + v[1] = (-pim1 + 5. * pi0 + 2. * pip1) / 6.; + v[2] = (2. * pi0 + 5. * pip1 - pip2) / 6.; + } + + if (dir_switch == 3) // RIGHT_PLUS + { + // Positive fluxes to compute primitive variables + // at the right cell boundary i+1/2: p^{R+}_{i+1/2} + pim2 = in[idx - stride]; + pim1 = in[idx]; + pi0 = in[idx + stride]; + pip1 = in[idx + 2 * stride]; + pip2 = in[idx + 3 * stride]; + + // ENO polynomials + v[0] = (-pim2 + 5. * pim1 + 2. * pi0) / 6.; + v[1] = (2. * pim1 + 5. * pi0 - pip1) / 6.; + v[2] = (11. * pi0 - 7. * pip1 + 2. * pip2) / 6.; + } + + // Smoothness indicators + beta[0] = (13. / 12.) * pow(pim2 - 2. * pim1 + pi0, 2) + + (1. / 4.) * pow(pim2 - 4. * pim1 + 3. * pi0, 2); + beta[1] = (13. / 12.) * pow(pim1 - 2. * pi0 + pip1, 2) + + (1. / 4.) * pow(pim1 - pip1, 2); + beta[2] = (13. / 12.) * pow(pi0 - 2. * pip1 + pip2, 2) + + (1. / 4.) * pow(3. * pi0 - 4. * pip1 + pip2, 2); + + // Weights + sum_alpha = 0.; + FOR(j) + { + alpha[j] = dd[j] / ((eW + beta[j]) * (eW + beta[j])); + sum_alpha += alpha[j]; + } + FOR(j) { weights[j] = alpha[j] / sum_alpha; } + + // primitive variables on the cell interface + data_t pFace = 0.; + FOR(j) { pFace += weights[j] * v[j]; } + return pFace; + } + + template