Skip to content
Snippets Groups Projects
Commit da40ef0e authored by Izzard, Robert Dr (Maths & Physics)'s avatar Izzard, Robert Dr (Maths & Physics)
Browse files

update logging for discs

parent 8bbfac9a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/disc/disc_logging.c
+ 5
− 2
View file @ da40ef0e
...@@ -196,7 +196,7 @@ void disc_logging(struct stardata_t * const stardata, ...@@ -196,7 +196,7 @@ void disc_logging(struct stardata_t * const stardata,
double Revap_in = disc->Revap_in; double Revap_in = disc->Revap_in;
/* output */ /* output */
Printf("DISC_EV %20.10e %20.10e %d %d %u %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n", Printf("DISC_EV %20.10e %20.10e %d %d %u %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
dtp, // 1 dtp, // 1
disc->lifetime/YEAR_LENGTH_IN_SECONDS, disc->lifetime/YEAR_LENGTH_IN_SECONDS,
stardata->star[0].stellar_type, stardata->star[0].stellar_type,
...@@ -248,7 +248,10 @@ void disc_logging(struct stardata_t * const stardata, ...@@ -248,7 +248,10 @@ void disc_logging(struct stardata_t * const stardata,
disc_total_gravitational_potential_energy(disc,b), disc_total_gravitational_potential_energy(disc,b),
disc_total_kinetic_energy(disc,b), disc_total_kinetic_energy(disc,b),
disc->sigma0, // 51 disc->sigma0, // 51
disc->Tvisc0 disc->Tvisc0,
disc->Rout / disc->Rin,
disc->Rin / Max(1e-50,b->separation),
disc->Rout / Max(1e-50,b->separation)
); );
} }
......
src/maths/GSL_integrator.c
+ 5
− 5
View file @ da40ef0e
#include "../binary_c.h" #include "../binary_c.h"
#include <gsl/gsl_integration.h> #include <gsl/gsl_integration.h>
double GSL_integrator(double lower, double GSL_integrator(const double lower,
double upper, const double upper,
double relative_tolerance, const double relative_tolerance,
double * error, double * const error,
int integrator, const int integrator,
double func(double x,void *params), double func(double x,void *params),
...) ...)
{ {
......
src/maths/maths_prototypes.h
+ 5
− 5
View file @ da40ef0e
...@@ -156,11 +156,11 @@ double polyalgorithm(const int s, ...@@ -156,11 +156,11 @@ double polyalgorithm(const int s,
Boolean inverse3( const double theMatrix [/*Y=*/3] [/*X=*/3], Boolean inverse3( const double theMatrix [/*Y=*/3] [/*X=*/3],
double theOutput [/*Y=*/3] [/*X=*/3] ); double theOutput [/*Y=*/3] [/*X=*/3] );
double GSL_integrator(double lower, double GSL_integrator(const double lower,
double upper, const double upper,
double tolerance, const double tolerance,
double * error, double * const error,
int integrator, const int integrator,
double func(double x,void *params), double func(double x,void *params),
...); ...);
void setup_GSL_handlers(struct stardata_t * RESTRICT const stardata); void setup_GSL_handlers(struct stardata_t * RESTRICT const stardata);
......
src/perl/scripts2/cbdiscs.pl
+ 336
− 172
View file @ da40ef0e
...@@ -105,7 +105,7 @@ $population->set( ...@@ -105,7 +105,7 @@ $population->set(
cbdisc_alpha => 1e-3, cbdisc_alpha => 1e-3,
cbdisc_kappa => 1e-2, cbdisc_kappa => 1e-2,
cbdisc_torquef => 1e-4, # 1e-4 cbdisc_torquef => 1e-4, # 1e-4
cbdisc_init_dM => 0.0, # 0 means do not adapt initial masses cbdisc_init_dM => 0.0, # 0 means do not adapt initial masses
cbdisc_init_dJdM => 0.5, cbdisc_init_dJdM => 0.5,
...@@ -148,15 +148,15 @@ if($ARGV[0] eq 'show_opts') ...@@ -148,15 +148,15 @@ if($ARGV[0] eq 'show_opts')
} }
$population->set( $population->set(
# slurm # slurm
slurm=>0, # use slurm if 1 (please set on command line) slurm=>0, # use slurm if 1 (please set on command line)
slurm_njobs=>4, # number of jobs slurm_njobs=>4, # number of jobs
slurm_dir=>'/users/ri0005/data/slurm', # working directory for scripts etc. slurm_dir=>'/users/ri0005/data/slurm', # working directory for scripts etc.
slurm_memory=>512, # RAM (MB) per job slurm_memory=>512, # RAM (MB) per job
slurm_postpone_join=>0, # if 1 do not join on slurm, join elsewhere slurm_postpone_join=>0, # if 1 do not join on slurm, join elsewhere
slurm_partition=>'all', # MUST be defined slurm_partition=>'all', # MUST be defined
slurm_jobname=>'binary_cbd', # not required but useful slurm_jobname=>'binary_cbd', # not required but useful
slurm_use_all_node_CPUs=>1, # required on eureka slurm_use_all_node_CPUs=>1, # required on eureka
); );
# scan command line arguments for extra options # scan command line arguments for extra options
...@@ -198,62 +198,63 @@ elsif($duplicity == 1) ...@@ -198,62 +198,63 @@ elsif($duplicity == 1)
q => $n, q => $n,
P => $n P => $n
}; };
my $mmin = 0.85; my $mmin = 0.85;
my $mmax = 6; my $mmax = 6;
$population->{_grid_options}{binary} = 1; $population->{_grid_options}{binary} = 1;
$population->add_grid_variable $population->add_grid_variable
( (
'name' => 'lnm1', 'name' => 'lnm1',
'longname' =>'Primary mass', 'longname' =>'Primary mass',
'range' =>[log($mmin),log($mmax)], 'range' =>[log($mmin),log($mmax)],
'resolution' => $resolution->{m1}, 'resolution' => $resolution->{m1},
'spacingfunc'=>"const(log($mmin),log($mmax),$resolution->{m1})", 'spacingfunc'=>"const(log($mmin),log($mmax),$resolution->{m1})",
'precode' =>'$m1=exp($lnm1);', 'precode' =>'$m1=exp($lnm1);',
'probdist' =>'Kroupa2001($m1)*$m1', 'probdist' =>'Kroupa2001($m1)*$m1',
'dphasevol' =>'$dlnm1', 'dphasevol' =>'$dlnm1',
); );
# q=M1/M2 distribution flat in q between 0.1/M1 and 1.0 # q=M1/M2 distribution flat in q between 0.1/M1 and 1.0
$population->add_grid_variable $population->add_grid_variable
( (
'condition' =>'$self->{_grid_options}{binary}==1', 'condition' =>'$self->{_grid_options}{binary}==1',
'name' =>'q', 'name' =>'q',
'longname' =>'Mass ratio', 'longname' =>'Mass ratio',
'range' =>['0.1/$m1',1.0], 'range' =>['0.1/$m1',1.0],
'resolution'=>$resolution->{q}, 'resolution'=>$resolution->{q},
'spacingfunc'=>"const(0.1/\$m1,1.0,$resolution->{q})", 'spacingfunc'=>"const(0.1/\$m1,1.0,$resolution->{q})",
'probdist' =>"flatsections\(\$q,\[ 'probdist' =>"flatsections\(\$q,\[
\{min=>0.1/\$m1,max=>0.8,height=>1.0\}, \{min=>0.1/\$m1,max=>0.8,height=>1.0\},
\{min=>0.8,max=>1.0,height=>1.0\}, \{min=>0.8,max=>1.0,height=>1.0\},
\]\)", \]\)",
precode =>'$m2=$q*$m1;', precode =>'$m2=$q*$m1;',
dphasevol =>'$dq', dphasevol =>'$dq',
); );
# orbital period Duquennoy and Mayor 1991 distribution # orbital period Duquennoy and Mayor 1991 distribution
$population->add_grid_variable $population->add_grid_variable
( (
'name' =>'logper', 'name' =>'logper',
'longname' =>'log(Orbital_Period)', 'longname' =>'log(Orbital_Period)',
# 'range' =>[-2.0,12.0], # 'range' =>[-2.0,12.0],
'range' =>[1.0,4.0], 'range' =>[1.0,4.0],
'resolution' =>$resolution->{P}, 'resolution' =>$resolution->{P},
'spacingfunc'=>"const(-1.0,10.0,$resolution->{P})", 'spacingfunc'=>"const(-1.0,10.0,$resolution->{P})",
'precode'=>" 'precode'=>"
\$per = 10.0 ** \$logper; \$per = 10.0 ** \$logper;
my \$eccentricity = 0.0; my \$eccentricity = 0.0;
\$sep=calc_sep_from_period(\$m1,\$m2,\$per) if(defined \$m1 && defined \$m2); \$sep=calc_sep_from_period(\$m1,\$m2,\$per) if(defined \$m1 && defined \$m2);
", ",
'postcode' => " 'postcode' => "
", ",
'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)", 'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)",
'dphasevol'=>'$dln10per' 'dphasevol'=>'$dln10per'
); );
} }
else else
{ {
...@@ -327,6 +328,214 @@ $population->require( ...@@ -327,6 +328,214 @@ $population->require(
'BISECT_FAIL_EXIT is off', 'BISECT_FAIL_EXIT is off',
); );
# list of variables in each data line, in order
# (note that dtp is the first data item, always)
#
# those that are associated with undef are not plotted
#
# those given an empty hash are plotted with no options
#
# otherwise a hash {...} can contain directives for each plot,
# e.g. whether data is binned, logged, multiplied or prepended
# with 'log'
my $vars = [
age => {
bin => 0.1,
log => 1
},
stellar_type0 => {},
stellar_type1 => {},
nzones => undef,
zone1_type => undef,
zone2_type => undef,
zone3_type => undef,
Rin => {
bin => 0.025,
log => 1,
prepend => 1,
},
Rout => {
bin => 0.025,
log => 1,
prepend => 1,
},
separation => {
bin => 0.025,
log => 1,
},
orbital_period => {
bin => 0.1, # years
multiplier => 365.25,
log => 1
},
Rism => undef,
Xray_domination => undef,
R0 => undef,
R1 => undef,
RL0 => undef,
RL1 => undef,
Ldisc =>
{
bin => 0.025,
log => 1,
prepend => 1,
},
Idisc => undef,
Jdisc => {
bin => 0.1,
log => 1,
prepend => 1,
},
Jbinary => {
bin => 0.025,
log => 1,
prepend => 1,
},
Mdisc => {
bin => 0.1,
log => 1,
prepend => 1,
},
Mbinary => {
bin => 0.025,
log => 1,
prepend => 1,
},
Mdot_viscous => {
bin => 0.1,
log => 1,
},
Mdot_global => {
bin => 0.1,
log => 1,
},
Mdot_inner_L2 => {
bin => 0.1,
log => 1,
},
Mdot_ISM => {
bin => 0.1,
log => 1,
},
Mdot_FUV => {
bin => 0.1,
log => 1,
},
Mdot_Xray => {
bin => 0.1,
log => 1,
},
eccentricity => {
bin => 0.025,
},
edot_resonance => {
bin => 0.1,
},
Tin => {
bin => 0.025,
log => 1,
prepend => 1,
},
Tout => {
bin => 0.025,
log => 1,
prepend => 1,
},
densityin => {
bin => 0.1,
log => 1,
prepend => 1,
},
densityout => {
bin => 0.1,
log => 1,
prepend => 1,
},
Pgrav_in => {
bin => 0.1,
log => 1,
prepend => 1,
},
Pgrav_out => {
bin => 0.1,
log => 1,
prepend => 1,
},
PISM => undef,
HRin => {
bin => 0.025,
log => 1,
prepend => 1,
},
HRout => {
bin => 0.025,
log => 1,
prepend => 1,
},
Teffmax => {
bin => 0.025,
log => 1,
prepend => 1,
},
Tin_over_Tstar => {
bin => 0.1,
},
Lstar => {
bin => 0.025,
prepend => 1,
log => 1,
},
LXstar => {
bin => 0.1,
prepend => 1,
log => 1,
},
Ldisc_over_Lstar => undef,
tvisc_in => {
bin => 0.1,
prepend => 1,
log => 1,
},
tvisc_out => {
bin => 0.1,
prepend => 1,
log => 1,
},
PE => undef,
KE => undef,
sigma0 => undef,
Tvisc0 => undef,
Rout_over_Rin => {
bin => 0.025,
log => 1,
},
Rin_over_separation => {
bin => 0.025,
log => 1,
},
Rout_over_separation => {
bin => 0.025,
log => 1,
},
];
# 2D plots : binning, logging etc. is as in 1D
my $twoD = {
'period_eccentricity' =>
{
x => 'orbital_period',
y => 'eccentricity',
},
'period_luminosity' =>
{
x => 'orbital_period',
y => 'logLstar',
},
};
# evolution the stellar population (this takes some time) # evolution the stellar population (this takes some time)
$population->evolve(); $population->evolve();
...@@ -346,13 +555,16 @@ sub parse_data ...@@ -346,13 +555,16 @@ sub parse_data
my $progenitor_stellar_type; my $progenitor_stellar_type;
my $disc_donor_stellar_type; my $disc_donor_stellar_type;
while(1) while(1)
{ {
# subsequent calls to tbse_line contain # subsequent calls to tbse_line contain
# (references to) arrays of data # (references to) arrays of data
my $la = $population->tbse_line(); my $la = $population->tbse_line();
# first element is the "header" line # first element is the "header" line
my $header = shift @$la; my $header = shift @$la;
...@@ -366,116 +578,78 @@ sub parse_data ...@@ -366,116 +578,78 @@ sub parse_data
#print "$header\n";# @$la\n", $population->{_grid_options}->{args},"\n"; #print "$header\n";# @$la\n", $population->{_grid_options}->{args},"\n";
if($header eq 'DISC_EV') if($header eq 'DISC_EV')
{ {
#print "DISC_EV @$la\n"; # first is always dtp : remove this
my ($dtp, my $dtp = shift @$la;
$t, #print "DISC EV : dtp = $dtp\n";
$stellar_type0, #print "ARRAY @$la\n";
$stellar_type1,
$nzones, # construct the $data hash
$zone1_type, my $data = {};
$zone2_type, my $c = 0;
$zone3_type, while(scalar @$la)
$Rin, {
$Rout, my $y = shift @$la;
$separation, # Rsun my $pre = '';
$orbital_period, # years my $k = $vars->[$c++]; # variable key
$Rism, my $h = $vars->[$c++]; # variable hash
$Xray_domination,
$R0, if(defined $h)
$R1, {
$RL0, # extra things to do...
$RL1,
$Ldisc, # multiplier
$Idisc, if(defined $h->{multiplier})
$Jdisc, {
$Jbinary, $y *= $h->{multiplier};
$Mdisc, }
$Mbinary,
$Mdot_viscous, # log10
$Mdot_global, if(defined $h->{log} && $h->{log})
$Mdot_inner_L2, {
$Mdot_ISM, $y = safelog10($y);
$Mdot_FUV,
$Mdot_Xray, # prepend 'log' to the data key
$eccentricity, if($h->{prepend})
$edot_resonance, {
$Tin, $pre = 'log';
$Tout, }
$density_in, }
$density_out,
$Pgrav_in, # bin data
$Pgrav_out, if(defined $h->{bin})
$PISM, {
$HRin, $y = $population->rebin($y,$h->{bin});
$HRout, }
$Teffmax,
$Tin_over_Tstar, # set the data hash : note that these data are unique
$Lstar, # there should be no existing data in $data->{$pre.$k}
$LXstar, $data->{$pre.$k}->{$y} = $dtp;
$Ldisc_over_Lstar, }
$tvisc_in, }
$tvisc_out,
$PE,
$KE,
$sigma0,
$Tvisc0) = @$la;
my $r; my $r;
foreach my $set ('all', foreach my $set ('all',
$progenitor_stellar_type) $progenitor_stellar_type)
{ {
$r = $results->{$set}; $r = $results->{$set};
# set the data
foreach my $x (keys %$data)
{
foreach my $y (keys %{$data->{$x}})
{
$results->{$set}->{$x}->{$y} += $data->{$x}->{$y};
}
}
# prebin some quantities
# we want orbital period in days
my $orbital_period_binned = $population->rebin(safelog10($orbital_period*365.25),0.1);
my $eccentricity_binned = $population->rebin($eccentricity,0.025);
my $logLstar_binned = $population->rebin(safelog10($Lstar),0.025);
# 1D distributions
$r->{Mdot_viscous}->{$population->rebin(safelog10($Mdot_viscous),0.1)} += $dtp;
$r->{Mdot_global}->{$population->rebin(safelog10($Mdot_global),0.1)} += $dtp;
$r->{Mdot_inner_L2}->{$population->rebin(safelog10($Mdot_inner_L2),0.1)} += $dtp;
$r->{Mdot_ISM}->{$population->rebin(safelog10($Mdot_ISM),0.1)} += $dtp;
$r->{Mdot_FUV}->{$population->rebin(safelog10($Mdot_FUV),0.1)} += $dtp;
$r->{Mdot_Xray}->{$population->rebin(safelog10($Mdot_Xray),0.1)} += $dtp;
$r->{logMdisc}->{$population->rebin(safelog10($Mdisc),0.1)} += $dtp;
$r->{logMbinary}->{$population->rebin(safelog10($Mbinary),0.025)} += $dtp;
$r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.025)} += $dtp;
$r->{logJdisc}->{$population->rebin(safelog10($Jdisc),0.1)} += $dtp;
$r->{logJbinary}->{$population->rebin(safelog10($Jbinary),0.025)} += $dtp;
$r->{age}->{$population->rebin(safelog10($t),0.1)} += $dtp;
$r->{separation}->{$population->rebin(safelog10($separation),0.025)} += $dtp;
$r->{orbital_period}->{$orbital_period_binned} += $dtp;
$r->{eccentricity}->{$eccentricity_binned} += $dtp;
$r->{edot_resonance}->{$population->rebin($edot_resonance,0.1)} += $dtp;
$r->{stellar_type0}->{$stellar_type0} += $dtp;
$r->{stellar_type1}->{$stellar_type1} += $dtp;
$r->{logRin}->{$population->rebin(safelog10($Rin),0.025)} += $dtp;
$r->{logRout}->{$population->rebin(safelog10($Rout),0.025)} += $dtp;
$r->{logRout_over_Rin}->{$population->rebin(safelog10($Rout/MAX(1e-50,$Rin)),0.025)} += $dtp;
$r->{logRin_over_separation}->{$population->rebin(safelog10($Rin/MAX(1e-50,$separation)),0.025)} += $dtp;
$r->{logRout_over_separation}->{$population->rebin(safelog10($Rout/MAX(1e-50,$separation)),0.025)} += $dtp;
$r->{logHRin}->{$population->rebin(safelog10($HRin),0.025)} += $dtp;
$r->{logHRout}->{$population->rebin(safelog10($HRout),0.025)} += $dtp;
$r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.1)} += $dtp;
$r->{logLstar}->{$logLstar_binned} += $dtp;
$r->{logLXstar}->{$population->rebin(safelog10($LXstar),0.1)} += $dtp;
$r->{logdensityin}->{$population->rebin(safelog10($density_in),0.1)} += $dtp;
$r->{logdensityout}->{$population->rebin(safelog10($density_out),0.1)} += $dtp;
$r->{logTin}->{$population->rebin(safelog10($Tin),0.025)} += $dtp;
$r->{logTout}->{$population->rebin(safelog10($Tout),0.025)} += $dtp;
$r->{logPgrav_in}->{$population->rebin(safelog10($Pgrav_in),0.1)} += $dtp;
$r->{logPgrav_out}->{$population->rebin(safelog10($Pgrav_out),0.1)} += $dtp;
$r->{logTeffmax}->{$population->rebin(safelog10($Teffmax),0.025)} += $dtp;
$r->{logtvisc_in}->{$population->rebin(safelog10($tvisc_in),0.1)} += $dtp;
# 2D plots # 2D plots
$r->{'2D_period_eccentricity'}->{$orbital_period_binned}->{$eccentricity_binned} += $dtp; foreach my $plot (keys %$twoD)
$r->{'2D_period_luminosity'}->{$orbital_period_binned}->{$logLstar_binned} += $dtp; {
my $x = (keys %{$data->{$twoD->{$plot}->{x}}})[0];
#print "2D $orbital_period_binned $eccentricity_binned $dtp\n"; my $y = (keys %{$data->{$twoD->{$plot}->{y}}})[0];
$r->{'2D_'.$plot}->{$x}->{$y} += $dtp;
}
} }
} }
elsif($header eq 'DISC_END') elsif($header eq 'DISC_END')
...@@ -499,18 +673,11 @@ sub parse_data ...@@ -499,18 +673,11 @@ sub parse_data
if(defined $systems{$x}) if(defined $systems{$x})
{ {
print "ERROR : double disc end ($systems{$x})\n$population->{_grid_options}->{args}\n"; print "ERROR : double disc end ($systems{$x})\n$population->{_grid_options}->{args}\n";
} }
$systems{$x}++; $systems{$x}++;
} }
printf "DISC_END Set prog type %s loglifetime = %g from %g += %g\n",
$progenitor_stellar_type,
$population->rebin(safelog10($lifetime),0.1),
$lifetime,
$p;
foreach my $set ('all', foreach my $set ('all',
$progenitor_stellar_type) $progenitor_stellar_type)
{ {
...@@ -537,7 +704,6 @@ sub parse_data ...@@ -537,7 +704,6 @@ sub parse_data
$p; $p;
print {$initparamspace_fp} $s; print {$initparamspace_fp} $s;
#print $s;
} }
} }
elsif($header eq 'DISC_START') elsif($header eq 'DISC_START')
...@@ -556,21 +722,19 @@ sub parse_data ...@@ -556,21 +722,19 @@ sub parse_data
$init_a, $init_a,
$init_e $init_e
) = @$la; ) = @$la;
print "Start donor type $donor_stellar_type (== $progenitor_stellar_type ?)\n";
foreach my $set ('all',
foreach my $set ('all',
$progenitor_stellar_type) $progenitor_stellar_type)
{ {
my $r = $results->{$set}; my $r = $results->{$set};
$r->{logstarttime}->{$population->rebin(safelog10($starttime),0.1)} += $p; $r->{logstarttime}->{$population->rebin(safelog10($starttime),0.1)} += $p;
$r->{initial_stellar_type0}->{$initial_stellar_type0} += $p; $r->{initial_stellar_type0}->{$initial_stellar_type0} += $p;
$r->{initial_stellar_type1}->{$initial_stellar_type1} += $p; $r->{initial_stellar_type1}->{$initial_stellar_type1} += $p;
$r->{overflower_stellar_type}->{$donor_stellar_type} += $p; $r->{overflower_stellar_type}->{$donor_stellar_type} += $p;
} }
# save for initial parameter space logging # save for initial parameter space logging
$disc_donor_stellar_type = $donor_stellar_type; $disc_donor_stellar_type = $donor_stellar_type;
#print "Donor type $disc_donor_stellar_type, prob $p\n";
} }
elsif($header eq 'DISC_COMENV') elsif($header eq 'DISC_COMENV')
{ {
......
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