On Sun, Jul 30, 2017 at 12:44 PM, Oleksandr Koshkarov <[email protected]> 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 <[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]> >>>> 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/
