[petsc-users] DMCreateMatrix with some dense row
Matthew Knepley
knepley at gmail.com
Tue Oct 27 20:48:14 CDT 2015
On Tue, Oct 27, 2015 at 8:18 PM, Sang pham van <pvsang002 at gmail.com> wrote:
> Hi Barry,
>
> I made my function for preallocating the DM matrix. I am using immersed
> boundary method to solve a problem of solid mechanics (3D, linear
> elasticity).
> At every solid nodes, I have 3 unknowns for displacement in 3 directions
> (nc=3), my DM has box stencil. Thus, at a node I have 3 equations, for each
> equation I preallocate 27*3 = 81 entries.
> When running the code, I got many of the following error when setting
> values into the preallocated matrix:
>
You might have memory corruption. Run with valgrind.
Thanks,
Matt
> [0]PETSC ERROR: --------------------- Error Message
> ------------------------------------
> [0]PETSC ERROR: Argument out of range!
> [0]PETSC ERROR: Local index 1496271904 too large 121500 (max) at 0!
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: Petsc Release Version 3.4.5, Jun, 29, 2014
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: ./sasMAIN on a arch-linux-cxx-opt named sang by sang Wed
> Oct 28 07:34:37 2015
> [0]PETSC ERROR: Libraries linked from
> /home/sang/petsc/petsc-3.4.5/arch-linux-cxx-opt/lib
> [0]PETSC ERROR: Configure run at Thu Sep 3 23:04:15 2015
> [0]PETSC ERROR: Configure options --download-mpich=1
> --with-scalar-type=real --with-clanguage=cxx --download-mumps=1
> --download-blacs=1 --download-parmetis=1 --download-scalapack=1
> --with-debugging=0 --download-hypre=1 --with-fc=gfortran --download-metis=1
> -download-cmake=1 --download-f-blas-lapack=1
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: ISLocalToGlobalMappingApply() line 444 in
> /home/sang/petsc/petsc-3.4.5/src/vec/is/utils/isltog.c
> [0]PETSC ERROR: MatSetValuesLocal() line 1968 in
> /home/sang/petsc/petsc-3.4.5/src/mat/interface/matrix.c
> [0]PETSC ERROR: MatSetValuesStencil() line 1339 in
> /home/sang/petsc/petsc-3.4.5/src/mat/interface/matrix.c
>
> Can you please explain to me what is possible reason for this error? it
> looks like there is problem with the mapping from LocaltoGobal index, but
> that part I just copy from the original code.
> One more question I have is that, when preallocating for a node, as you
> can see in my below code (based on original PETSc code), since I have 3
> equations, I preallocated 3 rows (in rows[l] array), each row has a number
> of columns index, all column indices of the three rows are stored in the
> cols[] array. At this point I don't understand how PETSc know in that
> cols[] arrays which column indices belonging to the first(second/third)
> row? Please help me to understand this.
>
> Thank you very much.
> Pham
>
> Below are my code for preallocating the DM matrix;
>
> PetscErrorCode sasMatVecPetsc::DMCreateMatrix_DA_3d_MPIAIJ_pvs(DM da,Mat
> pMat, sasSmesh* mesh,sasVector<int>& solidcells,sasVector<int>&
> isolidcellnearBI,sasVector<int>& solidcellnearBI,sasVector<int>&
> forcing_or_ghost_cells)
> {
> PetscErrorCode ierr;
> PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
> PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows =
> NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
> PetscInt
> istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
> MPI_Comm comm;
> PetscScalar *values;
> DMDABoundaryType bx,by,bz;
> ISLocalToGlobalMapping ltog,ltogb;
> DMDAStencilType st;
>
> PetscFunctionBegin;
>
> ierr =
> DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr);
> // nc = 3
> col = 2*s + 1;
> ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr);
> ierr =
> DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr);
> ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr);
> ierr =
> PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr);
> ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr);
> ierr =
> MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr);
> nCells = mesh->nCells;
> sasVector<int> cell_type(mesh->nCells,0);
> for(int i=0;i<solidcells.N;i++) cell_type[solidcells[i]] = 1;
> for(int i=0;i<forcing_or_ghost_cells.N;i++)
> cell_type[forcing_or_ghost_cells[i]] = 2;
> for(i = 0;i<nCells;i++)
> {
> slot = mesh->Cells_ijk[i].i - gxs + gnx*(mesh->Cells_ijk[i].j - gys) +
> gnx*gny*(mesh->Cells_ijk[i].k - gzs);
> cnt = 0;
> if(cell_type[i]==0) /// fluid cells
> {
> for (l=0; l<nc; l++)
> {
> cols[cnt++] = l + nc*slot;
> rows[l] = l + nc*slot;
> }
> }
>
> if(cell_type[i]==1) /// solid cells:
> {
> for (l=0; l<nc; l++)
> {
> for (int ll=0; ll<nc; ll++)
> for (ii=-1; ii<2; ii++) {
> for (jj=-1; jj<2; jj++) {
> for (kk=-1; kk<2; kk++) {
> cols[cnt++] = ll + nc*(slot + ii + gnx*jj + gnx*gny*kk);
> }
> }
> }
> rows[l] = l + nc*(slot);
> }
> }
> if(cell_type[i]!=2)
> ierr =
> MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr);
> }
> MatCreate(comm,&pMat);
>
> MatSetSizes(pMat,nc*gnx*gny*gnz,nc*gnx*gny*gnz,nc*gnx*gny*gnz,nc*gnx*gny*gnz);
> MatSetType(pMat,MATAIJ);
> ierr = MatSetBlockSize(pMat,nc);CHKERRQ(ierr);
> ierr = MatSeqAIJSetPreallocation(pMat,0,dnz);CHKERRQ(ierr);
> ierr = MatMPIAIJSetPreallocation(pMat,0,dnz,0,onz);CHKERRQ(ierr);
> ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
> ierr = MatSetLocalToGlobalMapping(pMat,ltog,ltog);CHKERRQ(ierr);
>
> ierr = PetscFree2(rows,cols);CHKERRQ(ierr);
>
> PetscFunctionReturn(0);
> }
>
> On Sun, Oct 25, 2015 at 12:46 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>> > On Oct 24, 2015, at 1:43 AM, Sang pham van <pvsang002 at gmail.com> wrote:
>> >
>> > Thank you, Dave,
>> >
>> > Can I just call My_DMCreateMatrix_DA_3d_MPIAIJ(DM da, Mat J, void*
>> my_context) when I need to create my matrix, and forget about
>> DMDASetGetMatrix()?
>>
>> You can do that if that is all you need. In some contexts when
>> something else in PETSc needs to create a specific matrix (like with
>> geometric multigrid) then that code calls DMCreateMatrix() which would then
>> use your creation routine.
>>
>> So, in summary you can start by just calling your own routine and
>> convert it to use DMSetGetMatrix() later if you need to.
>>
>> Barry
>>
>> >
>> > Thank you.
>> >
>> > Pham
>> >
>> > On Sat, Oct 24, 2015 at 1:33 PM, Dave May <dave.mayhem23 at gmail.com>
>> wrote:
>> > If you need to access a user defined context from within your
>> CreateMatrix function, you can attach it to the DM via PetscObjectCompose
>> >
>> >
>> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscObjectCompose.html
>> >
>> > If your context is not a petsc object, you can use
>> PetscContainerCreate(),
>> >
>> >
>> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscContainerCreate.html#PetscContainerCreate
>> >
>> > You would then set your user context pointer inside the container and
>> then use PetscObjectCompose() to attach the container to the DM
>> >
>> > Thanks,
>> > Dave
>> >
>> >
>> > On 24 October 2015 at 06:04, Sang pham van <pvsang002 at gmail.com> wrote:
>> > Hi Barry,
>> >
>> > The routine DMCreateMatrix_DA_3d_MPIAIJ has 2 input arguments (DM da,
>> Mat J), the function pointer in DMDASetGetMatrix() only accept function
>> with that two arguments.
>> > As you suggested, I am writing a routine (based on
>> DMCreateMatrix_DA_3d_MPIAIJ()) to preallocate the matrix in the way I wish,
>> to do that I think It needs to have one more input argument:
>> My_DMCreateMatrix_DA_3d_MPIAIJ(DM da, Mat J, void* my_context). But
>> DMDASetGetMatrix() will not accept pointer of this function. Please give a
>> suggestion to overcome this.
>> >
>> > Thank you.
>> >
>> > Pham
>> >
>> >
>> > On Fri, Sep 4, 2015 at 1:24 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>> >
>> > Look at DMCreateMatrix_DA_2d_MPIAIJ (or the 3d version if working
>> in 3d) You need to copy this routine and add whatever additional
>> preallocation information you need. Then call DMDASetGetMatrix() so that
>> the DM will use your routine to create the matrix for you.
>> >
>> > Barry
>> >
>> >
>> >
>> >
>> > > On Sep 3, 2015, at 11:28 AM, Sang pham van <pvsang002 at gmail.com>
>> wrote:
>> > >
>> > > Hi,
>> > >
>> > > I am using DMCreateMatrix to create matrix from a existed DM object
>> and defined stencil.
>> > > In my code, boundary nodes need to involve many inner nodes, thus
>> matrix rows corresponding to boundary nodes are almost dense. How can I
>> tell petsc that those rows need to be preallocated with more entries? I
>> don't want to use MatMPIAIJSetPreallocation() since advantages of DM might
>> be lost.
>> > >
>> > > Many thanks.
>> > >
>> > > Sam.
>> > >
>> > >
>> >
>> >
>> >
>> >
>>
>>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151027/dc639bbc/attachment.html>
More information about the petsc-users
mailing list