Haustiere addieren Matrixwerte

Question

Ich habe ein Problem mit PETSC. Ich habe einen Code in Matlab geschrieben und ich versuche, diesen Code mit der PETSC-Bibliothek nach C++ zu übersetzen. Ich bin eine Strömungssimulation für Mehrphasenströmungen zu schreiben und ich versuche, das äquivalent diesen Matlab Betrieb auf einfache Art und Weise zu tun:

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)) 

Gibt es eine Möglichkeit, dies zu tun, ohne die MATGETVALUES Funktion 7 mal anrufen zu müssen ?

Answer 1

Wie wäre es mit der Speicherung von u und v als Vektoren? Sie können von der verteilten Anordnung von Petsc interessiert sein, siehe DMDACreate2d() mit 2 Dofs u und v. Dann verwandelt die Funktion DMDAVecGetArrayDOF() einen Vektor von einer DMDA in ein Array. Am Ende könnte Ihre Zeile ut[i][j][U]=ut[i][j][U]+.....+u[i-1][j][V]) schreiben.

Das folgende Stück Code basiert auf this example zeigt, wie es gemacht werden kann:

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 

und lief durch:

mpirun -np 4 main3 -da_grid_x 100 -da_grid_y 80 

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); 
} 

Es wird von kompiliert werden kann