[petsc-users] A number of questions about DMDA with SNES and Quasi-Newton methods

[Jed Brown]

"zakaryah ." <zakaryah at gmail.com> 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 <jed at jedbrown.org> 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 ." <zakaryah at gmail.com> 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 . <zakaryah at gmail.com> 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?
>> >>
>> >>
>>