<div dir="ltr">Hi Matt,<div><br></div><div>Can you please answer my second question about the way PESCs understand the column indies? :</div><div><br class=""><span style="color:rgb(80,0,80)">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.</span><br></div><div><br></div><div>Thank you.</div><div> </div><div>Pham</div><div> </div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Oct 28, 2015 at 8:48 AM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span class="">On Tue, Oct 27, 2015 at 8:18 PM, Sang pham van <span dir="ltr"><<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>></span> wrote:<br></span><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">Hi Barry,<div><br></div><span class=""><div>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). </div><div>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. </div><div>When running the code, I got many of the following error when setting values into the preallocated matrix: </div></span></div></blockquote><div><br></div><div>You might have memory corruption. Run with valgrind.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div><div class="h5"><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div>[0]PETSC ERROR: --------------------- Error Message ------------------------------------</div><div>[0]PETSC ERROR: Argument out of range!</div><div>[0]PETSC ERROR: Local index 1496271904 too large 121500 (max) at 0!</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: Petsc Release Version 3.4.5, Jun, 29, 2014 </div><div>[0]PETSC ERROR: See docs/changes/index.html for recent updates.</div><div>[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.</div><div>[0]PETSC ERROR: See docs/index.html for manual pages.</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: ./sasMAIN on a arch-linux-cxx-opt named sang by sang Wed Oct 28 07:34:37 2015</div><div>[0]PETSC ERROR: Libraries linked from /home/sang/petsc/petsc-3.4.5/arch-linux-cxx-opt/lib</div><div>[0]PETSC ERROR: Configure run at Thu Sep 3 23:04:15 2015</div><div>[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</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: ISLocalToGlobalMappingApply() line 444 in /home/sang/petsc/petsc-3.4.5/src/vec/is/utils/isltog.c</div><div>[0]PETSC ERROR: MatSetValuesLocal() line 1968 in /home/sang/petsc/petsc-3.4.5/src/mat/interface/matrix.c</div><div>[0]PETSC ERROR: MatSetValuesStencil() line 1339 in /home/sang/petsc/petsc-3.4.5/src/mat/interface/matrix.c</div><div><br></div><div>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.</div><div>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.</div><div><br></div><div>Thank you very much.</div><div>Pham </div><div><br></div><div>Below are my code for preallocating the DM matrix; </div><div><br></div><div><div>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)</div><div>{</div><div> PetscErrorCode ierr;</div><div> PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;</div><div> PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;</div><div> PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;</div><div> MPI_Comm comm;</div><div> PetscScalar *values;</div><div> DMDABoundaryType bx,by,bz;</div><div> ISLocalToGlobalMapping ltog,ltogb;</div><div> DMDAStencilType st;</div><div><br></div><div> PetscFunctionBegin;</div><div><br></div><div> ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); // nc = 3</div><div> col = 2*s + 1;</div><div> ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr);<br></div><div> ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr);</div><div> ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr);</div><div> ierr = PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr);<br></div><div> ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr);</div><div> ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr);<br></div><div> nCells = mesh->nCells;</div><div> sasVector<int> cell_type(mesh->nCells,0);</div><div> for(int i=0;i<solidcells.N;i++) cell_type[solidcells[i]] = 1;</div><div> for(int i=0;i<forcing_or_ghost_cells.N;i++) cell_type[forcing_or_ghost_cells[i]] = 2;</div><div> for(i = 0;i<nCells;i++)</div><div> {</div><div> slot = mesh->Cells_ijk[i].i - gxs + gnx*(mesh->Cells_ijk[i].j - gys) + gnx*gny*(mesh->Cells_ijk[i].k - gzs);</div><div> cnt = 0;</div><div> if(cell_type[i]==0) /// fluid cells</div><div> {</div><div> for (l=0; l<nc; l++)</div><div> {</div><div> cols[cnt++] = l + nc*slot;</div><div> rows[l] = l + nc*slot;</div><div> }</div><div> }</div><div><br></div><div> if(cell_type[i]==1) /// solid cells: </div><div> {</div><div> for (l=0; l<nc; l++)</div><div> {</div><div> for (int ll=0; ll<nc; ll++)</div><div> for (ii=-1; ii<2; ii++) {</div><div> for (jj=-1; jj<2; jj++) {</div><div> for (kk=-1; kk<2; kk++) {</div><div> cols[cnt++] = ll + nc*(slot + ii + gnx*jj + gnx*gny*kk);</div><div> }</div><div> }</div><div> }</div><div> rows[l] = l + nc*(slot);</div><div> }</div><div> }</div><div> if(cell_type[i]!=2)</div><div> ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr);</div><div> }<br></div><div> MatCreate(comm,&pMat);<br></div><div> MatSetSizes(pMat,nc*gnx*gny*gnz,nc*gnx*gny*gnz,nc*gnx*gny*gnz,nc*gnx*gny*gnz);</div><div> MatSetType(pMat,MATAIJ);</div><div> ierr = MatSetBlockSize(pMat,nc);CHKERRQ(ierr);</div><div> ierr = MatSeqAIJSetPreallocation(pMat,0,dnz);CHKERRQ(ierr);</div><div> ierr = MatMPIAIJSetPreallocation(pMat,0,dnz,0,onz);CHKERRQ(ierr);</div><div> ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);</div><div> ierr = MatSetLocalToGlobalMapping(pMat,ltog,ltog);CHKERRQ(ierr);</div><div><br></div><div> ierr = PetscFree2(rows,cols);CHKERRQ(ierr);</div><div><br></div><div> PetscFunctionReturn(0);</div><div>}</div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Oct 25, 2015 at 12:46 AM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><span><br>
> On Oct 24, 2015, at 1:43 AM, Sang pham van <<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>> wrote:<br>
><br>
> Thank you, Dave,<br>
><br>
> 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()?<br>
<br>
</span> 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.<br>
<br>
So, in summary you can start by just calling your own routine and convert it to use DMSetGetMatrix() later if you need to.<br>
<span><font color="#888888"><br>
Barry<br>
</font></span><div><div><br>
><br>
> Thank you.<br>
><br>
> Pham<br>
><br>
> On Sat, Oct 24, 2015 at 1:33 PM, Dave May <<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>> wrote:<br>
> If you need to access a user defined context from within your CreateMatrix function, you can attach it to the DM via PetscObjectCompose<br>
><br>
> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscObjectCompose.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscObjectCompose.html</a><br>
><br>
> If your context is not a petsc object, you can use PetscContainerCreate(),<br>
><br>
> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscContainerCreate.html#PetscContainerCreate" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscContainerCreate.html#PetscContainerCreate</a><br>
><br>
> You would then set your user context pointer inside the container and then use PetscObjectCompose() to attach the container to the DM<br>
><br>
> Thanks,<br>
> Dave<br>
><br>
><br>
> On 24 October 2015 at 06:04, Sang pham van <<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>> wrote:<br>
> Hi Barry,<br>
><br>
> 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.<br>
> 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.<br>
><br>
> Thank you.<br>
><br>
> Pham<br>
><br>
><br>
> On Fri, Sep 4, 2015 at 1:24 AM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> 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.<br>
><br>
> Barry<br>
><br>
><br>
><br>
><br>
> > On Sep 3, 2015, at 11:28 AM, Sang pham van <<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>> wrote:<br>
> ><br>
> > Hi,<br>
> ><br>
> > I am using DMCreateMatrix to create matrix from a existed DM object and defined stencil.<br>
> > 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.<br>
> ><br>
> > Many thanks.<br>
> ><br>
> > Sam.<br>
> ><br>
> ><br>
><br>
><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>
</blockquote></div></div></div><span class=""><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</font></span></div></div>
</blockquote></div><br></div></div>