Commit 973588a0 authored by Phil Bull's avatar Phil Bull

BROKEN - Did initial modifications for variable bang time, not going well so far

parent 46fabc99
......@@ -13,6 +13,8 @@ int main(int argc, char* argv[]){
//if(argc<2){cerr << "ERROR: Not enough command-line arguments." << endl; return 1;}
//double rrr = atof(argv[1]);
// NOTE: Assumes that tB(r=0) = 0.
// Instantiate a model
/// Interesting one: CFLb mod(0.72, -5.75e-8, 1.0e-9); // (h, A, KR)
/// Another interesting one, with cusp!: CFLb mod(0.72, 0.5, 1.5e3); // (h, A, w)
......@@ -33,9 +35,6 @@ int main(int argc, char* argv[]){
//mod.output_fnz();
//mod.output_zin();
// TESTING
//mod.v_dipole(rrr);
// Finish timing the run
cerr << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
return 0;
......
......@@ -28,6 +28,170 @@ double Model::da1dot_da(double a, double m, double k){
}
int Model::get_max_j(double x){
/// Get the j value which will ensure that the grid will be
/// populated at least up to the maximum global time for a given x
int mid;
int a = 0;
int b = g_y.size()-1;
// Target y
double y = y_today - xTb(x);
// Get i<=>y by searching g_y
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(y <= g_y[mid]){
b = mid; // a <= y <= mid
}else{
a = mid; // mid < y <= b
} // end range check
} // end y-searching loop
return a; // Index of the maximum y value to calculate for this x
}
void Model::populate_timesample_vector(){
/// Populate the time sampling vector, which gives the local times
/// on the grid. Can use this to do variable time sampling
// FIXME: Basic (fixed) step-size definition for now
TBMIN = 0.0; // FIXME: Should set TBMIN (minimum value of Bang-time fn) in the defined models
double ystep = (y_today - (TBMIN/ts))/((double)YCELLS);
// Make YCELLS no. of samples
for(int i=0; i<YCELLS; i++){
g_y.push_back( y0 + (double)i*ystep );
} // end loop through cells
}
void Model::populate_grid(){
/// Populate a mesh with values for a1 using an RK4 ODE solution
/// to the LTB Friedmann equation
clock_t tclock = clock();
double x, y; int jmax;
//bool reached_asymptotic = false;
// Set-up arrays etc.
populate_timesample_vector();
// Initialise the GSL adaptive integrator
//const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd; // RK8 Prince-Dormand adaptive algorithm
const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45; // Runge-Kutta-Fehlberg (4,5) adaptive algorithm
gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1); // Allocate a 1D ODE array
gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-12, 0.0); // Set the relative/absolute error tolerances
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
int status; // Variable to hold status of GSL integration
// Make sure we start at the right place
double hh; double a1[1] = { A1_REC };
x = x0 - XSTEP;
ModelPointers p; p.mod = this;
// Loop through x values
for(int i=0; i<XCELLS+2 /*FIXME and !reached_asymptotic*/; i++){
// Set-up for integration
x += XSTEP; g_x.push_back(x); // Add the x value for this cell to the array
p.m = xM(x); p.k = xK(x);
a1[0] = A1_REC;
// Calculate the maximum j value to go up to
jmax = get_max_j(x);
// Create new temporary time slice vectors
// Add an extra 2 elements to each slice to avoid boundary problems when interpolating
vector<double> tmp_a1(jmax+2, 0.0);
vector<double> tmp_a1dot(jmax+2, 0.0);
// Define system of ODEs (ODE + Jacobian)
gsl_odeiv_system sys = {static_a1dot, a1dot_jac, 1, &p};
// Loop through sample points
y = y0; hh = 1e-16; // Make sure we reset hh
for(int j=0; j<jmax; j++){
// Integrate until the next sample point
while (y < g_y[j]){
status = gsl_odeiv_evolve_apply (e, c, s, &sys, &y, g_y[j], &hh, a1);
if (status != GSL_SUCCESS){ cerr << "Model::populate_grid(): Integration failed." << endl; break; }
} // end while loop
if(isnan(a1[0])){throw 111;}
// Add to mesh at this sample point
tmp_a1[j] = a1[0];
tmp_a1dot[j] = s_a1dot(a1[0], p.m, p.k)/ts;
} // end time steps
g_a1.push_back( tmp_a1 );
g_a1dot.push_back( tmp_a1dot );
} // end loop through x values
// Get the a2 values now the dust has settled
calculate_a2grid();
// Free all of the integration variables
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
cerr << "\tpopulate_grid(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
} // end initial integration function
void Model::calculate_a2grid(){
/// Calculate the a2 values on the grid now that we have a1, a1dot values
double x, y, r, t; int jmax;
double dr = XSTEP*rs; double invdr = 1./(2.*dr);
// TESTING
double a11, a12, a13, a14, a15, a16, a17, a18;
// Loop through x values
for(int i=0; i<XCELLS+2; i++){
x = g_x[i]; r = x*rs;
jmax = get_max_j(x); // Calculate the maximum j value to go up to
// Create new temporary time slice vectors
// Add an extra 2 elements to each slice to avoid boundary problems when interpolating
vector<double> tmp_a2(jmax+2, 0.0);
vector<double> tmp_a2dot(jmax+2, 0.0);
for(int j=0; j<jmax; j++){
y = g_y[j] + xTb(x);
t = y*ts;
if(i>1){
a11 = a1(r, t);
a12 = a1(r+dr, t);
a13 = a1(r-dr, t);
a15 = a1(r-2.*dr, t);
a16 = a1(r+2.*dr, t);
a17 = a1(r+3.*dr, t);
a14 = r*(a1(r+dr, t) - a1(r-dr, t))*invdr;
tmp_a2[j] = a11 + a14;
//tmp_a2[j] = a1(r, t) + r*(a1(r+dr, t) - a1(r-dr, t))*invdr;
tmp_a2dot[j] = a1dot(r, t) + r*(a1dot(r+dr, t) - a1dot(r-dr, t))*invdr;
if(tmp_a2[j] <= 0.0){cerr << "WARNING: Shell crossing occurred at r=" << x*rs << ", y/y_today=" << y/y_today << " (NumericGrid::populate_grid)" << endl; throw 99;} // end check
}else{
// a2, a2dot = a1, a1dot at r=0
tmp_a2[j] = a1(r, t);
tmp_a2dot[j] = a1dot(r, t);
}
} // end j loop y values
g_a2.push_back( tmp_a2 );
g_a2dot.push_back( tmp_a2dot );
} // end i loop through x values
}
/*
void Model::populate_grid(){
/// Populate a mesh with values for a1 using an RK4 ODE solution
/// to the LTB Friedmann equation
......@@ -132,3 +296,4 @@ void Model::populate_grid(){
cerr << "\tpopulate_grid(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
} // end initial integration function
*/
......@@ -20,18 +20,20 @@ double Model::get_array_value(double r, double t, int v){
int i, j, a, b, mid;
double x1, x2, y1, y2;
double xcorr = 0.0; double ycorr = 0.0;
a = 0; b = g_x.size() - 3; // Since XCELLS = xsize - 2
double x = r/rs;
double y = t/ts;
double y = t/ts - xTb(x);
double gxsize = g_x.size();
double gysize = g_y.size();
// Sanity checks on x, y
if(x < g_x[a]){
if(x < g_x[0]){
if(x < 0.0){
cerr << "WARNING: x went negative, taking fabs(x) [Model::get_array_value]" << endl;
x = fabs(x); // Try the +ve value of x
}else{
cerr << "ERROR: x=" << x << " outside calculated range: " << g_x[a] << " to " << g_x[b] << endl; throw 71;
cerr << "ERROR: x=" << x << " outside calculated range: " << g_x[0] << " to " << g_x[gxsize-1] << endl; throw 71;
} // end check whether x is +ve
}
// FIXME: Sometimes this gets thrown when compiled with -O3 (but seemingly doesn't happen when DEBUG flags used)
......@@ -42,11 +44,12 @@ double Model::get_array_value(double r, double t, int v){
}
// Check if x is too big, and if so, get the asymptotic value
if(x >= g_x[b-1]){
x = g_x[b-1]; // Set x to value where expn. vars. reach asymptotic values, it'll return the same results
if(x >= g_x[gxsize-1]){
x = g_x[gxsize-1]; // Set x to value where expn. vars. reach asymptotic values, it'll return the same results
}
// Get i<=>x by searching g_x
a = 0; b = gxsize - 3; // Since XCELLS = xsize - 2
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(x <= g_x[mid]){
......@@ -57,28 +60,43 @@ double Model::get_array_value(double r, double t, int v){
} // end x-searching loop
i = a;
// Get j<=>y by searching g_y
a = 0; b = gysize - 1;
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(y <= g_y[mid]){
b = mid; // a <= y <= mid
}else{
a = mid; // mid < y <= b
} // end range check
} // end y-searching loop
j = a;
// Get j<=>y by looking it up from the index (assumes fixed grid cell size in y)
j = (int)floor( ((y - y0) / (y_today - y0)) * (double)YCELLS);
if(j >= YCELLS){j = YCELLS - 1;} // Check we're not right at the top of the range
//j = (int)floor( ((y - y0) / (y_today - y0)) * (double)YCELLS);
//if(j >= YCELLS){j = YCELLS - 1;} // Check we're not right at the top of the range
// Check if x, y indices are safe to do forward interpolation
double gxsize = g_x.size();
if(i >= gxsize and j < YCELLS){
if(i >= gxsize and j < gysize){
// We're at the x boundary. FIXME: Check for validity
xcorr = g_x[i] - g_x[i-1]; // Correct for shift in cell used to calculate gradient
i -= 1; // We're at the x boundary, interpolate using a backwards step (fix the indices)
}else if(i < gxsize and j >= YCELLS){
ycorr = YSTEP; j -= 1; // We're at the y boundary
}else if(i >= gxsize and j >= YCELLS){
}else if(i < gxsize and j >= gysize){
ycorr = g_y[j] - g_y[j-1]; j -= 1; // We're at the y boundary
}else if(i >= gxsize and j >= gysize){
xcorr = g_x[i] - g_x[i-1];
ycorr = YSTEP; i -= 1; j-=1; // We're at the x and y boundary!
ycorr = g_y[j] - g_y[j-1]; i -= 1; j-=1; // We're at the x and y boundary!
}
// Get values for x1, x2, y1, y2
x1 = g_x[i]; x2 = g_x[i+1];
y1 = y0 + ((double)j * YSTEP); y2 = y1 + YSTEP;
y1 = g_y[j]; y2 = g_y[j+1];
//cerr << "x1=" << x1 << ", x2=" << x2 << ", y1=" << y1 << ", y2=" << y2 << endl;
//cerr << "x=" << x << ", y=" << y << endl;
// Do a simple (bilinear) interpolation
// FIXME: This is wrong! May not be [i, j+1]
double arr1, arr2, arr3, arr4;
switch(v){
case 0: arr1 = g_a1[i][j]; arr2 = g_a1[i][j+1];
......@@ -93,7 +111,25 @@ double Model::get_array_value(double r, double t, int v){
break;
}
return (1./(YSTEP*(x2 - x1))) * (
/*
cerr << "-------------\nV = " << v << endl;
cerr << " ca[i,j]=" << arr1 * (ycorr + y2 - y)*(xcorr + x2 - x)
<< ", ca[i,j+1]=" << arr2 * (ycorr + y - y1)*(xcorr + x2 - x)
<< ", ca[i+1,j]=" << arr3 * (ycorr + y2 - y)*(xcorr + x - x1)
<< ", ca[i+1,j+1]=" << arr4 * (ycorr + y - y1)*(xcorr + x - x1)
<< endl;
cerr << " a[i,j]=" << arr1
<< ", a[i,j+1]=" << arr2
<< ", a[i+1,j]=" << arr3
<< ", a[i+1,j+1]=" << arr4
<< endl;
cerr << "(ycorr + y2 - y)=" << (ycorr + y2 - y)
<< ", (xcorr + x2 - x)=" << (xcorr + x2 - x) << endl;
cerr << " ycorr=" << ycorr << ", xcorr=" << xcorr << endl;
*/
return (1./((y2 - y1)*(x2 - x1))) * (
arr1 * (ycorr + y2 - y)*(xcorr + x2 - x)
+ arr2 * (ycorr + y - y1)*(xcorr + x2 - x)
+ arr3 * (ycorr + y2 - y)*(xcorr + x - x1)
......
......@@ -9,6 +9,7 @@ protected:
vector<vector <double> > g_a2;
vector<vector <double> > g_a2dot;
vector<double> g_x;
vector<double> g_y;
vector<double> rval;
vector<double> tval;
......@@ -21,6 +22,10 @@ protected:
void generate_geodesic();
void hubble_time();
int get_max_j(double x);
void populate_timesample_vector();
void calculate_a2grid();
double z(double r);
double a1(double r, double t);
......@@ -62,6 +67,7 @@ public:
int XCELLS, YCELLS;
double XMAX, XMIN, XSTEP, YSTEP;
double TBMIN;
double a1_0, x0, y0, dy, y_today, t_today, r0, t0;
double dt;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment