-
Izzard, Robert Dr (Maths & Physics) authoredfix timesteps to correct for the wrong size of the arrays
Izzard, Robert Dr (Maths & Physics) authoredfix timesteps to correct for the wrong size of the arrays
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
evolution.c 8.07 KiB
#include "../binary_c.h"
/*
* Wrapper to evolve the system in time, starting at time t,
* and do:
*
* 1) stellar mass and angular momentum updates
* 2) stellar evolution
* 3) binary evolution
*
* leaving the stars at time t+dt.
*/
int evolution(struct stardata_t * RESTRICT stardata,
Boolean system_type,
Boolean update_time)
{
/*
* Apply stellar mass and angular momentum changes,
* e.g. wind loss, wind accretion, RLOF
*/
double dtm_in = stardata->model.dtm;
double dt_in = stardata->model.dt;
if((stardata->preferences->solver == SOLVER_RK2 ||
stardata->preferences->solver == SOLVER_RK4) &&
stardata->model.solver_step==0)
{
stardata->model.intermediate_step = TRUE;
}
else
{
stardata->model.intermediate_step = FALSE;
}
if(stardata->model.intermediate_step == TRUE &&
stardata->model.solver_step == 0)
{
stardata->model.dtm *= 0.5;
stardata->model.dt *= 0.5;
printf("DT *0.5 -> dtm = %g, dt = %g\n",
stardata->model.dtm,
stardata->model.dt);
}
printf("in evolution() model %d at t = %g, dt = %g, dtm = %g : solver_step %d (intermediate? %d)\n",
stardata->common.model_number,
stardata->model.time,
stardata->model.dt,
stardata->model.dtm,
stardata->model.solver_step,
stardata->model.intermediate_step
);
int stop = mass_angmom_and_evolutionary_changes(stardata,system_type);
int status = 0;
if(stop == 0)
{
/*
* When should we rejuv, when should we not?
* He-star fail suggests we should, but
* other systems suggest we should not.
*/
struct model_t * model = &stardata->model;
if(system_type==DETACHED_SYSTEM_CALL)
{
/*
* Detached system */
save_detached_stellar_types(stardata);
if(model->intpol==0)
{
model->tphys0 = model->time;
model->dtm0 = model->dtm;
}
}
else
{
/*
* RLOFing system.
*
* Rejuvenate the secondary and
* age the primary if they are on
* the main sequence.
*/
#ifdef BSE
rejuvenate_MS_secondary_and_age_primary(stardata);
#endif
}
/*
* Set the time to t+dt
*/
if(update_time == TRUE)
update_the_time(stardata);
Dprint("Advance the time: intpol=%d time was %30.22e (dt=%30.22e dtm=%30.22e) now %30.22e (tphys00 %30.22e, maxt=%g)\n",
model->intpol,
model->time-model->dtm,
model->dt,
model->dtm,
model->time,
model->tphys0,
model->max_evolution_time);
NANCHECK(model->dt);
NANCHECK(model->dtm);
/*
* Do stellar evolution over the timestep
*/
status = stellar_evolution(stardata,system_type);
}
/*
* Do binary stellar evolution over the timestep,
* also checks for supernovae and updates the orbit
* if required.
*/
binary_star_evolution(stardata,&status);
int retval;
if(stop!=0)
{
/*
* Return without updating the time:
* there has probably been a common envelope
* or a merger or something like that.
*/
retval = -stop;
}
else
{
/*
* Update stellar spins (by conserving angular momenta)
*/
calculate_spins(stardata);
if(system_type==DETACHED_SYSTEM_CALL)
{
if(stardata->model.sgl==FALSE)
{
/* in binaries, update the radius derivative */
adjust_radius_derivative(stardata);
}
else
{
/* in single stars, make Roche lobes huge */
make_roche_lobe_radii_huge(stardata);
}
if(stardata->model.supernova==TRUE) stardata->model.dtm = 0.0;
}
update_phase_start_times(stardata);
Boolean exhausted = check_for_time_exhaustion(stardata,stardata->model.intpol);
printf("exhausted? (evolution) %d\n",exhausted);
retval = exhausted == TRUE ? -STOP : status;
}
/*
* We may need to call this function again
* if using a non-Euler scheme
*/
if(retval == KEEP_LOOPING_RESULT)
{
#define Restore \
{ \
struct stardata_t ** stack = stardata->stardata_stack; \
evolution_restore_from_previous(stardata); \
stardata->stardata_stack = stack; \
}
#undef Restore
#define Restore /* */
if(stardata->model.solver_step>0)
{
/* do no calculations if we're an intermediate step */
return retval;
}
if(stardata->preferences->solver == SOLVER_RK2)
{
#warning TODO hand events properly in RK2, RK4
if(stardata->model.solver_step == 0)
{
/*
* RK2 second step and apply derivatives at time t + dt/2
* to obtain stars at t+dt
*/
stardata->model.solver_step = 1;
stardata->model.intermediate_step = FALSE;
int retval2 = evolution(stardata,system_type,UPDATE_TIME);
retval = retval2;
}
}
else if(stardata->preferences->solver == SOLVER_RK4)
{
/*
* we have k1 at t
* calculate k2 at t+dt/2
*/
printf("DT *0.5 -> dtm = %g, dt = %g\n",
stardata->model.dtm,
stardata->model.dt);
double t2 = stardata->model.time;
stardata->model.solver_step = 1;
evolution(stardata,system_type,DO_NOT_UPDATE_TIME);
/* calculate k3 at t+dt/2 (do not update the time!) */
Restore;
stardata->model.solver_step = 2;
stardata->model.time = t2;
evolution(stardata,system_type,UPDATE_TIME);
/* calculate k4 at t+dt */
Restore;
stardata->model.solver_step = 3;
stardata->model.intermediate_step = FALSE;
retval = evolution(stardata,system_type,DO_NOT_UPDATE_TIME);
}
else if(stardata->preferences->solver == SOLVER_PREDICTOR_CORRECTOR)
{
/*
* PEC method
* https://en.wikipedia.org/wiki/Predictor%E2%80%93corrector_method
*/
const int maxk = 10;
int i;
for(i=1;i<maxk;i++)
{
stardata->model.solver_step = i;
printf("PEC%d\n",i);
Restore;
evolution(stardata,system_type,DO_NOT_UPDATE_TIME);
}
}
/* go back to step 0 */
stardata->model.solver_step = 0;
/* restore timestep */
stardata->model.dt = dt_in;
stardata->model.dtm = dtm_in;
printf("DT RESTORE -> dtm = %g, dt = %g\n",
stardata->model.dtm,
stardata->model.dt
);
if(IS_ZERO(stardata->model.dt))
{
/* reset in case of repeat */
stardata->common.test = 0.0;
}
double absdiff = ABSDIFF(stardata->model.time*1e6,
stardata->common.test);
printf("POSTTEST %30.12e %30.12e : diff %g\n",
stardata->model.time*1e6,
stardata->common.test,
absdiff
);
if(absdiff > 1e-6)
{
printf("integration failed\n");
fflush(NULL);
_exit(0);
}
}
printf("evolution retval %d\n",retval);
return retval;
}