<div dir="ltr"><div dir="ltr">On Fri, Jan 24, 2020 at 9:25 AM Mark Cunningham <<a href="mailto:mark.cunningham@ariacoustics.com">mark.cunningham@ariacoustics.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>As a novice PETSc user, can I ask for some advice? I've found some test examples that do similar sorts of things but I would appreciate any suggestions.<br></div><div>I have a matrix generated from overset meshes with a structure like</div><div>A = |A0 B01 B0n|</div><div> |B01 A1 B1n|</div><div> : :</div><div> |B0n ... An| where the Ai are finite difference stencils and the Bij</div><div> represent interpolation between the grids. <br></div><div>1. If I create a DMDA composite object that defines each of the grid sizes through repeated calls to DMCreate1D (because we are omitting blanked points from the 3D mesh) and create the matrix</div><div>through DMCreateMatrix, can I still fill the matrix through MatSetValues using the global index?</div><div>(This is how the code operates now.)<br></div></div></blockquote><div><br></div><div>I cannot understand this question as is. First, I think you mean a DMComposite object, made from many DMDA objects. I do not understand using DMDACreate1D() here. However, DMComposite does not understand how to couple grids together. Thus DMCreateMatrix() will give back a block diagonal matrix, with no allocated nonzero between grids. You can put them in yourself, _or_ you can just insert nonzeros, incurring the overhead of reallocating (which can be large).</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"><div dir="ltr"><div></div><div><div>2. If I have to declare the matrix to be matnest and I have to create each of the subblocks individually, how do I cope with some of the Bij being all zero? Must there be a MatSetValues call for each block?</div></div></div></blockquote><div><br></div><div>You can leave some blocks as NULL in MatNest.</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"><div dir="ltr"><div>3. The principle reason at the moment for the change is to use diag(A0...An) as the precondition matrix and use ILU on the blocks. If I set the preconditioner to be bjacobi and then fetch the subksps, can I then set the subpcs to be ilu? <br></div></div></blockquote><div><br></div><div>Yes, if you want block ILU(0) you can get this using</div><div><br></div><div> -pc_type bjacobi -sub_pc_type ilu</div><div><br></div><div>I am not sure that BJACOBI makes the blocks match the MatNest layout. I will have to check. If not, this is easy</div><div>to do by hand by getting out the ISes for the blocks.</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"><div dir="ltr"><div></div><div>4. A further complication is we would like to have a boundary condition where the values of the boundary points are defined by an expansion in functions that satisfy the radiation condition. So, I would like to define an index set that identifies the boundary points and define a PCSHELL on A<br></div><div>that for the pcsetup phase will do a QR factorization of the auxiliary matrix: G, where we're solving <br></div><div>Gy = QRy =c and then on pcapply do the Ry = Q* c triangular solve that will enable me to update <br></div><div>the boundary points (y subset of x the solution vector.) How do I get the PCSHELL to work with the block jacobi strategy.</div></div></blockquote><div><br></div><div>I am not sure I would do things this way, but maybe I do not understand well enough. Suppose we magically have the boundary values.</div><div>Then we can just call</div><div><br></div><div> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatZeroRowsIS.html">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatZeroRowsIS.html</a></div><div><br></div><div>with your IS to set the rows of constrained unknowns to 1 and set the boundary values in the rhs b. So, it seems like you can just calculate</div><div>the boundary values using your QR, and then insert them using the function above.</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"><div dir="ltr"><div>Thanks for any suggestions.</div><div><br></div><div>Mark Cunningham<br></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><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.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>