From 0330e51ffc3f106f3234b555269b3bd37b3921e0 Mon Sep 17 00:00:00 2001 From: dinatraykova Date: Fri, 28 Jul 2023 17:29:56 +0200 Subject: [PATCH] Add fluid matter rhs class --- Examples/PerfectFluid/FluidCCZ4RHS.hpp | 118 +++++++++++++++++++ Examples/PerfectFluid/FluidCCZ4RHS.impl.hpp | 120 ++++++++++++++++++++ Examples/PerfectFluid/InitialFluidData.hpp | 7 +- Examples/PerfectFluid/PerfectFluid.hpp | 7 +- Examples/PerfectFluid/PerfectFluid.impl.hpp | 7 +- Examples/PerfectFluid/PerfectFluidLevel.cpp | 14 ++- Examples/PerfectFluid/WENODerivatives.hpp | 65 +++++++---- 7 files changed, 300 insertions(+), 38 deletions(-) create mode 100644 Examples/PerfectFluid/FluidCCZ4RHS.hpp create mode 100644 Examples/PerfectFluid/FluidCCZ4RHS.impl.hpp diff --git a/Examples/PerfectFluid/FluidCCZ4RHS.hpp b/Examples/PerfectFluid/FluidCCZ4RHS.hpp new file mode 100644 index 000000000..365d65381 --- /dev/null +++ b/Examples/PerfectFluid/FluidCCZ4RHS.hpp @@ -0,0 +1,118 @@ +/* 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 "WENODerivatives.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 "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 + 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/PerfectFluid/FluidCCZ4RHS.impl.hpp b/Examples/PerfectFluid/FluidCCZ4RHS.impl.hpp new file mode 100644 index 000000000..3b961feb2 --- /dev/null +++ b/Examples/PerfectFluid/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, d1, d1, d1); + + // 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/PerfectFluid/InitialFluidData.hpp b/Examples/PerfectFluid/InitialFluidData.hpp index d1943ee86..c9092b518 100644 --- a/Examples/PerfectFluid/InitialFluidData.hpp +++ b/Examples/PerfectFluid/InitialFluidData.hpp @@ -9,7 +9,7 @@ #include "ADMConformalVars.hpp" #include "Cell.hpp" #include "Coordinates.hpp" -#include "MatterCCZ4RHS.hpp" +#include "FluidCCZ4RHS.hpp" #include "PerfectFluid.hpp" #include "Tensor.hpp" #include "UserVariables.hpp" //This files needs NUM_VARS - total no. components @@ -53,9 +53,8 @@ class InitialFluidData // where am i? Coordinates coords(current_cell, m_dx, m_params.center); - // MetricVars metric_vars; - //m_background.compute_metric_background(metric_vars, coords); - + MetricVars metric_vars; + data_t x = coords.x; double y = coords.y; double z = coords.z; diff --git a/Examples/PerfectFluid/PerfectFluid.hpp b/Examples/PerfectFluid/PerfectFluid.hpp index 1d111abc4..4fd21779c 100644 --- a/Examples/PerfectFluid/PerfectFluid.hpp +++ b/Examples/PerfectFluid/PerfectFluid.hpp @@ -102,9 +102,10 @@ template class PerfectFluid 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 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 }; diff --git a/Examples/PerfectFluid/PerfectFluid.impl.hpp b/Examples/PerfectFluid/PerfectFluid.impl.hpp index 7eb0ee981..208111b5f 100644 --- a/Examples/PerfectFluid/PerfectFluid.impl.hpp +++ b/Examples/PerfectFluid/PerfectFluid.impl.hpp @@ -70,9 +70,10 @@ 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 diff2_vars_t> &d2, - const vars_t &advec) const + const vars_t> &lm, + const vars_t> &lp, + const vars_t> &rm, + const vars_t> &rp) const { using namespace TensorAlgebra; diff --git a/Examples/PerfectFluid/PerfectFluidLevel.cpp b/Examples/PerfectFluid/PerfectFluidLevel.cpp index 571f08061..ae8b773ca 100644 --- a/Examples/PerfectFluid/PerfectFluidLevel.cpp +++ b/Examples/PerfectFluid/PerfectFluidLevel.cpp @@ -13,7 +13,7 @@ #include "TraceARemoval.hpp" // For RHS update -#include "MatterCCZ4RHS.hpp" +#include "FluidCCZ4RHS.hpp" // For constraints calculation #include "NewMatterConstraints.hpp" @@ -54,6 +54,8 @@ void PerfectFluidLevel::initialData() // First set everything to zero then initial conditions for scalar field - // here a Kerr BH and a scalar field profile + // EoS eos(m_p.eos_params); + //PerfectFluidEoS perfect_fluid(eos); BoxLoops::loop( make_compute_pack(SetValue(0.), Minkowski(m_p.kerr_params, m_dx), InitialFluidData(m_p.initial_params, m_dx)), @@ -93,16 +95,18 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, PerfectFluidEoS perfect_fluid(eos); if (m_p.max_spatial_derivative_order == 4) { - MatterCCZ4RHS + 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); } else if (m_p.max_spatial_derivative_order == 6) { - MatterCCZ4RHS + 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); diff --git a/Examples/PerfectFluid/WENODerivatives.hpp b/Examples/PerfectFluid/WENODerivatives.hpp index 91b30fecf..0f266f886 100644 --- a/Examples/PerfectFluid/WENODerivatives.hpp +++ b/Examples/PerfectFluid/WENODerivatives.hpp @@ -14,8 +14,8 @@ class WENODerivatives { - public: - enum +public: + enum { LEFT_MINUS, LEFT_PLUS, @@ -23,22 +23,22 @@ class WENODerivatives RIGHT_PLUS, }; - // public: - WENODerivatives(double a_eW): m_eW(a_eW){} + //public: + WENODerivatives(double a_eW) : m_eW(a_eW){} template ALWAYS_INLINE data_t get_Pface(const double *in_ptr, const int idx, const int stride, int dir_switch) const { - double beta[3] = {0.,0.,0.}; + data_t beta[3] = {0.,0.,0.}; const double dd[3] = {3./10.,3./5.,1./10.}; //const double eW = 1.; - double alpha[3] = {0.,0.,0.}; - double sum_alpha; - double weights[3] = {0.,0.,0.}; - double v[3] = {0.,0.,0.}; - double u[3] = {0.,0.,0.}; - double pim2, pim1, pi0, pip1, pip2; + 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); @@ -63,11 +63,11 @@ class WENODerivatives { // Positive fluxes to compute primitive variables // at the left cell boundary i-1/2: p^{L+}_{i-1/2} - pim2 = [idx - 2.*stride]; - pim1 = [idx - stride]; - pi0 = [idx]; - pip1 = [idx + stride]; - pip2 = [idx + 2.*stride]; + 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.; @@ -96,11 +96,11 @@ class WENODerivatives { // Positive fluxes to compute primitive variables // at the right cell boundary i+1/2: p^{R+}_{i+1/2} - pim2 = [idx - stride]; - pim1 = [idx]; - pi0 = [idx + stride]; - pip1 = [idx + 2.*stride]; - pip2 = [idx + 3.*stride]; + 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.; @@ -120,7 +120,7 @@ class WENODerivatives sum_alpha = 0.; FOR(j) { - alpha[j] = dd[j]/((a_eW+beta[j])*(a_eW+beta[j])); + alpha[j] = dd[j]/((m_eW+beta[j])*(m_eW+beta[j])); sum_alpha += alpha[j]; } FOR(j) @@ -129,11 +129,30 @@ class WENODerivatives } // primitive variables on the cell interface - double pFace = 0.; + data_t pFace = 0.; FOR(j) {pFace += weights[j]*v[j];} return pFace; } + template