Re: [petsc-users] DMPlex with spring elements
On Mon, Sep 29, 2014 at 4:55 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi all I'm bumping this post because I have more questions related to the same problem. I am looping over the edges of my DMNetwork, then I obtain the vertices that make up each edge with DMNetworkGetConnectedNode(). Each of these vertices have two variables (or actually, two degrees of freedom for my problem). My intentions are to modify the solution vector entries that are affected by these variables in each vertex. I would call the function DMNetworkGetVariableOffset() to do this. What happens if one of the vertices is a ghost vertex? Can I still modify the solution vector? My problem is that the edge has information to provide to these nodes. This is unfortunately not explained well enough in our documentation. I will try to put a more substantial explanation in the manual. PETSc operates in two different spaces: local space and global space. In the global space, every degree of freedom (dof, or unknown) is owned by a single process and is a single entry which lives on that process in a global Vec. This is what solvers, like GMRES, deal with and also what is plotted for visualization. In the local space, each process has all of the dofs it owns AND all of those that it shares with other owning processes (ghosts). Thus in the local space a given dof can be entries on several processes in a local Vec. You usually use the local space to compute residuals/Jacobians. Our paradigm is: compute residual in the local space, so that you can set values for ghosts, and then call DMLocalToGlobal() to put those values into a global Vec (using some combination operator, like +) which can be handed to the solver. Thanks, Matt Thanks Miguel On Fri, Sep 26, 2014 at 12:33 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tut orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian
Re: [petsc-users] DMPlex with spring elements
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com Date: Mon, 29 Sep 2014 16:55:14 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements Hi all I'm bumping this post because I have more questions related to the same problem. I am looping over the edges of my DMNetwork, then I obtain the vertices that make up each edge with DMNetworkGetConnectedNode(). Each of these vertices have two variables (or actually, two degrees of freedom for my problem). My intentions are to modify the solution vector entries that are affected by these variables in each vertex. I would call the function DMNetworkGetVariableOffset() to do this. What happens if one of the vertices is a ghost vertex? Can I still modify the solution vector? My problem is that the edge has information to provide to these nodes. Sorry for the delay. I think you would want to use http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMNetworkIsG hostVertex.html and not modify the value for the ghost vertex. Shri Thanks Miguel On Fri, Sep 26, 2014 at 12:33 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example (http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tu t orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr
Re: [petsc-users] DMPlex with spring elements
On Tue, Sep 30, 2014 at 11:22 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: From: Miguel Angel Salazar de Troya salazardetr...@gmail.com Date: Mon, 29 Sep 2014 16:55:14 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements Hi all I'm bumping this post because I have more questions related to the same problem. I am looping over the edges of my DMNetwork, then I obtain the vertices that make up each edge with DMNetworkGetConnectedNode(). Each of these vertices have two variables (or actually, two degrees of freedom for my problem). My intentions are to modify the solution vector entries that are affected by these variables in each vertex. I would call the function DMNetworkGetVariableOffset() to do this. What happens if one of the vertices is a ghost vertex? Can I still modify the solution vector? My problem is that the edge has information to provide to these nodes. Sorry for the delay. I think you would want to use http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMNetworkIsG hostVertex.html and not modify the value for the ghost vertex. This depends on the discretization. In FD, you do not modify ghosts, but in FEM you do. Matt Shri Thanks Miguel On Fri, Sep 26, 2014 at 12:33 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tu t orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use
Re: [petsc-users] DMPlex with spring elements
No worries. Thanks for your responses. I'm assuming you suggested to use DMNetworkIsGhostVertex() and not modify its value for the case in which I were using the global vectors and not the local vectors, where it is possible, as Matt suggested. Miguel On Tue, Sep 30, 2014 at 11:22 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: From: Miguel Angel Salazar de Troya salazardetr...@gmail.com Date: Mon, 29 Sep 2014 16:55:14 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements Hi all I'm bumping this post because I have more questions related to the same problem. I am looping over the edges of my DMNetwork, then I obtain the vertices that make up each edge with DMNetworkGetConnectedNode(). Each of these vertices have two variables (or actually, two degrees of freedom for my problem). My intentions are to modify the solution vector entries that are affected by these variables in each vertex. I would call the function DMNetworkGetVariableOffset() to do this. What happens if one of the vertices is a ghost vertex? Can I still modify the solution vector? My problem is that the edge has information to provide to these nodes. Sorry for the delay. I think you would want to use http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMNetworkIsG hostVertex.html and not modify the value for the ghost vertex. Shri Thanks Miguel On Fri, Sep 26, 2014 at 12:33 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tu t orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts
Re: [petsc-users] DMPlex with spring elements
My reply was based on the assumption that the values set at a vertex do not have contributions from other vertices or edges. In which case, you would not want to set any value for the ghost vertex. For example, consider the following network v1e1 v2e2 v3 | | --- | that gets distributed on two processors as follows v1 e1v2' [0] | --- | v2e2 v3 [1] | - | Here, v2' is the ghost vertex on P0. Now, if you want to set f(v2) = 2, then you do not need to set anything for f(v2') and call DMLocalToGlobal() at the end with INSERT_VALUES flag (You can also use ADD_VALUES flag assuming that you have zeroed the local vector initially). On the other hand, if you have f(v2) = f(v1) + f(v3) then you would need to set f(v2') = f(v1) and f(v2) = f(v3) and call DMLocalToGlobal with the flag ADD_VALUES. Shri From: Miguel Angel Salazar de Troya salazardetr...@gmail.com Date: Tue, 30 Sep 2014 11:35:31 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements No worries. Thanks for your responses. I'm assuming you suggested to use DMNetworkIsGhostVertex() and not modify its value for the case in which I were using the global vectors and not the local vectors, where it is possible, as Matt suggested. Miguel On Tue, Sep 30, 2014 at 11:22 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: From: Miguel Angel Salazar de Troya salazardetr...@gmail.com Date: Mon, 29 Sep 2014 16:55:14 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements Hi all I'm bumping this post because I have more questions related to the same problem. I am looping over the edges of my DMNetwork, then I obtain the vertices that make up each edge with DMNetworkGetConnectedNode(). Each of these vertices have two variables (or actually, two degrees of freedom for my problem). My intentions are to modify the solution vector entries that are affected by these variables in each vertex. I would call the function DMNetworkGetVariableOffset() to do this. What happens if one of the vertices is a ghost vertex? Can I still modify the solution vector? My problem is that the edge has information to provide to these nodes. Sorry for the delay. I think you would want to use http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMNetworkIs G hostVertex.html and not modify the value for the ghost vertex. Shri Thanks Miguel On Fri, Sep 26, 2014 at 12:33 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example (http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/t u t orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works
Re: [petsc-users] DMPlex with spring elements
Hi all I'm bumping this post because I have more questions related to the same problem. I am looping over the edges of my DMNetwork, then I obtain the vertices that make up each edge with DMNetworkGetConnectedNode(). Each of these vertices have two variables (or actually, two degrees of freedom for my problem). My intentions are to modify the solution vector entries that are affected by these variables in each vertex. I would call the function DMNetworkGetVariableOffset() to do this. What happens if one of the vertices is a ghost vertex? Can I still modify the solution vector? My problem is that the edge has information to provide to these nodes. Thanks Miguel On Fri, Sep 26, 2014 at 12:33 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tut orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how
Re: [petsc-users] DMPlex with spring elements
Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 -- 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
That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tutorials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 -- 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 Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tutorials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 -- 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 -- 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
DMNetwork does not have the functionality to constraint a dof. The user has to do it by specifying the constraint equation. As Jed suggested, you could simply not set any values in the rows/columns (except the diagonal) corresponding to the constrained dof in the Jacobian evaluation to keep the matrix symmetric. Shri From: Matthew Knepley knep...@gmail.com Date: Thu, 25 Sep 2014 17:17:52 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data;Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses.Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- 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
Re: [petsc-users] DMPlex with spring elements
Yeah, but doesn't it only work with the local vectors localX and localF? Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tutorials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 -- 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 -- 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
Re: [petsc-users] DMPlex with spring elements
On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tutorials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 -- 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 -- What most experimenters take for granted
Re: [petsc-users] DMPlex with spring elements
Ok thanks. Miguel On Fri, Sep 26, 2014 at 10:28 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tutorials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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 -- 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
Re: [petsc-users] DMPlex with spring elements
What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example (http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tut orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data;Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses.Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you
Re: [petsc-users] DMPlex with spring elements
I understand. Thanks a lot. Miguel On Fri, Sep 26, 2014 at 10:53 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: What Matt is saying is that there are two interfaces in PETSc for setting the residual evaluation routine: i) SNESSetFunction takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzroutine(SNES snes, Vec X, Vec F, void* ctx); X and F are the global solution and residual vectors. To compute the global residual evaluation, typically one does -- (a) scattering X and F onto local vectors localX and localF (DMGlobalToLocal), (b) computing the local residual, and (c) gathering the localF in the global F (DMLocalToGlobal). This is what is done in the example. ii) DMSNESSetFunctionLocal takes in a function pointer for the residual evaluation routine that has the prototype PetscErrorCode xyzlocalroutine(DM, Vec localX, localF, void* ctx) In this case, the localX and localF get passed to the routine. So, you only have to do the local residual evaluation. PETSc does the LocalToGlobal gather to form the global residual. I chose to use SNESSetFunction in the example. You can use either of them. Shri From: Matthew Knepley knep...@gmail.com Date: Fri, 26 Sep 2014 10:28:26 -0500 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: Jed Brown j...@jedbrown.org, Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMPlex with spring elements On Fri, Sep 26, 2014 at 10:26 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Yeah, but doesn't it only work with the local vectors localX and localF? I am telling you what the interface for the functions is. You can do whatever you want inside. Matt Miguel On Fri, Sep 26, 2014 at 10:10 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 10:06 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That means that if we call SNESSetFunction() we don't build the residual vector in parallel? In the pflow example ( http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/src/snes/examples/tut orials/network/pflow/pf.c.html) the function FormFunction() (Input for SNESSetFunction() works with the local vectors. I don't understand this. FormFunction() in that link clearly takes in a global vector X and returns a global vector F. Inside, it converts them to local vectors. This is exactly what you would do for a function given to SNESSetFunction(). Matt Thanks Miguel On Fri, Sep 26, 2014 at 9:34 AM, Matthew Knepley knep...@gmail.com wrote: On Fri, Sep 26, 2014 at 9:31 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. I had another question about the DM and SNES and TS. There are similar routines to assign the residual and jacobian evaluation to both objects. For the SNES case are: DMSNESSetFunctionLocal DMSNESSetJacobianLocal What are the differences of these with: SNESSetFunction SNESSetJacobian SNESSetFunction() expects the user to construct the entire parallel residual vector. DMSNESSetFunctionLocal() expects the user to construct the local pieces of the residual, and then it automatically calls DMLocalToGlobal() to assembly the full residual. It also converts the input from global vectors to local vectors, and in the case of DMDA multidimensional arrays. Thanks, Matt and when should we use each? With Local, it is meant to evaluate the function/jacobian for the elements in the local processor? I could get the local edges in DMNetwork by calling DMNetworkGetEdgeRange? Miguel On Thu, Sep 25, 2014 at 5:17 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data;Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses.Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org
Re: [petsc-users] DMPlex with spring elements
Thanks. I think the term Component was confusing me, I thought it was related to the components of a field. I think this would be useful to me if I wanted to assign coordinates to the vertices, wouldn't it? Also, I was wondering how to set up dirichlet boundary conditions, basically fixing certain nodes position. Could I do it as the function *SetInitialValues *does it in the pflow example? These values are used to eliminate the zeroth-order energy modes of the stiffness matrix? Last question, in my case I have two degrees of freedom per node, when I grab the offset with DMNetworkVariableOffset, that's for the first degree of freedom in that node and the second degree of freedom would just be offset+1? Miguel On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: 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.com Date: Wed, 24 Sep 2014 17:38:11 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
Thanks. I think the term Component was confusing me, I thought it was related to the components of a field. I think this would be useful to me if I wanted to assign coordinates to the vertices, wouldn't it? Yes. You can put whatever data you want in the component data structure. Also, I was wondering how to set up dirichlet boundary conditions, basically fixing certain nodes position. You can add a component at each node with a field marking whether the node is a boundary node. Could I do it as the function SetInitialValues does it in the pflow example? No. You need to put in the component data structure before calling DMNetworkAddComponent() These values are used to eliminate the zeroth-order energy modes of the stiffness matrix? Last question, in my case I have two degrees of freedom per node, when I grab the offset with DMNetworkVariableOffset, that's for the first degree of freedom in that node and the second degree of freedom would just be offset+1? Yes. Shri Miguel On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: 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.com Date: Wed, 24 Sep 2014 17:38:11 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
Thanks. Once I have marked the nodes that are fixed nodes using the component data structure, how can I process it later? I mean, at what point does the solver know that those degrees of freedom are actually fixed and how I can tell it that they are fixed? Miguel On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Thanks. I think the term Component was confusing me, I thought it was related to the components of a field. I think this would be useful to me if I wanted to assign coordinates to the vertices, wouldn't it? Yes. You can put whatever data you want in the component data structure. Also, I was wondering how to set up dirichlet boundary conditions, basically fixing certain nodes position. You can add a component at each node with a field marking whether the node is a boundary node. Could I do it as the function SetInitialValues does it in the pflow example? No. You need to put in the component data structure before calling DMNetworkAddComponent() These values are used to eliminate the zeroth-order energy modes of the stiffness matrix? Last question, in my case I have two degrees of freedom per node, when I grab the offset with DMNetworkVariableOffset, that's for the first degree of freedom in that node and the second degree of freedom would just be offset+1? Yes. Shri Miguel On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: 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.com Date: Wed, 24 Sep 2014 17:38:11 -0500 To: Shri abhy...@mcs.anl.gov Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
The solver does not know anything about the boundary conditions. You would have
to specify it to the solver by describing the appropriate equations. For e.g.
in the power grid example, there is a part in the residual evaluation
if (bus-ide == REF_BUS || bus-ide == ISOLATED_BUS) {
farr[offset] = 0.0;
farr[offset+1] = 0.0;
break;
}
This sets the residual at the nodes marked with REF_BUS or ISOLATED_BUS to 0.0.
You can do something similar.
Shri
From: Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 10:52:16 -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. Once I have marked the nodes that are fixed nodes using the component
data structure, how can I process it later? I mean, at what point does the
solver know that those degrees of freedom are actually fixed and how I can tell
it that they are fixed?
Miguel
On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote:
Thanks. I think the term Component was confusing me, I thought it was
related to the components of a field. I think this would be useful to me
if I wanted to assign coordinates to the vertices, wouldn't it?
Yes. You can put whatever data you want in the component data structure.
Also, I was wondering how to set up dirichlet boundary conditions,
basically fixing certain nodes position.
You can add a component at each node with a field marking whether the node
is a boundary node.
Could I do it as the function SetInitialValues does it in the pflow
example?
No. You need to put in the component data structure before calling
DMNetworkAddComponent()
These values are used to eliminate the zeroth-order energy modes of the
stiffness matrix?
Last question, in my case I have two degrees of freedom per node, when I
grab the offset with DMNetworkVariableOffset, that's for the first degree
of freedom in that node and the second degree of freedom would just be
offset+1?
Yes.
Shri
Miguel
On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote:
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
Re: [petsc-users] DMPlex with spring elements
I see, and I guess I would have to assign a value of one to the diagonal
entry of that degree of freedom in the Jacobian right? Wouldn't this break
the symmetry of the Jacobian (in case it were symmetric)?
Thanks
Miguel
On Sep 25, 2014 11:32 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
The solver does not know anything about the boundary conditions. You
would have to specify it to the solver by describing the appropriate
equations. For e.g. in the power grid example, there is a part in the
residual evaluation
if (bus-ide == REF_BUS || bus-ide == ISOLATED_BUS) {
farr[offset] = 0.0;
farr[offset+1] = 0.0;
break;
}
This sets the residual at the nodes marked with REF_BUS or ISOLATED_BUS
to 0.0. You can do something similar.
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 10:52:16 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
Thanks. Once I have marked the nodes that are fixed nodes using the
component data structure, how can I process it later? I mean, at what point
does the solver know that those degrees of freedom are actually fixed and
how I can tell it that they are fixed?
Miguel
On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Thanks. I think the term Component was confusing me, I thought it was
related to the components of a field. I think this would be useful to me
if I wanted to assign coordinates to the vertices, wouldn't it?
Yes. You can put whatever data you want in the component data structure.
Also, I was wondering how to set up dirichlet boundary conditions,
basically fixing certain nodes position.
You can add a component at each node with a field marking whether the node
is a boundary node.
Could I do it as the function SetInitialValues does it in the pflow
example?
No. You need to put in the component data structure before calling
DMNetworkAddComponent()
These values are used to eliminate the zeroth-order energy modes of the
stiffness matrix?
Last question, in my case I have two degrees of freedom per node, when I
grab the offset with DMNetworkVariableOffset, that's for the first degree
of freedom in that node and the second degree of freedom would just be
offset+1?
Yes.
Shri
Miguel
On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
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.com
Date: Wed, 24 Sep 2014 17:38:11 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
From: Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 11:43:13 -0500
To: Shri abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
I see, and I guess I would have to assign a value of one to the diagonal entry
of that degree of freedom in the Jacobian right?
Yes.
Wouldn't this break the symmetry of the Jacobian (in case it were symmetric)?
No.
Shri
Thanks
Miguel
On Sep 25, 2014 11:32 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote:
The solver does not know anything about the boundary conditions. You would have
to specify it to the solver by describing the appropriate equations. For e.g.
in the power grid example, there is a part in the residual evaluation
if (bus-ide == REF_BUS || bus-ide == ISOLATED_BUS) {
farr[offset] = 0.0;
farr[offset+1] = 0.0;
break;
}
This sets the residual at the nodes marked with REF_BUS or ISOLATED_BUS to 0.0.
You can do something similar.
Shri
From: Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 10:52:16 -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. Once I have marked the nodes that are fixed nodes using the component
data structure, how can I process it later? I mean, at what point does the
solver know that those degrees of freedom are actually fixed and how I can tell
it that they are fixed?
Miguel
On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote:
Thanks. I think the term Component was confusing me, I thought it was
related to the components of a field. I think this would be useful to me
if I wanted to assign coordinates to the vertices, wouldn't it?
Yes. You can put whatever data you want in the component data structure.
Also, I was wondering how to set up dirichlet boundary conditions,
basically fixing certain nodes position.
You can add a component at each node with a field marking whether the node
is a boundary node.
Could I do it as the function SetInitialValues does it in the pflow
example?
No. You need to put in the component data structure before calling
DMNetworkAddComponent()
These values are used to eliminate the zeroth-order energy modes of the
stiffness matrix?
Last question, in my case I have two degrees of freedom per node, when I
grab the offset with DMNetworkVariableOffset, that's for the first degree
of freedom in that node and the second degree of freedom would just be
offset+1?
Yes.
Shri
Miguel
On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote:
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
Re: [petsc-users] DMPlex with spring elements
Why not? Wouldn't we have a row of zeros except for the diagonal term? The
column that corresponds to that degree of from doesn't have to be zero,
right?
Thanks
Miguel
On Sep 25, 2014 12:38 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 11:43:13 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
I see, and I guess I would have to assign a value of one to the diagonal
entry of that degree of freedom in the Jacobian right?
Yes.
Wouldn't this break the symmetry of the Jacobian (in case it were
symmetric)?
No.
Shri
Thanks
Miguel
On Sep 25, 2014 11:32 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
The solver does not know anything about the boundary conditions. You
would have to specify it to the solver by describing the appropriate
equations. For e.g. in the power grid example, there is a part in the
residual evaluation
if (bus-ide == REF_BUS || bus-ide == ISOLATED_BUS) {
farr[offset] = 0.0;
farr[offset+1] = 0.0;
break;
}
This sets the residual at the nodes marked with REF_BUS or ISOLATED_BUS
to 0.0. You can do something similar.
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 10:52:16 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
Thanks. Once I have marked the nodes that are fixed nodes using the
component data structure, how can I process it later? I mean, at what point
does the solver know that those degrees of freedom are actually fixed and
how I can tell it that they are fixed?
Miguel
On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Thanks. I think the term Component was confusing me, I thought it was
related to the components of a field. I think this would be useful to me
if I wanted to assign coordinates to the vertices, wouldn't it?
Yes. You can put whatever data you want in the component data structure.
Also, I was wondering how to set up dirichlet boundary conditions,
basically fixing certain nodes position.
You can add a component at each node with a field marking whether the
node
is a boundary node.
Could I do it as the function SetInitialValues does it in the pflow
example?
No. You need to put in the component data structure before calling
DMNetworkAddComponent()
These values are used to eliminate the zeroth-order energy modes of the
stiffness matrix?
Last question, in my case I have two degrees of freedom per node, when I
grab the offset with DMNetworkVariableOffset, that's for the first
degree
of freedom in that node and the second degree of freedom would just be
offset+1?
Yes.
Shri
Miguel
On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
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.com
Date: Wed, 24 Sep 2014 17:38:11 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
You are right. The Jacobian for the power grid application is indeed
non-symmetric. Is that a problem for your application?
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 12:57:53 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
Why not? Wouldn't we have a row of zeros except for the diagonal term?
The column that corresponds to that degree of from doesn't have to be
zero, right?
Thanks
Miguel
On Sep 25, 2014 12:38 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 11:43:13 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
I see, and I guess I would have to assign a value of one to the diagonal
entry of that degree of freedom in the Jacobian right?
Yes.
Wouldn't this break the symmetry of the Jacobian (in case it were
symmetric)?
No.
Shri
Thanks
Miguel
On Sep 25, 2014 11:32 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
The solver does not know anything about the boundary conditions. You
would have to specify it to the solver by describing the appropriate
equations. For e.g. in the power grid example, there is a part in the
residual evaluation
if (bus-ide == REF_BUS || bus-ide == ISOLATED_BUS) {
farr[offset] = 0.0;
farr[offset+1] = 0.0;
break;
}
This sets the residual at the nodes marked with REF_BUS or ISOLATED_BUS
to 0.0. You can do something similar.
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 10:52:16 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
Thanks. Once I have marked the nodes that are fixed nodes using the
component data structure, how can I process it later? I mean, at what
point does the solver know that those degrees of freedom are actually
fixed and how I can tell it that they
are fixed?
Miguel
On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Thanks. I think the term Component was confusing me, I thought it was
related to the components of a field. I think this would be useful to me
if I wanted to assign coordinates to the vertices, wouldn't it?
Yes. You can put whatever data you want in the component data structure.
Also, I was wondering how to set up dirichlet boundary conditions,
basically fixing certain nodes position.
You can add a component at each node with a field marking whether the node
is a boundary node.
Could I do it as the function SetInitialValues does it in the pflow
example?
No. You need to put in the component data structure before calling
DMNetworkAddComponent()
These values are used to eliminate the zeroth-order energy modes of the
stiffness matrix?
Last question, in my case I have two degrees of freedom per node, when I
grab the offset with DMNetworkVariableOffset, that's for the first degree
of freedom in that node and the second degree of freedom would just be
offset+1?
Yes.
Shri
Miguel
On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
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.com
Date: Wed, 24 Sep 2014 17:38:11 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
You are right. The Jacobian for the power grid application is indeed
non-symmetric. Is that a problem for your application?
If you need a symmetric Jacobian, you can use the BC facility in
PetscSection, which eliminates the
variables completely. This is how the FEM examples, like ex12, work.
Matt
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 12:57:53 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
Why not? Wouldn't we have a row of zeros except for the diagonal term?
The column that corresponds to that degree of from doesn't have to be
zero, right?
Thanks
Miguel
On Sep 25, 2014 12:38 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 11:43:13 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
I see, and I guess I would have to assign a value of one to the diagonal
entry of that degree of freedom in the Jacobian right?
Yes.
Wouldn't this break the symmetry of the Jacobian (in case it were
symmetric)?
No.
Shri
Thanks
Miguel
On Sep 25, 2014 11:32 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
The solver does not know anything about the boundary conditions. You
would have to specify it to the solver by describing the appropriate
equations. For e.g. in the power grid example, there is a part in the
residual evaluation
if (bus-ide == REF_BUS || bus-ide == ISOLATED_BUS) {
farr[offset] = 0.0;
farr[offset+1] = 0.0;
break;
}
This sets the residual at the nodes marked with REF_BUS or ISOLATED_BUS
to 0.0. You can do something similar.
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Thu, 25 Sep 2014 10:52:16 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMPlex with spring elements
Thanks. Once I have marked the nodes that are fixed nodes using the
component data structure, how can I process it later? I mean, at what
point does the solver know that those degrees of freedom are actually
fixed and how I can tell it that they
are fixed?
Miguel
On Thu, Sep 25, 2014 at 10:27 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Thanks. I think the term Component was confusing me, I thought it was
related to the components of a field. I think this would be useful to me
if I wanted to assign coordinates to the vertices, wouldn't it?
Yes. You can put whatever data you want in the component data structure.
Also, I was wondering how to set up dirichlet boundary conditions,
basically fixing certain nodes position.
You can add a component at each node with a field marking whether the node
is a boundary node.
Could I do it as the function SetInitialValues does it in the pflow
example?
No. You need to put in the component data structure before calling
DMNetworkAddComponent()
These values are used to eliminate the zeroth-order energy modes of the
stiffness matrix?
Last question, in my case I have two degrees of freedom per node, when I
grab the offset with DMNetworkVariableOffset, that's for the first degree
of freedom in that node and the second degree of freedom would just be
offset+1?
Yes.
Shri
Miguel
On Wed, Sep 24, 2014 at 9:52 PM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
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.com
Date: Wed, 24 Sep 2014 17:38:11 -0500
To: Shri abhy...@mcs.anl.gov
Cc: petsc-users@mcs.anl.gov 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.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
Re: [petsc-users] DMPlex with spring elements
Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). pgpkrj2MXOlP9.pgp Description: PGP signature
Re: [petsc-users] DMPlex with spring elements
If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 Thu, Sep 25, 2014 at 5:15 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. Would that be with PetscSectionSetConstraintDof ? For that I will need the PetscSection, DofSection, within DMNetwork, how can I obtain it? I could cast it to DM_Network from the dm, networkdm, declared in the main program, maybe something like this: DM_Network *network = (DM_Network*) networkdm-data; Then I would loop over the vertices and call PetscSectionSetConstraintDof if it's a boundary node (by checking the corresponding component) I admit to not completely understanding DMNetwork. However, it eventually builds a PetscSection for data layout, which you could get from DMGetDefaultSection(). The right thing to do is find where it builds the Section, and put in your BC there, but that sounds like it would entail coding. Thanks, Matt Thanks for your responses. Miguel On Thu, Sep 25, 2014 at 2:42 PM, Jed Brown j...@jedbrown.org wrote: Matthew Knepley knep...@gmail.com writes: On Thu, Sep 25, 2014 at 1:46 PM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: You are right. The Jacobian for the power grid application is indeed non-symmetric. Is that a problem for your application? If you need a symmetric Jacobian, you can use the BC facility in PetscSection, which eliminates the variables completely. This is how the FEM examples, like ex12, work. You can also use MatZeroRowsColumns() or do the equivalent transformation during assembly (my preference). -- *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 -- 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
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
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
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
[petsc-users] DMPlex with spring elements
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 Miguel -- *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 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 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
Re: [petsc-users] DMPlex with spring elements
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
