Re: [petsc-users] SNESSetJacobian
On 09/23/2014 06:11 PM, Barry Smith wrote: On Sep 23, 2014, at 9:50 AM, anton po...@uni-mainz.de wrote: Starting from version 3.5 the matrix parameters in SNESSetJacobian are no longer pointers, hence my question: What is the most appropriate place to call SNESSetJacobian if I need to change the Jacobian during solution? What about FormFunction? Could you please explain why you need to change the Mat? Our hope was that people would not need to change it. Note that you can change the type of a matrix at any time. So for example inside your FormJacobian you can have code like MatSetType(J,MATAIJ) this wipes out the old matrix data structure and gives you an empty matrix of the new type ready to be preallocated and then filled. Let us know what you need. How should a user switch from assembled to a matrix-free Jacobian (for example) within one run? Simplest is resetting SNES altogether, I guess. Anton Barry Thanks, Anton
Re: [petsc-users] SNESSetJacobian
On Wed, Sep 24, 2014 at 6:21 AM, anton po...@uni-mainz.de wrote: On 09/23/2014 06:11 PM, Barry Smith wrote: On Sep 23, 2014, at 9:50 AM, anton po...@uni-mainz.de wrote: Starting from version 3.5 the matrix parameters in SNESSetJacobian are no longer pointers, hence my question: What is the most appropriate place to call SNESSetJacobian if I need to change the Jacobian during solution? What about FormFunction? Could you please explain why you need to change the Mat? Our hope was that people would not need to change it. Note that you can change the type of a matrix at any time. So for example inside your FormJacobian you can have code like MatSetType(J,MATAIJ) this wipes out the old matrix data structure and gives you an empty matrix of the new type ready to be preallocated and then filled. Let us know what you need. How should a user switch from assembled to a matrix-free Jacobian (for example) within one run? Simplest is resetting SNES altogether, I guess. Set the type to MATSHELL and set your apply function. You still should not need to change the pointer, exactly as Barry says above. Matt Anton Barry Thanks, Anton -- 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] SNESSetJacobian
On Sep 24, 2014, at 5:21 AM, anton po...@uni-mainz.de wrote: On 09/23/2014 06:11 PM, Barry Smith wrote: On Sep 23, 2014, at 9:50 AM, anton po...@uni-mainz.de wrote: Starting from version 3.5 the matrix parameters in SNESSetJacobian are no longer pointers, hence my question: What is the most appropriate place to call SNESSetJacobian if I need to change the Jacobian during solution? What about FormFunction? Could you please explain why you need to change the Mat? Our hope was that people would not need to change it. Note that you can change the type of a matrix at any time. So for example inside your FormJacobian you can have code like MatSetType(J,MATAIJ) this wipes out the old matrix data structure and gives you an empty matrix of the new type ready to be preallocated and then filled. Let us know what you need. How should a user switch from assembled to a matrix-free Jacobian (for example) within one run? Simplest is resetting SNES altogether, I guess. So you want to run, say, three steps of Newton “matrix-free” and then four steps with an explicit matrix? Then I guess you could call SNESSetJacobian() inside a SNES monitor routine, or even in your compute Jacobian routine. Calling it within FormFunction would not be a good idea since FormFunction is called in a variety of places. Barry Anton Barry Thanks, Anton
Re: [petsc-users] DMPlex with spring elements
Thanks for your response. My discretization is based on spring elements. For the linear one dimensional case in which each spring has a coefficient k, their jacobian would be this two by two matrix. [ k-k ] [ -k k ] and the internal force [ k ( Ui - Uj) ] [ k ( Uj - Ui) ] where Ui and Uj are the node displacements (just one displacement per node because it's one dimensional) For the two dimensional case, assuming small deformations, we have a four-by-four matrix. Each node has two degrees of freedom. We obtain it by performing the outer product of the vector (t , -t) where t is the vector that connects both nodes in a spring. This is for the case of small deformations. I would need to assemble each spring contribution to the jacobian and the residual like they were finite elements. The springs share nodes, that's how they are connected. This example is just the linear case, I will have to implement a nonlinear case in a similar fashion. Seeing the DMNetwork example, I think it's what I need, although I don't know much of power electric grids and it's hard for me to understand what's going on. Do you have a good reference to be able to follow the code? For example, why are they adding components to the edges? 475: DMNetworkAddComponent http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/docs/manualpages/DM/DMNetworkAddComponent.html#DMNetworkAddComponent(networkdm,i,componentkey[0],pfdata.branch[i-eStart]); Miguel On Tue, Sep 23, 2014 at 11:13 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You may also want to take a look at the DMNetwork framework that can be used for general unstructured networks that don't use PDEs. Its description is given in the manual and an example is in src/snes/examples/tutorials/network/pflow. Shri From: Matthew Knepley knep...@gmail.com Date: Tue, 23 Sep 2014 22:40:52 -0400 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Tue, Sep 23, 2014 at 4:01 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi all I was wondering if it could be possible to build a model similar to the example snes/ex12.c, but with spring elements (for elasticity) instead of simplicial elements. Spring elements in a grid, therefore each element would have two nodes and each node two components. There would be more differences, because instead of calling the functions f0,f1,g0,g1,g2 and g3 to build the residual and the jacobian, I would call a routine that would build the residual vector and the jacobian matrix directly. I would not have shape functions whatsoever. My problem is discrete, I don't have a PDE and my equations are algebraic. What is the best way in petsc to solve this problem? Is there any example that I can follow? Thanks in advance Yes, ex12 is fairly specific to FEM. However, I think the right tools for what you want are DMPlex and PetscSection. Here is how I would proceed: 1) Make a DMPlex that encodes a simple network that you wish to simulate 2) Make a PetscSection that gets the data layout right. Its hard from the above for me to understand where you degrees of freedom actually are. This is usually the hard part. 3) Calculate the residual, so you can check an exact solution. Here you use the PetscSectionGetDof/Offset() for each mesh piece that you are interested in. Again, its hard to be more specific when I do not understand your discretization. Thanks, Matt Miguel -- Miguel Angel Salazar de Troya Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 tel:%28217%29%20550-2360 salaz...@illinois.edu -- 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 -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu
Re: [petsc-users] DMPlex with spring elements
On Wed, Sep 24, 2014 at 11:31 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks for your response. My discretization is based on spring elements.
For the linear one dimensional case in which each spring has a coefficient
k, their jacobian would be this two by two matrix.
[ k-k ]
[ -k k ]
and the internal force
[ k ( Ui - Uj) ]
[ k ( Uj - Ui) ]
where Ui and Uj are the node displacements (just one displacement per node
because it's one dimensional)
For the two dimensional case, assuming small deformations, we have a
four-by-four matrix. Each node has two degrees of freedom. We obtain it by
performing the outer product of the vector (t , -t) where t is the vector
that connects both nodes in a spring. This is for the case of small
deformations. I would need to assemble each spring contribution to the
jacobian and the residual like they were finite elements. The springs share
nodes, that's how they are connected. This example is just the linear case,
I will have to implement a nonlinear case in a similar fashion.
Seeing the DMNetwork example, I think it's what I need, although I don't
know much of power electric grids and it's hard for me to understand what's
going on. Do you have a good reference to be able to follow the code? For
example, why are they adding components to the edges?
Okay, I understand what you want now. I would recommend doing a simple
example by hand first. The PetscSection is very simple,
just 'dim' degrees of freedom for each vertex. It should look something like
DMPlexGetDepthStratum(dm, 0, vStart, vEnd);
PetscSectionCreate(comm, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetChart(s, vStart, vEnd);
for (v = vStart; v vEnd; ++v) {
PetscSectionSetDof(s, v, dim);
}
PetscSectionSetUp(s);
DMSetDefaultSection(dm, s);
PetscSectionDestroy(s);
Then when you form the residual, you loop over edges and add a contribution
to each vertex (I think, its still not
clear to me how your residual is defined)
VecGetArray(locF, residual);
DMPlexGetHeightStratum(dm, 0, cStart, cEnd);
for (c = cStart; c cEnd; ++c) {
const PetscInt *cone;
PetscInt offA, offB;
DMPlexGetCone(c, cone);
PetscSectionGetOffset(s, cone[0], offA);
PetscSectionGetOffset(s, cone[1], offB);
residual[offA] = contribution for vertex A;
residual[offB] = contribution for vertex B;
}
VecRestoreArray(locF, residual);
do LocalToGlobal unless PETSc is doing it for you
After that works, using DMNetwork should be much easier.
Thanks,
Matt
475: DMNetworkAddComponent
http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/docs/manualpages/DM/DMNetworkAddComponent.html#DMNetworkAddComponent(networkdm,i,componentkey[0],pfdata.branch[i-eStart]);
Miguel
On Tue, Sep 23, 2014 at 11:13 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
You may also want to take a look at the DMNetwork framework that can be
used for general unstructured networks that don't use PDEs. Its
description is given in the manual and an example is in
src/snes/examples/tutorials/network/pflow.
Shri
From: Matthew Knepley knep...@gmail.com
Date: Tue, 23 Sep 2014 22:40:52 -0400
To: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
On Tue, Sep 23, 2014 at 4:01 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Hi all
I was wondering if it could be possible to build a model similar to the
example snes/ex12.c, but with spring elements (for elasticity) instead of
simplicial elements. Spring elements in a grid, therefore each element
would have two nodes and each node two components. There would be more
differences, because instead of calling the functions f0,f1,g0,g1,g2 and
g3 to build the residual and the jacobian, I would call a routine that
would build the residual vector and the jacobian matrix directly. I would
not have shape functions whatsoever. My problem is discrete, I don't have
a PDE and my equations are algebraic. What is the best way in petsc to
solve this problem? Is there any example that I can follow? Thanks in
advance
Yes, ex12 is fairly specific to FEM. However, I think the right tools for
what you want are
DMPlex and PetscSection. Here is how I would proceed:
1) Make a DMPlex that encodes a simple network that you wish to
simulate
2) Make a PetscSection that gets the data layout right. Its hard from
the above
for me to understand where you degrees of freedom actually are.
This is usually
the hard part.
3) Calculate the residual, so you can check an exact solution. Here you
use the
PetscSectionGetDof/Offset() for each mesh piece that you are
interested in. Again,
its hard to be more specific when I do not understand your
discretization.
Thanks,
Matt
Miguel
--
Miguel Angel Salazar de Troya
Graduate Research Assistant
[petsc-users] How to set L1 norm as the converge test?
Hi all, For example, I am using LSQR or CG to solve Ax=b. When we turn on the monitor option of solving the linear system, we could see iteration number and residual for each iteration. I believe it's L2 norm of the objective function? If yes, is there a way that I could set it to L1 norm for solving the system? Thanks. -- Brian Yang U of Houston
Re: [petsc-users] How to set L1 norm as the converge test?
On Wed, Sep 24, 2014 at 2:08 PM, Brian Yang brianyang1...@gmail.com wrote: Hi all, For example, I am using LSQR or CG to solve Ax=b. When we turn on the monitor option of solving the linear system, we could see iteration number and residual for each iteration. I believe it's L2 norm of the objective function? If yes, is there a way that I could set it to L1 norm for solving the system? Yes, you can use http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetConvergenceTest.html and then use VecNorm() with NORM_1 instead. Thanks, Matt Thanks. -- Brian Yang U of Houston -- 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] How to set L1 norm as the converge test?
Brian, Do you wish to monitor the convergence of the L1 norm, or do you wish to use the L1 norm in the convergence test? Regardless you cannot change the LSQR or CG to be a different algorithm (with different convergence properties) with the L1 norm. Barry On Sep 24, 2014, at 2:08 PM, Brian Yang brianyang1...@gmail.com wrote: Hi all, For example, I am using LSQR or CG to solve Ax=b. When we turn on the monitor option of solving the linear system, we could see iteration number and residual for each iteration. I believe it's L2 norm of the objective function? If yes, is there a way that I could set it to L1 norm for solving the system? Thanks. -- Brian Yang U of Houston
Re: [petsc-users] How to set L1 norm as the converge test?
Thanks Mat and Barry, Yes I want to use the L1 norm in the convergence test. For the ||Ax-b||, I always want to calculate L1 norm instead of L2. Is Mat's way gonna work on this? On Wed, Sep 24, 2014 at 2:17 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Do you wish to monitor the convergence of the L1 norm, or do you wish to use the L1 norm in the convergence test? Regardless you cannot change the LSQR or CG to be a different algorithm (with different convergence properties) with the L1 norm. Barry On Sep 24, 2014, at 2:08 PM, Brian Yang brianyang1...@gmail.com wrote: Hi all, For example, I am using LSQR or CG to solve Ax=b. When we turn on the monitor option of solving the linear system, we could see iteration number and residual for each iteration. I believe it's L2 norm of the objective function? If yes, is there a way that I could set it to L1 norm for solving the system? Thanks. -- Brian Yang U of Houston -- Brian Yang U of Houston
Re: [petsc-users] How to set L1 norm as the converge test?
Brian, Your convergence test routine could call KSPBuildResidual() and then compute the 1-norm of the residual and make any decision it likes based on that norm. See KSPSetConvergenceTest() Barry On Sep 24, 2014, at 2:24 PM, Brian Yang brianyang1...@gmail.com wrote: Thanks Mat and Barry, Yes I want to use the L1 norm in the convergence test. For the ||Ax-b||, I always want to calculate L1 norm instead of L2. Is Mat's way gonna work on this? On Wed, Sep 24, 2014 at 2:17 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Do you wish to monitor the convergence of the L1 norm, or do you wish to use the L1 norm in the convergence test? Regardless you cannot change the LSQR or CG to be a different algorithm (with different convergence properties) with the L1 norm. Barry On Sep 24, 2014, at 2:08 PM, Brian Yang brianyang1...@gmail.com wrote: Hi all, For example, I am using LSQR or CG to solve Ax=b. When we turn on the monitor option of solving the linear system, we could see iteration number and residual for each iteration. I believe it's L2 norm of the objective function? If yes, is there a way that I could set it to L1 norm for solving the system? Thanks. -- Brian Yang U of Houston -- Brian Yang U of Houston
Re: [petsc-users] How to set L1 norm as the converge test?
On Wed, Sep 24, 2014 at 2:31 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Your convergence test routine could call KSPBuildResidual() and then compute the 1-norm of the residual and make any decision it likes based on that norm. See KSPSetConvergenceTest() Note that what Barry is saying is that the convergence theory for CG guarantees monotonicity in the energy (A) norm, but says nothing about L1 so you might get crap. Matt Barry On Sep 24, 2014, at 2:24 PM, Brian Yang brianyang1...@gmail.com wrote: Thanks Mat and Barry, Yes I want to use the L1 norm in the convergence test. For the ||Ax-b||, I always want to calculate L1 norm instead of L2. Is Mat's way gonna work on this? On Wed, Sep 24, 2014 at 2:17 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Do you wish to monitor the convergence of the L1 norm, or do you wish to use the L1 norm in the convergence test? Regardless you cannot change the LSQR or CG to be a different algorithm (with different convergence properties) with the L1 norm. Barry On Sep 24, 2014, at 2:08 PM, Brian Yang brianyang1...@gmail.com wrote: Hi all, For example, I am using LSQR or CG to solve Ax=b. When we turn on the monitor option of solving the linear system, we could see iteration number and residual for each iteration. I believe it's L2 norm of the objective function? If yes, is there a way that I could set it to L1 norm for solving the system? Thanks. -- Brian Yang U of Houston -- Brian Yang U of Houston -- 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] How to set L1 norm as the converge test?
Matthew Knepley knep...@gmail.com writes: On Wed, Sep 24, 2014 at 2:31 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Your convergence test routine could call KSPBuildResidual() and then compute the 1-norm of the residual and make any decision it likes based on that norm. See KSPSetConvergenceTest() Note that what Barry is saying is that the convergence theory for CG guarantees monotonicity in the energy (A) norm, but says nothing about L1 so you might get crap. Moreover, if you have an overdetermined linear system, but want to minimize the 1-norm of the residual (instead of the 2-norm that least squares minimizes), you are actually asking to solve a linear program and need to use an LP solver (a very different algorithm). If you merely watch the 1-norm of the residual, it might increase as the least squares solver converges. pgp9EAgFEgbFB.pgp Description: PGP signature
[petsc-users] Setting coordinates
All, I am trying to inset coordinate data into my program line 160 below comes back with a NULL pointer for coordinates. I did not expect this - what do I do, what does it mean? John 140 DM coordDA; 141 Veccoordinates; 142: DMDACoor3d ***coords; 143: PetscScalaru, ux, uy, uxx, uyy; 144: PetscReal D, K, hx, hy, hxdhy, hydhx; 145: PetscInt i,j, k; 150: D = user-D; 151: K = user-K; 152: hx = 1.0/(PetscReal)(info-mx-1); 153: hy = 1.0/(PetscReal)(info-my-1); 154: hxdhy = hx/hy; 155: hydhx = hy/hx; 156: /* 157: Compute function over the locally owned part of the grid 158: */ 159: DMGetCoordinateDA(info-da, coordDA); 160: DMGetCoordinatesLocal(info-da, coordinates); 161: DMDAVecGetArray(coordDA, coordinates, coords);
Re: [petsc-users] Setting coordinates
On Wed, Sep 24, 2014 at 3:03 PM, Alletto, John M john.m.alle...@lmco.com wrotAll, I am trying to inset coordinate data into my program line 160 below comes back with a NULL pointer for coordinates. I did not expect this - what do I do, what does it mean? Coordinates are not created by default since they take up a lot of memory. You either create the storage manually, or use something like this: http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMDASetUniformCoordinates.html Thanks, Matt John 140 DM coordDA; 141 Veccoordinates; 142: DMDACoor3d ***coords; 143: PetscScalaru, ux, uy, uxx, uyy; 144: PetscReal D, K, hx, hy, hxdhy, hydhx; 145: PetscInt i,j, k; 150: D = user-D; 151: K = user-K; 152: hx = 1.0/(PetscReal)(info-mx-1); 153: hy = 1.0/(PetscReal)(info-my-1); 154: hxdhy = hx/hy; 155: hydhx = hy/hx; 156: /* 157: Compute function over the locally owned part of the grid 158: */ 159: DMGetCoordinateDA(info-da, coordDA); 160: DMGetCoordinatesLocal(info-da, coordinates); 161: DMDAVecGetArray(coordDA, coordinates, coords); -- 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] How to set L1 norm as the converge test?
Thanks everyone, I will try KSPSetConvergenceTest to test. On Wed, Sep 24, 2014 at 2:54 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Wed, Sep 24, 2014 at 2:31 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Your convergence test routine could call KSPBuildResidual() and then compute the 1-norm of the residual and make any decision it likes based on that norm. See KSPSetConvergenceTest() Note that what Barry is saying is that the convergence theory for CG guarantees monotonicity in the energy (A) norm, but says nothing about L1 so you might get crap. Moreover, if you have an overdetermined linear system, but want to minimize the 1-norm of the residual (instead of the 2-norm that least squares minimizes), you are actually asking to solve a linear program and need to use an LP solver (a very different algorithm). If you merely watch the 1-norm of the residual, it might increase as the least squares solver converges. -- Brian Yang U of Houston
[petsc-users] I expected the Laplace run with a smaller delta to have a more accurate solution, yet they come out almost exactly the same.
I have set up a test for running a Laplacian solver with 2 sets of data. One has twice as many points as the other. Both cover the same range, I input the X, Y and Z variables set using SetCoordinates. I compare the results with an analytical model. I expected the Laplace run with a smaller delta to have a more accurate solution, yet they come out almost exactly the same. Any ideas?
[petsc-users] Generating xdmf from h5 file.
Hi, i was trying to use petsc_gen_xdmf.py to convert a h5 file to a xdmf file. The h5 file was generated by snes/ex12 which was run as, ex12 -dm_view hdf5:my.h5 When I do, petsc_gen_xdmf.py my.h5 I get the following error, File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 220, in module generateXdmf(sys.argv[1]) File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 208, in generateXdmf time= np.array(h5['time']).flatten() File /usr/lib/python2.7/dist-packages/h5py/_hl/group.py, line 153, in __getitem__ oid = h5o.open(self.id, self._e(name), lapl=self._lapl) File h5o.pyx, line 173, in h5py.h5o.open (h5py/h5o.c:3403) KeyError: unable to open object (Symbol table: Can't open object) I am not sure if the error is on my end. This is on Ubuntu 14.04 with the serial version of hdf5. I built petsc with --download-hdf5, is it necessary to use the same version of hdf5 to generate the xdmf file? Thanks Subramanya
Re: [petsc-users] I expected the Laplace run with a smaller delta to have a more accurate solution, yet they come out almost exactly the same.
On Wed, Sep 24, 2014 at 5:03 PM, Alletto, John M john.m.alle...@lmco.com wrote: I have set up a test for running a Laplacian solver with 2 sets of data. One has twice as many points as the other. Both cover the same range, I input the X, Y and Z variables set using SetCoordinates. I compare the results with an analytical model. I expected the Laplace run with a smaller delta to have a more accurate solution, yet they come out almost exactly the same. Any ideas? I think you may have the wrong idea of accuracy. If you are expecting to converge in the L_2 norm, you must do an integral to get the error, rather than just take the difference of vertex values (that is the l2 norm). You could be seeing superconvergence at the vertices, but I do not know what discretization you are using. 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
Re: [petsc-users] Generating xdmf from h5 file.
On Wed, Sep 24, 2014 at 5:08 PM, subramanya sadasiva pota...@outlook.com wrote: Hi, i was trying to use petsc_gen_xdmf.py to convert a h5 file to a xdmf file. The h5 file was generated by snes/ex12 which was run as, ex12 -dm_view hdf5:my.h5 When I do, petsc_gen_xdmf.py my.h5 I get the following error, File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 220, in module generateXdmf(sys.argv[1]) File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 208, in generateXdmf time= np.array(h5['time']).flatten() File /usr/lib/python2.7/dist-packages/h5py/_hl/group.py, line 153, in __getitem__ oid = h5o.open(self.id, self._e(name), lapl=self._lapl) File h5o.pyx, line 173, in h5py.h5o.open (h5py/h5o.c:3403) KeyError: unable to open object (Symbol table: Can't open object) I am not sure if the error is on my end. This is on Ubuntu 14.04 with the serial version of hdf5. I built petsc with --download-hdf5, is it necessary to use the same version of hdf5 to generate the xdmf file? That code is alpha, and mainly built for me to experiment with an application here, so it is not user-friendly. In your HDF5 file, there is no 'time' since you are not running a TS. This access to h5['time'] should just be protected, and an empty array should be put in if its not there. Matt Thanks Subramanya -- 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] Generating xdmf from h5 file.
Hi Matt, That did not help. Is there any other way to output the mesh to something that paraview can view? I tried outputting the file to a vtk file using ex12 -dm_view vtk:my.vtk:ascii_vtk which, I saw in another post on the forums, but that did not give me any output. (Sorry for sending this twice. but I noticed my reply did not go into the users forum) Subramanya Date: Wed, 24 Sep 2014 17:19:51 -0500 Subject: Re: [petsc-users] Generating xdmf from h5 file. From: knep...@gmail.com To: pota...@outlook.com CC: petsc-users@mcs.anl.gov On Wed, Sep 24, 2014 at 5:08 PM, subramanya sadasiva pota...@outlook.com wrote: Hi, i was trying to use petsc_gen_xdmf.py to convert a h5 file to a xdmf file. The h5 file was generated by snes/ex12 which was run as, ex12 -dm_view hdf5:my.h5 When I do, petsc_gen_xdmf.py my.h5 I get the following error, File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 220, in module generateXdmf(sys.argv[1]) File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 208, in generateXdmf time= np.array(h5['time']).flatten() File /usr/lib/python2.7/dist-packages/h5py/_hl/group.py, line 153, in __getitem__ oid = h5o.open(self.id, self._e(name), lapl=self._lapl) File h5o.pyx, line 173, in h5py.h5o.open (h5py/h5o.c:3403) KeyError: unable to open object (Symbol table: Can't open object) I am not sure if the error is on my end. This is on Ubuntu 14.04 with the serial version of hdf5. I built petsc with --download-hdf5, is it necessary to use the same version of hdf5 to generate the xdmf file? That code is alpha, and mainly built for me to experiment with an application here, so it is not user-friendly. In yourHDF5 file, there is no 'time' since you are not running a TS. This access to h5['time'] should just be protected, andan empty array should be put in if its not there. Matt Thanks Subramanya -- 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] Generating xdmf from h5 file.
On Wed, Sep 24, 2014 at 5:29 PM, subramanya sadasiva pota...@outlook.com wrote: Hi Matt, That did not help. That's not enough description to fix anything, and fixing it will require programming. Is there any other way to output the mesh to something that paraview can view? I tried outputting the file to a vtk file using ex12 -dm_view vtk:my.vtk:ascii_vtk which, I saw in another post on the forums, but that did not give me any output. This is mixing two different things. PETSc has a diagnostic ASCII vtk output, so the type would be ascii, not vtk, and format ascii_vtk . It also has a production VTU output, which is type vtk with format vtk_vtu. Thanks, Matt Subramanya -- Date: Wed, 24 Sep 2014 17:19:51 -0500 Subject: Re: [petsc-users] Generating xdmf from h5 file. From: knep...@gmail.com To: pota...@outlook.com CC: petsc-users@mcs.anl.gov On Wed, Sep 24, 2014 at 5:08 PM, subramanya sadasiva pota...@outlook.com wrote: Hi, i was trying to use petsc_gen_xdmf.py to convert a h5 file to a xdmf file. The h5 file was generated by snes/ex12 which was run as, ex12 -dm_view hdf5:my.h5 When I do, petsc_gen_xdmf.py my.h5 I get the following error, File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 220, in module generateXdmf(sys.argv[1]) File /home/ssadasiv/software/petsc/bin/pythonscripts/petsc_gen_xdmf.py, line 208, in generateXdmf time= np.array(h5['time']).flatten() File /usr/lib/python2.7/dist-packages/h5py/_hl/group.py, line 153, in __getitem__ oid = h5o.open(self.id, self._e(name), lapl=self._lapl) File h5o.pyx, line 173, in h5py.h5o.open (h5py/h5o.c:3403) KeyError: unable to open object (Symbol table: Can't open object) I am not sure if the error is on my end. This is on Ubuntu 14.04 with the serial version of hdf5. I built petsc with --download-hdf5, is it necessary to use the same version of hdf5 to generate the xdmf file? That code is alpha, and mainly built for me to experiment with an application here, so it is not user-friendly. In your HDF5 file, there is no 'time' since you are not running a TS. This access to h5['time'] should just be protected, and an empty array should be put in if its not there. Matt Thanks Subramanya -- 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] DMPlex with spring elements
On Wed, Sep 24, 2014 at 5:38 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks for your response. I'm attaching a pdf with a description of the model. The description of the PetscSection is necessary for the DMNetwork? It looks like DMNetwork does not use a PetscSection. It does internally. It takes the information you give it and makes one, much like ex12 takes the FEM information and makes a PetscSection. Thanks, Matt Miguel On Wed, Sep 24, 2014 at 1:43 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Thanks for your response. My discretization is based on spring elements. For the linear one dimensional case in which each spring has a coefficient k, their jacobian would be this two by two matrix. [ k-k ] [ -k k ] and the internal force [ k ( Ui - Uj) ] [ k ( Uj - Ui) ] where Ui and Uj are the node displacements (just one displacement per node because it's one dimensional) For the two dimensional case, assuming small deformations, we have a four-by-four matrix. Each node has two degrees of freedom. We obtain it by performing the outer product of the vector (t , -t) where t is the vector that connects both nodes in a spring. This is for the case of small deformations. I would need to assemble each spring contribution to the jacobian and the residual like they were finite elements. The springs share nodes, that's how they are connected. This example is just the linear case, I will have to implement a nonlinear case in a similar fashion. Seeing the DMNetwork example, I think it's what I need, although I don't know much of power electric grids and it's hard for me to understand what's going on. Do you have a good reference to be able to follow the code? Please see the attached document which has more description of DMNetwork and the equations for the power grid example. I don't have anything that describes how the power grid example is implemented. For example, why are they adding components to the edges? 475: DMNetworkAddComponent http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/docs/manualpages/DM/D MNetworkAddComponent.html#DMNetworkAddComponent(networkdm,i,componentkey[ 0],pfdata.branch[i-eStart]);Miguel Each edge or node can have several components (limited to 10) attached to it. The term components, taken from the circuit terminology, refers to the elements of a network. For example, a component could be a resistor, inductor, spring, or even edge/vertex weights (for graph problems). For code implementation, component is a data structure that holds the data needed for the residual, Jacobian, or any other function evaluation. In the case of power grid, there are 4 components: branches or transmission lines connecting nodes, buses or nodes, generators that are incident at a subset of the nodes, and loads that are also incident at a subset of the nodes. Each of the these components are defined by their data structures given in pf.h. DMNetwork is a wrapper class of DMPlex specifically for network applications that can be solely described using nodes, edges, and their associated components. If you have a PDE, or need FEM, or need other advanced features then DMPlex would be suitable. Please send us a write-up of your equations so that we can assist you better. Shri On Tue, Sep 23, 2014 at 11:13 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You may also want to take a look at the DMNetwork framework that can be used for general unstructured networks that don't use PDEs. Its description is given in the manual and an example is in src/snes/examples/tutorials/network/pflow. Shri From: Matthew Knepley knep...@gmail.com Date: Tue, 23 Sep 2014 22:40:52 -0400 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Tue, Sep 23, 2014 at 4:01 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi all I was wondering if it could be possible to build a model similar to the example snes/ex12.c, but with spring elements (for elasticity) instead of simplicial elements. Spring elements in a grid, therefore each element would have two nodes and each node two components. There would be more differences, because instead of calling the functions f0,f1,g0,g1,g2 and g3 to build the residual and the jacobian, I would call a routine that would build the residual vector and the jacobian matrix directly. I would not have shape functions whatsoever. My problem is discrete, I don't have a PDE and my equations are algebraic. What is the best way in petsc to solve this problem? Is there any example that I can follow? Thanks in advance Yes, ex12 is fairly specific to FEM. However, I think the right tools for what you want are DMPlex and PetscSection. Here is how I would proceed: 1) Make a DMPlex that
Re: [petsc-users] How to set L1 norm as the converge test?
For LSQR, http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/impls/lsqr/lsqr.c.html#KSPLSQR I checked the lsqr.c implementation, within routine KSPSolve_LSQR, I saw few places using VecNorm and before KSP_Monitor to output the rnorm, there were many intermediate updates for the iteration. I am just wondering whether I could replace NORM2 to NORM1 and remove related square (or square root) operations to achieve my goal? Thanks. On Wed, Sep 24, 2014 at 3:48 PM, Brian Yang brianyang1...@gmail.com wrote: Thanks everyone, I will try KSPSetConvergenceTest to test. On Wed, Sep 24, 2014 at 2:54 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Wed, Sep 24, 2014 at 2:31 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Your convergence test routine could call KSPBuildResidual() and then compute the 1-norm of the residual and make any decision it likes based on that norm. See KSPSetConvergenceTest() Note that what Barry is saying is that the convergence theory for CG guarantees monotonicity in the energy (A) norm, but says nothing about L1 so you might get crap. Moreover, if you have an overdetermined linear system, but want to minimize the 1-norm of the residual (instead of the 2-norm that least squares minimizes), you are actually asking to solve a linear program and need to use an LP solver (a very different algorithm). If you merely watch the 1-norm of the residual, it might increase as the least squares solver converges. -- Brian Yang U of Houston -- Brian Yang U of Houston
Re: [petsc-users] How to set L1 norm as the converge test?
On Wed, Sep 24, 2014 at 6:08 PM, Brian Yang brianyang1...@gmail.com wrote: For LSQR, http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/impls/lsqr/lsqr.c.html#KSPLSQR I checked the lsqr.c implementation, within routine KSPSolve_LSQR, I saw few places using VecNorm and before KSP_Monitor to output the rnorm, there were many intermediate updates for the iteration. I am just wondering whether I could replace NORM2 to NORM1 and remove related square (or square root) operations to achieve my goal? I think the depends on your goal. If it is the print the L1 norm out, then the monitor or convergence test is enough. If it is to have a mathematically consistent way of using the L1 norm, then no you will have to use a different solver. Thanks, Matt Thanks. On Wed, Sep 24, 2014 at 3:48 PM, Brian Yang brianyang1...@gmail.com wrote: Thanks everyone, I will try KSPSetConvergenceTest to test. On Wed, Sep 24, 2014 at 2:54 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Wed, Sep 24, 2014 at 2:31 PM, Barry Smith bsm...@mcs.anl.gov wrote: Brian, Your convergence test routine could call KSPBuildResidual() and then compute the 1-norm of the residual and make any decision it likes based on that norm. See KSPSetConvergenceTest() Note that what Barry is saying is that the convergence theory for CG guarantees monotonicity in the energy (A) norm, but says nothing about L1 so you might get crap. Moreover, if you have an overdetermined linear system, but want to minimize the 1-norm of the residual (instead of the 2-norm that least squares minimizes), you are actually asking to solve a linear program and need to use an LP solver (a very different algorithm). If you merely watch the 1-norm of the residual, it might increase as the least squares solver converges. -- Brian Yang U of Houston -- Brian Yang U of Houston -- 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] DMPlex with spring elements
If you have equations only at the nodes, with a part of it contributed by the edges (springs), then you can use DMNetwork. If you are planning to have equations for the beads in the future, or other higher layers, then DMPlex has better functionality to manage that. Shri From: Miguel Angel Salazar de Troya salazardetr...@gmail.commailto:salazardetr...@gmail.com Date: Wed, 24 Sep 2014 17:38:11 -0500 To: Shri abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements Thanks for your response. I'm attaching a pdf with a description of the model. The description of the PetscSection is necessary for the DMNetwork? It looks like DMNetwork does not use a PetscSection. Miguel On Wed, Sep 24, 2014 at 1:43 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote: Thanks for your response. My discretization is based on spring elements. For the linear one dimensional case in which each spring has a coefficient k, their jacobian would be this two by two matrix. [ k-k ] [ -k k ] and the internal force [ k ( Ui - Uj) ] [ k ( Uj - Ui) ] where Ui and Uj are the node displacements (just one displacement per node because it's one dimensional) For the two dimensional case, assuming small deformations, we have a four-by-four matrix. Each node has two degrees of freedom. We obtain it by performing the outer product of the vector (t , -t) where t is the vector that connects both nodes in a spring. This is for the case of small deformations. I would need to assemble each spring contribution to the jacobian and the residual like they were finite elements. The springs share nodes, that's how they are connected. This example is just the linear case, I will have to implement a nonlinear case in a similar fashion. Seeing the DMNetwork example, I think it's what I need, although I don't know much of power electric grids and it's hard for me to understand what's going on. Do you have a good reference to be able to follow the code? Please see the attached document which has more description of DMNetwork and the equations for the power grid example. I don't have anything that describes how the power grid example is implemented. For example, why are they adding components to the edges? 475: DMNetworkAddComponent http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/docs/manualpages/DM/D MNetworkAddComponent.html#DMNetworkAddComponent(networkdm,i,componentkey[ 0],pfdata.branch[i-eStart]);Miguel Each edge or node can have several components (limited to 10) attached to it. The term components, taken from the circuit terminology, refers to the elements of a network. For example, a component could be a resistor, inductor, spring, or even edge/vertex weights (for graph problems). For code implementation, component is a data structure that holds the data needed for the residual, Jacobian, or any other function evaluation. In the case of power grid, there are 4 components: branches or transmission lines connecting nodes, buses or nodes, generators that are incident at a subset of the nodes, and loads that are also incident at a subset of the nodes. Each of the these components are defined by their data structures given in pf.h. DMNetwork is a wrapper class of DMPlex specifically for network applications that can be solely described using nodes, edges, and their associated components. If you have a PDE, or need FEM, or need other advanced features then DMPlex would be suitable. Please send us a write-up of your equations so that we can assist you better. Shri On Tue, Sep 23, 2014 at 11:13 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote: You may also want to take a look at the DMNetwork framework that can be used for general unstructured networks that don't use PDEs. Its description is given in the manual and an example is in src/snes/examples/tutorials/network/pflow. Shri From: Matthew Knepley knep...@gmail.commailto:knep...@gmail.com Date: Tue, 23 Sep 2014 22:40:52 -0400 To: Miguel Angel Salazar de Troya salazardetr...@gmail.commailto:salazardetr...@gmail.com Cc: petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Tue, Sep 23, 2014 at 4:01 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.commailto:salazardetr...@gmail.com wrote: Hi all I was wondering if it could be possible to build a model similar to the example snes/ex12.c, but with spring elements (for elasticity) instead of simplicial elements. Spring elements in a grid, therefore each element would have two nodes and each node two components. There would be more differences, because instead of calling the functions f0,f1,g0,g1,g2 and g3 to build the residual and the jacobian, I would call a routine that would build
