[petsc-users] TS for explicit equations
Hello, I'trying to solve an ODE system using PETSc TS through its wrapper petsc4py, and I'm having problems with the setting. The problem is of the kind: u_t = G(t, u) with u being a vector with DoF from a FEM problem and G is computed by a series of matrix operations and manipulations, which are implemented for distributed domains. I set up the problem as explicit using the TSSetRHSFunction, but when I try to solve it using explicit solvers I got a weird behavior: the unknown variable (u above) is never updated (I checked it inside the RHS function I provided), despite having a non null evaluation of f = G(t, u). This happens for RK, SSP, SUNDIALS (adams), but strangely not for FE (where u is updated in every step). That is what disturbs me, because I'd like to use the more sophisticated solvers. I followed this steps in order to get the time evolution of u: # Create and Setup time stepper ts = PETSc.TS().create(comm=self.comm) # Set problem type ts.setProblemType(ts.ProblemType.NONLINEAR) # Set solver type ts.setType(ts.Type.RK) # Set initial time ts.setTime(sTime) # Set max time or max number of iterations ts.setDuration(eTime, max_steps=maxSteps) # Set how to treat the final step ts.setExactFinalTime(ts.ExactFinalTime.INTERPOLATE) # Allocate variables for internal computations of the RHS vector 'U' for the solution vector 'f' for the RHS evaluation # Set rhs function ts.setRHSFunction(self.evalRHSFunction, self.g) (where evalRHSFunction has the form f(self, ts, t, Vort, f)) # Set fuction to be run after each step ts.setPostStep(convergedStepFunction) # Set from other command line options ts.setFromOptions() # Set initial conditions in vector 'U' # Run TS ts.solve(U) Is there anything am I doing wrong? Thanks in advance, Alejandro
[petsc-users] Saving vector with hdf5 viewer
Hello, I am trying to save some field stored in a vector which has values associated with vertexes, edges and cells in a DMPlex. This vector was created (using petsc4py) from a petcs section, setting this as default section and then creating the vector from the DMPlex. The thing is that when I save that vector using the hdf5 viewer I found 2 entities inside it. One corresponding to the values for all DoF in the DMPlex (under the group fields), and another with values only for nodes (under the group: vertex_fields). That meaning a useless duplication of data, is it possible to set the viewer only to store one of those groups? Preferably the complete one, but in the future only the vertex one could be usefull. Thanks in advance, Alejandro
Re: [petsc-users] Distributing data along with DMPlex
Thanks for your answer. In principle, this is what I want, although I know that this doesn't scale. At this point I am reading the mesh in only one proc so I have all the properties there, that is why after that I need to distribute the mesh along with the data. Is there a better solution within this scheme (only one proc reads mesh data)? Regards, Alejandro On Tue, Dec 22, 2015 at 1:11 AM, Matthew Knepley wrote: > On Thu, Dec 17, 2015 at 1:27 PM, Justin Chang wrote: > >> So you would use something like DMPlexDistributeField() >> <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMPlexDistributeField.html> >> in >> that case. You have your original/global DM, PetscSection, and Vec of >> values. Then using the PetscSF that was created during DMPlexDistribute, >> you map the global stuff into the local stuff. >> > > Justin is correct. However, I will note that this strategy is not > scalable. Usually the data is much larger than the mesh, > so this runs out of memory much sooner. Are you sure that this is what you > want? > > Thanks, > > Matt > > >> On Thu, Dec 17, 2015 at 11:21 AM, Alejandro D Otero >> wrote: >> >>> Thanks, >>> The problem then is that after DMPlexDistribute the DMPlex 'points' are >>> renumbered. So if the values are related to each point in the original >>> numbering how do I set the values after the distribution. I know the >>> property stored in the vector related to the entities with the numbering of >>> the original mesh which I use to create the first DMPlex. >>> >>> Ideally for me, I would like to set the values in the vector before >>> DMPlexDistribute and get the vector components renumbered and redistributed >>> accordingly in a global vector. And then, get the local vector. >>> >>> Hope it could be more clear now. >>> Regards, >>> >>> Alejandro >>> >>> >>> >>> On Wed, Dec 16, 2015 at 7:01 PM, Justin Chang >>> wrote: >>> >>>> I think you would follow this order: >>>> >>>> 1*) create a DMPlex (depth, chart, etc) on rank 0. Other ranks have an >>>> empty DM >>>> >>>> 2) DMPlexDistribute() >>>> >>>> 3*) Create the PetscSection >>>> >>>> 4) DMCreateGlobalVector() >>>> >>>> 5) DMCreateLocalVector() >>>> >>>> Now you have a global vector and a local vector for your distributed >>>> DMPlex. The mapping/ghosting/etc of dofs is already taken care of. >>>> >>>> * if you're using standard Galerkin FE then in SNES examples 12 and 62 >>>> (and maybe others?) the first step is handled through the mesh generation >>>> functions and step 3 is handled through step 4 >>>> >>>> Thanks, >>>> Justin >>>> >>>> On Wednesday, December 16, 2015, Alejandro D Otero >>>> wrote: >>>> >>>>> Hi, I need some help understanding how to distribute data together >>>>> with a dmplex representing a FE mesh. >>>>> At the beginning I define the structure of the dmplex assigning >>>>> certain number of DoF to cells, edges and vertexes, in one process (the >>>>> dmplex in the rest is empty) >>>>> I create a petscsecton and I create an associated global vector with >>>>> the quantities I want to store. >>>>> Then I distribute the dmplex over all the processes. >>>>> * Although this does not perform well it is just a starting point. I >>>>> know it has to be improved. >>>>> >>>>> I would like to have the global vector distributed accordingly so that >>>>> each process has access to the corresponding local part with its DoF >>>>> (possibly adding some ghost values corresponding to the shared DoF not >>>>> taken care by it). >>>>> >>>>> Is there any 'correct' way to do that in PETSc? >>>>> >>>>> Thanks in advance, >>>>> >>>>> Alejandro >>>>> >>>> >>> >> > > > -- > 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 >
Re: [petsc-users] Distributing data along with DMPlex
Thanks, The problem then is that after DMPlexDistribute the DMPlex 'points' are renumbered. So if the values are related to each point in the original numbering how do I set the values after the distribution. I know the property stored in the vector related to the entities with the numbering of the original mesh which I use to create the first DMPlex. Ideally for me, I would like to set the values in the vector before DMPlexDistribute and get the vector components renumbered and redistributed accordingly in a global vector. And then, get the local vector. Hope it could be more clear now. Regards, Alejandro On Wed, Dec 16, 2015 at 7:01 PM, Justin Chang wrote: > I think you would follow this order: > > 1*) create a DMPlex (depth, chart, etc) on rank 0. Other ranks have an > empty DM > > 2) DMPlexDistribute() > > 3*) Create the PetscSection > > 4) DMCreateGlobalVector() > > 5) DMCreateLocalVector() > > Now you have a global vector and a local vector for your distributed > DMPlex. The mapping/ghosting/etc of dofs is already taken care of. > > * if you're using standard Galerkin FE then in SNES examples 12 and 62 > (and maybe others?) the first step is handled through the mesh generation > functions and step 3 is handled through step 4 > > Thanks, > Justin > > On Wednesday, December 16, 2015, Alejandro D Otero > wrote: > >> Hi, I need some help understanding how to distribute data together with a >> dmplex representing a FE mesh. >> At the beginning I define the structure of the dmplex assigning certain >> number of DoF to cells, edges and vertexes, in one process (the dmplex in >> the rest is empty) >> I create a petscsecton and I create an associated global vector with the >> quantities I want to store. >> Then I distribute the dmplex over all the processes. >> * Although this does not perform well it is just a starting point. I know >> it has to be improved. >> >> I would like to have the global vector distributed accordingly so that >> each process has access to the corresponding local part with its DoF >> (possibly adding some ghost values corresponding to the shared DoF not >> taken care by it). >> >> Is there any 'correct' way to do that in PETSc? >> >> Thanks in advance, >> >> Alejandro >> >
[petsc-users] Distributing data along with DMPlex
Hi, I need some help understanding how to distribute data together with a dmplex representing a FE mesh. At the beginning I define the structure of the dmplex assigning certain number of DoF to cells, edges and vertexes, in one process (the dmplex in the rest is empty) I create a petscsecton and I create an associated global vector with the quantities I want to store. Then I distribute the dmplex over all the processes. * Although this does not perform well it is just a starting point. I know it has to be improved. I would like to have the global vector distributed accordingly so that each process has access to the corresponding local part with its DoF (possibly adding some ghost values corresponding to the shared DoF not taken care by it). Is there any 'correct' way to do that in PETSc? Thanks in advance, Alejandro
Re: [petsc-users] [Posible SPAM] Re: DMPlex for high order elements
Ok, I'll think about it. Thank you! On Wed, Nov 18, 2015 at 4:29 PM, Jed Brown wrote: > Alejandro D Otero writes: > > > Hi, I am trying to define a spectral element mesh in the frame of DMPlex > > structure. > > To begin with I am thinking of quadrilateral 2d elements. Lets say I > have a > > 16 node element with 2 nodes on each edge besides the 2 corners and 4 > nodes > > inside the element. I'd like to treat all nodes as vertex in the DMPlex > > structure. > > That's not the intended way to use DMPlex. You associate an arbitrary > number of dofs with each topological entity (vertex, edge, face, cell). > So for your Q_3 elements in 2D, you would put 1 dof at vertices, 2 at > edges, and 4 at cell interiors. >
[petsc-users] DMPlex for high order elements
Hi, I am trying to define a spectral element mesh in the frame of DMPlex structure. To begin with I am thinking of quadrilateral 2d elements. Lets say I have a 16 node element with 2 nodes on each edge besides the 2 corners and 4 nodes inside the element. I'd like to treat all nodes as vertex in the DMPlex structure. I am wondering whether I can make internal edge nodes (those not in the extremes of the edge) depend on the corresponding edge and internal nodes (those not on the element edges) depend on the corresponding element? If this is not the right way of doing this from DMPlex philosophical point of view, could you please give me an example? Thanks in advance, best regards, Alejandro
[petsc-users] Question about Nested matrices
Hi all, I am trying to nest to square matrices to form a rectangular one with dims 2n x n (with n the size of the square matrices) I want the two matrices to have interlaced rows so I passed the following 2 list of indices to MatCreateNest: [[ 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48] [ 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49]] And we get the following error message: Global sizes (25,25) of nested submatrix (0,0) do not agree with space defined by index sets (50,25) The size of the 2 square matrices is 25 x 25. What am I doing wrong? Thanks in advance,
[petsc-users] Building a rectangular MPI matrix relating to PetscSections
Hi all, is it posible to build a rectangular matrix that relates two fields with (posible) different number of DoF per point defining the nonzero structure from 2 petscsections representing one field each? As an example, I want to build an operator that given a 2D velocity computes the corresponding vorticity by means of: W = Bcurl * V with V the vector with the velocity components at grid points and W the vorticity at the same points. Thanks in advance, regards, Alejandro
Re: [petsc-users] Nodes coordinates in distributed dmplex
Hi, I seems to be having problems with PETSc vector VTK viewer. The output file I get has some values which are different from what I actually have. I have a scalar field represented by a vector. If I print it on screen I got the following output. Vec Object: 2 MPI processes type: mpi Process [0] 0 0 0.25 0.25 0.5 0.5 0.75 0.75 1 1 Process [1] 0 0 0 0.25 0.25 0.25 0.5 0.5 0.5 0.75 0.75 0.75 1 1 1 But when I set the viewer to use VTK_VTU (in this case ASCII just for make it readable, but the same in native vtk format when read with paraview) I got: 0.00e+00 0.00e+00 2.50e-01 2.50e-01 5.00e-01 5.00e-01 7.155608e-01 7.198046e-01 7.892083e-01 1.00e+00 0.00e+00 2.248979e-01 2.439858e-01 2.788163e-01 2.50e-01 2.50e-01 4.997506e-01 5.042740e-01 5.003098e-01 7.50e-01 7.50e-01 7.50e-01 1.00e+00 1.00e+00 1.00e+00 The difference start in the 7th value on.Any clue of what could be happening? Thanks in advance, Alejandro
Re: [petsc-users] Nodes coordinates in distributed dmplex
Thanks Matt.
It is much clearer for me now.
Regards,
Alejandro
On Fri, Jul 10, 2015 at 6:14 PM, Matthew Knepley wrote:
> On Fri, Jul 10, 2015 at 6:48 AM, Alejandro D Otero
> wrote:
>
>> Hi Matt, thanks for your answer.
>>
>> I got a vector from getCoordinates(). How are their components indexed?
>> is (p * dim + d) with p: node, dim: of the problem, x_d coordinate of the
>> node, ok?
>> Which numbering for p? The local number of node, the number of point in
>> the DAG of the dm, the original number of node?
>>
>
> All data layouts for Plex are described by a PetscSection (which is how it
> should be in the rest of PETSc as well). The PetscSection is a map
> from mesh points to (# dof, offset) pairs. Thus, the way I get coordinates
> is the following:
>
> DMGetCoordinatesLocal(dm, &coordinates);
> DMGetCoordinateDM(dm, &cdm);
> DMGetDefaultSection(cdm, &cs);
> PetscSectionGetChart(cs, &pStart, &pEnd);
> VecGetArrayRead(coordinates, &coords);
> for (p = pStart; p < pEnd; ++p) {
> PetscInt dof, off;
>
> PetscSectionGetDof(cs, p, &dof);
> PetscSectionGetOffset(cs, p, &off);
> for (d = 0; d < dof; ++d)
> }
> VecRestoreArrayRead(coordinates, &coords);
>
>
>> I am trying a simple square mesh with 16 4-node square elements parted
>> into 2 process. Total of 25 nodes.
>> The distributed dm seems alright to me. Each process gets a dm with 8
>> elements an 15 nodes, which means that the 5 shared nodes are local to each
>> process. But one of the process gives negative values for the shared nodes.
>> How need them to be mapped to get the right number.
>>
>
> Where do you get negative point numbers? I encode off-process point
> numbers as -(remote point + 1).
>
>
>> It seems I'm using a wrong approach to this. May be I need to get the
>> coordinates in a somehow different way. I'm working on a from-scratch
>> implementation of a FEM code based on petsc4py. I want to code the problem
>> matrices assembly from elemental matrices. I've already done this
>> sequentially, but I got stuck when trying to compute elemental matrices in
>> parallel because I don't understand the way of obtaining the coordinates of
>> the nodes in for each element.
>>
>
> The way I get coordinates for the element c is
>
> PetscScalar *coords = NULL;
> PetscIntcsize;
>
> ierr = DMPlexVecGetClosure(dm, cs, coordinates, c, &csize,
> &coords);CHKERRQ(ierr);
> ierr = DMPlexVecRestoreClosure(dm, cs, coordinates, c, &csize,
> &coords);CHKERRQ(ierr);
>
> Thanks,
>
> Matt
>
>
>> Again, thanks for your help,
>>
>> A.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> On Thu, Jul 9, 2015 at 5:18 PM, Matthew Knepley
>> wrote:
>>
>>> On Thu, Jul 9, 2015 at 7:42 AM, Alejandro D Otero
>>> wrote:
>>>
>>>> Hi, sorry if this is an obvious question, but I cannot figure out how
>>>> to recover finite element nodes coordinates once I have distributed a mesh
>>>> stored as a dmplex. I am using petsc4py as interface to petsc rutines.
>>>>
>>>> I first created a dmplex using:
>>>> dm.createFromCellList()
>>>>
>>>> In a sequential run I got the coordinates with:
>>>> Coords = dm.getCoordinates()
>>>>
>>>> which gave a sequential vector with the coordinates of the mesh nodes.
>>>>
>>>> When I distribute the mesh with:
>>>> dm.distribute()
>>>>
>>>> each mpi process has it own dm but the indexing of the vector resulting
>>>> from getCoordinates() or getCoordinatesLocal() seems not consistent with
>>>> the local numbering of the cells and nodes.
>>>>
>>>
>>> When the mesh is distributed, the vertices are renumbered. Thus the
>>> coordinates you get out are
>>> for reordered local vertices, but they are consistent with the local
>>> topology (cells still contain the
>>> right vertices) and the overlap mapping (SF still connects the shared
>>> vertices).
>>>
>>> What do you need it to do?
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>
>>>> Which is the correct way of doing this in PETSc philosophy?
>>>>
>>>> Thanks in advance,
>>>> Alejandro
>>>>
>>>
>>>
>>>
>>> --
>>> 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
>
Re: [petsc-users] Nodes coordinates in distributed dmplex
Hi Matt, thanks for your answer. I got a vector from getCoordinates(). How are their components indexed? is (p * dim + d) with p: node, dim: of the problem, x_d coordinate of the node, ok? Which numbering for p? The local number of node, the number of point in the DAG of the dm, the original number of node? I am trying a simple square mesh with 16 4-node square elements parted into 2 process. Total of 25 nodes. The distributed dm seems alright to me. Each process gets a dm with 8 elements an 15 nodes, which means that the 5 shared nodes are local to each process. But one of the process gives negative values for the shared nodes. How need them to be mapped to get the right number. It seems I'm using a wrong approach to this. May be I need to get the coordinates in a somehow different way. I'm working on a from-scratch implementation of a FEM code based on petsc4py. I want to code the problem matrices assembly from elemental matrices. I've already done this sequentially, but I got stuck when trying to compute elemental matrices in parallel because I don't understand the way of obtaining the coordinates of the nodes in for each element. Again, thanks for your help, A. On Thu, Jul 9, 2015 at 5:18 PM, Matthew Knepley wrote: > On Thu, Jul 9, 2015 at 7:42 AM, Alejandro D Otero > wrote: > >> Hi, sorry if this is an obvious question, but I cannot figure out how to >> recover finite element nodes coordinates once I have distributed a mesh >> stored as a dmplex. I am using petsc4py as interface to petsc rutines. >> >> I first created a dmplex using: >> dm.createFromCellList() >> >> In a sequential run I got the coordinates with: >> Coords = dm.getCoordinates() >> >> which gave a sequential vector with the coordinates of the mesh nodes. >> >> When I distribute the mesh with: >> dm.distribute() >> >> each mpi process has it own dm but the indexing of the vector resulting >> from getCoordinates() or getCoordinatesLocal() seems not consistent with >> the local numbering of the cells and nodes. >> > > When the mesh is distributed, the vertices are renumbered. Thus the > coordinates you get out are > for reordered local vertices, but they are consistent with the local > topology (cells still contain the > right vertices) and the overlap mapping (SF still connects the shared > vertices). > > What do you need it to do? > > Thanks, > > Matt > > >> Which is the correct way of doing this in PETSc philosophy? >> >> Thanks in advance, >> Alejandro >> > > > > -- > 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 >
[petsc-users] Nodes coordinates in distributed dmplex
Hi, sorry if this is an obvious question, but I cannot figure out how to recover finite element nodes coordinates once I have distributed a mesh stored as a dmplex. I am using petsc4py as interface to petsc rutines. I first created a dmplex using: dm.createFromCellList() In a sequential run I got the coordinates with: Coords = dm.getCoordinates() which gave a sequential vector with the coordinates of the mesh nodes. When I distribute the mesh with: dm.distribute() each mpi process has it own dm but the indexing of the vector resulting from getCoordinates() or getCoordinatesLocal() seems not consistent with the local numbering of the cells and nodes. Which is the correct way of doing this in PETSc philosophy? Thanks in advance, Alejandro
