问题描述
Thisthe semi- form of a pair of --, whichamodel for the , ,
and loss of ozone and thein 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 = http://www.kingceram.com/post/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;realtypedx;//(XMAX-XMIN)/(MX-1)，x方向的变化率，常量realtypedy;// (YMAX-YMIN)/(MY-1)，y方向的变化率，常量realtypehdco;//KH/SQR(data->dx),水平方向上的扩散系数，公式（4）等式右边的第一项realtypehaco;//VEL/(TWO*data->dx),,水平方向上的对流系数，公式（4）等式右边的第二项realtypevdco;//(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= http://www.kingceram.com/post/(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