<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sat, Feb 25, 2017 at 10:12 AM, Gideon Simpson <span dir="ltr"><<a href="mailto:gideon.simpson@gmail.com" target="_blank">gideon.simpson@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><br><div>
<span class="m_4443323829976622494Apple-style-span" style="border-collapse:separate;color:rgb(0,0,0);font-family:Helvetica;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:normal;text-align:-webkit-auto;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">-gideon</span>

</div>
<br><div><blockquote type="cite"><div>On Feb 25, 2017, at 10:44 AM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:</div><br class="m_4443323829976622494Apple-interchange-newline"><div><div>Gideon Simpson <<a href="mailto:gideon.simpson@gmail.com" target="_blank">gideon.simpson@gmail.com</a>> writes:<br><br><blockquote type="cite">I’ve been continuing working on implementing a projection method problem, which, loosely, looks like the following.  Vec y contains the state vector for my system, y’ = f(y) which is solved with a TS, using, for now, rk4.<br><br>I have added to this a TSPostStep which does root finding on the nonlinear problem<br><br>f(lambda) = f(lambda; y, w) =  g(y + lambda * w) - g0<br></blockquote><br>You might want to treat this as a DAE.<br></div></div></blockquote><div><br></div>Yea, we discussed DAE previously.  For the sake of a reviewer, I’m trying to implement a classical method (explicit time step + projection), rather than DAE.</div><div><br><blockquote type="cite"><div><div><br><blockquote type="cite">y and w are vectors that are the size of the system (numbers of mesh points), and lambda is a scalar perturbation.  <br><br>Right now, the snes function for solving this amounts to:<br><br>PetscErrorCode projector(SNES snes, Vec lam, Vec f, void *ctx){<br><br>  AppCtx * user = (AppCtx *) ctx;<br>  const PetscScalar *lam_value;<br>  PetscScalar *f_value;<br><br>  VecGetArrayRead(lam, &lam_value);<br>  VecGetArray(f, &f_value);<br>  VecCopy(user->y, user->w);<br>  VecAXPY(user->w, lam_value[0], user->gradMy);<br><br>  f_value[0] = g(user->w) -user->g0;<br>  VecRestoreArrayRead(lam, &lam_value);<br>  VecRestoreArray(f, &f_value);<br><br>  return 0;<br>}<br><br>To get this all to work, I constructed the SNES and Vec lam with PETSC_COMM_SELF, effectively putting a copy of the nonlinear problem on each process.  Basically, the nonlinear problem is low dimensional, but it parametrically depends on the high dimensional, distributed, vectors y and w.<br><br>The one thing that bothers me about this is:  <br><br>1.  It would seem that each process is is doing all of these vector operations, which is entirely redundant, as the only thing that really needs to be shared amongst the processes is the value of lambda.  Is that correct? <br></blockquote><br>Only the SNES part is redundant, but it's super cheap because it's a<br>scalar problem.  The Vec operations are using global (parallel,<br>distributed) vectors, so they are not redundant.  Of course it is<br>critical for the SNES on every process to be configured identically and<br>deterministic so that the processes don't diverge.  And g(user->w) must<br>return the same value on each process (it probably finishes with an<br>MPI_Allreduce or similar).<br></div></div></blockquote><div><br></div><div>So within the root finding problem, when it calls, for instance, VecCopy (which is acting on distributed vectors), it’s not doing that once on process 0, once on process 1, once on process 2, etc.  Or it is, but you’re saying it’s too cheap to matter? </div></div></div></blockquote><div><br></div><div>No they are collective. That is why Jed said its crucial that the same branches are taken on all processes.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><div><div>Two things I’m mildly concerned about is that, since user->w changes at each step of the solver, by virtue of the lambda scalar changing,  if the SNES’s are running asynchronously, by virtue of them being PETSC_COMM_SELF constructs, couldn’t there be a collision?  </div></div></div></blockquote><div><br></div><div>Yes, as I say above. I think its much safer to run a COMM_WORLD SNES that just has nothing on most procs than to do the former</div><div>since you are already synchronizing everywhere.</div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><div><div>The operation g(user->w) does indeed conclude with an MPI_Allreduce, followed by DMDAVecRestoreArray and DMRestoreLocalVector.</div><div><br></div><div>Lastly, in the TSPostStep, it runs as:</div><div><div>  </div><div>SNESSolve(user->snes, NULL,user->lam);</div><div>VecGetArray(user->lam, &lam_value);</div><div>VecAXPY(y, lam_value[0], user->gradMy);</div><div>VecRestoreArray(user->lam, &lam_value);</div><div><br></div><div>Is there any communication issue with ensuring that all the SNES’s have finished on each process, before proceeding to do the vector operations?</div></div></div></div></blockquote><div><br></div><div>No.</div><div> <br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><div><div><div>I keep wondering if I should just write this as a global, with the data structures stored on process 0, just to avoid this kind of headache. </div></div></div></div></blockquote><div><br></div><div>Yes.</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><div><blockquote type="cite"><div><div><br><blockquote type="cite">2.  Is there an obvious way around this redundancy?<br><br>-gideon<br></blockquote></div></div></blockquote></div><br></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature">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></div>