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

Re: [petsc-dev] Using cuBlas as the vendor blas for PETSc
by Dave Nystrom 24 Feb '12

24 Feb '12

Hi Jack, Thanks for your comments. I had not thought of the idea of a wrapper. The idea of the overhead with small blocks is certainly worrisome. I have not really played around with blas much in a long time and so don't really have an idea of where the breakeven size might be. I might play around with this at some point just to satisfy my curiosity. Thanks again, Dave Jack Poulson writes: > Dave, > > That will probably not be a very good idea due to the overhead associated > with transferring data to and from the GPU being more expensive than the > computation itself for small problems. This issue can be somewhat avoided > by writing trivial wrappers for routines like dgemm which only run the > multiply on the GPU when the dimensions of the problem are above some > threshold, but this would require slightly more work than simply replacing > BLAS with CUBLAS. > > Jack > > On Fri, Feb 24, 2012 at 8:28 PM, Dave Nystrom <Dave.Nystrom(a)tachyonlogic.com > > wrote: > > > I was wondering if anyone had ever tried using cuBlas as a substitute for > > something like MKL with PETSc. I've been wondering if it would give better > > performance than MKL for my direct solves with cholmod even though the > > block > > sizes are small for cholmod i.e. 32x32 is the default I believe. If so, > > were > > there any tricky aspects to using cuBlas in this way? > > > > Thanks, > > > > Dave > >

1 0

Re: [petsc-dev] Using cuBlas as the vendor blas for PETSc
by Jack Poulson 24 Feb '12

24 Feb '12

Dave, That will probably not be a very good idea due to the overhead associated with transferring data to and from the GPU being more expensive than the computation itself for small problems. This issue can be somewhat avoided by writing trivial wrappers for routines like dgemm which only run the multiply on the GPU when the dimensions of the problem are above some threshold, but this would require slightly more work than simply replacing BLAS with CUBLAS. Jack On Fri, Feb 24, 2012 at 8:28 PM, Dave Nystrom <Dave.Nystrom(a)tachyonlogic.com > wrote: > I was wondering if anyone had ever tried using cuBlas as a substitute for > something like MKL with PETSc. I've been wondering if it would give better > performance than MKL for my direct solves with cholmod even though the > block > sizes are small for cholmod i.e. 32x32 is the default I believe. If so, > were > there any tricky aspects to using cuBlas in this way? > > Thanks, > > Dave >

1 0

Using cuBlas as the vendor blas for PETSc
by Dave Nystrom 24 Feb '12

24 Feb '12

I was wondering if anyone had ever tried using cuBlas as a substitute for something like MKL with PETSc. I've been wondering if it would give better performance than MKL for my direct solves with cholmod even though the block sizes are small for cholmod i.e. 32x32 is the default I believe. If so, were there any tricky aspects to using cuBlas in this way? Thanks, Dave

1 0

Re: [petsc-dev] Understanding Some Parallel Results with PETSc
by Dave Nystrom 24 Feb '12

24 Feb '12

That is what I did for the gpu case where I was running with 2 mpi processes per node since I had 2 gpus per node. But if that is not really the right way to be doing this, I'll be happy to learn of a better way. Note that for the CPU case in the runPetscProb_2 script I use a different binding. Anyway, if there is a better way to do this, I'll be happy to learn and try the better way. Also, this is recent work that has been done in a hurry because of trying to get ready for this new machine coming in. I probably need to go back and revisit the details to make sure I understand the binding I am really getting. Jed Brown writes: > Hmm, I'll look more carefully, but I think your binding is incorrect. > > numactl -l --cpunodebind=$OMPI_COMM_WORLD_LOCAL_RANK ex2 $KSP_ARGS > > NUMA node numbering is different from MPI ranks. > > On Fri, Feb 24, 2012 at 12:11, Nystrom, William D <wdn(a)lanl.gov> wrote: > > > Hi Jed, > > > > Attached is a gzipped tarball of the stuff I used to run these two test > > problems with > > numactl. Actually, I hacked them a bit this morning because I was running > > them in > > our test framework for doing acceptance testing of new systems. But the > > scripts in > > the tarball should give you all the info you need. There is a top level > > script called > > runPetsc that just invokes mpirun from openmpi and calls the wrapper > > scripts for using > > numactl. You could actually dispense with the top level script and just > > invoke the > > mpirun commands yourself. I include it as an easy way to document what I > > did. > > The runPetscProb_1 script runs petsc on the gpus using numactl to control > > the affinities > > of the gpus to the cpu numa nodes. The runPetscProb_2 script runs petsc > > on the cpus > > using numactl. Note that both of those wrapper scripts are using openmpi > > variables. > > I'm not sure how one would do the same thing with another flavor of mpi. > > But I imagine > > it is possible. Also, I'm not sure if there are other more elegant ways > > to run with numactl > > than using the wrapper script approach. Perhaps there is but this is what > > we have been > > doing. > > > > I've also included a Perl script called numa-maps that is useful for > > actually checking the > > affinities that you get while running in order to make sure that numactl > > is doing what > > you think it is doing. I'm not sure where this script comes from. I find > > it on some systems > > and not on others. > > > > I've also include logs with the output of cpuinfo, nvidia-smi and uname -a > > to answer any > > questions you had about the system I was running on. > > > > Finally, I've included runPetscProb_1.log and runPetscProb_2.log which > > contains the > > log_summary output for my latest runs on the gpu and cpu respectively. > > Using numactl > > reduced the runtime for the gpu case as well but not as much as for the > > cpu case. So > > the final result was that running the same problem while using all of the > > gpu resources > > on a node was about 2.5x times faster than using all of the cpu resources > > on the same > > number of nodes. > > > > Let me know if you need more any more info. I'm planning to use this > > stuff to help test > > a new gpu cluster that we have just started acceptance testing on. It has > > the same > > basic hardware as the testbed cluster for these results but has 308 > > nodes. That > > should be interesting and fun. > > > > > > Thanks, > > > > Dave > > > > -- > > Dave Nystrom > > LANL HPC-5 > > Phone: 505-667-7913 > > Email: wdn(a)lanl.gov > > Smail: Mail Stop B272 > > Group HPC-5 > > Los Alamos National Laboratory > > Los Alamos, NM 87545 > > > > ------------------------------ > > *From:* petsc-dev-bounces(a)mcs.anl.gov [petsc-dev-bounces(a)mcs.anl.gov] on > > behalf of Jed Brown [jedbrown(a)mcs.anl.gov] > > *Sent:* Thursday, February 23, 2012 10:43 PM > > > > *To:* For users of the development version of PETSc > > *Cc:* Dave Nystrom > > > > *Subject:* Re: [petsc-dev] Understanding Some Parallel Results with PETSc > > > > On Thu, Feb 23, 2012 at 23:41, Dave Nystrom <dnystrom1(a)comcast.net>wrote: > > > >> I could also send you my mpi/numactl command lines for gpu and cpu when I > >> am > >> back in the office. > >> > > > > Yes, please. > >

1 0

Re: [petsc-dev] --download-mpich --download-openmpi and mpicc
by Barry Smith 24 Feb '12

24 Feb '12

Sean hit the rather obvious nail on the head. Do to the poor design of the mpicc it usually doesn't make sense to have mpicc compile another MPI Satish or Matt could you fix this and put in Sean's error message :-) On Feb 21, 2012, at 11:31 AM, Sean Farley wrote: > I don't think this was the problem. We call updateCompilers() after the install, so it should override anything the user puts in. > > If I understand Barry correctly, it's that people who specify this: > > $ ./configure CC=mpicc --CXX=mpicxx --download-mpich > > clearly do not understand the --download-mpich option and should be told: > > ERROR: User needs to learn MPI; either use the cluster's custom built mpi XOR use --download-mpich > > because if configure && make succeed, then they will irrevocably use the mpiexec or mpicc in their $PATH and then use such great lines as "I have no idea what MPI is!" or "I just want to compile and run on my cluster! Why won't you let me do that?!"

1 0

Re: [petsc-dev] DMComplex and PCFieldSplit
by Barry Smith 24 Feb '12

24 Feb '12

cool On Feb 24, 2012, at 4:23 PM, Matthew Knepley wrote: > On Fri, Feb 24, 2012 at 4:16 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info. > > Found it. So you want PCFieldSplit to troll through those inner loops of > > /* Define fields */ > PetscSection section, sectionGlobal; > PetscInt *fieldSizes, **fieldIndices; > PetscInt nF, f, pStart, pEnd, p; > > ierr = DMComplexGetDefaultSection(pc->dm, &section);CHKERRQ(ierr); > ierr = DMComplexGetDefaultGlobalSection(pc->dm, &sectionGlobal);CHKERRQ(ierr); > ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); > ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt *,&fieldIndices);CHKERRQ(ierr); > ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); > for(f = 0; f < nF; ++f) { > fieldSizes[f] = 0; > } > for(p = pStart; p < pEnd; ++p) { > PetscInt gdof; > > ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); > if (gdof > 0) { > for(f = 0; f < nF; ++f) { > PetscInt fcdof; > > ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); > fieldSizes[f] += fcdof; > } > } > } > for(f = 0; f < nF; ++f) { > ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);CHKERRQ(ierr); > fieldSizes[f] = 0; > } > for(p = pStart; p < pEnd; ++p) { > PetscInt gdof, goff; > > ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); > if (gdof > 0) { > ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); > for(f = 0; f < nF; ++f) { > PetscInt fcdof, fc; > > ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); > for(fc = 0; fc < fcdof; ++fc, ++fieldSizes[f]) { > fieldIndices[f][fieldSizes[f]] = goff++; > } > } > } > } > for(f = 0; f < nF; ++f) { > IS field; > const char *fieldname; > > ierr = PetscSectionGetFieldName(section, f, &fieldname);CHKERRQ(ierr); > ierr = ISCreateGeneral(((PetscObject) pc)->comm, fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &field);CHKERRQ(ierr); > ierr = PetscPrintf(((PetscObject) pc)->comm, "Field %s\n", fieldname);CHKERRQ(ierr); > ierr = ISView(field, PETSC_NULL);CHKERRQ(ierr); > ierr = PCFieldSplitSetIS(pc, fieldname, field);CHKERRQ(ierr); > ierr = ISDestroy(&field);CHKERRQ(ierr); > } > ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr); > } > > every time instead of call DMGetFieldsIS(dm&is); That's insane, just take the crud above and call it DMGetFieldsIS_Complex() and we'll have PCFieldSplit independent of any DM guts. > > I don't understand what we are fighting about. > > Yes, I agree with that. I will move it. > > Matt > > > On Feb 24, 2012, at 3:49 PM, Matthew Knepley wrote: > > > > > It is the correct level of granularity, especially fr FEM which PETSc has almost none of right now > > because we make it unnaturally hard due to our ridiculous infatuation with Finite Difference. That kind > > of thinking colors everything. > > I never said that PetscSection doesn't belong anywhere. It just doesn't seem to belong anywhere in PC/KSP land, it is unnecessarily low for that. Sure it may be ok for the finite element assembly etc. > > Barry > > > > > > Regardless, say the concept of PetscSection as a way to keep the information about fields is fine. PetscSection is a better representation of an array of ISs. > > > > 1) Imbedding these constraint thingys in it seems completely wrong. Information about the constraint things should just be in another PetscSection not inside the PetscSection that tells you the field division. > > > > I disagree since you need this information in every situation where you ask about sizes, etc. Again, there is only ONE > > problem in PETSc that does this, written by me. It is often useful to try something before trying to figure out what is the > > best way to do it. > > > > 2) One could then have DMGetPetscSection(DM,section) that tells you the decomposition, but since the solvers need DMs for the pieces (in PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so different from my first email with PetscSection replacing my array of ISs. So why not fix IS to be a bit more dynamic and build there is no real need for PetscSection except as an array of ISs. > > > > I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info. > > > > Matt > > > > Jed is laughing his pants off at the reaction to his side comment in that email :-) > > > > > > Barry > > > > On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote: > > > > > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote: > > > > > > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote: > > > > > > > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > > > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote: > > > > > > > > > > > We need a richer field interface. Some way for a client of DM to discover vector and tensor structure. > > > > > > > > > > No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS *ises,DM *dms); and DMPushSubDMType(DM,"name of a type of splitting like u,v,p or velocity,p"); DMPopSubDMType(); and DMGetAvailableSubDMTypes(DM,char **); > > > > > > > > > > Do you need more than that? > > > > > > > > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index stuff into the DM interface. I think this is so > > > > > unwieldy and complicated. PetscSection is small, powerful, and does everything suggested and more. I see > > > > > no argument in support of stuff like the above. > > > > > > > > I have no idea what a PetscSection is and > > > > > > > > PetscSection - This is a mapping from DMMESH points to sets of values, which is > > > > our presentation of a fibre bundle. > > > > > > > > doesn't help explain it. PetscSection is small? It has more methods than DM, that is not small. > > > > > > > > Don't make me return the count. It has fewer, and definitely fewer concepts. Most of those are > > > > simple get/set. > > > > > > > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why does it return an array of integers? The universal way in PETSc would have it returning an IS. Wait a minute, it returns one int, the number of components? Then why is called GetFieldComponents and not GetNumFieldComponents() > > > > > > > > It could be changed to GetNumFieldComponents. > > > > > > > > And what is all this constraint stuff? PetscSectionGetConstraintDof() Huh? It looks inside the "section" and gets an array of ints from a section inside the section? > > > > > > > > Yes. It keeps track of constrained dofs. People, who through laziness or willful blindness, refuse to acknowledge that some > > > > problems make more sense when constrained dofs ar eliminated. > > > > > > > > Here's my guess: a PetscSection is way of "chopping up" a canonical vector into a bunch of fields. Fine, but then why does it have a bunch of other stuff in it (like constraints, what the hell are constraints)? And why isn't a PetscSection just an array of IS that define the fields? > > > > > > > > Constraints are just that. IS is designed to be immutable. Building on it internally is not a good idea for this reason. > > > > > > > > The reason to shove all the index stuff into the DM interface is because generally one doesn't just want to know the entries of a field but one also wants to create vectors for those field variables only, create matrices, know something about the mesh for that field so one can generate interpolations and coarsen for multigrid. One cannot do any of those things with PetscSection. > > > > > > > > Yes you can. > > > > > > Hmm, I guess I haven't pulled recently. I don't see a PetscSectionGetInterpolation() a PetscSectionCoarsen() a PetscSectionGetGlobalVector() etc etc. > > > > > > > > > > > Please explain in freshman linear algebra what a PetscSection is. > > > > > > > > Exactly what it says in the help. It is a mapping from mesh points (vertices, edges, faces, cells) to dofs. How can this be easier? > > > > > > Say I want an IS that represents all the cell center dofs, how do I obtain that from a PetscSection? I don't see the functionality for that but it seems the most basic kind of thing needed. > > > > > > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh points which are cells, height is like co-dimension */ > > > numCellDof = 0; > > > for(c = cStart; c < cEnd; ++c) { > > > PetscSectionGetDof(section, c, &dof); > > > numCellDof += dof; > > > } > > > PetscMalloc(numCellDof * sizeof(PetscInt), &indices); > > > numCellDof = 0; > > > for(c = cStart; c < cEnd; ++c) { > > > PetscSectionGetDof(section, c, &dof); > > > PetscSectionGetOffset(section, c, &off); > > > for(d = 0; d < dof; ++d, ++numCellDof) { > > > indices[numCellDof] = off+d; > > > } > > > } > > > > > > Matt > > > > > > > > > Barry > > > > > > > > > > > Matt > > > > > > > > > > > > Barry > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > Matt > > > > > > > > > > Barry > > > > > > > > > > Note that current DMDA's default subdms served up would be a DM for each field consisting of one of the degrees of freedom per node and it could also serve up DMs that represent any collections of those for example dof 0 and 1 per node, then dof 2 per node. We'll need to come up with some naming convention for these. > > > > > > > > > > > > > > > > Possibly also (eventually) a way to solve equations of state so that we can formulate reduced problems in non-conservative variables. > > > > > > > > > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley(a)gmail.com> wrote: > > > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > > > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a > > > > > > > > > > > > Matt, I want to get rid of these conditionals, not add more. We should have a DM base interface for getting fields on which to split, that common interface should _replace_ the DMComposite specialization. > > > > > > > > > > > > I agree. However, I won't put anything in until I do it by hand once. > > > > > > > > > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its not quite right. I am planning > > > > > > to put the IS method in PetscSection. I am hoping eventually DMDA uses PetscSection for > > > > > > layout. > > > > > > > > > > > > Matt > > > > > > > > > > > > -- > > > > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > > > > -- Norbert Wiener > > > > > > > > > > > > > > > > > > > > > > > > > -- > > > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > > > -- Norbert Wiener > > > > > > > > > > > > > > > > > > > > -- > > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > > -- Norbert Wiener > > > > > > > > > > > > > > > -- > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > -- Norbert Wiener > > > > > > > > > > -- > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > -- Norbert Wiener > > > > > -- > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > -- Norbert Wiener

1 0

Re: [petsc-dev] DMComplex and PCFieldSplit
by Matthew Knepley 24 Feb '12

24 Feb '12

On Fri, Feb 24, 2012 at 4:16 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > I have DMComplexGetDefaultSection(), which is the first thing I call > when constructing the FieldSplit info. > > Found it. So you want PCFieldSplit to troll through those inner loops of > > /* Define fields */ > PetscSection section, sectionGlobal; > PetscInt *fieldSizes, **fieldIndices; > PetscInt nF, f, pStart, pEnd, p; > > ierr = DMComplexGetDefaultSection(pc->dm, &section);CHKERRQ(ierr); > ierr = DMComplexGetDefaultGlobalSection(pc->dm, > &sectionGlobal);CHKERRQ(ierr); > ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); > ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt > *,&fieldIndices);CHKERRQ(ierr); > ierr = PetscSectionGetChart(sectionGlobal, &pStart, > &pEnd);CHKERRQ(ierr); > for(f = 0; f < nF; ++f) { > fieldSizes[f] = 0; > } > for(p = pStart; p < pEnd; ++p) { > PetscInt gdof; > > ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); > if (gdof > 0) { > for(f = 0; f < nF; ++f) { > PetscInt fcdof; > > ierr = PetscSectionGetFieldConstraintDof(section, p, f, > &fcdof);CHKERRQ(ierr); > fieldSizes[f] += fcdof; > } > } > } > for(f = 0; f < nF; ++f) { > ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt), > &fieldIndices[f]);CHKERRQ(ierr); > fieldSizes[f] = 0; > } > for(p = pStart; p < pEnd; ++p) { > PetscInt gdof, goff; > > ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); > if (gdof > 0) { > ierr = PetscSectionGetOffset(sectionGlobal, p, > &goff);CHKERRQ(ierr); > for(f = 0; f < nF; ++f) { > PetscInt fcdof, fc; > > ierr = PetscSectionGetFieldConstraintDof(section, p, f, > &fcdof);CHKERRQ(ierr); > for(fc = 0; fc < fcdof; ++fc, ++fieldSizes[f]) { > fieldIndices[f][fieldSizes[f]] = goff++; > } > } > } > } > for(f = 0; f < nF; ++f) { > IS field; > const char *fieldname; > > ierr = PetscSectionGetFieldName(section, f, > &fieldname);CHKERRQ(ierr); > ierr = ISCreateGeneral(((PetscObject) pc)->comm, fieldSizes[f], > fieldIndices[f], PETSC_OWN_POINTER, &field);CHKERRQ(ierr); > ierr = PetscPrintf(((PetscObject) pc)->comm, "Field %s\n", > fieldname);CHKERRQ(ierr); > ierr = ISView(field, PETSC_NULL);CHKERRQ(ierr); > ierr = PCFieldSplitSetIS(pc, fieldname, field);CHKERRQ(ierr); > ierr = ISDestroy(&field);CHKERRQ(ierr); > } > ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr); > } > > every time instead of call DMGetFieldsIS(dm&is); That's insane, just > take the crud above and call it DMGetFieldsIS_Complex() and we'll have > PCFieldSplit independent of any DM guts. > > I don't understand what we are fighting about. Yes, I agree with that. I will move it. Matt > > On Feb 24, 2012, at 3:49 PM, Matthew Knepley wrote: > > > > > It is the correct level of granularity, especially fr FEM which PETSc > has almost none of right now > > because we make it unnaturally hard due to our ridiculous infatuation > with Finite Difference. That kind > > of thinking colors everything. > > I never said that PetscSection doesn't belong anywhere. It just doesn't > seem to belong anywhere in PC/KSP land, it is unnecessarily low for that. > Sure it may be ok for the finite element assembly etc. > > Barry > > > > > > Regardless, say the concept of PetscSection as a way to keep the > information about fields is fine. PetscSection is a better representation > of an array of ISs. > > > > 1) Imbedding these constraint thingys in it seems completely wrong. > Information about the constraint things should just be in another > PetscSection not inside the PetscSection that tells you the field division. > > > > I disagree since you need this information in every situation where you > ask about sizes, etc. Again, there is only ONE > > problem in PETSc that does this, written by me. It is often useful to > try something before trying to figure out what is the > > best way to do it. > > > > 2) One could then have DMGetPetscSection(DM,section) that tells you the > decomposition, but since the solvers need DMs for the pieces (in > PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so > different from my first email with PetscSection replacing my array of ISs. > So why not fix IS to be a bit more dynamic and build there is no real need > for PetscSection except as an array of ISs. > > > > I have DMComplexGetDefaultSection(), which is the first thing I call > when constructing the FieldSplit info. > > > > Matt > > > > Jed is laughing his pants off at the reaction to his side comment in > that email :-) > > > > > > Barry > > > > On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote: > > > > > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > > > > > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote: > > > > > > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > > > > > > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote: > > > > > > > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > > > > > > > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote: > > > > > > > > > > > We need a richer field interface. Some way for a client of DM to > discover vector and tensor structure. > > > > > > > > > > No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS > *ises,DM *dms); and DMPushSubDMType(DM,"name of a type of splitting like > u,v,p or velocity,p"); DMPopSubDMType(); and > DMGetAvailableSubDMTypes(DM,char **); > > > > > > > > > > Do you need more than that? > > > > > > > > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index > stuff into the DM interface. I think this is so > > > > > unwieldy and complicated. PetscSection is small, powerful, and > does everything suggested and more. I see > > > > > no argument in support of stuff like the above. > > > > > > > > I have no idea what a PetscSection is and > > > > > > > > PetscSection - This is a mapping from DMMESH points to sets of > values, which is > > > > our presentation of a fibre bundle. > > > > > > > > doesn't help explain it. PetscSection is small? It has more > methods than DM, that is not small. > > > > > > > > Don't make me return the count. It has fewer, and definitely fewer > concepts. Most of those are > > > > simple get/set. > > > > > > > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why > does it return an array of integers? The universal way in PETSc would have > it returning an IS. Wait a minute, it returns one int, the number of > components? Then why is called GetFieldComponents and not > GetNumFieldComponents() > > > > > > > > It could be changed to GetNumFieldComponents. > > > > > > > > And what is all this constraint stuff? > PetscSectionGetConstraintDof() Huh? It looks inside the "section" and gets > an array of ints from a section inside the section? > > > > > > > > Yes. It keeps track of constrained dofs. People, who through > laziness or willful blindness, refuse to acknowledge that some > > > > problems make more sense when constrained dofs ar eliminated. > > > > > > > > Here's my guess: a PetscSection is way of "chopping up" a canonical > vector into a bunch of fields. Fine, but then why does it have a bunch of > other stuff in it (like constraints, what the hell are constraints)? And > why isn't a PetscSection just an array of IS that define the fields? > > > > > > > > Constraints are just that. IS is designed to be immutable. Building > on it internally is not a good idea for this reason. > > > > > > > > The reason to shove all the index stuff into the DM interface is > because generally one doesn't just want to know the entries of a field but > one also wants to create vectors for those field variables only, create > matrices, know something about the mesh for that field so one can generate > interpolations and coarsen for multigrid. One cannot do any of those things > with PetscSection. > > > > > > > > Yes you can. > > > > > > Hmm, I guess I haven't pulled recently. I don't see a > PetscSectionGetInterpolation() a PetscSectionCoarsen() a > PetscSectionGetGlobalVector() etc etc. > > > > > > > > > > > Please explain in freshman linear algebra what a PetscSection is. > > > > > > > > Exactly what it says in the help. It is a mapping from mesh points > (vertices, edges, faces, cells) to dofs. How can this be easier? > > > > > > Say I want an IS that represents all the cell center dofs, how do I > obtain that from a PetscSection? I don't see the functionality for that > but it seems the most basic kind of thing needed. > > > > > > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh > points which are cells, height is like co-dimension */ > > > numCellDof = 0; > > > for(c = cStart; c < cEnd; ++c) { > > > PetscSectionGetDof(section, c, &dof); > > > numCellDof += dof; > > > } > > > PetscMalloc(numCellDof * sizeof(PetscInt), &indices); > > > numCellDof = 0; > > > for(c = cStart; c < cEnd; ++c) { > > > PetscSectionGetDof(section, c, &dof); > > > PetscSectionGetOffset(section, c, &off); > > > for(d = 0; d < dof; ++d, ++numCellDof) { > > > indices[numCellDof] = off+d; > > > } > > > } > > > > > > Matt > > > > > > > > > Barry > > > > > > > > > > > Matt > > > > > > > > > > > > Barry > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > Matt > > > > > > > > > > Barry > > > > > > > > > > Note that current DMDA's default subdms served up would be a DM > for each field consisting of one of the degrees of freedom per node and it > could also serve up DMs that represent any collections of those for > example dof 0 and 1 per node, then dof 2 per node. We'll need to come up > with some naming convention for these. > > > > > > > > > > > > > > > > Possibly also (eventually) a way to solve equations of state so > that we can formulate reduced problems in non-conservative variables. > > > > > > > > > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley(a)gmail.com> > wrote: > > > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown(a)mcs.anl.gov> > wrote: > > > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a > > > > > > > > > > > > Matt, I want to get rid of these conditionals, not add more. We > should have a DM base interface for getting fields on which to split, that > common interface should _replace_ the DMComposite specialization. > > > > > > > > > > > > I agree. However, I won't put anything in until I do it by hand > once. > > > > > > > > > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its > not quite right. I am planning > > > > > > to put the IS method in PetscSection. I am hoping eventually > DMDA uses PetscSection for > > > > > > layout. > > > > > > > > > > > > Matt > > > > > > > > > > > > -- > > > > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > > > > -- Norbert Wiener > > > > > > > > > > > > > > > > > > > > > > > > > -- > > > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > > > -- Norbert Wiener > > > > > > > > > > > > > > > > > > > > -- > > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > > -- Norbert Wiener > > > > > > > > > > > > > > > -- > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > -- Norbert Wiener > > > > > > > > > > -- > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > -- Norbert Wiener > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0

Re: [petsc-dev] DMComplex and PCFieldSplit
by Barry Smith 24 Feb '12

24 Feb '12

> I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info. Found it. So you want PCFieldSplit to troll through those inner loops of /* Define fields */ PetscSection section, sectionGlobal; PetscInt *fieldSizes, **fieldIndices; PetscInt nF, f, pStart, pEnd, p; ierr = DMComplexGetDefaultSection(pc->dm, &section);CHKERRQ(ierr); ierr = DMComplexGetDefaultGlobalSection(pc->dm, &sectionGlobal);CHKERRQ(ierr); ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt *,&fieldIndices);CHKERRQ(ierr); ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); for(f = 0; f < nF; ++f) { fieldSizes[f] = 0; } for(p = pStart; p < pEnd; ++p) { PetscInt gdof; ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); if (gdof > 0) { for(f = 0; f < nF; ++f) { PetscInt fcdof; ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); fieldSizes[f] += fcdof; } } } for(f = 0; f < nF; ++f) { ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);CHKERRQ(ierr); fieldSizes[f] = 0; } for(p = pStart; p < pEnd; ++p) { PetscInt gdof, goff; ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); if (gdof > 0) { ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); for(f = 0; f < nF; ++f) { PetscInt fcdof, fc; ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); for(fc = 0; fc < fcdof; ++fc, ++fieldSizes[f]) { fieldIndices[f][fieldSizes[f]] = goff++; } } } } for(f = 0; f < nF; ++f) { IS field; const char *fieldname; ierr = PetscSectionGetFieldName(section, f, &fieldname);CHKERRQ(ierr); ierr = ISCreateGeneral(((PetscObject) pc)->comm, fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &field);CHKERRQ(ierr); ierr = PetscPrintf(((PetscObject) pc)->comm, "Field %s\n", fieldname);CHKERRQ(ierr); ierr = ISView(field, PETSC_NULL);CHKERRQ(ierr); ierr = PCFieldSplitSetIS(pc, fieldname, field);CHKERRQ(ierr); ierr = ISDestroy(&field);CHKERRQ(ierr); } ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr); } every time instead of call DMGetFieldsIS(dm&is); That's insane, just take the crud above and call it DMGetFieldsIS_Complex() and we'll have PCFieldSplit independent of any DM guts. I don't understand what we are fighting about. On Feb 24, 2012, at 3:49 PM, Matthew Knepley wrote: > > It is the correct level of granularity, especially fr FEM which PETSc has almost none of right now > because we make it unnaturally hard due to our ridiculous infatuation with Finite Difference. That kind > of thinking colors everything. I never said that PetscSection doesn't belong anywhere. It just doesn't seem to belong anywhere in PC/KSP land, it is unnecessarily low for that. Sure it may be ok for the finite element assembly etc. Barry > > Regardless, say the concept of PetscSection as a way to keep the information about fields is fine. PetscSection is a better representation of an array of ISs. > > 1) Imbedding these constraint thingys in it seems completely wrong. Information about the constraint things should just be in another PetscSection not inside the PetscSection that tells you the field division. > > I disagree since you need this information in every situation where you ask about sizes, etc. Again, there is only ONE > problem in PETSc that does this, written by me. It is often useful to try something before trying to figure out what is the > best way to do it. > > 2) One could then have DMGetPetscSection(DM,section) that tells you the decomposition, but since the solvers need DMs for the pieces (in PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so different from my first email with PetscSection replacing my array of ISs. So why not fix IS to be a bit more dynamic and build there is no real need for PetscSection except as an array of ISs. > > I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info. > > Matt > > Jed is laughing his pants off at the reaction to his side comment in that email :-) > > > Barry > > On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote: > > > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote: > > > > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote: > > > > > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote: > > > > > > > > > We need a richer field interface. Some way for a client of DM to discover vector and tensor structure. > > > > > > > > No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS *ises,DM *dms); and DMPushSubDMType(DM,"name of a type of splitting like u,v,p or velocity,p"); DMPopSubDMType(); and DMGetAvailableSubDMTypes(DM,char **); > > > > > > > > Do you need more than that? > > > > > > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index stuff into the DM interface. I think this is so > > > > unwieldy and complicated. PetscSection is small, powerful, and does everything suggested and more. I see > > > > no argument in support of stuff like the above. > > > > > > I have no idea what a PetscSection is and > > > > > > PetscSection - This is a mapping from DMMESH points to sets of values, which is > > > our presentation of a fibre bundle. > > > > > > doesn't help explain it. PetscSection is small? It has more methods than DM, that is not small. > > > > > > Don't make me return the count. It has fewer, and definitely fewer concepts. Most of those are > > > simple get/set. > > > > > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why does it return an array of integers? The universal way in PETSc would have it returning an IS. Wait a minute, it returns one int, the number of components? Then why is called GetFieldComponents and not GetNumFieldComponents() > > > > > > It could be changed to GetNumFieldComponents. > > > > > > And what is all this constraint stuff? PetscSectionGetConstraintDof() Huh? It looks inside the "section" and gets an array of ints from a section inside the section? > > > > > > Yes. It keeps track of constrained dofs. People, who through laziness or willful blindness, refuse to acknowledge that some > > > problems make more sense when constrained dofs ar eliminated. > > > > > > Here's my guess: a PetscSection is way of "chopping up" a canonical vector into a bunch of fields. Fine, but then why does it have a bunch of other stuff in it (like constraints, what the hell are constraints)? And why isn't a PetscSection just an array of IS that define the fields? > > > > > > Constraints are just that. IS is designed to be immutable. Building on it internally is not a good idea for this reason. > > > > > > The reason to shove all the index stuff into the DM interface is because generally one doesn't just want to know the entries of a field but one also wants to create vectors for those field variables only, create matrices, know something about the mesh for that field so one can generate interpolations and coarsen for multigrid. One cannot do any of those things with PetscSection. > > > > > > Yes you can. > > > > Hmm, I guess I haven't pulled recently. I don't see a PetscSectionGetInterpolation() a PetscSectionCoarsen() a PetscSectionGetGlobalVector() etc etc. > > > > > > > > Please explain in freshman linear algebra what a PetscSection is. > > > > > > Exactly what it says in the help. It is a mapping from mesh points (vertices, edges, faces, cells) to dofs. How can this be easier? > > > > Say I want an IS that represents all the cell center dofs, how do I obtain that from a PetscSection? I don't see the functionality for that but it seems the most basic kind of thing needed. > > > > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh points which are cells, height is like co-dimension */ > > numCellDof = 0; > > for(c = cStart; c < cEnd; ++c) { > > PetscSectionGetDof(section, c, &dof); > > numCellDof += dof; > > } > > PetscMalloc(numCellDof * sizeof(PetscInt), &indices); > > numCellDof = 0; > > for(c = cStart; c < cEnd; ++c) { > > PetscSectionGetDof(section, c, &dof); > > PetscSectionGetOffset(section, c, &off); > > for(d = 0; d < dof; ++d, ++numCellDof) { > > indices[numCellDof] = off+d; > > } > > } > > > > Matt > > > > > > Barry > > > > > > > > Matt > > > > > > > > > Barry > > > > > > > > > > > > > > > > > > > > > > > > > > Matt > > > > > > > > Barry > > > > > > > > Note that current DMDA's default subdms served up would be a DM for each field consisting of one of the degrees of freedom per node and it could also serve up DMs that represent any collections of those for example dof 0 and 1 per node, then dof 2 per node. We'll need to come up with some naming convention for these. > > > > > > > > > > > > > Possibly also (eventually) a way to solve equations of state so that we can formulate reduced problems in non-conservative variables. > > > > > > > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley(a)gmail.com> wrote: > > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a > > > > > > > > > > Matt, I want to get rid of these conditionals, not add more. We should have a DM base interface for getting fields on which to split, that common interface should _replace_ the DMComposite specialization. > > > > > > > > > > I agree. However, I won't put anything in until I do it by hand once. > > > > > > > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its not quite right. I am planning > > > > > to put the IS method in PetscSection. I am hoping eventually DMDA uses PetscSection for > > > > > layout. > > > > > > > > > > Matt > > > > > > > > > > -- > > > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > > > -- Norbert Wiener > > > > > > > > > > > > > > > > > > > > -- > > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > > -- Norbert Wiener > > > > > > > > > > > > > > > -- > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > -- Norbert Wiener > > > > > > > > > > -- > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > -- Norbert Wiener > > > > > -- > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > -- Norbert Wiener

1 0

Re: [petsc-dev] DMComplex and PCFieldSplit
by Matthew Knepley 24 Feb '12

24 Feb '12

On Fri, Feb 24, 2012 at 3:44 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > What the? Isn't the more nature interface > PetscSectionGetFieldIS(PetscSection, PetscInt fieldnumber, IS *is)? Why are > you providing an interface that provides NO aggregation at all? This way of > getting this information is completely unnatural for 99.9% of users. > It is the correct level of granularity, especially fr FEM which PETSc has almost none of right now because we make it unnaturally hard due to our ridiculous infatuation with Finite Difference. That kind of thinking colors everything. > Regardless, say the concept of PetscSection as a way to keep the > information about fields is fine. PetscSection is a better representation > of an array of ISs. > > 1) Imbedding these constraint thingys in it seems completely wrong. > Information about the constraint things should just be in another > PetscSection not inside the PetscSection that tells you the field division. > I disagree since you need this information in every situation where you ask about sizes, etc. Again, there is only ONE problem in PETSc that does this, written by me. It is often useful to try something before trying to figure out what is the best way to do it. > 2) One could then have DMGetPetscSection(DM,section) that tells you the > decomposition, but since the solvers need DMs for the pieces (in > PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so > different from my first email with PetscSection replacing my array of ISs. > So why not fix IS to be a bit more dynamic and build there is no real need > for PetscSection except as an array of ISs. > I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info. Matt > Jed is laughing his pants off at the reaction to his side comment in that > email :-) > > > Barry > > On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote: > > > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote: > > > > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > > > > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote: > > > > > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > > > > > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote: > > > > > > > > > We need a richer field interface. Some way for a client of DM to > discover vector and tensor structure. > > > > > > > > No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS > *ises,DM *dms); and DMPushSubDMType(DM,"name of a type of splitting like > u,v,p or velocity,p"); DMPopSubDMType(); and > DMGetAvailableSubDMTypes(DM,char **); > > > > > > > > Do you need more than that? > > > > > > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index > stuff into the DM interface. I think this is so > > > > unwieldy and complicated. PetscSection is small, powerful, and does > everything suggested and more. I see > > > > no argument in support of stuff like the above. > > > > > > I have no idea what a PetscSection is and > > > > > > PetscSection - This is a mapping from DMMESH points to sets of values, > which is > > > our presentation of a fibre bundle. > > > > > > doesn't help explain it. PetscSection is small? It has more methods > than DM, that is not small. > > > > > > Don't make me return the count. It has fewer, and definitely fewer > concepts. Most of those are > > > simple get/set. > > > > > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why > does it return an array of integers? The universal way in PETSc would have > it returning an IS. Wait a minute, it returns one int, the number of > components? Then why is called GetFieldComponents and not > GetNumFieldComponents() > > > > > > It could be changed to GetNumFieldComponents. > > > > > > And what is all this constraint stuff? PetscSectionGetConstraintDof() > Huh? It looks inside the "section" and gets an array of ints from a section > inside the section? > > > > > > Yes. It keeps track of constrained dofs. People, who through laziness > or willful blindness, refuse to acknowledge that some > > > problems make more sense when constrained dofs ar eliminated. > > > > > > Here's my guess: a PetscSection is way of "chopping up" a canonical > vector into a bunch of fields. Fine, but then why does it have a bunch of > other stuff in it (like constraints, what the hell are constraints)? And > why isn't a PetscSection just an array of IS that define the fields? > > > > > > Constraints are just that. IS is designed to be immutable. Building on > it internally is not a good idea for this reason. > > > > > > The reason to shove all the index stuff into the DM interface is > because generally one doesn't just want to know the entries of a field but > one also wants to create vectors for those field variables only, create > matrices, know something about the mesh for that field so one can generate > interpolations and coarsen for multigrid. One cannot do any of those things > with PetscSection. > > > > > > Yes you can. > > > > Hmm, I guess I haven't pulled recently. I don't see a > PetscSectionGetInterpolation() a PetscSectionCoarsen() a > PetscSectionGetGlobalVector() etc etc. > > > > > > > > Please explain in freshman linear algebra what a PetscSection is. > > > > > > Exactly what it says in the help. It is a mapping from mesh points > (vertices, edges, faces, cells) to dofs. How can this be easier? > > > > Say I want an IS that represents all the cell center dofs, how do I > obtain that from a PetscSection? I don't see the functionality for that > but it seems the most basic kind of thing needed. > > > > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh > points which are cells, height is like co-dimension */ > > numCellDof = 0; > > for(c = cStart; c < cEnd; ++c) { > > PetscSectionGetDof(section, c, &dof); > > numCellDof += dof; > > } > > PetscMalloc(numCellDof * sizeof(PetscInt), &indices); > > numCellDof = 0; > > for(c = cStart; c < cEnd; ++c) { > > PetscSectionGetDof(section, c, &dof); > > PetscSectionGetOffset(section, c, &off); > > for(d = 0; d < dof; ++d, ++numCellDof) { > > indices[numCellDof] = off+d; > > } > > } > > > > Matt > > > > > > Barry > > > > > > > > Matt > > > > > > > > > Barry > > > > > > > > > > > > > > > > > > > > > > > > > > Matt > > > > > > > > Barry > > > > > > > > Note that current DMDA's default subdms served up would be a DM for > each field consisting of one of the degrees of freedom per node and it > could also serve up DMs that represent any collections of those for > example dof 0 and 1 per node, then dof 2 per node. We'll need to come up > with some naming convention for these. > > > > > > > > > > > > > Possibly also (eventually) a way to solve equations of state so > that we can formulate reduced problems in non-conservative variables. > > > > > > > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley(a)gmail.com> > wrote: > > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown(a)mcs.anl.gov> > wrote: > > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a > > > > > > > > > > Matt, I want to get rid of these conditionals, not add more. We > should have a DM base interface for getting fields on which to split, that > common interface should _replace_ the DMComposite specialization. > > > > > > > > > > I agree. However, I won't put anything in until I do it by hand > once. > > > > > > > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its not > quite right. I am planning > > > > > to put the IS method in PetscSection. I am hoping eventually DMDA > uses PetscSection for > > > > > layout. > > > > > > > > > > Matt > > > > > > > > > > -- > > > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > > > -- Norbert Wiener > > > > > > > > > > > > > > > > > > > > -- > > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > > -- Norbert Wiener > > > > > > > > > > > > > > > -- > > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > > -- Norbert Wiener > > > > > > > > > > -- > > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > > -- Norbert Wiener > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0

Re: [petsc-dev] DMComplex and PCFieldSplit
by Barry Smith 24 Feb '12

24 Feb '12

What the? Isn't the more nature interface PetscSectionGetFieldIS(PetscSection, PetscInt fieldnumber, IS *is)? Why are you providing an interface that provides NO aggregation at all? This way of getting this information is completely unnatural for 99.9% of users. Regardless, say the concept of PetscSection as a way to keep the information about fields is fine. PetscSection is a better representation of an array of ISs. 1) Imbedding these constraint thingys in it seems completely wrong. Information about the constraint things should just be in another PetscSection not inside the PetscSection that tells you the field division. 2) One could then have DMGetPetscSection(DM,section) that tells you the decomposition, but since the solvers need DMs for the pieces (in PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so different from my first email with PetscSection replacing my array of ISs. So why not fix IS to be a bit more dynamic and build there is no real need for PetscSection except as an array of ISs. Jed is laughing his pants off at the reaction to his side comment in that email :-) Barry On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote: > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote: > > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote: > > > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote: > > > > > > > We need a richer field interface. Some way for a client of DM to discover vector and tensor structure. > > > > > > No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS *ises,DM *dms); and DMPushSubDMType(DM,"name of a type of splitting like u,v,p or velocity,p"); DMPopSubDMType(); and DMGetAvailableSubDMTypes(DM,char **); > > > > > > Do you need more than that? > > > > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index stuff into the DM interface. I think this is so > > > unwieldy and complicated. PetscSection is small, powerful, and does everything suggested and more. I see > > > no argument in support of stuff like the above. > > > > I have no idea what a PetscSection is and > > > > PetscSection - This is a mapping from DMMESH points to sets of values, which is > > our presentation of a fibre bundle. > > > > doesn't help explain it. PetscSection is small? It has more methods than DM, that is not small. > > > > Don't make me return the count. It has fewer, and definitely fewer concepts. Most of those are > > simple get/set. > > > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why does it return an array of integers? The universal way in PETSc would have it returning an IS. Wait a minute, it returns one int, the number of components? Then why is called GetFieldComponents and not GetNumFieldComponents() > > > > It could be changed to GetNumFieldComponents. > > > > And what is all this constraint stuff? PetscSectionGetConstraintDof() Huh? It looks inside the "section" and gets an array of ints from a section inside the section? > > > > Yes. It keeps track of constrained dofs. People, who through laziness or willful blindness, refuse to acknowledge that some > > problems make more sense when constrained dofs ar eliminated. > > > > Here's my guess: a PetscSection is way of "chopping up" a canonical vector into a bunch of fields. Fine, but then why does it have a bunch of other stuff in it (like constraints, what the hell are constraints)? And why isn't a PetscSection just an array of IS that define the fields? > > > > Constraints are just that. IS is designed to be immutable. Building on it internally is not a good idea for this reason. > > > > The reason to shove all the index stuff into the DM interface is because generally one doesn't just want to know the entries of a field but one also wants to create vectors for those field variables only, create matrices, know something about the mesh for that field so one can generate interpolations and coarsen for multigrid. One cannot do any of those things with PetscSection. > > > > Yes you can. > > Hmm, I guess I haven't pulled recently. I don't see a PetscSectionGetInterpolation() a PetscSectionCoarsen() a PetscSectionGetGlobalVector() etc etc. > > > > > Please explain in freshman linear algebra what a PetscSection is. > > > > Exactly what it says in the help. It is a mapping from mesh points (vertices, edges, faces, cells) to dofs. How can this be easier? > > Say I want an IS that represents all the cell center dofs, how do I obtain that from a PetscSection? I don't see the functionality for that but it seems the most basic kind of thing needed. > > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh points which are cells, height is like co-dimension */ > numCellDof = 0; > for(c = cStart; c < cEnd; ++c) { > PetscSectionGetDof(section, c, &dof); > numCellDof += dof; > } > PetscMalloc(numCellDof * sizeof(PetscInt), &indices); > numCellDof = 0; > for(c = cStart; c < cEnd; ++c) { > PetscSectionGetDof(section, c, &dof); > PetscSectionGetOffset(section, c, &off); > for(d = 0; d < dof; ++d, ++numCellDof) { > indices[numCellDof] = off+d; > } > } > > Matt > > > Barry > > > > > Matt > > > > > > Barry > > > > > > > > > > > > > > > > > > Matt > > > > > > Barry > > > > > > Note that current DMDA's default subdms served up would be a DM for each field consisting of one of the degrees of freedom per node and it could also serve up DMs that represent any collections of those for example dof 0 and 1 per node, then dof 2 per node. We'll need to come up with some naming convention for these. > > > > > > > > > > Possibly also (eventually) a way to solve equations of state so that we can formulate reduced problems in non-conservative variables. > > > > > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley(a)gmail.com> wrote: > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a > > > > > > > > Matt, I want to get rid of these conditionals, not add more. We should have a DM base interface for getting fields on which to split, that common interface should _replace_ the DMComposite specialization. > > > > > > > > I agree. However, I won't put anything in until I do it by hand once. > > > > > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its not quite right. I am planning > > > > to put the IS method in PetscSection. I am hoping eventually DMDA uses PetscSection for > > > > layout. > > > > > > > > Matt > > > > > > > > -- > > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > > -- Norbert Wiener > > > > > > > > > > > > > > > -- > > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > > -- Norbert Wiener > > > > > > > > > > -- > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > > -- Norbert Wiener > > > > > -- > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > -- Norbert Wiener

1 0