diff --git a/Examples/Cosmo/CosmoLevel.cpp b/Examples/Cosmo/CosmoLevel.cpp index 0b501e38b..9ec48ba82 100644 --- a/Examples/Cosmo/CosmoLevel.cpp +++ b/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 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( + 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 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( + 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 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 ¤t_state, const FArrayBox ¤t_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() { diff --git a/Examples/Cosmo/SimulationParameters.hpp b/Examples/Cosmo/SimulationParameters.hpp index 36e9d4867..ad3c169d6 100644 --- a/Examples/Cosmo/SimulationParameters.hpp +++ b/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 center_tag; InitialScalarData::params_t initial_params; Potential::params_t potential_params; KerrBH::params_t kerr_params; diff --git a/Examples/Cosmo/params_cheap.txt b/Examples/Cosmo/params_cheap.txt index e289e8dba..a94b5eaad 100644 --- a/Examples/Cosmo/params_cheap.txt +++ b/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 diff --git a/Examples/Cosmo/plot.py b/Examples/Cosmo/plot.py index 2c201e6ef..b918c485c 100644 --- a/Examples/Cosmo/plot.py +++ b/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() \ No newline at end of file diff --git a/Examples/Cosmo/plot_efold.pdf b/Examples/Cosmo/plot_efold.pdf index c98cbe567..e3b829054 100644 Binary files a/Examples/Cosmo/plot_efold.pdf and b/Examples/Cosmo/plot_efold.pdf differ diff --git a/Examples/Cosmo/plot_lineout.pdf b/Examples/Cosmo/plot_lineout.pdf index 97318f506..c4e9fa84c 100644 Binary files a/Examples/Cosmo/plot_lineout.pdf and b/Examples/Cosmo/plot_lineout.pdf differ diff --git a/Examples/Cosmo/plot_rho_mean.pdf b/Examples/Cosmo/plot_rho_mean.pdf index 950b3de94..a43c239bc 100644 Binary files a/Examples/Cosmo/plot_rho_mean.pdf and b/Examples/Cosmo/plot_rho_mean.pdf differ diff --git a/Source/TaggingCriteria/CosmoHamTaggingCriterion.hpp b/Source/TaggingCriteria/CosmoHamTaggingCriterion.hpp new file mode 100644 index 000000000..6089d14bf --- /dev/null +++ b/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 m_center_tag; + double m_rad; + double m_rho_mean; + + public: + CosmoHamTaggingCriterion(double dx, + std::array center_tag, + double rad, double rho_mean) + : m_dx(dx), m_center_tag(center_tag), m_rad(rad), + m_rho_mean(rho_mean){}; + + template void compute(Cell 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 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_ */