C++ CVode函数代码示例

int cvode_eval(solver_props *props, unsigned int modelid){
  cvode_mem *mem = props->mem;
  mem = &mem[modelid];

  // Stop the solver if the stop time has been reached
  props->running[modelid] = (props->time[modelid] + props->timestep) < props->stoptime;
  if(!props->running[modelid])
    return 0;

  // if a positive dt is specified, then we will have this function return after it reaches the next time point,
  // otherwise, it will just run one iteration and return
  if(props->timestep > 0) {
    // Reinitialize the function at each step
    if(CVodeReInit(mem->cvmem, props->time[modelid], mem->y0) != CV_SUCCESS) {
      PRINTF( "CVODE failed to reinitialize");
    }

    CDATAFORMAT stop_time = MIN(props->time[modelid] + props->timestep, props->stoptime);
    mem->first_iteration = TRUE;
    if(CVode(mem->cvmem, stop_time, mem->y0, &(props->next_time[modelid]), CV_NORMAL) != CV_SUCCESS){
      PRINTF( "CVODE failed to make a fixed step in model %d.\n", modelid);
      return 1;
    }
  }
  else {
    mem->first_iteration = TRUE;
    if(CVode(mem->cvmem, props->stoptime, mem->y0, &(props->next_time[modelid]), CV_ONE_STEP) != CV_SUCCESS){
      PRINTF( "CVODE failed to make a step in model %d.\n", modelid);
      return 1;
    }
  }

  return 0;
}

 void CVodeInt::integrate(double tout)
 {
   double t;
   int flag;
   flag = CVode(m_cvode_mem, tout, nv(m_y), &t, NORMAL);
   if (flag != SUCCESS) 
     throw CVodeErr(" CVode error encountered.");
 }

 double CVodesIntegrator::step(double tout)
 {
   double t;
   int flag;
   flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP);
   if (flag != CV_SUCCESS) 
     throw CVodesErr(" CVodes error encountered.");
   return t;
 }

void CVodesIntegrator::integrate(double tout)
{
    int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);
    if (flag != CV_SUCCESS) {
        throw CVodesErr("CVodes error encountered. Error code: " + int2str(flag) + "\n" + m_error_message +
                        "\nComponents with largest weighted error estimates:\n" + getErrorInfo(10));
    }
    m_sens_ok = false;
}

double CVodesIntegrator::step(double tout)
{
    int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_ONE_STEP);
    if (flag != CV_SUCCESS) {
        throw CanteraError("CVodesIntegrator::step",
            "CVodes error encountered. Error code: {}\n{}\n"
            "Components with largest weighted error estimates:\n{}",
            flag, m_error_message, getErrorInfo(10));

    }
    m_sens_ok = false;
    return m_time;
}

int ode_solver_solve(ode_solver* solver, const double t, double* y, double* tout){
  
  double lTout;
  
  /* Advance the solution */
  NV_DATA_S(solver->y) = y;
  int flag = CVode(solver->cvode_mem,t, solver->y, &lTout, CV_NORMAL);
  
  if(tout)
    *tout = lTout;
  
  return flag;
}

int SOLVER(cvode, eval, TARGET, SIMENGINE_STORAGE, cvode_mem *mem, unsigned int modelid) {
  // Stop the solver if the stop time has been reached
  mem->props->running[modelid] = mem->props->time[modelid] < mem->props->stoptime;
  if(!mem->props->running[modelid])
    return 0;

  mem[modelid].first_iteration = TRUE;
  if(CVode(mem[modelid].cvmem, mem[modelid].props->stoptime, ((N_Vector)(mem[modelid].y0)), &(mem[modelid].props->time[modelid]), CV_ONE_STEP) != CV_SUCCESS){
    fprintf(stderr, "CVODE failed to make a step in model %d.\n", modelid);
    return 1;
  }

  return 0;
}

real Solver::run(real tout, int &ncalls, real &rhstime)
{
#ifdef CHECK
  int msg_point = msg_stack.push("Running solver: solver::run(%e)", tout);
#endif

  MPI_Barrier(MPI_COMM_WORLD);

  rhs_wtime = 0.0;
  rhs_ncalls = 0;

  pre_Wtime = 0.0;
  pre_ncalls = 0.0;

  int flag = CVode(cvode_mem, tout, uvec, &simtime, CV_NORMAL);
  
  ncalls = rhs_ncalls;
  rhstime = rhs_wtime;

  // Copy variables
  load_vars(NV_DATA_P(uvec));

  // Call rhs function to get extra variables at this time
  real tstart = MPI_Wtime();
  (*func)(simtime);
  rhs_wtime += MPI_Wtime() - tstart;
  rhs_ncalls++;
  
  if(flag < 0) {
    output.write("ERROR CVODE solve failed at t = %e, flag = %d\n", simtime, flag);
    return -1.0;
  }
  
#ifdef CHECK
  msg_stack.pop(msg_point);
#endif

  return simtime;
}

int integrate(struct Integrator* integrator, double tout, double* t)
{
  if (integrator->em->nRates > 0)
  {
    /* need to integrate if we have any differential equations */
    int flag;
    /* Make sure we don't go past the specified end time - could run into
       trouble if we're almost reaching a threshold */
    flag = CVodeSetStopTime(integrator->cvode_mem,(realtype)tout);
    if (check_flag(&flag,"CVode",1)) return(ERR);
    flag = CVode(integrator->cvode_mem,tout,integrator->y,t,CV_NORMAL);
    if (check_flag(&flag,"CVode",1)) return(ERR);
    /* we also need to evaluate all the other variables that are not required
       to be updated during integration */
    integrator->em->evaluateVariables(*t);
  }
  else
  {
	/* no differential equations so just evaluate once */
	integrator->em->computeRates(tout);
	integrator->em->evaluateVariables(tout);
    *t = tout;
  }

  /*
   * Now that using CV_NORMAL_TSTOP this is no longer required?
   */
#ifdef OLD_CODE
  /* only the y array is gonna be at the desired tout, so we need to also
     update the full variables array */
  ud->BOUND[0] = *t;
  ud->methods->ComputeVariables(ud->BOUND,ud->RATES,ud->CONSTANTS,
    ud->VARIABLES);
#endif
  
  /* Make sure the outputs are up-to-date */
  integrator->em->getOutputs(*t);
  return(OK);
}

int CVAdataStore(CVadjMem ca_mem, CkpntMem ck_mem)
{
  CVodeMem cv_mem;
  DtpntMem *dt_mem;
  realtype t;
  long int i;
  int flag;

  cv_mem = ca_mem->cv_mem;
  dt_mem = ca_mem->dt_mem;

  /* Initialize cv_mem with data from ck_mem */
  flag = CVAckpntGet(cv_mem, ck_mem);
  if (flag != CV_SUCCESS) return(CV_REIFWD_FAIL);

  /* Set first structure in dt_mem[0] */
  dt_mem[0]->t = t0_;
  storePnt(cv_mem, dt_mem[0]);

  /* Run CVode to set following structures in dt_mem[i] */

  i = 1;
  do {
    flag = CVode(cv_mem, t1_, ytmp, &t, CV_ONE_STEP);
    if (flag < 0) return(CV_FWD_FAIL);
    dt_mem[i]->t = t;
    storePnt(cv_mem, dt_mem[i]);
    i++;
  } while (t<t1_);

  /* New data is now available */
  ckpntData = ck_mem;
  newData = TRUE;
  np  = i;

  return(CV_SUCCESS);
}

PetscErrorCode TSStep_Sundials(TS ts)
{
  TS_Sundials    *cvode = (TS_Sundials*)ts->data;
  PetscErrorCode ierr;
  PetscInt       flag;
  long int       its,nsteps;
  realtype       t,tout;
  PetscScalar    *y_data;
  void           *mem;

  PetscFunctionBegin;
  mem  = cvode->mem;
  tout = ts->max_time;
  ierr = VecGetArray(ts->vec_sol,&y_data);CHKERRQ(ierr);
  N_VSetArrayPointer((realtype*)y_data,cvode->y);
  ierr = VecRestoreArray(ts->vec_sol,NULL);CHKERRQ(ierr);

  ierr = TSPreStep(ts);CHKERRQ(ierr);

  /* We would like to call TSPreStep() when starting each step (including rejections) and TSPreStage() before each
   * stage solve, but CVode does not appear to support this. */
  if (cvode->monitorstep) flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
  else flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);

  if (flag) { /* display error message */
    switch (flag) {
      case CV_ILL_INPUT:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
        break;
      case CV_TOO_CLOSE:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
        break;
      case CV_TOO_MUCH_WORK: {
        PetscReal      tcur;
        ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr);
        ierr = CVodeGetCurrentTime(mem,&tcur);CHKERRQ(ierr);
        SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_WORK. At t=%G, nsteps %D exceeds mxstep %D. Increase '-ts_max_steps <>' or modify TSSetDuration()",tcur,nsteps,ts->max_steps);
      } break;
      case CV_TOO_MUCH_ACC:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
        break;
      case CV_ERR_FAILURE:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
        break;
      case CV_CONV_FAILURE:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
        break;
      case CV_LINIT_FAIL:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
        break;
      case CV_LSETUP_FAIL:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
        break;
      case CV_LSOLVE_FAIL:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
        break;
      case CV_RHSFUNC_FAIL:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
        break;
      case CV_FIRST_RHSFUNC_ERR:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
        break;
      case CV_REPTD_RHSFUNC_ERR:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
        break;
      case CV_UNREC_RHSFUNC_ERR:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
        break;
      case CV_RTFUNC_FAIL:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
        break;
      default:
        SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
    }
  }

  /* copy the solution from cvode->y to cvode->update and sol */
  ierr = VecPlaceArray(cvode->w1,y_data);CHKERRQ(ierr);
  ierr = VecCopy(cvode->w1,cvode->update);CHKERRQ(ierr);
  ierr = VecResetArray(cvode->w1);CHKERRQ(ierr);
  ierr = VecCopy(cvode->update,ts->vec_sol);CHKERRQ(ierr);
  ierr = CVodeGetNumNonlinSolvIters(mem,&its);CHKERRQ(ierr);
  ierr = CVSpilsGetNumLinIters(mem, &its);
  ts->snes_its = its; ts->ksp_its = its;

  ts->time_step = t - ts->ptime;
  ts->ptime     = t;
  ts->steps++;

  ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr);
  if (!cvode->monitorstep) ts->steps = nsteps;
  PetscFunctionReturn(0);
}

int main()
{
    realtype abstol=ATOL, reltol=RTOL, t, tout;
    N_Vector c;
    WebData wdata;
    void *cvode_mem;
    booleantype firstrun;
    int jpre, gstype, flag;
    int ns, mxns, iout;

    c = NULL;
    wdata = NULL;
    cvode_mem = NULL;

    /* Initializations */
    c = N_VNew_Serial(NEQ);
    if(check_flag((void *)c, "N_VNew_Serial", 0)) return(1);
    wdata = AllocUserData();
    if(check_flag((void *)wdata, "AllocUserData", 2)) return(1);
    InitUserData(wdata);
    ns = wdata->ns;
    mxns = wdata->mxns;

    /* Print problem description */
    PrintIntro();

    /* Loop over jpre and gstype (four cases) */
    for (jpre = PREC_LEFT; jpre <= PREC_RIGHT; jpre++) {
        for (gstype = MODIFIED_GS; gstype <= CLASSICAL_GS; gstype++) {

            /* Initialize c and print heading */
            CInit(c, wdata);
            PrintHeader(jpre, gstype);

            /* Call CVodeInit or CVodeReInit, then CVSpgmr to set up problem */

            firstrun = (jpre == PREC_LEFT) && (gstype == MODIFIED_GS);
            if (firstrun) {
                cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
                if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

                wdata->cvode_mem = cvode_mem;

                flag = CVodeSetUserData(cvode_mem, wdata);
                if(check_flag(&flag, "CVodeSetUserData", 1)) return(1);

                flag = CVodeInit(cvode_mem, f, T0, c);
                if(check_flag(&flag, "CVodeInit", 1)) return(1);

                flag = CVodeSStolerances(cvode_mem, reltol, abstol);
                if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);

                flag = CVSpgmr(cvode_mem, jpre, MAXL);
                if(check_flag(&flag, "CVSpgmr", 1)) return(1);

                flag = CVSpilsSetGSType(cvode_mem, gstype);
                if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);

                flag = CVSpilsSetEpsLin(cvode_mem, DELT);
                if(check_flag(&flag, "CVSpilsSetEpsLin", 1)) return(1);

                flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve);
                if(check_flag(&flag, "CVSpilsSetPreconditioner", 1)) return(1);

            } else {

                flag = CVodeReInit(cvode_mem, T0, c);
                if(check_flag(&flag, "CVodeReInit", 1)) return(1);

                flag = CVSpilsSetPrecType(cvode_mem, jpre);
                check_flag(&flag, "CVSpilsSetPrecType", 1);
                flag = CVSpilsSetGSType(cvode_mem, gstype);
                if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);

            }

            /* Print initial values */
            if (firstrun) PrintAllSpecies(c, ns, mxns, T0);

            /* Loop over output points, call CVode, print sample solution values. */
            tout = T1;
            for (iout = 1; iout <= NOUT; iout++) {
                flag = CVode(cvode_mem, tout, c, &t, CV_NORMAL);
                PrintOutput(cvode_mem, t);
                if (firstrun && (iout % 3 == 0)) PrintAllSpecies(c, ns, mxns, t);
                if(check_flag(&flag, "CVode", 1)) break;
                if (tout > RCONST(0.9)) tout += DTOUT;
                else tout *= TOUT_MULT;
            }

            /* Print final statistics, and loop for next case */
            PrintFinalStats(cvode_mem);

        }
    }

    /* Free all memory */
    CVodeFree(&cvode_mem);
    N_VDestroy_Serial(c);
    FreeUserData(wdata);
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  void *cvode_mem;
  UserData data;
  realtype t, tout;
  N_Vector y;
  int iout, flag, nthreads, nnz;

  realtype pbar[NS];
  int is; 
  N_Vector *yS;
  booleantype sensi, err_con;
  int sensi_meth;

  cvode_mem = NULL;
  data      = NULL;
  y         =  NULL;
  yS        = NULL;

  /* Process arguments */
  ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);

  /* User data structure */
  data = (UserData) malloc(sizeof *data);
  if (check_flag((void *)data, "malloc", 2)) return(1);
  data->p[0] = RCONST(0.04);
  data->p[1] = RCONST(1.0e4);
  data->p[2] = RCONST(3.0e7);

  /* Initial conditions */
  y = N_VNew_Serial(NEQ);
  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);

  Ith(y,1) = Y1;
  Ith(y,2) = Y2;
  Ith(y,3) = Y3;

  /* Call CVodeCreate to create the solver memory and specify the 
     Backward Differentiation Formula and the use of a Newton iteration */
  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

  /* Call CVodeInit to initialize the integrator memory and specify the
     user's right hand side function in y'=f(t,y), the initial time T0, and
     the initial dependent variable vector y. */
  flag = CVodeInit(cvode_mem, f, T0, y);
  if (check_flag(&flag, "CVodeInit", 1)) return(1);

  /* Call CVodeWFtolerances to specify a user-supplied function ewt that sets
     the multiplicative error weights W_i for use in the weighted RMS norm */
  flag = CVodeWFtolerances(cvode_mem, ewt);
  if (check_flag(&flag, "CVodeSetEwtFn", 1)) return(1);

  /* Attach user data */
  flag = CVodeSetUserData(cvode_mem, data);
  if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);

  /* Call CVKLU to specify the CVKLU sparse direct linear solver */
  nthreads = 1;                 /* no. of threads to use when factoring the system*/
  nnz = NEQ * NEQ;              /* max no. of nonzeros entries in the Jac */
  flag = CVSuperLUMT(cvode_mem, nthreads, NEQ, nnz);
  if (check_flag(&flag, "CVSuperLUMT", 1)) return(1);

  /* Set the Jacobian routine to Jac (user-supplied) */
  flag = CVSlsSetSparseJacFn(cvode_mem, Jac);
  if (check_flag(&flag, "CVSlsSetSparseJacFn", 1)) return(1);

  printf("\n3-species chemical kinetics problem\n");

  /* Sensitivity-related settings */
  if (sensi) {

    /* Set parameter scaling factor */
    pbar[0] = data->p[0];
    pbar[1] = data->p[1];
    pbar[2] = data->p[2];

    /* Set sensitivity initial conditions */
    yS = N_VCloneVectorArray_Serial(NS, y);
    if (check_flag((void *)yS, "N_VCloneVectorArray_Serial", 0)) return(1);
    for (is=0;is<NS;is++) N_VConst(ZERO, yS[is]);

    /* Call CVodeSensInit1 to activate forward sensitivity computations
       and allocate internal memory for COVEDS related to sensitivity
       calculations. Computes the right-hand sides of the sensitivity
       ODE, one at a time */
    flag = CVodeSensInit1(cvode_mem, NS, sensi_meth, fS, yS);
    if(check_flag(&flag, "CVodeSensInit", 1)) return(1);

    /* Call CVodeSensEEtolerances to estimate tolerances for sensitivity 
       variables based on the rolerances supplied for states variables and 
       the scaling factor pbar */
    flag = CVodeSensEEtolerances(cvode_mem);
    if(check_flag(&flag, "CVodeSensEEtolerances", 1)) return(1);

    /* Set sensitivity analysis optional inputs */
    /* Call CVodeSetSensErrCon to specify the error control strategy for 
       sensitivity variables */
    flag = CVodeSetSensErrCon(cvode_mem, err_con);
    if (check_flag(&flag, "CVodeSetSensErrCon", 1)) return(1);
//.........这里部分代码省略.........

int run_rate_state_sim(std::vector<std::vector<realtype> > &results, RSParams &params) {
	realtype		long_term_reltol, event_reltol, t, tout, tbase=0;
	N_Vector		y, long_term_abstol, event_abstol;
	unsigned int	i, n;
	int				flag, err_code;
	void			*long_term_cvode, *event_cvode, *current_cvode;
	
	// Create serial vector of length NEQ for I.C. and abstol
	y = N_VNew_Serial(params.num_eqs()*params.num_blocks());
	if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
	long_term_abstol = N_VNew_Serial(params.num_eqs()*params.num_blocks()); 
	if (check_flag((void *)long_term_abstol, "N_VNew_Serial", 0)) return(1);
	event_abstol = N_VNew_Serial(params.num_eqs()*params.num_blocks()); 
	if (check_flag((void *)event_abstol, "N_VNew_Serial", 0)) return(1);
	
	// Initialize y
	for (i=0;i<params.num_blocks();++i) {
		NV_Ith_S(y,i*params.num_eqs()+EQ_X) = params.init_val(i, EQ_X);
		NV_Ith_S(y,i*params.num_eqs()+EQ_V) = params.init_val(i, EQ_V);
		NV_Ith_S(y,i*params.num_eqs()+EQ_H) = params.init_val(i, EQ_H);
	}
	
	/* Initialize interactions */
	/*interaction = new realtype[NBLOCKS*NBLOCKS];
	double int_level = 1e-2;
	double dropoff = 1.1;
	for (i=0;i<NBLOCKS;++i) {
		for (n=0;n<NBLOCKS;++n) {
			interaction[i*NBLOCKS+n] = (i==n?(1.0-int_level):int_level);
		}
	}*/
	
	/* Set the scalar relative tolerance */
	long_term_reltol = RCONST(1.0e-12);
	event_reltol = RCONST(1.0e-12);
	/* Set the vector absolute tolerance */
	for (i=0;i<params.num_blocks();++i) {
		Xth(long_term_abstol,i) = RCONST(1.0e-12);
		Vth(long_term_abstol,i) = RCONST(1.0e-12);
		Hth(long_term_abstol,i) = RCONST(1.0e-12);
		Xth(event_abstol,i) = RCONST(1.0e-12);
		Vth(event_abstol,i) = RCONST(1.0e-12);
		Hth(event_abstol,i) = RCONST(1.0e-12);
	}
	
	/* Call CVodeCreate to create the solver memory and specify the 
	 * Backward Differentiation Formula and the use of a Newton iteration */
	long_term_cvode = CVodeCreate(CV_BDF, CV_NEWTON);
	if (check_flag((void *)long_term_cvode, "CVodeCreate", 0)) return(1);
	event_cvode = CVodeCreate(CV_BDF, CV_NEWTON);
	if (check_flag((void *)event_cvode, "CVodeCreate", 0)) return(1);
	
	// Turn off error messages
	//CVodeSetErrFile(long_term_cvode, NULL);
	//CVodeSetErrFile(event_cvode, NULL);
	
	/* Call CVodeInit to initialize the integrator memory and specify the
	 * user's right hand side function in y'=f(t,y), the inital time T0, and
	 * the initial dependent variable vector y. */
	flag = CVodeInit(long_term_cvode, func, T0, y);
	if (check_flag(&flag, "CVodeInit", 1)) return(1);
	flag = CVodeInit(event_cvode, func, T0, y);
	if (check_flag(&flag, "CVodeInit", 1)) return(1);
	
	/* Call CVodeSVtolerances to specify the scalar relative tolerance
	 * and vector absolute tolerances */
	flag = CVodeSVtolerances(long_term_cvode, long_term_reltol, long_term_abstol);
	if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
	flag = CVodeSVtolerances(event_cvode, event_reltol, event_abstol);
	if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
	
	/* Set the root finding function */
	//flag = CVodeRootInit(long_term_cvode, params.num_blocks(), vel_switch_finder);
	//flag = CVodeRootInit(event_cvode, params.num_blocks(), vel_switch_finder);
	//if (check_flag(&flag, "CVodeRootInit", 1)) return(1);
	
	/* Call CVDense to specify the CVDENSE dense linear solver */
	//flag = CVSpbcg(cvode_mem, PREC_NONE, 0);
	//if (check_flag(&flag, "CVSpbcg", 1)) return(1);
	//flag = CVSpgmr(cvode_mem, PREC_NONE, 0);
	//if (check_flag(&flag, "CVSpgmr", 1)) return(1);
	flag = CVDense(long_term_cvode, params.num_eqs()*params.num_blocks());
	if (check_flag(&flag, "CVDense", 1)) return(1);
	flag = CVDense(event_cvode, params.num_eqs()*params.num_blocks());
	if (check_flag(&flag, "CVDense", 1)) return(1);
	
	flag = CVodeSetUserData(long_term_cvode, &params);
	if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);
	flag = CVodeSetUserData(event_cvode, &params);
	if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);
	
	CVodeSetMaxNumSteps(long_term_cvode, 100000);
	CVodeSetMaxNumSteps(event_cvode, 100000);
	
	/* Set the Jacobian routine to Jac (user-supplied) */
	flag = CVDlsSetDenseJacFn(long_term_cvode, Jac);
	if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);
	flag = CVDlsSetDenseJacFn(event_cvode, Jac);
	if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);
	
//.........这里部分代码省略.........

int main(int narg, char **args)
{
    realtype reltol, t, tout;
    N_Vector state, abstol;
    void *cvode_mem;
    int flag, flagr;
    int rootsfound[NRF];
    int rootdir[] = {1,};

    FILE *pout;
    if(!(pout = fopen("results/iaf_v.dat", "w"))){
        fprintf(stderr, "Cannot open file results/iaf_v.dat. Are you trying to write to a non-existent directory? Exiting...\n");
        exit(1);
    }

    state = abstol = NULL;
    cvode_mem = NULL;

    state = N_VNew_Serial(NEQ);
    if (check_flag((void *)state, "N_VNew_Serial", 0)) return(1);
    abstol = N_VNew_Serial(NEQ); 
    if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);
    
    realtype reset = -0.07;
    realtype C = 3.2e-12;
    realtype thresh = -0.055;
    realtype gleak = 2e-10;
    realtype eleak = -0.053;
    realtype p[] = {reset, C, thresh, gleak, eleak, };

    realtype v = reset;
    NV_Ith_S(state, 0) = reset;

    reltol = RTOL;
    NV_Ith_S(abstol,0) = ATOL0;
 
    /* Allocations and initializations */
    cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
    if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
    
    flag = CVodeInit(cvode_mem, dstate_dt, T0, state);
    if (check_flag(&flag, "CVodeInit", 1)) return(1);
   
    flag = CVodeSetUserData(cvode_mem, p);
    if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);
   
    flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
    if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
   
        flag = CVodeRootInit(cvode_mem, NRF, root_functions);
    if (check_flag(&flag, "CVodeRootInit", 1)) return(1);

    CVodeSetRootDirection(cvode_mem, rootdir);
    if (check_flag(&flag, "CVodeSetRootDirection", 1)) return(1);
        
   
    flag = CVDense(cvode_mem, NEQ);
    if (check_flag(&flag, "CVDense", 1)) return(1);


    printf(" \n Integrating iaf \n\n");
    printf("#t v, \n");
    PrintOutput(pout, t, state);
   
    tout = DT;
    while(1) {
        flag = CVode(cvode_mem, tout, state, &t, CV_NORMAL);
        
        if(flag == CV_ROOT_RETURN) {
            /* Event detected */
            flagr = CVodeGetRootInfo(cvode_mem, rootsfound);
            if (check_flag(&flagr, "CVodeGetRootInfo", 1)) return(1);
            PrintRootInfo(t, state, rootsfound);
          
            if(rootsfound[0]){
                //condition_0
                v = NV_Ith_S(state, 0);
                NV_Ith_S(state, 0) = reset;

        }
            

        /* Restart integration with event-corrected state */
            flag = CVodeSetUserData(cvode_mem, p);
            if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);
        CVodeReInit(cvode_mem, t, state);
        //PrintRootInfo(t, state, rootsfound);
    }
        else
                {
            PrintOutput(pout, t, state);
            if(check_flag(&flag, "CVode", 1)) break;
            if(flag == CV_SUCCESS) {
                tout += DT;
            }
            if (t >= T1) break;
        }

    }

//.........这里部分代码省略.........