[petsc-users] PETSc for ISPH

[Matthew Knepley]

On Mon, Nov 28, 2016 at 2:18 PM, Massoud Rezavand <
Massoud.Rezavand at uibk.ac.at> wrote:

> Hi,
>
> Thanks.
>
> As you know, in SPH method, the calculations are done over the neighboring
> particles (j) that fall inside a support domain defined by a circle over
> the particle of interest (i). Since the Lagrangian nature the method, the
> number of neighboring particles are varying slightly over time, e.g. in a
> 2D domain this number is varying between 43 to 51 (in my experience).
>
> The number of nonzeros per row (A_ij) is equal to the number of
> neighboring particles and normally is not fixed over time, therefore, we
> put the elements dynamically at each time step and we have to calculate
> d_nz and o_nz at each time iteration.
>
> In order to preallocate the matrix, another way would be to calculate the
> number of neighboring particles and set that as the number of nonzeros per
> row.  Doing so, do you recommend to use :
>
> MatMPIAIJSetPreallocation()
>
> to preallocate A to achieve the best performance?
>

First, you should just do it in two passes. First count the nonzeros,
preallocate the matrix, then put in the values. I have never seen this
get even close to 1% of the total runtime with an implicit method. However,
if you measure it and its significant, then I would switch to
a MatShell where you just calculate the action of the operator.

   Matt


>
> Regards,
>
> Massoud
>
>
> On 11/28/2016 06:36 PM, Barry Smith wrote:
>
>> On Nov 28, 2016, at 10:30 AM, Massoud Rezavand <
>>> Massoud.Rezavand at uibk.ac.at> wrote:
>>>
>>> Dear Barry,
>>>
>>> You recommended me to directly use MatSetValues() and not to put the
>>> matrix in a parallel CSR matrix.
>>>
>>> In order to count the d_nz and o_nz I have to put the entries into a
>>> sequential CSR matrix
>>>
>>     If you don't know the number of nonzeros per row how are you going to
>> put the values into a sequential CSR format?
>> On the other hand if you can figure out the number of nonzeros per row
>> without creating the matrix how come you cannot figure out the d_nz and
>> o_nz?
>>
>>
>> and then do the MatMPIAIJSetPreallocation() and then do the MatSet
>>> Values().
>>>
>>      If you do put the values into a sequential CSR format, which it is
>> not clear to me is needed, then you can just call
>> MatCreateMPIAIJWithArrays() and skip the "MatMPIAIJSetPreallocation() and
>> then do the MatSet Values()"
>>
>>     Barry
>>
>>
>>
>> Does it effect the performance ?
>>>
>>>
>>> Regards,
>>>
>>> Massoud
>>>
>>> On 11/21/2016 08:10 PM, Barry Smith wrote:
>>>
>>>> On Nov 21, 2016, at 12:20 PM, Massoud Rezavand <
>>>>> Massoud.Rezavand at uibk.ac.at> wrote:
>>>>>
>>>>> Thank you very much.
>>>>>
>>>>> Yes I am developing the new 3D version in Parallel with the PPM (the
>>>>> new generation OpenFPM, not released yet) library which generates the
>>>>> particles and decomposes the domain.
>>>>>
>>>>> I don't have the parallel matrix generation yet. In the old version I
>>>>> had CSR format and a vector of knowns (b).
>>>>> So, should I use MatSetValuesStencil() ?
>>>>>
>>>>      MatSetValuesStencil is for finite differences on a structured
>>>> grid. I don't think it makes sense for your application.
>>>>
>>>>      You need to use MatMPIAIJSetPreallocation() and then
>>>> MatSetValues() to put the entries in.
>>>>
>>>> What do you recommend for creating the vector of knowns (b)?
>>>>>
>>>>     Just use VecCreateMPI()
>>>>
>>>>> On the other hand, due to the convergence issues for millions of
>>>>> particles in ISPH, I have to use a preconditioner. In a paper I saw they
>>>>> have used BoomerAMG from HYPRE. Do you have any recommendation?
>>>>>
>>>>     We have many to try, it is not clear that any would be particularly
>>>> good for SPH. Certainly try BoomerAMG
>>>>
>>>> I saw an example ( ex19.c) using BoomerAMG. Should I follow that?
>>>>>
>>>>>
>>>>> PS: regarding the unbalance sparsity in SPH, yes in contrast to the
>>>>> mesh-based methods, the A matrix in ISPH is changing over the time but the
>>>>> number of non-zeros is defined by the number of neighboring particles which
>>>>> in most cases is constant.
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Massoud
>>>>>
>>>>>
>>>>>
>>>>> On 11/21/2016 06:18 PM, Barry Smith wrote:
>>>>>
>>>>>> On Nov 21, 2016, at 10:33 AM, Massoud Rezavand <
>>>>>>> Massoud.Rezavand at uibk.ac.at>
>>>>>>>   wrote:
>>>>>>>
>>>>>>> Dear all,
>>>>>>>
>>>>>>> I am going to use PETSc in an Incompressible SPH code to solve the
>>>>>>> pressure Poisson equation as a linear system of equations.
>>>>>>>
>>>>>>> In my old sequential code I used the PCG method or the BiCGSTAB with
>>>>>>> jacobi preconditioner.
>>>>>>> I used to store the coefficient matrix (A) in CSR (AIJ) format and
>>>>>>> solve it.
>>>>>>>
>>>>>>> My question is that how should I import the CSR metrix and the known
>>>>>>> vector (b) into the library to be solved? Is there an example to show how
>>>>>>> to import and external system of eqs. into PETSc?
>>>>>>>
>>>>>>>     For sequential code it is straightforward.
>>>>>>
>>>>>>     If you already have the matrix in CSR format you can call
>>>>>> MatCreateSeqAIJWithArrays() to create the PETSc matrix without copying the
>>>>>> data. You can use VecCreateSeqWithArray() to provide the vector. Or you can
>>>>>> use VecPlaceArray() to use the array of vector values you provide.
>>>>>>
>>>>>>
>>>>>> In my case, the computational domain is decomposed by another
>>>>>>> library, so does it effect the performance of PETSc?
>>>>>>>
>>>>>>>     I read this to mean you want the new code to be parallel (but
>>>>>> the old one is sequential)?
>>>>>>
>>>>>>     If you don't currently have matrix generation in parallel I urge
>>>>>> you strongly to directly use MatSetValues() to generate your matrix, do not
>>>>>> first put the matrix entries into some kind of parallel CSR format. If you
>>>>>> already have the matrix in "parallel" CSR format you can use
>>>>>> MatCreateMPIAIJWithArrays() to copy the matrix over to CSR format.
>>>>>>
>>>>>>     It is my understanding that SPH can produce matrices with very
>>>>>> unbalance sparsity. It is important to take this into account if you wish
>>>>>> to run in parallel since if you end up with some processes having many more
>>>>>> nonzeros than other processes you will get very poor performance.
>>>>>>
>>>>>>
>>>>>>     Barry
>>>>>>
>>>>>>
>>>>>>
>>>>>> With the best regards,
>>>>>>> Massoud
>>>>>>>
>>>>>>>
>>>>>>>
>


-- 
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/20161128/b7bdcb3d/attachment-0001.html>