[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Bishesh Khanal
bisheshkh at gmail.com
Tue Nov 12 08:27:03 CST 2013
On Mon, Nov 11, 2013 at 11:39 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Mon, Nov 11, 2013 at 1:08 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>
>> On Mon, Nov 4, 2013 at 5:25 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>
>>> On Mon, Nov 4, 2013 at 4:54 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>
>>>> On Mon, Nov 4, 2013 at 4:49 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>
>>>>> On Fri, Nov 1, 2013 at 1:30 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>
>>>>>> On Fri, Nov 1, 2013 at 12:08 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>>>
>>>>>>> On Thu, Oct 31, 2013 at 11:34 PM, Matthew Knepley <knepley at gmail.com
>>>>>>> > wrote:
>>>>>>>
>>>>>>>> On Thu, Oct 31, 2013 at 5:02 PM, Bishesh Khanal <
>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> On Thu, Oct 31, 2013 at 7:43 PM, Matthew Knepley <
>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> 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
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Up to 4 (a), it is correct. There is a problem in the way
>>>>>>>>> DMDAGetCellPoint and PetscSectionGetOffset works in this case scenario. I
>>>>>>>>> will try to explain it below for the case of 4X4 grid.
>>>>>>>>> First case:
>>>>>>>>> If I set the computational domain to be all the points of the dmda
>>>>>>>>> grid, (i.e. isPosInDomain(..,i,j,..) returns true for all i and j in the
>>>>>>>>> dmda grid), the program runs fine and does not give any error.
>>>>>>>>>
>>>>>>>>> Second case:
>>>>>>>>> I want the computational domain to be some part of the whole grid.
>>>>>>>>> There is a problem in this case.
>>>>>>>>> The following test is in a case where,
>>>>>>>>> isPosInDomain(..,i,j,..) returns true ONLY for (i,j) pairs of
>>>>>>>>> (2,1) (1,2) (2,2) (3,2) and (3,2). The grid with its corresponding point
>>>>>>>>> number in the petscsection is shown below:
>>>>>>>>>
>>>>>>>>> 12 13 14 15
>>>>>>>>> 8 9 10 11
>>>>>>>>> 4 5 6 7
>>>>>>>>> 0 1 2 3
>>>>>>>>>
>>>>>>>>> of which 6, 9, 10, 11 and 14 correspond to the (i,j) points that
>>>>>>>>> returns true for isPosInDomain(..,i,j,..)
>>>>>>>>> MatView gives me the 16-row matrix with the star stencil non-zero
>>>>>>>>> structure as expected.
>>>>>>>>> The error I get is: new non-zero at (0,2) caused a malloc.
>>>>>>>>>
>>>>>>>>> This error is for the value of (i,j) = (2,1), i.e. point 6.
>>>>>>>>> The loop to set values in the matrix starts as:
>>>>>>>>> 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);
>>>>>>>>>
>>>>>>>>> Now, what I wanted from the beginning is to create the matrix
>>>>>>>>> containing the rows corresponding to only those points (i,j) which has true
>>>>>>>>> values for isPosInDomain(ctx,i,j), that means in this case 5 row matrix. In
>>>>>>>>> the current test, the first (i,j) that reaches DMDAGetCellPoint is (2,1),
>>>>>>>>> which gives sets point variable to 6.
>>>>>>>>> The line PetscSectionGetOffset(gs,point,&row) sets the value of
>>>>>>>>> row to 0.
>>>>>>>>> So I think this where the inconsistency lies.
>>>>>>>>> MatSetValues(A,1,&row,num,col,v,INSERT_VALUES) expects row
>>>>>>>>> variable and col[] array to correspond to (i,j) based on DMDA, but
>>>>>>>>> PetsSectionGetOffset gives result based on how we masked away some of the
>>>>>>>>> points from dmda grid. Or am I missing something very basic here ?
>>>>>>>>>
>>>>>>>>
>>>>>>>> You are missing something basic. The row is correctly identified as
>>>>>>>> 0, and that means the 2 is
>>>>>>>> point 10.
>>>>>>>>
>>>>>>> The problem must be preallocation. First, you should send the result
>>>>>>>> of MatView().
>>>>>>>>
>>>>>>>
>>>>>>> MatView shows that the preallocation and the matrix created is for
>>>>>>> the full dmda grid which is probably not what we want here, right ? Here is
>>>>>>> the matview result:
>>>>>>>
>>>>>>
>>>>>> Thats right. Did you call
>>>>>>
>>>>>> DMSetDefaultSection(da, s)
>>>>>>
>>>>>> after you made the new section and before you called DMGetMatrix()?
>>>>>>
>>>>> Do you mean before I called DMCreateMatrix() ? I called
>>>>> DMSetDefaultSection(da,s) as below:
>>>>>
>>>>> ...
>>>>> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
>>>>> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
>>>>> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
>>>>> // ierr = DMGetMatrix(dm,&A);CHKERRQ(ierr); //It seems DMGetMatrix
>>>>> is changed to DMCreateMatrix().
>>>>> ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);
>>>>> ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
>>>>>
>>>>
>>>> Yes, this is correct. I do exactly this in the example code and it
>>>> works fine. I suspect
>>>> you are not creating the section you think you are. You can send the
>>>> code if you want
>>>> and I will look at it.
>>>>
>>>
>>> Thanks, here is the code:
>>> https://github.com/bishesh/petscPoissonIrregular/blob/master/poissonIrregular.cxx
>>> It works for 2D and 3D without using PetscSection. (The choice of using
>>> Petsc_section is on line 68). Using PetsSection currently it implements
>>> only the 2D case with dof = 1. It works when the domain consists of all the
>>> cells of the dmda, it works (to check you can uncomment line 129 in
>>> isPosInDomain(). But it does not when I mask out some of the cells while
>>> creating the petscsection.
>>>
>>
>> Hi Matt, I'm afraid if my last email somehow escaped you. Could you
>> please let me know what could be possibly going wrong with my code linked
>> above ?
>>
>
> There is a bug in DMCreateMatrix() for DMDA when using a PetscSection. I
> am fixing it. I will mail when its done.
>
Thanks Matt, I've one more application where I think PetscSection might be
useful in stokes like equation with staggered grid; I will start another
thread for that while waiting for the bug to get fixed in this thread.
>
> Thanks,
>
> Matt
>
>
>>
>>
>>>
>>>> Matt
>>>>
>>>>
>>>>> This MatView on A shows (result below) that it still has 16 rows for
>>>>> all DMDA points instead of only those petscsection points for which I set a
>>>>> non-zero dof. Is it what it should not be doing ?
>>>>>
>>>>>> Mat Object: 1 MPI processes
>>>>>>>
>>>>>>> type: seqaij
>>>>>>>
>>>>>>> row 0: (0, 0) (1, 0) (4, 0)
>>>>>>>
>>>>>>> row 1: (0, 0) (1, 0) (2, 0) (5, 0)
>>>>>>>
>>>>>>> row 2: (1, 0) (2, 0) (3, 0) (6, 0)
>>>>>>>
>>>>>>> row 3: (2, 0) (3, 0) (7, 0)
>>>>>>>
>>>>>>> row 4: (0, 0) (4, 0) (5, 0) (8, 0)
>>>>>>>
>>>>>>> row 5: (1, 0) (4, 0) (5, 0) (6, 0) (9, 0)
>>>>>>>
>>>>>>> row 6: (2, 0) (5, 0) (6, 0) (7, 0) (10, 0)
>>>>>>>
>>>>>>> row 7: (3, 0) (6, 0) (7, 0) (11, 0)
>>>>>>>
>>>>>>> row 8: (4, 0) (8, 0) (9, 0) (12, 0)
>>>>>>>
>>>>>>> row 9: (5, 0) (8, 0) (9, 0) (10, 0) (13, 0)
>>>>>>>
>>>>>>> row 10: (6, 0) (9, 0) (10, 0) (11, 0) (14, 0)
>>>>>>>
>>>>>>> row 11: (7, 0) (10, 0) (11, 0) (15, 0)
>>>>>>>
>>>>>>> row 12: (8, 0) (12, 0) (13, 0)
>>>>>>>
>>>>>>> row 13: (9, 0) (12, 0) (13, 0) (14, 0)
>>>>>>>
>>>>>>> row 14: (10, 0) (13, 0) (14, 0) (15, 0)
>>>>>>>
>>>>>>> row 15: (11, 0) (14, 0) (15, 0)
>>>>>>>
>>>>>>> [0]PETSC ERROR: --------------------- Error Message
>>>>>>> ------------------------------------
>>>>>>>
>>>>>>> [0]PETSC ERROR: Argument out of range!
>>>>>>>
>>>>>>> [0]PETSC ERROR: New nonzero at (0,2) caused a malloc!
>>>>>>>
>>>>>>>
>>>>>>>
>>>>> For reference here is the code snippet again:
>>>>> ...
>>>>> 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);
>>>>> ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
>>>>>
>>>>> 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);
>>>>>
>>>>> ....
>>>>> The function computeMatrix2dSection:
>>>>>
>>>>> #undef __FUNCT__
>>>>> #define __FUNCT__ "computeMatrix2dSection"
>>>>> PetscErrorCode computeMatrix2dSection(KSP ksp, Mat A, Mat B,
>>>>> MatStructure *str, void *context) {
>>>>> PoissonModel *ctx = (PoissonModel*)context;
>>>>> PetscErrorCode ierr;
>>>>> DM da;
>>>>> DMDALocalInfo info;
>>>>> PetscInt row, col[5];
>>>>> PetscInt dof;
>>>>> PetscScalar v[5]; //array to store 5 point stencil for
>>>>> each row
>>>>> PetscInt num; //non-zero position in the current row
>>>>> // PetscScalar dScale = 1; //to scale Dirichlet
>>>>> identity rows if needed
>>>>> PetscSection gs; //Global Section that
>>>>> keeps the grid info and indices
>>>>> PetscInt point; //Current point of the
>>>>> petscSection
>>>>>
>>>>> PetscFunctionBeginUser;
>>>>> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
>>>>> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
>>>>> ierr = DMGetDefaultGlobalSection(da,&gs);CHKERRQ(ierr);
>>>>> ierr = DMCreateMatrix(da,&A);CHKERRQ(ierr);
>>>>> // ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);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 < ctx->mDof; ++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);
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>>
>>>>>>>> Matt
>>>>>>>>
>>>>>>>>
>>>>>>>>>
>>>>>>>>>> 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
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> 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/20131112/9ef74fa0/attachment-0001.html>
More information about the petsc-users
mailing list