<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
  </head>
  <body text="#000000" bgcolor="#FFFFFF">
    <p>Thanks for the references, and warning.</p>
    <p>I am trying to do exactly this<br>
    </p>
    <p><a class="moz-txt-link-freetext" href="http://www.sciencedirect.com/science/article/pii/S0021999115004738">http://www.sciencedirect.com/science/article/pii/S0021999115004738</a></p>
    <p>(I probably should have added the reference in the first place,
      sorry)</p>
    <p>the evolution formula is (15) in the article.</p>
    <p>So, I will do convolution only in 3D (while full system is 6D)
      because velocity space is not discretized with Fourier</p>
    <p>So As I understand, it will not be N^2 vs N logN. Moreover, I
      know that algorithm with FFT has significantly larger "constant".</p>
    <p>Please correct me if I am wrong here :)</p>
    <p>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).</p>
    <p>I hope it does make sense.</p>
    <p>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>
    <p>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.<br>
    </p>
    Best regards,<br>
    Oleksandr Koshkarov.<br>
    <br>
    <div class="moz-cite-prefix">On 07/30/2017 09:01 AM, Matthew Knepley
      wrote:<br>
    </div>
    <blockquote type="cite"
cite="mid:CAMYG4GnrLwxFkVaJ6x41uhEH=bk2jjF2MLdCV=N2JMoMZOsZgQ@mail.gmail.com">
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
      <div dir="ltr">
        <div class="gmail_extra">
          <div class="gmail_quote">On Sat, Jul 29, 2017 at 7:01 PM,
            Oleksandr Koshkarov <span dir="ltr"><<a
                href="mailto:olk548@mail.usask.ca" target="_blank"
                moz-do-not-send="true">olk548@mail.usask.ca</a>></span>
            wrote:<br>
            <blockquote class="gmail_quote" style="margin:0px 0px 0px
              0.8ex;border-left:1px solid
              rgb(204,204,204);padding-left:1ex">Thank you for the
              response,<br>
              <br>
              Now that you said about additional communication for fft,
              I want to ask another question:<br>
              <br>
              What if I disregard FFT and will compute convolution
              directly, I will rewrite it as vector-matrix-vector
              product<br>
            </blockquote>
            <div><br>
            </div>
            <div> I would caution you here. The algorithm is the most
              important thing to get right up front. In this case, the
              effect of</div>
            <div>dimension will be dominant, and anything that does not
              scale well with dimension will be thrown away, perhaps</div>
            <div>quite quickly, as you want to solve bigger problems.</div>
            <div><br>
            </div>
            <div>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</div>
            <div>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</div>
            <div>equal,</div>
            <div><br>
            </div>
            <div>   100^6 (log 100^6) = x^12</div>
            <div>   1e12 (6 2) = x^12</div>
            <div>   x ~ 10</div>
            <div><br>
            </div>
            <div>Past versions of PETSc had a higher dimensional DA
              construct (adda in prior versions), but no one used it so
              we removed it.</div>
            <div><br>
            </div>
            <div>I assume you are proposing something like</div>
            <div><br>
            </div>
            <div>  <a href="https://arxiv.org/abs/1306.4625"
                moz-do-not-send="true">https://arxiv.org/abs/1306.4625</a></div>
            <div><br>
            </div>
            <div>or maybe </div>
            <div><br>
            </div>
            <div>  <a href="https://arxiv.org/abs/1610.00397"
                moz-do-not-send="true">https://arxiv.org/abs/1610.00397</a></div>
            <div><br>
            </div>
            <div>The later is interesting since I think you could use a
              DMDA for the velocity and something like DMPlex for the
              space.</div>
            <div><br>
            </div>
            <div>  Thanks,</div>
            <div><br>
            </div>
            <div>    Matt</div>
            <div> </div>
            <blockquote class="gmail_quote" style="margin:0px 0px 0px
              0.8ex;border-left:1px solid
              rgb(204,204,204);padding-left:1ex">
              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).<br>
              <br>
              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?<br>
              <br>
              Best regards,<br>
              <br>
              Alex.
              <div class="gmail-m_6801210051757406332HOEnZb">
                <div class="gmail-m_6801210051757406332h5"><br>
                  <br>
                  <br>
                  On 07/29/2017 05:23 PM, Barry Smith wrote:<br>
                  <blockquote class="gmail_quote" style="margin:0px 0px
                    0px 0.8ex;border-left:1px solid
                    rgb(204,204,204);padding-left:1ex">
                       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.<br>
                    <br>
                       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.<br>
                    <br>
                        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.<br>
                    <br>
                        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.<br>
                    <br>
                         I wish I had better answers for managing the 6d
                    domain decomposition<br>
                    <br>
                        Barry<br>
                    <br>
                    <br>
                    <br>
                    <br>
                    <blockquote class="gmail_quote" style="margin:0px
                      0px 0px 0.8ex;border-left:1px solid
                      rgb(204,204,204);padding-left:1ex">
                      On Jul 29, 2017, at 5:06 PM, Oleksandr Koshkarov
                      <<a href="mailto:olk548@mail.usask.ca"
                        target="_blank" moz-do-not-send="true">olk548@mail.usask.ca</a>>
                      wrote:<br>
                      <br>
                      Dear All,<br>
                      <br>
                      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.<br>
                      <br>
                      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:<br>
                      <br>
                      (Equations are in latex format, the nice website
                      to see them is <a
                        href="https://www.codecogs.com/latex/eqneditor.php"
                        rel="noreferrer" target="_blank"
                        moz-do-not-send="true">https://www.codecogs.com/latex<wbr>/eqneditor.php</a>)<br>
                      <br>
                      \frac{dC_{m,m,p}}{dt} =\\<br>
                      \partial_x \left ( a_n C_{n+1,m,p} +b_n C_{n,m,p}
                      +c_n C_{n-1,m,p} \right ) \\<br>
                      + \partial_y \left ( a_m C_{n,m+1,p} +b_m
                      C_{n,m,p} +c_m C_{n,m-1,p}  \right ) \\<br>
                      + \partial_z \left ( a_p C_{n,m,p+1} +b_p
                      C_{n,m,p} +c_p C_{n,m,p-1}  \right ) \\<br>
                      + 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} \\<br>
                      + 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)<br>
                      <br>
                      <br>
                      where a,b,c,d,e,f are some constants which can
                      depend on n,m,p,<br>
                      <br>
                      n,m,p = 0...N,<br>
                      <br>
                      C_{n,m,p} = C_{n,m,p}(x,y,z),<br>
                      <br>
                      E_x = E_x(x,y,z), (same for E_y,B_x,...)<br>
                      <br>
                      and fields are described with equation of the sort
                      (linear pdes with source terms dependent on C):<br>
                      <br>
                      \frac{dE_x}{dt} = \partial_y B_z - \partial_z B_x
                      + (AC_{1,0,0} + BC_{0,0,0}) \\<br>
                      \frac{dB_x}{dt} = -\partial_y E_z + \partial_z E_x<br>
                      <br>
                      with A,B being some constants.<br>
                      <br>
                      I will need fully implicit Crank–Nicolson method,
                      so my plan is to use SNES or TS where<br>
                      <br>
                      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.<br>
                      <br>
                      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?<br>
                      <br>
                      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?)<br>
                      <br>
                      If you have some other suggestions, please feel
                      free to share<br>
                      <br>
                      Thanks and best regards,<br>
                      <br>
                      Oleksandr Koshkarov.<br>
                      <br>
                    </blockquote>
                  </blockquote>
                  <br>
                </div>
              </div>
            </blockquote>
          </div>
          <br>
          <br clear="all">
          <div><br>
          </div>
          -- <br>
          <div class="gmail-m_6801210051757406332gmail_signature">
            <div dir="ltr">
              <div>What most experimenters take for granted before they
                begin their experiments is infinitely more interesting
                than any results to which their experiments lead.<br>
                -- Norbert Wiener</div>
              <div><br>
              </div>
              <div><a href="http://www.caam.rice.edu/%7Emk51/"
                  target="_blank" moz-do-not-send="true">http://www.caam.rice.edu/~<wbr>mk51/</a><br>
              </div>
            </div>
          </div>
        </div>
      </div>
    </blockquote>
    <br>
  </body>
</html>