From 1df517d8ee18fb5881c615349a8068d1566ca03a Mon Sep 17 00:00:00 2001 From: Robert Izzard <r.izzard@surrey.ac.uk> Date: Tue, 11 Dec 2012 12:40:29 +0000 Subject: [PATCH] update trunk to my latest version --- cf_imf_plots.pl | 1 - configure | 224 ++---- doc/binary_grid-flexigrid.lyx | 359 ++++++++- grid-IMF.pl | 1 - grid-wrap.pl | 372 +++++++++ grid.pl | 1 - imftickfile | 1 - src/binary_parameters.h | 13 +- src/binary_prototypes.h | 3 +- .../common_envelope_evolution.c | 2 +- src/binary_structures.h | 5 +- src/claet/claet_low_mass_GB.c | 2 +- src/claet/claet_not_a_helium_star.c | 3 +- src/claet/claet_post_HG.c | 6 +- src/code_options.h | 6 +- src/free_memory.c | 2 +- src/hrdiag/hrdiag.c | 7 +- src/hrdiag/hrdiag_HeStar.c | 4 +- src/iteration/end_of_iteration.c | 27 +- src/iteration/init_star.c | 2 +- src/iteration/initialize_parameters.c | 8 +- src/iteration/iterate_logging.c | 27 +- src/iteration/to_loop_or_not_to_loop.c | 4 +- src/logging/close_log_files.c | 2 - src/logging/log_hr.c | 7 +- src/logging/open_log_files.c | 1 - src/logging/output_string_to_log.c | 6 +- src/logging/output_to_logfile.c | 2 + src/main.c | 2 + src/main.h | 3 +- src/main2.c | 68 +- src/main_prototypes.h | 1 + src/maths/gaussj.c | 22 +- src/maths/interpolate.c | 30 +- src/nucsyn/nucsyn_elemental_abundance.c | 11 +- src/nucsyn/nucsyn_gce_log.c | 756 +++++++++--------- src/nucsyn/nucsyn_hbb.c | 32 +- src/nucsyn/nucsyn_initial_abundances.c | 1 - src/nucsyn/nucsyn_log.c | 9 +- ..._sequence_convective_envelope_size_table.c | 4 +- src/nucsyn/nucsyn_mc1tp.c | 2 +- src/nucsyn/nucsyn_parameters.h | 13 +- src/nucsyn/nucsyn_preset_sigmav.c | 4 +- src/nucsyn/nucsyn_radioactive_decay.c | 2 +- src/nucsyn/nucsyn_set_1st_dup_abunds.c | 27 +- src/nucsyn/nucsyn_set_envelope_abunds.c | 2 - src/nucsyn/nucsyn_set_third_dredgeup_abunds.c | 8 +- src/nucsyn/nucsyn_short_log.c | 5 +- src/nucsyn/nucsyn_sn_core_collapse.c | 2 - src/nucsyn/nucsyn_third_dredge_up.c | 22 +- src/nucsyn/nucsyn_update_abundances.c | 2 +- .../nucsyn_update_abundances_with_profiles.c | 2 - src/perl/install_modules.pl | 22 +- src/perl/modules/binary_grid.pm | 532 ++++++++---- src/perl/modules/binary_grid/condor.pm | 503 ++++++++++++ src/perl/modules/rob_misc.pm | 7 +- src/perl/modules/spacing_functions.pm | 78 +- src/perl/modules/sync | 18 +- src/setup/parse_arguments.c | 234 +++--- src/setup/set_up_variables.c | 4 +- src/setup/version.c | 11 + .../calc_stellar_wind_mass_loss.c | 2 +- src/wind/wind_AGB.c | 5 +- src/wind/wind_WR.c | 2 +- src/zfuncs/menvf.c | 2 +- tbse | 22 +- 66 files changed, 2518 insertions(+), 1052 deletions(-) delete mode 120000 cf_imf_plots.pl delete mode 120000 grid-IMF.pl create mode 100755 grid-wrap.pl delete mode 120000 grid.pl delete mode 120000 imftickfile create mode 100644 src/perl/modules/binary_grid/condor.pm diff --git a/cf_imf_plots.pl b/cf_imf_plots.pl deleted file mode 120000 index 00c085c98..000000000 --- a/cf_imf_plots.pl +++ /dev/null @@ -1 +0,0 @@ -src/perl/scripts/cf_imf_plots.pl \ No newline at end of file diff --git a/configure b/configure index 33520e6b7..052b5247b 100755 --- a/configure +++ b/configure @@ -28,18 +28,19 @@ use strict; # ./configure cgi # # For DUMA -# ./configure - +# ./configure duma +# +# For condor support +# ./configure condor -use vars qw($acmlflags $arch $auto_opt $autoprofile $beowulf $binary $bits @c_files $cc $cflags $cgi $condor $coptflags $cpu %cpuinfo $debug $duma $extrainc $extralibs $hostname %icc_opts $ipod $ld $ldflags $ldpath $lib $make $makefilename $nolic $nprocessors $opt $os $owd $profile_arcs $profilepass $src_dir $static %suncc_opts $tbse_test $v $vb $global_warn); +use vars qw($acmlflags $arch $auto_opt $autoprofile $binary $bits @c_files $cc $cflags $cgi $condor $coptflags $cpu %cpuinfo $debug $duma $extrainc $extralibs $hostname %icc_opts $ipod $ld $ldflags $ldpath $lib $make $makefilename $nolic $nprocessors $opt $os $owd $profile_arcs $profilepass $src_dir $static %suncc_opts $tbse_test $v $vb $global_warn); defaults(); parse_args(); secondary_variables(); -#make_unrolled_loop_functions(); -if($autoprofile==1) +if($autoprofile) { if($cc=~/[ig]cc/) { @@ -90,12 +91,12 @@ sub defaults if($ENV{'CC'}) { $cc=$ENV{'CC'}; - $ld="ld"; # linker + $ld='ld'; # linker } else { - $cc="gcc"; # c-compiler - $ld="gcc"; + $cc='gcc'; # c-compiler + $ld='gcc'; } $ldflags=''; # linker flags e.g. '-dM -E -c' for a macro list @@ -109,7 +110,7 @@ sub defaults $nolic=("@ARGV"=~/nolic/); $v=get_version(); $extralibs=''; - $binary="binary_c"; # output executable + $binary='binary_c'; # output executable $make=find_make(); # make program $auto_opt=1; # automatically select optimisation parameters $ipod=0; @@ -124,7 +125,7 @@ sub defaults # specifiy the first files to be complied (the slowest) if you want to @c_files=get_cfiles(); - #$hostname=$ENV{'HOSTNAME'}; + if($hostname eq '') { # try bash @@ -141,11 +142,7 @@ sub defaults close FP; } } - $beowulf=0; - if($hostname eq 'krypton') - { - $beowulf=1; - } + # 32 or 64 bit? if($cpuinfo{'flags'}=~/ lm /) { @@ -165,17 +162,6 @@ sub defaults } } - - # use condor? - if(`condor -version 2>&1`=~/CondorVersion/) - { - $condor=1; - } - else - { - $condor=0; - } - } sub secondary_variables @@ -522,7 +508,7 @@ sub secondary_variables # debugging if($debug == 1) { - $cflags='-O0 -g -fPIC '; + $cflags='-O0 -g -fPIC -rdynamic '; if($coptflags=~/-pg/) { $coptflags='-pg'; @@ -535,12 +521,7 @@ sub secondary_variables # standard options $cflags.=" $tbse_test "; - - if($beowulf==1) - { - $cflags.=" -DBEOWULF "; - } - + # OS specific options if($os eq "Linux") { @@ -561,8 +542,9 @@ sub secondary_variables $cflags.=" -DCPUFREQ=".$cpuinfo{'cpu MHz'}.' '; } - #$cflags.=' -DCOMPARE_FORTRAN '; - if($cgi==1) + $cflags.=' -DSVN_REVISION=$(SVN_REVISION) -DSVN_URL="$(SVN_URL)" '; + + if($cgi) { $lib=~s/lib64/lib/go; $acmlflags=''; @@ -573,7 +555,7 @@ sub secondary_variables $extralibs=~s/-lacml//o; } - if($ipod==1) + if($ipod) { # ipod cross-compilation $cc='arm-elf-gcc'; @@ -590,40 +572,43 @@ sub secondary_variables #$cflags=$coptflags=$extralibs=''; #$coptflags=" -ffloat-store -fmath-errno -fno-unsafe-math-optimizations -fno-finite-math-only -ftrapping-math -fno-rounding-math -fsignaling-nans-O0 -g "; - if($duma==1) + if($duma) { $extralibs.=' -lduma '; $cflags.=' -DDUMA '; } - - print "Host: $hostname\nOS: $os\nCpu: $cpu\nArch $arch\nBits: $bits\ncc: $cc\nld: $ld\nCflags : $cflags\nExtra Cflags for optimization: $coptflags\nExtra libs: $extralibs\n"; - if($acmlflags) + if($condor) { - print "acml: $acmlflags\n"; + # code should be portable across the condor platform + $coptflags=~s/-msse\S+//g; + $cflags=~s/-mtune=\S+/-mtune=generic/g; + $coptflags=~s/-mtune=\S+/-mtune=generic/g; + $cflags=~s/-march=\S+//g; + $coptflags=~s/-march=\S+//g; + $coptflags=~s/-m(3dnow|popcnt)//g; } + + print "Host: $hostname\nOS: $os\nCpu: $cpu\nArch $arch\nBits: $bits\ncc: $cc\nld: $ld\nCflags : $cflags\nExtra Cflags for optimization: $coptflags\nExtra libs: $extralibs\n"; + print "acml: $acmlflags\n" if($acmlflags); } sub get_os { my $uname=`uname`; chomp $uname; - if($uname eq "Linux") + if($uname eq 'Linux' || $uname eq 'Darwin') { - return("Linux"); - } - elsif($uname eq "Darwin") - { - return("Darwin"); + return $uname; } elsif($uname=~/sunos/) { - return("SunOS"); + return 'SunOS'; } else { print "Warning : unknown OS detected!\n"; - return($uname); + return $uname; } } @@ -634,20 +619,15 @@ sub get_proccpu # on a Mac (for RS!) find a different way of getting the # CPU ID, frequency etc. my $hostinfo=`hostinfo`; - if($hostinfo=~/(\d+) processors/) - { - $nprocessors=$1; - } - else - { - $nprocessors=1; - } + $nprocessors = ($hostinfo=~/(\d+) processors/) ? $1 : 1; + if($hostinfo=~/Processor type: (.*)/) { $cpuinfo{'model name'}=$1; } my $sysctl=`sysctl -a`; + if($sysctl=~/hw.cpufrequency: (\S+)/) { $cpuinfo{'CPUFREQ'}=int($1/1000.0); @@ -664,10 +644,8 @@ sub get_proccpu while(<PROC>) { chomp; - if(/processor\s+\:/) - { - $nprocessors++; - } + $nprocessors++ if(/processor\s+\:/o); + if(/^(.*)\s+\:\s+(.*)$/) { my $k=$1; @@ -716,7 +694,7 @@ sub get_cpu { if($ipod==1) { - return('ipod'); + return 'ipod'; } else { @@ -725,26 +703,26 @@ sub get_cpu # intel chips if($cpuinfo{'model name'}=~/Pentium\(R\) M/) { - return('pentium-m'); + return 'pentium-m'; } elsif($cpuinfo{'cpu family'}>=15) { - return('pentium4'); + return 'pentium4'; } elsif($cpuinfo{'model name'}=~/Intel/) { if($cpuinfo{'model name'}=~/Core/) { - return('core2'); + return 'core2'; } else { - return('i686'); + return 'i686'; } } else { - return('i686'); + return 'i686'; } } elsif($cpuinfo{'model name'}=~/amd/i) @@ -752,34 +730,34 @@ sub get_cpu # amd chips if($cpuinfo{'model name'}=~/phenom\(tm\)/i) { - return('amdfam10'); + return 'amdfam10'; } if($cpuinfo{'model name'}=~/opteron/i) { - return('opteron'); + return 'opteron'; } elsif($cpuinfo{'model name'}=~/Athlon\(tm\) 64/) { - return('athlon64'); + return 'athlon64'; } elsif($cpuinfo{'model'}>=8) { - return('athlon-xp'); + return 'athlon-xp'; } elsif($cpuinfo{'model'}>=4) { - return('athlon-tbird'); + return 'athlon-tbird'; } else { # this is the default for modern machines - return('athlon'); + return 'athlon'; } } elsif($cpuinfo{'cpu'} =~/750/) { # ibook - return('750'); + return '750'; } elsif($os eq 'Darwin' && ($cpuinfo{'model name'}=~/ppc(\d+)/)) { @@ -788,15 +766,14 @@ sub get_cpu } else { - return('native'); + return 'native'; } } } sub get_arch { - my $arch=get_cpu(); - return($arch); + return get_cpu(); } sub parse_args @@ -849,6 +826,10 @@ sub parse_args { $debug=1; } + elsif($arg eq 'condor') + { + $condor=1; + } elsif($arg eq 'cgi') { $cgi=1; @@ -881,18 +862,14 @@ sub create_master_makefile if($ld) { print FP "LD := $ld\n"; - if($ldpath) - { - print FP "LD_LIBRARY_PATH := $ldpath\n"; - } - if($ldflags) - { - print FP "LDFLAGS := $ldflags\n"; - } + print FP "LD_LIBRARY_PATH := $ldpath\n"if($ldpath); + print FP "LDFLAGS := $ldflags\n"if($ldflags); } print FP "PROGRAM := $binary MAKE := $make +SVN_REVISION := \$(shell svn info|grep Revision|cut -d' ' -f2 ) +SVN_URL := \$(shell svn info|grep URL|cut -d' ' -f2|sed -e 's/^http\\:\\/\\///g') CFLAGS := $cflags $coptflags LIBDIRS := -L$lib -L. $static $acmlflags INCDIRS := -I. -I/usr/include $extrainc @@ -927,20 +904,12 @@ sub gcc_version } else { - return(''); + return ''; } } sub icc_version { - my $v=`icc --version`; - if($v=~/icc \(ICC\) (\d+\.\d+)/) - { - return $1; - } - else - { - return undef; - } + return `icc --version`=~/icc \(ICC\) (\d+\.\d+)/ ? $1 : undef; } sub get_cfiles @@ -970,9 +939,6 @@ sub get_cfiles } push(@c_files,sort{$a cmp $b}(keys(%c))); - # output - #map{print $_,"\n"} sort {$a<=>$b} @c_files; - chdir $owd; return(@c_files); } @@ -1041,66 +1007,6 @@ sub slurp return (do { local( $/ ) ; <$fh> } ); } -sub make_unrolled_loop_functions -{ - - # specialist hard-coded functions - - # get ISOTOPE_ARRAY_SIZE - my $n; - - my @x=grep {/^\#define ISOTOPE_ARRAY_SIZE/} `cd $src_dir;gcc -g3 -E -DNUCSYN nucsyn/nucsyn_isotopes.h 2>/dev/null`; - if($x[0]=~/\#define ISOTOPE_ARRAY_SIZE (\d+)/) - { - $n=$1; - } - else - { - print "Could not generate unrolled loop functions (ISOTOPE_ARRAY_SIZE undetermined)\n"; - return; - } - open(C,">$src_dir/nucsyn/nucsyn_unrolled_functions.h")|| - die("cannot open $src_dir/nucsyn/nucsyn_unrolled_function.h"); - - print C "\#ifdef NUCSYN\n"; - - # for nanchceks - print C "\#define NUCSYN_UNROLLED_NANCHECK if("; - for(my $i=0;$i<$n-1;$i++) - { - print C " (isnan(x[$i])!=0) || "; - } - print C "(isnan(x[",$n-2,"]!=0))) nancheck_oops(label,x);\n"; - - # for nucsyn_check_abundance_array - print C "\n\n#define NUCSYN_UNROLLED_LESS_THAN_ZERO_OR_NANCHECK {"; - for(my $i=0;$i<$n;$i++) - { - print C "if((X[$i]<0.0)||(isnan(X[$i])));"; - } - print C "}\n"; - - # for nucsyn_totalX - print C "\n\n#define NUCSYN_UNROLLED_TOTALX "; - for(my $i=0;$i<$n-1;$i++) - { - print C "X[$i]+"; - } - print C "X[",$n-1,"]\n\n"; - - # for nucsyn_molecular_weight - print C "\n\n#define NUCSYN_UNROLLED_MOLECULAR_WEIGHT "; - for(my $i=0;$i<$n;$i++) - { - print C "X[$i]*Y[$i]+"; - } - print C "X[",$n-1,"]*Y[",$n-1,"]\n\n"; - - print C "#endif // NUCSYN\n"; - - close C; -} - sub find_make { foreach my $f ('/usr/bin/make', @@ -1111,6 +1017,6 @@ sub find_make { return $f if(-s $f > 0); } - print "Could not locate make or gmake\n"; + print "Could not locate make or gmake : please check that GNU make is installed on your system\n"; exit; } diff --git a/doc/binary_grid-flexigrid.lyx b/doc/binary_grid-flexigrid.lyx index 690071602..8d70b04eb 100644 --- a/doc/binary_grid-flexigrid.lyx +++ b/doc/binary_grid-flexigrid.lyx @@ -4961,7 +4961,7 @@ name "sub:setup grid script" \end_layout \begin_layout Standard -You need to do two things: +You need to do a few things: \end_layout \begin_layout Enumerate @@ -4974,6 +4974,60 @@ hash by calling grid_defaults \emph default . + +\begin_inset Newline newline +\end_inset + +Note that you +\emph on +must +\emph default + link a hash in your script to +\family typewriter +$binary_grid::grid_options{results_hash} +\family default + e.g. + +\begin_inset Newline newline +\end_inset + + +\family typewriter +$binary_grid::grid_options{results_hash}=$results; +\family default + (if you have a global anonymous hash e.g. + +\family typewriter +$results={}; +\family default +) +\begin_inset Newline newline +\end_inset + +or +\begin_inset Newline newline +\end_inset + + +\family typewriter +$binary_grid::grid_options{results_hash}= +\backslash +%results; +\family default + (if you have a global hash e.g. + +\family typewriter +%results; +\family default +) +\begin_inset Newline newline +\end_inset + + +\begin_inset Newline newline +\end_inset + + \end_layout \begin_layout Enumerate @@ -7147,6 +7201,270 @@ Only the innermost loop contains the calls to push commands onto the thread queue. \end_layout +\begin_layout Subsection +Flexigrid and Condor +\end_layout + +\begin_layout Standard +\begin_inset CommandInset label +LatexCommand label +name "sub:Condor" + +\end_inset + +Warning: this is work in progress and subject to change! +\begin_inset Newline newline +\end_inset + + +\begin_inset Newline newline +\end_inset + +Flexigrid now supports the Condor distributed computing engine. + To use it, a few minor changes need to be made to your grid script. + Please see +\family typewriter +src/perl/scripts-flexigrid/grid-condor.pl +\family default + for an example. +\end_layout + +\begin_layout Itemize +Add a line (near the top of your Perl script) to use the condor module: +\begin_inset Newline newline +\end_inset + + +\family typewriter +use binary_grid::condor; +\end_layout + +\begin_layout Itemize +Replace the call to +\family typewriter +flexigrid(\SpecialChar \ldots{} +); +\family default + with +\family typewriter +condor_grid(); +\end_layout + +\begin_layout Itemize +At the beginning of your output() subroutine, add a line +\begin_inset Newline newline +\end_inset + + +\family typewriter +return if($binary_grid::grid_options{condor_command} eq 'run_flexigrid'); +\family default + +\begin_inset Newline newline +\end_inset + +(This suppresses output except for the final running job which has all the + data) +\end_layout + +\begin_layout Itemize +Put the following in your +\family typewriter +defaults() +\family default + subroutine: +\begin_inset Newline newline +\end_inset + + +\family typewriter +$binary_grid::grid_options{condor_dir}='some world readable directory'; +\family default + +\begin_inset Newline newline +\end_inset + +This sets the working directory for this Condor job: it should be readable + +\emph on +from every Condor machine! +\emph default + Usually this is an NFS mounted disk. + I have not tried to make Condor work without a universally-writable directory. +\end_layout + +\begin_layout Itemize +Every job you run should use a different output directory. +\end_layout + +\begin_layout Subsubsection* +Condor issues and FAQ +\end_layout + +\begin_layout Standard +(Some of these problems are AIfA specific, but may apply to you) +\end_layout + +\begin_layout Itemize +What are the requirements for the output directory? +\begin_inset Newline newline +\end_inset + +The output directory must be in a world-writable location, e.g. + an NFS mount (see the problem with +\family typewriter + /export +\family default + below). + Each time you run a grid script with a call to condor_grid, a new directory + must be used. + +\emph on +It is up to you to ensure this is the case! +\end_layout + +\begin_layout Itemize +How are the jobs numbered? +\begin_inset Newline newline +\end_inset + +The jobs are numbered +\emph on +x.y +\emph default + (where +\emph on +x +\emph default +is the job number, +\emph on +y +\emph default + is the total number of jobs). + These numbers just have to be unique, you are free to change them. +\end_layout + +\begin_layout Itemize +What goes into the output directory? +\begin_inset Newline newline +\end_inset + +Each job's output ( +\family typewriter +stdout +\family default + and +\family typewriter +stderr +\family default +) goes into the directories +\emph on +stdout +\emph default + and +\emph on +stderr +\emph default +. + +\begin_inset Newline newline +\end_inset + +Each job also has a log which goes into the +\emph on +log +\emph default + directory. + +\end_layout + +\begin_layout Itemize +How do I put extra parameters in the Condor submit scripts? +\begin_inset Newline newline +\end_inset + +You can use a hash called +\family typewriter +%binary_grid::grid_options{condor_options} +\family default + e.g. +\begin_inset Newline newline +\end_inset + + +\family typewriter +$binary_grid::grid_options{condor_options}{Requirements}='Machine!= +\begin_inset Quotes erd +\end_inset + +aibn73.astro.uni-bonn.de +\begin_inset Quotes erd +\end_inset + +'; +\end_layout + +\begin_layout Itemize +I have problems running my script with errors +\family typewriter +/export/\SpecialChar \ldots{} + +\family default + not found +\begin_inset Newline newline +\end_inset + +The problem here is the +\family typewriter +/export +\family default + is local to your PC: you should run your script from +\family typewriter +/vol/\SpecialChar \ldots{} + +\family default + (a world-readable NFS mount point) instead. + There is a function in the example condor scripts called +\family typewriter +fix_aifa_environment( +\family default +) ... + try using it. +\end_layout + +\begin_layout Itemize +Condor jobs are using all my CPU! +\begin_inset Newline newline +\end_inset + +They should be niced (nice -n 19) and preferably ioniced (ionice -c3) and + given SCHED_IDLE priority ( +\family typewriter +chrt -i -p 0 +\emph on +<pid> +\family default +\emph default +) as well. + The +\emph on +rob_misc +\emph default + Perl module is used to do this automatically (with a call to +\family typewriter +renice_me(); +\family default +) There are issues with the +\emph on +Linux +\emph default + kernel and +\emph on +cgroups +\emph default + which are pending. +\end_layout + \begin_layout Standard \begin_inset Newpage pagebreak \end_inset @@ -11791,6 +12109,22 @@ binary_c \begin_layout Description +\family typewriter +arg_checking +\family default + Sometimes, an argument to binary_c will not be recognised (e.g. + if you set something in a bse_option which is not compatible with your + binary_c, perhaps because you haven't enabled some feature in binary_c) + and the grid will give errors (and/or crash) which are rather cryptic. + Instead of this, you can set arg_checking to 1 to have strict checking + of each argument that is passed to binary_c. + This is, however, rather slow (because it requires a lot of I/O) so the + default is 0. + +\end_layout + +\begin_layout Description + \family typewriter args \family default @@ -11814,6 +12148,29 @@ binary \begin_layout Description +\family typewriter +cache_binary_c_output +\family default + If this is 1 (the default) then output from binary_c is cached in a local + array before being processed through +\family typewriter +tbse_line() +\family default + calls. + This should be faster, because it means less switching from one process + to another, but its use means that +\emph on +binary_c +\emph default + is out of sync with +\family typewriter +tbse_line() +\family default +, so debugging is more difficult. +\end_layout + +\begin_layout Description + \family typewriter colour \family default diff --git a/grid-IMF.pl b/grid-IMF.pl deleted file mode 120000 index 0bd133b93..000000000 --- a/grid-IMF.pl +++ /dev/null @@ -1 +0,0 @@ -src/perl/scripts/grid-IMF.pl \ No newline at end of file diff --git a/grid-wrap.pl b/grid-wrap.pl new file mode 100755 index 000000000..5148a0f46 --- /dev/null +++ b/grid-wrap.pl @@ -0,0 +1,372 @@ +#!/usr/bin/env perl +use 5.16.0; +use strict; +use Sort::Key qw(nsort); +use rob_misc; +use binary_grid; +use threads; +use threads::shared; +use Fcntl ':flock'; # import LOCK_* constants +use subs qw(logprint logprintf); +use Carp qw(confess); +use File::Path qw(make_path); +$|=1; +############################################################ + +my $outdir='/tmp/gridtest'; +my $nmachines=20; + +say "Run across (modulo=) ",scalar(@machines)," machines"; +my $modulo=scalar(@machines); + +my @pids; +my @datafiles; +my $offset=0; +my @spectral_type_strings=('O','B','A','F','G','K','M'); +my %spectral_type; revspechash(); +my @stellar_type_strings; +my @short_stellar_type_strings; +my %short_stellar_type_strings; +my $period_binwidth=0.1; # width in logP of histograms +my $log_fp; +defaults(); + +foreach my $machine (0..$nmachines) +{ + my $datafile="$outdir/results.$offset"; + my $cmd="./src/perl/scripts-flexigrid/grid-canberra.pl modulo=$modulo offset=$offset results_hash_dumpfile=$datafile"; + + say "Run $cmd"; + + my $forkr; + + # run the models : should outsource this to binary_grid::condor module + if(1){ + $forkr=fork(); + if($forkr==0) + { + print `$cmd 2>$outdir/grid.$offset.err > $outdir/grid.$offset.out`; + exit; + } + } + + push(@pids,$forkr); + push(@datafiles,$datafile); + + $offset++; +} + +# wait for children to finish +say "Waiting for children (pids: @pids)"; +map { waitpid $_,0; }@pids; + +say "Children finished\n"; + + +# put all the data together +my $results=merge_results_hash_dumps(@datafiles); + +# save ouput +output(); + +if(0) +{ + use Data::Dumper; + open($log_fp,">$outdir/datadump")||confess("cannot open $outdir/datadump"); + print {$log_fp} Dumper($results); + close $log_fp; +} + +say "Done\n"; +exit; + +############################################################ + + +sub output +{ + # open log file + + # check for an output directory, otherwise assume a default in /tmp + + # output binary_grid::grid_options/bse_options + open($log_fp,">$outdir/grid_options")||confess("cannot open $outdir/grid_options for output"); + print $log_fp "grid_options\n\n"; + map + { + print {$log_fp} $_.' : '.$binary_grid::grid_options{$_}."\n"; + }keys %binary_grid::grid_options; + close $log_fp; + open($log_fp,">$outdir/bse_options")||confess("cannot open $outdir/bse_options for output"); + print $log_fp "bse_options\n\n"; + map + { + print {$log_fp} $_.' : '.$binary_grid::bse_options{$_}."\n"; + }keys %binary_grid::bse_options; + close $log_fp; + + # output grid results + open ($log_fp,">$outdir/log") || confess ("cannot open $outdir/log for output"); + logprint "\n\n\nCANBERRA Grid results\n============================================================\n\n"; + logprint "SVN Revision ",`svn info|grep Revision|gawk \"\{print \\\$2\}\"`; + logprint "Run at ",scalar(localtime),"\n"; + + logprint "\n\nParameters:\n"; + logprint "Initial binary fraction : Raghavan 2010\n"; + foreach my $p ('alpha_ce','lambda_ce','lambda_ionisation','z','wd_sigma') + { + my $opt = ($p eq 'lambda_ce' && $binary_grid::bse_options{$p}==-1) + ? 'Fitted to detailed stellar models' + : + sprintf "%g",$binary_grid::bse_options{$p}; + + logprint sprintf "% -30s = %s\n",$p eq'z'?'Metallicity':$p,$opt; + } + + + my $dtp_wd=1e-30; # total number count of WDs (single or binary) + my %dtp_wd; + + { + my $k='number count'; + # get total WD number count + foreach my $duplicity (keys %{$results->{$k}}) + { + foreach my $wd_type (keys %{$results->{$k}->{$duplicity}}) + { + $dtp_wd += $$results{$k}->{$duplicity}->{$wd_type}; + $dtp_wd{$duplicity} += $$results{$k}->{$duplicity}->{$wd_type}; + } + } + } + + my $dtp_wdms=0.0; + { + # get total WD-MS number count + my $k='MSWD spectral types'; + foreach my $spectral_type (nsort keys %{$results->{$k}}) + { + $dtp_wdms += $$results{$k}{$spectral_type}; + } + } + + + foreach my $k (sort keys %$results) + { + next if($k eq 'thread_number'); + if($k eq 'number count') + { + logprint sprintf "WD binary fraction %3.2f %%\n",$dtp_wd{'Binary'}/$dtp_wd*100.0; + foreach my $duplicity (keys %{$results->{$k}}) + { + foreach my $wd_type (nsort keys %{$results->{$k}->{$duplicity}}) + { + logprint sprintf "% -6s % -30s %3.2f %%\n", + $duplicity, + $stellar_type_strings[$wd_type], + $$results{$k}->{$duplicity}->{$wd_type}/$dtp_wd*100.0; + } + } + } + elsif($k eq 'companion types') + { + my $test_tot=0.0; + logprint "\n\nWhite Dwarf - ? pairs (as a fraction of all White Dwarf Binaries)\n"; + foreach my $wd1 (nsort keys %{$results->{$k}}) + { + foreach my $wd2 (nsort keys %{$results->{$k}->{$wd1}}) + { + logprint sprintf "% -40s - % -40s % 5.2f %%\n", + $stellar_type_strings[$wd1], + $stellar_type_strings[$wd2], + $$results{$k}{$wd1}{$wd2} + /$dtp_wd{'Binary'}*100.0; + } + } + logprint "\n"; + } + elsif($k eq 'test') + { + # test initial binary fraction + my $test_tot=0.0; + logprint "\n\nTest numbers\n\n"; + foreach my $duplicity (keys %{$results->{$k}}) + { + $test_tot += $$results{$k}->{$duplicity}; + } + foreach my $duplicity (nsort keys %{$results->{$k}}) + { + logprint sprintf "MS $duplicity %3.2f %%\n", 100.0*$$results{$k}->{$duplicity}/$test_tot; + } + } + elsif($k eq 'test dist spectypes') + { + print "ZAMS Binary fraction by spectral type : ",keys %{$results->{$k}->{'Binary'}},"\n\n"; + foreach my $spec_type (@spectral_type_strings) + { + my $t=$$results{$k}->{'Single'}->{$spec_type}+$$results{$k}->{'Binary'}->{$spec_type}; + logprint sprintf "%s %5.2f %%\n",$spec_type,$$results{$k}->{'Binary'}->{$spec_type}/$t*100.0 if($t>0.0); + } + } + elsif($k eq 'test dist') + { + # test initial distributions + foreach my $duplicity (keys %{$results->{$k}}) + { + foreach my $dist (keys %{$results->{$k}->{$duplicity}}) + { + my $f="$outdir/init_".$duplicity.'_'.$dist.'.dat'; + logprint "Output distribution to $f\n"; + open(FP,">$f")||confess("can't open $f"); + foreach my $x (nsort keys %{$results->{$k}->{$duplicity}->{$dist}}) + { + printf FP "%g %g\n",$x,$$results{$k}{$duplicity}{$dist}{$x}; + } + close FP; + } + } + } + + elsif($k eq 'MSWD spectral types') + { + my $f="$outdir/MSWD_spectral_types.dat"; + open(FP,'>'.$f)||confess("cannot open $f for output"); + logprint "MSWD binaries: spectral types (see datafile $f)\n"; + pad_datahash($results->{$k},1); + foreach my $spectral_type (nsort keys %{$results->{$k}}) + { + printf FP "%d %g\n",$spectral_type,$$results{$k}{$spectral_type}/$dtp_wdms; + logprint sprintf "MSWD spectral type %d = %s = %g\n", + $spectral_type,$spectral_type_strings[$spectral_type],$$results{$k}{$spectral_type}/$dtp_wdms; + } + logprint "\n"; + close FP; + } + # misc period distributions + elsif($k =~/(\S+) period distribution/) + { + my $distname=$1; + my $f="$outdir/$distname".'_period_distribution.dat'; + #print "Pad datahash $results->{$k}\n"; + pad_datahash($results->{$k},$period_binwidth); + open(FP,'>'.$f)||confess("cannot open $f for output"); + logprint "$1 binaries: period distribution (see datafile $f)\n"; + foreach my $logP (nsort keys %{$results->{$k}}) + { + printf FP "%g %g\n",$logP,$$results{$k}{$logP}; + } + close FP; + } + } + + # compare to Holberg 2009 + { + my %obs=( + # Holberg 2009 + 'Single WDs'=>0.70, + 'WD +dM'=>0.12, + 'Sirius-Like'=>0.08, + 'WD + WD'=>0.10, + 'Total WDs'=>1.0, + ); + map{$obs{$_}*=100.0}keys %obs; + + my $format="\% -30s \% 5.2f \%\% (% 5.2f)\n"; + logprint "\n=============\nHolberg 2009 comparison\n\n"; + + # Single WD fraction + my $single_WD = $dtp_wd{'Single'}/$dtp_wd*100.0; + logprintf $format,'Single WDs',$single_WD,$obs{'Single WDs'}; + + # WD + dM fraction + my $WD_dM=$$results{'MSWD spectral types'}{$spectral_type{'M'}}/$dtp_wd*100.0; + logprintf $format,'WD +dM',$WD_dM,$obs{'WD +dM'}; + + # Sirius-like fraction + my $Sirius=0.0; + map{$Sirius+=$$results{'MSWD spectral types'}{$spectral_type{$_}}}('O','B','A','F','G','K'); + $Sirius = $Sirius/$dtp_wd * 100.0; + logprintf $format,'Sirius-Like',$Sirius,$obs{'Sirius-Like'}; + + # WD + WD fraction + my $WDWD=0.0; + foreach my $wd1 (10,11,12) + { + foreach my $wd2 (10,11,12) + { + $WDWD += $$results{'companion types'}{$wd1}{$wd2}; + } + } + $WDWD = $WDWD/$dtp_wd*100.0; + logprintf $format,'WD + WD',$WDWD,$obs{'WD + WD'}; + + # total (should be 100%!) + logprintf $format,'Total WDs',($single_WD + $WD_dM + $Sirius + $WDWD),$obs{'Total WDs'}; + } + + logprint "\n============================================================\n\n\n"; + + close $log_fp; +} + +sub logprint +{ + print {$log_fp} "@_"; + print "@_"; +} +sub logprintf +{ + my $format=shift @_; + printf {$log_fp} $format,@_; + printf $format,@_; +} + +sub revspechash +{ + for(my $i=0;$i<=$#spectral_type_strings;$i++) + { + $spectral_type{$spectral_type_strings[$i]}=$i; + } +} + +sub spec_to_float +{ + my $s=$_[0]; # spectral type to check + if($s=~/([A-Z])(\d)/) + { + return $spectral_type{$1}+0.1*$2; + } + else + { + print "Failed to convert spec type $_[1] to float\n"; + exit; + } +} + +sub spec_type_sorter +{ + return spec_to_float($a) <=> spec_to_float($b); +} + +sub pad_datahash +{ + # pad datahash $h so it has zeros where there's no data + my $h=shift; + my $d=shift; + my @k=nsort keys %$h; + my $xmax=$k[$#k]; # max value + my $x=$k[0]; # start (min) value + while($x<$xmax) + { + $$h{$x}=0.0 if(!defined($$h{$x})); + $x=sprintf '%g',$x+$d; + } +} + +sub defaults +{ + grid_defaults(); + @stellar_type_strings=@{(stellar_type_strings())[2]}; + @short_stellar_type_strings=@{(stellar_type_strings())[3]}; + %short_stellar_type_strings=%{(stellar_type_strings())[4]}; +} diff --git a/grid.pl b/grid.pl deleted file mode 120000 index 2a6f488b6..000000000 --- a/grid.pl +++ /dev/null @@ -1 +0,0 @@ -src/perl/scripts/grid.pl \ No newline at end of file diff --git a/imftickfile b/imftickfile deleted file mode 120000 index f27a91d62..000000000 --- a/imftickfile +++ /dev/null @@ -1 +0,0 @@ -./imfticks/imf_ticks-2GYr \ No newline at end of file diff --git a/src/binary_parameters.h b/src/binary_parameters.h index 851d48c11..c65c7c494 100644 --- a/src/binary_parameters.h +++ b/src/binary_parameters.h @@ -21,14 +21,13 @@ * except those which contain ONLY SETTINGS */ +#include <stdio.h> /* * code_options.h contains the machine-specific and build-specific coding * options, rather than options which should affect the behaviour of the code. * It should be read FIRST. */ -#include "code_options.h" -#include <stdio.h> #include <signal.h> #include <stdlib.h> #include <string.h> @@ -37,11 +36,12 @@ #include <fcntl.h> #include <sys/time.h> #include <sys/mman.h> +#include "code_options.h" /* * Make binary_c/nucsyn as similar to Amanda's models as possible */ -#define LIKE_AMANDA +//#define LIKE_AMANDA // Enable Axel's changes : WILL FAIL //#define AXEL @@ -57,7 +57,7 @@ * Alternative, simpler form from rob: use this instead of SMOOTH_AGB? * You will need ROBSMOOTHLOG to use testsmooth.pl */ -//#define ROBSMOOTH +#define ROBSMOOTH //#define ROBSMOOTHDEBUG //#define ROBSMOOTHLOG //#define ROBSMOOTH_STOP_AT_TPAGB @@ -689,7 +689,7 @@ #define HRDIAG_BIN_SIZE_LOGTIME 0.05 #define HRDIAG_START_LOG 5.0 #define HRDIAG_BIN_SIZE_TEFF 100.0 -#define HRDIAG_BIN_SIZE_LOGL 0.01 +#define HRDIAG_BIN_SIZE_LOGL 0.1 #define HRDIAG_BIN_SIZE_LOGG 0.5 // enable this for more time resolution at the AGB->WD transition @@ -700,6 +700,7 @@ #undef SINGLE_STAR_LIFETIMES #define SINGLE_STAR_LIFETIMES +#undef FILE_LOG #endif//HRDIAG @@ -1272,6 +1273,8 @@ #undef ROTATIONALLY_ENHANCED_MASS_LOSS #undef ADAPTIVE_RLOF #undef SELMA_MASS_EJECTA +#undef NANCHECKS +#undef NUCSYN #endif /* diff --git a/src/binary_prototypes.h b/src/binary_prototypes.h index 38ca8e9f2..7d2ce71fc 100644 --- a/src/binary_prototypes.h +++ b/src/binary_prototypes.h @@ -6,12 +6,13 @@ #ifndef BINARY_PROTOTYPES_H #define BINARY_PROTOTYPES_H +#include <stdio.h> +#include <stdlib.h> #include "code_options.h" #include "binary_parameters.h" #include "binary_macros.h" #include "binary_structures.h" -#include <stdlib.h> /* * Get the prototypes for each of the libraries (possibly diff --git a/src/binary_star_functions/common_envelope_evolution.c b/src/binary_star_functions/common_envelope_evolution.c index 4c355c64e..75792bc9b 100644 --- a/src/binary_star_functions/common_envelope_evolution.c +++ b/src/binary_star_functions/common_envelope_evolution.c @@ -134,7 +134,7 @@ int common_envelope_evolution (double *m01, #ifdef COMPARE_FORTRAN printf("IN COMENV\n"); #endif - Dprint("In comenv\n");fflush(stdout); + Dprint("In comenv\n"); #ifdef STELLAR_TYPE_LOG printf("STELLARTYPES COMENV\n"); #endif diff --git a/src/binary_structures.h b/src/binary_structures.h index 2564dc176..a201ac9cb 100644 --- a/src/binary_structures.h +++ b/src/binary_structures.h @@ -8,10 +8,10 @@ #ifndef BINARY_STRUCTURES_H #define BINARY_STRUCTURES_H - +#include <stdio.h> #include "binary_parameters.h" /* always put this first! */ #include "binary_macros.h" -#include <stdio.h> + #ifdef NUCSYN /* Grab the settings for the nucleosynthesis */ @@ -356,6 +356,7 @@ struct preferences_t { Boolean batchmode; int batch_submode; char batchmode_string[BATCHMODE_STRING_LENGTH]; + Boolean no_arg_checking; #endif #ifdef AXEL // overshooting TPAGB diff --git a/src/claet/claet_low_mass_GB.c b/src/claet/claet_low_mass_GB.c index bdeebfab6..bc366ca11 100644 --- a/src/claet/claet_low_mass_GB.c +++ b/src/claet/claet_low_mass_GB.c @@ -22,7 +22,7 @@ void claet_low_mass_GB(double GBp, luminosities[L_BGB] = lbgbf(mass,giant_branch_parameters); /* Set GB timescales */ - Dprint("Set GB timescales\n");fflush(stdout); + Dprint("Set GB timescales\n"); double y=iq/(GB[GB_EFFECTIVE_AH]*GB[GB_B]); double y2=ip/(GB[GB_EFFECTIVE_AH]*GB[GB_D]); diff --git a/src/claet/claet_not_a_helium_star.c b/src/claet/claet_not_a_helium_star.c index 0e32da8f2..50c47fd97 100644 --- a/src/claet/claet_not_a_helium_star.c +++ b/src/claet/claet_not_a_helium_star.c @@ -110,8 +110,7 @@ int claet_not_a_helium_star(const double mass, printf("CFSTARe branch mass=%g zpar3=%g\n",mass,metallicity_parameters[3]); #endif - Dprint("compare mass=%g vs metallicity_parameter[3]=%g\n",mass,metallicity_parameters[3]);fflush(stdout); - + Dprint("compare mass=%g vs metallicity_parameter[3]=%g\n",mass,metallicity_parameters[3]); if LESS_OR_EQUAL(mass,metallicity_parameters[3]) { diff --git a/src/claet/claet_post_HG.c b/src/claet/claet_post_HG.c index 975ffbe15..c82fc031f 100644 --- a/src/claet/claet_post_HG.c +++ b/src/claet/claet_post_HG.c @@ -29,7 +29,7 @@ int claet_post_HG(double tbagb, /* * Set the core mass at the BGB. */ - Dprint("Set BGB core mass\n");fflush(stdout); + Dprint("Set BGB core mass\n"); if LESS_OR_EQUAL(mass,metallicity_parameters[2]) { @@ -58,7 +58,7 @@ int claet_post_HG(double tbagb, /* * FAGB time parameters */ - Dprint("Set FAGB parameters\n");fflush(stdout); + Dprint("Set FAGB parameters\n"); ilmx=1.0/luminosities[L_LMX]; ilbagb=1.0/luminosities[L_BAGB]; @@ -77,7 +77,7 @@ int claet_post_HG(double tbagb, timescales[7],timescales[8],timescales[9]); #endif - Dprint("Find LTP etc\n");fflush(stdout); + Dprint("Find LTP etc\n"); /* * Now to find Ltp and ttp using Mc,He,tp diff --git a/src/code_options.h b/src/code_options.h index 90cd2b3a0..3692437a2 100644 --- a/src/code_options.h +++ b/src/code_options.h @@ -306,7 +306,7 @@ * Enable to generate extra code for NaN checking * WARNING: very slow but sometimes necessary */ -#define NANCHECKS +//#define NANCHECKS /* * Perhaps unroll loops? (speeds up nanchecking!) @@ -364,4 +364,8 @@ #define printf(...) printf(__VA_ARGS__) //binary_c_printf(__VA_ARGS__) + + +//#define fflush(...) binary_c_flush(__VA_ARGS__) + #endif //CODE_OPTIONS_H diff --git a/src/free_memory.c b/src/free_memory.c index f548c6955..c49f2e1cc 100644 --- a/src/free_memory.c +++ b/src/free_memory.c @@ -14,7 +14,7 @@ void free_memory(struct stardata_t *stardata) #endif /* MEMMAP */ #ifdef DEBUG - Dprint("Freeing memory\n");fflush(stdout);fflush(stderr); + Dprint("Freeing memory\n"); #endif int kt; for(kt=0;kt<KTYPE_ARRAY_SIZE;kt++) diff --git a/src/hrdiag/hrdiag.c b/src/hrdiag/hrdiag.c index 9a41b6952..2654bff98 100644 --- a/src/hrdiag/hrdiag.c +++ b/src/hrdiag/hrdiag.c @@ -549,9 +549,11 @@ void hrdiag(double * mass, /** star mass in solar units**/ { fprintf(stdout,"OOPS GB to HG in one step!\n"); fprintf(stderr,"OOPS GB to HG in one step!\n"); - Dprint("OOPS GB to HG in one step!\n"); + Dprint("OOPS GB to HG in one step!\n"); +#ifdef FILE_LOG if(stardata->model.log_fp!=NULL) fprintf(stardata->model.log_fp,"OOPS GB to HG in one step!\n"); +#endif exit(0); } @@ -563,11 +565,12 @@ void hrdiag(double * mass, /** star mass in solar units**/ stardata->model.dtm,stardata->model.dt,stardata->model.intpol); Dprint("OOPS CHeB to HG in one step! dtm=%g dt=%g intpol=%dx\n", stardata->model.dtm,stardata->model.dt,stardata->model.intpol); +#ifdef FILE_LOG if(stardata->model.log_fp!=NULL) fprintf(stardata->model.log_fp,"OOPS CHeB to HG in one step! dtm=%g dt=%g intpol=%dx\n", stardata->model.dtm,stardata->model.dt,stardata->model.intpol ); - +#endif exit(0); } } diff --git a/src/hrdiag/hrdiag_HeStar.c b/src/hrdiag/hrdiag_HeStar.c index ef168cdae..3c6351627 100644 --- a/src/hrdiag/hrdiag_HeStar.c +++ b/src/hrdiag/hrdiag_HeStar.c @@ -86,7 +86,7 @@ void hrdiag_HeStar(Stellar_type *stellar_type, *r = *rx*(1.0+am*(tau-POW6(tau))); *rg = *rx; -#ifdef NUCSYN_WR +#if (defined(NUCSYN)&&defined(NUCSYN_WR)) if(!RLOF_CALL) { /* @@ -171,7 +171,7 @@ void hrdiag_HeStar(Stellar_type *stellar_type, Dprint("core_mass=%g, mt=%g, mtc = %g, mcmax=%g (from mass=%g)\n", *core_mass,*mt,*mtc,mcmax,*mass); -#ifdef NUCSYN_WR +#if (defined(NUCSYN)&&defined(NUCSYN_WR)) if(!RLOF_CALL) star->he_mcmax=mcmax; // save for use in nucsyn_WR #endif diff --git a/src/iteration/end_of_iteration.c b/src/iteration/end_of_iteration.c index 102de298d..f6d8da00f 100644 --- a/src/iteration/end_of_iteration.c +++ b/src/iteration/end_of_iteration.c @@ -10,22 +10,15 @@ void end_of_iteration(struct stardata_t *stardata) { -#ifdef FILE_LOG - if(stardata->model.log_fp==NULL) - { - printf("NULL stardata ? fp=%p stardata=%p (when it should be defined!)\n",stardata->model.log_fp,stardata); - return; - } -#endif if(stardata->model.com==TRUE) - { + { stardata->model.com = FALSE; - } + } #if (defined(LOGGING)&&defined(FILE_LOG)) else { - if(stardata->model.coel==TRUE) + if(stardata->model.coel==TRUE) { output_to_logfile(stardata->model.log_fp, stardata->model.model_time, @@ -38,7 +31,7 @@ void end_of_iteration(struct stardata_t *stardata) 0.0, -1.0, stardata->preferences->label[COELESCE_LABEL],stardata -); + ); } else if(BOTH_MASSLESS) @@ -79,8 +72,8 @@ void end_of_iteration(struct stardata_t *stardata) char tempstring[STRING_LENGTH]; snprintf(tempstring,STRING_LENGTH, - "Probability : %.10g\n", - stardata->model.probability); + "Probability : %.10g\n", + stardata->model.probability); output_string_to_log(stardata, tempstring); } @@ -92,7 +85,7 @@ void end_of_iteration(struct stardata_t *stardata) if((stardata->star[1].stellar_type==11)&& (stardata->star[2].stellar_type<10)&& (stardata->model.sgl==0)) - { + { double *N=malloc(ISOTOPE_MEMSIZE); X_to_N(stardata->common.imnuc, 10.0,N,stardata->star[2].Xsurf, @@ -118,10 +111,10 @@ void end_of_iteration(struct stardata_t *stardata) XF19,XFe56), nucsyn_square_multibracket(stardata->star[2].Xsurf, stardata->common.Xsolar, - (Isotope[]){XC12,XC13,XN14,XN15},4,(Isotope[]){XFe56},1) - ); + (Isotope[]){XC12,XC13,XN14,XN15},4,(Isotope[]){XFe56},1) + ); - } + } #endif iterate_logging(stardata); diff --git a/src/iteration/init_star.c b/src/iteration/init_star.c index b9e5e1634..c0dbc451f 100644 --- a/src/iteration/init_star.c +++ b/src/iteration/init_star.c @@ -283,5 +283,5 @@ void init_star(struct star_t * RESTRICT star, star->stellar_type=1; star->last_stellar_type=1; - Dprint("Done for this star\n");fflush(stdout); + Dprint("Done for this star\n"); } diff --git a/src/iteration/initialize_parameters.c b/src/iteration/initialize_parameters.c index 0c44fc429..53b78d56d 100644 --- a/src/iteration/initialize_parameters.c +++ b/src/iteration/initialize_parameters.c @@ -195,13 +195,17 @@ void initialize_parameters(struct stardata_t * RESTRICT stardata) #ifdef IDUM_LOG printf("IDUM=%ld\n",common->idum); +#ifdef FILE_LOG if(model->log_fp!=NULL) { fprintf(model->log_fp,"IDUM=%ld\n",common->idum); + } +#endif + #ifdef IDUM_LOG_STREAM - fprintf(IDUM_LOG_STREAM,"IDUM=%ld\n",common->idum); + if(IDUM_LOG_STREAM!=NULL) fprintf(IDUM_LOG_STREAM,"IDUM=%ld\n",common->idum); #endif - } + #endif model->kmin = 1; diff --git a/src/iteration/iterate_logging.c b/src/iteration/iterate_logging.c index 478a77def..ced4408f5 100644 --- a/src/iteration/iterate_logging.c +++ b/src/iteration/iterate_logging.c @@ -31,6 +31,15 @@ void iterate_logging(struct stardata_t * stardata) /* single star condition */ Boolean single = system_is_single(stardata); +//#define MIMF +#ifdef MIMF + if(stardata->star[1].stellar_type==TPAGB) + { + printf("MCWD %g %g\n", + stardata->star[1].pms_mass, + stardata->star[1].core_mass); + } +#endif #ifdef TPAGBTRACKS if(first_timestep) @@ -365,8 +374,6 @@ void iterate_logging(struct stardata_t * stardata) X[XO18] ); - fflush(stdout); - nucsyn_burn(T,dt,N,&(stardata->star[1]),stardata); N_to_X(stardata->common.mnuc, @@ -2375,8 +2382,6 @@ void iterate_logging(struct stardata_t * stardata) output_to_detailed_compact_logfile(stardata); #endif - fflush(stdout); - } static Boolean system_is_single(struct stardata_t * stardata) @@ -2384,7 +2389,7 @@ static Boolean system_is_single(struct stardata_t * stardata) /* * Is a system single (TRUE) or binary (FALSE) */ - return + Boolean single = // sgl flag stardata->model.sgl || stardata->star[1].stellar_type==MASSLESS_REMNANT || @@ -2396,4 +2401,16 @@ static Boolean system_is_single(struct stardata_t * stardata) stardata->common.separation < 0.0 || stardata->common.eccentricity < 0.0 ; + + Dprint("system_is_single sgl %d : MLR1 %d : MLR2 %d : M2<0.01 %d : P>long %d : a<0 %d : e<0 %d : single ? %d\n", + stardata->model.sgl, + stardata->star[1].stellar_type==MASSLESS_REMNANT , + stardata->star[2].stellar_type==MASSLESS_REMNANT , + stardata->star[2].mass < 0.010001 , + stardata->common.orbital_period > 1e7*YEAR_LENGTH_IN_DAYS , + stardata->common.separation < 0.0 , + stardata->common.eccentricity < 0.0, + single); + + return single; } diff --git a/src/iteration/to_loop_or_not_to_loop.c b/src/iteration/to_loop_or_not_to_loop.c index 0bfd77451..d35941cbc 100644 --- a/src/iteration/to_loop_or_not_to_loop.c +++ b/src/iteration/to_loop_or_not_to_loop.c @@ -14,9 +14,7 @@ Boolean to_loop_or_not_to_loop(struct stardata_t * RESTRICT stardata) /* that is the question! */ { - Dprint("Determine whether to loop again\n");fflush(stdout); - - + Dprint("Determine whether to loop again\n"); if(stardata->model.change==TRUE) { stardata->model.change = FALSE; diff --git a/src/logging/close_log_files.c b/src/logging/close_log_files.c index fb251d4fb..1d90e1f3e 100644 --- a/src/logging/close_log_files.c +++ b/src/logging/close_log_files.c @@ -30,8 +30,6 @@ void close_log_files(FILE *(*log_fp), fflush(NULL); exit(LOG_FILE_CLOSE_ERROR); } - fflush(stderr); - fflush(stdout); *log_fp=NULL; #endif/* FILE_LOG */ diff --git a/src/logging/log_hr.c b/src/logging/log_hr.c index 92c74186d..a36ed74ad 100644 --- a/src/logging/log_hr.c +++ b/src/logging/log_hr.c @@ -28,8 +28,6 @@ static void out_hr_log(double t, double fill, Boolean use2); - - void log_hr(struct stardata_t * RESTRICT stardata) { static double prevlog=0.0; @@ -216,7 +214,6 @@ void log_hr(struct stardata_t * RESTRICT stardata) /* add up the probability in nuclear-burning systems */ p_burning = ((star1->stellar_type < HeWD)||((use2)&&(star2->stellar_type < HeWD))) ? 1.0 : 0.0; - double dtp; if(zams) { logt=0.0; @@ -260,8 +257,8 @@ void log_hr(struct stardata_t * RESTRICT stardata) if(f1>TINY) { - dtp = f1 * dt * stardata->model.probability; #ifdef SAMPLE_DEBUG + double dtp = f1 * dt * stardata->model.probability; printf("dtp1=%g, f1=%g\n",dtp,f1); #endif out_hr_log(stardata->model.model_time, @@ -415,10 +412,12 @@ static void out_hr_log(double t, } printf("\n"); + /* if(star1->stellar_type<10) { printf("HRDEXACT %g %g\n",TEFFSTAR(star1),star1->luminosity); } + */ } } diff --git a/src/logging/open_log_files.c b/src/logging/open_log_files.c index c01a5c55d..16fceb653 100644 --- a/src/logging/open_log_files.c +++ b/src/logging/open_log_files.c @@ -101,7 +101,6 @@ void open_log_files(FILE *(*log_fp), #if (DEBUG==1) fflush(NULL); #endif - fflush(stdout); #endif/*FILE_LOG*/ #endif/*LOGGING*/ diff --git a/src/logging/output_string_to_log.c b/src/logging/output_string_to_log.c index d883044e4..20555b336 100644 --- a/src/logging/output_string_to_log.c +++ b/src/logging/output_string_to_log.c @@ -14,11 +14,7 @@ void output_string_to_log(struct stardata_t * RESTRICT stardata, #ifdef FILE_LOG if(stardata->model.log_fp!=NULL) { - fprintf( - stardata->model.log_fp, - "%s", - c); - + fprintf(stardata->model.log_fp,"%s",c); #ifdef FILE_LOG_FLUSH fflush(stardata->model.log_fp); #endif diff --git a/src/logging/output_to_logfile.c b/src/logging/output_to_logfile.c index af13686ff..459223900 100644 --- a/src/logging/output_to_logfile.c +++ b/src/logging/output_to_logfile.c @@ -42,6 +42,8 @@ void output_to_logfile(FILE *fp, struct stardata_t * RESTRICT stardata ) { + if(fp==NULL) return; // do nothing if no log file is open + #ifdef LOGGING #ifdef FILE_LOG char logstring[LOGSTRINGLENGTH]; diff --git a/src/main.c b/src/main.c index c9e69fd3d..7d0208a7d 100644 --- a/src/main.c +++ b/src/main.c @@ -75,5 +75,7 @@ int main ( int argc, char * * argv ) ,repeat #endif ); + + fflush(stdout); return(NORMAL_EXIT); } diff --git a/src/main.h b/src/main.h index 48514dafa..a571d01a7 100644 --- a/src/main.h +++ b/src/main.h @@ -1,5 +1,5 @@ -#include "code_options.h" + /** header files **/ #include <stdio.h> @@ -13,6 +13,7 @@ #include <sys/time.h> #include <sys/mman.h> #include <time.h> +#include "code_options.h" /* binary_c/nucsyn's header files */ #include "misc/misc.h" diff --git a/src/main2.c b/src/main2.c index 8ef087d14..5893646c1 100644 --- a/src/main2.c +++ b/src/main2.c @@ -21,6 +21,8 @@ #include <sys/time.h> #include <sys/mman.h> #include <time.h> +#include <execinfo.h> + inline void catch_signal(void) { @@ -451,7 +453,7 @@ void batchmode_loop(struct preferences_t * preferences, #ifdef CPUFREQ /* calculate total run time */ double runtime=((float)ticks)/(1e6*((float)CPUFREQ)); - printf("Tick count %ld, ~ runtime = %g",ticks,runtime); + printf("Tick count %ld, ~ runtime = %g\n",ticks,runtime); #else printf("Tick count %ld (CPUFREQ undefined)\n",ticks); #endif // CPUFREQ @@ -462,13 +464,12 @@ void batchmode_loop(struct preferences_t * preferences, first = FALSE; #endif// RDTSC_TIMER_PER_SYSTEM - fflush(stdout); #endif//RDTSC_TIMER printf("fin\n"); - fflush(stdout); + fflush(stdout); // flush is REQUIRED here } } - else + else { while((batchargs[++b]=strtok(NULL," "))&&(b<98)){} /* @@ -481,7 +482,7 @@ void batchmode_loop(struct preferences_t * preferences, */ { parse_arguments(0,b,batchargs,stardata); - fflush(stdout); + //fflush(stdout); } } } @@ -536,3 +537,60 @@ void final_timing_output(ticks start_tick #endif// RDTSC_TIMER } + +static void print_trace (void); +int binary_c_flush(FILE *stream) +{ + /* wrapper around fflush to catch and log flushes */ + fprintf(stderr,"fflush() at >>"); + print_trace(); + fprintf(stderr,"<< to stream %p which is ",stream); + + if(stream==NULL) + { + fprintf(stderr,"NULL"); + } + else if(stream==stdout) + { + fprintf(stderr,"stdout"); + } + else if(stream==stderr) + { + fprintf(stderr,"stderr"); + } +#ifdef FILE_LOG + else if(stream==global_stardata_p->model.log_fp) + { + fprintf(stderr,"file_log"); + } +#endif + else + { + fprintf(stderr,"unknown/private\n"); + } + fprintf(stderr,"\n"); + + /* undef fflush to use library call */ +#undef fflush + fflush(stderr); + return fflush(stream); +} + +/* Obtain a backtrace and print it to stdout. */ +static void print_trace (void) +{ + void *array[10]; + size_t size; + char **strings; + size_t i; + + size = backtrace (array, 10); + strings = backtrace_symbols (array, size); + //fprintf (stderr,"Obtained %zd stack frames.\n", size); + for (i = 0; i < size; i++) + { + fprintf (stderr,"%s > ", strings[i]); + } + free (strings); +} + diff --git a/src/main_prototypes.h b/src/main_prototypes.h index 3491f485f..5728204cb 100644 --- a/src/main_prototypes.h +++ b/src/main_prototypes.h @@ -31,5 +31,6 @@ void final_timing_output(ticks start_tick ); int buffer_stack_pop(char *c); +int binary_c_flush(FILE *stream); #endif diff --git a/src/maths/gaussj.c b/src/maths/gaussj.c index fd73c16d8..971480d78 100644 --- a/src/maths/gaussj.c +++ b/src/maths/gaussj.c @@ -18,13 +18,13 @@ void gaussj(double ** RESTRICT a,const int n,double ** RESTRICT b,const int m) unsigned int icol=0,irow=0; const unsigned int nn=n+1,mm=m+1; - Dprint("GAUSSJ n=%d m=%d\n",n,m);fflush(stdout); + Dprint("GAUSSJ n=%d m=%d\n",n,m); indxc=ivector(1,n); indxr=ivector(1,n); ipiv=ivector(1,n); for (j=1;j<nn;j++) ipiv[j]=0; - Dprint("Main loop\n");fflush(stdout); + Dprint("Main loop\n"); for (i=1;i<nn;i++) { @@ -48,23 +48,23 @@ void gaussj(double ** RESTRICT a,const int n,double ** RESTRICT b,const int m) } } ++(ipiv[icol]); - Dprint("Pivot at icol %d (irow=%d)\n",icol,irow);fflush(stdout); + Dprint("Pivot at icol %d (irow=%d)\n",icol,irow); if (irow != icol) { - Dprint("Swap a\n");fflush(stdout); + Dprint("Swap a\n"); for (l=1;l<nn;l++) SWAP(a[irow][l],a[icol][l]); - Dprint("Swap b[irow=%d][1] to b[icol=%d][1]\n",irow,icol);fflush(stdout); + Dprint("Swap b[irow=%d][1] to b[icol=%d][1]\n",irow,icol); for (l=1;l<mm;l++) SWAP(b[irow][l],b[icol][l]); - Dprint("ok\n");fflush(stdout); + Dprint("ok\n"); } - Dprint("Set indexes...\n");fflush(stdout); + Dprint("Set indexes...\n"); indxr[i]=irow; indxc[i]=icol; - Dprint("Done\n");fflush(stdout); + Dprint("Done\n"); if (IS_ZERO(a[icol][icol])) { @@ -72,7 +72,7 @@ void gaussj(double ** RESTRICT a,const int n,double ** RESTRICT b,const int m) exit(0); } pivinv=1.0/a[icol][icol]; - Dprint("Set pivinv\n");fflush(stdout); + Dprint("Set pivinv\n"); a[icol][icol]=1.0; for (l=1;l<nn;l++) a[icol][l] *= pivinv; @@ -87,7 +87,7 @@ void gaussj(double ** RESTRICT a,const int n,double ** RESTRICT b,const int m) } } } - Dprint("Loop 2\n");fflush(stdout); + Dprint("Loop 2\n"); for (l=n;l>0;l--) { if (indxr[l] != indxc[l]) @@ -104,7 +104,7 @@ void gaussj(double ** RESTRICT a,const int n,double ** RESTRICT b,const int m) free_ivector(ipiv,1); free_ivector(indxr,1); free_ivector(indxc,1); - Dprint("Exit gaussj\n");fflush(stdout); + Dprint("Exit gaussj\n"); } diff --git a/src/maths/interpolate.c b/src/maths/interpolate.c index 5c6868e3f..edd5dc433 100644 --- a/src/maths/interpolate.c +++ b/src/maths/interpolate.c @@ -191,7 +191,11 @@ static size_t interpolate_mem_allocated=0; #undef TINY #define TINY 1e-20 +#ifdef INTERPOLATE_DEBUG #define FLUSH fflush(stdout); +#else +#define FLUSH /* */ +#endif #ifdef NEWCACHE /* @@ -343,7 +347,7 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table printf("Interpolate: check cache (table_id=%d, cache_array=%p,results_array=%p)\n", table_id,cache_array,results_array); #endif - fflush(stdout); + FLUSH; #endif #ifdef INTERPOLATE_CACHE @@ -512,7 +516,7 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table * Setup / (re)allocate memory */ #ifdef INTERPOLATE_DEBUG - printf("Interpolate: memory allocation\n");fflush(stdout); + printf("Interpolate: memory allocation\n");FLUSH; #endif #ifdef INTERPOLATE_USE_REALLOC /* realloc a single piece of memory */ @@ -848,17 +852,17 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table k+=lnl; #ifdef INTERPOLATE_DEBUG - printf("Line %d : ",i);fflush(stdout); + printf("Line %d : ",i);FLUSH; for(j=0;j<n;j++) { - printf("% 3.3e ",*(int_table+i*lnl+j));fflush(stdout); + printf("% 3.3e ",*(int_table+i*lnl+j));FLUSH; } - printf(" | ");fflush(stdout); + printf(" | ");FLUSH; for(j=n;j<lnl;j++) { - printf("% 3.3e ",*(int_table+i*lnl+j));fflush(stdout); + printf("% 3.3e ",*(int_table+i*lnl+j));FLUSH; } - printf(" %d/%d\n",i,max_itc_size-1);fflush(stdout); + printf(" %d/%d\n",i,max_itc_size-1);FLUSH; #endif } @@ -885,7 +889,7 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table */ #ifdef INTERPOLATE_DEBUG printf("Interpolate dim=%d f=%g\n",dim,f[dim]); - fflush(stdout); + FLUSH; #endif #ifdef POINTER_ARITHMETIC v=*(f+dim); @@ -949,7 +953,7 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table } #ifdef INTERPOLATE_DEBUG - printf("memcopy results\n");fflush(stdout); + printf("memcopy results\n");FLUSH; #endif /* @@ -970,7 +974,7 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table { printf("% 3.3e ",*(int_table+j)); } - printf("\n");fflush(stdout); + printf("\n");FLUSH; #endif #ifdef INTERPOLATE_CACHE @@ -1009,17 +1013,17 @@ void interpolate(const double * RESTRICT table, // (pointer to) the data table /* free alloc'd memory */ #ifdef INTERPOLATE_DEBUG fprintf(stdout,"Free int_table=%p\n",int_table); - fflush(stdout); + FLUSH; #endif SAFE_FREE(int_table); #ifdef INTERPOLATE_DEBUG fprintf(stdout,"Free f = %p\n",f); - fflush(stdout); + FLUSH; #endif #ifdef INTERPOLATE_DEBUG fprintf(stderr,"Free sum = %p\n",sum); - fflush(stdout); + FLUSH; #endif SAFE_FREE(sum); #endif diff --git a/src/nucsyn/nucsyn_elemental_abundance.c b/src/nucsyn/nucsyn_elemental_abundance.c index f1a35b6f9..37f3c21a3 100644 --- a/src/nucsyn/nucsyn_elemental_abundance.c +++ b/src/nucsyn/nucsyn_elemental_abundance.c @@ -11,13 +11,14 @@ Abundance nucsyn_elemental_abundance(const char * element, Abundance Xelem=0.0; // calculated mass fraction (returned) - printf("Find abundance of element %s in %p (%g,%g,%g,%g,%g) XH1=%g XHe4=%g XFe56=%g\n",element, + /*printf("Find abundance of element %s in %p (%g,%g,%g,%g,%g) XH1=%g XHe4=%g XFe56=%g\n",element, X,X[0],X[1],X[2],X[3],X[4],X[XH1],X[XHe4],X[XFe56]); + */ /* calcualte the atomic number of the element */ Atomic_number Z=element_to_atomic_number(element); - printf("Z=%d\n",Z); + //printf("Z=%d\n",Z); /* * ID isotopes of this element @@ -27,13 +28,12 @@ Abundance nucsyn_elemental_abundance(const char * element, { if(stardata->common.atomic_number[i] == Z) { - printf("Add isotope number %d (abundance %g)\n",i,X[i]); + //printf("Add isotope number %d (abundance %g)\n",i,X[i]); Xelem += X[i]; } } - printf("Final abundance of %s is %g\n",element,Xelem); - + //printf("Final abundance of %s is %g\n",element,Xelem); return Xelem; } @@ -57,3 +57,4 @@ Atomic_number element_to_atomic_number(const char * element) exit(0); } #endif + diff --git a/src/nucsyn/nucsyn_gce_log.c b/src/nucsyn/nucsyn_gce_log.c index e1b4ff916..6016484de 100644 --- a/src/nucsyn/nucsyn_gce_log.c +++ b/src/nucsyn/nucsyn_gce_log.c @@ -4,77 +4,78 @@ #include "nucsyn_parameters.h" #include <stdio.h> +#define FLUSH /* fflush(stdout); */ #ifdef NUCSYN_GCE void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) { - /* log stuff for 2d plots in GCE, as well as output for the - * "present day" i.e. observations for NOW - */ + /* log stuff for 2d plots in GCE, as well as output for the + * "present day" i.e. observations for NOW + */ #ifdef NUCSYN_STELLAR_POPULATIONS_ENSEMBLE - Dprint("Enter GCE log\n"); + Dprint("Enter GCE log\n"); - double f; - double ensemble_weight[NUCSYN_ENSEMBLE_N]; - static double ensemble_store[NUMBER_OF_STARS+1][NUCSYN_ENSEMBLE_N]; - double teff,l,lacc,ncno,mu; - static double lastensembleouttime; - static double ftot=0.0; - unsigned int k; - Star_number i; - static Boolean rate_mask[NUCSYN_ENSEMBLE_N]; - Boolean reset_time=FALSE; + double f; + double ensemble_weight[NUCSYN_ENSEMBLE_N]; + static double ensemble_store[NUMBER_OF_STARS+1][NUCSYN_ENSEMBLE_N]; + double teff,l,lacc,ncno,mu; + static double lastensembleouttime; + static double ftot=0.0; + unsigned int k; + Star_number i; + static Boolean rate_mask[NUCSYN_ENSEMBLE_N]; + Boolean reset_time=FALSE; #endif // NUCSYN_STELLAR_POPULATIONS_ENSEMBLE #ifdef NUCSYN_PRESENT_DAY_ABUNDS - static double prevtpresentday; + static double prevtpresentday; #endif - double dt,dt_since_last; - static double prevt; - /* - * Calculate timestep, force to be zero for first call. - */ - if(stardata->model.model_time<0.01) + double dt,dt_since_last; + static double prevt; + /* + * Calculate timestep, force to be zero for first call. + */ + if(stardata->model.model_time<0.01) { - //printf("reset prevt\n"); - prevt=0.0; - dt=0.0; + //printf("reset prevt\n"); + prevt=0.0; + dt=0.0; #ifdef NUCSYN_PRESENT_DAY_ABUNDS - prevtpresentday=0.0; + prevtpresentday=0.0; #endif #ifdef NUCSYN_STELLAR_POPULATIONS_ENSEMBLE - for(k=0;k<NUCSYN_ENSEMBLE_N;k++) + for(k=0;k<NUCSYN_ENSEMBLE_N;k++) { - for(i=1;i<NUMBER_OF_STARS;i++) + for(i=1;i<NUMBER_OF_STARS;i++) { - ensemble_store[i][k]=0.0; + ensemble_store[i][k]=0.0; } - rate_mask[k]=FALSE; + rate_mask[k]=FALSE; } - /* some ensemble variables are rates : treat these differently */ - rate_mask[NUCSYN_ENSEMBLE_SN_IA_He]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_IA_ELD]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_IA_CHAND]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_AIC]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_IA_He_Coel]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_IA_CHAND_Coel]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_IBC]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_II]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_SN_IIa]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_GRB_NS_NS]=TRUE; - rate_mask[NUCSYN_ENSEMBLE_GRB_COLLAPSAR]=TRUE; + /* some ensemble variables are rates : treat these differently */ + rate_mask[NUCSYN_ENSEMBLE_SN_IA_He]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_IA_ELD]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_IA_CHAND]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_AIC]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_IA_He_Coel]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_IA_CHAND_Coel]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_IBC]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_II]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_SN_IIa]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_GRB_NS_NS]=TRUE; + rate_mask[NUCSYN_ENSEMBLE_GRB_COLLAPSAR]=TRUE; #endif } - else + else { - dt=stardata->model.model_time-prevt; - //printf("Calc dt=%g from t=%g prevt=%g\n",dt,stardata->model.model_time,prevt); + dt=stardata->model.model_time-prevt; + //printf("Calc dt=%g from t=%g prevt=%g\n",dt,stardata->model.model_time,prevt); } @@ -87,7 +88,7 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) for(i=1;i<=NUMBER_OF_STARS;i++) - { + { Dprint("Star %d",i); @@ -95,15 +96,15 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) * First set all detections to false. */ for(k=0;k<NUCSYN_ENSEMBLE_N;k++) - { + { ensemble_weight[k]=0.0; - } + } /* * Stars with zams mass 0.01 Msun are obviously pointless objects */ if(stardata->star[i].zams_mass > 0.02) - { + { Dprint("HRD"); /* * calculate effective temperature, log L @@ -136,11 +137,11 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) ((stardata->star[1].stellar_type==MASSLESS_REMNANT)&& (stardata->star[2].stellar_type!=MASSLESS_REMNANT)&& (i==2))) - { + { if((stardata->star[3-i].stellar_type==0)&& (stardata->star[3-i].mass<0.011)) - { + { /* * Probably a single star: or at least it would be * observed as such. Note: companion must be stellar type 0 @@ -149,14 +150,14 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) ensemble_weight[NUCSYN_ENSEMBLE_NUMBER_OF_SINGLE_STARS]=1; ensemble_weight[NUCSYN_ENSEMBLE_LUMINOSITY_OF_SINGLE_STARS]=stardata->star[i].luminosity; ensemble_weight[NUCSYN_ENSEMBLE_MASS_IN_SINGLE_STARS]=stardata->star[i].mass; - } + } else - { + { /* * Probably a binary star */ if(stardata->star[3-i].stellar_type==MASSLESS_REMNANT) - { + { /* * Merger: will be a single star AND a merged binary */ @@ -167,9 +168,9 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) ensemble_weight[NUCSYN_ENSEMBLE_LUMINOSITY_OF_SINGLE_STARS]=stardata->star[i].luminosity; ensemble_weight[NUCSYN_ENSEMBLE_MASS_IN_SINGLE_STARS]=stardata->star[i].mass; - } + } else if(stardata->common.separation > TINY) - { + { /* * Still a binary! */ @@ -177,9 +178,9 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) ensemble_weight[NUCSYN_ENSEMBLE_LUMINOSITY_OF_BINARY_STARS]=stardata->star[1].luminosity+stardata->star[2].luminosity; ensemble_weight[NUCSYN_ENSEMBLE_MASS_IN_BINARY_STARS]=stardata->star[1].mass+stardata->star[2].mass; - } + } else if(stardata->common.separation < TINY) - { + { /* * Disrupted binary: will appear as two single stars so count twice! */ @@ -190,11 +191,11 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) ensemble_weight[NUCSYN_ENSEMBLE_LUMINOSITY_OF_SINGLE_STARS]=stardata->star[1].luminosity+stardata->star[2].luminosity; ensemble_weight[NUCSYN_ENSEMBLE_MASS_IN_SINGLE_STARS]=stardata->star[1].mass+stardata->star[2].mass; - } + } - } + } - } + } /* * Save the number of each stellar type and the luminosity in each stellar type @@ -208,59 +209,59 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) * Check for WR stars */ if(stardata->star[i].stellar_type<10) - { + { // avoid WDs, which would qualify as WO stars! switch(nucsyn_WR_type(&(stardata->star[i]))) - { - case WR_WNL: + { + case WR_WNL: // WNL star ensemble_weight[NUCSYN_ENSEMBLE_WNL]=1; ensemble_weight[NUCSYN_ENSEMBLE_WNL_LUM]=stardata->star[i].luminosity; break; - case WR_WNE: + case WR_WNE: // WNE star ensemble_weight[NUCSYN_ENSEMBLE_WNE]=1; ensemble_weight[NUCSYN_ENSEMBLE_WNE_LUM]=stardata->star[i].luminosity; break; - case WR_WC: + case WR_WC: // WC star ensemble_weight[NUCSYN_ENSEMBLE_WC]=1; ensemble_weight[NUCSYN_ENSEMBLE_WC_LUM]=stardata->star[i].luminosity; break; - case WR_WO: + case WR_WO: // WO star ensemble_weight[NUCSYN_ENSEMBLE_WO]=1; ensemble_weight[NUCSYN_ENSEMBLE_WO_LUM]=stardata->star[i].luminosity; break; - } - } + } + } Dprint("B/RSG checks"); /* * Check for BSG or RSG */ if(stardata->star[i].Xsurf[XH1]>0.4) - { + { if(l>4.9) - { + { if(teff<3.66) - { + { ensemble_weight[NUCSYN_ENSEMBLE_RSG]=1; ensemble_weight[NUCSYN_ENSEMBLE_RSG_LUM]=stardata->star[i].luminosity; - } + } else if(teff>3.9) - { + { ensemble_weight[NUCSYN_ENSEMBLE_BSG]=1; ensemble_weight[NUCSYN_ENSEMBLE_BSG_LUM]=stardata->star[i].luminosity; - } + } else - { + { ensemble_weight[NUCSYN_ENSEMBLE_ISG]=1; ensemble_weight[NUCSYN_ENSEMBLE_ISG_LUM]=stardata->star[i].luminosity; - } - } - } + } + } + } Dprint("SN checks"); /* @@ -268,84 +269,84 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) * when calculating populations */ if(stardata->star[i].sn_last_time!=-1) - { + { switch(stardata->star[i].sn_last_time) - { - case SN_IA_He: + { + case SN_IA_He: ensemble_weight[NUCSYN_ENSEMBLE_SN_IA_He]=1; break; - case SN_IA_ELD: + case SN_IA_ELD: ensemble_weight[NUCSYN_ENSEMBLE_SN_IA_ELD]=1; break; - case SN_IA_CHAND: + case SN_IA_CHAND: ensemble_weight[NUCSYN_ENSEMBLE_SN_IA_CHAND]=1; break; - case SN_AIC: + case SN_AIC: ensemble_weight[NUCSYN_ENSEMBLE_SN_AIC]=1; break; - case SN_IA_He_Coel: + case SN_IA_He_Coel: ensemble_weight[NUCSYN_ENSEMBLE_SN_IA_He_Coel]=1; break; - case SN_IA_CHAND_Coel: + case SN_IA_CHAND_Coel: ensemble_weight[NUCSYN_ENSEMBLE_SN_IA_CHAND_Coel]=1; break; - case SN_IBC: + case SN_IBC: ensemble_weight[NUCSYN_ENSEMBLE_SN_IBC]=1; break; - case SN_II: + case SN_II: ensemble_weight[NUCSYN_ENSEMBLE_SN_II]=1; break; - case SN_IIa: + case SN_IIa: ensemble_weight[NUCSYN_ENSEMBLE_SN_IIa]=1; break; - case SN_GRB_NS_NS: + case SN_GRB_NS_NS: ensemble_weight[NUCSYN_ENSEMBLE_GRB_NS_NS]=1; break; - case SN_GRB_COLLAPSAR: + case SN_GRB_COLLAPSAR: /* * A collapsar is also a SN 1bc, so set SN_IBC as well */ ensemble_weight[NUCSYN_ENSEMBLE_GRB_COLLAPSAR]=1; ensemble_weight[NUCSYN_ENSEMBLE_SN_IBC]=1; break; - } + } // reset sn_last_time stardata->star[i].sn_last_time=-1; - } + } Dprint("Xray binary checks"); /* * X-ray binaries: luminosity is calculated in log_xray_binary */ if(stardata->star[i].Xray_luminosity > TINY) - { + { if(ON_MAIN_SEQUENCE(stardata->star[3-i].stellar_type)&& (stardata->star[3-i].mass<=2.0)) - { + { /* * LMXB */ ensemble_weight[NUCSYN_ENSEMBLE_XRAY_LMXB]=1; ensemble_weight[NUCSYN_ENSEMBLE_XRAY_LMXB_LUM]=stardata->star[i].Xray_luminosity; - } + } else if(WHITE_DWARF(stardata->star[3-i].stellar_type)) - { + { /* * WDXB */ ensemble_weight[NUCSYN_ENSEMBLE_XRAY_WDXB]=1; ensemble_weight[NUCSYN_ENSEMBLE_XRAY_WDXB_LUM]=stardata->star[i].Xray_luminosity; - } + } else - { + { /* * HMXB */ ensemble_weight[NUCSYN_ENSEMBLE_XRAY_HMXB]=1; ensemble_weight[NUCSYN_ENSEMBLE_XRAY_HMXB_LUM]=stardata->star[i].Xray_luminosity; //printf("XRAY HMXB Detect Lx=%g\n",stardata->star[i].Xray_luminosity); - } - } + } + } Dprint("Spectral types"); /* @@ -355,157 +356,159 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) */ if((stardata->star[i].stellar_type>2)&& (stardata->star[i].stellar_type<10)) - { + { if(l>4.9) - { + { // supergiant? if(teff>log10(26000)) - { + { // O ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_O]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_O_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(9730)) - { + { // B ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_B]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_B_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(7700)) - { + { // A ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_A]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_A_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(6550)) - { + { // F ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_F]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_F_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(4420)) - { + { // G ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_G]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_G_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(3650)) - { + { // K ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K_LUM]=stardata->star[i].luminosity; - } + } else - { + { // M ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M_LUM]=stardata->star[i].luminosity; - } - } + } + } else - { + { // giant if(teff>log10(29000)) - { + { // O ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_O]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_O_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(8100)) - { + { // B ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_B]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_B_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(6470)) - { + { // A ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_A]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_A_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(5150)) - { + { // F ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_F]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_F_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(3950)) - { + { // G ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_G]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_G_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(3620)) - { + { // K ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K_LUM]=stardata->star[i].luminosity; - } + } else - { + { // M ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M_LUM]=stardata->star[i].luminosity; - } - } - } + } + } + } else - { + { // everything else is a dwarf if(teff>log10(30000)) - { + { // O ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_O]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_O_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(9520)) - { + { // B ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_B]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_B_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(6030)) - { + { // A ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_A]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_A_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(5250)) - { + { // F ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_F]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_F_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(3850)) - { + { // G ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_G]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_G_LUM]=stardata->star[i].luminosity; - } + } else if(teff>log10(3370)) - { + { // K ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K_LUM]=stardata->star[i].luminosity; - } + } else - { + { // M ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M_LUM]=stardata->star[i].luminosity; - } - } + } + } + + //ncno=(stardata->star[i].Xsurf[XC12]/12.0)/ + //(TINY+stardata->star[i].Xsurf[XO16]/16.0); + ncno=nucsyn_elemental_abundance("C",stardata->star[i].Xsurf,stardata)/ + nucsyn_elemental_abundance("O",stardata->star[i].Xsurf,stardata); - ncno=(stardata->star[i].Xsurf[XC12]/12.0)/ - (TINY+stardata->star[i].Xsurf[XO16]/16.0); - if(ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_M]==1) - { + { /* * M-type stars can be C stars */ if(ncno > 1.05) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_C]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_C_LUM]=stardata->star[i].luminosity; @@ -514,29 +517,29 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) */ if(((stardata->star[i].Xsurf[XC12]/12.0)/ (stardata->star[i].Xsurf[XC13]/13.0))<12.0) - { - ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_J]=1; - ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_J_LUM]=stardata->star[i].luminosity; - } + { + ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_J]=1; + ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_J_LUM]=stardata->star[i].luminosity; + } - } + } else if(ncno > 0.95) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_S]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_S_LUM]=stardata->star[i].luminosity; - } - } + } + } else if(ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_K]==1) - { + { /* * K-type stars which are C-rich are assumed to be R-stars */ if(ncno>1.0) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_R]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_R_LUM]=stardata->star[i].luminosity; - } - } + } + } /* * Red G/K type giants which are barium rich @@ -548,20 +551,21 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) (nucsyn_elemental_abundance("Ba",stardata->star[i].Xsurf,stardata)/ nucsyn_elemental_abundance("Ba",stardata->common.XZAMS,stardata)>2.0) ) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_Ba]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_Ba_LUM]=stardata->star[i].luminosity; - } + } /* * Look for Tc-rich stars which are GB stars? */ if((stardata->star[i].stellar_type==GIANT_BRANCH)&& - (stardata->star[i].Xsurf[XTc]/stardata->common.XZAMS[XTc]>2.0)) - { + ((nucsyn_elemental_abundance("Tc",stardata->star[i].Xsurf,stardata)/ + nucsyn_elemental_abundance("Tc",stardata->common.XZAMS,stardata)>2.0))) + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_Tc]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_Tc_LUM]=stardata->star[i].luminosity; - } + } /* * Look for extrinsic dwarf carbon stars: @@ -571,98 +575,98 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) if((stardata->star[i].stellar_type<=6)&& (ncno>1.0)&& (WHITE_DWARF(stardata->star[3-i].stellar_type))) - { + { if(ON_MAIN_SEQUENCE(stardata->star[i].stellar_type)) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_EXTRINSIC_DWARF_C]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_EXTRINSIC_DWARF_C_LUM]=stardata->star[i].luminosity; - } + } else - { + { /* * C-rich pre-TPAGB stars are "extrinsic" */ if(stardata->star[i].num_thermal_pulses_since_mcmin<0.5) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_EXTRINSIC_GIANT_C]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_EXTRINSIC_GIANT_C_LUM]=stardata->star[i].luminosity; - } - } - } + } + } + } /* * Classify WR stars by Hurley Pols Tout 2002 prescription: * Assume WDs cannot be WR stars (dubious?) */ if((stardata->star[i].stellar_type>1)&& (stardata->star[i].stellar_type<10)) - { + { mu = ((stardata->star[i].mass-stardata->star[i].core_mass)/stardata->star[i].mass)*MIN(5.0,MAX(1.20,pow((stardata->star[i].luminosity/LUM0),KAP))); if(mu<1.0) - { + { ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_WR]=1; ensemble_weight[NUCSYN_ENSEMBLE_SPECTRAL_TYPE_WR_LUM]=stardata->star[i].luminosity; - } - } + } + } /* * Luminous blue variable */ if((stardata->star[i].luminosity>6e5)&& (1e-5*stardata->star[i].radius*sqrt(stardata->star[i].luminosity)>1.0)) - { + { ensemble_weight[NUCSYN_ENSEMBLE_LBV]=1; ensemble_weight[NUCSYN_ENSEMBLE_LBV_LUM]=stardata->star[i].luminosity; - } + } Dprint("Special binaries"); /* * Special binary stellar systems */ if(i==1) - { + { /* * Check for double degenerates (only check star 1) */ if(WHITE_DWARF(stardata->star[i].stellar_type)) - { + { /* * Primary is a white dwarf */ if(WHITE_DWARF(stardata->star[3-i].stellar_type)) - { + { /* * WD WD system */ ensemble_weight[NUCSYN_ENSEMBLE_WDWD_DD]=1; - } + } else if(stardata->star[3-i].stellar_type==NS) - { + { ensemble_weight[NUCSYN_ENSEMBLE_WDNS_DD]=1; - } + } else if(stardata->star[3-i].stellar_type==BH) - { + { ensemble_weight[NUCSYN_ENSEMBLE_WDBH_DD]=1; - } - } + } + } else if(stardata->star[i].stellar_type==NS) - { + { if(stardata->star[3-i].stellar_type==NS) - { + { ensemble_weight[NUCSYN_ENSEMBLE_NSNS_DD]=1; - } + } else if(stardata->star[3-i].stellar_type==BH) - { + { ensemble_weight[NUCSYN_ENSEMBLE_NSBH_DD]=1; - } - } + } + } else if(stardata->star[i].stellar_type==BH) - { + { if(stardata->star[3-i].stellar_type==BH) - { + { ensemble_weight[NUCSYN_ENSEMBLE_BHBH_DD]=1; - } - } - } + } + } + } /* * Cataclysmic variables @@ -671,40 +675,40 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) if(WHITE_DWARF(stardata->star[i].stellar_type)&& (stardata->star[3-i].radius > stardata->star[3-i].rol)&& (stardata->star[3-i].stellar_type<10)) - { + { /* * CV */ if(stardata->star[3-i].stellar_type<=1) - { + { /* * Donor on MS: classical CV */ ensemble_weight[NUCSYN_ENSEMBLE_CLASSICAL_CV]=1; - } + } else if(stardata->star[3-i].stellar_type==2) - { + { /* * Donor is HG star: GK Per CV */ ensemble_weight[NUCSYN_ENSEMBLE_GK_PER_CV]=1; - } + } else if((stardata->star[3-i].stellar_type>=3)&& (stardata->star[3-i].stellar_type<=6)) - { + { /* * Donor is helium-burning: Symbiotic star */ ensemble_weight[NUCSYN_ENSEMBLE_SYMBIOTIC_CV]=1; - } + } else if(NAKED_HELIUM_STAR(stardata->star[3-i].stellar_type)) - { + { /* * Dubious check for SDB star! */ ensemble_weight[NUCSYN_ENSEMBLE_SDB]=1; - } - } + } + } /* @@ -713,55 +717,55 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) Dprint("Algols"); if(ON_MAIN_SEQUENCE(stardata->star[i].stellar_type)&& (stardata->star[3-i].radius > stardata->star[3-i].rol)) - { + { if(stardata->star[3-i].mass/stardata->star[i].mass > 1.0) - { + { ensemble_weight[NUCSYN_ENSEMBLE_PRE_ALGOL]=1; - } + } else - { + { if(ON_MAIN_SEQUENCE(stardata->star[3-i].stellar_type)) - { + { ensemble_weight[NUCSYN_ENSEMBLE_MS_ALGOL]=1; - } + } else - { + { if(stardata->star[i].mass<=1.25) - { + { ensemble_weight[NUCSYN_ENSEMBLE_COLD_ALGOL]=1; - } + } else - { + { ensemble_weight[NUCSYN_ENSEMBLE_HOT_ALGOL]=1; - } - } - } - } + } + } + } + } Dprint("Blue stragglers"); /* * Blue stragglers (set in acquire_stellar_parameters...) */ if((stardata->model.bss==TRUE)&&(i==1)) - { + { ensemble_weight[NUCSYN_ENSEMBLE_BLUE_STRAGGLER]=1; - } + } Dprint("Low-mass WDs"); /* * Low-mass white dwarfs */ if(stardata->star[i].mass<0.5) - { + { if(stardata->star[i].stellar_type==HeWD) - { + { ensemble_weight[NUCSYN_ENSEMBLE_LOW_MASS_He_WD]=1; - } + } else if(stardata->star[i].stellar_type==COWD) - { + { ensemble_weight[NUCSYN_ENSEMBLE_LOW_MASS_CO_WD]=1; - } - } + } + } Dprint("Symbiotic stars"); /* @@ -769,28 +773,28 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) */ if(stardata->star[i].stellar_type<=12) - { + { lacc=3.14E+07*stardata->star[i].mass*stardata->star[i].dmt/stardata->star[i].radius; if((lacc>10.0)||(lacc/stardata->star[3-i].luminosity>0.01)) - { + { if(stardata->star[3-i].stellar_type==6) - { + { /* * D-type symbiotic */ ensemble_weight[NUCSYN_ENSEMBLE_D_TYPE_SYMBIOTIC]=1; - } + } else if(stardata->star[3-i].stellar_type<6) - { + { /* * S-type symbiotic */ ensemble_weight[NUCSYN_ENSEMBLE_S_TYPE_SYMBIOTIC]=1; - } - } - } + } + } + } /* * Output @@ -799,27 +803,27 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) #ifdef NUCSYN_STELLAR_POPULATIONS_ENSEMBLE_SPARSE if(stardata->model.model_time<0.01) - { + { /* * First timestep */ lastensembleouttime=0.0; - } + } /* add up the ensemble store */ for(k=0;k<NUCSYN_ENSEMBLE_N;k++) - { + { if(rate_mask[k]==TRUE) - { + { /* this is a rate: just count */ ensemble_store[i][k]+=ensemble_weight[k]; - } + } else - { + { /* weight by timestep */ ensemble_store[i][k]+=ensemble_weight[k]*dt; - } - } + } + } //printf("ADD star %d weight=%g dt=%g t=%g\n", //i,ensemble_weight[23],dt,stardata->model.model_time); @@ -832,77 +836,77 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) /* output (if > the next output time) */ if(MORE_OR_EQUAL(dt_since_last,stardata->preferences->yields_dt)) - { - printf("STELLARPOPENSEMBLE%d__ %6.6e ", - i, - stardata->model.model_time); - //dt_since_last, - //stardata->model.probability); + { + printf("STELLARPOPENSEMBLE%d__ %6.6e ", + i, + stardata->model.model_time); + //dt_since_last, + //stardata->model.probability); - /* - * We want to output the value of x taken as a - * time-averaged mean over the time bin. - * - * The values which are in ensemble_store are - * - * sum(x * delta t) - * - * where delta t is the each *code* timestep - * - * x is the variable in question (=1 for simple counts, but - * possibly a true stellar variable) - * - * We must now modulate by the probability to get the true - * 'dtp' contribution to each time bin, then divide by the - * time since the last log to get the time-average. - * - * Rates are treated identically except that they are counted as - * - * sum = # of events - * - * i.e. there is no 'delta t' term, hence the rate is output directly. - */ - - - double fac=stardata->model.probability/dt_since_last; - - for(k=0;k<NUCSYN_ENSEMBLE_N;k++) - { - if(rate_mask[k]==TRUE) - { - printf("%g ",fac*ensemble_store[i][k]); - } - else - { - //printf("%d=%g ",k,fac*ensemble_store[i][k]); - printf("%g ",fac*ensemble_store[i][k]); - } - - /* reset */ - ensemble_store[i][k]=0.0; - } - printf("\n"); - reset_time=TRUE; - } + /* + * We want to output the value of x taken as a + * time-averaged mean over the time bin. + * + * The values which are in ensemble_store are + * + * sum(x * delta t) + * + * where delta t is the each *code* timestep + * + * x is the variable in question (=1 for simple counts, but + * possibly a true stellar variable) + * + * We must now modulate by the probability to get the true + * 'dtp' contribution to each time bin, then divide by the + * time since the last log to get the time-average. + * + * Rates are treated identically except that they are counted as + * + * sum = # of events + * + * i.e. there is no 'delta t' term, hence the rate is output directly. + */ + + + double fac=stardata->model.probability/dt_since_last; + + for(k=0;k<NUCSYN_ENSEMBLE_N;k++) + { + if(rate_mask[k]==TRUE) + { + printf("%g ",fac*ensemble_store[i][k]); + } + else + { + //printf("%d=%g ",k,fac*ensemble_store[i][k]); + printf("%g ",fac*ensemble_store[i][k]); + } + + /* reset */ + ensemble_store[i][k]=0.0; + } + printf("\n"); + reset_time=TRUE; + } #else /* output now */ printf("STELLARPOPENSEMBLE%d__ %6.6e %3.3e ", i, stardata->model.model_time, stardata->model.probability); - fflush(stdout); + FLUSH; for(k=0;k<NUCSYN_ENSEMBLE_N;k++) - { + { //printf("%d=%g ",k,ensemble_weight[k]); printf("%g ",ensemble_weight[k]); // nb use %g because the weight might be an integer - } + } printf("\n"); #endif - } - } // i loop over stars + } + } // i loop over stars if(reset_time==TRUE) lastensembleouttime=stardata->model.model_time; @@ -912,82 +916,82 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata) #ifdef NUCSYN_PRESENT_DAY_ABUNDS Dprint("Present day abundances\n"); - if(stardata->preferences->present_day > 0) - { - - if((stardata->preferences->present_day_dumped==FALSE)&& - (MORE_OR_EQUAL(stardata->model.model_time, - stardata->preferences->present_day-stardata->preferences->present_day_period))) - { - if(stardata->model.model_time>stardata->preferences->present_day) - { - // set the flag to say we are past the present day - // so we should stop output - stardata->preferences->present_day_dumped=TRUE; - } + if(stardata->preferences->present_day > 0) + { + + if((stardata->preferences->present_day_dumped==FALSE)&& + (MORE_OR_EQUAL(stardata->model.model_time, + stardata->preferences->present_day-stardata->preferences->present_day_period))) + { + if(stardata->model.model_time>stardata->preferences->present_day) + { + // set the flag to say we are past the present day + // so we should stop output + stardata->preferences->present_day_dumped=TRUE; + } - /*printf("Check %g > %g > %g\n", - stardata->preferences->present_day, - stardata->model.model_time, - stardata->preferences->present_day-stardata->preferences->present_day_period); - */ + /*printf("Check %g > %g > %g\n", + stardata->preferences->present_day, + stardata->model.model_time, + stardata->preferences->present_day-stardata->preferences->present_day_period); + */ - /* - * Time-weighting factor: this is the fraction of the present_day_period - * during this timestep - */ - f=stardata->model.probability* - (MIN(stardata->model.model_time,stardata->preferences->present_day) - -MAX(stardata->preferences->present_day-stardata->preferences->present_day_period, - prevtpresentday))/stardata->preferences->present_day_period; - f=MIN(1.0,f); // just in case! + /* + * Time-weighting factor: this is the fraction of the present_day_period + * during this timestep + */ + f=stardata->model.probability* + (MIN(stardata->model.model_time,stardata->preferences->present_day) + -MAX(stardata->preferences->present_day-stardata->preferences->present_day_period, + prevtpresentday))/stardata->preferences->present_day_period; + f=MIN(1.0,f); // just in case! - prevtpresentday=stardata->model.model_time; + prevtpresentday=stardata->model.model_time; - ftot+=f; + ftot+=f; - /* - printf("t=%g : Weight f=%g (ftot %g) from p=%g dt=%g / per=%g\n", - stardata->model.model_time, - f,ftot, - stardata->model.probability, - dt, - stardata->preferences->present_day_period); - */ - // dump output for NOW for each star - if(f>TINY) - { - for(i=1;i<=NUMBER_OF_STARS;i++) - { - printf("GCE_PRESENT_DAY%d%g %d %g %g %g %g %g ", - i, - f, - stardata->star[i].stellar_type, - stardata->star[i].luminosity, - TEFF(i), - log10(GRAVITATIONAL_CONSTANT* - (M_SUN*stardata->star[i].mass)/ - POW2(R_SUN*stardata->star[i].radius)) - ,stardata->star[i].radius - ); - for(k=0;k<GCE_ISOTOPES;k++) - { - printf("%0+3e ",stardata->star[i].Xsurf[k]); - } - printf("\n"); - } - } - } - } + /* + printf("t=%g : Weight f=%g (ftot %g) from p=%g dt=%g / per=%g\n", + stardata->model.model_time, + f,ftot, + stardata->model.probability, + dt, + stardata->preferences->present_day_period); + */ + // dump output for NOW for each star + if(f>TINY) + { + for(i=1;i<=NUMBER_OF_STARS;i++) + { + printf("GCE_PRESENT_DAY%d%g %d %g %g %g %g %g ", + i, + f, + stardata->star[i].stellar_type, + stardata->star[i].luminosity, + TEFF(i), + log10(GRAVITATIONAL_CONSTANT* + (M_SUN*stardata->star[i].mass)/ + POW2(R_SUN*stardata->star[i].radius)) + ,stardata->star[i].radius + ); + for(k=0;k<GCE_ISOTOPES;k++) + { + printf("%0+3e ",stardata->star[i].Xsurf[k]); + } + printf("\n"); + } + } + } + } #endif // NUCSYN_PRESENT_DAY_ABUNDS - Dprint("Set prevt\n"); + Dprint("Set prevt\n"); - prevt=stardata->model.model_time; - //printf("Set prevt %g\n",prevt); + prevt=stardata->model.model_time; + //printf("Set prevt %g\n",prevt); } #endif // NUCSYN_GCE diff --git a/src/nucsyn/nucsyn_hbb.c b/src/nucsyn/nucsyn_hbb.c index a7271ab3e..41dd55221 100644 --- a/src/nucsyn/nucsyn_hbb.c +++ b/src/nucsyn/nucsyn_hbb.c @@ -8,7 +8,7 @@ #define TIMESTEP (1e6*stardata->model.dtm) //#define NUCSYN_TPAGB_HBB_CLAYTON_TABLE #ifdef NUCSYN_TPAGB_HBB_CLAYTON_TABLE - void nucsyn_clayton_table(double *mnuc); +void nucsyn_clayton_table(double *mnuc); #endif #if (DEBUG==1) @@ -24,7 +24,7 @@ void nucsyn_hbb(const double T, /* (log10) temperature */ Abundance * RESTRICT Xhbb, /* incoming mass fractions */ struct star_t *star, struct stardata_t * stardata - ) + ) { const Nuclear_mass *mnuc = stardata->common.mnuc; const Nuclear_mass *imnuc = stardata->common.imnuc; @@ -61,13 +61,13 @@ void nucsyn_hbb(const double T, /* (log10) temperature */ #ifdef NANCHECKS for(i=0;i<ISOTOPE_ARRAY_SIZE;i++) - { + { if(isnan(Nhbb[i])!=0) - { + { printf("nucsyn_hbb (preburn) isotope %d () is %g\n",i,Nhbb[i]); exit(0); - } - } + } + } #endif /* do the burn ( N(t) -> N(t+dt) ) */ @@ -75,13 +75,13 @@ void nucsyn_hbb(const double T, /* (log10) temperature */ #ifdef NANCHECKS for(i=0;i<ISOTOPE_ARRAY_SIZE;i++) - { + { if(isnan(Nhbb[i])!=0) - { + { printf("nucsyn_hbb (postburn) isotope %d is %g\n",i,Nhbb[i]); exit(0); - } - } + } + } #endif /* covert number density to mass fractions */ @@ -95,23 +95,23 @@ void nucsyn_hbb(const double T, /* (log10) temperature */ */ Xtot=XXsum(Xhbb); if(fabs(Xtot-1.0)>0.1) - { + { fflush(stdout); fprintf(stdout,"Xtot (%g) differs from 1.0 by > 1%% in nucsyn_hbb\n",Xtot); for(i=0;i<ISOTOPE_ARRAY_SIZE;i++) - { + { printf("X%d=%g ",i,Xhbb[i]); - } + } printf("\n"); fflush(stdout); fflush(stderr); exit(HBB_MASS_FRACTION_ERROR); - } + } Xtot=1.0/Xtot; for(i=0;i<ISOTOPE_ARRAY_SIZE;i++) - { + { Xhbb[i]*=Xtot; - } + } #endif return; diff --git a/src/nucsyn/nucsyn_initial_abundances.c b/src/nucsyn/nucsyn_initial_abundances.c index f52fa24c8..debb65290 100644 --- a/src/nucsyn/nucsyn_initial_abundances.c +++ b/src/nucsyn/nucsyn_initial_abundances.c @@ -603,7 +603,6 @@ void nucsyn_initial_abundances(Abundance * X, Dprint("m=%g\n",m); X[XH1]+=(1-m); - Dprint("Abundance mixture from Z=%g (Xtot now %g star's Z=%g)\n",z,XXsum(X),XXsum(X)-X[XH1]-X[XHe4]); Dprint("H %g He %g C12 %g C13 %g N14 %g O16 %g\nLOGS : H %g He %g C12 %g C13 %g N14 %g O16 %g\n", diff --git a/src/nucsyn/nucsyn_log.c b/src/nucsyn/nucsyn_log.c index 897267be7..3f5e1067d 100644 --- a/src/nucsyn/nucsyn_log.c +++ b/src/nucsyn/nucsyn_log.c @@ -32,6 +32,8 @@ /* only log STAR_TYPE */ #define STAR_TYPE TPAGB +#define FLUSH /* fflush(OUTPUT); */ + #ifdef NUCSYN_LOG /* Function to output details of nucleosynthesis to a log */ void nucsyn_log(const double t, @@ -81,7 +83,7 @@ void nucsyn_log(const double t, fprintf(OUTPUT,"%s%-14.14s%-14.14s%-14.14s", HEADER, "#Time","Mass","Mc"); - fflush(OUTPUT); + FLUSH; count=1; for(i=0;i<ISOTOPE_ARRAY_SIZE;i++) @@ -166,9 +168,8 @@ void nucsyn_log(const double t, ); /* number ratios this time! */ - fprintf(OUTPUT,"%8e %8e ",NRATIO(XC12,XC13),NCORATIO); - fprintf(OUTPUT,"\n"); /* trailing newline */ - fflush(OUTPUT); + fprintf(OUTPUT,"%8e %8e \n",NRATIO(XC12,XC13),NCORATIO); + FLUSH; } diff --git a/src/nucsyn/nucsyn_main_sequence_convective_envelope_size_table.c b/src/nucsyn/nucsyn_main_sequence_convective_envelope_size_table.c index 41511e5ad..0ae70722c 100644 --- a/src/nucsyn/nucsyn_main_sequence_convective_envelope_size_table.c +++ b/src/nucsyn/nucsyn_main_sequence_convective_envelope_size_table.c @@ -62,7 +62,7 @@ double lookup_MS_conv_size(double Z,double m_ZAMS,double frac_age,struct stardat while(m>0.0) { - Dprint(" Check at m=%g\n",m);fflush(stdout); + Dprint(" Check at m=%g\n",m); nucsyn_temperature_density_and_convection(m_ZAMS, m, @@ -73,7 +73,7 @@ double lookup_MS_conv_size(double Z,double m_ZAMS,double frac_age,struct stardat stardata, MAIN_SEQUENCE); - Dprint(" -> Convective = %d\n",convective); fflush(stdout); + Dprint(" -> Convective = %d\n",convective); if(convective==FALSE) return(MAX(0.0,m_ZAMS-(m+0.5*dm))); m-=dm; diff --git a/src/nucsyn/nucsyn_mc1tp.c b/src/nucsyn/nucsyn_mc1tp.c index f746cba77..d0e682af7 100644 --- a/src/nucsyn/nucsyn_mc1tp.c +++ b/src/nucsyn/nucsyn_mc1tp.c @@ -249,7 +249,7 @@ double nucsyn_mc1tp (double m, interpolate(data,1,1,12,coeffs,x,TABLE_POELARENDS_MC1TP,TRUE); } - Dprint("Interpolate for mc1tp m=%g mc1tp=%g\n",m,x[0]);fflush(stdout); + Dprint("Interpolate for mc1tp m=%g mc1tp=%g\n",m,x[0]); /* phase in */ y=AJ_SMOOTH(m,y,x[0]); diff --git a/src/nucsyn/nucsyn_parameters.h b/src/nucsyn/nucsyn_parameters.h index 7db4d5d22..baeaa337c 100644 --- a/src/nucsyn/nucsyn_parameters.h +++ b/src/nucsyn/nucsyn_parameters.h @@ -20,8 +20,6 @@ * nucleosynthesis routines. */ - - /* First dredge up */ /* Activate first dredge up : you want this, it uses fits to Amanda's @@ -191,7 +189,6 @@ #define USE_TABULAR_INTERSHELL_ABUNDANCES_KARAKAS_2002 //#define USE_TABULAR_INTERSHELL_ABUNDANCES_KARAKAS_2012 - // perhaps fudge the Ne22 in the 5msun z=1e-4 model (only if // you use tabular intershell abundances) @@ -262,7 +259,6 @@ /* HBB in TPAGB : you want this */ #define NUCSYN_TPAGB_HBB - /* * How would you like to burn? * Choose NUCSYN_ANAL_BURN for the old-fashioned, but reliable, @@ -715,7 +711,7 @@ * probably best to leave it and turn logging on and off below for * each specific case. */ -#define NUCSYN_LOGGING +//#define NUCSYN_LOGGING // Used to calibrate the free parameters for HBB to amanda's models. // You probably don't want it because I have done the calibration and @@ -843,7 +839,7 @@ //#define NUCSYN_SHORT_LOG /* normal log : use with cf_amanda.pl and cf_models-varMZ.pl */ -//#define NUCSYN_LOG +#define NUCSYN_LOG /* special log for testing globular cluster abundance ratios */ //#define NUCSYN_Globular_Cluster_LOG @@ -1399,11 +1395,8 @@ #undef NUCSYN_THIRD_DREDGE_UP_TABULAR_NCAL #undef NUCSYN_THIRD_DREDGE_UP_HYDROGEN_SHELL -#undef NUCSYN_THIRD_DREDGE_UP_MULTIPLIERS +#undef NUCSYN_THIRD_DREDGE_UP_MULTIPLIER -#undef MINIMUM_ENVELOPE_MASS_FOR_THIRD_DREDGEUP -#define MINIMUM_ENVELOPE_MASS_FOR_THIRD_DREDGEUP (0.5) - #undef USE_TABULAR_INTERSHELL_ABUNDANCES #undef NUCSYN_TPAGB diff --git a/src/nucsyn/nucsyn_preset_sigmav.c b/src/nucsyn/nucsyn_preset_sigmav.c index e4965584d..ae54ad42b 100644 --- a/src/nucsyn/nucsyn_preset_sigmav.c +++ b/src/nucsyn/nucsyn_preset_sigmav.c @@ -95,14 +95,14 @@ void print_sigmav_table(struct stardata_t * RESTRICT stardata) printf("Sigmav table\n"); for(i=0;i<NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION;i++) { - printf("Sigmav table line %d at %p: ",i,stardata->preferences->sigmav+i*sigmav_size);fflush(stdout); + printf("Sigmav table line %d at %p: ",i,stardata->preferences->sigmav+i*sigmav_size); for(j=0;j<6;j++) { printf("%d=%g ",j,*(stardata->preferences->sigmav+i*sigmav_size+j)); fflush(stdout); } printf("\n"); - fflush(stdout); + } } #endif diff --git a/src/nucsyn/nucsyn_radioactive_decay.c b/src/nucsyn/nucsyn_radioactive_decay.c index b61b67b5d..4952cf90f 100644 --- a/src/nucsyn/nucsyn_radioactive_decay.c +++ b/src/nucsyn/nucsyn_radioactive_decay.c @@ -237,7 +237,7 @@ void nucsyn_radioactive_decay(struct stardata_t * RESTRICT stardata) #endif } - Dprint("Done radioactive decay\n");fflush(stdout); + Dprint("Done radioactive decay\n"); } diff --git a/src/nucsyn/nucsyn_set_1st_dup_abunds.c b/src/nucsyn/nucsyn_set_1st_dup_abunds.c index 41645d7f2..1f05d2b53 100644 --- a/src/nucsyn/nucsyn_set_1st_dup_abunds.c +++ b/src/nucsyn/nucsyn_set_1st_dup_abunds.c @@ -494,7 +494,15 @@ void nucsyn_first_dup_apply(double m,double z, Dprint("# Richard's table (modified) r array : DH1=%g DHe4=%g DC12=%g DC13=%g DN14=%g DO16=%g\n", r[38],r[39],r[40],r[5],r[41],r[42]); - X[XH2]+=r[0];X[XHe3]+=r[1];X[XLi7]+=r[2];X[XBe7]+=r[3];X[XB11]+=r[4];X[XC13]+=r[5];X[XC14]+=r[6];X[XN15]+=r[7];X[XO17]+=r[8];X[XO18]+=r[9];X[XF19]+=r[10];X[XNe21]+=r[11];X[XNe22]+=r[12];X[XNa22]+=r[13];X[XNa23]+=r[14];X[XMg24]+=r[15];X[XMg25]+=r[16];X[XMg26]+=r[17];X[XAl26]+=r[18]+r[19]; X[XAl27]+=r[20];X[XSi28]+=r[21];X[XSi29]+=r[22];X[XSi30]+=r[23];X[XP31]+=r[24];X[XS32]+=r[25];X[XS33]+=r[26];X[XS34]+=r[27];X[XFe56]+=r[28];X[XFe57]+=r[29];X[XFe58]+=r[30];X[XFe59]+=r[31];X[XFe60]+=r[32];X[XCo59]+=r[33];X[XNi58]+=r[34];X[XNi59]+=r[35];X[XNi60]+=r[36];X[XNi61]+=r[37];X[XH1]+=r[38];X[XHe4]+=r[39];X[XC12]+=r[40];X[XN14]+=r[41];X[XO16]+=r[42];X[XNe20]+=r[43]; + X[XH2]+=r[0];X[XHe3]+=r[1];X[XLi7]+=r[2];X[XBe7]+=r[3];X[XB11]+=r[4];X[XC13]+=r[5]; +//X[XC14]+=r[6]; +X[XN15]+=r[7];X[XO17]+=r[8];X[XO18]+=r[9];X[XF19]+=r[10];X[XNe21]+=r[11];X[XNe22]+=r[12];X[XNa22]+=r[13];X[XNa23]+=r[14];X[XMg24]+=r[15];X[XMg25]+=r[16];X[XMg26]+=r[17];X[XAl26]+=r[18]+r[19]; X[XAl27]+=r[20];X[XSi28]+=r[21];X[XSi29]+=r[22];X[XSi30]+=r[23];X[XP31]+=r[24];X[XS32]+=r[25];X[XS33]+=r[26];X[XS34]+=r[27];X[XFe56]+=r[28];X[XFe57]+=r[29];X[XFe58]+=r[30]; + +//X[XFe59]+=r[31];X[XFe60]+=r[32]; + +X[XCo59]+=r[33];X[XNi58]+=r[34]; +//X[XNi59]+=r[35]; +X[XNi60]+=r[36];X[XNi61]+=r[37];X[XH1]+=r[38];X[XHe4]+=r[39];X[XC12]+=r[40];X[XN14]+=r[41];X[XO16]+=r[42];X[XNe20]+=r[43]; Dprint("# Richard's table (modified) final X : H1=%g He4=%g C12=%g C13=%g N14=%g O16=%g\n", X[XH1],X[XHe4],X[XC12],X[XC13],X[XN14],X[XO16]); @@ -507,7 +515,22 @@ void nucsyn_first_dup_apply(double m,double z, Dprint("# Amanda : DH1=%g DHe4=%g DC12=%g DC13=%g DN14=%g DO16=%g\n", r[0],r[3],r[7],r[8],r[11],r[15]); - X[XH1]+=r[0]; X[XH2]+=r[1]; X[XHe3]+=r[2]; X[XHe4]+=r[3];X[XLi7]+=r[4]; X[XLi7]=MAX(0.0,X[XLi7]); X[XBe7]+=r[5]; X[XBe7]=MAX(0.0,X[XBe7]); X[XC12]+=r[7]; X[XC13]+=r[8]; X[XN13]+=r[10]; X[XN14]+=r[11]; X[XN15]+=r[12]; X[XO15]+=r[14]; X[XO16]+=r[15]; X[XO17]+=r[16]; X[XO18]+=r[17]; X[XF17]+=r[19]; X[XF18]+=r[20]; X[XF19]+=r[21]; X[XNe20]+=r[24]; X[XNe21]+=r[25]; X[XNe22]+=r[26]; X[XNa22]+=r[28]; X[XNa23]+=r[29]; X[XNa24]+=r[30]; X[XMg24]+=r[32]; X[XMg25]+=r[33]; X[XMg26]+=r[34]; X[XAl26]+=r[37]; X[XAl27]+=r[39]; X[XAl28]+=r[40]; X[XSi28]+=r[42]; X[XSi29]+=r[43]; X[XSi30]+=r[44]; X[XP30]+=r[49]; X[XP31]+=r[50]; X[XS32]+=r[54]; X[XS33]+=r[55]; X[XS34]+=r[56]; X[XS35]+=r[57]; X[XFe56]+=r[58]; X[XFe57]+=r[59]; X[XFe58]+=r[60]; X[XFe59]+=r[61]; X[XFe60]+=r[62]; X[XCo59]+=r[64]; X[XCo60]+=r[65]; X[XCo61]+=r[66]; X[XNi58]+=r[67]; X[XNi59]+=r[68]; X[XNi60]+=r[69]; X[XNi61]+=r[70]; + X[XH1]+=r[0]; X[XH2]+=r[1]; X[XHe3]+=r[2]; X[XHe4]+=r[3];X[XLi7]+=r[4]; X[XLi7]=MAX(0.0,X[XLi7]); X[XBe7]+=r[5]; X[XBe7]=MAX(0.0,X[XBe7]); X[XC12]+=r[7]; X[XC13]+=r[8]; X[XN13]+=r[10]; X[XN14]+=r[11]; X[XN15]+=r[12]; X[XO15]+=r[14]; X[XO16]+=r[15]; X[XO17]+=r[16]; X[XO18]+=r[17]; X[XF17]+=r[19]; X[XF18]+=r[20]; X[XF19]+=r[21]; X[XNe20]+=r[24]; X[XNe21]+=r[25]; X[XNe22]+=r[26]; X[XNa22]+=r[28]; X[XNa23]+=r[29]; +// X[XNa24]+=r[30]; + + X[XMg24]+=r[32]; X[XMg25]+=r[33]; X[XMg26]+=r[34]; X[XAl26]+=r[37]; X[XAl27]+=r[39]; +//X[XAl28]+=r[40]; + + X[XSi28]+=r[42]; X[XSi29]+=r[43]; X[XSi30]+=r[44]; X[XP30]+=r[49]; X[XP31]+=r[50]; X[XS32]+=r[54]; X[XS33]+=r[55]; X[XS34]+=r[56]; +//X[XS35]+=r[57]; + X[XFe56]+=r[58]; X[XFe57]+=r[59]; X[XFe58]+=r[60]; +//X[XFe59]+=r[61]; +// X[XFe60]+=r[62]; + X[XCo59]+=r[64]; +// X[XCo60]+=r[65]; X[XCo61]+=r[66]; + X[XNi58]+=r[67]; +//X[XNi59]+=r[68]; + X[XNi60]+=r[69]; X[XNi61]+=r[70]; } #endif // NUCSYN_FIRST_DREDGE_UP_AMANDAS_TABLE diff --git a/src/nucsyn/nucsyn_set_envelope_abunds.c b/src/nucsyn/nucsyn_set_envelope_abunds.c index 84816c3a7..3216fe4f1 100644 --- a/src/nucsyn/nucsyn_set_envelope_abunds.c +++ b/src/nucsyn/nucsyn_set_envelope_abunds.c @@ -28,8 +28,6 @@ void nucsyn_set_envelope_abunds(struct star_t * RESTRICT star, memcpy(star->profileX[i],X,ISOTOPE_MEMSIZE); } } - - } #endif//NUCSYN_PROFILES diff --git a/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c b/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c index 804082635..7ebbc48d7 100644 --- a/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c +++ b/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c @@ -29,6 +29,8 @@ void nucsyn_set_third_dredgeup_abunds(Abundance * RESTRICT dup_material, * set isotopic abundances of dredgeup material according to Amanda's paper * (Karakas et al. 2002). The array dup_material contains the abundances, *star * is provided so that we can extend these to be time dependent + * + * Table format is : Z,M,Npulse */ static const double intershell_table[NUCSYN_INTERSHELL_TABLE_LINES*NUCSYN_INTERSHELL_TABLE_LINE_LENGTH]={NUCSYN_INTERSHELL_TABLE}; double x[3]; // input data @@ -96,8 +98,8 @@ void nucsyn_set_third_dredgeup_abunds(Abundance * RESTRICT dup_material, dup_material[XSi30]=intershell[15]; dup_material[XP31]=intershell[16]; - - Dprint("PULSE %g x=%g %g %g set dup material He4=%g C12=%g Ne20=%g Ne22=%g Mg24=%g Mg25=%g Mg26=%g\n", + /* + printf("PULSE %g x=%g %g %g set dup material He4=%g C12=%g Ne20=%g Ne22=%g Mg24=%g Mg25=%g Mg26=%g\n", star->num_thermal_pulses, x[0],x[1],x[2], dup_material[XHe4], @@ -107,7 +109,7 @@ void nucsyn_set_third_dredgeup_abunds(Abundance * RESTRICT dup_material, dup_material[XMg24], dup_material[XMg25], dup_material[XMg26]); - + */ #ifdef NUCSYN_STPAGB diff --git a/src/nucsyn/nucsyn_short_log.c b/src/nucsyn/nucsyn_short_log.c index fe0927cb0..25a46472f 100644 --- a/src/nucsyn/nucsyn_short_log.c +++ b/src/nucsyn/nucsyn_short_log.c @@ -284,9 +284,7 @@ void nucsyn_short_log(struct stardata_t * RESTRICT stardata) stardata->star[i].mass ); -#if (DEBUG==1) - fflush(stdout); -#endif + } /* s-process */ @@ -400,7 +398,6 @@ void nucsyn_short_log(struct stardata_t * RESTRICT stardata) S_HS(stardata->star[i].Xsurf), S_LS(stardata->star[i].Xsurf), S_HS_LS(stardata->star[i].Xsurf)); - fflush(stdout); } #endif diff --git a/src/nucsyn/nucsyn_sn_core_collapse.c b/src/nucsyn/nucsyn_sn_core_collapse.c index a4eb87a97..3d3a2a9d0 100644 --- a/src/nucsyn/nucsyn_sn_core_collapse.c +++ b/src/nucsyn/nucsyn_sn_core_collapse.c @@ -436,14 +436,12 @@ void nucsyn_sn_core_collapse(Abundance *X,//space calloc'd previously Dprint("Use Chieffi and Limongi 2004 yields\n"); #endif - fflush(stdout); Dprint("SN mco=%g menv=%g : explosive yields (M=%g) : C=%g O=%g Fe56=%g (Xtot=%g)\n", mco,menv,mco+menv, dm_core*X[XC12]/Xtot, dm_core*X[XO16]/Xtot, dm_core*X[XFe56]/Xtot, Xtot); - fflush(stdout); #endif nucsyn_calc_yields(stardata, diff --git a/src/nucsyn/nucsyn_third_dredge_up.c b/src/nucsyn/nucsyn_third_dredge_up.c index f7810a9ce..83ceffd5a 100644 --- a/src/nucsyn/nucsyn_third_dredge_up.c +++ b/src/nucsyn/nucsyn_third_dredge_up.c @@ -8,7 +8,6 @@ #ifdef NUCSYN_THIRD_DREDGE_UP_AMANDA - /* growth rate of lambda (from Onno) */ #ifdef NUCSYN_THIRD_DREDGE_UP_TABULAR_NCAL /* use tabular fit for NCAL - use this */ @@ -41,10 +40,6 @@ else\ #define MIN_MASS_4_3DUP nucsyn_mcmin_for_dredgeup_rs(*mt,stardata->common.metallicity,stardata->common.delta_mcmin) #endif // NUCSYN_THIRD_DREDGE_UP_RICHARD -/* Conditions for third dredge up */ - -#define TDUP_CONDITIONS ((*mc >= MAX(star->mc0,MIN_MASS_4_3DUP))&&(menv>MINIMUM_ENVELOPE_MASS_FOR_THIRD_DREDGEUP)) - /************************************************************/ /****************** Third dredge up *************************/ /************************************************************/ @@ -199,8 +194,18 @@ double nucsyn_third_dredge_up(struct star_t *star, /* calculate lambda (3DUP efficiency) - this is returned at the end */ CALCLAMBDA; - if(TDUP_CONDITIONS) - { + /* + printf("lambda=%g lammax=%g ncal=%g L=%g dmc_prev_pulse=%g\n",lambda,lammax,ncal, + star->luminosity,star->dmc_prev_pulse); + printf("check mc=%g>%g and menv=%g>%g (%g)\n", + *mc,MAX(star->mc0,MIN_MASS_4_3DUP), + menv,MINIMUM_ENVELOPE_MASS_FOR_THIRD_DREDGEUP, + stardata->preferences->minimum_envelope_mass_for_third_dredgeup); + */ + + if((*mc >= MAX(star->mc0,MIN_MASS_4_3DUP))&& + (menv>MINIMUM_ENVELOPE_MASS_FOR_THIRD_DREDGEUP)) + { /* Third dredgeup has occurred */ star->num_thermal_pulses_since_mcmin++; @@ -226,8 +231,7 @@ double nucsyn_third_dredge_up(struct star_t *star, memset(dup_material,0,ISOTOPE_MEMSIZE); nucsyn_set_third_dredgeup_abunds(dup_material,star,stardata->common.metallicity,stardata,*mc); - - Dprint("DUPB mix dredge up material masses: intershell %g (H=%g He=%g C=%g), envelope: %g (H=%g He=%g C=%g)\n", + Dprint("DUPB mix dredge up material masses: intershell %g (H=%g He=%g C=%g), envelope: %g (H=%g He=%g C=%g)\n", dm_dredge_up, dup_material[XH1], dup_material[XHe4], diff --git a/src/nucsyn/nucsyn_update_abundances.c b/src/nucsyn/nucsyn_update_abundances.c index 68181efbe..820ae5d4b 100644 --- a/src/nucsyn/nucsyn_update_abundances.c +++ b/src/nucsyn/nucsyn_update_abundances.c @@ -133,7 +133,7 @@ void nucsyn_update_abundances(struct stardata_t * RESTRICT stardata) nXenv[k]=malloc(ISOTOPE_MEMSIZE); #ifdef NUCSYN_PROFILES - printf("Update : profile alloced\n");fflush(stdout); + printf("Update : profile alloced\n"); #endif #ifdef ALLOC_CHECKS diff --git a/src/nucsyn/nucsyn_update_abundances_with_profiles.c b/src/nucsyn/nucsyn_update_abundances_with_profiles.c index 69001f708..84b09df72 100644 --- a/src/nucsyn/nucsyn_update_abundances_with_profiles.c +++ b/src/nucsyn/nucsyn_update_abundances_with_profiles.c @@ -154,9 +154,7 @@ void nucsyn_update_abundances_with_profiles(struct stardata_t * RESTRICT stardat #ifdef NUCSYN_PROFILES_DEBUG printf("Donor: H=%g C=%g; ",Xdonor[XH1],Xdonor[XC12]); - fflush(stdout); printf("Accretor: H=%g C=%g; ",Xaccretor[XH1],Xaccretor[XC12]); - fflush(stdout); #endif // Mix the winds in-system (to X) diff --git a/src/perl/install_modules.pl b/src/perl/install_modules.pl index 8e24d7c95..54307aaf5 100755 --- a/src/perl/install_modules.pl +++ b/src/perl/install_modules.pl @@ -6,14 +6,20 @@ use 5.16.0; -foreach my $mod (#'local::lib', - 'common::sense', - 'Compress::Zlib', +foreach my $mod ( 'Carp', 'Clone', - #'Config', + 'common::sense', + 'Compress::Zlib', + 'Cwd', + 'Data::Dumper', + 'Data::Serializer', 'Fcntl', + 'File::Basename', + 'FindBin', 'IO::File', + 'IO::Handle', + 'IO::Interactive', 'IPC::Open3', 'Sort::Key', 'Sys::Hostname', @@ -29,13 +35,7 @@ foreach my $mod (#'local::lib', 'constant', 'List::Util', 'Memoize', - #'vars', - #'subs', - #'strict' - ) -{ - install($mod); -} + ){ install($mod); } my $pwd=$ENV{PWD}; diff --git a/src/perl/modules/binary_grid.pm b/src/perl/modules/binary_grid.pm index e02411786..9c32ca1de 100644 --- a/src/perl/modules/binary_grid.pm +++ b/src/perl/modules/binary_grid.pm @@ -36,7 +36,7 @@ use common::sense; # many variables gone, added snapshot ability, new Perl 5.16 # syntax introduced where appropriate -our $VERSION = '0.15'; +our $VERSION = '0.16'; $|=1; # autoflush please (for log output) @@ -47,11 +47,13 @@ use Carp qw(confess); use Clone qw(clone); use Config; # for signals use Fcntl; # for file (un)blocking +use Fcntl qw(:flock SEEK_END); # import LOCK_* and SEEK_END constants use IO::File; use IPC::Open3; use Sort::Key qw(nsort); use Sys::Hostname; use Time::HiRes qw(sleep gettimeofday tv_interval); +use Cwd; ## threads modules and friends use threads; @@ -76,7 +78,7 @@ require Exporter; our @ISA = qw(Exporter); # Functions to export -my @funcs= qw(&flexigrid &grid_defaults &parse_grid_args &setup_binary_grid &tbse_kill &tbse_line &isotope_list &isotope_hash &source_list &tbse_launch &tbse_land &tbse &nuclear_mass_list &nuclear_mass_hash &get_from_header &show_options &options_string &binary_grid_mem_usage &get_bsegridvar &set_bsegridvar &ensemble_list &stellar_type_strings +my @funcs= qw(&flexigrid &grid_defaults &parse_grid_args &setup_binary_grid &tbse_kill &tbse_line &isotope_list &isotope_hash &source_list &tbse_launch &tbse_land &tbse &nuclear_mass_list &nuclear_mass_hash &get_from_header &show_options &options_string &binary_grid_mem_usage &get_bsegridvar &set_bsegridvar &ensemble_list &stellar_type_strings &merge_results_hash_dumps ); our %EXPORT_TAGS = ( 'all' => [ @funcs ] ); our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } ); @@ -115,6 +117,8 @@ my %threadinfo; # Subroutines follow ############################################################ + +my $joincount=0; use subs qw(print_binary_in binary_out verbosity); sub verbosity @@ -319,9 +323,27 @@ sub grid_defaults 'snapshot_file'=>$tmp.'/binary_c-snapshot', + # if arg_checking is on, each argument is checked for validity: + # this involves a lot of I/O and is slow: don't do it unless you have crashes + 'arg_checking'=>0, + + ######################################## + # infomration about the running grid script + ######################################## + 'working_directory'=>cwd(), # the starting directory + 'perlscript'=>$0, # the name of the perlscript + 'perlscript_arguments'=>join(' ',@ARGV), # arguments as a string + 'perl_executable'=>$^X, # the perl executable + 'command_line'=>join(' ',$0,@ARGV), # full command line + ######################################## # GRID ######################################## + 'rungrid'=>1, # usually run the grid, but can be 0 + # to skip it (e.g. for condor runs) + + 'merge_datafiles'=>'', + 'binary'=>0, # set to 0 for single stars, 1 for binaries # where is binary_c? need this to get the values of some counters @@ -362,12 +384,29 @@ sub grid_defaults modulo=>1, offset=>0, + # data dump file : undef by default (do nothing) + results_hash_dumpfile => '', + + # put binary_c's output in a cache and process it afterward: + # might be quicker... ? (at least should prevent some CPU + # switching from one task to another) + 'cache_binary_c_output'=>1, + + ######################################## + # Condor stuff + ######################################## + 'condor_command'=>'',# condor command e.g. "run_flexigrid", "join_datafiles" + 'condor_dir'=>'', # working directory containing e.g. scripts, output, logs (e.g. should be NFS available to all) + 'condor_njobs'=>'', # number of scripts + 'condor_jobid'=>'', # condor job id + ######################################## # Code, Timeouts, Signals ######################################## 'timeout'=>30, # timeout in seconds 'timeout_vb'=>0, # no timeout logging 'tvb'=>0, # no thread logging + 'nfs_sleep'=>1, # time to wait for NFS to catch up with file accesses # flexigrid checks the timeouts every # flexigrid_timeout_check_interval seconds @@ -614,7 +653,7 @@ sub options_string }keys %ENV; # and the binary_c defined options - my $cmd=$binary_grid::grid_options{'prog'}.' --version'; + my $cmd=binary_c_prog().' --version'; my @r=`$cmd`; map { @@ -624,7 +663,6 @@ sub options_string return join('',@s); } - sub setup_sig_handlers { if(! $binary_grid::grid_options{'no signals'}) @@ -679,7 +717,7 @@ sub shutdown_binary_grid stop_flexigrid_threads(); # if this fails, KILL - sleep 1; + sleep 1 if($binary_grid::grid_options{condor_jobid} eq ''); kill_flexigrid_binary_c_pids(); # kill main process with given signal @@ -689,14 +727,14 @@ sub shutdown_binary_grid if(fork()==0) { say "Child waiting to clean up with signal $sig..."; - sleep 5; + sleep 5 if($binary_grid::grid_options{condor_jobid} eq ''); print "Everything died because:\n",color('red bold'),$flexigrid_shared{error},"\n",(defined $flexigrid_shared{global_error_string} ? $flexigrid_shared{global_error_string} : ''), color('reset') if(defined($flexigrid_shared{error})); kill $sig => $pid; exit; } } - sleep 1; + sleep 1 if($binary_grid::grid_options{condor_jobid} eq ''); map{ printf "Send SIGSTOP to thread %d\n",$_->tid if(verbosity); @@ -704,7 +742,7 @@ sub shutdown_binary_grid }(threads->list(threads::running)); say 'Threads killed'; - sleep 1; + sleep 1 if($binary_grid::grid_options{condor_jobid} eq ''); say 'Try tbse_land, tbse_kill'; tbse_land(); @@ -774,7 +812,7 @@ sub tbse_restart say STDERR 'Trying to restart binary_c and then continue'; tbse_land(); tbse_kill(); - sleep 1; + sleep 1 if($binary_grid::grid_options{condor_jobid} eq ''); tbse_launch(); say STDERR 'Restart complete'; } @@ -831,6 +869,11 @@ sub tbse_kill } +sub binary_c_prog +{ + return -e $binary_grid::grid_options{'rootpath'}.'/'.$binary_grid::grid_options{'prog'} ? $binary_grid::grid_options{'rootpath'}.'/'.$binary_grid::grid_options{'prog'} : $binary_grid::grid_options{'prog'}; +} + sub tbse_launch { ############################################################ @@ -840,7 +883,7 @@ sub tbse_launch # choose executable for binary_c : usually it is specified, but if not, # try some defaults - my $prog= -e $binary_grid::grid_options{'rootpath'}.'/'.$binary_grid::grid_options{'prog'} ? $binary_grid::grid_options{'rootpath'}.'/'.$binary_grid::grid_options{'prog'} : $binary_grid::grid_options{'prog'}; + my $prog=binary_c_prog(); my $command; # command to be set up @@ -893,7 +936,7 @@ sub tbse_launch # logging my $s=sprintf "Thread %d :binary_c pid %d \n", - $threadinfo{binary_c_pid},$threadinfo{thread_number}; + $threadinfo{thread_number},$threadinfo{binary_c_pid}; thread_log($s,$threadinfo{thread_number}); @@ -902,7 +945,7 @@ sub tbse_launch say '# Bi-directional binary_c/nucsyn pipe open for business' if(verbosity); # sleep to allow threads to kick in - sleep $binary_grid::grid_options{'thread_presleep'} if(defined($binary_grid::grid_options{'thread_presleep'})); + sleep $binary_grid::grid_options{'thread_presleep'} if(defined($binary_grid::grid_options{'thread_presleep'}) && ($binary_grid::grid_options{condor_jobid} eq '')); return $threadinfo{binary_c_pid}; } @@ -922,7 +965,7 @@ sub tbse # in as arguments. # It sends output the given parser function. ############################################################ - + # set timeout in case of failure/pause/lockup if($binary_grid::grid_options{'timeout'} && @@ -1001,11 +1044,10 @@ sub tbse # save the previous alive time for each thread $threadinfo{thread_prev_alive}=time(); - vbout(3, "About to log args\n"); - # log arguments to a file if($binary_grid::grid_options{'log_args'}) { + vbout(3, "About to log args\n"); my $f = $binary_grid::grid_options{'log_args_dir'}."/binary_c-args.thread".$threadinfo{thread_number}; @@ -1038,65 +1080,78 @@ sub tbse flush_binary_out(); blocking(); # enforce blocking - - # send arguments one by one, check that they work - map - { - my $arg='--'.$_; - vbout(3,"Send arg $arg .. \n"); - # send arg - print_binary_in $arg."\n"; - vbout(3,"sent $arg : checking response\n"); - - # get response, was it ok? - my ($x1,$x2,$x); - - # 3 attempts to get responses: each time we sleep for a second - # to let the process catch up. That should be enough time, although - # of course on some machines perhaps not... but then the binary_out - # call should block. Why it doesn't (sometimes) I don't know. - my $failcount=0; - while($failcount<3) + if($binary_grid::grid_options{arg_checking}) + { + # send arguments one by one, check that they work + map { - # try to get a line of output - $x1=binary_out; - - # ignore 'set batchmode to 1' - $x2=binary_out if(defined($x1) && ($x1 eq "set batchmode to 1\n")); + my $arg='--'.$_; + vbout(3,"Send arg $arg .. \n"); - # choose which to use: x2 if we have it (is newer) - $x= $x2 // $x1; + # send arg + print_binary_in $arg."\n"; + vbout(3,"sent $arg : checking response\n"); - if(defined($x)) + # get response, was it ok? + my ($x1,$x2,$x); + + # 3 attempts to get responses: each time we sleep for a second + # to let the process catch up. That should be enough time, although + # of course on some machines perhaps not... but then the binary_out + # call should block. Why it doesn't (sometimes) I don't know. + my $failcount=0; + while($failcount<3) { - last; + # try to get a line of output + $x1=binary_out; + + # ignore 'set batchmode to 1' + $x2=binary_out if(defined($x1) && ($x1 eq "set batchmode to 1\n")); + + # choose which to use: x2 if we have it (is newer) + $x= $x2 // $x1; + + if(defined($x)) + { + last; + } + else + { + # sleep and try again (warn after first failure, let that slip by into the ether) + say STDERR "warning: thread $threadinfo{thread_number} sleeping (awaiting confirmation of argument '$arg' from binary_c: x1='$x1', x2='$x2')" if($failcount>0); + $failcount++; + sleep 1 if($binary_grid::grid_options{condor_jobid} eq ''); + } } - else + + if((!(defined($x) && $x=~/arg ok/o)) || ($failcount>=3)) { - # sleep and try again (warn after first failure, let that slip by into the ether) - say STDERR "warning: thread $threadinfo{thread_number} sleeping (awaiting confirmation of argument '$arg' from binary_c: x1='$x1', x2='$x2')" if($failcount>0); - $failcount++; - sleep 1; - } - } + # error detected : close the pipes before damage is done + say STDERR "ARG $arg failed (output $x)"; - if((!(defined($x) && $x=~/arg ok/o)) || ($failcount>=3)) - { - # error detected : close the pipes before damage is done - say STDERR "ARG $arg failed (output $x)"; + close BINARY_IN; + close BINARY_OUT; + close BINARY_ERR; - close BINARY_IN; - close BINARY_OUT; - close BINARY_ERR; + waitpid $threadinfo{binary_c_pid},0; - waitpid $threadinfo{binary_c_pid},0; + return("Argument to binary_c '$arg' failed (received x='$x', x1='$x1', x2='$x2' from binary_c instead) : check that it is correct in binary_c or in your grid script, (\%bse_options) and also check that binary_c is running ok. (failcount=$failcount)\n"); + } + vbout(3,"response to '$arg' checked and ok\n"); + + }grep {$_ ne ''} split(/--/o,$args); - return("Argument to binary_c '$arg' failed (received x='$x', x1='$x1', x2='$x2' from binary_c instead) : check that it is correct in binary_c or in your grid script, (\%bse_options) and also check that binary_c is running ok. (failcount=$failcount)\n"); - } - vbout(3,"response to '$arg' checked and ok\n"); - - }grep {$_ ne ''} split(/--/o,$args); + } + else + { + # no arg checking: just send them + vbout(3,"Send args $args .. \n"); + map + { + print_binary_in '--'.$_."\n"; + }grep {$_ ne ''} ('no_arg_checking',split(/--/o,$args)); + } vbout(3,"binary_c go\n"); @@ -1234,17 +1289,43 @@ sub tbse_line { ############################################################ # get a line from the binary_c/nucsyn process, check - # for some special cases and then return it + # for some special cases and then return it. + # Note that the line returned is guaranteed to be chomped. ############################################################ # (perhaps) flush BINARY_ERR stream before we do anything flush_binary_err() if($binary_grid::grid_options{'always flush binary_err'}); - my $x=binary_out; # the data line + my $x; + state @cache; + + if($binary_grid::grid_options{cache_binary_c_output}) + { + # if we should cache output, save it all now + if($#cache==-1) + { + while($x=binary_out) + { + push(@cache,$x); + last if($x eq "fin\n"); + } + chomp @cache; + vbout(3,"Filled binary_c cache with $#cache lines\n"); + } + # shift off the first cache line + $x=shift(@cache); + } + else + { + # no caching : get output directly + $x=binary_out; # the data line + chomp $x; + } + if(defined($x)) { - return 'fin' if($x=~/^Tick count/o); - + #return 'fin' if($x=~/^Tick count/o); + # detect if the thread has become stuck in an # infinite loop of doom @@ -1278,7 +1359,6 @@ sub tbse_line } # then just return the data to the grid script - chomp $x; return $x; } else @@ -1331,10 +1411,10 @@ sub increment_probtot sub increment_count { # increment counter: NB this MUST be a function in binary_grid - print "INC COUNT\n"; $binary_grid::flexigrid{count}++; } + sub vb1out { ############################################################ @@ -1370,27 +1450,34 @@ sub vb1out } chomp $etf; + state $hostname = hostname(); + my $pstring=' '; { my $vb1string=' % '.(length($ngrid)+1).'d/%d %s%s%s %s -> '; # output to the terminal (stdout) - my $s=sprintf " %s$vb1string%s%s % 5.1f %% complete%s ETA=% 7.1f $units tpr=%2.2e %sETF=%s %s %s pid=$$ mem:grd=%d fam=%d (kids:%3.2f%%) %s", + my $s=sprintf " %s$vb1string%s%s % 5.1f %% complete%s ETA=% 7.1f $units tpr=%2.2e %sETF=%s %s %s pid=$$ mem:grd=%d fam=%d (kids:%3.2f%%) %s %s", color("yellow bold on_black"), $binary_grid::flexigrid{count},$ngrid, color('green'),thetime(), color('yellow'),$t, - color("blue on_white"),$pstring,100.0*$binary_grid::flexigrid{count}/$ngrid, - color("bold red on_black"),$eta,$tpr, + color('blue on_white'),$pstring,100.0*$binary_grid::flexigrid{count}/$ngrid, + color('bold red on_black'),$eta,$tpr, color('blue'),$etf, color('magenta'), defined($binary_grid::grid_options{'binary_c_pid'}) ? 'pid='.$binary_grid::grid_options{'binary_c_pid'}.' ' : ' ', binary_grid_mem_usage(), binary_grid_family_mem_usage(), 100.0*frac_mem_usage(), + $binary_grid::grid_options{condor_jobid} eq '' ? '' : color('magenta')."\@$hostname (job=$binary_grid::grid_options{condor_jobid}) ", color('reset'); + $s=~s/weight=1//; # trivial if true + $s=~s/\(kids\:0.00\%\)//; + $s=~s/\s+/ /g; # save space + thread_log($s,$binary_grid::threadinfo{thread_number}); thread_vb($s) if verbosity(3); } @@ -1670,89 +1757,135 @@ sub flexigrid # Takes one argument: the number of threads # ############################################################ - - # if there is no flexigrid hash in existence (it may have been loaded from a snapshot) - # then make a new one - if(!defined($flexigrid{id})) - { - # set up new flexigrid hash (all threads have a copy of this, although - # they must not change it) - %flexigrid=( - count=>0, - global_error_string=>undef, - id=>random_string(8), # random string to ID the flexigrid - maxq=>$binary_grid::grid_options{maxq_per_thread}*$_[0], # total maximum queue length before we pause - modulo=>$binary_grid::grid_options{modulo}, # run modulo n - monitor_files=>[], - nextlogtime=>0, - nthreads=>$_[0], # number of shreads (passed in) - offset=>$binary_grid::grid_options{offset}, # start at model offset (0-based, so first model is zero) - probtot=>0.0, - resolution=>{ - shift=>0, - previous=>0, - n=>{} # per-variable resolution - }, - results_hash=>$binary_grid::grid_options{results_hash}, # set in your grid script - t0=>[gettimeofday], # flexigrid start time - thread_q=>undef, - threads=>undef, # array of threads objects - tstart=>time() # flexigrid start time - ); - convert_functions_to_pointers(); - } - $flexigrid_shared{error}=0; - if(!defined($flexigrid{results_hash})) + if($binary_grid::grid_options{rungrid}==1) { - confess("results_hash is undefined in flexigrid: please make sure \$binary_grid::grid_options\{results_hash\} points to an anonymous array"); - } + # if there is no flexigrid hash in existence (it may have been loaded from a snapshot) + # then make a new one + if(!defined($flexigrid{id})) + { + # set up new flexigrid hash (all threads have a copy of this, although + # they must not change it) + %flexigrid=( + count=>0, + global_error_string=>undef, + id=>random_string(8), # random string to ID the flexigrid + maxq=>$binary_grid::grid_options{maxq_per_thread}*($_[0]//1), # total maximum queue length before we pause + modulo=>$binary_grid::grid_options{modulo}, # run modulo n + monitor_files=>[], + nextlogtime=>0, + nthreads=>($_[0]//1), # number of shreads (passed in) + offset=>$binary_grid::grid_options{offset}, # start at model offset (0-based, so first model is zero) + probtot=>0.0, + resolution=>{ + shift=>0, + previous=>0, + n=>{} # per-variable resolution + }, + results_hash=>$binary_grid::grid_options{results_hash}, # set in your grid script + t0=>[gettimeofday], # flexigrid start time + thread_q=>undef, + threads=>undef, # array of threads objects + tstart=>time() # flexigrid start time + ); + convert_functions_to_pointers(); + } + $flexigrid_shared{error}=0; + $binary_grid::flexigrid{count}=0; + if(!defined($flexigrid{results_hash})) + { + confess("results_hash is undefined in flexigrid: please make sure \$binary_grid::grid_options\{results_hash\} points to an anonymous array"); + } - # erase arrays - @flexigrid_binary_c_pids=(); + # erase arrays + @flexigrid_binary_c_pids=(); - my $starcount; - - # check/set threads' stack size - set_threads_stacksize(); - - # launch threads - vbout(0,"Launch $flexigrid{nthreads} threads\n"); + my $starcount; + + # check/set threads' stack size + set_threads_stacksize(); + + # launch threads + vbout(0,"Launch $flexigrid{nthreads} threads\n"); - ($flexigrid{threads},$flexigrid{thread_q})= - launch_flexigrid_threads($flexigrid{nthreads}); + ($flexigrid{threads},$flexigrid{thread_q})= + launch_flexigrid_threads($flexigrid{nthreads}); - # set grid option in case external programs want to know the number of threads - $binary_grid::grid_options{'nthreads'}=$flexigrid{nthreads}; + # set grid option in case external programs want to know the number of threads + $binary_grid::grid_options{'nthreads'}=$flexigrid{nthreads}; - vbout(1,"Make flexigrid:\n"); - - # make flexigrid code on the fly while we wait for the threads to start - my $gridcode=make_flexigrid(); - - vbout(1,"Run flexigrid\n"); + vbout(1,"Make flexigrid:\n"); + + # make flexigrid code on the fly while we wait for the threads to start + my $gridcode=make_flexigrid(); + + vbout(1,"Run flexigrid\n"); - # set up timeout - if($binary_grid::grid_options{'timeout'} && - (!$binary_grid::grid_options{'disable signal'}{'ALRM'})) - { - Time::HiRes::alarm($binary_grid::grid_options{'flexigrid_timeout_check_interval'}); + # set up timeout + if($binary_grid::grid_options{'timeout'} && + (!$binary_grid::grid_options{'disable signal'}{'ALRM'})) + { + Time::HiRes::alarm($binary_grid::grid_options{'flexigrid_timeout_check_interval'}); + } + + # run the grid to send commands to the threads + my $gridcode_retval = eval $gridcode; + + # check for a return value (undef means error!) + flexigrid_eval_error($@,$gridcode) if(!defined($gridcode_retval)); + + vbout(1, "Grid finished: joining threads... (be patient)\n"); + + # join threads to get the data + my $local_threadinfo=join_flexigrid_threads($flexigrid{thread_q},$flexigrid{threads}); + + vbout(1, "DUMP ? results_hash_dumpfile=$binary_grid::grid_options{results_hash_dumpfile}\n"); + # dump results hash if required + dump_results_hash() + if(defined($binary_grid::grid_options{results_hash_dumpfile})&& + ($binary_grid::grid_options{results_hash_dumpfile} ne '')); + + # timing stats + printf "Flexigrid took %g seconds\n",tv_interval($flexigrid{t0},[gettimeofday]); + + if($binary_grid::grid_options{condor_jobid} ne '') + { + # end of a condor run : call the condor_job_hook function in the condor module + use binary_grid::condor; + binary_grid::condor::condor_job_hook(); + } } + else + { + # rungrid is zero: do something else - # run the grid to send commands to the threads - my $gridcode_retval = eval $gridcode; + # merge a set of datafiles, the result goes into results_hash, + # and is dumped if a results_hash_dumpfile is specified + if($binary_grid::grid_options{merge_datafiles}) + { + # the merge_datafiles variable should be a comma-separated list of filenames + # (without any spaces in it!) + my @datafiles=split(',',$binary_grid::grid_options{'merge_datafiles'}); + + say "Merge datafiles: @datafiles"; - # check for a return value (undef means error!) - flexigrid_eval_error($@,$gridcode) if(!defined($gridcode_retval)); + # now do the merge + $binary_grid::grid_options{results_hash}=merge_results_hash_dumps(@datafiles); - vbout(1, "Grid finished: joining threads... (be patient)\n"); - - # join threads to get the data - my $local_threadinfo=join_flexigrid_threads($flexigrid{thread_q},$flexigrid{threads}); - - # timing stats - printf "Flexigrid took %g seconds\n",tv_interval($flexigrid{t0},[gettimeofday]); + print "MERGED : $binary_grid::grid_options{results_hash}, dump to $binary_grid::grid_options{results_hash_dumpfile}\n"; + # save if dumpfile name is given + dump_results_hash() if($binary_grid::grid_options{results_hash_dumpfile} ne ''); + } + + # load the results_hash from a file (that was previously dumped) + if(defined($binary_grid::grid_options{load_results_from_file})) + { + # fill the results hash from data in a file : this allows us to process + # the data with a single grid script + $binary_grid::grid_options{results_hash}=load_results_hash_file($binary_grid::grid_options{'load_results_from_file'}); + } + } return; } @@ -1781,7 +1914,12 @@ sub wait_for_queue { # wait for $go to be signalled lock $flexigrid_lock; - cond_wait($flexigrid_lock); +# cond_wait($flexigrid_lock); + + # hmm cond_wait just locks condor: instead, we use a 1 second timeout + # so don't rely on this function for an empty queue slot! You must check + # that thread_q->pending() < $flexigrid{maxq} also. + cond_timedwait($flexigrid_lock,1); } sub flexigrid_timeout_check @@ -1805,9 +1943,11 @@ sub join_flexigrid_thread # You do NOT have to use this function, however it is recommended vbout(1,"Join hashes global results_hash=$binary_grid::grid_options{results_hash} and local contribution h=$_[0]\n"); - + arithmetically_add_hashes($binary_grid::grid_options{results_hash},$_[0]); + $joincount++; + return @_; } @@ -1905,7 +2045,6 @@ sub run_flexigrid_thread $binary_grid::grid_options{'running args'}=$args; my $err=tbse($args,$threadinfo{h}); # run a star - if(defined($err)) { say "binary_c error detected : $err"; @@ -2531,9 +2670,12 @@ TRAILING_BLOCK # final indentation $gridcode=indent_gridcode($gridcode); - # verbose output - output_gridcode_with_linenumbers($gridcode) ; - output_gridcode($gridcode); + if(!$binary_grid::grid_options{no_gridcode_dump}) + { + # verbose output + output_gridcode_with_linenumbers($gridcode) ; + output_gridcode($gridcode); + } return $gridcode; } @@ -2653,28 +2795,29 @@ sub run_binary_c { # Queue binary_c job with the given arguments - # but do nothing if modulo/offset combo is wrong - return if($binary_grid::flexigrid{count} < $binary_grid::grid_options{offset} || - ($binary_grid::flexigrid{count}+$binary_grid::grid_options{offset}) % $binary_grid::grid_options{modulo} != 0 - ); + # but do nothing if modulo/offset combo is wrong + return if(($binary_grid::flexigrid{count}+$binary_grid::grid_options{offset}) % $binary_grid::grid_options{modulo} != 0); # make arg string my $thread_q=shift @_; my $args=make_binary_c_argstring(@_); - vbout(2,"ENQUEUE $args on thread_q=$thread_q (currently ".$thread_q->pending()." pending, maxq=$flexigrid{maxq}\n"); + vbout(2,"thread_q: send to enqueue loop $args on thread_q=$thread_q (currently ".$thread_q->pending()." pending, maxq=$flexigrid{maxq}\n"); # wait for a queue spot my $brk=0; while($brk==0) { + #vbout(4,'thread_q pending pre-wait : pending=',$thread_q->pending()); if($thread_q->pending()<$flexigrid{maxq}) { + vbout(2,'thread_q: enqueue '.$args); $thread_q->enqueue($args); $brk=1; } else { + #vbout(4,'thread_q: wait_for_queue'); wait_for_queue(); } } @@ -2770,6 +2913,7 @@ sub join_flexigrid_threads my $f=$binary_grid::grid_options{'threads_join_function_pointer'}; vbout(2,"Join thread $j\n"); + my ($r,$threadinfo)=&$f(($j->join)); vbout(2,"Thread $j has been joined\n"); @@ -3025,21 +3169,16 @@ sub blocking sub isotope_list { - # obtain a list of isotopes from the binary_c code - check_grid_defaults(); - - my $dir=$binary_grid::grid_options{'srcpath'}.'/nucsyn/'; - my $file='nucsyn_isotopes.h'; - open(BINARY_GRID_ISOTOPE_LIST,"<$dir/$file")||confess("cannot find $dir/$file in binary_grid::isotope_list()"); + open(my $binc,binary_c_prog().' --version 2>&1 |')||confess("cannot open binary_c prog ".binary_c_prog()); my @isotopes; - - while(<BINARY_GRID_ISOTOPE_LIST>) + while(<$binc>) { - next if /(DECAY|TRACE)/o; - $isotopes[$2]=$1 if(/\#define X(\S+)\s+(\d+)/o); + if(/Isotope (\d+) is (\S+)/) + { + $isotopes[$1]=$2; + } } - - close BINARY_GRID_ISOTOPE_LIST; + close $binc; return @isotopes; } @@ -3141,6 +3280,8 @@ sub save_snapshot $binary_grid::flexigrid{actions}{save_snapshot}=0; } + + sub load_snapshot { ############################################################ @@ -3189,7 +3330,7 @@ sub load_snapshot # inherit the count and probtot, but not the rest of the flexigrid $binary_grid::flexigrid{count} = $$h{flexigrid}{count}; -$binary_grid::flexigrid{probtot} = $$h{flexigrid}{probtot}; + $binary_grid::flexigrid{probtot} = $$h{flexigrid}{probtot}; say "Set initial count to $binary_grid::flexigrid{count} from $$h{flexigrid}{count}, probtot to $binary_grid::flexigrid{probtot}"; # set the grid variables @@ -3233,6 +3374,65 @@ sub okdir return (-d $_[0]) ? $_[0] : undef; } + +sub dump_results_hash +{ + ############################################################ + # save the results_hash into a file, to be called after + # all threads are joined + ############################################################ + + # make data object + my $o=serializer_object(); + my $f=$_[0] // $binary_grid::grid_options{results_hash_dumpfile}; + say "dump_results_hash: data::serializer object $o, save to $f"; + + # dump to file : include some information other than simply the results + # just in case we want it + $o->store($binary_grid::grid_options{results_hash},$f) + ||confess("cannot write results_hash to $f"); + + #human-readable version + $f=$binary_grid::grid_options{results_hash_clean_dumpfile}; + if(defined($f) && $f ne '') + { + open(my $log_fp,'>'.$f)||confess("cannot open $f.datadump"); + print {$log_fp} Dumper($binary_grid::grid_options{results_hash}); + close $log_fp; + } +} + + + +sub load_results_hash_file +{ + # load results hash file (name passed in), return only the results_hash part + my $o=serializer_object(); + return $o->retrieve($_[0]); +} + + +sub merge_results_hash_dumps +{ + # merge results hash dumps with filenames + # given in @_ + + my $results_hash=$binary_grid::grid_options{results_hash} // {}; # returned + + # add the data for each file + foreach my $file (sort @_) + { + print "Merge data for file $file ... "; + my $h=load_results_hash_file($file); + print "add $h to $results_hash ..."; + arithmetically_add_hashes($results_hash,$h); + say 'done'; + } + + # return the completed results hash + return $results_hash; +} + 1; diff --git a/src/perl/modules/binary_grid/condor.pm b/src/perl/modules/binary_grid/condor.pm new file mode 100644 index 000000000..8ec5af848 --- /dev/null +++ b/src/perl/modules/binary_grid/condor.pm @@ -0,0 +1,503 @@ +package binary_grid::condor; +$|=1; +use 5.16.0; +use strict; +use IO::Handle; +use common::sense; + +######### binary_grid module for binary_c/nucsyn +######### >>> condor extensions <<< +# +# Written by Robert Izzard (robert.izzard@astro.uni-bonn.de +# or via gmail: rob.izzard) (c) 2000-2012 and onwards into the 21st century. + +# version number should match the binary_grid version +our $VERSION = '0.16'; + +$|=1; # autoflush please (for log output) + +# Standard modules: +# if you don't have them then check http://www.cpan.org/ or, better, +# use cpanm or packages to install them. +use Carp qw(confess); +use Clone qw(clone); +use Config; # for signals +use Fcntl; # for file (un)blocking +use IO::File; +use IPC::Open3; +use Sort::Key qw(nsort); +use Sys::Hostname; +use Time::HiRes qw(sleep gettimeofday tv_interval); +use File::Basename; +use rob_misc qw(mkdirhier); +use FindBin qw($Bin $Script $RealBin); +use lib "$Bin/../lib"; +use IO::Interactive; + +# module setup +require Exporter; +our @ISA = qw(Exporter); + +# Functions to export +my @funcs= qw( &condor_grid &condor_job_hook); + +our %EXPORT_TAGS = ( 'all' => [ @funcs ] ); +our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } ); +our @EXPORT = @funcs; +use vars qw(@ISA @EXPORT @EXPORT_OK); + +############################################################ +# Module-global variables +############################################################ + +############################################################ +# Subroutines follow +############################################################ + +# verbose logging +my $vb=1; # 0 = no logging, 1= logging +my $vbmirror=(! IO::Interactive::is_interactive()); # mirror stderr if not a terminal + +sub vb +{ + return if($vb==0); + + # set autoflush first time (just in case) + state $inited=0; + if($inited==0) + { + $inited=1; + STDOUT->autoflush(1); + STDERR->autoflush(1); + } + + # join args to make string, force newline at the end (but just one) + my $s="@_"; + chomp $s; + $s .= "\n"; + + # output + print STDOUT 'condor_grid(vb,stdout): ',$s; + print STDERR 'condor_grid(vb,stdout): ',$s if($vbmirror); +} + +sub condor_grid +{ + # grid wrapper for condor : should be called exactly like + # flexigrid() but instead distributes the work to condor jobs + vb("Bin $Bin : script $Script : RealBin $RealBin : perl executable $^X"); + + # disable timeouts and non-stdout/stderr logging + $binary_grid::grid_options{timeout}=0; + $binary_grid::grid_options{thread_presleep}=0; + $binary_grid::grid_options{nfs_sleep}=0; + $binary_grid::grid_options{log_args}=0; + $binary_grid::grid_options{tvb}=0; + $binary_grid::grid_options{no_gridcode_dump}=1; + $binary_grid::bse_options{log_filename}='/dev/null'; + + vb("condor_command is $binary_grid::grid_options{condor_command}"); + + if($binary_grid::grid_options{condor_command} eq 'run_flexigrid') + { + # in a job : just run flexigrid on one CPU + vb("call flexigrid(1)"); + return binary_grid::flexigrid(1); + } + elsif($binary_grid::grid_options{condor_command} eq 'join_datafiles') + { + vb("Joining data"); + # join data : don't run the grid + $binary_grid::grid_options{rungrid}=0; + # but do call flexigrid to join the data + return binary_grid::flexigrid(1); + } + else + { + # no condor command : we're in the main controller script + # so we have to set up the condor scripts + vb("In main controller thread: set up scripts"); + + # set up working directory + vb("Set workingdir()"); + condor_workingdir(); + + # how many scripts? default to 24 (max!) if nothing is given + my $njobs=$_[0] // $binary_grid::grid_options{condor_njobs} || 24; + say "Use $njobs condor jobs\n"; + vb("Use $njobs condor jobs"); + + # clean up joining files + vb("Unlink joining files"); + unlink($binary_grid::grid_options{condor_dir}.'/status/joining'); + + + # make scripts + vb("Make condor scripts"); + foreach my $n (1..$njobs) + { + my $opts={ + n=>$n, + njobs=>$njobs, + dir=>$binary_grid::grid_options{condor_dir}, + cmdline=>$binary_grid::grid_options{command_line}, + wd=>$binary_grid::grid_options{working_directory}, + }; + + my $scriptfile=make_condor_script($opts); + + vb("made script $n : submit"); + condor_submit($opts,$scriptfile); + } + + vb("Scripts submitted : all done"); + exit; # testing + } +} + + +sub condor_workingdir +{ + # define a working directory for condor : this is + $binary_grid::grid_options{condor_dir} = + $binary_grid::grid_options{condor_dir} || + ($binary_grid::grid_options{tmp}//'/tmp').'/condor/'.(rand()=~/\d\.(\d+)/)[0]; + + say "Condor working directory: $binary_grid::grid_options{condor_dir})"; + + # make directory tree + mkdirhier($binary_grid::grid_options{condor_dir}); + foreach ('scripts','stdout','stderr','results','logs','status') + { + mkdirhier($binary_grid::grid_options{condor_dir}.'/'.$_); + } + + return $binary_grid::grid_options{condor_dir}; +} + +sub make_condor_script +{ + say "Make condor_script $_[0]{n}"; + write_condor_script(condor_script_data(@_),@_); +} + +sub write_condor_script +{ + # write condor script to a file + my $scriptdata=shift @_; # data is always first arg + my $scriptfile=$binary_grid::grid_options{condor_dir}.'/scripts/'.jobid($_[0]); + open(my $fp,'>'.$scriptfile)||confess("can't open condor script $scriptfile for writing"); + print {$fp} $scriptdata; + close $fp; + say "Script written to $scriptfile"; + return $scriptfile; +} + +sub set_condor_job_status +{ + # set condor job status file (args=jobid,status string) + my $jobid=$_[0]; + my $status=$_[1]; + my $statusfile=$binary_grid::grid_options{condor_dir}.'/status/'.$jobid; + open(my $fp,'>'.$statusfile)||confess("cannot open status file $statusfile"); + print {$fp} $status; + close $fp; +} + +sub check_and_merge_condor_jobs +{ + # if all condor jobs have finished, merge the data and return it, + # otherwise return undef + if(check_condor_jobs_done(@_)) + { + return (merge_condor_jobs(@_)); + } + else + { + return (undef,undef); + } +} + +sub check_condor_jobs_done +{ + # check that all jobs in $path are finished, in which case return 1 + # otherwise return 0 + my $path=$_[0]; + my $njobs=$_[1]; + foreach (1..$njobs) + { + my $jobid=$_.'.'.$njobs; + my $f="$path/status/$jobid"; + if(open(my $fp,'<'.$f)) + { + if(!(<$fp> =~ /finished/)) + { + say "Job $jobid has not finished\n"; + return 0; + } + else + { + say "Job $jobid has finished\n"; + } + } + else + { + print "WARNING : Probable error : could not find status file $f - it should at least exist!\n"; + vb("WARNING : Probable error : could not find status file $f - it should at least exist!\n"); + return 0; + } + } + + say "All condor jobs have finished\n"; + + # no error, all jobs are finished : return 1 + return 1; +} + +sub merge_condor_jobs +{ + # merge all job data in $path + my $path=$_[0]; + my $njobs=$_[1]; + + # make file list + my @f; + foreach (1..$njobs) + { + my $jobid=$_.'.'.$njobs; + push(@f,"$path/results/$jobid"); + } + + # save + + # save data to file + my $f="$path/results/joined"; + binary_grid::dump_results_hash($f); + + # return data and filename + return (binary_grid::merge_results_hash_dumps(@_),$f); +} + + +sub condor_submit +{ + # submit a condor script + my $opts=$_[0]; + my $scriptfile=$_[1]; + vb( "Submit $scriptfile\n"); + `condor_submit $scriptfile`; + set_condor_job_status(jobid($opts),'submitted'); +} + +sub jobid +{ + # job identifier number : usually m.n where + # m is the job number and n the number of jobs in this group + return $_[0]{n}.'.'.$_[0]{njobs}; +} + +sub condor_script_data +{ + # return condor script data + + # the perl script + my $script=$_[0]{cmdline}; + $script=~s/condor_njobs=\d+//g; + $script=~s/\s+/ /g; $script=~s/\s+$//; + + # job id + my $id=jobid($_[0]); + + # set up environment + my @environment; + foreach ('PATH','PERL5LIB') + { + push(@environment, $_.'='.$ENV{$_}) if(defined($ENV{$_})); + } + my $environment = join(' ',@environment); + + my $s= " + ################################################# + # + # Condor script to run a binary_grid flexigrid + # + # script $_[0]{n} of $_[0]{njobs} + # + ################################################# + + executable = $^X + arguments = $script run_flexigrid=1 offset=".($_[0]{n}-1)." modulo=$_[0]{njobs} rungrid=1 results_hash_dumpfile=$_[0]{dir}/results/$id condor_command=run_flexigrid condor_jobid=$id condor_njobs=$_[0]{njobs} condor_dir=$_[0]{dir} vb=1 + environment = \"$environment\" + universe = vanilla + output = $_[0]{dir}/stdout/$id + error = $_[0]{dir}/stderr/$id + log = $_[0]{dir}/logs/$id + initialdir = $_[0]{wd} + remote_initialdir = $_[0]{wd} + stream_output = True + stream_error = True + +"; + + # optional condor parameters are given as condor_* in grid_options + foreach my $opt (grep {defined($binary_grid::grid_options{'condor_'.$_})} + ('request_memory') + ) + { + $s .= " $opt = ".$binary_grid::grid_options{'condor_'.$_}."\n"; + } + + $s.= " + + queue +"; + + return $s; +} + +sub datafiles_string +{ + # make a string of the datafiles, comma-separated + my @f; + foreach (1..$binary_grid::grid_options{condor_njobs}) + { + my $jobid=$_.'.'.$binary_grid::grid_options{condor_njobs}; + push(@f,$binary_grid::grid_options{condor_dir}.'/results/'.$jobid); + } + return join(',',@f); +} + +sub condor_rerun_command +{ + # make condor rerun command + my $datafiles=datafiles_string(); + + say "Parent: condor jobs are done, call grid again to join datafiles : $datafiles"; + my $cmd=$binary_grid::grid_options{command_line}; + + # remove unwanted arguments + $cmd=~s/condor_jobid\S+//o; + $cmd=~s/condor_njobs\S+//o; + $cmd=~s/rungrid=1//o; + $cmd=~s/run_flexigrid=1//o; + $cmd=~s/offset=\d+//o; + $cmd=~s/modulo=\d+//o; + $cmd=~s/results_hash_dumpfile=\S+//o; + $cmd=~s/condor_command=\S+//o; + + $cmd=join(' ', + $cmd, + 'rungrid=0', + 'merge_datafiles='.$datafiles, + 'results_hash_dumpfile='.$binary_grid::grid_options{condor_dir}.'/results/joined', + 'condor_command=join_datafiles', + '2>'.$binary_grid::grid_options{condor_dir}.'/stderr/join', + '>'.$binary_grid::grid_options{condor_dir}.'/stdout/join', + 'condor_jobid=join', + ); + + $cmd=~s/\s+/ /go; + return $cmd; +} + +sub condor_check_joining_file +{ + # check for joining file : + if(-f $binary_grid::grid_options{condor_dir}.'/status/joining') + { + # joining file exists : a different process is joining + say "Some other process is joining : I shall not"; + exit; + } + else + { + # touch the "joining" file + open(my $touch,'>'.$binary_grid::grid_options{condor_dir}.'/status/joining'); + print $touch $$; + close $touch; + } + + sleep $binary_grid::grid_options{nfs_sleep}; # wait for the touch to take effect (prevent NFS issues) + + # check joining file has our PID + if(-f $binary_grid::grid_options{condor_dir}.'/status/joining') + { + # joining file exists : is it ours? + open(my $fp,'<'.$binary_grid::grid_options{condor_dir}.'/status/joining'); + my $n=<$fp>; + chomp $n; + if($n!=$$) + { + say "Joining file belongs to a different process : I am doing nothing\n"; + exit; + } + } + else + { + confess("joining file has gone missing!\n"); + } + return condor_rerun_command(); +} + + +sub condor_job_hook +{ + # subroutine called after a condor-flexigrid job has finished + + # find job id + my $jobid=$binary_grid::grid_options{condor_jobid}; + + # set job status to finished + set_condor_job_status($jobid,'finished'); + + # save parent pid and fork child process + my $parent_pid=$$; + my $child_pid=fork(); + + if($child_pid!=0) + { + # parent process + say "condor_check : Parent process $$"; + + # free memory (no longer required, just annoying waste of space) + $binary_grid::grid_options{results_hash}=undef; + + # wait for the child to finish + say "Parent: waitpid for child : $child_pid"; + my $r=waitpid $child_pid,0; + + say "Parent : waitpid returned $r\nParent: check_condor_jobs_done"; + + if(check_condor_jobs_done($binary_grid::grid_options{condor_dir}, + $binary_grid::grid_options{condor_njobs})) + { + # check if the joining file exists, if yes: exits (from condor.pm) + # if no, returns the command to rerun the grid to join all the data + my $cmd=binary_grid::condor::condor_check_joining_file(); + say "Join condor data sets with exec of:\n$cmd"; + `$cmd`; + } + else + { + say 'Parent : condor jobs not finished'; + } + say 'Parent: exit'; + + # exit : do not even return! + exit; + } + else + { + # child process: just return and finish off as per usual + say "condor_check : Child process $$ continuing until output() and end."; + return; + } +} + +############################################################ +############################################################ + +1; + +__END__ + diff --git a/src/perl/modules/rob_misc.pm b/src/perl/modules/rob_misc.pm index 5fa5ddfef..733fa4257 100644 --- a/src/perl/modules/rob_misc.pm +++ b/src/perl/modules/rob_misc.pm @@ -122,6 +122,7 @@ sub trem sub log10 { # log10 function: why doesn't perl have one? + confess("cannot take log of '$_[0]' (zero, nan or inf?)\n")if($_[0]==0.0 || $_[0]=~/(?:nan|inf)/oi); return(log($_[0])/value_of_log10); } @@ -865,7 +866,7 @@ sub get_from_C_header return 'DEFINED'; } }@res; - return ''; + return undef; } sub ssleep @@ -1289,7 +1290,7 @@ sub arithmetically_add_hashes # pass in hash references my $h1=$_[0]; my $h2=$_[1]; - + if(Hash::RobMerge::get_behavior() ne 'SumHashes') { # use the behaviour given in the Hash::Merge documentation @@ -1318,7 +1319,7 @@ sub arithmetically_add_hashes # add the two hashes %$h1 = %{ Hash::RobMerge::merge($h1,$h2) }; - + return $h1; } diff --git a/src/perl/modules/spacing_functions.pm b/src/perl/modules/spacing_functions.pm index edce32899..3fd4593fa 100644 --- a/src/perl/modules/spacing_functions.pm +++ b/src/perl/modules/spacing_functions.pm @@ -72,7 +72,7 @@ sub const_dt my $function=$_[1]; # what to do? # override verbosity ? - $vb=$$opts{vb}//0; + $vb=$$opts{vb}//$vb; print "Const_dt called : opts=$opts, function=$function : const_dt_mass_array length = $#const_dt_mass_array (was $const_dt_mass_array_initial_length)\n"if($vb); @@ -80,11 +80,13 @@ sub const_dt # make required data if it hasn't been made or if the reset call is given if($#const_dt_mass_array==-1 || $function eq 'reset') { + return undef if(defined($const_dt_mass_array_initial_length)); + no warnings; my $optsstring=join '',sort keys %$opts, sort values %$opts; use warnings; - binary_grid::flexigrid_grid_resolution_shift()if($function eq 'reset'); + binary_grid::flexigrid_grid_resolution_shift() if($function eq 'reset'); if(defined($const_dt_results{$optsstring})) { @@ -135,8 +137,9 @@ sub const_dt } elsif($function eq 'next') { - # subsequent calls simply remove masses from the mass_array - my $m=shift @const_dt_mass_array; + # subsequent calls simply remove masses from the mass_array, or use the final + # value if no more are available + my $m = $#const_dt_mass_array==0 ? $const_dt_mass_array[0] : shift @const_dt_mass_array; my $dlnm=$const_dt_mass_spacing{$m}; print STDERR "Shift next mass off the list: M=$m > dlnM=$dlnm\n"if($vb>=2); return $dlnm; @@ -155,7 +158,7 @@ sub adaptive_masshash # input: my $opts=$_[0]; # options hash - + my $vb=1; # output: my $results=$_[1]; # results hash (see const_dt() for definition) @@ -186,11 +189,16 @@ sub adaptive_masshash my $t0; # start time my $tmin; # minimum time my @masses; - - if($vb) + + print "Called by ",(caller(2))[3],"\n"; + +# if($vb) { - printf "adaptive_masshash : range %g to %g \n", - $$opts{mmin},$$opts{mmax}; + printf "adaptive_masshash : opts:\n"; + map + { + print " $_ = $$opts{$_}\n"; + }keys %$opts; print "adaptive_masshash : call make_stellar_lifetime_table($opts) \n"; } @@ -242,11 +250,12 @@ sub adaptive_masshash # calculate time and duration of the helium flash $tflash=stellar_lifetime($mflash,$results); - $dtflash = $logtimes ? $dtflash=10.0 : log10(10.0**$tflash+10.0)*0.5; + $dtflash = $logtimes ? log10(10.0**$tflash+10.0)*0.5 : 10.0; # calculate nearest time bin $tflash=$dt*(int($tflash/$dt)+0.5); $tflash=log10($tflash) if($logtimes); + printf "Extra He-flash resolution at t=%g +- %g\n",$tflash,$dtflash; } # include low-mass stars which would not make it @@ -282,52 +291,47 @@ sub adaptive_masshash ################################### # make main part of the mass list # ################################### - my $prevm; - + # we are required to start at a later time # than $t0 because there is contribution to (e.g) hertzsprung gap # stars at late times which is exaggerated because dlnm will be too # large... now we can assume these have ~20% of the MS lifetime (seems # to be ok...) - + $vb=2; my $tt0=$t0/0.8; printf "Add normal stars from time %g to %g\n",($logtimes ? (10.0**$tmin,10.0**$tt0) : ($tmin,$tt0)) if($vb); for(my $t=$tt0;$t>$tmin;$t-=$dt) { - if(($$opts{extra_flash_resolution})&&(abs($t-$dtflash)<100)) - { - print "Skip flash at $t\n"; - next; - } - my $tt = $logtimes ? 10.0**$t : $t; + my $tt = $logtimes ? 10.0**$t : $t; # the time # from the lifetime, calculate the mass my $m=stellar_mass_from_lifetime($tt,$results); - print "From time $tt > mass = $m\n"if($vb>=2); - $prevm=$m; + print "From time $tt > mass = $m (<mbreak=$mbreak)\n"if($vb>=2); - if($m<$mbreak) - { - push(@masses,$m); - } - else - { - last; - } - $prevm=$m; + last if ($m>=$mbreak); + push(@masses,$m); } # replace bins at $tflash +/ $dt with higher resolution version if($$opts{extra_flash_resolution}) { printf "Increase He flash resolution\n"if($vb); + # increase flash resolution between $tmin and $tmax my $tmin=$tflash-$dtflash+0.5*$dt; my $tmax=$tflash+$dtflash-0.5*$dt; my $hedt=$dt*0.1; $tmin+=0.5*$hedt; $tmax-=0.5*$hedt; + # calculate masses corresponding to $tmin and $tmax + my $mmin=stellar_mass_from_lifetime($logtimes?10.0**$tmin:$tmin,$results); + my $mmax=stellar_mass_from_lifetime($logtimes?10.0**$tmin:$tmin,$results); + + # and remove masses in that region from the masses array + @masses=grep {$_<$mmin || $_>$mmax} @masses; + + # now add extra flash resolution for(my $t=$tmin;$t<=$tmax;$t+=$hedt) { push(@masses,stellar_mass_from_lifetime($logtimes?10.0**$t:$t,$results)); @@ -439,7 +443,7 @@ sub stellar_lifetime { my $m=$_[0]; my $results=$_[1]; - print "Stellar lifetime from $m, $results\n"; + #print "Stellar lifetime from $m, $results\n"; return 10.0**(generic_interpolation_wrapper(log10($m),$$results{evoltimes_list},$$results{evoltimes_hash})); } @@ -467,13 +471,19 @@ sub make_stellar_lifetime_table printf "Making stellar lifetime table (nm=$$opts{stellar_lifetime_table_nm}) from M=%g Msun to %g Msun \n", exp($lnmmin),exp($lnmmax); - # back up hashes + # back up hashes (restored on exit) my %backup_grid_options=%binary_grid::grid_options; my %backup_bse_options=%binary_grid::bse_options; # force single stars, long lived $binary_grid::grid_options{binary}=0; + # force modulo 1 and offset 0 so we run all the stars + $binary_grid::grid_options{modulo}=1; + $binary_grid::grid_options{offset}=0; + + # force single thread only + $$opts{nthreads}//=1; # force evolution time and timestep $binary_grid::grid_options{max_evolution_time}=2.0*$$opts{max_evolution_time}; @@ -514,6 +524,7 @@ sub make_stellar_lifetime_table my $dd=1+int($nm/$$opts{nthreads}); my $low=0; my $nt=0; + print "Running models to make mass grid : please wait...\n"; while($low<$nm) { my $high=$low+$dd; @@ -555,7 +566,7 @@ sub make_stellar_lifetime_table sleep($$opts{thread_sleep}) ; } } - + printf "Done models for mass grid, currently %d threads running\n",scalar(threads->list(threads::running)); print "Finished threads\n"if($vb); @{$$results{evoltimes_list}}=nsort (keys %{$$results{evoltimes_hash}}); @@ -625,6 +636,7 @@ sub lifetime_thread foreach my $lnm (@$mass_splice) { my $args=binary_grid::make_binary_c_argstring(exp($lnm),@alwaysargs); + #printf "RUN M=%g\n",exp($lnm) if($binary_grid::grid_options{vb}); binary_grid::tbse($args,$h); } diff --git a/src/perl/modules/sync b/src/perl/modules/sync index 02fc47f18..000d3f100 100755 --- a/src/perl/modules/sync +++ b/src/perl/modules/sync @@ -7,7 +7,8 @@ my $root=$ENV{HOME}.'/progs/perl/modules'; foreach my $mod ('misc/rob_misc/lib/rob_misc.pm', 'astro/cosmology/lib/cosmology.pm', 'astro/binary_grid/lib/binary_grid.pm', - 'astro/binary_grid-deprecated/lib/binary_grid-deprecated.pm', + 'astro/binary_grid/lib/binary_grid/condor.pm', + #'astro/binary_grid-deprecated/lib/binary_grid-deprecated.pm', 'astro/IMF/lib/IMF.pm', 'astro/distribution_functions/lib/distribution_functions.pm', 'astro/spacing_functions/lib/spacing_functions.pm', @@ -17,10 +18,17 @@ foreach my $mod ('misc/rob_misc/lib/rob_misc.pm', 'maths/Interpolation/lib/Interpolation.pm' ) { - print "Module $mod\n"; - `cp $root/$mod .`; - my $m=$mod;$m=~s/.*\///; - `svn add $m 2>\&1`; + + my $path=$mod; + $path=~s/.*lib\///; + + if($path=~m!(.*)/!) + { + `mkdir -p $1;svn add $1 2>/dev/null`; + } + + print "Module $mod (copy to path $path)\n"; + `cp $root/$mod $path; svn add $path 2>\&1`; } # special case: RobMerge diff --git a/src/setup/parse_arguments.c b/src/setup/parse_arguments.c index b9633ea4a..5fcfe0fc9 100644 --- a/src/setup/parse_arguments.c +++ b/src/setup/parse_arguments.c @@ -26,7 +26,7 @@ #define POSSIBLE_ARGUMENTS "--version, --M_1 <m1>, --M_2 <m2>, --max_evolution_time <time>, --orbital_period <period>, --stellar_type_1 <0-15>, --stellar_type_2 <0-15>, --metallicity <1e-4-0.03>, --eccentricity <0-1>, --eddfac <number>, {these are obselete --monte_carlo_kicks, --hbbfax a b c d e, --testagb, --m1limits <m1min m1max m1n>, --m2limits <m2min m2max m2n>, --seplimits <sepmin sepmax sepn> --snfile <sn_file>, --polyrfac <number>, --alpha_ce <number>, --lambda_ce <number>, --wr_wind_fac <number>, --BH_prescription {0,1}, --snfile <log filename>, --noNeNaMgAl, --dtfac <number>, --nothirdddup, --epsnov <number>, --rlof_gamma <number>, --sn_sigma <number>, --gb_reimers_eta <number>, --bb <number>, --acc2 <number>, --CE_duration <number>, --overshooting_tpagb <number> (if AXEL is defined), --hbbtfac <number>, --initial_abunds_only, --c13_eff <number>, --gbwindfac <number>, --lambda_min <number>, --delta_mcmin <number>, --gbwind {0,1}, --tpagbwind {0,1,2,3}, --superwind_mira_switchon <number>, --tpagb_reimers_eta <number>, --stpagbwind {0,1,2,3}, --stpagb_lambda <number>, --nostpagb, --stpagb_superwind_mira_switchon <number>, --stpagb_reimers_eta <number>, --wr_wind {0,1,2}, --nucreacmult <reaction number> <multiplier>, --hbbfax <5 numbers>, -hbb3fax <3 numbers>, -bbtempmultiplier <number>, --init_abund <isotope number> <initial abundance> --comenv_prescription <0,1> --hachisu_disk_wind <factor>. Doubtless there are others!" - +#define BATCHFLUSH /* */ void parse_arguments(int start, int argc, @@ -47,31 +47,16 @@ void parse_arguments(int start, /** Take a pointer to argv[i] so we don't have to dereference **/ /** the array each time - also to the next arg**/ current_argument=argv[i]; - if(i+1<argc) - { - next_argument=argv[i+1]; - } - else - { - next_argument=argv[i]; - } + next_argument = (i+1<argc) ? argv[i+1] : argv[i]; - /* - Dprint("Argument %d is \"%s\", next argument is %s\n", - i,current_argument,next_argument); - fflush(stdout); - */ + Dprint("Argument %d is \"%s\", next argument is %s\n", + i,current_argument,next_argument); + #ifdef BATCHMODE /* Batch mode arguments */ - /* why are the continue statements in here? */ if(preferences->batchmode>0) { - if(preferences->batchmode==2) - { - printf("Parse argument %s\n",current_argument); - fflush(stdout); - } ARGIS("echo") { preferences->batchmode=2; @@ -79,118 +64,118 @@ void parse_arguments(int start, continue; } else ARGIS("go") - { - int iterate_return_value = iterate(stardata); - printf("iterate returned %d\n",iterate_return_value); - if(preferences->batch_submode==1) - { - printf("fin\n"); - } - continue; - } + { + int iterate_return_value = iterate(stardata); + printf("iterate returned %d\n",iterate_return_value); + if(preferences->batch_submode==1) + { + printf("fin\n"); + } + continue; + } else ARGIS("noecho") - { - printf("set echo off\n");fflush(stdout); - preferences->batchmode=1; - continue; - } + { + printf("set echo off\n");fflush(stdout); + preferences->batchmode=1; + continue; + } else ARGIS("fin") - { - printf("fin\n"); - continue; - } + { + printf("fin\n"); + continue; + } else ARGIS("flush") - { - fflush(stdout); - printf("fin\n"); - fflush(stdout); - continue; - } + { + BATCHFLUSH; + printf("fin\n"); + fflush(stdout); + continue; + } else ARGIS("noechonow") - { - preferences->batchmode=1; - continue; - } + { + preferences->batchmode=1; + continue; + } else ARGIS("bye") - { - if(preferences->batchmode==2) printf("Bye then!\n"); - exit_binaryc(stardata,0); - } + { + if(preferences->batchmode==2) printf("Bye then!\n"); + exit_binaryc(stardata,0); + } else ARGIS("reset_prefs") - { - // reset preferences - if(preferences->batchmode==2) printf("Reset preferences\n"); - reset_preferences(stardata); - continue; - } + { + // reset preferences + if(preferences->batchmode==2) printf("Reset preferences\n"); + reset_preferences(stardata); + continue; + } else ARGIS("reset_stars") - { - // reset preferences - if(preferences->batchmode==2) - { - printf("Reset the stars\n"); - fflush(stdout); - } - init_star(&(stardata->star[1]), - stardata->common.metallicity); - init_star(&(stardata->star[2]), - stardata->common.metallicity); - init_model(&(stardata->model)); - init_common(stardata); - stardata->star[1].starnum=1; - stardata->star[2].starnum=2; - continue; - } + { + // reset preferences + if(preferences->batchmode==2) + { + printf("Reset the stars\n"); + BATCHFLUSH; + } + init_star(&(stardata->star[1]), + stardata->common.metallicity); + init_star(&(stardata->star[2]), + stardata->common.metallicity); + init_model(&(stardata->model)); + init_common(stardata); + stardata->star[1].starnum=1; + stardata->star[2].starnum=2; + continue; + } else ARGIS("defaults") - { - // set just enough defaults to make the code run - if(preferences->batchmode==2) - { - printf("Applying default stellar parameters\n"); - fflush(stdout); - } - stardata->star[1].mass=6.0; - stardata->star[1].zams_mass=0.0; - stardata->star[1].pms_mass=0.0; - stardata->star[1].initial_mass=0.0; - stardata->star[2].mass=3.0; - stardata->star[2].zams_mass=0.0; - stardata->star[2].pms_mass=0.0; - stardata->star[2].initial_mass=0.0; + { + // set just enough defaults to make the code run + if(preferences->batchmode==2) + { + printf("Applying default stellar parameters\n"); + BATCHFLUSH; + } + stardata->star[1].mass=6.0; + stardata->star[1].zams_mass=0.0; + stardata->star[1].pms_mass=0.0; + stardata->star[1].initial_mass=0.0; + stardata->star[2].mass=3.0; + stardata->star[2].zams_mass=0.0; + stardata->star[2].pms_mass=0.0; + stardata->star[2].initial_mass=0.0; #ifdef NUCSYN - stardata->star[1].hezamsmetallicity=0.0; - stardata->star[2].hezamsmetallicity=0.0; - stardata->star[1].num_thermal_pulses=-1.0; - stardata->star[2].num_thermal_pulses=-1.0; - stardata->star[1].temp_mult=1.0; - stardata->star[2].temp_mult=1.0; - + stardata->star[1].hezamsmetallicity=0.0; + stardata->star[2].hezamsmetallicity=0.0; + stardata->star[1].num_thermal_pulses=-1.0; + stardata->star[2].num_thermal_pulses=-1.0; + stardata->star[1].temp_mult=1.0; + stardata->star[2].temp_mult=1.0; + #endif //NUCSYN - stardata->star[1].radius=0.0; - stardata->star[2].radius=0.0; - stardata->common.orbital_period=1e3; - stardata->common.separation=calculate_separation(stardata); - stardata->star[1].stellar_type=1; - stardata->star[2].stellar_type=1; - stardata->common.eccentricity=0; - stardata->common.metallicity=0.02; - stardata->model.max_evolution_time=16000; - preferences->monte_carlo_sn_kicks=TRUE; - continue; - } + stardata->star[1].radius=0.0; + stardata->star[2].radius=0.0; + stardata->common.orbital_period=1e3; + stardata->common.separation=calculate_separation(stardata); + stardata->star[1].stellar_type=1; + stardata->star[2].stellar_type=1; + stardata->common.eccentricity=0; + stardata->common.metallicity=0.02; + stardata->model.max_evolution_time=16000; + preferences->monte_carlo_sn_kicks=TRUE; + continue; + } else ARGIS("status") - { + { #ifdef STARDATA_STATUS - stardata_status(stardata,0); - stardata_status(stardata,1); - stardata_status(stardata,2); + stardata_status(stardata,0); + stardata_status(stardata,1); + stardata_status(stardata,2); #endif - continue; - } - fflush(stdout); + continue; + } + BATCHFLUSH; } - #endif // BATCHMODE + if((argc<2)&&(preferences->batchmode==0)) { printf("You must provide an argument! Try one or more of %s\n",POSSIBLE_ARGUMENTS); @@ -323,7 +308,11 @@ void parse_arguments(int start, } BREAK_OR_ERROR; case 'n': - ARGIS("--nucsyn_metallicity") + ARGIS("--no_arg_checking") + { + stardata->preferences->no_arg_checking=1; + } + EARGIS("--nucsyn_metallicity") { /** set the GCE metallicity **/ #ifdef NUCSYN @@ -1056,6 +1045,10 @@ void parse_arguments(int start, stardata->model.agbwindfac=atof(next_argument); ++i; } + EARGIS("--arg_checking") + { + stardata->preferences->no_arg_checking=0; + } BREAK_OR_ERROR; case 'h': @@ -1328,8 +1321,9 @@ void parse_arguments(int start, fflush(stdout); } /* switch */ - //fprintf(stderr,"CHeck %d %d\n",preferences->batchmode,preferences->batch_submode); - if(preferences->batchmode==1 && preferences->batch_submode==0) + if(preferences->no_arg_checking==0 && + preferences->batch_submode==0 && + preferences->batchmode==1 ) { printf("arg %s %s\n", error==FALSE ? "ok" : "fail", diff --git a/src/setup/set_up_variables.c b/src/setup/set_up_variables.c index 65fc79281..414fc75f8 100644 --- a/src/setup/set_up_variables.c +++ b/src/setup/set_up_variables.c @@ -29,7 +29,7 @@ void set_up_variables( ) { Dprint("set_up_variables preferences=%p stardata->preferences=%p\n", - preferences,stardata->preferences);fflush(stdout); + preferences,stardata->preferences); if(preferences==NULL) { @@ -66,7 +66,7 @@ void set_up_variables( /* set preferences to defaults */ reset_preferences(stardata); - Dprint("done\nargc: %d\n",argc);fflush(stdout); + Dprint("done\nargc: %d\n",argc); if(argc<=1) { diff --git a/src/setup/version.c b/src/setup/version.c index 1ab62c6ff..8ff55c3d4 100644 --- a/src/setup/version.c +++ b/src/setup/version.c @@ -25,6 +25,17 @@ void version(struct stardata_t * RESTRICT stardata) #ifdef AXEL struct preferences_t *preferences=stardata->preferences; #endif +#ifdef SVN_REVISION + printf("SVN revision %d\n",SVN_REVISION); +#else + printf("SVN revision unknown\n"); +#endif +#ifdef SVN_URL + printf("SVN URL %s\n",STRINGOF(SVN_URL)); +#else + printf("SVN URL unknown\n"); +#endif + #ifdef __VERSION__ printf("Compiler version %s\n",__VERSION__); #ifdef __OPTIMIZE__ diff --git a/src/single_star_functions/calc_stellar_wind_mass_loss.c b/src/single_star_functions/calc_stellar_wind_mass_loss.c index 04fbb760c..dc8f9f8fd 100644 --- a/src/single_star_functions/calc_stellar_wind_mass_loss.c +++ b/src/single_star_functions/calc_stellar_wind_mass_loss.c @@ -192,7 +192,7 @@ double calc_stellar_wind_mass_loss(const Stellar_type stellar_type, double rvin,rmva,rmvb,rmjj,rmwr,rmjj2,COHe ; const double eldridge_power=0.5; - Dprint("Getting WR type - stellar type =%d\n",star->stellar_type);fflush(stdout); + Dprint("Getting WR type - stellar type =%d\n",star->stellar_type); wr_type = nucsyn_WR_type(star); Dprint("WR type=%d\n",wr_type); diff --git a/src/wind/wind_AGB.c b/src/wind/wind_AGB.c index 6df938c5a..81ca5aab4 100644 --- a/src/wind/wind_AGB.c +++ b/src/wind/wind_AGB.c @@ -222,13 +222,14 @@ static double wind_TPAGB(WIND_PROTOTYPE_ARGS) star->Xsurf[XC13]/13.0)/ (star->Xsurf[XN14]/14.0+ star->Xsurf[XN14]/15.0); - + printf("M=%g Mc=%g C/O=%g\n",star->mass,star->core_mass,r); + /* * We need more C than O, otherwise CO forms * and more C than N, otherwise CN forms * ASSUME THESE HAVE LOW OPACITY compared to C-C CCC etc. */ - if((r>0.99)&&(rn>0.99)&&(star->num_thermal_pulses>10)) + if((r>0.99))//&&(rn>0.99)&&(star->num_thermal_pulses>10)) { // nC/nO > 0.99 mdot_TPAGB=1e-4; diff --git a/src/wind/wind_WR.c b/src/wind/wind_WR.c index ae0d16735..135b0a623 100644 --- a/src/wind/wind_WR.c +++ b/src/wind/wind_WR.c @@ -51,7 +51,7 @@ double wind_WR(WIND_PROTOTYPE_ARGS) */ (stellar_type!=6)&&(star->zams_mass>8)) { - Dprint("Getting WR type - stellar type =%d\n",stellar_type);fflush(stdout); + Dprint("Getting WR type - stellar type =%d\n",stellar_type); int wr_type = nucsyn_WR_type(star); Dprint("WR type=%d\n",wr_type); diff --git a/src/zfuncs/menvf.c b/src/zfuncs/menvf.c index 1889c7881..e846e6607 100644 --- a/src/zfuncs/menvf.c +++ b/src/zfuncs/menvf.c @@ -46,7 +46,7 @@ CONST_FUNC double menvf(const int kw, result = m - mc; } result = MAX(result,1.0E-10); - Dprint("%f\n",result);fflush(stdout); + Dprint("%f\n",result); return(result); } diff --git a/tbse b/tbse index df621ea70..88933230e 100755 --- a/tbse +++ b/tbse @@ -7,7 +7,7 @@ # the terminal. # perhaps be nice? -NICELEVEL=19 +NICELEVEL=0 # make sure the path has the current directory in it export PATH=$PATH:. @@ -28,16 +28,20 @@ fi #REPEAT=500; RANDOM_SYSTEMS=" --random_systems " REPEAT=1 +# set LOG_FILENAME to /dev/null for no logging +#LOG_FILENAME="/tmp/c_log2.dat" +LOG_FILENAME="/dev/null"; + # stellar types (1=Main sequence default) ST1=1 ST2=1 # masses -M1=1.23 -M2=0.5 +M1=1.0 +M2=0.1 # orbital period and eccentricity -PER=10 +PER=100000000000000 ECC=0 # maximum evolution time @@ -46,7 +50,7 @@ EVOLTIME=13700 YIELDS_DT=1 # metallicity for stellar evolution, must be in the range 1e-4<=Z<=0.03 -Z=0.008 +Z=0.0001 JORB_LOSS=0 # orb ang mom loss prescription 0=Tout, 1=lw, 2=lw hybrid for fast wind WIND_DJORB_FAC=1.0 # multiplies Tout's Jorbdot @@ -121,9 +125,9 @@ CEMP_CFE_MINIMUM=1.0 CEMP_LOGG_MAXIMUM=4.0 # dredge up fudge for low-Z stars #DELTA_MCMIN=-0.1 -DELTA_MCMIN=0 +DELTA_MCMIN=0.0 #LAMBDA_MIN=0.5 -LAMBDA_MIN=0 +LAMBDA_MIN=0.0 # Z=0.008 -> 0.5 # Z=0.004 -> 0.65 # Z=0.0001 -> 0.8 @@ -135,7 +139,7 @@ LAMBDA_MIN=0 # minimum age for cemps (12000 Myr) CEMP_MINIMUM_AGE=10000 # minimum envelope mass for third DUP (0.5=Straniero, 0.0=max DUP) -MINIMUM_ENV_FOR_3DUP=0.5 +MINIMUM_ENV_FOR_3DUP=0.01 # Shift in VW93 relation # +ve shortens the tpagb lifetime, -ve reduces it (0.0) VW93_SHIFT=0.0 @@ -319,7 +323,7 @@ INIT_ABUND_MIX=1 ############################################################ # Now set up the argument string -ARGS=" --stellar_type_1 $ST1 --stellar_type_2 $ST2 --M_1 $M1 --M_2 $M2 --eccentricity $ECC --orbital_period $PER --max_evolution_time $EVOLTIME --metallicity $Z --monte_carlo_kicks --wr_wind $WRWIND --wr_wind_fac $WRWINDFAC --BH_prescription $BHP --c13_eff $C13_EFF --delta_mcmin $DELTA_MCMIN --lambda_min $LAMBDA_MIN --dtfac $DTFAC --gbwind $GBWIND --tpagbwind $TPAGBWIND --gb_reimers_eta $GB_ETA --probability $PROB --eddfac $EDDFAC --acc2 $ACC2 --hbbtfac $HBBTFAC --superwind_mira_switchon $TPAGB_SUPERWIND_MIRA --tpagb_reimers_eta $TPAGB_REIMERS_ETA --sn_sigma $SN_SIGMA --wd_sigma $WD_SIGMA $INIT_ABUNDS --stpagb_superwind_mira_switchon $STPAGB_SUPERWIND_MIRA --stpagb_reimers_eta $STPAGB_REIMERS_ETA --stpagbwind $STPAGBWIND --stpagb_lambda $STPAGB_LAMBDA --stpagbcoremass $STPAGB_COREMASS --CE_duration $CE_DURATION $NUCMULTS --alpha_ce $ALPHA_CE --lambda_ce $LAMBDA_CE $NUCSYN_ROTATION_EFFECTS --present_day 13680 --present_day_period 1 --escape_velocity 1000000000 --escape_fraction 0 --gce_metallicity $GCE_Z --minimum_timestep $MINDT --maximum_timestep $MAXDT --epsnov $EPSNOV --comenv_ns_accretion_fraction $CE_NS_ACC_FRAC --comenv_ns_accretion_mass $CE_NS_ACC_MASS --bb $BB --log_filename /tmp/c_log2.dat --prestpagb_wind $PRESTPAGB_WIND --stpagb_phi $STPAGB_PHI --rlof_gamma $RLOF_GAMMA $NO_PRODUCTION --lambda_multiplier $LAMBDA_MULT $HACH $COMENV $ELD --CEMP_cfe_minimum $CEMP_CFE_MINIMUM $NO_THERMOHALINE --CEMP_logg_maximum $CEMP_LOGG_MAXIMUM --CEMP_minimum_age $CEMP_MINIMUM_AGE --minimum_envelope_mass_for_third_dredgeup $MINIMUM_ENV_FOR_3DUP --vw93_mira_shift $VW93_SHIFT --vw93_multiplier $VW93_MULTIPLIER --settling $SETTLING --comenv_ms_accretion_mass $COMENV_MS_ACCRETION_MASS --yields_dt $YIELDS_DT $HE_IG_EXTRA --nucsyn_metallicity $NUCSYN_Z --rotationally_enhanced_mass_loss $ROTATIONALLY_ENHANCED_MASS_LOSS --rotationally_enhanced_fudge_power $ROTATIONALLY_ENHANCED_FUDGE_POWER --RLOF_timestep_modulator $RLOF_TIMESTEP_MODULATOR --timestep_modulator $TIMESTEP_MODULATOR --vrot_logging_dt $VROT_LOGGING_DT --RLOF_f $RLOF_f --RLOF_mdot_factor $RLOF_MDOT_FACTOR --vrot1 $VROT1 --vrot2 $VROT2 --RLOF_method $RLOF_METHOD --repeat $REPEAT $RANDOM_SYSTEMS --tidal_strength_factor $TIDAL_STRENGTH_FACTOR --wind_mass_loss $WIND_MASS_LOSS --merger_angular_momentum_factor $MERGER_ANGULAR_MOMENTUM_FACTOR $COMENV --mattsson_Orich_tpagbwind $MATTSON_ORICH_TPAGBWIND --initial_abundance_mix $INIT_ABUND_MIX --jorb_loss $JORB_LOSS --lw $LW --wind_djorb_fac $WIND_DJORB_FAC --lambda_ionisation $LAMBDA_ION --wd_kick_direction $WD_KICK_DIRECTION --wd_kick_when $WD_KICK_WHEN --wd_kick_pulse_number $WD_KICK_PULSE_NUMBER --E2_prescription $E2 --starsearch_timestep $STARSEARCH_DT --sn_kick_distribution $SN_KICK_DISTRIBUTION --qcrit_GB_method $QCRIT_GB_METHOD" +ARGS=" --stellar_type_1 $ST1 --stellar_type_2 $ST2 --M_1 $M1 --M_2 $M2 --eccentricity $ECC --orbital_period $PER --max_evolution_time $EVOLTIME --metallicity $Z --monte_carlo_kicks --wr_wind $WRWIND --wr_wind_fac $WRWINDFAC --BH_prescription $BHP --c13_eff $C13_EFF --delta_mcmin $DELTA_MCMIN --lambda_min $LAMBDA_MIN --dtfac $DTFAC --gbwind $GBWIND --tpagbwind $TPAGBWIND --gb_reimers_eta $GB_ETA --probability $PROB --eddfac $EDDFAC --acc2 $ACC2 --hbbtfac $HBBTFAC --superwind_mira_switchon $TPAGB_SUPERWIND_MIRA --tpagb_reimers_eta $TPAGB_REIMERS_ETA --sn_sigma $SN_SIGMA --wd_sigma $WD_SIGMA $INIT_ABUNDS --stpagb_superwind_mira_switchon $STPAGB_SUPERWIND_MIRA --stpagb_reimers_eta $STPAGB_REIMERS_ETA --stpagbwind $STPAGBWIND --stpagb_lambda $STPAGB_LAMBDA --stpagbcoremass $STPAGB_COREMASS --CE_duration $CE_DURATION $NUCMULTS --alpha_ce $ALPHA_CE --lambda_ce $LAMBDA_CE $NUCSYN_ROTATION_EFFECTS --present_day 13680 --present_day_period 1 --escape_velocity 1000000000 --escape_fraction 0 --gce_metallicity $GCE_Z --minimum_timestep $MINDT --maximum_timestep $MAXDT --epsnov $EPSNOV --comenv_ns_accretion_fraction $CE_NS_ACC_FRAC --comenv_ns_accretion_mass $CE_NS_ACC_MASS --bb $BB --log_filename $LOG_FILENAME --prestpagb_wind $PRESTPAGB_WIND --stpagb_phi $STPAGB_PHI --rlof_gamma $RLOF_GAMMA $NO_PRODUCTION --lambda_multiplier $LAMBDA_MULT $HACH $COMENV $ELD --CEMP_cfe_minimum $CEMP_CFE_MINIMUM $NO_THERMOHALINE --CEMP_logg_maximum $CEMP_LOGG_MAXIMUM --CEMP_minimum_age $CEMP_MINIMUM_AGE --minimum_envelope_mass_for_third_dredgeup $MINIMUM_ENV_FOR_3DUP --vw93_mira_shift $VW93_SHIFT --vw93_multiplier $VW93_MULTIPLIER --settling $SETTLING --comenv_ms_accretion_mass $COMENV_MS_ACCRETION_MASS --yields_dt $YIELDS_DT $HE_IG_EXTRA --nucsyn_metallicity $NUCSYN_Z --rotationally_enhanced_mass_loss $ROTATIONALLY_ENHANCED_MASS_LOSS --rotationally_enhanced_fudge_power $ROTATIONALLY_ENHANCED_FUDGE_POWER --RLOF_timestep_modulator $RLOF_TIMESTEP_MODULATOR --timestep_modulator $TIMESTEP_MODULATOR --vrot_logging_dt $VROT_LOGGING_DT --RLOF_f $RLOF_f --RLOF_mdot_factor $RLOF_MDOT_FACTOR --vrot1 $VROT1 --vrot2 $VROT2 --RLOF_method $RLOF_METHOD --repeat $REPEAT $RANDOM_SYSTEMS --tidal_strength_factor $TIDAL_STRENGTH_FACTOR --wind_mass_loss $WIND_MASS_LOSS --merger_angular_momentum_factor $MERGER_ANGULAR_MOMENTUM_FACTOR $COMENV --mattsson_Orich_tpagbwind $MATTSON_ORICH_TPAGBWIND --initial_abundance_mix $INIT_ABUND_MIX --jorb_loss $JORB_LOSS --lw $LW --wind_djorb_fac $WIND_DJORB_FAC --lambda_ionisation $LAMBDA_ION --wd_kick_direction $WD_KICK_DIRECTION --wd_kick_when $WD_KICK_WHEN --wd_kick_pulse_number $WD_KICK_PULSE_NUMBER --E2_prescription $E2 --starsearch_timestep $STARSEARCH_DT --sn_kick_distribution $SN_KICK_DISTRIBUTION --qcrit_GB_method $QCRIT_GB_METHOD" #echo $ARGS;exit ionice -c3 nice -n $NICELEVEL $BIN_C $ARGS; exit -- GitLab