<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><br></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<br><br></div><div>Timothée<br></div></div>