[petsc-users] help me to solve 6-dim PDE with PETSc with advice

[Matthew Knepley]

On Sun, Jul 30, 2017 at 12:44 PM, Oleksandr Koshkarov <olk548 at mail.usask.ca>
wrote:

> 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 think you are wrong here. First, I do not understand why you think the
asymptotics have changed. Second, the constant for
FFT is about 6, which would be really hard to beat with a quadrature.

Last, this discretization looks fully coupled to me. I cannot see how you
would parallelize this right off the bat. This is a problem
with spectral discretizations, and the main motivator for spectral elements.

  Thanks,

    Matt

> 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 <olk548 at mail.usask.ca
> > 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 <olk548 at mail.usask.ca>
>>>> 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)
>>>>
>>>> \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/
>
>
>


-- 
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/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170730/9b6471ad/attachment-0001.html>