Petsc Добавление значений матрицы
У меня проблема с PETSC. Я написал код в Matlab, и я пытаюсь перевести этот код на C ++ с помощью библиотеки PETSC. Я пишу симуляцию гидродинамики для многофазных потоков, и я пытаюсь сделать эквивалент этой операции Matlab простым способом:
ut(i,j)=u(i,j)+(u(i+1,j)+u(i,j))^2-(u(i,j)+u(i-1,j))*(v(i, j) + v(i-1, j))
Есть ли способ сделать это, не вызывая функцию MATGETVALUES 7 раз?
Решение
Как насчет хранения u v как векторы? Вас может заинтересовать распределенный массив Petsc, см. DMDACreate2d() используя 2 степеней свободы u а также v, Тогда функция DMDAVecGetArrayDOF() превращает вектор из DMDA в массив. В конце концов, ваша строка может написать ut[i][j][U]=ut[i][j][U]+.....+u[i-1][j][V]),
Следующий кусок кода на основе этот пример показывает, как это можно сделать:
static char help[] = "Just a piece of code. See ts/examples/tutorials/ex12.c\n";
#include <petscdm.h>
#include <petscdmda.h>
#undef __FUNCT__
#define __FUNCT__ "main"int main(int argc,char **argv)
{
Vec x,r,xloc,rloc;
PetscErrorCode ierr;
DM da;
PetscScalar ***xarray,***rarray;
PetscInt xm,ym,xs,ys,i,j,Mx,My;
PetscInitialize(&argc,&argv,(char*)0,help);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create distributed array (DMDA) to manage parallel grid and vectors
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
2,1,NULL,NULL,&da);CHKERRQ(ierr);
ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&r);CHKERRQ(ierr);ierr = DMGetLocalVector(da,&xloc);CHKERRQ(ierr);
ierr = DMGetLocalVector(da,&rloc);CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(da,x,INSERT_VALUES,xloc);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(da,x,INSERT_VALUES,xloc);CHKERRQ(ierr);
ierr = DMDAVecGetArrayDOF(da,xloc,&xarray);CHKERRQ(ierr);
ierr = DMDAVecGetArrayDOF(da,rloc,&rarray);CHKERRQ(ierr);
// Get global size
ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
//Get local grid boundaries
ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
//boundary conditions... Your job !
continue;
}
// ut(i,j)=u(i,j)+(u(i+1,j)+u(i,j))^2-(u(i,j)+u(i-1,j))*(v(i, j) + v(i-1, j))
rarray[j][i][0]=xarray[j][i][0]+pow(xarray[j][i+1][0]+xarray[j][i][0],2)-(xarray[j][i][0]+xarray[j][i-1][0])*(xarray[j][i][1]+xarray[j][i-1][1]);
// v component
rarray[j][i][1]=xarray[j][i][1]+pow(xarray[j+1][i][1]+xarray[j][i][1],2)-(xarray[j][i][1]+xarray[j-1][i][1])*(xarray[j][i][0]+xarray[j-1][i][0]);
}
}
// the arrays are placed back in the vector.
ierr = DMDAVecRestoreArrayDOF(da,xloc,&xarray);CHKERRQ(ierr);
ierr = DMDAVecRestoreArrayDOF(da,rloc,&rarray);CHKERRQ(ierr);
// local vector rloc is inserted in global vector.
ierr = DMLocalToGlobalBegin(da,rloc,INSERT_VALUES,r);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(da,rloc,INSERT_VALUES,r);CHKERRQ(ierr);
// local vector are destroyed.
ierr = DMRestoreLocalVector(da,&xloc);CHKERRQ(ierr);
ierr = DMRestoreLocalVector(da,&rloc);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = PetscFinalize();
PetscFunctionReturn(0);
}
Может быть скомпилировано:
mpicc main3.c -o main3 -O3 -I/.../petsc-3.6.4/include -I/.../petsc-3.6.4/linux-gnu-c-debug/include -L/.../petsc-3.6.4/linux-gnu-c-debug/lib -lpetsc
И побежал мимо:
mpirun -np 4 main3 -da_grid_x 100 -da_grid_y 80
Другие решения
Других решений пока нет …