[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Matthew Knepley
knepley at gmail.com
Mon Nov 11 16:39:20 CST 2013
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
>
>
>>
>>> 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/20131111/3ef5c410/attachment-0001.html>
More information about the petsc-users
mailing list