-
Izzard, Robert Dr (Maths & Physics) authoredadd pthread warnings
Izzard, Robert Dr (Maths & Physics) authoredadd pthread warnings
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
evolution_split.c 14.67 KiB
#include "../binary_c.h"
No_empty_translation_unit_warning;
/*
* Evolution run splitting algorithm
*
* Required EVOLUTION_SPLITTING to be defined (in binary_c_parameters.h)
* and
* stardata->preferences->evolution_splitting to be non-zero
*
* 1) Detect the end of an evolutionary run, if we have previously
* split, restore the previous evolutionary sequence and go again
* if required.
*
* The max number of splits is set on the command line
* in the preference evolution_splitting_maxdepth, up to a
* hard-wired maximum EVOLUTION_SPLITTING_HARD_MAX_DEPTH.
* Obviously, too much splitting will cause the population
* to grind to a halt.
*
* Returns EVOLUTION_SPLIT_RESTART in this case.
*
* 2) Detect when a new split should occur. This is controlled by
* the splittest function (below) which returns the number of runs
* to be a forked at each split.
*
* Currently, this is implemented for supernovae, but it could
* work at any point in the evolutionary run (e.g. common envelope
* evolution).
*
* Returns EVOLUTION_SPLIT_NEW_SPLIT in this case.
*
* 3) Returns EVOLUTION_SPLIT_CONTINUE otherwise.
*
* When a new split is made, a splitinfo struct is created which
* contains the number of runs remaining, and a copy of the stardata
* at the point of the split. This is then iterated over until the number
* of runs left is zero.
*
* Note that at the end of evolution, evolution_cleanup is called.
* Currently this does NOT free memory or close log files, these are
* allowed to continue, which is probably what you want.
*
*/
#ifdef EVOLUTION_SPLITTING
#define In_splitdepth(STARDATA) \
((STARDATA)->preferences->current_splitdepth < \
Min((STARDATA)->preferences->evolution_splitting_maxdepth, \
EVOLUTION_SPLITTING_HARD_MAX_DEPTH))
// split test function
static int splittest(struct stardata_t * RESTRICT const stardata);
static void modulate_stardata(struct stardata_t * RESTRICT const stardata,
const double f);
// debugging output if required
#define SDEBUG Dprint
//#define SDEBUG(...) if(stardata->model.time > 0) { printf("split : ");printf(__VA_ARGS__); }
#ifdef NUCSYN
#if defined NUCSYN_YIELD_COMPRESSION || defined NUCSYN_LOG_BINARY_DX_YIELDS
#define SAVE_YIELDS
#endif#endif//NUCSYN
#endif // EVOLUTION_SPLITTING
int evolution_split(struct stardata_t * RESTRICT const stardata,
int * RESTRICT const status)
{
int retval = EVOLUTION_SPLIT_CONTINUE;
#ifdef EVOLUTION_SPLITTING
if(stardata->preferences->evolution_splitting)
{
#define Break_evolution_Q \
(abs(*status) == EVOLUTION_BREAK || \
abs(*status) == EVOLUTION_STOP)
SDEBUG("splitting is on at t=%g : status = %d cf %d %d (break? %s): stellar types %d %d : splittest = %d, doingsplit = %d, In_splitdepth = %d (current depth %d vs Min(%d,%d))\n",
stardata->model.time,
*status,
EVOLUTION_BREAK,
EVOLUTION_STOP,
Yesno(Break_evolution_Q),
stardata->star[0].stellar_type,
stardata->star[1].stellar_type,
splittest(stardata),
stardata->model.doingsplit,
In_splitdepth(stardata),
stardata->preferences->current_splitdepth,
stardata->preferences->evolution_splitting_maxdepth,
EVOLUTION_SPLITTING_HARD_MAX_DEPTH
);
if(Break_evolution_Q)
{
/*
* flush yield logs of previous evolution, except
* the very last time which is handled by evolve_system_*
*/
if(!(stardata->preferences->current_splitdepth==1 &&
stardata->preferences->splitinfo[stardata->preferences->current_splitdepth - 1]->count==0))
{
/* do logging and cleanup */
SDEBUG("Split do logging and cleanup\n");
Evolution_logging;
evolution_cleanup(stardata,FALSE);
}
while(stardata->preferences->current_splitdepth > 0)
{
struct splitinfo_t * p = stardata->preferences->splitinfo[stardata->preferences->current_splitdepth - 1];
/*
* Reached end of evolution
*
* Check if we need to rerun from a previous
* splitpoint
*/
SDEBUG("SPLIT status is evolution_break : currently working on split-info struct %d at %p, of which count = %d\n",
stardata->preferences->current_splitdepth-1,
(void*)p,
p->count);
if(p->count > 0)
{
/*
* yes, need to, restore almost everything
* and lower the count
*
* NB save random number seed information: this
* should not be reset otherwise all the post-split stars * will be identical to the pre-split stars
*/
long random_seed = stardata->common.random_seed;
struct binary_c_random_buffer_t random_buffer = stardata->common.random_buffer;
p->count--;
erase_events(stardata);
free_previous_stardatas(stardata);
copy_stardata(p->stardata,
stardata,
COPY_STARDATA_PREVIOUS_COPY,
COPY_STARDATA_PERSISTENT_FROM_POINTER);
stardata->common.random_seed = random_seed;
stardata->common.random_buffer = random_buffer;
*status = EVOLUTION_CONTINUE;
SDEBUG("SPLIT restarting depth=%d, count=%d at t=%g, stellar types %d %d\n",
stardata->preferences->current_splitdepth - 1,
p->count,
stardata->model.time,
stardata->star[0].stellar_type,
stardata->star[1].stellar_type
);
Append_logstring(LOG_SPLIT_EVOLUTION,
"Evolution restart @ t=%g, %d.%d",
stardata->model.time,
stardata->preferences->current_splitdepth-1,
p->count);
#ifdef LINUX
reset_binary_c_timeout();
#endif
retval = EVOLUTION_SPLIT_RESTART;
break;
}
else
{
stardata->preferences->current_splitdepth--;
SDEBUG("SPLIT last done (count=%d) : next depth is %d\n",
p->count,
stardata->preferences->current_splitdepth-1);
/* clean up memory */
struct stardata_t * s = stardata->preferences->splitinfo[stardata->preferences->current_splitdepth]->stardata;
free_previous_stardatas(s);
Safe_free(s);
Safe_free(stardata->preferences->splitinfo[stardata->preferences->current_splitdepth]);
retval = EVOLUTION_SPLIT_CONTINUE;
}
}
stardata->model.split_last_time = FALSE;
}
else if(stardata->model.doingsplit == FALSE &&
In_splitdepth(stardata))
{
/*
* check for a new split case in this timestep
*/
int split = splittest(stardata);
SDEBUG("Splittest : %d\n",split);
if(split>=1)
{
/*
* Found a new split case
*/ SDEBUG("SPLIT AT t = %g : current_split_depth = %d\n",
stardata->model.time,
stardata->preferences->current_splitdepth);
/* add new struct */
struct splitinfo_t * p
= stardata->preferences->splitinfo[stardata->preferences->current_splitdepth]
= Malloc(sizeof(struct splitinfo_t));
SDEBUG("malloced new split-info struct at %p\n",(void*)p);
p->depth = stardata->preferences->current_splitdepth;
/* should be rerun this many times */
p->count = split - 1;
/*
* allocate space
*/
p->stardata = new_stardata(stardata->preferences);
/*
* Copy previous stardata structure to p->stardata
*/
copy_stardata(
stardata->previous_stardata,
p->stardata,
COPY_STARDATA_PREVIOUS_COPY,
COPY_STARDATA_PERSISTENT_FROM_POINTER
);
/*
* Keep preferences
*/
p->stardata->preferences = stardata->preferences;
SDEBUG("starting post-split %d at t=%g\n",
p->count,
p->stardata->model.time);
/*
* modulate probability, yields, ensemble etc. which
* have been saved up until now
*/
modulate_stardata(p->stardata,1.0/(double)split);
SDEBUG("Rerunning from time stardata = %g, p->stardata = %g; restore stellar types to %d %d\n",
stardata->model.time,
p->stardata->model.time,
p->stardata->star[0].stellar_type,
p->stardata->star[1].stellar_type
);
/* set the split_last_time flag to avoid a double split */
p->stardata->model.split_last_time = TRUE;
/* set the doingsplit number */
p->stardata->model.doingsplit = TRUE;
/* next split should be one level up */
p->stardata->preferences->current_splitdepth++;
SDEBUG("SPLIT set current_splitdepth=%d\n",stardata->preferences->current_splitdepth);
/*
* Free stardata's previous_stardatas, events etc.
* so we can replace stardata
*/
free_previous_stardatas(stardata);
erase_events(stardata);
/*
* put back in stardata for restart
*/
copy_stardata(p->stardata,
stardata,
COPY_STARDATA_PREVIOUS_COPY,
COPY_STARDATA_PERSISTENT_FROM_POINTER);
stardata->model.doingsplit = TRUE;
Append_logstring(LOG_SPLIT_EVOLUTION,
"Evolution split and restart at t=%g, %d.%d",
stardata->model.time,
stardata->preferences->current_splitdepth-1,
p->count);
#ifdef LINUX
reset_binary_c_timeout();
#endif
/* return and restart this timestep */
retval = EVOLUTION_SPLIT_NEW_SPLIT;
}
}
}
#endif // EVOLUTION_SPLITTING
/*
* We're not doing a split if we're just
* continuing
*/
if(retval == EVOLUTION_SPLIT_CONTINUE)
{
stardata->model.doingsplit = FALSE;
}
return retval;
}
#ifdef EVOLUTION_SPLITTING
static int splittest(struct stardata_t * RESTRICT const stardata)
{
struct stardata_t * previous_stardata = stardata->previous_stardata;
int split = 0;
Star_number i;
/*
* Detect new split cases after stellar evolution has
* been done in this system. You can use stardata and
* stardata->previous_stardata to do this.
*
* Returns a number saying with which resolution the split
* should act. e.g. returning 10 means the evolution is split into
* 10 new sub-runs.
*
* Returns zero if there should be no split.
*/
Starloop(i)
{
if(!split)
{
// Supernova
SDEBUG("check SN %d : type %d\n",
stardata->model.supernova,
stardata->star[i].stellar_type);
if(stardata->model.supernova || ((stardata->star[i].stellar_type == BLACK_HOLE ||
stardata->star[i].stellar_type == NEUTRON_STAR) &&
(previous_stardata->star[i].stellar_type<=HeGB)))
{
split = stardata->preferences->evolution_splitting_sn_n;
}
}
}
return split;
}
static void modulate_stardata(struct stardata_t * RESTRICT const stardata,
const double f)
{
/*
* Function to "modulate" a stardata, i.e. multiply some integrals
* by a factor f to represent an evolution split into 1/f runs.
*
* Not all things will want to be modulated, e.g. stellar structure
* variables. However, some - like integrated yields which are output
* on a timestep-by-timestep basis, or at the end of evolution
* (which is really just the sum of the timestep-by-timestep basis data)
* will need to be modulated, as will the probability. The reason
* is that the subsequent evolution will be repeated 1/f times,
* and we don't want to overcount.
*/
/* probability should be reduced for rerun */
stardata->model.probability *= f;
/* various arrays also need modulating */
#ifdef NUCSYN
{
Star_number k;
#ifdef NUCSYN_ID_SOURCES
Source_number j;
#endif
#ifdef STELLAR_POPULATIONS_ENSEMBLE
Ensemble_index i;
#endif
Starloop(k)
{
SETstar(k);
Xmult(star->Xyield,f);
#ifdef LOG_BINARY_MPYIELDS
Xmult(star->mpyield,f);
#endif
#ifdef SAVE_YIELDS
if(k==0) Xmult(stardata->model.oldyields,f);
#ifdef NUCSYN_ID_SOURCES
Sourceloop(j)
{
Xmult(&(stardata->model.oldyieldsources[k][j][0]),f);
}
#endif //NUCSYN_ID_SOURCES
#endif // NUCSYN_YIELD_COMPRESSION || NUCSYN_LOG_BINARY_DX_YIELDS
#ifdef NUCSYN_ID_SOURCES
Sourceloop(j)
{
Xmult(&(star->Xsource[j][0]),f);
}
#endif
#ifdef STELLAR_POPULATIONS_ENSEMBLE
for(i=0;i<ENSEMBLE_NUMBER;i++)
{
stardata->model.ensemble_store[k][i] *= f;
}
#endif
}
}
#endif //NUCSYN
}
#endif // EVOLUTION_SPLITTING