CVODE的cvDiurnal_kry例子学习

原创

已于 2022-03-14 10:40:11 修改 · 765 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#SUNDIALS

于 2022-03-14 05:01:50 首次发布

本文详细介绍了一种程序，用于解决大气层中氧气和臭氧的对流-扩散方程的半离散化形式，涉及非线性日间动力学、水平对流和非均匀垂直扩散。通过10x10网格划分，展示了从初始分布计算到完整求解的过程，包括关键代码片段和数据结构。

问题描述

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

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

PDE方程

在这里插入图片描述

c1,c2分别是两种成分的微分方程
变量：水平方向的坐标x,垂直方向的坐标y，以及时间t,
K_h:水平扩散系数，常量=4.0e-6
V：对流系数，常量=0.001
K_v:表达式=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)

最低0.47元/天解锁文章