<div dir="ltr">I wanted to revive this thread and move it to petsc-dev. This problem seems to be harder than I realized.<div><br><div>Suppose MatMult inside KSPSolve() inside SNESSolve() cannot compute a valid output vector. </div><div>For example, <span style="font-size:13.1999998092651px">it's a MatMFFD and as part of its function evaluation it  has to evaluate an implicitly-defined </span></div><div><span style="font-size:13.1999998092651px">constitutive model (e.g., solve an equation of state) and this inner solve diverges </span></div><div><span style="font-size:13.1999998092651px">(e.g., the time step is too big).  I want to be able to abort the linear </span></div><div><span style="font-size:13.1999998092651px">solve and the nonlinear solve, return a suitable "converged" reason and let the user retry, maybe with a </span></div><div><span style="font-size:13.1999998092651px">different timestep size.  This is for a hand-rolled time stepper, but TS would face similar issues. </span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px">Based on the previous thread here </span><span style="font-size:13.1999998092651px"><a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2014-August/022597.html">http://lists.mcs.anl.gov/pipermail/petsc-users/2014-August/022597.html</a></span></div><div><span style="font-size:13.1999998092651px">I tried marking the result of MatMult as "invalid" and let it propagate up to KSPSolve() where it can be handled.</span></div><div><span style="font-size:13.1999998092651px">This quickly gets out of control, since the invalid Vec isn't returned to the KSP immediately.  It could be a work</span></div><div><span style="font-size:13.1999998092651px">vector, which is fed into PCApply() along various code paths, depending on the side of the preconditioner, whether it's a</span></div><div><span style="font-size:13.1999998092651px">transpose solve, etc.  Each of these transformations (e.g., PCApply()) would then have to check the validity of</span></div><div><span style="font-size:13.1999998092651px">the input argument, clear its error condition and set it on the output argument, etc.  Very error-prone and fragile.</span></div><div><span style="font-size:13.1999998092651px">Not to mention the large amount of code to sift through.</span></div><div><br></div><div><div><span style="font-size:13.1999998092651px">This is a general problem of exception handling -- we want to "unwind" the stack to the point where the problem should </span></div><div><span style="font-size:13.1999998092651px">be handled, but there doesn't seem to a good way to do it.  We also want to be able to clear all of the error conditions</span></div><div><span style="font-size:13.1999998092651px">on the way up (e.g., mark vectors as valid again, but not too early), otherwise we leave the solver in an invalid state.</span></div></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div>Instead of passing an exception condition up the stack I could try storing that condition in one of the more globally-visible </div><div>objects (e.g., the Mat), but if the error occurs inside the evaluation of the residual that's being differenced, it doesn't really </div><div>have access to the Mat.  This probably raises various thread safety issues as well.</div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px">Using SNESSetFunctionDomainError() doesn't seem to be a solution: a MatMFFD created with MatCreateSNESMF()</span></div><div><span style="font-size:13.1999998092651px">has a pointer to SNES, but the function evaluation code actually has no clue about that. More generally, I don't</span></div><div><span style="font-size:13.1999998092651px">know whether we want to wait for the linear solve to complete before handling this exception: it is unnecessary, </span></div><div><span style="font-size:13.1999998092651px">it might be an expensive linear solve </span><span style="font-size:13.1999998092651px">and the result of such a KSPSolve() is probably undefined and might blow up in </span></div><div><span style="font-size:13.1999998092651px">unexpected ways.  I suppose if there is a way to get a hold of SNES, each subsequent MatMult_MFFD has to check</span></div><div><span style="font-size:13.1999998092651px">whether the domain error is set and return early in that case?  We would still have to wait for the linear solve to grind</span></div><div><span style="font-size:13.1999998092651px">through the rest of its iterations.    I don't know, however, if there is a good way to guarantee </span><span style="font-size:13.1999998092651px">that linear solver will get </span></div><div><span style="font-size:13.1999998092651px">through this quickly and without unintended consequences. </span><span style="font-size:13.1999998092651px">Should MatMFFD also get a hold of the KSP and set a flag </span></div><div><span style="font-size:13.1999998092651px">there to abort?  I still don't know what the intervening code (e.g., the various PCApply()) will do before the KSP has a</span></div><div><span style="font-size:13.1999998092651px">chance to deal with this.</span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px">I'm now thinking that setting some vector entries to NaN might be a good solution: I hope </span><span style="font-size:13.1999998092651px">this NaN will propagate all the </span></div><div><span style="font-size:13.1999998092651px">way up through the subsequent arithmetic operations (does the IEEE floating-point arithmetic guarantees?), this "error</span></div><div><span style="font-size:13.1999998092651px">condition" gets automatically cleared the next time the vector is recomputed, since its values are reset.  Finally, I want</span></div><div><span style="font-size:13.1999998092651px">this exception to be detected globally but without incurring an extra reduction every time the residual is evaluated,</span></div><div><span style="font-size:13.1999998092651px">and NaN will be show up in the norm that (most) KSPs would compute anyway.  That way KSP could diverge with a</span></div><div><span style="font-size:13.1999998092651px">KSP_DIVERGED_NAN or a similar reason and the user would have an option to retry.  The problem with this approach</span></div><div><span style="font-size:13.1999998092651px">is that VecValidEntries() in MatMult() and PCApply() will throw an error before this can work, so I'm trying to think about</span></div><div><span style="font-size:13.1999998092651px">good ways of turning it off.  Any ideas about how to do this?</span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px">Incidentally, I realized that I don't understand how SNESFunctionDomainError can be handled gracefully in the current </span></div><div><span style="font-size:13.1999998092651px">set up: it's not set or checked collectively, so there isn't a good way to abort and retry across the whole comm, is there?</span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px">Dmitry.</span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><br></div><div><br></div><div><br></div><div><span style="font-size:13.1999998092651px"><br></span></div><div><br><br><div class="gmail_quote">On Sun Aug 31 2014 at 10:12:53 PM Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Dmitry Karpeyev <<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>> writes:<br>
<br>
> Handling this at the KSP level (I actually think the Mat level is more<br>
> appropriate, since the operator, not the solver, knows its domain),<br>
<br>
We are dynamically discovering the domain, but I don't think it's<br>
appropriate for Mat to refuse to evaluate any more matrix actions until<br>
some action is taken at the MatMFFD/SNESMF level.  Marking the Vec<br>
invalid is fine, but some action needs to be taken and if Derek wants<br>
the SNES to skip further evaluations, we need to propagate the<br>
information up the stack somehow.<br>
</blockquote></div></div></div></div>