[petsc-users] DMPlex with spring elements

[Abhyankar, Shrirang G.]

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 <salazardetroya at gmail.com>
Date:  Thu, 25 Sep 2014 12:57:53 -0500
To:  Shri <abhyshr at mcs.anl.gov>
Cc:  <petsc-users at 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." <abhyshr at mcs.anl.gov>
>wrote:
>
>
>
>From: Miguel Angel Salazar de Troya <salazardetroya at gmail.com>
>Date: Thu, 25 Sep 2014 11:43:13 -0500
>To: Shri <abhyshr at mcs.anl.gov>
>Cc: <petsc-users at 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." <abhyshr at 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 <salazardetroya at gmail.com>
>Date: Thu, 25 Sep 2014 10:52:16 -0500
>To: Shri <abhyshr at mcs.anl.gov>
>Cc: "petsc-users at mcs.anl.gov" <petsc-users at 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.
><abhyshr at 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.
>><abhyshr at 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 <salazardetroya at gmail.com>
>>Date: Wed, 24 Sep 2014 17:38:11 -0500
>>To: Shri <abhyshr at mcs.anl.gov>
>>Cc: "petsc-users at mcs.anl.gov" <petsc-users at 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.
>><abhyshr at 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,componentke
>>>y
>>>[
>>>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.
>>><abhyshr at 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 <knepley at gmail.com>
>>>Date:  Tue, 23 Sep 2014 22:40:52 -0400
>>>To:  Miguel Angel Salazar de Troya <salazardetroya at gmail.com>
>>>Cc:  "petsc-users at mcs.anl.gov" <petsc-users at 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
>>>><salazardetroya at 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>
>>>>salaza11 at 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
>>>salaza11 at illinois.edu
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>--
>>Miguel Angel Salazar de Troya
>>Graduate Research Assistant
>>Department of Mechanical Science and Engineering
>>University of Illinois at Urbana-Champaign
>>(217) 550-2360
>>salaza11 at illinois.edu
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>--
>>Miguel Angel Salazar de Troya
>>Graduate Research Assistant
>>Department of Mechanical Science and Engineering
>>University of Illinois at Urbana-Champaign
>>(217) 550-2360
>>salaza11 at illinois.edu
>
>
>
>
>
>
>
>
>
>-- 
>Miguel Angel Salazar de Troya
>Graduate Research Assistant
>Department of Mechanical Science and Engineering
>University of Illinois at Urbana-Champaign
>(217) 550-2360
>salaza11 at illinois.edu