C++ N_VNew_Serial函数代码示例

void RefSol(realtype tout, N_Vector yref)
{
  void *cpode_mem;
  N_Vector yy, yp;
  realtype tol, t, th, thd;
  int flag;
  

  yy = N_VNew_Serial(2);
  yp = N_VNew_Serial(2);
  Ith(yy,1) = 0.0;  /* theta */
  Ith(yy,2) = 0.0;  /* thetad */
  tol = TOL_REF;

  cpode_mem = CPodeCreate(CP_EXPL, CP_BDF, CP_NEWTON);  
  flag = CPodeSetMaxNumSteps(cpode_mem, 100000);
  flag = CPodeInit(cpode_mem, (void *)fref, NULL, 0.0, yy, yp, CP_SS, tol, &tol);
  flag = CPDense(cpode_mem, 2);

  flag = CPodeSetStopTime(cpode_mem, tout);
  flag = CPode(cpode_mem, tout, &t, yy, yp, CP_NORMAL_TSTOP);
  th  = Ith(yy,1);
  thd = Ith(yy,2);
  Ith(yref,1) = cos(th);
  Ith(yref,2) = sin(th);
  Ith(yref,3) = -thd*sin(th);
  Ith(yref,4) =  thd*cos(th);
  
  N_VDestroy_Serial(yy);
  N_VDestroy_Serial(yp);
  CPodeFree(&cpode_mem);

  return;
}

void GetSol(void *cpode_mem, N_Vector yy0, realtype tol, 
            realtype tout, booleantype proj, N_Vector yref)
{
  N_Vector yy, yp;
  realtype t, x, y, xd, yd, g;
  int flag;
  long int nst, nfe, nsetups, nje, nfeLS, ncfn, netf;

  if (proj) {
    printf(" YES   ");
    CPodeSetProjFrequency(cpode_mem, 1);
  } else {
    CPodeSetProjFrequency(cpode_mem, 0);
    printf(" NO    ");
  }

  yy = N_VNew_Serial(4);
  yp = N_VNew_Serial(4);

  flag = CPodeReInit(cpode_mem, (void *)f, NULL, 0.0, yy0, NULL, CP_SS, tol, &tol);

  flag = CPode(cpode_mem, tout, &t, yy, yp, CP_NORMAL_TSTOP);

  x  = Ith(yy,1);
  y  = Ith(yy,2);
  g = ABS(x*x + y*y - 1.0);

  N_VLinearSum(1.0, yy, -1.0, yref, yy);

  N_VAbs(yy, yy);

  x  = Ith(yy,1);
  y  = Ith(yy,2);  
  xd = Ith(yy,3);
  yd = Ith(yy,4);


  printf("%9.2e  %9.2e  %9.2e  %9.2e  |  %9.2e  |",  
         Ith(yy,1),Ith(yy,2),Ith(yy,3),Ith(yy,4),g);

  CPodeGetNumSteps(cpode_mem, &nst);
  CPodeGetNumFctEvals(cpode_mem, &nfe);
  CPodeGetNumLinSolvSetups(cpode_mem, &nsetups);
  CPodeGetNumErrTestFails(cpode_mem, &netf);
  CPodeGetNumNonlinSolvConvFails(cpode_mem, &ncfn);

  CPDlsGetNumJacEvals(cpode_mem, &nje);
  CPDlsGetNumFctEvals(cpode_mem, &nfeLS);


  printf(" %6ld   %6ld+%-4ld  %4ld (%3ld)  |  %3ld  %3ld\n",
         nst, nfe, nfeLS, nsetups, nje, ncfn, netf);

  N_VDestroy_Serial(yy);
  N_VDestroy_Serial(yp);

  return;
}

int main(void)
{
  realtype dx, dy, reltol, abstol, t, tout, umax;
  N_Vector u, up;
  UserData data;
  void *cvode_mem;
  int iout, flag;
  long int nst;

  u = NULL;
  data = NULL;
  cvode_mem = NULL;

  u  = N_VNew_Serial(NEQ);
  up = N_VNew_Serial(NEQ);

  reltol = ZERO;
  abstol = ATOL;

  data = (UserData) malloc(sizeof *data);
  dx = data->dx = XMAX/(MX+1);
  dy = data->dy = YMAX/(MY+1);
  data->hdcoef = ONE/(dx*dx);
  data->hacoef = HALF/(TWO*dx);
  data->vdcoef = ONE/(dy*dy);

  SetIC(u, data);

  cvode_mem = CPodeCreate(CP_EXPL, CP_BDF, CP_NEWTON);
  flag = CPodeInit(cvode_mem, (void *)f, data, T0, u, NULL, CP_SS, reltol, &abstol);

  flag = CPLapackBand(cvode_mem, NEQ, MY, MY);
  flag = CPDlsSetJacFn(cvode_mem, (void *)Jac, data);

  /* In loop over output points: call CPode, print results, test for errors */
  umax = N_VMaxNorm(u);
  PrintHeader(reltol, abstol, umax);
  for(iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {
    flag = CPode(cvode_mem, tout, &t, u, up, CP_NORMAL);
    umax = N_VMaxNorm(u);
    flag = CPodeGetNumSteps(cvode_mem, &nst);
    PrintOutput(t, umax, nst);
  }

  PrintFinalStats(cvode_mem);

  N_VDestroy_Serial(u);
  CPodeFree(&cvode_mem);
  free(data);

  return(0);
}

	void OpenSMOKE_CVODE_Sundials<T>::MemoryAllocation(const int n)
	{	
		this->n_		=	n;					// Number of equations

		this->y0_ 	= new double[this->n_];
		this->y_ 	= new double[this->n_];
 
		y0Sundials_ = N_VNew_Serial(this->n_);
		if (check_flag((void *)y0Sundials_, std::string("N_VNew_Serial"), 0))	exit(-1);

		ySundials_ = N_VNew_Serial(this->n_);
		if (check_flag((void *)ySundials_, std::string("N_VNew_Serial"), 0))	exit(-1);
	}

  void SundialsInterface::init_memory(void* mem) const {
    Integrator::init_memory(mem);
    auto m = static_cast<SundialsMemory*>(mem);

    // Allocate n-vectors
    m->xz = N_VNew_Serial(nx_+nz_);
    m->q = N_VNew_Serial(nq_);
    m->rxz = N_VNew_Serial(nrx_+nrz_);
    m->rq = N_VNew_Serial(nrq_);

    // Allocate memory
    m->p.resize(np_);
    m->rp.resize(nrp_);
  }

static WebData AllocUserData(void)
{
  int i, ngrp = NGRP, ns = NS;
  WebData wdata;

  wdata = (WebData) malloc(sizeof *wdata);
  for(i=0; i < ngrp; i++) {
    (wdata->P)[i] = newDenseMat(ns, ns);
    (wdata->pivot)[i] = newIndexArray(ns);
  }
  wdata->rewt  = N_VNew_Serial(NEQ+1);
  wdata->vtemp = N_VNew_Serial(NEQ+1);

  return(wdata);
}

void CVodesIntegrator::sensInit(double t0, FuncEval& func)
{
    m_np = func.nparams();
    m_sens_ok = false;

    N_Vector y = N_VNew_Serial(static_cast<sd_size_t>(func.neq()));
    m_yS = N_VCloneVectorArray_Serial(static_cast<sd_size_t>(m_np), y);
    for (size_t n = 0; n < m_np; n++) {
        N_VConst(0.0, m_yS[n]);
    }
    N_VDestroy_Serial(y);

    int flag = CVodeSensInit(m_cvode_mem, static_cast<sd_size_t>(m_np),
                             CV_STAGGERED, CVSensRhsFn(0), m_yS);

    if (flag != CV_SUCCESS) {
        throw CanteraError("CVodesIntegrator::sensInit", "Error in CVodeSensInit");
    }
    vector_fp atol(m_np);
    for (size_t n = 0; n < m_np; n++) {
        // This scaling factor is tuned so that reaction and species enthalpy
        // sensitivities can be computed simultaneously with the same abstol.
        atol[n] = m_abstolsens / func.m_paramScales[n];
    }
    flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data());
}

N_Vector Cvode::nvnew(long int n) {
#if PARANEURON
	if (use_partrans_) {
		if (net_cvode_instance->use_long_double_) {
			return N_VNew_NrnParallelLD(0, n, global_neq_);
		}else{
			return N_VNew_Parallel(0, n, global_neq_);
		}
	}
#endif
	if (nctd_ > 1) {
		assert(n == neq_);
		if (!nthsizes_) {
			nthsizes_ = new long int[nrn_nthread];
			for (int i = 0; i < nrn_nthread; ++i) {
				nthsizes_[i] = ctd_[i].nvsize_;
			}
		}
#if 1
		int sum = 0;
		for (int i=0; i < nctd_; ++i) { sum += nthsizes_[i];}
		assert(sum == neq_);
#endif
		if (net_cvode_instance->use_long_double_) {
			return N_VNew_NrnThreadLD(n, nctd_, nthsizes_);
		}else{
			return N_VNew_NrnThread(n, nctd_, nthsizes_);
		}
	}
	if (net_cvode_instance->use_long_double_) {
		return N_VNew_NrnSerialLD(n);
	}else{
		return N_VNew_Serial(n);
	}
}

void CVodesIntegrator::sensInit(double t0, FuncEval& func)
{
    m_np = func.nparams();
    size_t nv = func.neq();
    m_sens_ok = false;

    doublereal* data;
    N_Vector y;
    y = N_VNew_Serial(static_cast<sd_size_t>(nv));
    m_yS = N_VCloneVectorArray_Serial(static_cast<sd_size_t>(m_np), y);
    for (size_t n = 0; n < m_np; n++) {
        data = NV_DATA_S(m_yS[n]);
        for (size_t j = 0; j < nv; j++) {
            data[j] =0.0;
        }
    }

    int flag = CVodeSensInit(m_cvode_mem, static_cast<sd_size_t>(m_np),
                             CV_STAGGERED, CVSensRhsFn(0), m_yS);

    if (flag != CV_SUCCESS) {
        throw CVodesErr("Error in CVodeSensMalloc");
    }
    vector_fp atol(m_np, m_abstolsens);
    double rtol = m_reltolsens;
    flag = CVodeSensSStolerances(m_cvode_mem, rtol, atol.data());

}

  void SundialsInterface::init_memory(void* mem) const {
    Integrator::init_memory(mem);
    auto m = static_cast<SundialsMemory*>(mem);

    // Allocate n-vectors
    m->xz = N_VNew_Serial(nx_+nz_);
    m->q = N_VNew_Serial(nq_);
    m->rxz = N_VNew_Serial(nrx_+nrz_);
    m->rq = N_VNew_Serial(nrq_);

    // Reset linear solvers
    linsolF_.reset(get_function("jacF").sparsity_out(0));
    if (nrx_>0) {
      linsolB_.reset(get_function("jacB").sparsity_out(0));
    }
  }

int ViCaRS::init(SimEquations *eqns) {
	unsigned int	num_local, i;
	int				flag;
	
	_eqns = eqns;
	_num_equations = eqns->num_equations();
	
	//flag = fill_greens_matrix();
	//if (flag) return flag;
	
	num_local = (unsigned int)_bdata.size();
	_vars = N_VNew_Serial(_num_global_blocks*_num_equations);
	if (_vars == NULL) return 1;
	
	// Init CVodes structures
	flag = _eqns->init(this);
	if (flag) return flag;
    
	for (i=0;i<_solvers.size();++i) _solvers[i]->init_solver(this);
	
	_cur_solver = 0;
	_cur_time = 0;
	
	return 0;
}

static void Problem2(void)
{
  void *fct;
  void *cpode_mem;
  N_Vector y, yp;
  realtype reltol=RTOL, abstol=ATOL, t, tout, erm, hu;
  int flag, iout, qu;

  y = NULL;
  yp = NULL;
  cpode_mem = NULL;

  y = N_VNew_Serial(P2_NEQ);
  N_VConst(ZERO, y);
  NV_Ith_S(y,0) = ONE;

  yp = N_VNew_Serial(P2_NEQ);

  if (ODE == CP_EXPL) {
    fct = (void *)f2;
  } else {
    fct = (void *)res2;
    f2(P2_T0, y, yp, NULL);
  }

  cpode_mem = CPodeCreate(ODE, CP_ADAMS, CP_FUNCTIONAL);
  /*  flag = CPodeSetInitStep(cpode_mem, 2.0e-9);*/
  flag = CPodeInit(cpode_mem, fct, NULL, P2_T0, y, yp, CP_SS, reltol, &abstol);

  printf("\n      t        max.err      qu     hu \n");
  for(iout=1, tout=P2_T1; iout <= P2_NOUT; iout++, tout*=P2_TOUT_MULT) {
    flag = CPode(cpode_mem, tout, &t, y, yp, CP_NORMAL);
    if (flag != CP_SUCCESS) break;
    erm = MaxError(y, t);
    flag = CPodeGetLastOrder(cpode_mem, &qu);
    flag = CPodeGetLastStep(cpode_mem, &hu);
    printf("%10.3f  %12.4le   %2d   %12.4le\n", t, erm, qu, hu);
  }

  PrintFinalStats(cpode_mem);

  CPodeFree(&cpode_mem);
  N_VDestroy_Serial(y);
  N_VDestroy_Serial(yp);

  return;
}

int main()
{
  void *fct;
  void *cpode_mem;
  N_Vector y, yp;
  realtype reltol=RTOL, abstol=ATOL, t, tout, hu;
  int flag, iout, qu;

  y = NULL;
  yp = NULL;
  cpode_mem = NULL;

  y = N_VNew_Serial(P1_NEQ);
  NV_Ith_S(y,0) = TWO;
  NV_Ith_S(y,1) = ZERO;

  yp = N_VNew_Serial(P1_NEQ);
  
  if (ODE == CP_EXPL) {
    fct = (void *)f;
  } else {
    fct = (void *)res;
    f(P1_T0, y, yp, NULL);
  }

  cpode_mem = CPodeCreate(ODE, CP_ADAMS, CP_FUNCTIONAL);
  /*  flag = CPodeSetInitStep(cpode_mem, 4.0e-9);*/
  flag = CPodeInit(cpode_mem, fct, NULL, P1_T0, y, yp, CP_SS, reltol, &abstol);

  printf("\n     t           x              xdot         qu     hu \n");
  for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {
    flag = CPode(cpode_mem, tout, &t, y, yp, CP_NORMAL);
    if (flag != CP_SUCCESS)  break;
    flag = CPodeGetLastOrder(cpode_mem, &qu);
    flag = CPodeGetLastStep(cpode_mem, &hu);    
    printf("%10.5f    %12.5le   %12.5le   %2d    %6.4le\n", t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);
  }
  
  PrintFinalStats(cpode_mem);

  CPodeFree(&cpode_mem);
  N_VDestroy_Serial(y);
  N_VDestroy_Serial(yp);

  return 0;
}

int main()
{
  void *cpode_mem;
  N_Vector yref, yy0;
  realtype tol, tout;
  int i, flag;


  tout = 30.0;
 
  /* Get reference solution */
  yref = N_VNew_Serial(4);
  RefSol(tout, yref);

  /* Initialize solver */
  tol = TOL;
  cpode_mem = CPodeCreate(CP_EXPL, CP_BDF, CP_NEWTON);  
  yy0 = N_VNew_Serial(4);
  Ith(yy0,1) = 1.0;  /* x */
  Ith(yy0,2) = 0.0;  /* y */
  Ith(yy0,3) = 0.0;  /* xd */
  Ith(yy0,4) = 0.0;  /* yd */
  flag = CPodeInit(cpode_mem, (void *)f, NULL, 0.0, yy0, NULL, CP_SS, tol, &tol);
  flag = CPodeSetMaxNumSteps(cpode_mem, 50000);
  flag = CPodeSetStopTime(cpode_mem, tout);
  flag = CPodeProjDefine(cpode_mem, proj, NULL);
  flag = CPDense(cpode_mem, 4);

  for (i=0;i<5;i++) {

    printf("\n\n%.2e\n", tol);
    GetSol(cpode_mem, yy0, tol, tout, TRUE, yref);
    GetSol(cpode_mem, yy0, tol, tout, FALSE, yref);
    tol /= 10.0;
  }

  N_VDestroy_Serial(yref);
  CPodeFree(&cpode_mem);

  return(0);
}

 void CVodesIntegrator::setTolerances(double reltol, int n, double* abstol) {
   m_itol = CV_SV;
   m_nabs = n;
   if (n != m_neq) {
     if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
     m_abstol = reinterpret_cast<void*>(N_VNew_Serial(n));
   }
   for (int i=0; i<n; i++) {
     NV_Ith_S(nv(m_abstol), i) = abstol[i];
   }
   m_reltol = reltol; 
 }