<div dir="ltr"><div class="gmail_default" style="font-size:small">Thanks Jed, that seems like it will work.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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.  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?</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Oct 24, 2017 at 11:01 PM, Jed Brown <span dir="ltr"><<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hmm, this is a use case we would like to support, but I'm not sure how<br>
to do it.  I think at this point that your local submatrix will only<br>
have a stencil set if you were using MatNest, which I don't recommend<br>
(whatever you do, don't depend on it).  But since you have the DMDA for<br>
your J_hh block, you can call MatSetStencil() yourself.<br>
<br>
  ierr = DMDAGetGhostCorners(da,&<wbr>starts[0],&starts[1],&starts[<wbr>2],&dims[0],&dims[1],&dims[2])<wbr>;CHKERRQ(ierr);<br>
  ierr = MatSetStencil(A,dim,dims,<wbr>starts,dof);CHKERRQ(ierr);<br>
<br>
I don't know of code you could use to translate the stencil coordinates<br>
to local indices (it's just a lexicographic ordering so simple for you<br>
to write), but you can use MatSetValuesLocal on the J_hr/J_rh blocks.<br>
<div class="HOEnZb"><div class="h5"><br>
"zakaryah ." <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> writes:<br>
<br>
> Well I made a little progress by considering SNES ex28.c.  In the Jacobian<br>
> routine, I call DMCompositeGetLocalISs, then use the IS to call<br>
> MatGetLocalSubMatrix.  I call these J_rr, J_rh, J_hr, and J_hh, where r<br>
> represents the redundant variables and h represents the displacements.  I<br>
> assume I can call MatSetValuesStencil on J_hh, as before, and MatSetValue<br>
> on J_rr (which is 1x1).  I'm guessing that J_rr, J_rh, and J_hr can only be<br>
> set on the processor which owns the redundant variable - is this correct?<br>
> How do I determine the ordering for J_hr and J_rh?<br>
><br>
> On Tue, Oct 24, 2017 at 2:45 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> wrote:<br>
><br>
>> I see - I use a local variable to add up the terms on each process, then<br>
>> call MPI_Reduce within the function on process 0, which owns the redundant<br>
>> variable.<br>
>><br>
>> I have one more question - for the calculation of the Jacobian, my life is<br>
>> made much much easier by using MatSetValuesStencil.  However, the matrix<br>
>> which the SNES Jacobian receives as argument is the "full" matrix,<br>
>> containing the DMDA variables (displacements), plus the redundant<br>
>> variable.  How do I access the submatrix corresponding just to the DMDA?<br>
>> If I can do that, then I can call MatSetValuesStencil on the submatrix.  Is<br>
>> this the right approach?  I'm not sure how to set the elements of the<br>
>> Jacobian which correspond to the redundant variable, either - i.e., how do<br>
>> I get the ordering?<br>
>><br>
>><br>
</div></div></blockquote></div><br></div>