问题描述

问题描述

This the semi- form of a pair of -- , which a model for the , ,

and loss of ozone and the in the upper . The , and , and .

问题研究的是，两种成分的对流-扩散运动方程的微分方程的半离散化形式的求解。

比如大气层中，氧气，臭氧的转换，产生，损耗，可以用上面的形式来建立一个简单的模型。

该问题包含了:

,水平的对流和扩散非均匀的垂直方向的扩散 PDE方程

c1,c2分别是两种成分的微分方程变量：水平方向的坐标x,垂直方向的坐标y，以及时间t,Kh:水平扩散系数，常量=4.0e-6V：对流系数，常量=0.001Kv:表达式=10−8*exp(y/5)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,q4s = 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,yc1 = IJKth(udata,1,jx,jy);c2 = IJKth(udata,2,jx,jy);//计算R1,R2qq1 = 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;//R1rkin2 = 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);
}