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

Re: [petsc-dev] What needs to be done before a release?
by Barry Smith 29 May '12

29 May '12

On May 29, 2012, at 3:29 AM, Blaise Bourdin wrote: > > On May 29, 2012, at 1:18 AM, Barry Smith wrote: > >> >> Blaise, >> >> How about the attached? My goal is to have the usage be as similar to the other PetscOptions... calls and too not require the Fortran user to be C aware that is, they should not have to append null characters or declare special C type declaration. In other words, the user writes Fortran77 code :-).[see attached file: ex1f.F90] >> >> I don't think it should bleed, but then I am no longer a Fortran programmer. > > Once a fortran programmer, always a fortran programmer... It's not like there is a lot to forget. > > Now I understand what you want. I agree with your goal of making the call standard and transparent. There is one issue with this approach for a function like PetscBagRegisterEnum (as far as I understand the C code), where you don't copy the string array in the bag item structure, but keep a reference to it (l. 80 in bag.c). After Carray and list1 are deallocated, the bag item->list points to some nonsense. I lost my mind trying to find a workaround last week. > > Would you be OK with requesting that each of the strings in list have some padding at the end, and be modified in the wrapper (adding the C_NULL_CHAR)? The alternative would be to modify PetscBagRegisterEnum, in order to store the content of list instead of just a reference, which may not be something you want to do right now. Actually. I think PetscBagRegisterEnum() should be modified to copy the list, it is an unsafe function even from C now. Implicitly I assumed that list was declared as a global static variable (hence the current PetscBagRegisterEnum() is ok) but I should not have made that assumption. Barry > > Blaise > > >> >> Barry >> >> >> >> On May 27, 2012, at 10:16 AM, Blaise Bourdin wrote: >> >>> Barry, >>> >>> If I frame your email, will you sign it for me? ;) >>> >>> I am traveling so I did not have much time to do it, even in my sleep, but Tim's solution is essentially right. >>> >>> Have a look at the attached code. I had to write my own interface to PetscOptionsGetEnum, but I have to do the equivalent of FIXCHAR by myself (basically appending //C_NULL_CHAR). I checked in the debugger, and the list argument that is passed to PetscOptionsGetEList is the same as the one passed from C. Is this what you want? >>> >>> Do you want to do the FIXCHAR business in the C or fortran side? I suspect that you would prefer not having to append the C_NULL_CHAR in fortran, but I can't figure how to do that without allocating a new array of strings, and causing a memory leak. What do you think? >>> >>> If this is not what you want, could you meet me half way and write the fortran interface you would like to use, in which case I will try to write the missing fortran parts. >>> >>> Blaise >>> <PetscF2003.F90> >>>> >>>> Tim, >>>> >>>> Thanks. But what I want is all bundled up a routine that takes in the Fortran array of character strings and allocates and copies over them to the C ptr business. Likely I could figure it out myself with but Blaise is such a Fortran wizard he can do it optimally in his sleep :-) >>>> >>>> Barry >>>> >>>> On May 25, 2012, at 8:14 PM, Tim Gallagher wrote: >>>> >>>>> For what it's worth, I answered a question about how to pass string arrays awhile ago on StackOverflow: >>>>> >>>>> http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-b… >>>>> >>>>> Tim >>>>> >>>>> ----- Original Message ----- >>>>> From: "Barry Smith" <bsmith(a)mcs.anl.gov> >>>>> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> >>>>> Sent: Friday, May 25, 2012 8:57:52 PM >>>>> Subject: Re: [petsc-dev] What needs to be done before a release? >>>>> >>>>> Blaise >>>>> >>>>> If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. >>>>> >>>>> Barry >>>>> >>>>> On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: >>>>> >>>>>> Barry, >>>>>> >>>>>>>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>>>>>>> is a bit tricky since it involves passing fortran arrays of strings. >>>>>>> >>>>>>> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). >>>>>> >>>>>> I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. >>>>>> >>>>>> Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. >>>>>> >>>>>> Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? >>>>>> >>>>>> Blaise >>>>>> >>>>>> PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 >>>>>> >>>>>> <PetscF2003Binding.F90> >>>>>> -- >>>>>> Department of Mathematics and Center for Computation & Technology >>>>>> Louisiana State University, Baton Rouge, LA 70803, USA >>>>>> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>> >>>> >>> >>> -- >>> Department of Mathematics and Center for Computation & Technology >>> Louisiana State University, Baton Rouge, LA 70803, USA >>> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >>> >>> >>> >>> >>> >>> >>> >> >> <ex1f.F90> > > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin > > > > > > >

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Dmitry Karpeev 29 May '12

29 May '12

On Mon, May 28, 2012 at 5:22 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Mon, May 28, 2012 at 5:10 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: > >> On Mon, May 28, 2012 at 5:02 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: >> >>> This is an enormous amount of new code and changes to put in during a >>> code freeze, but my test case works now. >> >> We can revert, if it starts causing problems. The extra code is due to >> the API chance, hence several new interface functions. >> >> >>> A couple comments from looking at the patch: >>> >>> I'm not sure I like the direction of the link in DMSetEmbedding(DM dm, >>> IS embedding, DM ambientdm). Normally the composite DM knows how to fit its >>> pieces together, but the sub-DM does not have "upward links" to the global >>> problem. >>> >> This is to be able to combine domain and field splits: sometimes a subDM >> has to know how to "relativize" itself with respect to another subDM. >> > > I don't believe it has to. The current approach in PETSc has always been > that we store links going the other direction. I don't like the > inconsistency and I think the model is much more sustainable to hold > downward links only. > I removed the embedding code, so now there is just about the minimum necessary to resolve the issue of field/subdomain decomposition mixing: http://petsc.cs.iit.edu/petsc/petsc-dev/rev/b207eabfcac1 Little of PETSc's code depends on this, but external DM impls (e.g., DM_libMesh) need this API to serve up decompositions correctly. > > What algorithm can only be implemented with these upward links? > The "upward" links are still necessary at least in my experience, but they can be hidden in the impl of the decomposition DM. The difference with other PETSc code (I assume you are alluding to DMComposite, which only needs "downward" links to its constituent pieces) is that of putting things together versus taking them apart: DMComposite inductively _defines_ a global numbering by ordering the constituent pieces with an existing relative numbering. Those pieces don't need to know whether/how their parent composite DM might be embedded into something else. A DM wrapping a piece of an existing mesh, on the other hand, already has the global numbering of that piece, but may need to compute the relative numbering inside another, bigger mesh piece. To do that, the smaller DM needs to know how the ambient DM is embedded into the overall mesh. For example, say I have a problem with variables (u,v,p) numbered globally by libMesh. I first want to split off p decomposing the system into (u,v) and (p), and then further decompose the first piece into (u) and (v). I only have the global numberings of the u, v and p fields to work with. In order to know how u is embedded into (u,v) I need to know (u,v) --> (u,v,p) and then make (u) --> (u,v,p) factor through that. The problem is that there can in principle be an arbitrary number of various decompositions, so storing all of the "downward" links is impractical. In contrast, DMComposite only deals with a single decomposition. Of course, the decomposition DM obtained via DMCreateField/DomainDecomposition() is supposed to resolve that issue by serving as a proxy for a particular collection of "downward" or "upward" links in an impl-defined way. Incidentally, another use of embedding for relativization is in composing active set DMs with (other) field/domain decompositions: if I am only solving a subsystem defined by (in)active variables, want to use PCFIELDSPLIT or PC(G)ASM with it, and want to use the splits/subdomains defined by the original DM, I would have to recalculate the ISs defining the original subdomains, as I'm now only interested in their intersections with the smaller (in)active space. > > >> >>> There are dead links in your documentation because you reference >>> functions that don't exist (were renamed), e.g. DMCreateDecomposition() in >>> the DMSetEmbedding() man page. There may be more. You should generate the >>> docs and check the man pages. >>> >> These should be gone now. > >>> How does your interface support "Neumann" subdomains? >>> >> >> I think this issue is orthogonal to the definition of the decomposition, >> since there has to be additional support at the PC >> and Mat level to assemble the Neumann subproblems. At this point you >> would have to do some combination of PCASMSetModifySubmatrices() with >> nullspace removal. >> > > Sure, you have to assemble to a different format, but what interface are > we going to use. In the linear case, we can, in principle, do everything > with matrices. For the nonlinear case, we need DM APIs. So in the nonlinear > context, how are we going to build Neumann problems? > I think the "PETSc-facing" API is more or less clear: DMFunction, DMJacobian on each of the decomposition pieces define the subproblem. The issue with the user-facing API is how to make defining the subproblems relatively easy? DMDAComputeFunctionLocal()? In general you'd want something like the Moose API: the user provides a quadrature-point kernel (maybe that's too fine a granularity; an element-loop instead?), and the DM traverses the subdomain and puts the kernel invocations together? Dmitry. > > >> Your DMCreateDomainDecomposition() seems to provide nominally overlapping >>> index sets, but how do we determine what part is owned? (RASM, especially, >>> needs both an overlapping decomposition and a (non-overlapping) partition. >>> >> My intent was to have these be nonoverlapping and expand them with >> another call (e.g., DMDomainDecompositionIncreaseOverlap(). >> > > But with separate functions, we have to jump through extra hoops to use a > user-specified decomposition (e.g. that carefully selects Lagrange > multipliers or pressure variables to maintain inf-sup stability). > > >> The current API plays better with GASM, since that PC assumes that the >> basic subdomains are nonoverlapping. That fix will have to wait for the >> next release, I guess. >> > > PCGASMSetSubdomains() also has two IS arguments, so I don't know what you > mean. > > > We need to work out an interface for nonlinear DD (like ASPIN) this > summer, but I don't think it belongs in this release. I'm not so wild about > putting new interfaces into this release that will be changed immediately > after. >

1 0

Re: [petsc-dev] What needs to be done before a release?
by Blaise Bourdin 29 May '12

29 May '12

On May 29, 2012, at 1:18 AM, Barry Smith wrote: > > Blaise, > > How about the attached? My goal is to have the usage be as similar to the other PetscOptions... calls and too not require the Fortran user to be C aware that is, they should not have to append null characters or declare special C type declaration. In other words, the user writes Fortran77 code :-).[see attached file: ex1f.F90] > > I don't think it should bleed, but then I am no longer a Fortran programmer. Once a fortran programmer, always a fortran programmer... It's not like there is a lot to forget. Now I understand what you want. I agree with your goal of making the call standard and transparent. There is one issue with this approach for a function like PetscBagRegisterEnum (as far as I understand the C code), where you don't copy the string array in the bag item structure, but keep a reference to it (l. 80 in bag.c). After Carray and list1 are deallocated, the bag item->list points to some nonsense. I lost my mind trying to find a workaround last week. Would you be OK with requesting that each of the strings in list have some padding at the end, and be modified in the wrapper (adding the C_NULL_CHAR)? The alternative would be to modify PetscBagRegisterEnum, in order to store the content of list instead of just a reference, which may not be something you want to do right now. Blaise > > Barry > > > > On May 27, 2012, at 10:16 AM, Blaise Bourdin wrote: > >> Barry, >> >> If I frame your email, will you sign it for me? ;) >> >> I am traveling so I did not have much time to do it, even in my sleep, but Tim's solution is essentially right. >> >> Have a look at the attached code. I had to write my own interface to PetscOptionsGetEnum, but I have to do the equivalent of FIXCHAR by myself (basically appending //C_NULL_CHAR). I checked in the debugger, and the list argument that is passed to PetscOptionsGetEList is the same as the one passed from C. Is this what you want? >> >> Do you want to do the FIXCHAR business in the C or fortran side? I suspect that you would prefer not having to append the C_NULL_CHAR in fortran, but I can't figure how to do that without allocating a new array of strings, and causing a memory leak. What do you think? >> >> If this is not what you want, could you meet me half way and write the fortran interface you would like to use, in which case I will try to write the missing fortran parts. >> >> Blaise >> <PetscF2003.F90> >>> >>> Tim, >>> >>> Thanks. But what I want is all bundled up a routine that takes in the Fortran array of character strings and allocates and copies over them to the C ptr business. Likely I could figure it out myself with but Blaise is such a Fortran wizard he can do it optimally in his sleep :-) >>> >>> Barry >>> >>> On May 25, 2012, at 8:14 PM, Tim Gallagher wrote: >>> >>>> For what it's worth, I answered a question about how to pass string arrays awhile ago on StackOverflow: >>>> >>>> http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-b… >>>> >>>> Tim >>>> >>>> ----- Original Message ----- >>>> From: "Barry Smith" <bsmith(a)mcs.anl.gov> >>>> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> >>>> Sent: Friday, May 25, 2012 8:57:52 PM >>>> Subject: Re: [petsc-dev] What needs to be done before a release? >>>> >>>> Blaise >>>> >>>> If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. >>>> >>>> Barry >>>> >>>> On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: >>>> >>>>> Barry, >>>>> >>>>>>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>>>>>> is a bit tricky since it involves passing fortran arrays of strings. >>>>>> >>>>>> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). >>>>> >>>>> I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. >>>>> >>>>> Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. >>>>> >>>>> Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? >>>>> >>>>> Blaise >>>>> >>>>> PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 >>>>> >>>>> <PetscF2003Binding.F90> >>>>> -- >>>>> Department of Mathematics and Center for Computation & Technology >>>>> Louisiana State University, Baton Rouge, LA 70803, USA >>>>> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>> >>> >> >> -- >> Department of Mathematics and Center for Computation & Technology >> Louisiana State University, Baton Rouge, LA 70803, USA >> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >> >> >> >> >> >> >> > > <ex1f.F90> -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Re: [petsc-dev] What needs to be done before a release?
by Barry Smith 28 May '12

28 May '12

Blaise, How about the attached? My goal is to have the usage be as similar to the other PetscOptions... calls and too not require the Fortran user to be C aware that is, they should not have to append null characters or declare special C type declaration. In other words, the user writes Fortran77 code :-).[see attached file: ex1f.F90] I don't think it should bleed, but then I am no longer a Fortran programmer. Barry On May 27, 2012, at 10:16 AM, Blaise Bourdin wrote: > Barry, > > If I frame your email, will you sign it for me? ;) > > I am traveling so I did not have much time to do it, even in my sleep, but Tim's solution is essentially right. > > Have a look at the attached code. I had to write my own interface to PetscOptionsGetEnum, but I have to do the equivalent of FIXCHAR by myself (basically appending //C_NULL_CHAR). I checked in the debugger, and the list argument that is passed to PetscOptionsGetEList is the same as the one passed from C. Is this what you want? > > Do you want to do the FIXCHAR business in the C or fortran side? I suspect that you would prefer not having to append the C_NULL_CHAR in fortran, but I can't figure how to do that without allocating a new array of strings, and causing a memory leak. What do you think? > > If this is not what you want, could you meet me half way and write the fortran interface you would like to use, in which case I will try to write the missing fortran parts. > > Blaise > <PetscF2003.F90> >> >> Tim, >> >> Thanks. But what I want is all bundled up a routine that takes in the Fortran array of character strings and allocates and copies over them to the C ptr business. Likely I could figure it out myself with but Blaise is such a Fortran wizard he can do it optimally in his sleep :-) >> >> Barry >> >> On May 25, 2012, at 8:14 PM, Tim Gallagher wrote: >> >>> For what it's worth, I answered a question about how to pass string arrays awhile ago on StackOverflow: >>> >>> http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-b… >>> >>> Tim >>> >>> ----- Original Message ----- >>> From: "Barry Smith" <bsmith(a)mcs.anl.gov> >>> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> >>> Sent: Friday, May 25, 2012 8:57:52 PM >>> Subject: Re: [petsc-dev] What needs to be done before a release? >>> >>> Blaise >>> >>> If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. >>> >>> Barry >>> >>> On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: >>> >>>> Barry, >>>> >>>>>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>>>>> is a bit tricky since it involves passing fortran arrays of strings. >>>>> >>>>> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). >>>> >>>> I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. >>>> >>>> Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. >>>> >>>> Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? >>>> >>>> Blaise >>>> >>>> PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 >>>> >>>> <PetscF2003Binding.F90> >>>> -- >>>> Department of Mathematics and Center for Computation & Technology >>>> Louisiana State University, Baton Rouge, LA 70803, USA >>>> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>> >> > > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin > > > > > > >

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Jed Brown 28 May '12

28 May '12

On Mon, May 28, 2012 at 5:10 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > On Mon, May 28, 2012 at 5:02 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> This is an enormous amount of new code and changes to put in during a >> code freeze, but my test case works now. > > We can revert, if it starts causing problems. The extra code is due to > the API chance, hence several new interface functions. > > >> A couple comments from looking at the patch: >> >> I'm not sure I like the direction of the link in DMSetEmbedding(DM dm, IS >> embedding, DM ambientdm). Normally the composite DM knows how to fit its >> pieces together, but the sub-DM does not have "upward links" to the global >> problem. >> > This is to be able to combine domain and field splits: sometimes a subDM > has to know how to "relativize" itself with respect to another subDM. > I don't believe it has to. The current approach in PETSc has always been that we store links going the other direction. I don't like the inconsistency and I think the model is much more sustainable to hold downward links only. What algorithm can only be implemented with these upward links? > >> There are dead links in your documentation because you reference >> functions that don't exist (were renamed), e.g. DMCreateDecomposition() in >> the DMSetEmbedding() man page. There may be more. You should generate the >> docs and check the man pages. >> > >> How does your interface support "Neumann" subdomains? >> > > I think this issue is orthogonal to the definition of the decomposition, > since there has to be additional support at the PC > and Mat level to assemble the Neumann subproblems. At this point you > would have to do some combination of PCASMSetModifySubmatrices() with > nullspace removal. > Sure, you have to assemble to a different format, but what interface are we going to use. In the linear case, we can, in principle, do everything with matrices. For the nonlinear case, we need DM APIs. So in the nonlinear context, how are we going to build Neumann problems? > Your DMCreateDomainDecomposition() seems to provide nominally overlapping >> index sets, but how do we determine what part is owned? (RASM, especially, >> needs both an overlapping decomposition and a (non-overlapping) partition. >> > My intent was to have these be nonoverlapping and expand them with another > call (e.g., DMDomainDecompositionIncreaseOverlap(). > But with separate functions, we have to jump through extra hoops to use a user-specified decomposition (e.g. that carefully selects Lagrange multipliers or pressure variables to maintain inf-sup stability). > The current API plays better with GASM, since that PC assumes that the > basic subdomains are nonoverlapping. That fix will have to wait for the > next release, I guess. > PCGASMSetSubdomains() also has two IS arguments, so I don't know what you mean. We need to work out an interface for nonlinear DD (like ASPIN) this summer, but I don't think it belongs in this release. I'm not so wild about putting new interfaces into this release that will be changed immediately after.

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Dmitry Karpeev 28 May '12

28 May '12

On Mon, May 28, 2012 at 5:02 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > This is an enormous amount of new code and changes to put in during a code > freeze, but my test case works now. We can revert, if it starts causing problems. The extra code is due to the API chance, hence several new interface functions. > A couple comments from looking at the patch: > > I'm not sure I like the direction of the link in DMSetEmbedding(DM dm, IS > embedding, DM ambientdm). Normally the composite DM knows how to fit its > pieces together, but the sub-DM does not have "upward links" to the global > problem. > This is to be able to combine domain and field splits: sometimes a subDM has to know how to "relativize" itself with respect to another subDM. > > There are dead links in your documentation because you reference functions > that don't exist (were renamed), e.g. DMCreateDecomposition() in the > DMSetEmbedding() man page. There may be more. You should generate the docs > and check the man pages. > > How does your interface support "Neumann" subdomains? > I think this issue is orthogonal to the definition of the decomposition, since there has to be additional support at the PC and Mat level to assemble the Neumann subproblems. At this point you would have to do some combination of PCASMSetModifySubmatrices() with nullspace removal. > > Your DMCreateDomainDecomposition() seems to provide nominally overlapping > index sets, but how do we determine what part is owned? (RASM, especially, > needs both an overlapping decomposition and a (non-overlapping) partition. > My intent was to have these be nonoverlapping and expand them with another call (e.g., DMDomainDecompositionIncreaseOverlap(). The current API plays better with GASM, since that PC assumes that the basic subdomains are nonoverlapping. That fix will have to wait for the next release, I guess. Dmitry. > > On Mon, May 28, 2012 at 4:32 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: > >> I refactored the code to get rid of field/domain decomposition mixing. >> The use of DM in ASM is reenabled here: >> http://petsc.cs.iit.edu/petsc/petsc-dev/rev/2e2ee6a01d1a >> I didn't see any errors related to it in alltest, but I'll keep an eye on >> the nightlies. >> >> Dmitry. >> >> >> On Sun, May 27, 2012 at 9:40 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: >> >>> In particular, read through my first email in this thread. I explain how >>> your implementation leads to several instances of either the wrong >>> algorithm or a crash. >>> >>> I pushed SNES runex48_5, it is a basic test of PCASM using AIJ, BAIJ, >>> and SBAIJ. >>> >>> If you revert http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd, >>> the AIJ test will take much more iterations and the other two will crash. >>> >>> >>> On Sun, May 27, 2012 at 9:02 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: >>> >>>> ./ex48 -pc_type asm >>>> >>>> does the wrong algorithm. >>>> On May 27, 2012 8:56 PM, "Dmitry Karpeev" <karpeev(a)mcs.anl.gov> wrote: >>>> >>>>> I think I have a fix and want to reenable DM to serve up >>>>> decompositions to PCASM. >>>>> Was there a particular test that was failing? I don't see any >>>>> decomposition-related errors in alltest. >>>>> Dmitry. >>>>> >>>>> On Sat, May 26, 2012 at 10:38 PM, Jed Brown <jedbrown(a)mcs.anl.gov>wrote: >>>>> >>>>>> On Sat, May 26, 2012 at 10:24 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov >>>>>> > wrote: >>>>>> >>>>>>> How about having a DM-specific function for that designation? >>>>>>>> >>>>>>> DMXXXGetFieldSplits()/DMXXXGetDomainDecompositions(), for example? >>>>>>> >>>>>> >>>>>> For the C interface, maybe configure a named set with DMXxxSet... and >>>>>> tell PCASM to use it. For the options database interface, maybe have the >>>>>> option processed by DMXxx code instead of by PCASM. I don't know, this >>>>>> might be just as bad. >>>>>> >>>>>> >>>>>>> I think that would work, if the implementation defined only a small >>>>>>> number of decompositions. >>>>>>> If instead there is a short grammar defining well-formed strings >>>>>>> naming such decompositions, I'm not sure how to >>>>>>> communicate that to the user, other than in documentation. I have >>>>>>> DMLibMeshGetVars() to return all vars. FieldSplits >>>>>>> can be defined as arbitrary groupings of those. Each split will >>>>>>> then advertise only the vars in the groupings that defines the split. >>>>>>> >>>>>> >>>>>> I still don't know what to do about a nonlinear change of basis. >>>>>> Advertising var names seems reasonable to me. >>>>>> >>>>>> Note that ASM domains are essentially vertex-oriented which is >>>>>> different from element-centered domains (natural for non-overlapping >>>>>> methods). >>>>>> >>>>> >>>>> >>> >> >

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Jed Brown 28 May '12

28 May '12

This is an enormous amount of new code and changes to put in during a code freeze, but my test case works now. A couple comments from looking at the patch: I'm not sure I like the direction of the link in DMSetEmbedding(DM dm, IS embedding, DM ambientdm). Normally the composite DM knows how to fit its pieces together, but the sub-DM does not have "upward links" to the global problem. There are dead links in your documentation because you reference functions that don't exist (were renamed), e.g. DMCreateDecomposition() in the DMSetEmbedding() man page. There may be more. You should generate the docs and check the man pages. How does your interface support "Neumann" subdomains? Your DMCreateDomainDecomposition() seems to provide nominally overlapping index sets, but how do we determine what part is owned? (RASM, especially, needs both an overlapping decomposition and a (non-overlapping) partition. On Mon, May 28, 2012 at 4:32 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > I refactored the code to get rid of field/domain decomposition mixing. > The use of DM in ASM is reenabled here: > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/2e2ee6a01d1a > I didn't see any errors related to it in alltest, but I'll keep an eye on > the nightlies. > > Dmitry. > > > On Sun, May 27, 2012 at 9:40 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> In particular, read through my first email in this thread. I explain how >> your implementation leads to several instances of either the wrong >> algorithm or a crash. >> >> I pushed SNES runex48_5, it is a basic test of PCASM using AIJ, BAIJ, and >> SBAIJ. >> >> If you revert http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd, >> the AIJ test will take much more iterations and the other two will crash. >> >> >> On Sun, May 27, 2012 at 9:02 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: >> >>> ./ex48 -pc_type asm >>> >>> does the wrong algorithm. >>> On May 27, 2012 8:56 PM, "Dmitry Karpeev" <karpeev(a)mcs.anl.gov> wrote: >>> >>>> I think I have a fix and want to reenable DM to serve up decompositions >>>> to PCASM. >>>> Was there a particular test that was failing? I don't see any >>>> decomposition-related errors in alltest. >>>> Dmitry. >>>> >>>> On Sat, May 26, 2012 at 10:38 PM, Jed Brown <jedbrown(a)mcs.anl.gov>wrote: >>>> >>>>> On Sat, May 26, 2012 at 10:24 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: >>>>> >>>>>> How about having a DM-specific function for that designation? >>>>>>> >>>>>> DMXXXGetFieldSplits()/DMXXXGetDomainDecompositions(), for example? >>>>>> >>>>> >>>>> For the C interface, maybe configure a named set with DMXxxSet... and >>>>> tell PCASM to use it. For the options database interface, maybe have the >>>>> option processed by DMXxx code instead of by PCASM. I don't know, this >>>>> might be just as bad. >>>>> >>>>> >>>>>> I think that would work, if the implementation defined only a small >>>>>> number of decompositions. >>>>>> If instead there is a short grammar defining well-formed strings >>>>>> naming such decompositions, I'm not sure how to >>>>>> communicate that to the user, other than in documentation. I have >>>>>> DMLibMeshGetVars() to return all vars. FieldSplits >>>>>> can be defined as arbitrary groupings of those. Each split will then >>>>>> advertise only the vars in the groupings that defines the split. >>>>>> >>>>> >>>>> I still don't know what to do about a nonlinear change of basis. >>>>> Advertising var names seems reasonable to me. >>>>> >>>>> Note that ASM domains are essentially vertex-oriented which is >>>>> different from element-centered domains (natural for non-overlapping >>>>> methods). >>>>> >>>> >>>> >> >

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Dmitry Karpeev 28 May '12

28 May '12

I refactored the code to get rid of field/domain decomposition mixing. The use of DM in ASM is reenabled here: http://petsc.cs.iit.edu/petsc/petsc-dev/rev/2e2ee6a01d1a I didn't see any errors related to it in alltest, but I'll keep an eye on the nightlies. Dmitry. On Sun, May 27, 2012 at 9:40 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > In particular, read through my first email in this thread. I explain how > your implementation leads to several instances of either the wrong > algorithm or a crash. > > I pushed SNES runex48_5, it is a basic test of PCASM using AIJ, BAIJ, and > SBAIJ. > > If you revert http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd, > the AIJ test will take much more iterations and the other two will crash. > > > On Sun, May 27, 2012 at 9:02 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> ./ex48 -pc_type asm >> >> does the wrong algorithm. >> On May 27, 2012 8:56 PM, "Dmitry Karpeev" <karpeev(a)mcs.anl.gov> wrote: >> >>> I think I have a fix and want to reenable DM to serve up decompositions >>> to PCASM. >>> Was there a particular test that was failing? I don't see any >>> decomposition-related errors in alltest. >>> Dmitry. >>> >>> On Sat, May 26, 2012 at 10:38 PM, Jed Brown <jedbrown(a)mcs.anl.gov>wrote: >>> >>>> On Sat, May 26, 2012 at 10:24 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: >>>> >>>>> How about having a DM-specific function for that designation? >>>>>> >>>>> DMXXXGetFieldSplits()/DMXXXGetDomainDecompositions(), for example? >>>>> >>>> >>>> For the C interface, maybe configure a named set with DMXxxSet... and >>>> tell PCASM to use it. For the options database interface, maybe have the >>>> option processed by DMXxx code instead of by PCASM. I don't know, this >>>> might be just as bad. >>>> >>>> >>>>> I think that would work, if the implementation defined only a small >>>>> number of decompositions. >>>>> If instead there is a short grammar defining well-formed strings >>>>> naming such decompositions, I'm not sure how to >>>>> communicate that to the user, other than in documentation. I have >>>>> DMLibMeshGetVars() to return all vars. FieldSplits >>>>> can be defined as arbitrary groupings of those. Each split will then >>>>> advertise only the vars in the groupings that defines the split. >>>>> >>>> >>>> I still don't know what to do about a nonlinear change of basis. >>>> Advertising var names seems reasonable to me. >>>> >>>> Note that ASM domains are essentially vertex-oriented which is >>>> different from element-centered domains (natural for non-overlapping >>>> methods). >>>> >>> >>> >

1 0

Re: [petsc-dev] What needs to be done before a release?
by Blaise Bourdin 27 May '12

27 May '12

Barry, If I frame your email, will you sign it for me? ;) I am traveling so I did not have much time to do it, even in my sleep, but Tim's solution is essentially right. Have a look at the attached code. I had to write my own interface to PetscOptionsGetEnum, but I have to do the equivalent of FIXCHAR by myself (basically appending //C_NULL_CHAR). I checked in the debugger, and the list argument that is passed to PetscOptionsGetEList is the same as the one passed from C. Is this what you want? Do you want to do the FIXCHAR business in the C or fortran side? I suspect that you would prefer not having to append the C_NULL_CHAR in fortran, but I can't figure how to do that without allocating a new array of strings, and causing a memory leak. What do you think? If this is not what you want, could you meet me half way and write the fortran interface you would like to use, in which case I will try to write the missing fortran parts. Blaise > > Tim, > > Thanks. But what I want is all bundled up a routine that takes in the Fortran array of character strings and allocates and copies over them to the C ptr business. Likely I could figure it out myself with but Blaise is such a Fortran wizard he can do it optimally in his sleep :-) > > Barry > > On May 25, 2012, at 8:14 PM, Tim Gallagher wrote: > >> For what it's worth, I answered a question about how to pass string arrays awhile ago on StackOverflow: >> >> http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-b… >> >> Tim >> >> ----- Original Message ----- >> From: "Barry Smith" <bsmith(a)mcs.anl.gov> >> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> >> Sent: Friday, May 25, 2012 8:57:52 PM >> Subject: Re: [petsc-dev] What needs to be done before a release? >> >> Blaise >> >> If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. >> >> Barry >> >> On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: >> >>> Barry, >>> >>>>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>>>> is a bit tricky since it involves passing fortran arrays of strings. >>>> >>>> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). >>> >>> I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. >>> >>> Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. >>> >>> Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? >>> >>> Blaise >>> >>> PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 >>> >>> <PetscF2003Binding.F90> >>> -- >>> Department of Mathematics and Center for Computation & Technology >>> Louisiana State University, Baton Rouge, LA 70803, USA >>> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >>> >>> >>> >>> >>> >>> >>> >> > -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Jed Brown 26 May '12

26 May '12

On Sat, May 26, 2012 at 10:24 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: > How about having a DM-specific function for that designation? >> > DMXXXGetFieldSplits()/DMXXXGetDomainDecompositions(), for example? > For the C interface, maybe configure a named set with DMXxxSet... and tell PCASM to use it. For the options database interface, maybe have the option processed by DMXxx code instead of by PCASM. I don't know, this might be just as bad. > I think that would work, if the implementation defined only a small number > of decompositions. > If instead there is a short grammar defining well-formed strings naming > such decompositions, I'm not sure how to > communicate that to the user, other than in documentation. I have > DMLibMeshGetVars() to return all vars. FieldSplits > can be defined as arbitrary groupings of those. Each split will then > advertise only the vars in the groupings that defines the split. > I still don't know what to do about a nonlinear change of basis. Advertising var names seems reasonable to me. Note that ASM domains are essentially vertex-oriented which is different from element-centered domains (natural for non-overlapping methods).

1 0