size_integrals.tpp

Go to the documentation of this file.

//             _xm 
// Calculate _/ w(z,t)u(z)dz
//         xlow
// implementation from orig plant model 
// ----
// I += (x_hi - x_lo) * (f_hi + f_lo);
// x_hi = x_lo;
// f_hi = f_lo;
// if (x_lo < xlow) break;
// ---- 
template<typename wFunc>
double Solver::integrate_wudx_above(wFunc w, double t, double xlow, int species_id){

    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");
    }
}





//             _xm 
// Calculate _/ w(z,t)u(z,t)dz
//         xb
// This will eventually be replaced with a call to integrate_wudx_above
template<typename wFunc>
double Solver::integrate_x(wFunc w, double t, int species_id){
    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");
    }
}