Re: [petsc-users] How to create a local to global mapping and construct matrix correctly

[Barry Smith]

> On Sep 16, 2016, at 9:00 PM, Ji Zhang  wrote:
> 
> Thanks for your previous suggestion and the construction from little m to big 
> M have accomplished. 
> 
> For a MPI program, a arbitrary matrix is been shorted in different cups (i.e. 
> 3), and each cup only contain part of then. So I think the matrix have two 
> kinds of indexes, a local one indicate the location of values at the 
> corresponding cup, and the global one indicate the location at the whole 
> matrix. I would like to know the relation between them and find the way to 
> shift the index from one to another. 

   This depends on what the little matrices are that you are putting into the 
large matrix. For PDE type problems you can look at the PETSc KSP and SNES 
tutorial examples but for your problem I don't know.
> 
> I have one more question, the function VecGetArray() only return a pointer to 
> the local data array. What should I do if I need a pointer to the whole data 
> array? 

  VecScatterCreateToAll() but note this is not scalable for very large problems
> 
> Wayne
> 
> On Sat, Sep 17, 2016 at 9:00 AM, Barry Smith  wrote:
> 
> > On Sep 16, 2016, at 7:52 PM, Ji Zhang  wrote:
> >
> > Sorry. What I mean is that, for example, I have a matrix
> > [a1, a2, a3]
> > mij = [b1, b2, b3] ,
> > [c1, c2, c3]
> > and using 3 cups. Thus, mij in cpu 2 is
> > mij_2 = [b1, b2, b3] .
> >
> > The local index of element b1 is (1, 1) and it's global index is (2, 1). 
> > How can I get the global index from the local index, and local index from 
> > global index?
> 
>That is something your code needs to generate and deal with, it is not 
> something PETSc can do for you directly. You are defining the little m's and 
> the big M and deciding where to put the little m's into the big ends. 
> PETSc/we have no idea what the little m's represent in terms of the big M and 
> where they would belong, that is completely the business of your application.
> 
>Barry
> 
> 
> 
> >
> > Thanks.
> > 2016-09-17
> > Best,
> > Regards,
> > Zhang Ji
> > Beijing Computational Science Research Center
> > E-mail: got...@gmail.com
> >
> >
> >
> >
> > Wayne
> >
> > On Sat, Sep 17, 2016 at 2:24 AM, Barry Smith  wrote:
> >
> >"Gives wrong answers" is not very informative. What answer do you expect 
> > and what answer do you get?
> >
> >   Note that each process is looping over mSizes?
> >
> > for i in range(len(mSizes)):
> > for j in range(len(mSizes)):
> >
> >   Is this what you want? It doesn't seem likely that you want all processes 
> > to generate all information in the matrix. Each process should be doing a 
> > subset of the generation.
> >
> >Barry
> >
> > > On Sep 16, 2016, at 11:03 AM, Ji Zhang  wrote:
> > >
> > > Dear all,
> > >
> > > I have a number of small 'mpidense' matrices mij, and I want to construct 
> > > them to a big 'mpidense' matrix M like this:
> > >  [  m11  m12  m13  ]
> > > M =  |  m21  m22  m23  |   ,
> > >  [  m31  m32  m33  ]
> > >
> > > And a short demo is below. I'm using python, but their grammar are 
> > > similar.
> > > import numpy as np
> > > from petsc4py import PETSc
> > > import sys, petsc4py
> > >
> > >
> > > petsc4py.init(sys.argv)
> > > mSizes = (2, 2)
> > > mij = []
> > >
> > > # create sub-matrices mij
> > > for i in range(len(mSizes)):
> > > for j in range(len(mSizes)):
> > > temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > > temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
> > > temp_m.setType('mpidense')
> > > temp_m.setFromOptions()
> > > temp_m.setUp()
> > > temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
> > > temp_m.assemble()
> > > temp_m.view()
> > > mij.append(temp_m)
> > >
> > > # Now we have four sub-matrices. I would like to construct them into a 
> > > big matrix M.
> > > M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > > M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes
> > > M.setType('mpidense')
> > > M.setFromOptions()
> > > M.setUp()
> > > mLocations = np.insert(np.cumsum(mSizes), 0, 0)# mLocations = [0, 
> > > mSizes]
> > > for i in range(len(mSizes)):
> > > for j in range(len(mSizes)):
> > > temp_m = mij[i*len(mSizes)+j].getDenseArray()
> > > for k in range(temp_m.shape[0]):
> > > M.setValues(mLocations[i]+k, 
> > > np.arange(mLocations[j],mLocations[j+1],dtype='int32'), temp_m[k, :])
> > > M.assemble()
> > > M.view()
> > > The code works well in a single cup, but give wrong answer for 2 and more 
> > > cores.
> > >
> > > Thanks.
> > > 2016-09-17
> > > Best,
> > > Regards,
> > > Zhang Ji
> > > Beijing Computational Science Research Center
> > > E-mail: got...@gmail.com
> > >
> > >
> >
> >
> 
> 

Re: [petsc-users] How to create a local to global mapping and construct matrix correctly

[Ji Zhang]

Thanks for your previous suggestion and the construction from little m to
big M have accomplished.

For a MPI program, a arbitrary matrix is been shorted in different cups
(i.e. 3), and each cup only contain part of then. So I think the matrix
have two kinds of indexes, a local one indicate the location of values at
the corresponding cup, and the global one indicate the location at the
whole matrix. I would like to know the relation between them and find the
way to shift the index from one to another.

I have one more question, the function VecGetArray() only return a pointer
to the local data array. What should I do if I need a pointer to the whole
data array?

Wayne

On Sat, Sep 17, 2016 at 9:00 AM, Barry Smith  wrote:

>
> > On Sep 16, 2016, at 7:52 PM, Ji Zhang  wrote:
> >
> > Sorry. What I mean is that, for example, I have a matrix
> > [a1, a2, a3]
> > mij = [b1, b2, b3] ,
> > [c1, c2, c3]
> > and using 3 cups. Thus, mij in cpu 2 is
> > mij_2 = [b1, b2, b3] .
> >
> > The local index of element b1 is (1, 1) and it's global index is (2, 1).
> How can I get the global index from the local index, and local index from
> global index?
>
>That is something your code needs to generate and deal with, it is not
> something PETSc can do for you directly. You are defining the little m's
> and the big M and deciding where to put the little m's into the big ends.
> PETSc/we have no idea what the little m's represent in terms of the big M
> and where they would belong, that is completely the business of your
> application.
>
>Barry
>
>
>
> >
> > Thanks.
> > 2016-09-17
> > Best,
> > Regards,
> > Zhang Ji
> > Beijing Computational Science Research Center
> > E-mail: got...@gmail.com
> >
> >
> >
> >
> > Wayne
> >
> > On Sat, Sep 17, 2016 at 2:24 AM, Barry Smith  wrote:
> >
> >"Gives wrong answers" is not very informative. What answer do you
> expect and what answer do you get?
> >
> >   Note that each process is looping over mSizes?
> >
> > for i in range(len(mSizes)):
> > for j in range(len(mSizes)):
> >
> >   Is this what you want? It doesn't seem likely that you want all
> processes to generate all information in the matrix. Each process should be
> doing a subset of the generation.
> >
> >Barry
> >
> > > On Sep 16, 2016, at 11:03 AM, Ji Zhang  wrote:
> > >
> > > Dear all,
> > >
> > > I have a number of small 'mpidense' matrices mij, and I want to
> construct them to a big 'mpidense' matrix M like this:
> > >  [  m11  m12  m13  ]
> > > M =  |  m21  m22  m23  |   ,
> > >  [  m31  m32  m33  ]
> > >
> > > And a short demo is below. I'm using python, but their grammar are
> similar.
> > > import numpy as np
> > > from petsc4py import PETSc
> > > import sys, petsc4py
> > >
> > >
> > > petsc4py.init(sys.argv)
> > > mSizes = (2, 2)
> > > mij = []
> > >
> > > # create sub-matrices mij
> > > for i in range(len(mSizes)):
> > > for j in range(len(mSizes)):
> > > temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > > temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
> > > temp_m.setType('mpidense')
> > > temp_m.setFromOptions()
> > > temp_m.setUp()
> > > temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
> > > temp_m.assemble()
> > > temp_m.view()
> > > mij.append(temp_m)
> > >
> > > # Now we have four sub-matrices. I would like to construct them into a
> big matrix M.
> > > M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > > M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes
> > > M.setType('mpidense')
> > > M.setFromOptions()
> > > M.setUp()
> > > mLocations = np.insert(np.cumsum(mSizes), 0, 0)# mLocations = [0,
> mSizes]
> > > for i in range(len(mSizes)):
> > > for j in range(len(mSizes)):
> > > temp_m = mij[i*len(mSizes)+j].getDenseArray()
> > > for k in range(temp_m.shape[0]):
> > > M.setValues(mLocations[i]+k, np.arange(mLocations[j],
> mLocations[j+1],dtype='int32'), temp_m[k, :])
> > > M.assemble()
> > > M.view()
> > > The code works well in a single cup, but give wrong answer for 2 and
> more cores.
> > >
> > > Thanks.
> > > 2016-09-17
> > > Best,
> > > Regards,
> > > Zhang Ji
> > > Beijing Computational Science Research Center
> > > E-mail: got...@gmail.com
> > >
> > >
> >
> >
>
>

Re: [petsc-users] How to create a local to global mapping and construct matrix correctly

[Barry Smith]

> On Sep 16, 2016, at 7:52 PM, Ji Zhang  wrote:
> 
> Sorry. What I mean is that, for example, I have a matrix 
> [a1, a2, a3]
> mij = [b1, b2, b3] , 
> [c1, c2, c3]
> and using 3 cups. Thus, mij in cpu 2 is 
> mij_2 = [b1, b2, b3] .
> 
> The local index of element b1 is (1, 1) and it's global index is (2, 1). How 
> can I get the global index from the local index, and local index from global 
> index? 

   That is something your code needs to generate and deal with, it is not 
something PETSc can do for you directly. You are defining the little m's and 
the big M and deciding where to put the little m's into the big ends. PETSc/we 
have no idea what the little m's represent in terms of the big M and where they 
would belong, that is completely the business of your application.

   Barry



> 
> Thanks. 
> 2016-09-17
> Best, 
> Regards,
> Zhang Ji 
> Beijing Computational Science Research Center 
> E-mail: got...@gmail.com
> 
> 
> 
> 
> Wayne
> 
> On Sat, Sep 17, 2016 at 2:24 AM, Barry Smith  wrote:
> 
>"Gives wrong answers" is not very informative. What answer do you expect 
> and what answer do you get?
> 
>   Note that each process is looping over mSizes?
> 
> for i in range(len(mSizes)):
> for j in range(len(mSizes)):
> 
>   Is this what you want? It doesn't seem likely that you want all processes 
> to generate all information in the matrix. Each process should be doing a 
> subset of the generation.
> 
>Barry
> 
> > On Sep 16, 2016, at 11:03 AM, Ji Zhang  wrote:
> >
> > Dear all,
> >
> > I have a number of small 'mpidense' matrices mij, and I want to construct 
> > them to a big 'mpidense' matrix M like this:
> >  [  m11  m12  m13  ]
> > M =  |  m21  m22  m23  |   ,
> >  [  m31  m32  m33  ]
> >
> > And a short demo is below. I'm using python, but their grammar are similar.
> > import numpy as np
> > from petsc4py import PETSc
> > import sys, petsc4py
> >
> >
> > petsc4py.init(sys.argv)
> > mSizes = (2, 2)
> > mij = []
> >
> > # create sub-matrices mij
> > for i in range(len(mSizes)):
> > for j in range(len(mSizes)):
> > temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
> > temp_m.setType('mpidense')
> > temp_m.setFromOptions()
> > temp_m.setUp()
> > temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
> > temp_m.assemble()
> > temp_m.view()
> > mij.append(temp_m)
> >
> > # Now we have four sub-matrices. I would like to construct them into a big 
> > matrix M.
> > M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes
> > M.setType('mpidense')
> > M.setFromOptions()
> > M.setUp()
> > mLocations = np.insert(np.cumsum(mSizes), 0, 0)# mLocations = [0, 
> > mSizes]
> > for i in range(len(mSizes)):
> > for j in range(len(mSizes)):
> > temp_m = mij[i*len(mSizes)+j].getDenseArray()
> > for k in range(temp_m.shape[0]):
> > M.setValues(mLocations[i]+k, 
> > np.arange(mLocations[j],mLocations[j+1],dtype='int32'), temp_m[k, :])
> > M.assemble()
> > M.view()
> > The code works well in a single cup, but give wrong answer for 2 and more 
> > cores.
> >
> > Thanks.
> > 2016-09-17
> > Best,
> > Regards,
> > Zhang Ji
> > Beijing Computational Science Research Center
> > E-mail: got...@gmail.com
> >
> >
> 
> 

Re: [petsc-users] How to create a local to global mapping and construct matrix correctly

[Ji Zhang]

Sorry. What I mean is that, for example, I have a matrix
[a1, a2, a3]
mij = [b1, b2, b3] ,
[c1, c2, c3]
and using 3 cups. Thus, mij in cpu 2 is
mij_2 = [b1, b2, b3] .

The local index of element b1 is (1, 1) and it's global index is (2, 1).
How can I get the global index from the local index, and local index from
global index?

Thanks.
2016-09-17
Best,
Regards,
Zhang Ji
Beijing Computational Science Research Center
E-mail: got...@gmail.com




Wayne

On Sat, Sep 17, 2016 at 2:24 AM, Barry Smith  wrote:

>
>"Gives wrong answers" is not very informative. What answer do you
> expect and what answer do you get?
>
>   Note that each process is looping over mSizes?
>
> for i in range(len(mSizes)):
> for j in range(len(mSizes)):
>
>   Is this what you want? It doesn't seem likely that you want all
> processes to generate all information in the matrix. Each process should be
> doing a subset of the generation.
>
>Barry
>
> > On Sep 16, 2016, at 11:03 AM, Ji Zhang  wrote:
> >
> > Dear all,
> >
> > I have a number of small 'mpidense' matrices mij, and I want to
> construct them to a big 'mpidense' matrix M like this:
> >  [  m11  m12  m13  ]
> > M =  |  m21  m22  m23  |   ,
> >  [  m31  m32  m33  ]
> >
> > And a short demo is below. I'm using python, but their grammar are
> similar.
> > import numpy as np
> > from petsc4py import PETSc
> > import sys, petsc4py
> >
> >
> > petsc4py.init(sys.argv)
> > mSizes = (2, 2)
> > mij = []
> >
> > # create sub-matrices mij
> > for i in range(len(mSizes)):
> > for j in range(len(mSizes)):
> > temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
> > temp_m.setType('mpidense')
> > temp_m.setFromOptions()
> > temp_m.setUp()
> > temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
> > temp_m.assemble()
> > temp_m.view()
> > mij.append(temp_m)
> >
> > # Now we have four sub-matrices. I would like to construct them into a
> big matrix M.
> > M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> > M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes
> > M.setType('mpidense')
> > M.setFromOptions()
> > M.setUp()
> > mLocations = np.insert(np.cumsum(mSizes), 0, 0)# mLocations = [0,
> mSizes]
> > for i in range(len(mSizes)):
> > for j in range(len(mSizes)):
> > temp_m = mij[i*len(mSizes)+j].getDenseArray()
> > for k in range(temp_m.shape[0]):
> > M.setValues(mLocations[i]+k, np.arange(mLocations[j],
> mLocations[j+1],dtype='int32'), temp_m[k, :])
> > M.assemble()
> > M.view()
> > The code works well in a single cup, but give wrong answer for 2 and
> more cores.
> >
> > Thanks.
> > 2016-09-17
> > Best,
> > Regards,
> > Zhang Ji
> > Beijing Computational Science Research Center
> > E-mail: got...@gmail.com
> >
> >
>
>

Re: [petsc-users] fieldsplit preconditioner for indefinite matrix

[Barry Smith]

> On Sep 16, 2016, at 6:09 PM, Hoang Giang Bui  wrote:
> 
> Hi Barry
> 
> You are right, using MatCreateAIJ() eliminates the first issue. Previously I 
> ran the mpi code with one process so A,B,C,D is all MPIAIJ
> 
> And how about the second issue, this error will always be thrown if A11 is 
> nonzero, which is my case?
> 
> Nevertheless, I would like to report my simple finding: I changed the part 
> around line 552 to 

  I'm sorry what file are you talking about? What version of PETSc? What other 
lines of  code are around 552? I can't figure out where you are doing this.

   Barry

> 
> if (D) {
> ierr = MatAXPY(*S, -1.0, D, SUBSET_NONZERO_PATTERN);CHKERRQ(ierr);
>  }
> 
> I could get ex42 works with
> 
> ierr = KSPSetOperators(ksp_S,A,A);CHKERRQ(ierr);
> 
> parameters:
> mpirun -np 1 ex42 \
> -stokes_ksp_monitor \
> -stokes_ksp_type fgmres \
> -stokes_pc_type fieldsplit \
> -stokes_pc_fieldsplit_type schur \
> -stokes_pc_fieldsplit_schur_fact_type full \
> -stokes_pc_fieldsplit_schur_precondition full \
> -stokes_fieldsplit_u_ksp_type preonly \
> -stokes_fieldsplit_u_pc_type lu \
> -stokes_fieldsplit_u_pc_factor_mat_solver_package mumps \
> -stokes_fieldsplit_p_ksp_type gmres \
> -stokes_fieldsplit_p_ksp_monitor_true_residual \
> -stokes_fieldsplit_p_ksp_max_it 300 \
> -stokes_fieldsplit_p_ksp_rtol 1.0e-12 \
> -stokes_fieldsplit_p_ksp_gmres_restart 300 \
> -stokes_fieldsplit_p_ksp_gmres_modifiedgramschmidt \
> -stokes_fieldsplit_p_pc_type lu \
> -stokes_fieldsplit_p_pc_factor_mat_solver_package mumps \
> 
> Output:
>   Residual norms for stokes_ solve.
>   0 KSP Residual norm 1.327791371202e-02 
> Residual norms for stokes_fieldsplit_p_ solve.
> 0 KSP preconditioned resid norm 1.651372938841e+02 true resid norm 
> 5.775755720828e-02 ||r(i)||/||b|| 1.e+00
> 1 KSP preconditioned resid norm 1.172753353368e+00 true resid norm 
> 2.072348962892e-05 ||r(i)||/||b|| 3.588013522487e-04
> 2 KSP preconditioned resid norm 3.931379526610e-13 true resid norm 
> 1.878299731917e-16 ||r(i)||/||b|| 3.252041503665e-15
>   1 KSP Residual norm 3.385960118582e-17 
> 
> inner convergence is much better although 2 iterations (:-( ??
> 
> I also obtain the same convergence behavior for the problem with A11!=0
> 
> Please suggest if this makes sense, or I did something wrong.
> 
> Giang
> 
> On Fri, Sep 16, 2016 at 8:31 PM, Barry Smith  wrote:
> 
>Why is your C matrix an MPIAIJ matrix on one process? In general we 
> recommend creating a SeqAIJ matrix for one process and MPIAIJ for multiple. 
> You can use MatCreateAIJ() and it will always create the correct one.
> 
>We could change the code as you suggest but I want to make sure that is 
> the best solution in your case.
> 
>Barry
> 
> 
> 
> > On Sep 16, 2016, at 3:31 AM, Hoang Giang Bui  wrote:
> >
> > Hi Matt
> >
> > I believed at line 523, src/ksp/ksp/utils/schurm.c
> >
> > ierr = MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
> >
> > in my test case C is MPIAIJ and AinvB is SEQAIJ, hence it throws the error.
> >
> > In fact I guess there are two issues with it
> > line 521, ierr = MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, 
> > );CHKERRQ(ierr);
> > shall we convert this to type of C matrix to ensure compatibility ?
> >
> > line 552, if(norm > PETSC_MACHINE_EPSILON) 
> > SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet 
> > implemented for Schur complements with non-vanishing D");
> > with this the Schur complement with A11!=0 will be aborted
> >
> > Giang
> >
> > On Thu, Sep 15, 2016 at 4:28 PM, Matthew Knepley  wrote:
> > On Thu, Sep 15, 2016 at 9:07 AM, Hoang Giang Bui  wrote:
> > Hi Matt
> >
> > Thanks for the comment. After looking carefully into the manual again, the 
> > key take away is that with selfp there is no option to compute the exact 
> > Schur, there are only two options to approximate the inv(A00) for selfp, 
> > which are lump and diag (diag by default). I misunderstood this previously.
> >
> > There is online manual entry mentioned about PC_FIELDSPLIT_SCHUR_PRE_FULL, 
> > which is not documented elsewhere in the offline manual. I tried to access 
> > that by setting
> > -pc_fieldsplit_schur_precondition full
> >
> > Yep, I wrote that specifically for testing, but its very slow so I did not 
> > document it to prevent people from complaining.
> >
> > but it gives the error
> >
> > [0]PETSC ERROR: - Error Message 
> > --
> > [0]PETSC ERROR: Arguments are incompatible
> > [0]PETSC ERROR: MatMatMult requires A, mpiaij, to be compatible with B, 
> > seqaij
> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
> > trouble shooting.
> > [0]PETSC ERROR: Petsc Release Version 3.7.3, Jul, 24, 2016
> > [0]PETSC ERROR: python on a arch-linux2-c-opt named bermuda by hbui Thu Sep 

Re: [petsc-users] Question about memory usage in Multigrid preconditioner

[Jed Brown]

Barry Smith  writes:

>> On Sep 15, 2016, at 1:10 PM, Dave May  wrote:
>> 
>> 
>> 
>> On Thursday, 15 September 2016, Barry Smith  wrote:
>> 
>>Should we have some simple selection of default algorithms based on 
>> problem size/number of processes? For example if using more than 1000 
>> processes then use scalable version etc?  How would we decide on the 
>> parameter values?
>> 
>> I don't like the idea of having "smart" selection by default as it's 
>> terribly annoying for the user when they try and understand the performance 
>> characteristics of a given method when they do a strong/weak scaling test. 
>> If such a smart selection strategy was adopted, the details of it should be 
>> made abundantly clear to the user.
>> 
>> These algs are dependent on many some factors, thus making the smart 
>> selection for all use cases hard / impossible.
>> 
>> I would be happy with unifying the three inplemtationa with three different 
>> options AND having these implantation options documented in the man page. 
>> Maybe even the man page should advise users which to use in particular 
>> circumstances (I think there is something similar on the VecScatter page).
>> 
>> I have these as suggestions for unifying the options names using bools 
>> 
>> -matptap_explicit_transpose
>> -matptap_symbolic_transpose_dense
>> -matptap_symbolic_transpose
>> 
>> Or maybe enums is more clear
>> -matptap_impl {explicit_pt,symbolic_pt_dense,symbolic_pt}
>> 
>> which are equivalent to these options
>> 1) the current default
>> 2) -matrap 0
>> 3) -matrap 0 -matptap_scalable
>> 
>> Maybe there could be a fourth option
>> -matptap_dynamic_selection
>> which chooses the most appropriate alg given machine info, problem size, 
>> partition size, At least if the user explicitly chooses the 
>> dynamic_selection mode, they wouldn't be surprised if there were any bumps 
>> appearing in any scaling study they conducted.
>
>I like the idea of enum types with the final enum type being "dynamically 
> select one for me".

I also like enums and "-matptap_impl auto" (which could be the default).


signature.asc
Description: PGP signature

Re: [petsc-users] fieldsplit preconditioner for indefinite matrix

[Barry Smith]

   Why is your C matrix an MPIAIJ matrix on one process? In general we 
recommend creating a SeqAIJ matrix for one process and MPIAIJ for multiple. You 
can use MatCreateAIJ() and it will always create the correct one.

   We could change the code as you suggest but I want to make sure that is the 
best solution in your case.

   Barry



> On Sep 16, 2016, at 3:31 AM, Hoang Giang Bui  wrote:
> 
> Hi Matt
> 
> I believed at line 523, src/ksp/ksp/utils/schurm.c
> 
> ierr = MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
> 
> in my test case C is MPIAIJ and AinvB is SEQAIJ, hence it throws the error.
> 
> In fact I guess there are two issues with it
> line 521, ierr = MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, 
> );CHKERRQ(ierr);
> shall we convert this to type of C matrix to ensure compatibility ?
> 
> line 552, if(norm > PETSC_MACHINE_EPSILON) 
> SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented 
> for Schur complements with non-vanishing D");
> with this the Schur complement with A11!=0 will be aborted
> 
> Giang
> 
> On Thu, Sep 15, 2016 at 4:28 PM, Matthew Knepley  wrote:
> On Thu, Sep 15, 2016 at 9:07 AM, Hoang Giang Bui  wrote:
> Hi Matt
> 
> Thanks for the comment. After looking carefully into the manual again, the 
> key take away is that with selfp there is no option to compute the exact 
> Schur, there are only two options to approximate the inv(A00) for selfp, 
> which are lump and diag (diag by default). I misunderstood this previously.
> 
> There is online manual entry mentioned about PC_FIELDSPLIT_SCHUR_PRE_FULL, 
> which is not documented elsewhere in the offline manual. I tried to access 
> that by setting
> -pc_fieldsplit_schur_precondition full
> 
> Yep, I wrote that specifically for testing, but its very slow so I did not 
> document it to prevent people from complaining.
>  
> but it gives the error
> 
> [0]PETSC ERROR: - Error Message 
> --
> [0]PETSC ERROR: Arguments are incompatible
> [0]PETSC ERROR: MatMatMult requires A, mpiaij, to be compatible with B, seqaij
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
> trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.7.3, Jul, 24, 2016 
> [0]PETSC ERROR: python on a arch-linux2-c-opt named bermuda by hbui Thu Sep 
> 15 15:46:56 2016
> [0]PETSC ERROR: Configure options --with-shared-libraries --with-debugging=0 
> --with-pic --download-fblaslapack=yes --download-suitesparse 
> --download-ptscotch=yes --download-metis=yes --download-parmetis=yes 
> --download-scalapack=yes --download-mumps=yes --download-hypre=yes 
> --download-ml=yes --download-pastix=yes --with-mpi-dir=/opt/openmpi-1.10.1 
> --prefix=/home/hbui/opt/petsc-3.7.3
> [0]PETSC ERROR: #1 MatMatMult() line 9514 in 
> /home/hbui/sw/petsc-3.7.3/src/mat/interface/matrix.c
> [0]PETSC ERROR: #2 MatSchurComplementComputeExplicitOperator() line 526 in 
> /home/hbui/sw/petsc-3.7.3/src/ksp/ksp/utils/schurm.c
> [0]PETSC ERROR: #3 PCSetUp_FieldSplit() line 792 in 
> /home/hbui/sw/petsc-3.7.3/src/ksp/pc/impls/fieldsplit/fieldsplit.c
> [0]PETSC ERROR: #4 PCSetUp() line 968 in 
> /home/hbui/sw/petsc-3.7.3/src/ksp/pc/interface/precon.c
> [0]PETSC ERROR: #5 KSPSetUp() line 390 in 
> /home/hbui/sw/petsc-3.7.3/src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: #6 KSPSolve() line 599 in 
> /home/hbui/sw/petsc-3.7.3/src/ksp/ksp/interface/itfunc.c
> 
> Please excuse me to insist on forming the exact Schur complement, but as you 
> said, I would like to track down what creates problem in my code by starting 
> from a very exact but ineffective solution.
> 
> Sure, I understand. I do not understand how A can be MPI and B can be Seq. Do 
> you know how that happens?
> 
>   Thanks,
> 
>  Matt
>  
> Giang
> 
> On Thu, Sep 15, 2016 at 2:56 PM, Matthew Knepley  wrote:
> On Thu, Sep 15, 2016 at 4:11 AM, Hoang Giang Bui  wrote:
> Dear Barry
> 
> Thanks for the clarification. I got exactly what you said if the code changed 
> to
> ierr = KSPSetOperators(ksp_S,B,B);CHKERRQ(ierr);
>   Residual norms for stokes_ solve.
>   0 KSP Residual norm 1.327791371202e-02 
> Residual norms for stokes_fieldsplit_p_ solve.
> 0 KSP preconditioned resid norm 0.e+00 true resid norm 
> 0.e+00 ||r(i)||/||b||   -nan
>   1 KSP Residual norm 3.997711925708e-17 
> 
> but I guess we solve a different problem if B is used for the linear system.
> 
> in addition, changed to
> ierr = KSPSetOperators(ksp_S,A,A);CHKERRQ(ierr);
> also works but inner iteration converged not in one iteration
> 
>   Residual norms for stokes_ solve.
>   0 KSP Residual norm 1.327791371202e-02 
> Residual norms for stokes_fieldsplit_p_ solve.
> 0 KSP preconditioned resid norm 5.308049264070e+02 true resid norm 
> 5.775755720828e-02 ||r(i)||/||b|| 

Re: [petsc-users] How to create a local to global mapping and construct matrix correctly

[Barry Smith]

   "Gives wrong answers" is not very informative. What answer do you expect and 
what answer do you get?

  Note that each process is looping over mSizes?

for i in range(len(mSizes)):
for j in range(len(mSizes)):

  Is this what you want? It doesn't seem likely that you want all processes to 
generate all information in the matrix. Each process should be doing a subset 
of the generation.

   Barry

> On Sep 16, 2016, at 11:03 AM, Ji Zhang  wrote:
> 
> Dear all, 
> 
> I have a number of small 'mpidense' matrices mij, and I want to construct 
> them to a big 'mpidense' matrix M like this:
>  [  m11  m12  m13  ]
> M =  |  m21  m22  m23  |   ,
>  [  m31  m32  m33  ]
> 
> And a short demo is below. I'm using python, but their grammar are similar. 
> import numpy as np
> from petsc4py import PETSc
> import sys, petsc4py
> 
> 
> petsc4py.init(sys.argv)
> mSizes = (2, 2)
> mij = []
> 
> # create sub-matrices mij
> for i in range(len(mSizes)):
> for j in range(len(mSizes)):
> temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
> temp_m.setType('mpidense')
> temp_m.setFromOptions()
> temp_m.setUp()
> temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
> temp_m.assemble()
> temp_m.view()
> mij.append(temp_m)
> 
> # Now we have four sub-matrices. I would like to construct them into a big 
> matrix M.
> M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes
> M.setType('mpidense')
> M.setFromOptions()
> M.setUp()
> mLocations = np.insert(np.cumsum(mSizes), 0, 0)# mLocations = [0, mSizes]
> for i in range(len(mSizes)):
> for j in range(len(mSizes)):
> temp_m = mij[i*len(mSizes)+j].getDenseArray()
> for k in range(temp_m.shape[0]):
> M.setValues(mLocations[i]+k, 
> np.arange(mLocations[j],mLocations[j+1],dtype='int32'), temp_m[k, :])
> M.assemble()
> M.view()
> The code works well in a single cup, but give wrong answer for 2 and more 
> cores. 
> 
> Thanks. 
> 2016-09-17
> Best, 
> Regards,
> Zhang Ji 
> Beijing Computational Science Research Center 
> E-mail: got...@gmail.com
> 
> 

Re: [petsc-users] Question about PETScSF usage in DMNetwork/DMPlex

[Adrian Maldonado]

Just one addition about one thing I've noticed.

The section:

PetscSection Object: 2 MPI processes
  type not yet set
Process 0:
  (   0) dim  1 offset   0
  (   1) dim  1 offset   1
  (   2) dim  1 offset   2
  (   3) dim  1 offset   3
  (   4) dim -2 offset  -8
  (   5) dim -2 offset  -9
  (   6) dim -2 offset -10
Process 1:
  (   0) dim  1 offset   4
  (   1) dim  1 offset   5
  (   2) dim  1 offset   6
  (   3) dim  1 offset   7
  (   4) dim  1 offset   8
  (   5) dim  1 offset   9

For the ghost values 4, 5, 6... is encoding the ghost values as rank = -(-2
+ 1) and offset = -(-8 + 1) ?

On Fri, Sep 16, 2016 at 11:36 AM, Adrian Maldonado 
wrote:

> Hi,
>
> I am trying to understand some of the data structures DMPlex/DMNetwork
> creates and the relationship among them.
>
> As an example, I have an small test circuit (/src/ksp/ksp/examples/
> tutorials/network/ex1.c).
>
> This is a graph that consists on 6 edges and 4 vertices, each one of those
> having one degree of freedom.  When ran with two processors, each rank will
> own 3 edges. Rank 0 will own one vertex (3 ghost) and Rank 1 will own 3
> vertices.
>
> These are some data structures for this problem. I am getting these data
> structures inside DMNetworkDistribute
> 
>
> DM Object: Parallel Mesh 2 MPI processes
>   type: plex
> Parallel Mesh in 1 dimensions:
>   0-cells: 4 3
>   1-cells: 3 3
> Labels:
>   depth: 2 strata of sizes (4, 3)
>
> This, as  I understand, is printing a tree with all the vertices and edges
> in each processor (owned and ghost).
>
> PetscSection Object: 2 MPI processes
>   type not yet set
> Process 0:
>   (   0) dim  1 offset   0
>   (   1) dim  1 offset   1
>   (   2) dim  1 offset   2
>   (   3) dim  1 offset   3
>   (   4) dim -2 offset  -8
>   (   5) dim -2 offset  -9
>   (   6) dim -2 offset -10
> Process 1:
>   (   0) dim  1 offset   4
>   (   1) dim  1 offset   5
>   (   2) dim  1 offset   6
>   (   3) dim  1 offset   7
>   (   4) dim  1 offset   8
>   (   5) dim  1 offset   9
>
> This is a global PETSc section that gives me the global numbering for the
> owned points and (garbage?) negative values for ghost.
>
> Until here everything is good. But then I print the PetscSF that is
> created by  'DMPlexDistribute'. This I do not understand:
>
> PetscSF Object: Migration SF 2 MPI processes
>   type: basic
> sort=rank-order
>   [0] Number of roots=10, leaves=7, remote ranks=1
>   [0] 0 <- (0,0)
>   [0] 1 <- (0,1)
>   [0] 2 <- (0,3)
>   [0] 3 <- (0,6)
>   [0] 4 <- (0,7)
>   [0] 5 <- (0,8)
>   [0] 6 <- (0,9)
>   [1] Number of roots=0, leaves=6, remote ranks=1
>   [1] 0 <- (0,2)
>   [1] 1 <- (0,4)
>   [1] 2 <- (0,5)
>   [1] 3 <- (0,7)
>   [1] 4 <- (0,8)
>   [1] 5 <- (0,9)
>   [0] Roots referenced by my leaves, by rank
>   [0] 0: 7 edges
>   [0]0 <- 0
>   [0]1 <- 1
>   [0]2 <- 3
>   [0]3 <- 6
>   [0]4 <- 7
>   [0]5 <- 8
>   [0]6 <- 9
>   [1] Roots referenced by my leaves, by rank
>   [1] 0: 6 edges
>   [1]0 <- 2
>   [1]1 <- 4
>   [1]2 <- 5
>   [1]3 <- 7
>   [1]4 <- 8
>   [1]5 <- 9
>
> I understand that SF is a data structure that saves references to pieces
> of data that are now owned by the process (https://arxiv.org/pdf/1506.
> 06194v1.pdf, page 4).
>
> Since the only ghost nodes appear in rank 0 (three ghost vertices) I would
> expect something like:
> *rank 0:*
>   4 - (1, 3)  (to read: point 4 is owned by rank 1 and is rank's 1 point 3)
>   etc...
> *rank 1:*
>   nothing
>
> Is my intuition correct? If so, what does the star forest that I get from
> DMPlexDistribute mean? I am printing the wrong thing?
>
> Thank you
>
> --
> D. Adrian Maldonado, PhD Candidate
> Electrical & Computer Engineering Dept.
> Illinois Institute of Technology
> 3301 S. Dearborn Street, Chicago, IL 60616
>



-- 
D. Adrian Maldonado, PhD Candidate
Electrical & Computer Engineering Dept.
Illinois Institute of Technology
3301 S. Dearborn Street, Chicago, IL 60616

Re: [petsc-users] Question about memory usage in Multigrid preconditioner

[Hengjie Wang]

Hi Dave,

I add both options and test it by solving the poisson eqn in  a 1024 
cube with 32^3 cores. This test used to give the OOM error. Now it runs 
well.

I attach the ksp_view and log_view's output in case you want to know.
I also test my original code with those petsc options by simulating a 
decaying turbulence in a 1024 cube. It also works.  I am going to test 
the code on a larger scale. If there is any problem then, I will let you 
know.

This really helps a lot. Thank you so much.

Regards,
Frank


On 9/15/2016 3:35 AM, Dave May wrote:

HI all,

I the only unexpected memory usage I can see is associated with the 
call to MatPtAP().

Here is something you can try immediately.
Run your code with the additional options
  -matrap 0 -matptap_scalable

I didn't realize this before, but the default behaviour of MatPtAP in 
parallel is actually to to explicitly form the transpose of P (e.g. 
assemble R = P^T) and then compute R.A.P.

You don't want to do this. The option -matrap 0 resolves this issue.

The implementation of P^T.A.P has two variants.
The scalable implementation (with respect to memory usage) is selected 
via the second option -matptap_scalable.


Try it out - I see a significant memory reduction using these options 
for particular mesh sizes / partitions.


I've attached a cleaned up version of the code you sent me.
There were a number of memory leaks and other issues.
The main points being
  * You should call DMDAVecGetArrayF90() before VecAssembly{Begin,End}
  * You should call PetscFinalize(), otherwise the option -log_summary 
(-log_view) will not display anything once the program has completed.



Thanks,
  Dave


On 15 September 2016 at 08:03, Hengjie Wang > wrote:


Hi Dave,

Sorry, I should have put more comment to explain the code.
The number of process in each dimension is the same: Px = Py=Pz=P.
So is the domain size.
So if the you want to run the code for a  512^3 grid points on
16^3 cores, you need to set "-N 512 -P 16" in the command line.
I add more comments and also fix an error in the attached code. (
The error only effects the accuracy of solution but not the memory
usage. )

Thank you.
Frank


On 9/14/2016 9:05 PM, Dave May wrote:



On Thursday, 15 September 2016, Dave May > wrote:



On Thursday, 15 September 2016, frank  wrote:

Hi,

I write a simple code to re-produce the error. I hope
this can help to diagnose the problem.
The code just solves a 3d poisson equation.


Why is the stencil width a runtime parameter?? And why is the
default value 2? For 7-pnt FD Laplace, you only need
a stencil width of 1.

Was this choice made to mimic something in the
real application code?


Please ignore - I misunderstood your usage of the param set by -P


I run the code on a 1024^3 mesh. The process partition is
32 * 32 * 32. That's when I re-produce the OOM error.
Each core has about 2G memory.
I also run the code on a 512^3 mesh with 16 * 16 * 16
processes. The ksp solver works fine.
I attached the code, ksp_view_pre's output and my petsc
option file.

Thank you.
Frank

On 09/09/2016 06:38 PM, Hengjie Wang wrote:

Hi Barry,

I checked. On the supercomputer, I had the option
"-ksp_view_pre" but it is not in file I sent you. I am
sorry for the confusion.

Regards,
Frank

On Friday, September 9, 2016, Barry Smith
 wrote:


> On Sep 9, 2016, at 3:11 PM, frank
 wrote:
>
> Hi Barry,
>
> I think the first KSP view output is from
-ksp_view_pre. Before I submitted the test, I was
not sure whether there would be OOM error or not. So
I added both -ksp_view_pre and -ksp_view.

  But the options file you sent specifically does
NOT list the -ksp_view_pre so how could it be from that?

   Sorry to be pedantic but I've spent too much time
in the past trying to debug from incorrect
information and want to make sure that the
information I have is correct before thinking.
Please recheck exactly what happened. Rerun with the
exact input file you emailed if that is needed.

   Barry

>
> Frank
>
>
> On 09/09/2016 12:38 PM, Barry Smith wrote:
>>   Why does ksp_view2.txt have two KSP views in it
  

Re: [petsc-users] 2D vector in 3D dmda

[Barry Smith]

> On Sep 16, 2016, at 9:29 AM, Aurelien PONTE  wrote:
> 
> Hi,
> 
> I've started using petsc4py in order to solve a 3D problem (inversion of 
> elliptic operator).
> I would like to store 2D metric terms describing the grid


 What do you mean by 2D metric terms describing the grid?

 Do you want to store something like a little dense 2d array for each grid 
point? If so create another 3D DA with a dof = the product of the dimensions of 
the little dense 2d array and then store the little dense 2d arrays at in a 
global vector obtained from that DA. 

Or is the grid uniform in one dimension and not uniform in the other two 
and hence you want to store the information about the non-uniformness in only a 
2d array so as to not "waste" the redundant information in the third direction? 
Then I recommend just "waste" the redundant information in the third dimension; 
it is trivial compared to all the data you need to solve the problem.

Or do you mean something else?

  Barry
   
> I am working on but don't know
> how to do that given my domain is tiled in 3D directions:
> 
> self.da = PETSc.DMDA().create([self.grid.Nx, self.grid.Ny, self.grid.Nz],
>  stencil_width=2)
> 
> I create my 3D vectors with, for example:
> 
> self.Q = self.da.createGlobalVec()
> 
> What am I supposed to do for a 2D vector?
> Is it a bad idea?
> 
> thanks
> 
> aurelien
> 

[petsc-users] Question about PETScSF usage in DMNetwork/DMPlex

[Adrian Maldonado]

Hi,

I am trying to understand some of the data structures DMPlex/DMNetwork
creates and the relationship among them.

As an example, I have an small test circuit
(/src/ksp/ksp/examples/tutorials/network/ex1.c).

This is a graph that consists on 6 edges and 4 vertices, each one of those
having one degree of freedom.  When ran with two processors, each rank will
own 3 edges. Rank 0 will own one vertex (3 ghost) and Rank 1 will own 3
vertices.

These are some data structures for this problem. I am getting these data
structures inside DMNetworkDistribute


DM Object: Parallel Mesh 2 MPI processes
  type: plex
Parallel Mesh in 1 dimensions:
  0-cells: 4 3
  1-cells: 3 3
Labels:
  depth: 2 strata of sizes (4, 3)

This, as  I understand, is printing a tree with all the vertices and edges
in each processor (owned and ghost).

PetscSection Object: 2 MPI processes
  type not yet set
Process 0:
  (   0) dim  1 offset   0
  (   1) dim  1 offset   1
  (   2) dim  1 offset   2
  (   3) dim  1 offset   3
  (   4) dim -2 offset  -8
  (   5) dim -2 offset  -9
  (   6) dim -2 offset -10
Process 1:
  (   0) dim  1 offset   4
  (   1) dim  1 offset   5
  (   2) dim  1 offset   6
  (   3) dim  1 offset   7
  (   4) dim  1 offset   8
  (   5) dim  1 offset   9

This is a global PETSc section that gives me the global numbering for the
owned points and (garbage?) negative values for ghost.

Until here everything is good. But then I print the PetscSF that is created
by  'DMPlexDistribute'. This I do not understand:

PetscSF Object: Migration SF 2 MPI processes
  type: basic
sort=rank-order
  [0] Number of roots=10, leaves=7, remote ranks=1
  [0] 0 <- (0,0)
  [0] 1 <- (0,1)
  [0] 2 <- (0,3)
  [0] 3 <- (0,6)
  [0] 4 <- (0,7)
  [0] 5 <- (0,8)
  [0] 6 <- (0,9)
  [1] Number of roots=0, leaves=6, remote ranks=1
  [1] 0 <- (0,2)
  [1] 1 <- (0,4)
  [1] 2 <- (0,5)
  [1] 3 <- (0,7)
  [1] 4 <- (0,8)
  [1] 5 <- (0,9)
  [0] Roots referenced by my leaves, by rank
  [0] 0: 7 edges
  [0]0 <- 0
  [0]1 <- 1
  [0]2 <- 3
  [0]3 <- 6
  [0]4 <- 7
  [0]5 <- 8
  [0]6 <- 9
  [1] Roots referenced by my leaves, by rank
  [1] 0: 6 edges
  [1]0 <- 2
  [1]1 <- 4
  [1]2 <- 5
  [1]3 <- 7
  [1]4 <- 8
  [1]5 <- 9

I understand that SF is a data structure that saves references to pieces of
data that are now owned by the process (
https://arxiv.org/pdf/1506.06194v1.pdf, page 4).

Since the only ghost nodes appear in rank 0 (three ghost vertices) I would
expect something like:
*rank 0:*
  4 - (1, 3)  (to read: point 4 is owned by rank 1 and is rank's 1 point 3)
  etc...
*rank 1:*
  nothing

Is my intuition correct? If so, what does the star forest that I get from
DMPlexDistribute mean? I am printing the wrong thing?

Thank you

-- 
D. Adrian Maldonado, PhD Candidate
Electrical & Computer Engineering Dept.
Illinois Institute of Technology
3301 S. Dearborn Street, Chicago, IL 60616

[petsc-users] How to create a local to global mapping and construct matrix correctly

[Ji Zhang]

Dear all,

I have a number of small 'mpidense' matrices mij, and I want to construct
them to a big 'mpidense' matrix M like this:

 [  m11  m12  m13  ]
M =  |  m21  m22  m23  |   ,
 [  m31  m32  m33  ]

And a short demo is below. I'm using python, but their grammar are similar.

import numpy as np
from petsc4py import PETSc
import sys, petsc4py


petsc4py.init(sys.argv)
mSizes = (2, 2)
mij = []

# create sub-matrices mij
for i in range(len(mSizes)):
for j in range(len(mSizes)):
temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
temp_m.setType('mpidense')
temp_m.setFromOptions()
temp_m.setUp()
temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
temp_m.assemble()
temp_m.view()
mij.append(temp_m)

# Now we have four sub-matrices. I would like to construct them into a
big matrix M.
M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes
M.setType('mpidense')
M.setFromOptions()
M.setUp()
mLocations = np.insert(np.cumsum(mSizes), 0, 0)# mLocations = [0, mSizes]
for i in range(len(mSizes)):
for j in range(len(mSizes)):
temp_m = mij[i*len(mSizes)+j].getDenseArray()
for k in range(temp_m.shape[0]):
M.setValues(mLocations[i]+k,
np.arange(mLocations[j],mLocations[j+1],dtype='int32'), temp_m[k, :])
M.assemble()
M.view()

The code works well in a single cup, but give wrong answer for 2 and more
cores.

Thanks.
2016-09-17
Best,
Regards,
Zhang Ji
Beijing Computational Science Research Center
E-mail: got...@gmail.com

[petsc-users] 2D vector in 3D dmda

[Aurelien PONTE]

Hi,

I've started using petsc4py in order to solve a 3D problem (inversion of 
elliptic operator).
I would like to store 2D metric terms describing the grid I am working 
on but don't know

how to do that given my domain is tiled in 3D directions:

self.da = PETSc.DMDA().create([self.grid.Nx, self.grid.Ny, self.grid.Nz],
  stencil_width=2)

I create my 3D vectors with, for example:

self.Q = self.da.createGlobalVec()

What am I supposed to do for a 2D vector?
Is it a bad idea?

thanks

aurelien