"zakaryah ." <[email protected]> writes: > Thanks Jed, that seems like it will work. > > I still have a problem, which I suspect is pretty simple to solve. To > initialize the solve, I need to run the SNES once without the control > parameter. My composite DM is dac, the DMDA inside dac is dah, the > redundant field is dab, and I have the vectors set up for the solve and the > residual. But I don't know how to setup the matrices for the Jacobian, in > the main loop, i.e. in the scope from which SNESSolve is called. I think I > have the matrix set up properly for the composite, and the evaluation > routine for the composite system extracts the submatrices as we discussed. > But I'm not sure how to setup a global matrix for the first solve, which is > in a smaller space. I've gotten a lot out of SNES ex28 but there the > smaller solve and the full solve never run in the same launch.
Why not use a different SNES? Or just make that extra equation trivial u_r = 0? > I guess I can create two matrices for the Jacobians, one like > DMCreateMatrix(dac,&J) for the full system and > DMCreateMatrix(dah,&J_hh) for the initial solve. Is there a more > efficient way - can I just create J and somehow extract the global > submatrix J_hh? > > On Tue, Oct 24, 2017 at 11:01 PM, Jed Brown <[email protected]> wrote: > >> Hmm, this is a use case we would like to support, but I'm not sure how >> to do it. I think at this point that your local submatrix will only >> have a stencil set if you were using MatNest, which I don't recommend >> (whatever you do, don't depend on it). But since you have the DMDA for >> your J_hh block, you can call MatSetStencil() yourself. >> >> ierr = DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[ >> 2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr); >> ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr); >> >> I don't know of code you could use to translate the stencil coordinates >> to local indices (it's just a lexicographic ordering so simple for you >> to write), but you can use MatSetValuesLocal on the J_hr/J_rh blocks. >> >> "zakaryah ." <[email protected]> writes: >> >> > Well I made a little progress by considering SNES ex28.c. In the >> Jacobian >> > routine, I call DMCompositeGetLocalISs, then use the IS to call >> > MatGetLocalSubMatrix. I call these J_rr, J_rh, J_hr, and J_hh, where r >> > represents the redundant variables and h represents the displacements. I >> > assume I can call MatSetValuesStencil on J_hh, as before, and MatSetValue >> > on J_rr (which is 1x1). I'm guessing that J_rr, J_rh, and J_hr can only >> be >> > set on the processor which owns the redundant variable - is this correct? >> > How do I determine the ordering for J_hr and J_rh? >> > >> > On Tue, Oct 24, 2017 at 2:45 PM, zakaryah . <[email protected]> wrote: >> > >> >> I see - I use a local variable to add up the terms on each process, then >> >> call MPI_Reduce within the function on process 0, which owns the >> redundant >> >> variable. >> >> >> >> I have one more question - for the calculation of the Jacobian, my life >> is >> >> made much much easier by using MatSetValuesStencil. However, the matrix >> >> which the SNES Jacobian receives as argument is the "full" matrix, >> >> containing the DMDA variables (displacements), plus the redundant >> >> variable. How do I access the submatrix corresponding just to the DMDA? >> >> If I can do that, then I can call MatSetValuesStencil on the >> submatrix. Is >> >> this the right approach? I'm not sure how to set the elements of the >> >> Jacobian which correspond to the redundant variable, either - i.e., how >> do >> >> I get the ordering? >> >> >> >> >>
