On Fri, Jan 24, 2020 at 9:25 AM Mark Cunningham < [email protected]> wrote:
> 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. > I have a matrix generated from overset meshes with a structure like > A = |A0 B01 B0n| > |B01 A1 B1n| > : : > |B0n ... An| where the Ai are finite difference stencils > and the Bij > represent interpolation between the grids. > 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 > through DMCreateMatrix, can I still fill the matrix through MatSetValues > using the global index? > (This is how the code operates now.) > 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). > 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? > You can leave some blocks as NULL in MatNest. > 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? > Yes, if you want block ILU(0) you can get this using -pc_type bjacobi -sub_pc_type ilu I am not sure that BJACOBI makes the blocks match the MatNest layout. I will have to check. If not, this is easy to do by hand by getting out the ISes for the blocks. > 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 > that for the pcsetup phase will do a QR factorization of the auxiliary > matrix: G, where we're solving > Gy = QRy =c and then on pcapply do the Ry = Q* c triangular solve that > will enable me to update > the boundary points (y subset of x the solution vector.) How do I get the > PCSHELL to work with the block jacobi strategy. > 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. Then we can just call https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatZeroRowsIS.html 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 the boundary values using your QR, and then insert them using the function above. Thanks, Matt > Thanks for any suggestions. > > Mark Cunningham > -- 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 https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
