February 2012 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Jed Brown 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 14:23, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > That's why encapsulating the related functions and their linearizations, > one (R,J) pair per DM, is helpful. Okay, but those DM interfaces are actually different. So I have a bunch of DMs running around, but how do I figure out which functions I need to implement to use with TS? Since the arguments are different, I can't reuse the same function pointers/dispatch/registration.

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Jed Brown 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 14:20, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > Just out of curiosity: how many transient problems have we seen where the > transient residual doesn't split like R(U,U_t) = R_0(U_t) + R_1(U)? There > is the issue of explicit dependence on time, but if that could be > addressed, the problem would naturally split and, as a side benefit, there > would be a natural "mass" matrix coming out of the linearization of R_0. > Moving mesh problems often have terms involving both U_t and U. The IMEX interface is g(t,X,X_t) = f(t,X) which lets us split the implicit and explicit parts. Having a separate mass matrix is normally quite wasteful and it's confusing for people doing matrix-free methods where it is natural to keep both together.

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Dmitry Karpeev 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 2:21 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Tue, Feb 28, 2012 at 14:16, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > >> Okay. I would still refrain from calling these things by the class of >> the object consuming them. DMFormFunction() and DMFormIFunction() seem >> fine, since these have a meaning independent of TS (although the DM has to >> then encode a transient problem). >> > > The IFunction is a particular way of representing a (part of) a transient > problem. Its arguments include the time and the time derivative. We have to > make that clear somehow and I think naming it after TS is helpful. > Otherwise I'll eventually just have a bunch of names (function, residual, > rhsfunction, objective, gradient, jacobian, ijacobian, hessian) and not > know which ones go together. > That's why encapsulating the related functions and their linearizations, one (R,J) pair per DM, is helpful. Dmitry. > > >> But I don't like DMGetSNESFunction() any more than I like >> DMGetKSPJacobian() and DMGetPCJacobian(). >> >> By the way, should we start using Residual instead of Function? >> > > Yes. >

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Jed Brown 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 14:16, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > Okay. I would still refrain from calling these things by the class of the > object consuming them. DMFormFunction() and DMFormIFunction() seem fine, > since these have a meaning independent of TS (although the DM has to then > encode a transient problem). > The IFunction is a particular way of representing a (part of) a transient problem. Its arguments include the time and the time derivative. We have to make that clear somehow and I think naming it after TS is helpful. Otherwise I'll eventually just have a bunch of names (function, residual, rhsfunction, objective, gradient, jacobian, ijacobian, hessian) and not know which ones go together. > But I don't like DMGetSNESFunction() any more than I like > DMGetKSPJacobian() and DMGetPCJacobian(). > > By the way, should we start using Residual instead of Function? > Yes.

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Dmitry Karpeev 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 2:06 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Tue, Feb 28, 2012 at 13:52, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > >> I don't understand the need for multiple function evaluation routines and >> their explicit naming. >> I think the situation here is completely analogous to the situation with >> matrices: I thought we agreed that each DM would be capable of generating a >> single matrix, the "Jacobian", rather than multiple matrices like J and >> J_pre, and that the different necessary matrices (again, e.g., J and J_pre) >> would instead be produced by their respective DMs -- the KSP DM and the PC >> DM, so to speak. These subDMs would be pulled out of the SNES DM and >> passed to the corresponding KSP and PC objects. Thus, a DM encapsulates a >> single "problem", and if we want to get a different problem -- the "exact" >> linearization of the original, for example -- we extract the KSP DM; >> likewise, PC DM is an approximate linearization, and so on. >> >> I don't think TS is any different -- we can have the TSI DM (maybe under >> a better name), which could have DMFormFunction encapsulating the IFunction >> (and, incidentally, a linearization encapsulating the mass matrix), etc. >> > > I'm concerned that this becomes insanity. A single TS might need > RHSFunction, IFunction, and IJacobian. They might also partition terms > differently for the functions and the Jacobian. Having different DMs seems > awfully confusing. > > In short -- I'm for each DM having one function (DMFormFunction or, >> better, DMFormResidual) and one matrix (DMFormJacobian) -- *a* >> linearization of the residual. All other functions and matrices are >> obtained from the corresponding subDMs. >> > > The interface is actually different for these (compare TSIFunction, a > transient residual, to SNESFunction) > Just out of curiosity: how many transient problems have we seen where the transient residual doesn't split like R(U,U_t) = R_0(U_t) + R_1(U)? There is the issue of explicit dependence on time, but if that could be addressed, the problem would naturally split and, as a side benefit, there would be a natural "mass" matrix coming out of the linearization of R_0. Dmitry.

1 0

Re: [petsc-dev] OpenMP/Vec
by Jed Brown 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 10:41, Lawrence Mitchell <lawrence.mitchell(a)ed.ac.uk > wrote: > It's possible that if one asks the gcc folk nicely enough that the > overhead might go to zero. Alternately, if it turns out that overheads are > unavoidable, and we've gone down the macro-ification route anyway (see > below) we could explicitly make a runtime decision based on numthreads in > PETSc. So that: > > PetscOMPParallelFor(..., numthreads) > for ( ... ) > ; > > would turn into > if ( numthreads == 1 ) > for ( ... ) > ; > else > #pragma omp parallel for > for ( ... ) > ; > > The branch should hopefully be cheaper than the thread sync overhead. > The branch is trivially cheap, the growth in binary size (having a threaded and non-threaded) is possibly an issue. > As far as attaching the threading information to the object goes, it would > be nice to do so staying within the sequential implementation (rather than > copying all the code as is done for pthreads right now). I'd envisage > something like this: > > PetscPragmaOMPObject(vec, parallel) > { > PetscPragmaOMP(for...) > for ( i = PetscVecLowerBound(vec); i < PetscVecUpperBound(vec); i++ ) > ...; > } > > Where one has: > > #if defined(PETSC_HAVE_OPENMP) > #define PetscVecLowerBound(vec) compute_lower_bound > #define PetscVecUpperBound(vec) compute_upper_bound > #else > #define PetscVecLowerBound(vec) 0 > #define PetscVecUpperBound(vec) vec->map->n > #endif > > Or is this too ugly for words? > The reason for my magic was simply to avoid calling omp_get_num_threads() and omp_get_thread_num() many times to make this determination. Note that OpenMP implementations create functions (run by the threads) that capture the local state and call these functions to find their work unit. Thus the generated code is quite similar to PETSc's pthreads versions. > So to take a single, representative, kernel we move from something like: > > PetscErrorCode VecConjugate_Seq(Vec xin) > { > PetscScalar *x; > PetscInt n = xin->map->n; > PetscInt i; > PetscErrorCode ierr; > > PetscFunctionBegin; > ierr = VecGetArray(xin,&x);CHKERRQ(**ierr); > #pragma omp parallel for schedule(static) default(none) private(i) > shared(n, x) > for (i=0; i<n; i++) x[i] = PetscConj(x[i]); > ierr = VecRestoreArray(xin,&x);**CHKERRQ(ierr); > PetscFunctionReturn(0); > } > > To something like this: > > PetscErrorCode VecConjugate_Seq(Vec xin) > { > PetscScalar *x; > PetscInt i; > PetscErrorCode ierr; > PetscFunctionBegin; > ierr = VecGetArray(xin,&x);CHKERRQ(**ierr); > PetscPragmaOMPObject(xin, parallel private(i) shared(xin, x)) > { > for (i=PetscVecLowerBound(xin); i<PetscVecUpperBound(xin); i++) > x[i] = PetscConj(x[i]); > } > ierr = VecRestoreArray(xin,&x);**CHKERRQ(ierr); > PetscFunctionReturn(0); > } > > I prefer inlining the loop bounds into macro/function calls at point of > use, rather than relying on magic names from a macro definition, but I'm > not completely wedded to the idea. > Not a big deal, we could optimize either variant.

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Dmitry Karpeev 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 2:06 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Tue, Feb 28, 2012 at 13:52, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > >> I don't understand the need for multiple function evaluation routines and >> their explicit naming. >> I think the situation here is completely analogous to the situation with >> matrices: I thought we agreed that each DM would be capable of generating a >> single matrix, the "Jacobian", rather than multiple matrices like J and >> J_pre, and that the different necessary matrices (again, e.g., J and J_pre) >> would instead be produced by their respective DMs -- the KSP DM and the PC >> DM, so to speak. These subDMs would be pulled out of the SNES DM and >> passed to the corresponding KSP and PC objects. Thus, a DM encapsulates a >> single "problem", and if we want to get a different problem -- the "exact" >> linearization of the original, for example -- we extract the KSP DM; >> likewise, PC DM is an approximate linearization, and so on. >> >> I don't think TS is any different -- we can have the TSI DM (maybe under >> a better name), which could have DMFormFunction encapsulating the IFunction >> (and, incidentally, a linearization encapsulating the mass matrix), etc. >> > > I'm concerned that this becomes insanity. A single TS might need > RHSFunction, IFunction, and IJacobian. They might also partition terms > differently for the functions and the Jacobian. Having different DMs seems > awfully confusing. > > In short -- I'm for each DM having one function (DMFormFunction or, >> better, DMFormResidual) and one matrix (DMFormJacobian) -- *a* >> linearization of the residual. All other functions and matrices are >> obtained from the corresponding subDMs. >> > > The interface is actually different for these (compare TSIFunction, a > transient residual, to SNESFunction) > Okay. I would still refrain from calling these things by the class of the object consuming them. DMFormFunction() and DMFormIFunction() seem fine, since these have a meaning independent of TS (although the DM has to then encode a transient problem). But I don't like DMGetSNESFunction() any more than I like DMGetKSPJacobian() and DMGetPCJacobian(). By the way, should we start using Residual instead of Function? Dmitry.

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Jed Brown 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 13:52, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > I don't understand the need for multiple function evaluation routines and > their explicit naming. > I think the situation here is completely analogous to the situation with > matrices: I thought we agreed that each DM would be capable of generating a > single matrix, the "Jacobian", rather than multiple matrices like J and > J_pre, and that the different necessary matrices (again, e.g., J and J_pre) > would instead be produced by their respective DMs -- the KSP DM and the PC > DM, so to speak. These subDMs would be pulled out of the SNES DM and > passed to the corresponding KSP and PC objects. Thus, a DM encapsulates a > single "problem", and if we want to get a different problem -- the "exact" > linearization of the original, for example -- we extract the KSP DM; > likewise, PC DM is an approximate linearization, and so on. > > I don't think TS is any different -- we can have the TSI DM (maybe under a > better name), which could have DMFormFunction encapsulating the IFunction > (and, incidentally, a linearization encapsulating the mass matrix), etc. > I'm concerned that this becomes insanity. A single TS might need RHSFunction, IFunction, and IJacobian. They might also partition terms differently for the functions and the Jacobian. Having different DMs seems awfully confusing. In short -- I'm for each DM having one function (DMFormFunction or, better, > DMFormResidual) and one matrix (DMFormJacobian) -- *a* linearization of the > residual. All other functions and matrices are obtained from the > corresponding subDMs. > The interface is actually different for these (compare TSIFunction, a transient residual, to SNESFunction)

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Dmitry Karpeev 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 1:24 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Tue, Feb 28, 2012 at 12:13, Matthew Knepley <knepley(a)gmail.com> wrote: > >> But this is what is causing the trouble. Vecs need to hold a DM to tell >> them how to write themselves, but DMs need >> those Vecs to pass to sovlers. These are very different jobs, and thus >> cause loops. >> > > Most possible loops are managed by deterministic ordering of destruction. > Work vectors are counted. If the only problem we need to solve is reference > count cycle analysis for PetscObjectCompose, we can do that. > > > It sounds like at present, nobody has complained loudly about DM becoming > the one-stop place for defining the problem to be solved (at least the > nonlinear problem). Is that true? So do we move DMSetFunction() to > DMSetSNESFunction() and add DMSetTSIFunction(), etc? > I don't understand the need for multiple function evaluation routines and their explicit naming. I think the situation here is completely analogous to the situation with matrices: I thought we agreed that each DM would be capable of generating a single matrix, the "Jacobian", rather than multiple matrices like J and J_pre, and that the different necessary matrices (again, e.g., J and J_pre) would instead be produced by their respective DMs -- the KSP DM and the PC DM, so to speak. These subDMs would be pulled out of the SNES DM and passed to the corresponding KSP and PC objects. Thus, a DM encapsulates a single "problem", and if we want to get a different problem -- the "exact" linearization of the original, for example -- we extract the KSP DM; likewise, PC DM is an approximate linearization, and so on. I don't think TS is any different -- we can have the TSI DM (maybe under a better name), which could have DMFormFunction encapsulating the IFunction (and, incidentally, a linearization encapsulating the mass matrix), etc. The situation with splitting the problem into subdomains or subfields -- with DMs encapsulating these subproblems -- is completely analogous. And if we really wanted to eliminate the dependence on SNES, KSP, TS, etc. -- even the logical dependence in the names of the functions calls like DMGetKSPDM -- we could resort to object composition. In short -- I'm for each DM having one function (DMFormFunction or, better, DMFormResidual) and one matrix (DMFormJacobian) -- *a* linearization of the residual. All other functions and matrices are obtained from the corresponding subDMs. Dmitry. > To avoid the dependency loop, those methods would actually live in > src/snes and src/ts, but their argument would be a DM. The actual > evaluation context for SNES would be composed with the DM, so it wouldn't > actually reside in _p_DM. > > Do we preserve SNESSetFunction() which would do SNESGetDM() and > DMSetSNESFunction(). If no DM had been set, SNESGetDM() would create a > dummy DM just as a place to attach its containers. > > The DMComputeSNESFunction() would always be called by > SNESComputeFunction(). If the DM was just the dummy, it would call straight > through to the user's implementation, even forwarding the SNES along as the > context. If the user had called DMSetSNESFunction() (or > DMDASetLocalSNESFunction(), etc), that version would be called. (This is > similar to current logic, but dispatch _always_ goes through some DM, even > if it's just a dummy to hold function pointers.) >

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Jed Brown 28 Feb '12

28 Feb '12

On Tue, Feb 28, 2012 at 12:13, Matthew Knepley <knepley(a)gmail.com> wrote: > But this is what is causing the trouble. Vecs need to hold a DM to tell > them how to write themselves, but DMs need > those Vecs to pass to sovlers. These are very different jobs, and thus > cause loops. > Most possible loops are managed by deterministic ordering of destruction. Work vectors are counted. If the only problem we need to solve is reference count cycle analysis for PetscObjectCompose, we can do that. It sounds like at present, nobody has complained loudly about DM becoming the one-stop place for defining the problem to be solved (at least the nonlinear problem). Is that true? So do we move DMSetFunction() to DMSetSNESFunction() and add DMSetTSIFunction(), etc? To avoid the dependency loop, those methods would actually live in src/snes and src/ts, but their argument would be a DM. The actual evaluation context for SNES would be composed with the DM, so it wouldn't actually reside in _p_DM. Do we preserve SNESSetFunction() which would do SNESGetDM() and DMSetSNESFunction(). If no DM had been set, SNESGetDM() would create a dummy DM just as a place to attach its containers. The DMComputeSNESFunction() would always be called by SNESComputeFunction(). If the DM was just the dummy, it would call straight through to the user's implementation, even forwarding the SNES along as the context. If the user had called DMSetSNESFunction() (or DMDASetLocalSNESFunction(), etc), that version would be called. (This is similar to current logic, but dispatch _always_ goes through some DM, even if it's just a dummy to hold function pointers.)

1 0