[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Matthew Knepley
knepley at gmail.com
Thu Oct 31 13:43:05 CDT 2013
On Thu, Oct 31, 2013 at 1:25 PM, Bishesh Khanal <bisheshkh at gmail.com> wrote:
> I did not call DMCreateMatrix() before when using just dmda (and it is
> working). In any case, I added a call to DMCreateMatrix() now but still
> there are problems. Here are the code snippets and the error I get:
>
Then you were allowing it to be called automatically by PETSc, and it would
have been this time as well.
> //Setting up the section and linking it to DM:
> ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr);
> ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
> ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr);
> for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {
> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {
> PetscInt point;
> if(isPosInDomain(&testPoisson,i,j,0)) {
> ierr = DMDAGetCellPoint(dm, i, j, 0, &point);CHKERRQ(ierr);
> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); //
> No. of dofs associated with the point.
> }
> }
> }
> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
> ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);
>
> //Set up KSP:
> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
> ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);
> ierr =
> KSPSetComputeOperators(ksp,computeMatrix2dSection,(void*)&testPoisson);CHKERRQ(ierr);
> ierr =
> KSPSetComputeRHS(ksp,computeRhs2dSection,(void*)&testPoisson);CHKERRQ(ierr);
> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
> ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
>
> ------------------------------------------------------
> The computeMatrix2dSection function has:
>
> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
> ierr = DMGetDefaultGlobalSection(da,&gs);CHKERRQ(ierr);
> for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {
> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {
> if (isPosInDomain(ctx,i,j)) {
> ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);
> ierr = PetscSectionGetOffset(gs,point,&row);
> ierr = PetscSectionGetDof(gs,point,&dof);
> for(PetscInt cDof = 0; cDof < dof; ++cDof) {
> num = 0;
> row+=cDof;
> col[num] = row; //(i,j) position
> v[num++] = -4;
> if(isPosInDomain(ctx,i+1,j)) {
> ierr =
> DMDAGetCellPoint(da,i+1,j,0,&point);CHKERRQ(ierr);
> ierr = PetscSectionGetOffset(gs,point,&col[num]);
> col[num] += cDof;
> v[num++] = 1;
> }
> if(isPosInDomain(ctx,i-1,j)) {
> ierr =
> DMDAGetCellPoint(da,i-1,j,0,&point);CHKERRQ(ierr);
> ierr = PetscSectionGetOffset(gs,point,&col[num]);
> col[num] += cDof;
> v[num++] = 1;
> }
> if(isPosInDomain(ctx,i,j+1)) {
> ierr =
> DMDAGetCellPoint(da,i,j+1,0,&point);CHKERRQ(ierr);
> ierr = PetscSectionGetOffset(gs,point,&col[num]);
> col[num] += cDof;
> v[num++] = 1;
> }
> if(isPosInDomain(ctx,i,j-1)) {
> ierr =
> DMDAGetCellPoint(da,i,j-1,0,&point);CHKERRQ(ierr);
> ierr = PetscSectionGetOffset(gs,point,&col[num]);
> col[num] += cDof;
> v[num++] = 1;
> }
> ierr =
> MatSetValues(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
> }
> }
> }
> }
> ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>
> But this is not working. For e.g. for a 6X6 grid size I get the following
> error:
>
Okay, here is how to debug this:
0) Go to a single scalar equations to make things easier
1) Use MatView(A, PETSC_VIEWER_STDOUT_WORLD), which should output a 36
row matrix with
the preallocated nonzero pattern, filled with zeros
2) Make sure this is the pattern you want
3) Run in the debugger with -start_in_debugger
4) When you get the error, see
a) If the (i, j) is one that should be setting a value
b) Why this (i, j) was not preallocated
Thanks,
Matt
> [0]PETSC ERROR: --------------------- Error Message
> ------------------------------------
> [0]PETSC ERROR: Argument out of range!
> [0]PETSC ERROR: New nonzero at (0,6) caused a malloc!
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: Petsc Development GIT revision:
> 61e5e40bb2c5bf2423e94b71a15fef47e411b0da GIT Date: 2013-10-25 21:47:45
> -0500
> [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: build/poissonIrregular on a arch-linux2-cxx-debug named
> edwards by bkhanal Thu Oct 31 19:23:58 2013
> [0]PETSC ERROR: Libraries linked from
> /home/bkhanal/Documents/softwares/petsc/arch-linux2-cxx-debug/lib
> [0]PETSC ERROR: Configure run at Sat Oct 26 16:35:15 2013
> [0]PETSC ERROR: Configure options --download-mpich
> -download-f-blas-lapack=1 --download-metis --download-parmetis
> --download-superlu_dist --download-scalapack --download-mumps
> --download-hypre --with-clanguage=cxx
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: MatSetValues_SeqAIJ() line 413 in
> src/mat/impls/aij/seq/aij.c
> [0]PETSC ERROR: MatSetValues() line 1130 in src/mat/interface/matrix.c
> [0]PETSC ERROR: computeMatrix2dSection() line 318 in
> /user/bkhanal/home/works/poissonIrregular/poissonIrregular.cxx
> [0]PETSC ERROR: KSPSetUp() line 228 in src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: KSPSolve() line 399 in src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: main() line 598 in
> /user/bkhanal/home/works/poissonIrregular/poissonIrregular.cxx
> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0
> [unset]: aborting job:
> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0
>
>
>> Matt
>>
>>
>>> However, now since I do not want the rows in the matrix corresponding to
>>> the points in the grid which do not lie in the computational domain. This
>>> means the size of the matrix will not necessarily equal the total number of
>>> cells in DMDA. So in the following code:
>>>
>>> ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
>>> ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr);
>>> ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
>>> ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr);
>>> for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {
>>> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {
>>> PetscInt point;
>>> if(isPosInDomain(&testPoisson,i,j,0)) {
>>> ierr = DMDAGetCellPoint(dm, i, j, 0,
>>> &point);CHKERRQ(ierr);
>>> ierr = PetscSectionSetDof(s, point, testPoisson.mDof);
>>> // No. of dofs associated with the point.
>>> }
>>>
>>> }
>>> }
>>> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
>>> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
>>> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
>>>
>>> should I myself compute proper nC (i.e. total no. of points for which I
>>> want the matrix to have a corresponding row) before calling
>>> PetscSectionSetChart() or is the masking out of the points inside the loop
>>> with if(isPosInDomain(&testPoisson,i,j,0)) sufficient ?
>>> And, when you write:
>>> PetscSectionGetOffset(gs, p, &ind);
>>> row = ind < 0 ? -(ind+1) : ind;
>>>
>>> it seems you allow the possibility of getting a -ve ind when using
>>> PetscSectionGetOffset. When would an offset to a particular dof and point?
>>>
>>>
>>>
>>>
>>>>
>>>>> for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {
>>>>> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {
>>>>> num = 0;
>>>>> /*ierr =
>>>>> DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);*/ /*Get the PetscPoint*/
>>>>> /*But I did now understand how I would emulate the
>>>>> row.c and col.c info with PetscSection*/
>>>>> row.i = i; row.j = j;
>>>>>
>>>>
>>>> Here you would call
>>>>
>>>> ierr = DMDAGetCellPoint(da, i, j, 0, &cell);CHKERRQ(ierr);
>>>> ierr = PetscSectionGetOffset(da, cell, &row);CHKERRQ(ierr);
>>>> ierr = PetscSectionGetDof(da, cell, &dof);CHKERRQ(ierr);
>>>>
>>>> This means that this cell has global rows = [row, row+dof).
>>>>
>>>>
>>>>> col[num].i = i; col[num].j = j;
>>>>> if (isPosInDomain(ctx,i,j)) { /*put the operator
>>>>> for only the values for only the points lying in the domain*/
>>>>> v[num++] = -4;
>>>>> if(isPosInDomain(ctx,i+1,j)) {
>>>>> col[num].i = i+1; col[num].j = j;
>>>>> v[num++] = 1;
>>>>> }
>>>>> if(isPosInDomain(ctx,i-1,j)) {
>>>>> col[num].i = i-1; col[num].j = j;
>>>>> v[num++] = 1;
>>>>> }
>>>>> if(isPosInDomain(ctx,i,j+1)) {
>>>>> col[num].i = i; col[num].j = j+1;
>>>>> v[num++] = 1;
>>>>> }
>>>>> if(isPosInDomain(ctx,i,j-1)) {
>>>>> col[num].i = i; col[num].j = j-1;
>>>>> v[num++] = 1;
>>>>> }
>>>>> } else {
>>>>> v[num++] = dScale; /*set Dirichlet identity rows
>>>>> for points in the rectangle but outside the computational domain*/
>>>>> }
>>>>>
>>>>
>>>> You do the same thing for the columns, and then call
>>>>
>>>> ierr = MatSetValues(A, dof, rows, num*dof, cols, v,
>>>> INSERT_VALUES);CHKERRQ(ierr);
>>>>
>>>>
>>>>> ierr =
>>>>> MatSetValuesStencil(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
>>>>> }
>>>>> }
>>>>> }
>>>>>
>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> Matt
>>>>>>
>>>>>>
>>>>>>>
>>>>>>>> Matt
>>>>>>>>
>>>>>>>>
>>>>>>>>>
>>>>>>>>>> Matt
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>> Thanks,
>>>>>>>>>>>>
>>>>>>>>>>>> Matt
>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>> If not, just put the identity into
>>>>>>>>>>>>>>>> the rows you do not use on the full cube. It will not hurt
>>>>>>>>>>>>>>>> scalability or convergence.
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> In the case of Poisson with Dirichlet condition this might
>>>>>>>>>>>>>>> be the case. But is it always true that having identity rows in the system
>>>>>>>>>>>>>>> matrix will not hurt convergence ? I thought otherwise for the following
>>>>>>>>>>>>>>> reasons:
>>>>>>>>>>>>>>> 1) Having read Jed's answer here :
>>>>>>>>>>>>>>> http://scicomp.stackexchange.com/questions/3426/why-is-pinning-a-point-to-remove-a-null-space-bad/3427#3427
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Jed is talking about a constraint on a the pressure at a
>>>>>>>>>>>>>> point. This is just decoupling these unknowns from the rest
>>>>>>>>>>>>>> of the problem.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> 2) Some observation I am getting (but I am still doing more
>>>>>>>>>>>>>>> experiments to confirm) while solving my staggered-grid 3D stokes flow with
>>>>>>>>>>>>>>> schur complement and using -pc_type gamg for A00 matrix. Putting the
>>>>>>>>>>>>>>> identity rows for dirichlet boundaries and for ghost cells seemed to have
>>>>>>>>>>>>>>> effects on its convergence. I'm hoping once I know how to use PetscSection,
>>>>>>>>>>>>>>> I can get rid of using ghost cells method for the staggered grid and get
>>>>>>>>>>>>>>> rid of the identity rows too.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> It can change the exact iteration, but it does not make the
>>>>>>>>>>>>>> matrix conditioning worse.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Anyway please provide me with some pointers so that I can
>>>>>>>>>>>>>>> start trying with petscsection on top of a dmda, in the beginning for
>>>>>>>>>>>>>>> non-staggered case.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Thanks,
>>>>>>>>>>>>>>> Bishesh
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Thanks,
>>>>>>>>>>>>>>>>> Bishesh
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> --
>>>>>>>>>>>>>>>> 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
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> --
>>>>>>>>>>>>>> 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
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> --
>>>>>>>>>>>> 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
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> 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
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> 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
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> 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
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>> --
>>>> 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
>>>>
>>>
>>>
>>
>>
>> --
>> 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
>>
>
>
--
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/20131031/cb57b64d/attachment-0001.html>
More information about the petsc-users
mailing list