[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Bishesh Khanal
bisheshkh at gmail.com
Thu Oct 31 13:25:19 CDT 2013
On Thu, Oct 31, 2013 at 6:13 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Thu, Oct 31, 2013 at 12:02 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>
>> On Tue, Oct 29, 2013 at 5:37 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>
>>> On Tue, Oct 29, 2013 at 5:31 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>
>>>> On Mon, Oct 28, 2013 at 9:00 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>
>>>>> On Mon, Oct 28, 2013 at 2:51 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>>
>>>>>> On Mon, Oct 28, 2013 at 6:19 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>>
>>>>>>> On Mon, Oct 28, 2013 at 10:13 AM, Bishesh Khanal <
>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>
>>>>>>>> On Mon, Oct 28, 2013 at 3:49 PM, Matthew Knepley <knepley at gmail.com
>>>>>>>> > wrote:
>>>>>>>>
>>>>>>>>> On Mon, Oct 28, 2013 at 9:06 AM, Bishesh Khanal <
>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> On Mon, Oct 28, 2013 at 1:40 PM, Matthew Knepley <
>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>
>>>>>>>>>>> On Mon, Oct 28, 2013 at 5:30 AM, Bishesh Khanal <
>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> On Fri, Oct 25, 2013 at 10:21 PM, Matthew Knepley <
>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>> On Fri, Oct 25, 2013 at 2:55 PM, Bishesh Khanal <
>>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Fri, Oct 25, 2013 at 8:18 PM, Matthew Knepley <
>>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> On Fri, Oct 25, 2013 at 12:09 PM, Bishesh Khanal <
>>>>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Dear all,
>>>>>>>>>>>>>>>> I would like to know if some of the petsc objects that I
>>>>>>>>>>>>>>>> have not used so far (IS, DMPlex, PetscSection) could be useful in the
>>>>>>>>>>>>>>>> following case (of irregular domains):
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Let's say that I have a 3D binary image (a cube).
>>>>>>>>>>>>>>>> The binary information of the image partitions the cube
>>>>>>>>>>>>>>>> into a computational domain and non-computational domain.
>>>>>>>>>>>>>>>> I must solve a pde (say a Poisson equation) only on the
>>>>>>>>>>>>>>>> computational domains (e.g: two isolated spheres within the cube). I'm
>>>>>>>>>>>>>>>> using finite difference and say a dirichlet boundary condition
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> I know that I can create a dmda that will let me access the
>>>>>>>>>>>>>>>> information from this 3D binary image, get all the coefficients, rhs values
>>>>>>>>>>>>>>>> etc using the natural indexing (i,j,k).
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Now, I would like to create a matrix corresponding to the
>>>>>>>>>>>>>>>> laplace operator (e.g. with standard 7 pt. stencil), and the corresponding
>>>>>>>>>>>>>>>> RHS that takes care of the dirchlet values too.
>>>>>>>>>>>>>>>> But in this matrix it should have the rows corresponding to
>>>>>>>>>>>>>>>> the nodes only on the computational domain. It would be nice if I can
>>>>>>>>>>>>>>>> easily (using (i,j,k) indexing) put on the rhs dirichlet values
>>>>>>>>>>>>>>>> corresponding to the boundary points.
>>>>>>>>>>>>>>>> Then, once the system is solved, put the values of the
>>>>>>>>>>>>>>>> solution back to the corresponding positions in the binary image.
>>>>>>>>>>>>>>>> Later, I might have to extend this for the staggered grid
>>>>>>>>>>>>>>>> case too.
>>>>>>>>>>>>>>>> So is petscsection or dmplex suitable for this so that I
>>>>>>>>>>>>>>>> can set up the matrix with something like DMCreateMatrix ? Or what would
>>>>>>>>>>>>>>>> you suggest as a suitable approach to this problem ?
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> I have looked at the manual and that led me to search for a
>>>>>>>>>>>>>>>> simpler examples in petsc src directories. But most of the ones I
>>>>>>>>>>>>>>>> encountered are with FEM (and I'm not familiar at all with FEM, so these
>>>>>>>>>>>>>>>> examples serve more as a distraction with FEM jargon!)
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> It sounds like the right solution for this is to use
>>>>>>>>>>>>>>> PetscSection on top of DMDA. I am working on this, but it is really
>>>>>>>>>>>>>>> alpha code. If you feel comfortable with that level of
>>>>>>>>>>>>>>> development, we can help you.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Thanks, with the (short) experience of using Petsc so far and
>>>>>>>>>>>>>> being familiar with the awesomeness (quick and helpful replies) of this
>>>>>>>>>>>>>> mailing list, I would like to give it a try. Please give me some pointers
>>>>>>>>>>>>>> to get going for the example case I mentioned above. A simple example of
>>>>>>>>>>>>>> using PetscSection along with DMDA for finite volume (No FEM) would be
>>>>>>>>>>>>>> great I think.
>>>>>>>>>>>>>> Just a note: I'm currently using the petsc3.4.3 and have not
>>>>>>>>>>>>>> used the development version before.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> Okay,
>>>>>>>>>>>>>
>>>>>>>>>>>>> 1) clone the repository using Git and build the 'next' branch.
>>>>>>>>>>>>>
>>>>>>>>>>>>> 2) then we will need to create a PetscSection that puts
>>>>>>>>>>>>> unknowns where you want them
>>>>>>>>>>>>>
>>>>>>>>>>>>> 3) Setup the solver as usual
>>>>>>>>>>>>>
>>>>>>>>>>>>> You can do 1) an 3) before we do 2).
>>>>>>>>>>>>>
>>>>>>>>>>>>> I've done 1) and 3). I have one .cxx file that solves the
>>>>>>>>>>>> system using DMDA (putting identity into the rows corresponding to the
>>>>>>>>>>>> cells that are not used).
>>>>>>>>>>>> Please let me know what I should do now.
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Okay, now write a loop to setup the PetscSection. I have the
>>>>>>>>>>> DMDA partitioning cells, so you would have
>>>>>>>>>>> unknowns in cells. The code should look like this:
>>>>>>>>>>>
>>>>>>>>>>> PetscSectionCreate(comm, &s);
>>>>>>>>>>> DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
>>>>>>>>>>> PetscSectionSetChart(s, 0, nC);
>>>>>>>>>>> for (k = zs; k < zs+zm; ++k) {
>>>>>>>>>>> for (j = ys; j < ys+ym; ++j) {
>>>>>>>>>>> for (i = xs; i < xs+xm; ++i) {
>>>>>>>>>>> PetscInt point;
>>>>>>>>>>>
>>>>>>>>>>> DMDAGetCellPoint(dm, i, j, k, &point);
>>>>>>>>>>> PetscSectionSetDof(s, point, dof); // You know how many
>>>>>>>>>>> dof are on each vertex
>>>>>>>>>>> }
>>>>>>>>>>> }
>>>>>>>>>>> }
>>>>>>>>>>> PetscSectionSetUp(s);
>>>>>>>>>>> DMSetDefaultSection(dm, s);
>>>>>>>>>>> PetscSectionDestroy(&s);
>>>>>>>>>>>
>>>>>>>>>>> I will merge the necessary stuff into 'next' to make this work.
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> I have put an example without the PetscSection here:
>>>>>>>>>> https://github.com/bishesh/petscPoissonIrregular/blob/master/poissonIrregular.cxx
>>>>>>>>>> From the code you have written above, how do we let PetscSection
>>>>>>>>>> select only those cells that lie in the computational domain ? Is it that
>>>>>>>>>> for every "point" inside the above loop, we check whether it lies in the
>>>>>>>>>> domain or not, e.g using the function isPosInDomain(...) in the .cxx file I
>>>>>>>>>> linked to?
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> 1) Give me permission to comment on the source (I am 'knepley')
>>>>>>>>>
>>>>>>>>> 2) You mask out the (i,j,k) that you do not want in that loop
>>>>>>>>>
>>>>>>>>
>>>>>>>> Done.
>>>>>>>> I mask it out using isPosInDomain() function::
>>>>>>>> if(isPosInDomain(&testPoisson,i,j,k)) {
>>>>>>>> ierr = DMDAGetCellPoint(dm, i, j, k, &point);CHKERRQ(
>>>>>>>> ierr);
>>>>>>>> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); //
>>>>>>>> You know how many dof are on each vertex
>>>>>>>> }
>>>>>>>>
>>>>>>>> And please let me know when I can rebuild the 'next' branch of
>>>>>>>> petsc so that DMDAGetCellPoint can be used. Currently compiler complains as
>>>>>>>> it cannot find it.
>>>>>>>>
>>>>>>>
>>>>>>> Pushed.
>>>>>>>
>>>>>>
>>>>>> Pulled it thanks, but compiler was complaining that it didn't find
>>>>>> DMDAGetCellPoint.
>>>>>> Doing grep -R 'DMDAGetCellPoint' include/
>>>>>> returned nothing although the implementation is there in dalocal.c
>>>>>> So I added the declaration to petscdmda.h file just above the
>>>>>> declaration of DMDAGetNumCells.
>>>>>>
>>>>>
>>>>> I will add the declaration.
>>>>>
>>>>>
>>>>>> I have also added you as a collaborator in the github project
>>>>>> repository so that you can edit and add comments to the file I linked
>>>>>> above. My computeMatrix still uses the DMDA instead of the PetscSection, so
>>>>>> is it where the changes need to be made ?
>>>>>>
>>>>>
>>>>> Instead of MatSetValuesStencil(), you would go back to using
>>>>> MatSetValues(), and calculate the indices using PetscSection.
>>>>> You get the section
>>>>>
>>>>> DMGetDefaultGlobalSection(dm, &gs)
>>>>>
>>>>> and then
>>>>>
>>>>> DMDAGetCellPoint(dm, i, j, k, &p);
>>>>> PetscSectionGetOffset(gs, p, &ind);
>>>>> row = ind < 0 ? -(ind+1) : ind;
>>>>>
>>>>>
>>>>
>>>> I'm sorry but did not understood. For e.g., in 2D, for each point p, I
>>>> have 2 dof. Does PetsSectionGetOffset give me the two offsets for the dof
>>>> sth like ind[0] = 0 and ind[1] = 1 ?
>>>> How would you change the following code that sets 2D laplacian stencil
>>>> for 2 dof to use petscsection and Matsetvalues instead ?
>>>>
>>>
>>> Forget about multiple degrees of freedom until your scalar problem works.
>>>
>>>
>>>> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
>>>> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
>>>> ierr = DMGetDefaultGlobalSection(dm,&gs);CHKERRQ(ierr); /*Here I
>>>> get the PetscSection*/
>>>> for(PetscInt dof = 0; dof < 2; ++dof) {
>>>> row.c = dof;
>>>> for(PetscInt k = 0; k < 5; ++k) {
>>>> col[k].c = dof;
>>>> }
>>>>
>>>
>>> I can't imagine that this is what you want. You have two non-interacting
>>> problems.
>>>
>>
>> I just wanted a very basic example to use PetscSection with multiple
>> dofs. In this example I'm solving Lap(u) = f where Lap is a component wise
>> laplacian and u is a vector. So if I want to apply the Lap component-wise
>> independently I could solve a system thrice for fx, fy and fz; but here
>> just to demonstrate, I wanted to have all fx, fy and fz line up in a single
>> vector.
>>
>> I did as you suggeted to write the values in the matrix, thanks. But now
>> is there a way to let Petsc know automatically the size of the matrix it
>> needs to create for the linear system based on the PetsSection I created?
>> When I used only DMDA without Petscsection, this was possible, so I could
>> just use:
>> ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);
>> ierr =
>> KSPSetComputeOperators(ksp,computeMatrix2d,(void*)&testPoisson);CHKERRQ(ierr);
>> ierr =
>> KSPSetComputeRHS(ksp,computeRhs2d,(void*)&testPoisson);CHKERRQ(ierr);
>> where in my computeMatrix2d function I would set the non-zero values in
>> the matrix.
>>
>
> DMCreateMatrix(), just as before.
>
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:
//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:
[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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131031/0d83aaca/attachment-0001.html>
More information about the petsc-users
mailing list