petsc 结构化数据使用

来源:互联网 发布:js时间滑动选择插件 编辑:程序博客网 时间:2024/05/19 18:37

       前文讲到petsc两类基本数据对象(mat, vec)。而petsc的这两类数据对象默认都是按一维存储和创建的。实际pde问题,往往是2维或者3维的。故petsc内建一些映射函数,将默认的一维数据映射成逻辑上的多维数组。本文主要介绍向规则的2D数组映射的函数。


        petsc中逻辑上的多维数组称为 "DA"  ( distributed array), 创建一个2D DA:

 DACreate2d(MPI_Comm comm, DAPeriodicType wrap, DAStencilType stencil_type, PetscInt M, PetscInt N, PetscINt m, PetscInt n, PetscInt dof, PetscInt s, PetscInt *lx,  PetscInt *ly, DA *da) /*      comm MPI通信容器         wrap  周期边界         stencil_type  差分模式(box, star)         M , N  x, y方向整体空间节点数         m, n , x, y方向cpu节点数         dof      空间节点自由度数         s       差分模式的厚度(即,how many 层数 空间节点)         *lx, *ly   cpu节点所管辖的部分节点数组         *da,  所创建DA对象  *

       获取该DA对象的全局及局部信息(即分配到cpu节点上的单元信息)采用一 DAGetLocalInfo (DA da, DALocalInfo *info) ,其中

typedef struct {  PetscInt dim, dof, sw ;   PetscInt mx, my, mz; /* global number of grid points in x, y, z */  PetscInt xs, ys, zs; /* indices to the start of grid points in current cpu, excluiding ghost points(gp) */  PetscInt xm, ym, zm;/* local number of grid points in current cpu, no GP */  PetscInt gxs, gys, gzs; /* similar to xs, but including GP */  PetscInt gxm, gym, gzm; /* similar to xm, but including GP */  DAPeriodicType  pt;   DAStencilType  st;  DA      da;} DALocalInfo ;

       DA对象相当于一个逻辑层,用于管理其所有空间节点以及该通信器里的所有cpu节点。与之对应的物理层,就是DAVec 或者 DAMat。 将数据通过逻辑层管理的好处是直观。 比如,将一个计算域离散成n * n 单元 ,采用DA表述就是我们头脑中的 n * n 的样子,而在物理层仍然是一维线,不直观。

       但是,数据实际上仍然是存在物理层的。因此对数据操作还是要再回到VEC, MAT。 (抽象出来的DA一方面是更直观了,但是又不能直接使用

       正确地讲,DA就是规则网格信息。DAVectGetArray函数用于建立列向量,且该列向量元素个数等于DA网格节点数×自由度数。DAGetMatrix用于建立矩阵A,且矩阵宽度等于DA两个维度之积。

 DACreateGlobalVector(DA da, Vec *global_vec) DACreateLocalVector(DA da, Vec *local_vec) /* 分别声明了一个全局向量和一个本地向量(包含GP),均基于DA对象。*/  DAVecGetArray (DA da, Vec vec, void* array) ;  /* 赋值需要对物理层的VEC进行 */  DAVecRestoreArray(DA da, Vec vec, void* array) ;  /* 赋值完成,再回到逻辑层DA ×/  

   举例:全局(整体)解向量 u(x,y ) =0 , v(x, y) = 1。首先创建DA,(造了空房子);然后通过DA创建GlobalVector,(将空房子与一个GLOBALVECT关联起来 ); 然后获取一个指向该DA(已包含将要赋值的GLOBALVECT)的指针,再通过该指针完成对元素的访问/修改。

DA da;DALocalInfo info;PetscInt i,j;Vec solv;typedef struct  {  PEtscScalar u, v;} Field ;Field **g;DACreate2d(...);DACreateGlobalVector(da, &solv);DAGetLocalInfo(da, &info);DAVecGetArray(da, solv, &g);for(i=info.xs; i < info.xs + info.xm; i++ ){  for(j=info.ys; j<info.ys + info.ym; j++){    g[j][i].u = 0;    g[j][i].v = 1;  } }DAVecRestoreArray(da, solv, &g);

 另外,更新GP采用:

DAGlobalToLocalBegin(DA da, Vec global, InsertMode mode, Vec local);DAGlobalToLocalend(DA da, Vec global, InsertMode mode, Vec local); 
更新GP是不同cpu节点之间的通信,耗时较长。故调用这两个函数之间允许做其他本地操作。

 类似的,DA 到MAT 的映射:

DAGetMatrix(DA da, MatType mtype, Mat *mat) /* create Matrix based on DA array */MatSetValuesStencil(Mat mat, PetscInt m, const MatStencil row[], PetscINt n, const MatStencil col[], const PetscScalar val[], InsertMode mode) /*    m, n :: number of rows and colums;   row[], col[] :: grid indices(associated with stencil structure(pattern));  and components indices;   val ::  values to insert   petsc默认按行划分,col[]实际上就存当前行中不为0的元素。  *//* 对MAT赋值完之后,类似VEC,需要再返映射到“逻辑层”*/  MatAssemblyBegin(Mat mat, MatAssemblyType type);  MatAssemblyEnd(Mat mat, MatAssemblyType type);

举例:一维laplace差分离散

DA da;DALocalInfo info;Mat M;PetscInt i,N=25, num_cols;MatStructure row, col[3];PetscScalar val[3], h=1.00/(N-1);DACreate1d(...);DAGetLocalInfo(da, &info);DAGetMatrix(da, MATMPIAIJ, &M);for(i=info.xs; i < info.xs + info.xm; i++ ){  row.i = i ; /* grid space indices */  row.c = 0 ; /* first component */  if(i ==0 || i = N -1 ) /* B.C. */ {    col[0].i = i ;    col[0].c = 0 ; /* only one component  */    val[0] = 1.00 ;     num_cols = 1;  } else {    col[0].i = i -1;/* left point */    col[0].c = 0;    val[0] = 1/h/h;    col[1].i = i;    col[1].c = 0;    val[1] = -2/h/h;    col[2].i = i+1;    col[2].c = 0;    val[2] = 1/h/h;        num_cols = 3;  }  MatSetValuesStencil(M,1,&row, num_cols, col, val, INSERT_VALUES); }
 MatAssemblyBegin(Mat mat, MatAssemblyType type);  MatAssemblyEnd(Mat mat, MatAssemblyType type);


总结:  规则计算域的数据采用DA管理较直观 (映射成DA);但DA对象不能直接访问/赋值VEC/MAT数据。故需将DA对象返映射回VECT/MAT。然后对于VECT,使用函数DAVecGetArray返回指向该VECT的指针,并由该指针访问/赋值VECT的元素;对于MAT,使用函数MatSetValuesStencil实现访问。

DA管理规则计算域,关联此计算域的方程组A * x = b 的矩阵A,列向量b 都可由DA获取。


际上,DA对象很绕圈子很麻烦。但这正是petsc的设计原则之一:petsc函数对用户数据不可见,所以只好通过指针访问。


0 0