Add CosmoHamTagging

Examples/Cosmo/CosmoLevel.cpp

@@ -18,7 +18,7 @@
 #include "NewMatterConstraints.hpp"
 
 // For tag cells
-#include "FixedGridsTaggingCriterion.hpp"
+#include "CosmoHamTaggingCriterion.hpp"
 
 // Problem specific includes
 #include "AMRReductions.hpp"
@@ -78,6 +78,25 @@ void CosmoLevel::initialData()
     InitialK<ScalarFieldWithPotential> my_initial_K(scalar_field, m_dx,
                                                     m_p.G_Newton);
     BoxLoops::loop(my_initial_K, m_state_new, m_state_new, EXCLUDE_GHOST_CELLS);
+    // Calculate constraints and some diagnostics as we need it in tagging
+    // criterion
+    BoxLoops::loop(MatterConstraints<ScalarFieldWithPotential>(
+                       scalar_field, m_dx, m_p.G_Newton, c_Ham,
+                       Interval(c_Mom, c_Mom), c_Ham_abs_sum,
+                       Interval(c_Mom_abs_sum, c_Mom_abs_sum)),
+                   m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS);
+    CosmoDiagnostics<ScalarFieldWithPotential> cosmo_diagnostics(
+        scalar_field, m_dx, m_p.G_Newton);
+    BoxLoops::loop(cosmo_diagnostics, m_state_new, m_state_diagnostics,
+                   EXCLUDE_GHOST_CELLS);
+    // Assign initial rho_mean here
+    // from rho = 1/2 m^2 phi^2;
+    // phi0 = A sin(2 pi n x /L);
+    // mean of phi0^2 = A^2 sin^2 = 0.5*A^2;
+    // m = m_mode
+    m_cosmo_amr.set_rho_mean(
+        0.5 * m_p.scalar_field_mode * m_p.scalar_field_mode *
+        (0.5 * m_p.initial_params.amplitude * m_p.initial_params.amplitude));
 }
 
 #ifdef CH_USE_HDF5
@@ -141,17 +160,28 @@ void CosmoLevel::specificUpdateODE(GRLevelData &a_soln,
 
 void CosmoLevel::preTagCells()
 {
-    // we don't need any ghosts filled for the fixed grids tagging criterion
-    // used here so don't fill any
+    // Pre tagging - fill ghost cells and calculate Ham terms
+    fillAllGhosts();
+    Potential potential(m_p.potential_params, m_p.L, m_p.scalar_field_mode);
+    ScalarFieldWithPotential scalar_field(potential);
+    BoxLoops::loop(MatterConstraints<ScalarFieldWithPotential>(
+                       scalar_field, m_dx, m_p.G_Newton, c_Ham,
+                       Interval(c_Mom, c_Mom), c_Ham_abs_sum,
+                       Interval(c_Mom_abs_sum, c_Mom_abs_sum)),
+                   m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS);
+    CosmoDiagnostics<ScalarFieldWithPotential> cosmo_diagnostics(
+        scalar_field, m_dx, m_p.G_Newton);
+    BoxLoops::loop(cosmo_diagnostics, m_state_new, m_state_diagnostics,
+                   EXCLUDE_GHOST_CELLS);
 }
 
 void CosmoLevel::computeTaggingCriterion(
     FArrayBox &tagging_criterion, const FArrayBox &current_state,
     const FArrayBox &current_state_diagnostics)
 {
-    BoxLoops::loop(
-        FixedGridsTaggingCriterion(m_dx, m_level, 2.0 * m_p.L, m_p.center),
-        current_state, tagging_criterion);
+    BoxLoops::loop(CosmoHamTaggingCriterion(m_dx, m_p.center_tag, m_p.rad,
+                                            m_cosmo_amr.get_rho_mean()),
+                   current_state_diagnostics, tagging_criterion);
 }
 void CosmoLevel::specificPostTimeStep()
 {

Examples/Cosmo/SimulationParameters.hpp

@@ -36,8 +36,13 @@ class SimulationParameters : public SimulationParametersBase
         pp.load("scalar_mass", potential_params.scalar_mass, 0.1);
         pp.load("scalar_field_mode", scalar_field_mode, 1.0);
 
+        // Lineout params
         pp.load("lineout_num_points", lineout_num_points, 10);
 
+        // Tagging params
+        pp.load("center_tag", center_tag, center);
+        pp.load("rad", rad, L);
+
 #ifdef USE_AHFINDER
         double AH_guess =
             8. * initial_params.amplitude * initial_params.amplitude;
@@ -59,8 +64,9 @@ class SimulationParameters : public SimulationParametersBase
 
     // Initial data for matter and potential and BH
     double G_Newton;
-    double scalar_field_mode;
+    double scalar_field_mode, rad;
     int lineout_num_points;
+    std::array<double, CH_SPACEDIM> center_tag;
     InitialScalarData::params_t initial_params;
     Potential::params_t potential_params;
     KerrBH::params_t kerr_params;

Examples/Cosmo/params_cheap.txt

@@ -52,6 +52,10 @@ scalar_field_mode = 1.0
 # lineout params
 lineout_num_points = 10
 
+# Tagging
+# center_tag = 0 0 0
+# rad = 7
+
 #################################################
 # Grid parameters
 
@@ -66,15 +70,15 @@ N_full = 16
 L_full = 1
 
 # Maximum number of times you can regrid above coarsest level
-max_level = 0 # There are (max_level+1) grids, so min is zero
+max_level = 1 # There are (max_level+1) grids, so min is zero
 
 # Frequency of regridding at each level and thresholds on the tagging
 # Need one for each level except the top one, ie max_level items
 # Generally you do not need to regrid frequently on every level
 # in this example turn off regridding on all levels
 # Level Regridding: 0   1   2   3   4   5
-regrid_interval   = 0   0   0   0   0   0
-# regrid_threshold  = 0.5
+regrid_interval   = 10   0   0   0   0   0
+regrid_threshold  = 500
 
 regrid_interval = 0 0 0 0 0 0
 

Examples/Cosmo/plot.py

@@ -12,8 +12,15 @@
 plt.rcParams.update({'figure.figsize'    :  '6, 4.2'})
 plt.rcParams.update({'figure.autolayout': True})
 
-data = np.loadtxt('cheap_ex_data/data_out.dat')
-rho_dat = np.loadtxt('cheap_ex_data/lineout.dat')
+# Read header
+with open('cheap_data/lineout.dat', 'r') as file:
+    header = file.readline().strip()
+coord_interp = header.split()
+coord_interp = coord_interp[2:]
+
+# Read data
+data = np.loadtxt('cheap_data/data_out.dat')
+rho_dat = np.loadtxt('cheap_data/lineout.dat')
 t_data = data[:,0]
 chi_mean = data[:,3]
 rho_mean = data[:,4]
@@ -47,11 +54,13 @@
 dt_multiplier = 0.25
 dt = dx*dt_multiplier # dt = dx * dt_multiplier = 0.5 * 0.25
 t = 1/dt
-num_t_step = 25 # lineout every time = num_t_step
+num_t_step = 50 # lineout every time = num_t_step
 t_plot = np.arange(0,len(rho_interp),t*num_t_step)
 
+x_tick = np.arange(len(coord_interp))
+
 for it_plot, t_plot in enumerate(t_plot):
-    plt.plot(rho_interp[int(t_plot)],color=cm.Reds(8./(1+it_plot*10.),1.),
+    plt.plot(x_tick, rho_interp[int(t_plot)],color=cm.Reds(8./(1+it_plot*10.),1.),
          label='t = '+ str(t_plot*dt), marker = "o")
 
 plt.title('lineout of ' + r'$\rho$' + ' along x-axis')
@@ -62,5 +71,6 @@
 plt.ylabel(r'$\rho$', fontsize=14)
 #plt.yscale('log')
 #plt.ylim(0.,1.3e-3)
+plt.xticks(ticks=x_tick ,labels=coord_interp, rotation='vertical', fontsize=10)
 plt.savefig('plot_lineout.pdf',dpi=256, bbox_inches='tight',pad_inches = 0.1)
 plt.close()

Source/TaggingCriteria/CosmoHamTaggingCriterion.hpp

@@ -0,0 +1,52 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef COSMOHAMTAGGINGCRITERION_HPP_
+#define COSMOHAMTAGGINGCRITERION_HPP_
+
+#include "Cell.hpp"
+#include "CosmoAMR.hpp"
+#include "DimensionDefinitions.hpp"
+#include "FourthOrderDerivatives.hpp"
+#include "Tensor.hpp"
+
+class CosmoHamTaggingCriterion
+{
+  protected:
+    const double m_dx;
+    std::array<double, CH_SPACEDIM> m_center_tag;
+    double m_rad;
+    double m_rho_mean;
+
+  public:
+    CosmoHamTaggingCriterion(double dx,
+                             std::array<double, CH_SPACEDIM> center_tag,
+                             double rad, double rho_mean)
+        : m_dx(dx), m_center_tag(center_tag), m_rad(rad),
+          m_rho_mean(rho_mean){};
+
+    template <class data_t> void compute(Cell<data_t> current_cell) const
+    {
+        auto Ham_abs_sum = current_cell.load_vars(c_Ham_abs_sum);
+        auto sqrt_gamma = current_cell.load_vars(c_sqrt_gamma);
+
+        const Coordinates<data_t> coords(current_cell, m_dx, m_center_tag);
+        const data_t r = coords.get_radius();
+
+        // Divide Ham_abs_sum by rho_mean
+        data_t criterion = Ham_abs_sum / m_rho_mean * sqrt_gamma * m_dx;
+        auto regrid = simd_compare_gt(r, m_rad);
+
+        // data_t criterion = 0.0;
+        // auto regrid = simd_compare_lt(r, m_rad);
+
+        criterion = simd_conditional(regrid, 0.0, criterion);
+
+        // Write back into the flattened Chombo box
+        current_cell.store_vars(criterion, 0);
+    }
+};
+
+#endif /* COSMOHAMTAGGINGCRITERION_HPP_ */