Solver Class Reference

#include <solver.h>

Collaboration diagram for Solver:

Public Member Functions
	Solver (PSPM_SolverType _method, std::string ode_method="rk45ck")
void	addSystemVariables (int _s)
void	addSpecies (int _J, double _xb, double _xm, bool log_breaks, Species_Base *_mod, int n_extra_vars, double input_birth_flux=-1)
void	addSpecies (std::vector< double > xbreaks, Species_Base *_mod, int n_extra_vars, double input_birth_flux=-1)
int	n_species ()
void	setEnvironment (EnvironmentBase *_env)
void	resetState (double t0=0)
void	resizeStateFromSpecies ()
void	initialize ()
void	copyStateToCohorts (std::vector< double >::iterator state_begin)
void	copyCohortsToState ()
double	maxSize ()
void	updateEnv (double t, std::vector< double >::iterator S, std::vector< double >::iterator dSdt)
void	calcRates_FMU (double t, std::vector< double >::iterator S, std::vector< double >::iterator dSdt)
void	calcRates_iFMU (double t, std::vector< double >::iterator S, std::vector< double >::iterator dSdt)
void	stepU_iFMU (double t, std::vector< double > &S, std::vector< double > &dSdt, double dt)
void	calcRates_EBT (double t, std::vector< double >::iterator S, std::vector< double >::iterator dSdt)
void	addCohort_EBT ()
void	removeDeadCohorts_EBT ()
void	calcRates_CM (double t, std::vector< double >::iterator S, std::vector< double >::iterator dSdt)
void	addCohort_CM ()
void	removeCohort_CM ()
template<typename AfterStepFunc >
void	step_to (double tstop, AfterStepFunc &afterStep_user)
void	step_to (double tstop)
double	calcSpeciesBirthFlux (int k, double t)
std::vector< double >	newborns_out (double t)
std::vector< double >	u0_out (double t)
void	print ()
template<typename wFunc >
double	integrate_x (wFunc w, double t, int species_id)
	Computes the integral \[I = \int_{x_b}^{x_m} w(z,t)u(z)dz\] for the specified species. For details, see integrate_wudx_above. More...
template<typename wFunc >
double	integrate_wudx_above (wFunc w, double t, double xlow, int species_id)
	Computes the partial integral \[I = \int_{x_{low}}^{x_m} w(z,t)u(z)dz\] for the specified species. More...
std::vector< double >	getDensitySpecies_EBT (int k, std::vector< double > breaks)

Public Attributes
OdeSolver	odeStepper
EnvironmentBase *	env
double	current_time
std::vector< double >	state
std::vector< double >	rates
std::vector< Species_Base * >	species_vec
struct {
double   ode_eps = 1e-6
double   ode_initial_step_size = 1e-6
double   convergence_eps = 1e-6
double   cm_grad_dx = 1e-6
bool   update_cohorts = true
int   max_cohorts = 500
double   ebt_ucut = 1e-10
double   ebt_grad_dx = 1e-6
double   ode_rk4_stepsize = 0.1
double   ode_ifmu_stepsize = 0.1
bool   ifmu_centered_grids = true
bool   integral_interpolate = true
}	control
bool	use_log_densities = true

Private Attributes
PSPM_SolverType	method
int	n_statevars_internal = 0
int	n_statevars_system = 0

Detailed Description

Definition at line 15 of file solver.h.

Constructor & Destructor Documentation

Solver()

Solver::Solver	(	PSPM_SolverType	_method,
		std::string	ode_method = "rk45ck"
	)

Definition at line 37 of file solver.cpp.

                                                         : odeStepper(ode_method, 0, 1e-6, 1e-6) {
    method = _method;
}

Member Function Documentation

addCohort_CM()

void Solver::addCohort_CM

(

)

Definition at line 58 of file cm.cpp.

                         {
    cout << ".......... add cohorts ...............\n";
    updateEnv(current_time, state.begin(), rates.begin());
    
    for (int s = 0; s< species_vec.size(); ++s){
        auto spp = species_vec[s];
        
        // calculate u for boundary cohort
        double gb = spp->growthRate(-1, spp->xb, current_time, env);
        double pe = spp->establishmentProbability(current_time, env);
        if (spp->birth_flux_in < 0){    
            double birthFlux = calcSpeciesBirthFlux(s,current_time) * pe;
            spp->set_ub(birthFlux/(gb+1e-20));
        }
        else{
            spp->calc_boundary_u(gb, pe);  // init density of boundary cohort
        }

        spp->initBoundaryCohort(current_time, env); // init extra state variables and birth time of the boundary cohort
        spp->addCohort();   // introduce copy of boundary cohort in system
    }
    
    resizeStateFromSpecies();
    copyCohortsToState();
}

Here is the caller graph for this function:

addCohort_EBT()

void Solver::addCohort_EBT

(

)

Definition at line 70 of file ebt.cpp.

                          {
    
    for (auto spp : species_vec){
        // 1. internalize the pi0-cohort (this cohort's birth time would have been already set when it was inserted)
        // - get pi0, N0 from last cohort
        double   pi0  =  spp->getX(spp->J-1);
        double   N0   =  spp->getU(spp->J-1);

        // - update the recently internalized pi0-cohort with actual x0 value
        double x0 = spp->xb + pi0/(N0+1e-12);
        spp->setX(spp->J-1, x0);
        spp->setU(spp->J-1, N0);

        // 2. insert a new cohort (copy of boundary cohort, to be the new pi0-cohort)
        spp->initBoundaryCohort(current_time, env); // update initial extra state and birth-time of boundary cohort
        spp->addCohort(); // introduce copy of boundary cohort into species

        // 3. set x,u of the new cohort to 0,0, thus marking it as the new pi0-cohort
        spp->setX(spp->J-1, 0); 
        spp->setU(spp->J-1, 0);
    }

    // 4. reset state from cohorts
    resizeStateFromSpecies(); 
    copyCohortsToState();
    
}

Here is the caller graph for this function:

addSpecies() [1/2]

void Solver::addSpecies	(	int	_J,
		double	_xb,
		double	_xm,
		bool	log_breaks,
		Species_Base *	_mod,
		int	n_extra_vars,
		double	input_birth_flux = -1
	)

Definition at line 106 of file solver.cpp.

                                                                         {
    vector<double> xbreaks;
    if (log_breaks) xbreaks = logseq(_xb, _xm, _J+1);
    else            xbreaks =    seq(_xb, _xm, _J+1);

    addSpecies(xbreaks, s, n_extra_vars, input_birth_flux);
}

Here is the call graph for this function:

addSpecies() [2/2]

void Solver::addSpecies	(	std::vector< double >	xbreaks,
		Species_Base *	_mod,
		int	n_extra_vars,
		double	input_birth_flux = -1
	)

Definition at line 42 of file solver.cpp.

                                                                                                            {
    s->set_inputBirthFlux(input_birth_flux);
    s->n_extra_statevars = n_extra_vars;

    if (method == SOLVER_FMU || method == SOLVER_IFMU){
        std::sort(xbreaks.begin(), xbreaks.end(), std::less<double>());  // sort cohorts ascending for FMU
        s->xb = *xbreaks.begin();
    }
    else {
        std::sort(xbreaks.begin(), xbreaks.end(), std::greater<double>());  // sort cohorts descending for CM, EBT
        s->xb = *xbreaks.rbegin();
    }

    int J;
    if      (method == SOLVER_FMU)  J = xbreaks.size()-1;   
    else if (method == SOLVER_IFMU) J = xbreaks.size()-1;   
    else if (method == SOLVER_MMU)  J = xbreaks.size()-1;  
    else if (method == SOLVER_CM )  J = xbreaks.size();
    else if (method == SOLVER_EBT)  J = xbreaks.size();
    else    throw std::runtime_error("Unsupported method");

    s->x = xbreaks;
    s->resize(J);

    // FMU has only 1 internal state variable (x), reset have 2 (x,u)
    bool cond = (method == SOLVER_FMU || method == SOLVER_IFMU);
    n_statevars_internal = (cond)? 1:2;
    
    int state_size = s->J*(n_statevars_internal + n_extra_vars);    // x/u and extra vars

    //s->start_index = state.size();    // New species will be appended to the end of state vector
    state.resize(state.size()+state_size);  // This will resize state for all species additions, but this in only initialization so its okay.
    rates.resize(rates.size()+state_size);

    species_vec.push_back(s);

    if (method == SOLVER_FMU){  
        s->X.resize(s->J);
        for (size_t i=0; i<s->J; ++i) s->X[i] = (s->x[i]+s->x[i+1])/2.0;
        
        s->h.resize(s->J);  // This will be used only by FMU
        for (size_t i=0; i<s->J; ++i) s->h[i] = s->x[i+1] - s->x[i];    
    }

    if (method == SOLVER_IFMU){
        s->X.resize(s->J);
        s->h.resize(s->J);  // This will be used only by FMU
        
        if (control.ifmu_centered_grids){ // grids are labelled by the grid center  
            for (size_t i=0; i<s->J; ++i) s->X[i] = (s->x[i]+s->x[i+1])/2.0;
            for (size_t i=0; i<s->J; ++i) s->h[i] = s->x[i+1] - s->x[i];    
            // for centered grids, also set xb to X[0] rather than x[0]. 
            s->xb = s->X[0];
        }
        else{  // grids are labelled by the lower edge (this ensures that recruitment happens exactly at xb)
            for (size_t i=0; i<s->J; ++i) s->X[i] = s->x[i];
            for (size_t i=0; i<s->J; ++i) s->h[i] = s->x[i+1] - s->x[i];    
            s->xb = s->X[0];
        }
    }

}

Here is the call graph for this function:

addSystemVariables()

void Solver::addSystemVariables

(

int

_s

)

Definition at line 116 of file solver.cpp.

                                     {
    n_statevars_system = _s;
    state.resize(state.size() + _s);
    rates.resize(state.size() + _s);
}

calcRates_CM()

void Solver::calcRates_CM	(	double	t,
		std::vector< double >::iterator	S,
		std::vector< double >::iterator	dSdt
	)

Definition at line 9 of file cm.cpp.

                                                                                        {

    //auto ss = species_vec[0]; 
    //cout << "svec: " << t << " " << S[0] << " " << S[1] << " " << S[2] << " " << S[3] << "\n";
    //cout << "coho: " << t << " " << ss->getX(0) << " " << ss->getU(0) << " " << ss->getX(1) << " " << ss->getU(1) << "\n";
    
    vector<double>::iterator its = S    + n_statevars_system; // Skip system variables
    vector<double>::iterator itr = dSdt + n_statevars_system;

    for (int s = 0; s<species_vec.size(); ++s){
        Species_Base* spp = species_vec[s];
        cout << "calcRates(species = " << s << ")\n";
        
        double pe = spp->establishmentProbability(t, env);
        double gb = spp->growthRate(-1, spp->xb, t, env);
        if (spp->birth_flux_in < 0){    
            double birthFlux = calcSpeciesBirthFlux(s,t) * pe;
            spp->set_ub(birthFlux/(gb+1e-20));
        }
        else{
            spp->calc_boundary_u(gb, pe); // this will set u in boundary cohort
        }

        
        for (int i=0; i<spp->J; ++i){
            double x = spp->getX(i);
            vector<double> g_gx = spp->growthRateGradient(i, x, t, env, control.cm_grad_dx);  // FIXME: x can go
            double dx =  g_gx[0];
            double du = -(spp->mortalityRate(i, x, t, env) + g_gx[1]); 
            if (!use_log_densities) du *= spp->getU(i);
        
            *itr++ = dx;
            *itr++ = du;    
            
            its+=2;
        }

        
        if (spp->n_extra_statevars > 0){
            auto itr_prev = itr;
            spp->getExtraRates(itr);
            assert(distance(itr_prev, itr) == spp->n_extra_statevars*spp->J); 
            its += spp->n_extra_statevars*spp->J;   
        }
        //cout << "---\n";
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

calcRates_EBT()

void Solver::calcRates_EBT	(	double	t,
		std::vector< double >::iterator	S,
		std::vector< double >::iterator	dSdt
	)

Definition at line 8 of file ebt.cpp.

                                                                                         {

    vector<double>::iterator its = S    + n_statevars_system; // Skip system variables
    vector<double>::iterator itr = dSdt + n_statevars_system;

    for (int s = 0; s < species_vec.size(); ++s){
        Species_Base * spp = species_vec[s];

        double   pi0  =  spp->getX(spp->J-1);    // last cohort is pi0, N0
        double   N0   =  spp->getU(spp->J-1);
        //std::cout << "pi = " << pi0 << ", N0 = " << N0 << "\n";

        std::vector<double> g_gx = spp->growthRateGradient(-1, spp->xb, t, env, control.ebt_grad_dx);
        std::vector<double> m_mx = spp->mortalityRateGradient(-1, spp->xb, t, env, control.ebt_grad_dx);
        //std::cout << "g = " << g_gx[0] << ", gx = " << g_gx[1] << "\n";
        //std::cout << "m = " << m_mx[0] << ", mx = " << m_mx[1] << "\n";

        double mb = m_mx[0], mortGrad = m_mx[1], gb = g_gx[0], growthGrad = g_gx[1];    
        
        double birthFlux;
        double pe = spp->establishmentProbability(t, env);
        if (spp->birth_flux_in < 0){    
            birthFlux = calcSpeciesBirthFlux(s,t) * pe;
        }
        else{
            double u0 = spp->calc_boundary_u(gb, pe);
            birthFlux = u0*gb;
        }

        for (int i=0; i<spp->J-1; ++i){ // go down to the second last cohort (exclude boundary cohort)
            double dx = spp->growthRate(i, spp->getX(i), t, env);                           // dx/dt
            double du = -spp->mortalityRate(i, spp->getX(i), t, env) * spp->getU(i);        // du/dt
            //std::cout << "S/C = " << s << "/" << i << " " << spp->getX(i) << " " << dx << " " << du << "\n";
            *itr++ = dx;
            *itr++ = du;
            its += 2;
        }

        // dpi0/dt and dN0/dt
        //std::cout << "S/C = " << s << "/" << "b" << " | pi0/N0 = " << pi0 << " " << N0;
        //if (pi0 <= 0) pi0 = 1e-40;
        //if (N0  <= 0) N0  = 1e-40;
        double dN0  = -mb*N0 - mortGrad*pi0 + birthFlux;
        double dpi0 = gb*N0 + growthGrad*pi0 - mb*pi0;
        //std::cout << " | dpi0/dN0 = " << dpi0 << " " << dN0 << " | mx/gx/mb/gb = " << m_mx[1] << " " << g_gx[1] << " " << m_mx[0] << " " << g_gx[0] << "\n";
        *itr++ = dpi0;
        *itr++ = dN0;
        its += 2;

        if (spp->n_extra_statevars > 0){
            auto itr_prev = itr;
            spp->getExtraRates(itr);
            assert(distance(itr_prev, itr) == spp->n_extra_statevars*spp->J); 
            its += spp->n_extra_statevars*spp->J;   
        }
        
    }

}

Here is the call graph for this function:

Here is the caller graph for this function:

calcRates_FMU()

void Solver::calcRates_FMU	(	double	t,
		std::vector< double >::iterator	S,
		std::vector< double >::iterator	dSdt
	)

Definition at line 11 of file fmu.cpp.

                                                                                         {

    vector<double>::iterator its = S    + n_statevars_system; // Skip system variables
    vector<double>::iterator itr = dSdt + n_statevars_system;
    
    for (int s = 0; s<species_vec.size(); ++s){
        Species_Base* spp = species_vec[s];
        
        // [S S S u u u u u a b c a b c a b c a b c a b c] <--- full SV for species is this
        //        ^ its, itr
        double *U = &(*its); // Since FMU only has U in state, start of species is actually U
        double *dUdt = &(*itr); 
        int J = spp->J;         // xsize of species.
        vector<double> &x = spp->x;
        vector<double> &X = spp->X;
        vector<double> &h = spp->h;

        vector <double> growthArray(J+1);
        growthArray[0] = spp->growthRate(-1, x[0], t, env); // growth rate of boundary cohort
        for (int i=1; i<J+1; ++i) growthArray[i] = spp->growthRateOffset(i-1, x[i], t, env);

        vector <double> u(J+1);
        
        // i=0
        double pe = spp->establishmentProbability(t, env);
        if (spp->birth_flux_in < 0){    
            double birthFlux = calcSpeciesBirthFlux(s,t) * pe;
            u[0] = birthFlux/(growthArray[0]+1e-12); // Q: is this correct? or g(X0,env)? - A: It is g(xb,env) - growthArray[] indexes x, so (0) is xb. Hence correct. 
        }
        else{
            u[0] = spp->calc_boundary_u(growthArray[0], pe);
        }

        // i=1 (calc u1 assuming linear u(x) in first interval)
        u[1] = 2*U[0]-u[0];  // NOTE: for g(x) < 0 this can be calculated with upwind scheme 
        
        for (int i=2; i<J-1; ++i){ // dU[i] ~ u[i+1] <-- U[i],U[i-1], u[i] <-- U[i-1],U[i-2]
            if(growthArray[i] >=0){
                double rMinus = ((U[i]-U[i-1])/(x[i]-x[i-1]))/((U[i-1]-U[i-2]+1e-12)/(x[i-1]-x[i-2]));
                u[i] = U[i-1] + phi(rMinus)*(U[i-1]-U[i-2])*(x[i]-x[i-1])/(x[i+1]-x[i-1]); 
            }   
            else{
                double rPlus  = ((U[i]-U[i-1])/(x[i]-x[i-1]))/((U[i+1]-U[i]+1e-12)/(x[i+1]-x[i]));
                u[i] = U[i] - phi(rPlus)*(U[i+1]-U[i])*(x[i+1]-x[i])/(x[i+2]-x[i]); 
            }
        }
        
        u[J-1] = 2*U[J-2] - u[J-2]; // NOTE: for g(x) > 0 This can be calc with upwind scheme
        u[J] = 2*U[J-1] - u[J-1];

        // [S S S u u u u u a b c a b c a b c a b c a b c] <--- full SV for species is this
        //        ^ itr, its. Therefore, advance 1 at a time while setting dUdt, and advance its by 3*5 after.
        for (int i=0; i<spp->J; ++i){ // dU[i] ~ u[i+1] <-- U[i],U[i-1], u[i] <-- U[i-1],U[i-2]
            dUdt[i] = -spp->mortalityRate(i,X[i], t, env)*U[i] - (growthArray[i+1]*u[i+1] - growthArray[i]*u[i])/h[i];
            ++itr; ++its;
        }
        if (spp->n_extra_statevars > 0){
            auto itr_prev = itr;
            spp->getExtraRates(itr);
            assert(distance(itr_prev, itr) == spp->n_extra_statevars*spp->J);
            its += spp->n_extra_statevars*spp->J;   
        }
            

    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

calcRates_iFMU()

void Solver::calcRates_iFMU	(	double	t,
		std::vector< double >::iterator	S,
		std::vector< double >::iterator	dSdt
	)

Definition at line 5 of file ifmu.cpp.

                                                                                          {
    vector<double>::iterator its = S    + n_statevars_system; // Skip system variables
    vector<double>::iterator itr = dSdt + n_statevars_system;
    
    for (int s = 0; s<species_vec.size(); ++s){
        auto spp = species_vec[s];
        
        its += spp->J;// skip u 
        for (int i=0; i<spp->J; ++i) *itr++ = 0; // set du/dt to 0 
    
        if (spp->n_extra_statevars > 0){
            auto itr_prev = itr;
            spp->getExtraRates(itr); // TODO/FIXME: Does calc of extra rates need t and env?
            assert(distance(itr_prev, itr) == spp->n_extra_statevars*spp->J);
            its += spp->n_extra_statevars*spp->J;   
        }
    }

}

Here is the caller graph for this function:

calcSpeciesBirthFlux()

double Solver::calcSpeciesBirthFlux	(	int	k,
		double	t
	)

Definition at line 352 of file solver.cpp.

                                                  {
    std::cout << "calc birthFlux...\n";
    auto spp = species_vec[k];  
    auto newborns_production = [this, spp](int i, double _t){
        double b1 = spp->birthRate(i, spp->getX(i), _t, env);
        return b1;  
    }; 
    double birthFlux = integrate_x(newborns_production, t, k);
    return birthFlux;   
}

Here is the call graph for this function:

Here is the caller graph for this function:

copyCohortsToState()

void Solver::copyCohortsToState

(

)

Definition at line 297 of file solver.cpp.

                               {
    std::cout << "state <--- cohorts\n";
    vector<double>::iterator it = state.begin() + n_statevars_system; // no need to copy system state
    
    for (int k=0; k<species_vec.size(); ++k){
        Species_Base* s = species_vec[k];
        
        if (method == SOLVER_FMU || method == SOLVER_IFMU){
            for (size_t i=0; i<s->J; ++i){
                *it++ = s->getU(i);
            }
        }
        if (method == SOLVER_CM){
            for (size_t i=0; i<s->J; ++i){
                double X = s->getX(i); 
                double U = s->getU(i);
                U = (use_log_densities)? log(U) : U;    // log(u) in state 
                *it++ = X;  // set x to state
                *it++ = U;  // set u to state
            }
        }
        if (method == SOLVER_EBT){
            // x, u for boundary and internal cohorts
            for (size_t i=0; i<s->J; ++i){
                double X = s->getX(i); 
                double U = s->getU(i);
                *it++ = X; 
                *it++ = U;
            }
        }

        if (s->n_extra_statevars > 0){  // If extra state variables have been requested, initialize them
            auto it_prev = it;
            s->copyCohortsExtraToState(it);
            assert(distance(it_prev, it) == s->n_extra_statevars*s->J); 
        }
    }
}

Here is the call graph for this function:

copyStateToCohorts()

void Solver::copyStateToCohorts

(

std::vector< double >::iterator

state_begin

)

Definition at line 256 of file solver.cpp.

                                                                    {
    std::cout << "state ---> cohorts\n";
    std::vector<double>::iterator it = state_begin + n_statevars_system; // no need to copy system state
    
    for (int k=0; k<species_vec.size(); ++k){
        Species_Base* s = species_vec[k];
        
        s->set_xb(s->xb); // Important: "touch" boundary cohort to trigger a precompute. Note that setSize() triggers it.
        if (method == SOLVER_FMU || method == SOLVER_IFMU){
            for (size_t i=0; i<s->J; ++i){
                s->setU(i,*it++);
            }
        }
        if (method == SOLVER_CM){
            for (size_t i=0; i<s->J; ++i){
                double X = *it++;   // get x from state
                double U = *it++;   // get u from state
                U = (use_log_densities)? exp(U) : U;    // u in cohorts 
                s->setX(i,X); 
                s->setU(i,U);
            }
        }
        if (method == SOLVER_EBT){
            // x, u for boundary and internal cohorts
            for (size_t i=0; i<s->J; ++i){
                double X = *it++; 
                double U = *it++;
                s->setX(i,X); 
                s->setU(i,U);
            }
        }

        if (s->n_extra_statevars > 0){  // If extra state variables have been requested, initialize them
            auto it_prev = it;
            s->copyExtraStateToCohorts(it);
            assert(distance(it_prev, it) == s->n_extra_statevars*s->J); 
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

getDensitySpecies_EBT()

std::vector< double > Solver::getDensitySpecies_EBT	(	int	k,
		std::vector< double >	breaks
	)

Definition at line 430 of file solver.cpp.

                                                                           {
    auto spp = species_vec[k];

    if (method == SOLVER_EBT){
        double xm = spp->getX(0)+1e-6;

        vector<point> points(breaks.size()-1);

        // assuming breaks are sorted ascending
        // cohorts are sorted descending
        int current_interval = breaks.size()-2;
        for (int i=0; i<spp->J; ++i){ // loop over all cohorts except boundary cohort
            double x = spp->getX(i);
            double N = spp->getU(i);
            if (i == spp->J-1) x = spp->xb + x/(N+1e-12); // real-ize x if it's boundary cohort
            
            while(breaks[current_interval]>x) --current_interval; // decrement interval until x fits
            //cout << current_interval << ", x = " << x << "(" << N << ") in [" << breaks[current_interval] << ", " << breaks[current_interval+1] << "]\n"; cout.flush();
            if (N>0){
                points[current_interval].abund += N;
                points[current_interval].count += 1;
                points[current_interval].xmean += N*x;
            }
        }

        
        for (int i=0; i<points.size(); ++i) if (points[i].count>0) points[i].xmean /= points[i].abund;
        
        // remove 0-count points
        auto pred = [this](const point& p) -> bool {
            return p.count == 0; 
        };

        auto pend = std::remove_if(points.begin(), points.end(), pred);
        points.erase(pend, points.end());

        //cout << "mean x and abund (removed):\n";
        //for (int i=0; i<points.size(); ++i) cout << i << "\t" << points[i].count << "\t" << points[i].xmean << "\t" << points[i].abund << "\n";   
        //cout << "--\n";

        if (points.size() > 2){
            vector<double> h(points.size());
            h[0] = (points[1].xmean+points[0].xmean)/2 - spp->xb;
            for (int i=1; i<h.size()-1; ++i) h[i] = (points[i+1].xmean - points[i-1].xmean)/2;
            h[h.size()-1] = xm - (points[h.size()-1].xmean+points[h.size()-2].xmean)/2;

            vector <double> xx, uu;
            xx.reserve(points.size());
            uu.reserve(points.size());
            for (int i=0; i<points.size(); ++i){
                xx.push_back(points[i].xmean);
                uu.push_back(points[i].abund / h[i]);
            }
            
            Spline spl;
            spl.splineType = Spline::LINEAR; //Spline::CONSTRAINED_CUBIC;
            spl.extrapolate = Spline::ZERO;
            spl.set_points(xx, uu);
            
            vector <double> dens;
            dens.reserve(points.size());
            for (int i=0; i<breaks.size(); ++i){
                dens.push_back(spl.eval(breaks[i]));            
            }
            
            return dens;
        }
        else return vector<double>(breaks.size(), 0);
    }
    else {
        throw std::runtime_error("This function can only be called for the EBT solver");
    }

}

Here is the call graph for this function:

initialize()

void Solver::initialize

(

)

Definition at line 194 of file solver.cpp.

                       {
    
    vector<double>::iterator it = state.begin() + n_statevars_system; // TODO: replace with init_sState() 
    
    for (int k=0; k<species_vec.size(); ++k){
        Species_Base* s = species_vec[k];
    
        // Boundary cohort is not in state, but used as a reference.    
        s->set_xb(s->xb); // set x of boundary cohort 
        s->set_ub(0);     // set initial density of boundary cohort to 0.
        
        // set birth time for each cohort to current_time
        for (int i=0; i<s->J; ++i) s->set_birthTime(i, current_time);

        // set x, u for all cohorts
        if (method == SOLVER_FMU || method == SOLVER_IFMU){
            for (size_t i=0; i<s->J; ++i){
                double X = s->X[i];         // for FMU, X is the midpoint of the cell edges
                double U = s->init_density(i, X, env); 
                s->setX(i,X); 
                s->setU(i,U);
                *it++ = U;      // u in state (only)
            }
        }
        if (method == SOLVER_CM){
            for (size_t i=0; i<s->J; ++i){
                double X = s->x[i];
                double U = s->init_density(i, X, env); 
                s->setX(i,X); 
                s->setU(i,U);
                *it++ = X;                                  // x in state
                *it++ = (use_log_densities)? log(U) : U;    // u in state 
            }
        }
        if (method == SOLVER_EBT){
            // x, u for internal cohorts in state and it cohorts
            for (size_t i=0; i<s->J-1; ++i){
                double X = (s->x[i]+s->x[i+1])/2.0;         
                double U = s->init_density(i, X, env)*(s->x[i]-s->x[i+1]); 
                s->setX(i,X); 
                s->setU(i,U);
                *it++ = X;  // x in state
                *it++ = U;  // u in state 
            }
            // set pi0, N0 as x, u for the last cohort. This scheme allows using this last cohort with xb+pi0 in integrals etc 
            *it++ = 0; *it++ = 0;
            s->setX(s->J-1,0); // TODO: should this be set to xb for init_state and set to 0 again later? maybe not, as init_state is not expected to be x dependent
            s->setU(s->J-1,0);      
        }

        // initialize extra state for each cohort and copy it to state
        s->initAndCopyExtraState(current_time, env, it);
        //if (s->n_extra_statevars > 0){  // FIXME: maybe redundant
            //auto it_prev = it;
            //s->init_ExtraState(it);  // this also inits the extra state of boundary cohort, but without advancing the iterator
            //assert(distance(it_prev, it) == s->n_extra_statevars*s->J); 
        //}

    }
}

Here is the call graph for this function:

integrate_wudx_above()

template<typename wFunc >

double Solver::integrate_wudx_above	(	wFunc	w,
		double	t,
		double	xlow,
		int	species_id
	)

Computes the partial integral

\[I = \int_{x_{low}}^{x_m} w(z,t)u(z)dz\]

for the specified species.

Parameters

w	A function or function-object of the form w(int i, double t) that returns a double. It can be a lambda. This function should access the 'i'th cohort and compute the weight from the cohort's properties. The function w should be able access to the Solver in order to access cohorts.
t	The current time (corresponding to the current physiological state). This will be passed to w to allow the weights to be a direct function of time.
xlow	The lower limit of the integral
species_id	The id of the species for which the integral should be computed

the computation of this integral depends on the solver method. For different solvers, the integral is defined as follows:

FMU: \(\quad I = \sum_{i=i_0}^J h_i w_i u_i\)

EBT: \(\quad I = \sum_{i=i_0}^J w_i N_i\), with \(x_0 = x_b + \pi_0/N_0\)

CM : \(\quad I = \sum_{i=i_0}^J h_i (w_{i+1}u_{i+1}+w_i u_i)/2\)

where \(i_0 = \text{argmax}(x_i \le x_{low})\), \(h_i = x_{i+1}-x_i\), and \(w_i = w(x_i)\).

If interpolation is turned on, \(h_{i_0}=x_{i_0+1}-x_{low}\), whereas \(u(x_{low})\) is set to \(u(x_{i_0})\) in FMU and calculated by bilinear interpolation in CM (See Figure). Interpolation does not play a role in EBT.

In the CM method, the density of the boundary cohort is obtained from the boundary condition of the PDE: \(u_b=B/g(x_b)\), where \(B\) is the input flux of newborns. In the current implementation of the CM method, \(B\) must be set to a constant. In future implementations which may allow real-time calculation of \(B\), \(u_b\) must be calculated recurvisely by solving \(u_b = B(u_b)/g(x_b)\), where \(B(u_b) = \int_{x_b}^{x_m} f(x)u(x)dx\).

Definition at line 44 of file size_integrals.tpp.

                                                                                 {

    Species_Base* spp = species_vec[species_id];

    // cohorts are inserted at the end, hence sorted in descending order
    // FIXME: should there be an optional "sort_cohorts" control parameter? Maybe some freak models are such that cohorts dont remain in sorted order?
    if (method == SOLVER_CM){
        // integrate using trapezoidal rule 
        // Note, new cohorts are inserted at the end, so x will be descending
        bool integration_completed = false;
        double I = 0;
        double u_hi = spp->getU(0); 
        double x_hi = spp->getX(0);
        double f_hi = w(0, t)*u_hi;
        
        for (int i=1; i<spp->J; ++i){  // going from high ----> low
            double u_lo = spp->getU(i); 
            double x_lo = spp->getX(i);
            double f_lo = w(i, t)*u_lo;
    
            if (x_lo <= xlow){  // check if xlow falls within the current interval. 
                // * if interpolation is on, stop exactly at xlow, else take the full interval
                double f = (control.integral_interpolate)? f_lo + (f_hi-f_lo)/(x_hi-x_lo)*(xlow - x_lo) : f_lo;  
                double x = (control.integral_interpolate)? xlow : x_lo;
                I += (x_hi - x) * (f_hi + f);
                integration_completed = true;
            }
            else{
                I += (x_hi - x_lo) * (f_hi + f_lo);
            }
        
            x_hi = x_lo;
            f_hi = f_lo;
            if (integration_completed) break;
        }
        
        //if (integration_completed) assert(x_hi < xlow);
        //if (!integration_completed) assert(x_hi >= xlow || spp.J == 1);
        
        // if integration has not completed, continue to the boundary interval  
        if (spp->J == 1 || !integration_completed){  
            // boundary at xb
            double u0 = spp->get_boundary_u();
            double x_lo = spp->xb;
            if (x_hi > x_lo){ 
                double f_lo =  w(-1, t)*u0; // -1 is boundary cohort
                double f = (control.integral_interpolate)? f_lo + (f_hi-f_lo)/(x_hi-x_lo)*(xlow - x_lo) : f_lo;  
                double x = (control.integral_interpolate)? xlow : x_lo;
                I += (x_hi-x)*(f_hi+f);
            }
        }
        
        return I*0.5;
    }

    else if (method == SOLVER_FMU || method == SOLVER_IFMU){
        //if (xlow < spp->xb) throw std::runtime_error("integral lower bound must be >= xb");
        if (xlow > spp->x[spp->J]) return 0;  // if xlow is above the maximum size, integral is 0

        // integrate using midpoint quadrature rule
        double I=0;
        for (int i=spp->J-1; i>=0; --i){  // in FMU, cohorts are sorted ascending
            if (spp->x[i] <= xlow){ // check if last interval is reached 
                double h = (control.integral_interpolate)?  spp->x[i+1]-xlow  :  spp->h[i];
                I += h * w(i,t) * spp->getU(i);
                break;
            }
            else{
                I += spp->h[i] * w(i,t) * spp->getU(i);
            }
        }
        
        return I;
    }
    
    else if (method == SOLVER_EBT){
        // set up cohorts to integrate
        // backup pi0, N0 from last (youngest) cohort <-- cohorts are sorted descending
        double   pi0  =  spp->getX(spp->J-1);
        double   N0   =  spp->getU(spp->J-1);

        // real-ize pi0-cohort with actual x0 value
        double x0 = spp->xb + pi0/(N0+1e-12);
        spp->setX(spp->J-1, x0);

        // calculate integral
        double I = 0;
        for (int i=0; i<spp->J; ++i){  // in EBT, cohorts are sorted descending
            if (spp->getX(i) < xlow) break; // if X == xlow, we still include it in the intgral
            else I += w(i, t)*spp->getU(i);
        }
        
        // restore the original pi0-cohort
        spp->setX(spp->J-1, pi0);

        return I;
    }
    
    else{
        throw std::runtime_error("Unsupported solver method");
    }
}

Here is the call graph for this function:

integrate_x()

template<typename wFunc >

double Solver::integrate_x	(	wFunc	w,
		double	t,
		int	species_id
	)

Computes the integral

\[I = \int_{x_b}^{x_m} w(z,t)u(z)dz\]

for the specified species. For details, see integrate_wudx_above.

Definition at line 156 of file size_integrals.tpp.

                                                           {
    Species_Base* spp = species_vec[species_id];

    if (method == SOLVER_FMU || method == SOLVER_IFMU){
        // integrate using midpoint quadrature rule
        double I=0;
        for (unsigned int i=0; i<spp->J; ++i){
            I += spp->h[i]*w(i, t)*spp->getU(i);  // TODO: Replace with std::transform after profiling
        }
        return I;
    }
    
    else if (method == SOLVER_EBT){
        // integrate using EBT rule (sum over cohorts)
        
        // backup pi0, N0 from last cohort 
        double   pi0  =  spp->getX(spp->J-1);
        double   N0   =  spp->getU(spp->J-1);

        // real-ize pi0-cohort with actual x0 value
        double x0 = spp->xb + pi0/(N0+1e-12);
        spp->setX(spp->J-1, x0);

        // calculate integral
        double I = 0;
        for (int i=0; i<spp->J; ++i) I += w(i, t)*spp->getU(i);
        
        // restore the original pi0-cohort
        spp->setX(spp->J-1, pi0);

        return I;
    }
    
    if (method == SOLVER_CM){
        // integrate using trapezoidal rule 
        // Note, new cohorts are inserted at the beginning, so x will be ascending
        double I = 0;
        double u_hi = spp->getU(0); 
        double x_hi = spp->getX(0);
        double f_hi = w(0, t)*u_hi;

        for (int i=1; i<spp->J; ++i){
            double u_lo = spp->getU(i); 
            double x_lo = spp->getX(i);
            double f_lo = w(i, t)*u_lo;
    
            I += (x_hi - x_lo) * (f_hi + f_lo);
            x_hi = x_lo;
            f_hi = f_lo;
        }
        
        // boundary at xb
        double u0 = spp->get_boundary_u();
        double x_lo = spp->xb;
        double f_lo =  w(-1, t)*u0; // -1 is boundary cohort
        I += (x_hi-x_lo)*(f_hi+f_lo);
        
        return I*0.5;
    }
    
    else{
        throw std::runtime_error("Unsupported solver method");
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

maxSize()

double Solver::maxSize

(

)

Definition at line 158 of file solver.cpp.

                      {
    double maxsize = 0;
    for (auto spp : species_vec) maxsize = std::max(maxsize, spp->get_maxSize());
    return maxsize;
}

n_species()

int Solver::n_species

(

)

Definition at line 153 of file solver.cpp.

                     {
    return species_vec.size();
}

newborns_out()

vector< double > Solver::newborns_out

(

double

t

)

Definition at line 364 of file solver.cpp.

                                           {  // TODO: make recompute env optional
    // update Environment from latest state
    //copyStateToCohorts(state.begin());  // not needed here because this is done in afterStep or upon cohorts update
    //env->computeEnv(t, this, state.begin(), rates.begin());
    updateEnv(t, state.begin(), rates.begin()); // this will trigger a precompute
//  for (int k = 0; k<species_vec.size(); ++k) preComputeSpecies(k,t);  

    vector<double> b_out;
    for (int k=0; k<species_vec.size(); ++k){   
        // calculate birthflux
        // [solved] wudx doesnt work here. Why?? - works now, no idea why it was not working earlier!
        //auto newborns_production = [this, k](int i, double t){
            //double z = species_vec[k]->getX(i);
            //double b = species_vec[k]->birthRate(i,z,t,env);  
            //return b; 
        //}; 
        //double birthFlux = integrate_x(newborns_production, t, k);
        //double birthFlux = integrate_wudx_above(newborns_production, t, 0, k);
        double birthFlux = calcSpeciesBirthFlux(k, t);
        b_out.push_back(birthFlux);
    }
    return b_out;
}

Here is the call graph for this function:

Here is the caller graph for this function:

print()

void Solver::print

(

)

Definition at line 167 of file solver.cpp.

                  {
    std::cout << ">> SOLVER \n";
    string types[] = {"FMU", "MMU", "CM", "EBT", "Implicit FMU"};
    std::cout << "+ Type: " << types[method] << std::endl;

    std::cout << "+ State size = " << state.size() << "\n";
    std::cout << "+ Rates size = " << rates.size() << "\n";
    std::cout << "+ Species:\n";
    std::cout.flush();
    for (int i=0; i<species_vec.size(); ++i) {
        std::cout << "Sp (" << i << "):\n";
        species_vec[i]->print();
    }
    
}

removeCohort_CM()

void Solver::removeCohort_CM

(

)

Definition at line 85 of file cm.cpp.

                            {

    for (auto& spp : species_vec){
        if (spp->J > control.max_cohorts) // remove a cohort if number of cohorts in the species exceeds a threshold
            spp->removeDensestCohort();
    }
    
    resizeStateFromSpecies();
    copyCohortsToState();
    cout << "........................................\n";
    
}

Here is the caller graph for this function:

removeDeadCohorts_EBT()

void Solver::removeDeadCohorts_EBT

(

)

Definition at line 99 of file ebt.cpp.

                                  {

    for (auto spp : species_vec){
        spp->removeDeadCohorts(control.ebt_ucut);   
    }

    resizeStateFromSpecies();
    copyCohortsToState();

}

Here is the caller graph for this function:

resetState()

void Solver::resetState

(

double

t0 = 0

)

Definition at line 123 of file solver.cpp.

                                {  // FIXME: This is currently redundant, and needs to be improved with reset of both state and cohorts for a true reset of state
    current_time = t0;
    odeStepper.reset(t0, control.ode_eps, control.ode_eps); // = RKCK45<vector<double>> (0, control.ode_eps, control.ode_initial_step_size);  // this is a cheap operation, but this will empty the internal containers, which will then be (automatically) resized at next 1st ODE step. Maybe add a reset function to the ODE stepper? 

    std::fill(state.begin(), state.end(), 0); 
    std::fill(rates.begin(), rates.end(), -999); // DEBUG
}

Here is the call graph for this function:

resizeStateFromSpecies()

void Solver::resizeStateFromSpecies

(

)

Definition at line 132 of file solver.cpp.

                                   {
    int state_size_new = n_statevars_system;
    for (auto& spp : species_vec){
        state_size_new += spp->J*(n_statevars_internal + spp->n_extra_statevars);
    }
    
    state.resize(state_size_new, -999);
    rates.resize(state_size_new, -999);
}

setEnvironment()

void Solver::setEnvironment

(

EnvironmentBase *

_env

)

Definition at line 143 of file solver.cpp.

                                                 {
    env = _env;
}

step_to() [1/2]

template<typename AfterStepFunc >

void Solver::step_to	(	double	tstop,
		AfterStepFunc &	afterStep_user
	)

Parameters

tstop	The time until which the ODE solver should be stepped. The stepper will stop exactly at tstop, for which the final step size is truncated if necessary.
afterStep_user	A function of the form f(double t) to be called after every successful ODE step.

Note

1. current_time is updated by the ODE solver at every (internal) step

2. After the last ODE step, the state vector is updated but cohorts still hold an intenal ODE state (y+k5*h etc). normally, this will not be a problem because state will be copied to cohorts before rates call of the next timestep. But since add/remove cohort after step_to will rewrite the state from cohorts, the updated state vector will be lost. To avoid this, we should ensure that the state is copied to cohorts after every successful ODE step (or at least after completion of step_to). This is achieved by the afterStep function

Definition at line 15 of file solver.tpp.

                                                               {
    // do nothing if tstop is <= current_time
    if (tstop <= current_time) return;
    
    auto after_step = [this, afterStep_user](double t, std::vector<double>::iterator S){
        std::cout << "After step: t = " << t << "\n";
        copyStateToCohorts(S);
        afterStep_user(t);
    };

    
    if (method == SOLVER_FMU){  
        auto derivs = [this](double t, std::vector<double>::iterator S, std::vector<double>::iterator dSdt, void* params){
            copyStateToCohorts(S); 
            updateEnv(t, S, dSdt);
            calcRates_FMU(t, S, dSdt);
        };
        
        // integrate
        odeStepper.step_to(tstop, current_time, state, derivs, after_step); // rk4_stepsize is only used if method is "rk4"
    }
    
    
    if (method == SOLVER_IFMU){ 
        while (current_time < tstop){
            double dt = std::min(control.ode_ifmu_stepsize, tstop-current_time);
            
            //copyStateToCohorts(state.begin()); // not needed here because it is called by the odestepper below
            updateEnv(current_time, state.begin(), rates.begin());
            std::vector<double> rates_prev(rates.begin(), rates.begin()+n_statevars_system);  // save system variable rates
            
            // use implicit stepper to advance u
            stepU_iFMU(current_time, state, rates, dt);
            // current_time += dt; // not needed here, as current time is advanced by the ODE stepper below.
            
            if (n_statevars_system > 0){
                copyStateToCohorts(state.begin());   // copy updated u to cohorts 
                updateEnv(current_time, state.begin(), rates.begin());  // recompute env with updated u
                // FIXME: use fully implicit stepper here?
                for (int i=0; i<n_statevars_system; ++i){
                    state[i] += (rates_prev[i]+rates[i])/2*dt;  // use average of old and updated rates for stepping system vars
                }
            }

            // use the ODE-stepper for other state variables
            // this stepper is called even if there are no extra state variables, so copyStateToCohorts is accomplished here
            auto derivs = [this](double t, std::vector<double>::iterator S, std::vector<double>::iterator dSdt, void* params){
                copyStateToCohorts(S);
                // precompute and env computation is not needed here, because it depends on x and u, which are not updated by the solver.
                calcRates_iFMU(t,S,dSdt);
            };
            // this step below will do afterstep. FIXME: But what if there are no extra state variables? 
            odeStepper.step_to(current_time+dt, current_time, state, derivs, after_step); // rk4_stepsize is only used if method is "rk4"
    
        }

    }
    
    if (method == SOLVER_EBT){
        auto derivs = [this](double t, std::vector<double>::iterator S, std::vector<double>::iterator dSdt, void* params){
            copyStateToCohorts(S);
            updateEnv(t, S, dSdt);
            calcRates_EBT(t, S, dSdt);
        };
        
        // integrate 
        odeStepper.step_to(tstop, current_time, state, derivs, after_step); // rk4_stepsize is only used if method is "rk4"
        
        // update cohorts
        removeDeadCohorts_EBT();
        addCohort_EBT();  // Add new cohort if N0 > 0. Add after removing dead ones otherwise this will also be removed. 
    }
    
    
    if (method == SOLVER_CM){
        auto derivs = [this](double t, std::vector<double>::iterator S, std::vector<double>::iterator dSdt, void* params){
            std::cout << "derivs()\n";
            copyStateToCohorts(S);  // this triggers precompute
            updateEnv(t, S, dSdt);
            calcRates_CM(t, S, dSdt);
        };
        
        // integrate 
        odeStepper.step_to(tstop, current_time, state, derivs, after_step); // rk4_stepsize is only used if method is "rk4"
        
        // update cohorts
        if (control.update_cohorts){
            addCohort_CM();     // add before so that it becomes boundary cohort and first internal cohort can be (potentially) removed
            removeCohort_CM();
        }
        //env->computeEnv(current_time, this); // is required here IF rescaleEnv is used in derivs
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

step_to() [2/2]

void Solver::step_to

(

double

tstop

)

Definition at line 338 of file solver.cpp.

                                {
    auto func = [](double t){};
    step_to(tstop, func);
}

Here is the call graph for this function:

stepU_iFMU()

void Solver::stepU_iFMU	(	double	t,
		std::vector< double > &	S,
		std::vector< double > &	dSdt,
		double	dt
	)

Definition at line 25 of file ifmu.cpp.

                                                                                   {
        
    vector<double>::iterator its = S.begin() + n_statevars_system; // Skip system variables
    
    // 1. Take implicit step for U
    for (int s = 0; s<species_vec.size(); ++s){
        Species_Base* spp = species_vec[s];
        
        // [S S S u u u u u a b c a b c a b c a b c a b c] <--- full SV for species is this
        //        ^ its, itr
        double *U = &(*its); // Since FMU only has U in state, start of species is actually U
        int J = spp->J;         // xsize of species.
        vector<double> &h = spp->h;
        
        vector <double> growthArray(J);
        for (int i=0; i<J; ++i) growthArray[i] = spp->growthRate(i, spp->getX(i), t, env);
        
        double birthFlux;
        double pe = spp->establishmentProbability(t, env);
        if (spp->birth_flux_in < 0){    
            birthFlux = calcSpeciesBirthFlux(s,t) * pe;
        }
        else{
            birthFlux = spp->calc_boundary_u(growthArray[0], pe)*growthArray[0];
        }
        
        //cout << t << "\t" << birthFlux/growthArray[0] << " -> " << calcSpeciesBirthFlux(s,t)/growthArray[0] << "\n";
        double B0  = 1 + dt/h[0]*growthArray[0] + dt*spp->mortalityRate(0, spp->getX(0), t, env);
//      std::cout << "B0 = " << B0 << endl; 
        double C0 = spp->getU(0) + dt/h[0]*birthFlux;
        U[0] = C0/B0;
//      std::cout << "U0 = " << U[0] << endl; 

        for (int w = 1; w < J; ++w){
            double Aw = -growthArray[w-1]*dt/h[w];
            double Bw  = 1 + dt/h[w]*growthArray[w] + dt*spp->mortalityRate(w, spp->getX(w), t, env);
            double Cw = spp->getU(w);

            U[w] = (Cw - Aw*U[w-1])/Bw;
        }
        
        its += J*(1+spp->n_extra_statevars);

    }
    
}

Here is the call graph for this function:

Here is the caller graph for this function:

u0_out()

vector< double > Solver::u0_out

(

double

t

)

Definition at line 391 of file solver.cpp.

                                     {
    vector <double> u0out;
    vector <double> newbornsout = newborns_out(t);
    for (int k=0; k < species_vec.size(); ++k){
        //species_vec[k]->preCompute(-1, t, env); // not req because precomputeAllCohorts called in newborns_out() precomputes BC too.  
        u0out.push_back(newbornsout[k]/species_vec[k]->growthRate(-1, species_vec[k]->xb, t, env));
    }
    return u0out;
}

Here is the call graph for this function:

updateEnv()

void Solver::updateEnv	(	double	t,
		std::vector< double >::iterator	S,
		std::vector< double >::iterator	dSdt
	)

Definition at line 344 of file solver.cpp.

                                                                                           {
    std::cout << "update Env...\n";
    for (auto spp : species_vec) spp->triggerPreCompute();
    env->computeEnv(t, this, S, dSdt);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

cm_grad_dx

double Solver::cm_grad_dx = 1e-6

Definition at line 37 of file solver.h.

control

struct { ... } Solver::control

convergence_eps

double Solver::convergence_eps = 1e-6

Definition at line 36 of file solver.h.

current_time

double Solver::current_time

Definition at line 27 of file solver.h.

ebt_grad_dx

double Solver::ebt_grad_dx = 1e-6

Definition at line 41 of file solver.h.

ebt_ucut

double Solver::ebt_ucut = 1e-10

Definition at line 40 of file solver.h.

env

EnvironmentBase* Solver::env

Definition at line 24 of file solver.h.

ifmu_centered_grids

bool Solver::ifmu_centered_grids = true

Definition at line 45 of file solver.h.

integral_interpolate

bool Solver::integral_interpolate = true

Definition at line 46 of file solver.h.

max_cohorts

int Solver::max_cohorts = 500

Definition at line 39 of file solver.h.

method

PSPM_SolverType Solver::method

private

Definition at line 17 of file solver.h.

n_statevars_internal

int Solver::n_statevars_internal = 0

private

Definition at line 19 of file solver.h.

n_statevars_system

int Solver::n_statevars_system = 0

private

Definition at line 20 of file solver.h.

ode_eps

double Solver::ode_eps = 1e-6

Definition at line 34 of file solver.h.

ode_ifmu_stepsize

double Solver::ode_ifmu_stepsize = 0.1

Definition at line 44 of file solver.h.

ode_initial_step_size

double Solver::ode_initial_step_size = 1e-6

Definition at line 35 of file solver.h.

ode_rk4_stepsize

double Solver::ode_rk4_stepsize = 0.1

Definition at line 43 of file solver.h.

odeStepper

OdeSolver Solver::odeStepper

Definition at line 23 of file solver.h.

rates

std::vector<double> Solver::rates

Definition at line 29 of file solver.h.

species_vec

std::vector<Species_Base*> Solver::species_vec

Definition at line 31 of file solver.h.

state

std::vector<double> Solver::state

Definition at line 28 of file solver.h.

update_cohorts

bool Solver::update_cohorts = true

Definition at line 38 of file solver.h.

use_log_densities

bool Solver::use_log_densities = true

Definition at line 49 of file solver.h.

---

The documentation for this class was generated from the following files:

• include/solver.h
• src/cm.cpp
• src/ebt.cpp
• src/fmu.cpp
• src/ifmu.cpp
• src/size_integrals.tpp
• src/solver.cpp
• src/solver.tpp