Add boosted bh proca example

Examples/BoostedBHProca/BoostedBHProcaLevel.cpp

@@ -0,0 +1,164 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+// General includes common to most GR problems
+#include "BoostedBHProcaLevel.hpp"
+#include "AMRReductions.hpp"
+#include "BoxLoops.hpp"
+#include "ComputePack.hpp"
+#include "NanCheck.hpp"
+#include "SetValue.hpp"
+#include "SmallDataIO.hpp"
+
+// For RHS update
+#include "BoostedBH.hpp"
+#include "MatterEvolution.hpp"
+
+// For tag cells
+#include "FixedGridsTaggingCriterion.hpp"
+
+// Problem specific includes
+#include "EnergyConservation.hpp"
+#include "ExcisionDiagnostics.hpp"
+#include "ExcisionEvolution.hpp"
+#include "FluxExtraction.hpp"
+#include "InitialProcaData.hpp"
+#include "LinearMomConservation.hpp"
+#include "ProcaConstraint.hpp"
+#include "ProcaField.hpp"
+
+// Initial data for field and metric variables
+void BoostedBHProcaLevel::initialData()
+{
+    CH_TIME("BoostedBHProcaLevel::initialData");
+    if (m_verbosity)
+        pout() << "BoostedBHProcaLevel::initialData " << m_level << endl;
+
+    // First set everything to zero, then set the value of the conformal factor
+    // This is just for the diagnostics
+    SetValue set_zero(0.0);
+    BoostedBH boosted_bh(m_p.bg_params, m_dx); // just calculates chi
+    auto compute_pack = make_compute_pack(set_zero, boosted_bh);
+    BoxLoops::loop(compute_pack, m_state_diagnostics, m_state_diagnostics,
+                   SKIP_GHOST_CELLS);
+
+    // Now set the actual evolution variables
+    InitialProcaData initial_pf(m_p.proca_amplitude, m_p.proca_mass, m_p.center,
+                                m_p.bg_params, m_dx);
+    BoxLoops::loop(initial_pf, m_state_new, m_state_new, FILL_GHOST_CELLS);
+
+    // now the gauss constraint
+    fillAllGhosts();
+    ProcaConstraint enforce_constraint(m_p.center, m_p.bg_params,
+                                       m_p.proca_mass, m_dx);
+    BoxLoops::loop(enforce_constraint, m_state_new, m_state_new,
+                   EXCLUDE_GHOST_CELLS);
+
+    // excise evolution vars within horizon, turn off simd vectorisation
+    BoxLoops::loop(
+        ExcisionEvolution<ProcaField, BoostedBH>(m_dx, m_p.center, boosted_bh),
+        m_state_new, m_state_new, SKIP_GHOST_CELLS, disable_simd());
+}
+
+void BoostedBHProcaLevel::specificPostTimeStep()
+{
+    // Check for nans on every level
+    if (m_p.nan_check)
+        BoxLoops::loop(NanCheck(), m_state_new, m_state_new, SKIP_GHOST_CELLS,
+                       disable_simd());
+
+    // At any level, but after the timestep on the minimum extraction level
+    int min_level = m_p.extraction_params.min_extraction_level();
+    bool calculate_diagnostics = at_level_timestep_multiple(min_level);
+    if (calculate_diagnostics)
+    {
+        fillAllGhosts();
+        ProcaField proca_field(m_p.proca_mass, m_p.proca_damping);
+        BoostedBH boosted_bh(m_p.bg_params, m_dx);
+        EnergyConservation<ProcaField, BoostedBH> energies(
+            proca_field, boosted_bh, m_dx, m_p.center);
+        int direction = 0; // we want the x direction for the momentum
+        LinearMomConservation<ProcaField, BoostedBH> linear_momenta(
+            proca_field, boosted_bh, direction, m_dx, m_p.center);
+        BoxLoops::loop(make_compute_pack(energies, linear_momenta), m_state_new,
+                       m_state_diagnostics, SKIP_GHOST_CELLS);
+
+        // excise within/outside specified radii, no simd
+        BoxLoops::loop(
+            ExcisionDiagnostics<ProcaField, BoostedBH>(
+                m_dx, m_p.center, boosted_bh, m_p.inner_r, m_p.outer_r),
+            m_state_diagnostics, m_state_diagnostics, SKIP_GHOST_CELLS,
+            disable_simd());
+    }
+
+    // write out the integral after each timestep on the min_level
+    if (m_p.activate_extraction == 1)
+    {
+        if (m_level == min_level)
+        {
+            bool first_step = (m_time == m_dt);
+            // integrate the densities and write to a file
+            AMRReductions<VariableType::diagnostic> amr_reductions(m_gr_amr);
+            double rhoEnergy_sum = amr_reductions.sum(c_rhoEnergy);
+            double rhoLinMom_sum = amr_reductions.sum(c_rhoLinMom);
+            double sourceLinMom_sum = amr_reductions.sum(c_sourceLinMom);
+
+            SmallDataIO integral_file(m_p.data_path + "EnergyIntegrals", m_dt,
+                                      m_time, m_restart_time,
+                                      SmallDataIO::APPEND, first_step);
+            // remove any duplicate data if this is post restart
+            integral_file.remove_duplicate_time_data();
+
+            std::vector<double> data_for_writing = {
+                rhoEnergy_sum, rhoLinMom_sum, sourceLinMom_sum};
+
+            // write data
+            if (first_step)
+            {
+                integral_file.write_header_line({"Energy density.",
+                                                 "Lin. Mom. density",
+                                                 "Lin. Mom. source"});
+            }
+            integral_file.write_time_data_line(data_for_writing);
+
+            // Now refresh the interpolator and do the interpolation
+            // only fill the actual ghost cells needed to save time
+            bool fill_ghosts = false;
+            m_gr_amr.m_interpolator->refresh(fill_ghosts);
+            m_gr_amr.fill_multilevel_ghosts(
+                VariableType::diagnostic, Interval(c_fluxEnergy, c_fluxLinMom));
+            FluxExtraction my_extraction(m_p.extraction_params, m_dt, m_time,
+                                         m_restart_time);
+            my_extraction.execute_query(m_gr_amr.m_interpolator, m_p.data_path);
+        }
+    }
+}
+
+// Things to do in RHS update, at each RK4 step
+void BoostedBHProcaLevel::specificEvalRHS(GRLevelData &a_soln,
+                                          GRLevelData &a_rhs,
+                                          const double a_time)
+{
+    // Calculate right hand side with matter_t = ProcaField
+    // and background_t = BoostedBH
+    ProcaField proca_field(m_p.proca_mass, m_p.proca_damping);
+    BoostedBH boosted_bh(m_p.bg_params, m_dx);
+    MatterEvolution<ProcaField, BoostedBH> my_evolution(
+        proca_field, boosted_bh, m_p.sigma, m_dx, m_p.center);
+    BoxLoops::loop(my_evolution, a_soln, a_rhs, SKIP_GHOST_CELLS);
+
+    // Do excision within horizon, don't use vectorisation
+    BoxLoops::loop(
+        ExcisionEvolution<ProcaField, BoostedBH>(m_dx, m_p.center, boosted_bh),
+        a_soln, a_rhs, SKIP_GHOST_CELLS, disable_simd());
+}
+
+// Note that for the fixed grids this only happens on the initial timestep
+void BoostedBHProcaLevel::computeTaggingCriterion(
+    FArrayBox &tagging_criterion, const FArrayBox &current_state)
+{
+    BoxLoops::loop(FixedGridsTaggingCriterion(m_dx, m_level, m_p.L, m_p.center),
+                   current_state, tagging_criterion);
+}

Examples/BoostedBHProca/BoostedBHProcaLevel.hpp

@@ -0,0 +1,43 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef BOOSTEDBHPROCALEVEL_HPP_
+#define BOOSTEDBHPROCALEVEL_HPP_
+
+#include "DefaultLevelFactory.hpp"
+#include "GRAMRLevel.hpp"
+// Problem specific includes
+#include "ProcaField.hpp"
+
+//!  A class for the evolution of a proca field, minimally coupled to gravity
+/*!
+     The class takes some initial data for a proca field (variables phi and Pi)
+     and evolves it on a fixed metric background.
+*/
+class BoostedBHProcaLevel : public GRAMRLevel
+{
+    friend class DefaultLevelFactory<BoostedBHProcaLevel>;
+    // Inherit the contructors from GRAMRLevel
+    using GRAMRLevel::GRAMRLevel;
+
+    // Typedef for proca field
+    typedef ProcaField ProcaField;
+
+    //! Initialize data for the field and metric variables
+    virtual void initialData();
+
+    //! RHS routines used at each RK4 step
+    virtual void specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
+                                 const double a_time);
+
+    // to do post each time step on every level
+    virtual void specificPostTimeStep();
+
+    //! Tell Chombo how to tag cells for regridding
+    virtual void computeTaggingCriterion(FArrayBox &tagging_criterion,
+                                         const FArrayBox &current_state);
+};
+
+#endif /* BOOSTEDBHPROCALEVEL_HPP_ */

Examples/BoostedBHProca/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
+enum
+{
+    c_chi,
+    c_rhoEnergy,
+    c_fluxEnergy,
+    c_rhoLinMom,
+    c_fluxLinMom,
+    c_sourceLinMom,
+
+    NUM_DIAGNOSTIC_VARS
+};
+
+namespace DiagnosticVariables
+{
+static const std::array<std::string, NUM_DIAGNOSTIC_VARS> variable_names = {
+    "chi",       "rhoEnergy",  "fluxEnergy",
+    "rhoLinMom", "fluxLinMom", "sourceLinMom"};
+}
+
+#endif /* DIAGNOSTICVARIABLES_HPP */

Examples/BoostedBHProca/ExcisionDiagnostics.hpp

@@ -0,0 +1,50 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef EXCISIONDIAGNOSTICS_HPP_
+#define EXCISIONDIAGNOSTICS_HPP_
+
+#include "Cell.hpp"
+#include "Coordinates.hpp"
+#include "GRInterval.hpp"
+#include "Tensor.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components
+#include "VarsTools.hpp"
+
+//! Does excision for the diagnostic variables - setting them to zero
+//! outside of the integration volume
+template <class matter_t, class background_t> class ExcisionDiagnostics
+{
+  protected:
+    const double m_dx;                              //!< The grid spacing
+    const std::array<double, CH_SPACEDIM> m_center; //!< The BH center
+    const FourthOrderDerivatives m_deriv;
+    const background_t m_background;
+    const double m_inner_r;
+    const double m_outer_r;
+
+  public:
+    ExcisionDiagnostics(const double a_dx,
+                        const std::array<double, CH_SPACEDIM> a_center,
+                        background_t a_background, const double a_inner_r,
+                        const double a_outer_r)
+        : m_dx(a_dx), m_deriv(m_dx), m_center(a_center),
+          m_background(a_background), m_inner_r(a_inner_r), m_outer_r(a_outer_r)
+    {
+    }
+
+    void compute(const Cell<double> current_cell) const
+    {
+        const Coordinates<double> coords(current_cell, m_dx, m_center);
+        if (coords.get_radius() < m_inner_r || coords.get_radius() > m_outer_r)
+        {
+            current_cell.store_vars(0.0, c_rhoLinMom);
+            current_cell.store_vars(0.0, c_rhoEnergy);
+            current_cell.store_vars(0.0, c_sourceLinMom);
+        }
+    }
+};
+
+#endif /* EXCISIONDIAGNOSTICS_HPP_ */

Examples/BoostedBHProca/FluxExtraction.hpp

@@ -0,0 +1,80 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef FLUXEXTRACTION_HPP_
+#define FLUXEXTRACTION_HPP_
+
+#include "SphericalExtraction.hpp"
+//!  The class allows extraction of the values of the force components on
+//!  spheroidal shells at specified radii, and integration over those shells
+/*!
+   The class allows the user to extract data from the grid for the force
+   components over spheroidal shells at specified radii. The values may then be
+   written to an output file, or integrated across the surfaces.
+*/
+class FluxExtraction : public SphericalExtraction
+{
+  public:
+    //! The constructor
+    FluxExtraction(const spherical_extraction_params_t &a_params, double a_dt,
+                   double a_time, bool a_first_step,
+                   double a_restart_time = 0.0)
+        : SphericalExtraction(a_params, a_dt, a_time, a_first_step,
+                              a_restart_time)
+    {
+        add_var(c_fluxEnergy, VariableType::diagnostic);
+        add_var(c_fluxLinMom, VariableType::diagnostic);
+    }
+
+    //! The old constructor which assumes it is called in specificPostTimeStep
+    //! so the first time step is when m_time == m_dt
+    FluxExtraction(const spherical_extraction_params_t &a_params, double a_dt,
+                   double a_time, double a_restart_time = 0.0)
+        : FluxExtraction(a_params, a_dt, a_time, (a_dt == a_time),
+                         a_restart_time)
+    {
+    }
+
+    // the references of the vars as used in the integrator
+    enum M_VARS
+    {
+        m_fluxEnergy,
+        m_fluxLinMom,
+        NUM_EXTRACTION_COMPS
+    };
+
+    //! Execute the query
+    void execute_query(AMRInterpolator<Lagrange<4>> *a_interpolator,
+                       std::string a_datapath)
+    {
+        // extract the values of the Flux scalars on the spheres
+        extract(a_interpolator);
+
+        // this would write out the values at every point on the sphere
+        if (m_params.write_extraction)
+        {
+            write_extraction(a_datapath + "FluxExtractionOut_");
+        }
+
+        // Setup to integrate fluxes
+        std::vector<std::vector<double>> force_integrals(NUM_EXTRACTION_COMPS);
+        add_var_integrand(m_fluxEnergy, force_integrals[m_fluxEnergy],
+                          IntegrationMethod::simpson);
+        add_var_integrand(m_fluxLinMom, force_integrals[m_fluxLinMom],
+                          IntegrationMethod::simpson);
+
+        // do the integration over the surface
+        integrate();
+
+        // write the integrals
+        std::vector<std::string> labels(NUM_EXTRACTION_COMPS);
+        labels[m_fluxEnergy] = "Energy Flux";
+        labels[m_fluxLinMom] = "Lin. Mom. Flux";
+        std::string filename = a_datapath + "FluxIntegrals";
+        write_integrals(filename, force_integrals, labels);
+    }
+};
+
+#endif /* FLUXEXTRACTION_HPP_ */

Examples/BoostedBHProca/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 := Set locally (e.g. export GRCHOMBO_SOURCE=path/to/GRChombo/Source in bash)
+
+GRDZHADZHA_SOURCE = ../../Source
+
+ebase := Main_BoostedBHProca
+
+LibNames := AMRTimeDependent AMRTools BoxTools
+
+src_dirs := $(GRCHOMBO_SOURCE)/utils \
+            $(GRCHOMBO_SOURCE)/simd  \
+            $(GRCHOMBO_SOURCE)/BoxUtils  \
+            $(GRCHOMBO_SOURCE)/TaggingCriteria  \
+            $(GRCHOMBO_SOURCE)/GRChomboCore  \
+            $(GRCHOMBO_SOURCE)/AMRInterpolator  \
+            $(GRDZHADZHA_SOURCE)/Background  \
+            $(GRDZHADZHA_SOURCE)/Matter
+
+include $(CHOMBO_HOME)/mk/Make.test