Matt, Okay that makes a lot of sense now, thank you so much. Sorry I have one more question then I will leave you all alone :)
>From what I understand about creating DMPlex, it seems that if there are multiple processes, only the root rank populates the DM structure. I see that DMPlexCreateBoxMesh as well as the GMSH/Exodus/CGNS readers have rank 0 read the mesh files and populate the sieve/mesh points within the DMPlex (assuming that the DM was created with PETSC_COMM_WORLD). However, when I look at the DMPlexCreateFromDAG and DMPlexCreateFromCellList functions, I don't see anywhere where only rank 0 handles the creation of the DM structure. I do notice that DMPlexCreateFromDAG already requires the DM structure to be created. So I guess my question is, is it possible to simply have rank 0 invoke the DMPlexCreateFromDAG function? Because that way only the root process reads the mesh data file and creates/distributes the DMPlex among the remaining processors. Thanks, Justin On Sat, Nov 1, 2014 at 1:30 PM, Matthew Knepley <[email protected]> wrote: > On Sat, Nov 1, 2014 at 12:20 PM, Justin Chang <[email protected]> wrote: > >> Thanks for the response. On a related note, >> >> In the source code of SNES ex12, where is the PetscSection actually being >> created? I can't seem to find anywhere in the code or its routines where >> the creation has been explicitly called. Because when the >> discretization/problem is define, I imagine a PetscSection for the fields >> and constraints has to be invoked somewhere. It seems to me as soon as >> DMPlexAddBoundary() is called it creates the PetscSectionConstraints but I >> can't seem to really confirm this within the source of that function. The >> reason I want to know this is because if I really were to declare my own >> constraints, I want to make sure that the PetscSection for the field >> variables has already been set up. >> > > Much of the mechanism has been moved into the library. The reason I did > this was to make things like > nonlinear multigrid automatic, rather than making the user define a > different section on each level. Let > me explain what is happening. When we call DMGetDefaultSection(), if the > section is not present it > will call the createdefaultsection() member function: > > > https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/interface/dm.c?at=master#cl-2977 > > which for Plex is here > > > https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plex.c?at=master#cl-5601 > > It uses the discretization information from PetscDS and the BC information > from DMPlexAddBoundary() to build the inputs for > DMPlexCreateSection() > > > https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plex.c?at=master#cl-5739 > > I realize that call is still limited. For example, it constrains all dofs > on a point, but sometimes you only want to constrain some. In > the most general case, you would just build the Section manually, using > the same tools I use in DMPlexCreateSection(): > > > https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plex.c?at=master#cl-3204 > > Thanks, > > Matt > > >> Thanks, >> Justin >> >> >> Sent from my iPhone >> >> On Nov 1, 2014, at 11:40 AM, Matthew Knepley <[email protected]> wrote: >> >> On Fri, Oct 31, 2014 at 4:15 PM, Justin Chang <[email protected]> >> wrote: >> >>> Matt, >>> >>> Thanks for the response. One more question: >>> >>> If I go with approach 2, will manually setting the constraints/indices >>> and its values be compatible with the DMPlexSNESComputeResidual/JacobianFEM >>> routines? When I look at the source code of those routines it seems the >>> constrained values are added to the local solution vectors via >>> DMPlexInsertBoundaryValuesFEM. If I choose not to use DMPlexAddBoundary and >>> wish to manually declare my boundary values, i assume I should call >>> DMPlexProjectField with mode INSERT_BC_VALUES? >>> >> >> Yes, exactly. I use DMPlexProjectFunctionLabelLocal(), >> >> >> https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plexfem.c?at=master#cl-576 >> >> Thanks, >> >> Matt >> >> >>> Thanks, >>> Justin >>> >>> Sent from my iPhone >>> >>> On Oct 30, 2014, at 4:24 PM, Matthew Knepley <[email protected]> wrote: >>> >>> On Thu, Oct 30, 2014 at 4:16 PM, Justin Chang <[email protected]> >>> wrote: >>> >>>> Matt, thanks for the quick response. >>>> >>>> What about for (dirichlet) boundary conditions? Would it be possible to >>>> do something similar for those, like using those PetscSectionSetConstraint >>>> functions? >>>> >>> >>> Yes. There are generally two ways of handling Dirichlet conditions: >>> >>> 1) Replace those rows of the Jacobian with the identity and put the >>> boundary value in the rhs. I find >>> this cumbersome for nonlinear solves. >>> >>> 2) Remove these unknowns from the system using >>> PetscSectionSetConstraintDof() and ConstratinIndices(). >>> >>> The DMPlexAddBoundary() functions are automating this processing by >>> marking boundaries using DMLabels, >>> and then constraining dofs on those boundaries. >>> >>> Thanks, >>> >>> Matt >>> >>> >>>> Thanks, >>>> Justin >>>> >>>> On Thu, Oct 30, 2014 at 3:35 PM, Matthew Knepley <[email protected]> >>>> wrote: >>>> >>>>> On Thu, Oct 30, 2014 at 3:29 PM, Justin Chang <[email protected]> >>>>> wrote: >>>>> >>>>>> Hi all, >>>>>> >>>>>> So I am writing an FEM code where it reads input data (e.g., >>>>>> auxiliary coefficients, source terms, etc) from a text file. I have >>>>>> preprocessed the data text file so that each vertex point has its >>>>>> corresponding data. For instance, if my source term for a diffusion >>>>>> problem >>>>>> has a sin or cos function of the coordinates, then this data is already >>>>>> tabulated and simply needs to be fed into my PETSc program. The data text >>>>>> file containing both the auxiliary data and the coordinate/connectivities >>>>>> will also be used to provide the input arguments for the >>>>>> DMPlexCreateFromDAG() function. >>>>>> >>>>>> Normally, I would hard-code the sin/cos functions into the source >>>>>> code itself, but i want my program to be "universal" in that it can take >>>>>> as >>>>>> input any user-defined function for not only these auxiliaries but for >>>>>> other things like boundary conditions. I see that PETSc has a lot of >>>>>> examples on how to read data into vectors and matrices, but I guess my >>>>>> question is is there a way to project data from a text file into a vector >>>>>> that has already been created with a defined DM structure? >>>>>> >>>>> >>>>> If you have the functions available, I think it is far more universal >>>>> to use the function itself, since then you can be independent >>>>> of mesh and discretization when specifying input and BC. >>>>> >>>>> However, if you want to read it in, and you guarantee that it matches >>>>> the mesh and discretization, I think the easiest thing to do is >>>>> demand that it come in the same order as the vertices and use >>>>> >>>>> VecGetArray(V, &a); >>>>> for (v = vStart, i = 0; v < vEnd; ++v) { >>>>> PetscSectionGetDof(s, v, &dof); >>>>> PetscSectionGetOffset(s, v, &off); >>>>> for (d = 0; d < dof; ++d) a[off+d] = text_data[i++]; >>>>> } >>>>> VecRestoreArray(V, &a); >>>>> >>>>> Matt >>>>> >>>>> >>>>>> Thanks, >>>>>> Justin >>>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> 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 >
