<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Sep 4, 2015 at 8:20 PM, Timothée Nicolas <span dir="ltr"><<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>I see, good. Then there should be no problem. Do you also have a suggestion for questions 2. and 3. ?<br></div></div></blockquote><div><br></div><div>2. Makes little difference. Get the solver converging and scaling first.</div><div><br></div><div>3. The DMDA allocates the matrix correctly for the stencil you specify. You should not have the error</div><div> if your insertions match the specification.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div></div>Timothée<br></div><div class="gmail_extra"><br><div class="gmail_quote">2015-09-05 3:42 GMT+09:00 Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Fri, Sep 4, 2015 at 11:20 AM, Timothée Nicolas <span dir="ltr"><<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div><div><div>Hi,<br><br>(short version at the end) <br><br>I am starting to write the matrices for the preconditioner of my application. It's a physics-based preconditioner (L Chacon Journal of Comput. Phys. 2002, 2003 and Physics of Plasmas 2010 if you want to know more), so I must compute the linearized operators as matrices. I then plan to use KSP Multigrid to invert the resulting linear systems.<br></div><br></div>My question is, what is the most efficient way of setting the matrices for linearized operators in the context of a 3D grid handled with DMDA and several degrees of freedom (8 in my case) ? Take for instance a term like curl(v x B) in the induction equation, which, linearized, gives curl(v0 x . ) for the magnetic part. I think the logical approach is to use vectors of the form (0,...,0,1,0,...,0) and to compute the result of the operator, and store it in the matrix by columns using MatSetValuesStencil. The matrix is created with DMCreateMatrix.<br><br>Since all the vectors have mostly zeros, by writing by hand what parts you need to compute, the number of evaluations for each vector is only the number of nonzeros in a column, typically 21 in my case (7 point stencil * 3 directions), so the total number of evaluations of the operator should be reasonable and efficient that way.<br><br>It gives promising results, however I am having trouble handling the regions at the interface between the processors. Because if your vector has its only non zero value on a grid point which is on the boundary between two processors, then it means that one of the non-zero value in the result vector will be living on a ghost point. And to compute the operator at this ghost point, you need an *extra* ghost point, which makes the whole thing quite cumbersome. You could of course increase the stencil_width of the DMDA by one, but it does not seem the good approach.<br></div></div></div></div></div></blockquote><div><br></div></span><div>Why not? The only reason that the stencil width exists is to specify how much data from other processors</div><div>is necessary to compute the residual/Jacobian. It sounds like two ghosts are needed.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><span><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div></div>Instead, I was wondering if it would be possible to write the matrix by column vectors. So I would still create vectors of the form (0,...,0,1,0,...,0) but in Petsc Vec type, then compute the result directly in Petsc Vec, and put this vector in the matrix as a column. But somehow it does not seem really efficient, unless there is a way Petsc knows that the vector you input has mostly zeros, as well as the result vector.<br><br></div>So, in short<br><br></div>1. How would you compute a matrix representing an operator in a 3D system with several degrees of freedom in the most efficient way ?<br><br></div>And two other question related to this matter:<br><div><br></div><div>2. When I do the Multigrid inversion, I will momentarily ignore 3 degrees of freedom (the velocity). What is the more efficient between (i) putting 1 everywhere on the diagonal in the velocity part or (ii) somehow create a new DMDA with 8-3 = 5 degrees of freedom to handle the problem separately.<br><br></div><div>3. When I use DMCreateMatrix, do I have to preallocate the matrix ? I have tried, but got errors like<br>[11]PETSC ERROR: New nonzero at (296,296) caused a malloc<br>Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check<br><br></div><div>which could not be prevented by using the advised MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)<br><br></div><div>Sorry for the long question<br></div><div><br></div><div>Best<span><font color="#888888"><br><br></font></span></div><span><font color="#888888"><div>Timothée<br></div></font></span></div>
</blockquote></span></div><span><font color="#888888"><br><br clear="all"><span class="HOEnZb"><font color="#888888"><div><br></div>-- <br><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>
</font></span></font></span></div></div>
</blockquote></div><br></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">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></div>