Thanks for the references, and warning.

I am trying to do exactly this

http://www.sciencedirect.com/science/article/pii/S0021999115004738

(I probably should have added the reference in the first place, sorry)

the evolution formula is (15) in the article.

So, I will do convolution only in 3D (while full system is 6D) because velocity space is not discretized with Fourier

So As I understand, it will not be N^2 vs N logN. Moreover, I know that algorithm with FFT has significantly larger "constant".

Please correct me if I am wrong here :)

I know for sure that FFT is a must at the end, but I thought I can start with matrix product to get more familiar with PETSc as it will be simpler to implement (At least I think so).

I hope it does make sense.

Another question about convolution via matrix product. Is there better way to do it other then to construct NxN Toeplitz-like matrix each time step on the "fly" from state vector and then multiply it to state vector?

p.s. Probably I should do FFT from the start, but I am a little bit intimidated by the problem taking into account that I do not have PETSc experience.

Best regards,
Oleksandr Koshkarov.

On 07/30/2017 09:01 AM, Matthew Knepley wrote:
On Sat, Jul 29, 2017 at 7:01 PM, Oleksandr Koshkarov <[email protected] <mailto:[email protected]>> wrote:

    Thank you for the response,

    Now that you said about additional communication for fft, I want
    to ask another question:

    What if I disregard FFT and will compute convolution directly, I
    will rewrite it as vector-matrix-vector product


I would caution you here. The algorithm is the most important thing to get right up front. In this case, the effect of dimension will be dominant, and anything that does not scale well with dimension will be thrown away, perhaps
quite quickly, as you want to solve bigger problems.

Small calculation: Take a 6-dim box which is 'n' on a side, so that N = n^6. FFT is N log N, whereas direct evaluation of a convolution is N^2. So if we could do n = 100 points on a side with FFT, you can do, assuming the constants are about
equal,

   100^6 (log 100^6) = x^12
   1e12 (6 2) = x^12
   x ~ 10

Past versions of PETSc had a higher dimensional DA construct (adda in prior versions), but no one used it so we removed it.

I assume you are proposing something like

https://arxiv.org/abs/1306.4625

or maybe

https://arxiv.org/abs/1610.00397

The later is interesting since I think you could use a DMDA for the velocity and something like DMPlex for the space.

  Thanks,

    Matt

    x^T A x (where T - is transpose, x - is my state vector, and A is
    a matrix which represents convolution and probably it will be not
    very sparse).

    Yes - I will do more work (alghorithm will be slower), but will I
    win if we take the parallelization into account - less
    communication + probably more scalable algorithm?

    Best regards,

    Alex.



    On 07/29/2017 05:23 PM, Barry Smith wrote:

           DMDA provides a simple way to do domain decomposition of
        structured grid meshes in the x, y, and z direction, that is
        it makes it easy to efficiently chop up vectors in all 3
        dimensions, allowing very large problems easily. Unfortunately
        it only goes up to 3 dimensions and attempts to produce
        versions for higher dimensions have failed.

           If you would like to do domain decomposition across all six
        of your dimensions (which I think you really need to do to
        solve large problems) we don't have a tool that helps with
        doing the domain decomposition. To make matters more
        complicated, likely the 3d mpi-fftw that need to be applied
        require rejiggering the data across the nodes to get it into a
        form where the ffts can be applied efficiently.

            I don't think there is necessarily anything theoretically
        difficult about managing the six dimensional domain
        decomposition I just don't know of any open source software
        out there that helps to do this. Maybe some PETSc users are
        aware of such software and can chime in.

            Aside from these two issues :-) PETSc would be useful for
        you since you could use our time integrators and nonlinear
        solvers. Pulling out subvectors to process on can be done with
        with either VecScatter or VecStrideScatter, VecStrideGather
        etc depending on the layout of the unknowns, i.e. how the
        domain decomposition is done.

             I wish I had better answers for managing the 6d domain
        decomposition

            Barry




            On Jul 29, 2017, at 5:06 PM, Oleksandr Koshkarov
            <[email protected] <mailto:[email protected]>> wrote:

            Dear All,

            I am a new PETSc user and I am still in the middle of the
            manual (I have finally settled to choose PETSc as my main
            numerical library) so I am sorry beforehand if my
            questions would be naive.

            I am trying to solve 6+1 dimensional Vlasov equation with
            spectral methods. More precisely, I will try to solve
            half-discretized equations of the form (approximate form)
            with pseudospectral Fourier method:

            (Equations are in latex format, the nice website to see
            them is https://www.codecogs.com/latex/eqneditor.php
            <https://www.codecogs.com/latex/eqneditor.php>)

            \frac{dC_{m,m,p}}{dt} =\\
            \partial_x \left ( a_n C_{n+1,m,p} +b_n C_{n,m,p} +c_n
            C_{n-1,m,p} \right ) \\
            + \partial_y \left ( a_m C_{n,m+1,p} +b_m C_{n,m,p} +c_m
            C_{n,m-1,p}  \right ) \\
            + \partial_z \left ( a_p C_{n,m,p+1} +b_p C_{n,m,p} +c_p
            C_{n,m,p-1}  \right ) \\
            + d_n E_x C_{n-1,m,p} + d_m E_x C_{n,m-1,p} + d_p E_x
            C_{n,m,p-1} \\
            + B_x (e_{m,p} C_{n,m-1,p-1} + f_{m,p}C_{n,m-1,p+1} +
            \dots) + B_y (\dots) + B_z (\dots)


            where a,b,c,d,e,f are some constants which can depend on
            n,m,p,

            n,m,p = 0...N,

            C_{n,m,p} = C_{n,m,p}(x,y,z),

            E_x = E_x(x,y,z), (same for E_y,B_x,...)

            and fields are described with equation of the sort (linear
            pdes with source terms dependent on C):

            \frac{dE_x}{dt} = \partial_y B_z - \partial_z B_x +
            (AC_{1,0,0} + BC_{0,0,0}) \\
            \frac{dB_x}{dt} = -\partial_y E_z + \partial_z E_x

            with A,B being some constants.

            I will need fully implicit Crank–Nicolson method, so my
            plan is to use SNES or TS where

            I will use one MPI PETSc vector which describes the state
            of the system (all, C, E, B),  then I can evolve linear
            part with just matrix multiplication.

            The first question is, should I use something like DMDA?
            if so, it is only 3dimensional but I need 6 dimensional
            vectots? Will it be faster then matrix multiplication?

            Then to do nonlinear part I will need 3d mpi-fftw to
            compute convolutions. The problem is, how do I extract
            subvectors from full big vector state? (what is the best
            approach?)

            If you have some other suggestions, please feel free to share

            Thanks and best regards,

            Oleksandr Koshkarov.





--
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

http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/%7Emk51/>

Reply via email to