CVODE的cvDiurnal_kry例子学习

问题描述

This program solves the semi-discretized form of a pair of kinetics-advection-diffusion partial differential equations, which represent a simplified model for the transport, production,
and loss of ozone and the oxygen singlet in the upper atmosphere. The problem includes nonlinear diurnal kinetics, horizontal advection and diffusion, and nonuniform vertical diffusion.

cvDiurnal_kry问题研究的是，两种成分的对流-扩散运动方程的微分方程的半离散化形式的求解。
比如大气层中，氧气，臭氧的转换，产生，损耗，可以用上面的形式来建立一个简单的模型。
该问题包含了:

1. diurnal kinetics,
2. 水平的对流和扩散
3. 非均匀的垂直方向的扩散

PDE方程

1. c1,c2分别是两种成分的微分方程
2. 变量：水平方向的坐标x,垂直方向的坐标y，以及时间t,
3. K_h:水平扩散系数，常量=4.0e-6
4. V：对流系数，常量=0.001
5. K_v:表达式=10−8*exp(y/5)
6. R1,R2是反应项

x代表了空间的横坐标
y代表了空间的纵坐标
t=0的时候，气体成分在空间各个位置的分布用(6)描述

我们将空间划分成10x10的格子，在每个格子中，对问题进行求解，因为成分有2种，一共有10x10个
切片，因此该问题的总变量数就是100x2=200.

---

代码实现

u:计算每个网格中,t=0时，(6)的c1,c2的值

u = N_VNew_Serial(NEQ);
 SetInitialProfiles(u, data->dx, data->dy);
static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy)
{
  int jx, jy;
  realtype x, y, cx, cy;
  realtype *udata;

  /* Set pointer to data array in vector u. */

  udata = N_VGetArrayPointer(u);

  /* Load initial profiles of c1 and c2 into u vector */
  //计算公式(6)时，取格子里起始点+变化率作为本格子本格子的起始点
  for (jy=0; jy < MY; jy++) {
    y = YMIN + jy*dy;
    cy = SQR(RCONST(0.1)*(y - YMID));
    cy = ONE - cy + RCONST(0.5)*SQR(cy);
    for (jx=0; jx < MX; jx++) {
      x = XMIN + jx*dx;
      cx = SQR(RCONST(0.1)*(x - XMID));
      cx = ONE - cx + RCONST(0.5)*SQR(cx);
      IJKth(udata,1,jx,jy) = C1_SCALE*cx*cy;
      IJKth(udata,2,jx,jy) = C2_SCALE*cx*cy;
    }
  }
}

存储数据的结构体：

typedef struct {
  realtype **P[MX][MY];
  realtype **Jbd[MX][MY];
  sunindextype *pivot[MX][MY];
 
  realtype q4;
  realtype om;
  realtype  dx;//(XMAX-XMIN)/(MX-1)，x方向的变化率，常量
  realtype  dy;// (YMAX-YMIN)/(MY-1)，y方向的变化率，常量
  realtype  hdco;//KH/SQR(data->dx),水平方向上的扩散系数，公式（4）等式右边的第一项
  realtype  haco;//VEL/(TWO*data->dx),,水平方向上的对流系数，公式（4）等式右边的第二项
  realtype  vdco;//(ONE/SQR(data->dy))*KV0,垂直方向上的扩散系数，公式（4）等式右边的第三项
} *UserData;

微分方程的实现

/* f routine. Compute RHS function f(t,u). */

static int f(realtype t, N_Vector u, N_Vector udot, void *user_data)
{
  realtype q3, c1, c2, c1dn, c2dn, c1up, c2up, c1lt, c2lt;
  realtype c1rt, c2rt, cydn, cyup, hord1, hord2, horad1, horad2;
  realtype qq1, qq2, qq3, qq4, rkin1, rkin2, s, vertd1, vertd2, ydn, yup;
  realtype q4coef, dely, verdco, hordco, horaco;
  realtype *udata, *dudata;
  int jx, jy, idn, iup, ileft, iright;
  UserData data;

  data   = (UserData) user_data;
  udata  = N_VGetArrayPointer(u);
  dudata = N_VGetArrayPointer(udot);

  /* Set diurnal rate coefficients. */
// 先计算sin(wt),根据值来求q3,q4
  s = sin(data->om*t);
  if (s > ZERO) {
    q3 = exp(-A3/s);
    data->q4 = exp(-A4/s);
  } else {
    q3 = ZERO;
    data->q4 = ZERO;
  }

  /* Make local copies of problem variables, for efficiency. */

  q4coef = data->q4;
  dely = data->dy;
  verdco = data->vdco;
  hordco  = data->hdco;
  horaco  = data->haco;

  /* Loop over all grid points. */
//对每个格子进行计算，MY=10，表示垂直方向的格子划分
  for (jy=0; jy < MY; jy++) {

    /* Set vertical diffusion coefficients at jy +- 1/2 */
    //垂直方向上，每个格子的范围是[jy-1/2,jy+1/2]
    ydn = YMIN + (jy - RCONST(0.5))*dely;
    yup = ydn + dely;
    
    cydn = verdco*exp(RCONST(0.2)*ydn);
    cyup = verdco*exp(RCONST(0.2)*yup);
    
    //idn：当前格子的前一个格子的偏移
    //iup:当前格式的下一个格子的偏移
    idn = (jy == 0) ? 1 : -1;
    iup = (jy == MY-1) ? -1 : 1;
    for (jx=0; jx < MX; jx++) {

      /* Extract c1 and c2, and set kinetic rate terms. */

	  //取出对应格子的空间起始点x,y
      c1 = IJKth(udata,1,jx,jy);
      c2 = IJKth(udata,2,jx,jy);

	  //计算R1,R2
      qq1 = Q1*c1*C3;//R1，R2的第一项
      qq2 = Q2*c1*c2;//R1,R2的第二项
      qq3 = q3*C3;//R1的第三项
      qq4 = q4coef*c2;//R1的第四项，也是R2的第三项
      rkin1 = -qq1 - qq2 + TWO*qq3 + qq4;//R1
      rkin2 = qq1 - qq2 - qq4;//R2

      /* Set vertical diffusion terms. */
       //取出c1,c2在前一个格子，后一个格子的值
      c1dn = IJKth(udata,1,jx,jy+idn);
      c2dn = IJKth(udata,2,jx,jy+idn);
      c1up = IJKth(udata,1,jx,jy+iup);
      c2up = IJKth(udata,2,jx,jy+iup);
      vertd1 = cyup*(c1up - c1) - cydn*(c1 - c1dn);
      vertd2 = cyup*(c2up - c2) - cydn*(c2 - c2dn);

      /* Set horizontal diffusion and advection terms. */

      ileft = (jx == 0) ? 1 : -1;
      iright =(jx == MX-1) ? -1 : 1;
      c1lt = IJKth(udata,1,jx+ileft,jy);
      c2lt = IJKth(udata,2,jx+ileft,jy);
      c1rt = IJKth(udata,1,jx+iright,jy);
      c2rt = IJKth(udata,2,jx+iright,jy);
      hord1 = hordco*(c1rt - TWO*c1 + c1lt);
      hord2 = hordco*(c2rt - TWO*c2 + c2lt);
      horad1 = horaco*(c1rt - c1lt);
      horad2 = horaco*(c2rt - c2lt);

      /* Load all terms into udot. */

      IJKth(dudata, 1, jx, jy) = vertd1 + hord1 + horad1 + rkin1;
      IJKth(dudata, 2, jx, jy) = vertd2 + hord2 + horad2 + rkin2;
    }
  }

  return(0);
}