[petsc-users] Parallelism of the Mat.convert() function
Hello,
The following code returns a different answer depending on how many
processors I use. With one processor, the last MPIAIJ matrix is correctly
formed:
row 0: (0, 1.) (1, 2.) (2, 3.) (3, 4.) (4, -1.)
row 1: (0, 5.) (1, 6.) (2, 7.) (3, 8.) (4, -2.)
row 2: (0, 9.) (1, 10.) (2, 11.) (3, 12.) (4, -3.)
row 3: (0, 13.) (1, 14.) (2, 15.) (3, 16.) (4, -4.)
row 4: (0, 17.) (1, 18.) (2, 19.) (3, 20.) (4, -5.)
row 5: (0, 21.) (1, 22.) (2, 23.) (3, 24.) (4, -6.)
row 6: (0, 25.) (1, 26.) (2, 27.) (3, 28.) (4, -7.)
With two processors though, the column matrix is placed in between:
row 0: (0, 1.) (1, 2.) (2, -1.) (3, 3.) (4, 4.)
row 1: (0, 5.) (1, 6.) (2, -2.) (3, 7.) (4, 8.)
row 2: (0, 9.) (1, 10.) (2, -3.) (3, 11.) (4, 12.)
row 3: (0, 13.) (1, 14.) (2, -4.) (3, 15.) (4, 16.)
row 4: (0, 17.) (1, 18.) (2, -5.) (3, 19.) (4, 20.)
row 5: (0, 21.) (1, 22.) (2, -6.) (3, 23.) (4, 24.)
row 6: (0, 25.) (1, 26.) (2, -7.) (3, 27.) (4, 28.)
Am I not building the nested matrix correctly, perhaps? I am using the
Firedrake PETSc fork. Can you reproduce it?
Thanks,
Miguel
```python
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import COMM_WORLD
input_array = np.array(
[
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
[13, 14, 15, 16],
[17, 18, 19, 20],
[21, 22, 23, 24],
[25, 26, 27, 28],
],
dtype=np.float64,
)
n_11_global_rows, n_11_global_cols = input_array.shape
size = ((None, n_11_global_rows), (None, n_11_global_cols))
mat = PETSc.Mat().createAIJ(size=size, comm=COMM_WORLD)
mat.setUp()
mat.setValues(range(n_11_global_rows), range(n_11_global_cols), input_array)
mat.assemblyBegin()
mat.assemblyEnd()
mat.view()
input_array = np.array([[-1, -2, -3, -4, -5, -6, -7]], dtype=np.float64)
global_rows, global_cols = input_array.T.shape
size = ((None, global_rows), (None, global_cols))
mat_2 = PETSc.Mat().createAIJ(size=size, comm=COMM_WORLD)
mat_2.setUp()
mat_2.setValues(range(global_rows), range(global_cols), input_array)
mat_2.assemblyBegin()
mat_2.assemblyEnd()
N = PETSc.Mat().createNest([[mat, mat_2]], comm=COMM_WORLD)
N.assemblyBegin()
N.assemblyEnd()
PETSc.Sys.Print(f"N sizes: {N.getSize()}")
N.convert("mpiaij").view()
```
--
* Miguel Angel Salazar de Troya *
Head of Software Engineering
EPFL Innovation Park Building C
1015 Lausanne
Email: miguel.sala...@corintis.com
Website:
https://urldefense.us/v3/__http://www.corintis.com__;!!G_uCfscf7eWS!bC7zvBxQx0RuDXxzlOgxr_PdSp5N9ZdzjgTPmjG_ZU5WbNvHboZHFBhZksYgyDF2nO1IRXABTx5zmJLaL2NK_EYg2deym84H$
Re: [petsc-users] Error handling in petsc4py
Hello,
Is there any way to get the PETSc error codes in the python interface? The
test I provided below is just a simple example that I know will run out of
memory.
Miguel
On Wed, Nov 15, 2023 at 10:00 AM Miguel Angel Salazar de Troya <
miguel.sala...@corintis.com> wrote:
> Hello,
>
> The following simple petsc4py snippet runs out of memory, but I would
> like to handle it from python with the usual try-except. Is there any way
> to do so? How can I get the PETSc error codes in the python interface?
>
> Thanks
>
> from petsc4py import PETSc
> import sys, petsc4py
> petsc4py.init(sys.argv)
> try:
> m, n = 100, 100
> A = PETSc.Mat().createAIJ([m, n], nnz=1e6)
>
> A.assemblyBegin()
> A.assemblyEnd()
> except Exception as e:
> print(f"An error occurred: {e}")
>
> An error occurred: error code 55
> [0] MatSeqAIJSetPreallocation() at
> /Users/miguel/repos/firedrake-glacierware/src/petsc/src/mat/impls/aij/seq/aij.c:3942
> [0] MatSeqAIJSetPreallocation_SeqAIJ() at
> /Users/miguel/repos/firedrake-glacierware/src/petsc/src/mat/impls/aij/seq/aij.c:4008
> [0] PetscMallocA() at
> /Users/miguel/repos/firedrake-glacierware/src/petsc/src/sys/memory/mal.c:408
> [0] PetscMallocAlign() at
> /Users/miguel/repos/firedrake-glacierware/src/petsc/src/sys/memory/mal.c:53
> [0] Out of memory. Allocated: 0, Used by process: 59752448
> [0] Memory requested 18446744064984991744
>
>
>
[petsc-users] Error handling in petsc4py
Hello,
The following simple petsc4py snippet runs out of memory, but I would like
to handle it from python with the usual try-except. Is there any way to do
so? How can I get the PETSc error codes in the python interface?
Thanks
from petsc4py import PETSc
import sys, petsc4py
petsc4py.init(sys.argv)
try:
m, n = 100, 100
A = PETSc.Mat().createAIJ([m, n], nnz=1e6)
A.assemblyBegin()
A.assemblyEnd()
except Exception as e:
print(f"An error occurred: {e}")
An error occurred: error code 55
[0] MatSeqAIJSetPreallocation() at
/Users/miguel/repos/firedrake-glacierware/src/petsc/src/mat/impls/aij/seq/aij.c:3942
[0] MatSeqAIJSetPreallocation_SeqAIJ() at
/Users/miguel/repos/firedrake-glacierware/src/petsc/src/mat/impls/aij/seq/aij.c:4008
[0] PetscMallocA() at
/Users/miguel/repos/firedrake-glacierware/src/petsc/src/sys/memory/mal.c:408
[0] PetscMallocAlign() at
/Users/miguel/repos/firedrake-glacierware/src/petsc/src/sys/memory/mal.c:53
[0] Out of memory. Allocated: 0, Used by process: 59752448
[0] Memory requested 18446744064984991744
[petsc-users] Structured (DMDA) vs Unstructured (DMPlex) meshes
Hello, I am trying to understand if I should make the effort to make my code use structured meshes instead of unstructured ones. My domain is cartesian so that is the first check for structured meshes. However, the problem size I am looking at is ~20 million degrees of freedom. My understanding is that for this problem size, most of the time is spent on the solver. In this case, do structured meshes still have an advantage? Can they run Krylov methods faster than when using structured meshes? What about other solvers and preconditioners? Thanks, Miguel
[petsc-users] Segregated solvers in PETSc
Hello, I am solving the Navier-Stokes equation and an advection-diffusion equation to model the temperature. They are fully coupled because the viscosity is temperature dependent. I plan to solve the fully-coupled problem with a segregated approach: I first solve the Navier-Stokes equation for a fixed temperature and feed the velocity to the thermal equation, then use that temperature back in the Navier-Stokes equation to solve for the velocity again until I reach convergence. If I assemble the residual and jacobian for the fully coupled system with the proper fields section for the fieldsplit preconditioner (I am using Firedrake), is there a way to tell PETSc to solve the problem with a segregated approach? Thanks, Miguel
Re: [petsc-users] TSAdjoint not working correctly when using TSThetaSetEndpoint(true)
Thanks! On Tue, Jun 2, 2020, 20:51 Zhang, Hong via petsc-users < petsc-users@mcs.anl.gov> wrote: > I checked the results again and can confirm that it is a bug on our side. > A merge request has been created to fix it: > https://gitlab.com/petsc/petsc/-/merge_requests/2830 > > Thanks, > Hong (Mr.) > > On Jun 2, 2020, at 6:26 PM, Salazar De Troya, Miguel < > salazardet...@llnl.gov> wrote: > > Hong, > > That is not the correct result, however, we can obtain the correct result > if we comment ts.setThetaEndpoint(True). Why is it so? > > Thanks > Miguel > > *From: *"Zhang, Hong" > *Date: *Tuesday, June 2, 2020 at 3:31 PM > *To: *"Salazar De Troya, Miguel" > *Cc: *"petsc-users@mcs.anl.gov" > *Subject: *Re: [petsc-users] TSAdjoint not working correctly when using > TSThetaSetEndpoint(true) > > Miguel, > > After I uncommented ts.setThetaEndpoint(True) and > removed ts.setProblemType(False.ProblemType.LINEAR) in your code, I got the > following result: > > hongzhang@Hongs-MacBook-Pro$ python3 petsc_question_2.py > Cost function > Vec Object: 1 MPI processes > type: seq > 127.781 > Exact value: 127.78112 > > Numerical sensitivity w.r.t. a > 12.778122049945212 > Real sensitivity w.r.t a > 12.7781121978613 > > Numerical sensitivity w.r.t. b > 83.8907580310942 > Real sensitivity w.r.t b > 211.67168296791954 > > -ts_type cn gives the same result. What was the problem you encountered? > > Thanks, > Hong (Mr.) > > > On Jun 2, 2020, at 2:36 PM, Salazar De Troya, Miguel via petsc-users < > petsc-users@mcs.anl.gov> wrote: > > Hello, > > I am attaching a simple example that uses TSAdjoint to calculate the > sensitivity of a 1 degree of freedom ODE. When using theta methods, it > returns the right sensitivity as given by the analytical solution. When I > set `ts.setThetaEndpoint(True)`, it does not work. Is there a theoretical > reason why? I was using this option because using the options > > ts.setType(ts.Type.THETA) > ts.setTheta(0.5) > ts.setThetaEndpoint(True) > > is equivalent to using ts.Type.CN > > Thanks > Miguel > > Miguel A. Salazar de Troya > Postdoctoral Researcher, Lawrence Livermore National Laboratory > B141 > Rm: 1085-5 > Ph: 1(925) 422-6411 > > > >
[petsc-users] Passing a context to TSPreStep()
Hi all I would like to record the information about the failed steps in an adaptive time stepping algorithm. I have noticed that I could call TSPreStep() and get the current time step, but I want to store it in an external data structure. Now, for TSMonitor(), one of the arguments is a context that I can use to store some results. Is there a similar thing for TSPreStep() How could I implement it? Thanks 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] Discrete adjoint and adaptive time stepping
If we calculate the gradient using the discrete adjoint without differentiating the controller, and then calculate the same gradient using finite difference (allowing the time steps to freely change), how different these results are? Miguel On Wed, Apr 29, 2015 at 2:12 PM, Barry Smith bsm...@mcs.anl.gov wrote: On Apr 29, 2015, at 1:51 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: If so, how is it done? We just keep the history of the time-step sizes and then use those time-step sizes when doing the backward integration. Seems simple to me, am I missing something? Barry, if you're using this for optimization, you might want the gradient to be exactly consistent with the objective functional. But for that, you would need to differentiate the controller, which is non-smooth in practice because the number of time steps can change and stages could be rejected (solver failure). Ahh, yes for multiple forward runs yup. One approach would be to save the timestep sequence and have the controller use that in subsequent *forward* runs. If the dynamical system behaves similarly for those steps, it would be okay to use the same timestep sequence. Presumably if that single set of dt (from the first run) is not sufficient for some later runs one could possibly use the union of the dt of several runs for all the runs. (that is run adaptively and inconsistently several runs to determine where dt needs to be controlled and then use the various smaller of the dt at the different time regions for a full set of consistent runs). Of course if the various smaller of the dt requires a tiny dt for all time steps then you are not getting an advantage of adaptive time-stepping, but ok. The idea of actually propagating the gradients through the time-step controller seems IMHO to be absurd; I won't even put it on our game plan until we have many more things done and much more practical experience. Barry -- *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
[petsc-users] Discrete adjoint and adaptive time stepping
Hi I was wondering if your implementation to calculate sensitivities with discrete adjoints will account for the adaptive time stepping. The time step history depends on the parameters, so it should be included in the derivation of the discrete adjoint. Is this calculated in your implementation? If so, how is it done? 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
[petsc-users] FSAL property in TS
Hi I realized that even though there's a boolean named FSAL (First Same As Last) inside the TS object, it seems this property is not implemented yet. I just wanted to double check this is like this and also ask if it will be implemented soon. Thanks 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] FSAL property in TS
Sorry I should have been more specific. I was talking about the explicit Runge-Kutta method. Thanks Miguel On Mon, Apr 13, 2015 at 9:29 AM, Jed Brown j...@jedbrown.org wrote: Miguel Angel Salazar de Troya salazardetr...@gmail.com writes: Hi I realized that even though there's a boolean named FSAL (First Same As Last) inside the TS object, It's not in TS, but is in some implementations. it seems this property is not implemented yet. What do you mean implemented? Some of the methods test their schemes for this property. -- *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
[petsc-users] Partitioning in DMNetwork
Hi In DMNetwork we have edges and vertices and we can add an arbitrary number of variables to both edges and vertices. When the partitioning is carried out, is this done according to the number of edges and vertices or the variables? Say for example, that we assigned a very large number of variables to the first edge, a number much greater than the total number of edges or vertices. After the partitioning, the processor that contains that edge with many variables would have a very large portion of the global vector, wouldn't it? This case is hypothetical and something to avoid, but, could it happen? Would the partitioning be made in some other way to avoid this load unbalance? Thanks 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] DMNetworkGetEdgeRange() in parallel
I modified the DMNetwork example to include the new DM with the modified
section. It has the same problems. Please find attached the code to this
email.
Thanks
On Tue, Feb 24, 2015 at 6:49 PM, Matthew Knepley knep...@gmail.com wrote:
On Tue, Feb 24, 2015 at 6:42 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I implemented the code as agreed, but I don't get the results I expected.
When I create the vector with DMCreateGlobalVector(), I obtain a vector
with a layout similar to the original DMNetwork, instead of the cloned
network with the new PetscSection. The code is as follows:
DMClone(dm, dmEdge);
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dmEdge, eStart, eEnd)
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
DMSetDefaultSection(dmEdge s);
DMCreateGlobalVector(dmEdge, globalVec);
When I get into DMCreateGlobalVector(dmEdge, globalVec) in the debugger,
in the function DMCreateSubDM_Section_Private() I call
PetscSectionView() on the section
I have no idea why you would be in DMCreateSubDM().
Just view globalVec. If the code is as above, it will give you a vector
with that layout. If not
it should be trivial to make a small code and send it. I do this
everywhere is PETSc, so the
basic mechanism certainly works.
Thanks,
Matt
obtained by DMGetDefaultGlobalSection(dm, sectionGlobal), and I obtain
a PetscSection nothing like the one I see when I call PetscSectionView()
on the PetscSection I created above. Does this have anything to do? I
tried to compare this strange PetscSection with the one from the original
DMNetwork, I call DMGetDefaultGlobalSection(dm, sectionGlobal) before
the first line of the snippet above and I get this error message.
0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: DM must have a default PetscSection in order to create a
global PetscSection
Thanks in advance
Miguel
On Mon, Feb 23, 2015 at 3:24 PM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 2:15 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks a lot, the partition should be done before setting up the
section, right?
The partition will be automatic. All you have to do is make the local
section. The DM is already partitioned,
and the Section will inherit that.
Matt
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the
code slower? I'm using the global vector in a TS, using one of the
explicit
RK schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the
global
vector had the edge variables, it would be a much larger vector, and all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because it
might require less modifications in my code. However, I don't know how I
could do this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the
section already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange().
For each edge, I extract or modify its corresponding value in a global
petsc vector. Therefore that vector must have as many components as
edges
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Now it works, thanks a lot!
Miguel
On Wed, Feb 25, 2015 at 10:12 AM, Matthew Knepley knep...@gmail.com wrote:
On Wed, Feb 25, 2015 at 9:59 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Miguel,
I'm a bit tied up today. I'll try to debug this issue tomorrow and get
back to you.
The problem is the way that DMNetwork is using the default section. When
DMCreateGlobalVec() is called,
it uses the default section for its included network-plex, but
DMSetDefaultSection() sets the section associated
with network. This is inconsistent.
I am sending the code which I hacked to show the correct result by using
the private header.
I think that
1) The underlying Plex should be exposed by DMNetworkGetPlex()
2) The default section business should be made consistent
Thanks,
Matt
Thanks,
Shri
From: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Date: Wed, 25 Feb 2015 08:48:10 -0600
To: Matthew Knepley knep...@gmail.com
Cc: Shri abhy...@mcs.anl.gov, petsc-users@mcs.anl.gov
petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
I modified the DMNetwork example to include the new DM with the
modified section. It has the same problems. Please find attached the code
to this email.
Thanks
On Tue, Feb 24, 2015 at 6:49 PM, Matthew Knepley knep...@gmail.com
wrote:
On Tue, Feb 24, 2015 at 6:42 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I implemented the code as agreed, but I don't get the results I
expected. When I create the vector with DMCreateGlobalVector(), I obtain a
vector with a layout similar to the original DMNetwork, instead of the
cloned network with the new PetscSection. The code is as follows:
DMClone(dm, dmEdge);
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dmEdge, eStart, eEnd)
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
DMSetDefaultSection(dmEdge s);
DMCreateGlobalVector(dmEdge, globalVec);
When I get into DMCreateGlobalVector(dmEdge, globalVec) in the
debugger, in the function DMCreateSubDM_Section_Private() I call
PetscSectionView() on the section
I have no idea why you would be in DMCreateSubDM().
Just view globalVec. If the code is as above, it will give you a
vector with that layout. If not
it should be trivial to make a small code and send it. I do this
everywhere is PETSc, so the
basic mechanism certainly works.
Thanks,
Matt
obtained by DMGetDefaultGlobalSection(dm, sectionGlobal), and I
obtain a PetscSection nothing like the one I see when I call
PetscSectionView()
on the PetscSection I created above. Does this have anything to do? I
tried to compare this strange PetscSection with the one from the original
DMNetwork, I call DMGetDefaultGlobalSection(dm, sectionGlobal) before
the first line of the snippet above and I get this error message.
0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: DM must have a default PetscSection in order to create
a global PetscSection
Thanks in advance
Miguel
On Mon, Feb 23, 2015 at 3:24 PM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 2:15 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks a lot, the partition should be done before setting up the
section, right?
The partition will be automatic. All you have to do is make the
local section. The DM is already partitioned,
and the Section will inherit that.
Matt
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the
code slower? I'm using the global vector in a TS, using one of the
explicit
RK schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the
global
vector had the edge variables, it would be a much larger vector, and
all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because
it might require less modifications in my code. However, I don't know
how I
could do this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
I implemented the code as agreed, but I don't get the results I expected.
When I create the vector with DMCreateGlobalVector(), I obtain a vector
with a layout similar to the original DMNetwork, instead of the cloned
network with the new PetscSection. The code is as follows:
DMClone(dm, dmEdge);
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dmEdge, eStart, eEnd)
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
DMSetDefaultSection(dmEdge s);
DMCreateGlobalVector(dmEdge, globalVec);
When I get into DMCreateGlobalVector(dmEdge, globalVec) in the debugger,
in the function DMCreateSubDM_Section_Private() I call PetscSectionView()
on the section obtained by DMGetDefaultGlobalSection(dm, sectionGlobal),
and I obtain a PetscSection nothing like the one I see when I call
PetscSectionView()
on the PetscSection I created above. Does this have anything to do? I tried
to compare this strange PetscSection with the one from the original
DMNetwork, I call DMGetDefaultGlobalSection(dm, sectionGlobal) before the
first line of the snippet above and I get this error message.
0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: DM must have a default PetscSection in order to create a
global PetscSection
Thanks in advance
Miguel
On Mon, Feb 23, 2015 at 3:24 PM, Matthew Knepley knep...@gmail.com wrote:
On Mon, Feb 23, 2015 at 2:15 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks a lot, the partition should be done before setting up the section,
right?
The partition will be automatic. All you have to do is make the local
section. The DM is already partitioned,
and the Section will inherit that.
Matt
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the
code slower? I'm using the global vector in a TS, using one of the explicit
RK schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the global
vector had the edge variables, it would be a much larger vector, and all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because it
might require less modifications in my code. However, I don't know how I
could do this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the
section already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange().
For each edge, I extract or modify its corresponding value in a global
petsc vector. Therefore that vector must have as many components as edges
there are in the network. To extract the value in the vector, I use
VecGetArray() and a variable counter that is incremented in each
iteration.
The array that I obtain in VecGetArray() has to be the same size
than the edge range. That variable counter starts as 0, so if the array
that I obtained in VecGetArray() is x_array, x_array[0] must be the
component in the global vector that corresponds with the start edge given
in DMNetworkGetEdgeRange()
I need that global petsc vector because I will use it in other
operations, it's not just data. Sorry for the confusion. Thanks in
advance.
This sounds like an assembly operation. The usual paradigm is to
compute in the local space
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Thanks, that will help me. Now what I would like to have is the following: if I have two processors and ten edges, the partitioning results in the first processor having the edges 0-4 and the second processor, the edges 5-9. I also have a global vector with as many components as edges, 10. How can I partition it so the first processor also has the 0-4 components and the second, the 5-9 components of the vector? Miguel On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks 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 -- *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] DMNetworkGetEdgeRange() in parallel
I'm iterating through local edges given in DMNetworkGetEdgeRange(). For each edge, I extract or modify its corresponding value in a global petsc vector. Therefore that vector must have as many components as edges there are in the network. To extract the value in the vector, I use VecGetArray() and a variable counter that is incremented in each iteration. The array that I obtain in VecGetArray() has to be the same size than the edge range. That variable counter starts as 0, so if the array that I obtained in VecGetArray() is x_array, x_array[0] must be the component in the global vector that corresponds with the start edge given in DMNetworkGetEdgeRange() I need that global petsc vector because I will use it in other operations, it's not just data. Sorry for the confusion. Thanks in advance. Miguel On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com wrote: On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks, that will help me. Now what I would like to have is the following: if I have two processors and ten edges, the partitioning results in the first processor having the edges 0-4 and the second processor, the edges 5-9. I also have a global vector with as many components as edges, 10. How can I partition it so the first processor also has the 0-4 components and the second, the 5-9 components of the vector? I think it would help to know what you want to accomplish. This is how you are proposing to do it.' If you just want to put data on edges, DMNetwork has a facility for that already. Thanks, Matt Miguel On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks 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 -- *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 -- 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
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Thanks a lot, the partition should be done before setting up the section,
right?
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley knep...@gmail.com wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the code
slower? I'm using the global vector in a TS, using one of the explicit RK
schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the global
vector had the edge variables, it would be a much larger vector, and all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because it might
require less modifications in my code. However, I don't know how I could do
this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the section
already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange().
For each edge, I extract or modify its corresponding value in a global
petsc vector. Therefore that vector must have as many components as edges
there are in the network. To extract the value in the vector, I use
VecGetArray() and a variable counter that is incremented in each iteration.
The array that I obtain in VecGetArray() has to be the same size than
the edge range. That variable counter starts as 0, so if the array that I
obtained in VecGetArray() is x_array, x_array[0] must be the component
in the global vector that corresponds with the start edge given in
DMNetworkGetEdgeRange()
I need that global petsc vector because I will use it in other
operations, it's not just data. Sorry for the confusion. Thanks in advance.
This sounds like an assembly operation. The usual paradigm is to compute
in the local space, and then communicate to get to the global space. So you
would make a PetscSection that had 1 (or some) unknowns on each cell (edge)
and then you can use DMCreateGlobal/LocalVector() and DMLocalToGlobal() to
do this.
Does that make sense?
Thanks,
Matt
Miguel
On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks, that will help me. Now what I would like to have is the
following: if I have two processors and ten edges, the partitioning
results
in the first processor having the edges 0-4 and the second processor, the
edges 5-9. I also have a global vector with as many components as edges,
10. How can I partition it so the first processor also has the 0-4
components and the second, the 5-9 components of the vector?
I think it would help to know what you want to accomplish. This is how
you are proposing to do it.'
If you just want to put data on edges, DMNetwork has a facility for
that already.
Thanks,
Matt
Miguel
On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Miguel,
One possible way is to store the global numbering of any
edge/vertex in the component attached to it. Once the mesh gets
partitioned, the components are also distributed so you can easily
retrieve
the global number of any edge/vertex by accessing its component. This is
what is done in the DMNetwork example pf.c although the global
numbering is
not used for anything.
Shri
From: Matthew Knepley knep...@gmail.com
Date: Mon, 23 Feb 2015 07:54:34 -0600
To: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks. Once I obtain
[petsc-users] DMNetworkGetEdgeRange() in parallel
Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? Thanks 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] DMNetworkGetEdgeRange() in parallel
Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks 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 -- *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
[petsc-users] DMNetworkSetSizes called by all processors
Hi I was checking this example: http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/network/pflow/pf.c.html In line 443, the data is read by processor 0 and DMNetworkSetSizes is called in 461 by all processors, since it is collective. However, the data read by processor 0 has not been broadcasted and numEdges and numVertices and have still a value of zero for the other processors. Does this mean that the other processors create a dummy DMNetwork and wait to receive their partition once the partitioning is called in line 508? Is this a standard procedure to create meshes in parallel? Thanks 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] Changing TSAdapt
Thanks a lot for the fix. Miguel On Thu, Dec 11, 2014 at 10:58 PM, Barry Smith bsm...@mcs.anl.gov wrote: Yikes, seven days out and not yet put into maint! Meanwhile Miguel has to waste a whole day because no one told him about the fix. I won't appologize for fixing it again. Since a fix that is not accessible is not a fix. Barry On Dec 11, 2014, at 10:50 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: Miguel, Thanks for reporting this, you have found a bug in our code. When we changed the adapt type we did not zero out the function pointers for the old basic adaptor hence they were improperly called when the object was finally destroyed at the end. Lisandro fixed this here. https://bitbucket.org/petsc/petsc/commits/40813bd21acd4c08b8080bc6cc1eef9949a22ac8 -- *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
[petsc-users] Changing TSAdapt
Hi I'm trying to use the same TS twice, the first time using the basic TSAdaptType, then I change it to none like this: TSAdapt adapt; TSGetAdapt(ts,adapt); TSAdaptSetType(adapt,none); However, when I destroy the TS, I get this error: 0x7605c4f2 in VecDestroy (v=0x28) at /home/miguel/petsc/src/vec/vec/interface/vector.c:423 423 if (!*v) PetscFunctionReturn(0); (gdb) backtrace #0 0x7605c4f2 in VecDestroy (v=0x28) at /home/miguel/petsc/src/vec/vec/interface/vector.c:423 #1 0x76f330a5 in TSAdaptDestroy_Basic (adapt=0xfdacd0) at /home/miguel/petsc/src/ts/adapt/impls/basic/adaptbasic.c:66 #2 0x76f2c433 in TSAdaptDestroy (adapt=0xfccbc8) at /home/miguel/petsc/src/ts/adapt/interface/tsadapt.c:238 #3 0x76f03093 in TSDestroy (ts=0x7fffdd80) at /home/miguel/petsc/src/ts/interface/ts.c:1906 It's trying to destroy the TSAdaptDestroy_Basic, but I think it was already destroyed when I changed the TSAdaptType to none, is this true? How can I effectively change the TSAdaptType without having this error? Thanks 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
[petsc-users] PetscFunctionBegin, -malloc_dump and C++ classes with PETSc objects
Hi I'm implementing a problem using the TS. Most of my functions are methods inside of a class, except for the callbacks (to form the RHS and the TS monitor), which are outside of the class, although in the same .C file where the class methods are implemented. For these callbacks I followed the network example: https://bitbucket.org/petsc/petsc/src/a614f7369d93d476173b8fc6bf2463276dcbdb3a/src/snes/examples/tutorials/network/pflow/pf.c?at=master Therefore, the callbacks have the PetscFunctionBegin at the beginning and PetscFunctionReturn(0) at the end. My problems come when I run the program with -malloc_dump and I get a lot of unfreed memory. Inspecting the output I see that the line of my code where the memory is allocated corresponds with the line when PetscFunctionBegin is called. Later in the file, I see that the function DMGetLocalVector() is called within a petsc internal routine (at the file dmget.c). I also call this routine in my callback methods few lines after PetscFunctionBegin. The procedure that I follow to use the local vectors is as the one in the network example. For vectors that I want to modify this is: ierr = DMGetLocalVector(networkdm,localX);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); ierr = VecGetArray(localX,xarr);CHKERRQ(ierr); Modify values in xarr ierr = VecRestoreArray(localX,xarr);CHKERRQ(ierr); ierr = DMLocalToGlobalBegin(networkdm,localX,INSERT_VALUES,X);CHKERRQ(ierr); ierr = DMLocalToGlobalEnd(networkdm,localX,INSERT_VALUES,X);CHKERRQ(ierr); ierr = DMRestoreLocalVector(networkdm,localX);CHKERRQ(ierr); One last thing that I think it might be a issue here is how I destroy the petsc objects. I create the petsc objects within a class. For instance, the class has a petsc vector that later passes to the TS object to get the solution. To destroy the petsc objects, I use the class destructor, where at the end I call PetscFinalize() Inside the class I pass the callbacks to the TS routines that need them (e.g. TSSetRHSFunction() ) I can compile the code and run it, but many memory allocations are not freed. What can be the issue here? Do you know of an example using C++ classes to implement PETSc methods? 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] PetscFunctionBegin, -malloc_dump and C++ classes with PETSc objects
PetscFinalize() is inside of the class destructor (the last line of the destructor), so when the object goes out of scope, the class destructor is called and PetscFinalize() as well. Is it better to have PetscFinalize() outside of the destructor and call the destructor explicitly before? Miguel On Tue, Nov 18, 2014 at 11:32 AM, Barry Smith bsm...@mcs.anl.gov wrote: On Nov 18, 2014, at 11:19 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I'm implementing a problem using the TS. Most of my functions are methods inside of a class, except for the callbacks (to form the RHS and the TS monitor), which are outside of the class, although in the same .C file where the class methods are implemented. For these callbacks I followed the network example: https://bitbucket.org/petsc/petsc/src/a614f7369d93d476173b8fc6bf2463276dcbdb3a/src/snes/examples/tutorials/network/pflow/pf.c?at=master Therefore, the callbacks have the PetscFunctionBegin at the beginning and PetscFunctionReturn(0) at the end. My problems come when I run the program with -malloc_dump and I get a lot of unfreed memory. Inspecting the output I see that the line of my code where the memory is allocated corresponds with the line when PetscFunctionBegin is called. This is normal. We cannot register the exact line the memory allocated, only the location of the PETScFunctionBegin; Later in the file, I see that the function DMGetLocalVector() is called within a petsc internal routine (at the file dmget.c). I also call this routine in my callback methods few lines after PetscFunctionBegin. The procedure that I follow to use the local vectors is as the one in the network example. For vectors that I want to modify this is: ierr = DMGetLocalVector(networkdm,localX);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); ierr = VecGetArray(localX,xarr);CHKERRQ(ierr); Modify values in xarr ierr = VecRestoreArray(localX,xarr);CHKERRQ(ierr); ierr = DMLocalToGlobalBegin(networkdm,localX,INSERT_VALUES,X);CHKERRQ(ierr); ierr = DMLocalToGlobalEnd(networkdm,localX,INSERT_VALUES,X);CHKERRQ(ierr); ierr = DMRestoreLocalVector(networkdm,localX);CHKERRQ(ierr); One last thing that I think it might be a issue here is how I destroy the petsc objects. I create the petsc objects within a class. For instance, the class has a petsc vector that later passes to the TS object to get the solution. To destroy the petsc objects, I use the class destructor, where at the end I call PetscFinalize() Inside the class I pass the callbacks to the TS routines that need them (e.g. TSSetRHSFunction() ) I can compile the code and run it, but many memory allocations are not freed. What can be the issue here? Do you know of an example using C++ classes to implement PETSc methods? Thanks in advance. Do you call the class destructor yourself explicitly before the PetscFinalize()? You need to, otherwise the class may not be destroyed until after PetscFinalize() and hence the PETSc objects won't be freed when you call PetscFinalize(). You also need to make sure that you destroy ALL PETSc objects, if you miss even one PETSc object, since the objects have references to each other it may be that many PETSc objects do not get freed and hence -malloc_dump shows many objects still alive. Barry 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 -- *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] Adaptive controllers in TS
Sorry to be pushy, but could anyone help me on this? Thanks Miguel On Thu, Oct 23, 2014 at 10:46 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I decided to implement the continuous adjoint because it is clearer to me how to do it and the gradient will converge anyways. I was trying to interpolate the solution, but in the adjoint analysis I need to do it once the simulation has finished. I have saved all the solutions for each time step. My idea was to reset the TS at the time step immediately previous to the time step in which I want to interpolate to: interpolate_timestep : where I want to interpolate to previous_timestep : last time step in our time step history which is smaller that the interpolate_timestep Current_Sol : Solution at the time step previous_timestep TSSetTime(ts,previous_timestep); TSSetSolution(ts,Current_Sol); TSSetRetainStages(ts,PETSC_TRUE); TSInterpolate(ts,interpolate_timestep,X); Vector X should have a close value to Current_Sol, but it's nowhere close. I also compare it to the solution of the next time step after previous_timestep (therefore interpolate_timestep is between these guys) and it's not close either. I've read that TSInterpolate() has to be extended to support continuous adjoints. Thanks Miguel On Tue, Oct 21, 2014 at 6:04 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That might be a reasonable argument, but I'm not sure. This is one of the papers that explains it. I'm re-reading to see if I skipped any details: http://www.sciencedirect.com/science/article/pii/S0377042709006062 Miguel On Mon, Oct 20, 2014 at 11:42 PM, Jed Brown j...@jedbrown.org wrote: Miguel Angel Salazar de Troya salazardetr...@gmail.com writes: Thanks for your response I'm struggling with this problem because the literature is not clear for me on how to calculate the discrete adjoint with adaptive time stepping algorithms. They cover the details when automatic differentiation tools are used. They mention that because the time step depend on the solution, it also depends on the parameter. Hence, there are terms that represent the derivative of the time step w.r.t. the parameters. What it is confusing is that they mention these terms must be removed. I don't understand this. I'm planning to hard-code the discrete adjoint problem (and use the TS if possible). Are they suggesting that the time step sizes for a given run should be frozen (at least locally) so that you have consistent gradients for a while? -- *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 -- *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 -- *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] Adaptive controllers in TS
I decided to implement the continuous adjoint because it is clearer to me how to do it and the gradient will converge anyways. I was trying to interpolate the solution, but in the adjoint analysis I need to do it once the simulation has finished. I have saved all the solutions for each time step. My idea was to reset the TS at the time step immediately previous to the time step in which I want to interpolate to: interpolate_timestep : where I want to interpolate to previous_timestep : last time step in our time step history which is smaller that the interpolate_timestep Current_Sol : Solution at the time step previous_timestep TSSetTime(ts,previous_timestep); TSSetSolution(ts,Current_Sol); TSSetRetainStages(ts,PETSC_TRUE); TSInterpolate(ts,interpolate_timestep,X); Vector X should have a close value to Current_Sol, but it's nowhere close. I also compare it to the solution of the next time step after previous_timestep (therefore interpolate_timestep is between these guys) and it's not close either. I've read that TSInterpolate() has to be extended to support continuous adjoints. Thanks Miguel On Tue, Oct 21, 2014 at 6:04 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: That might be a reasonable argument, but I'm not sure. This is one of the papers that explains it. I'm re-reading to see if I skipped any details: http://www.sciencedirect.com/science/article/pii/S0377042709006062 Miguel On Mon, Oct 20, 2014 at 11:42 PM, Jed Brown j...@jedbrown.org wrote: Miguel Angel Salazar de Troya salazardetr...@gmail.com writes: Thanks for your response I'm struggling with this problem because the literature is not clear for me on how to calculate the discrete adjoint with adaptive time stepping algorithms. They cover the details when automatic differentiation tools are used. They mention that because the time step depend on the solution, it also depends on the parameter. Hence, there are terms that represent the derivative of the time step w.r.t. the parameters. What it is confusing is that they mention these terms must be removed. I don't understand this. I'm planning to hard-code the discrete adjoint problem (and use the TS if possible). Are they suggesting that the time step sizes for a given run should be frozen (at least locally) so that you have consistent gradients for a while? -- *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 -- *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] Cloning a DMNetwork
Thanks. Miguel On Wed, Oct 22, 2014 at 10:16 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, I've added DMClone_Network to shri/clone-dmnetwork and next branch. Not that DMClone() creates a shallow copy of DMNetwork. If you change the values of the components in the cloned DM then those will also be reflected in the original DM (since they share the address space). Shri -Original Message- From: Shri abhy...@mcs.anl.gov Date: Wed, 22 Oct 2014 03:44:35 + To: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] Cloning a DMNetwork There is no clone method implemented for DMNetwork yet. I'll work on it tomorrow. Thanks, Shri From: Miguel Angel Salazar de Troya salazardetr...@gmail.com Date: Tue, 21 Oct 2014 18:10:42 -0500 To: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: [petsc-users] Cloning a DMNetwork Hi all I have a DMNetwork and I'm interested in cloning it to have another network with the same components. Can I do this with DMClone()? I have noticed that the components are not copied. Is there any function that takes care of this or it has to be hard coded? Thanks 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 -- *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] Adaptive controllers in TS
That might be a reasonable argument, but I'm not sure. This is one of the papers that explains it. I'm re-reading to see if I skipped any details: http://www.sciencedirect.com/science/article/pii/S0377042709006062 Miguel On Mon, Oct 20, 2014 at 11:42 PM, Jed Brown j...@jedbrown.org wrote: Miguel Angel Salazar de Troya salazardetr...@gmail.com writes: Thanks for your response I'm struggling with this problem because the literature is not clear for me on how to calculate the discrete adjoint with adaptive time stepping algorithms. They cover the details when automatic differentiation tools are used. They mention that because the time step depend on the solution, it also depends on the parameter. Hence, there are terms that represent the derivative of the time step w.r.t. the parameters. What it is confusing is that they mention these terms must be removed. I don't understand this. I'm planning to hard-code the discrete adjoint problem (and use the TS if possible). Are they suggesting that the time step sizes for a given run should be frozen (at least locally) so that you have consistent gradients for a while? -- *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
[petsc-users] Cloning a DMNetwork
Hi all I have a DMNetwork and I'm interested in cloning it to have another network with the same components. Can I do this with DMClone()? I have noticed that the components are not copied. Is there any function that takes care of this or it has to be hard coded? Thanks 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
[petsc-users] Adaptive controllers in TS
Hi all I'm trying to find out what is the adaptive controller scheme used in the Runge-Kutta time stepping algorithm. Basically I want to know the function that determines the next time step. I see that the next time step is set with the function TSAdaptChoose(), which calls the pointer function (*choose)() within the structure _TSAdaptOps, but where is this set? I need to know this for an adjoint analysis to calculate gradients. I was wondering if there is an example of such analysis with the TS structure. I've seen that the RK source file includes a funciton to interpolate the solution, something that is necessary in an adjoint analysis, hence my question about the example. Thanks 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] Adaptive controllers in TS
Thanks for your response I'm struggling with this problem because the literature is not clear for me on how to calculate the discrete adjoint with adaptive time stepping algorithms. They cover the details when automatic differentiation tools are used. They mention that because the time step depend on the solution, it also depends on the parameter. Hence, there are terms that represent the derivative of the time step w.r.t. the parameters. What it is confusing is that they mention these terms must be removed. I don't understand this. I'm planning to hard-code the discrete adjoint problem (and use the TS if possible). I know this might not be the place to ask this. Maybe Hong can help me here. Miguel On Sun, Oct 19, 2014 at 1:41 PM, Jed Brown j...@jedbrown.org wrote: Miguel Angel Salazar de Troya salazardetr...@gmail.com writes: Hi all I'm trying to find out what is the adaptive controller scheme used in the Runge-Kutta time stepping algorithm. Basically I want to know the function that determines the next time step. I see that the next time step is set with the function TSAdaptChoose(), which calls the pointer function (*choose)() within the structure _TSAdaptOps, but where is this set? via TSAdaptSetType. The default will be TSAdaptChoose_Basic (adaptbasic.c). You can find this information easily by running in a debugger. I need to know this for an adjoint analysis to calculate gradients. I was wondering if there is an example of such analysis with the TS structure. I've seen that the RK source file includes a funciton to interpolate the solution, something that is necessary in an adjoint analysis, hence my question about the example. Interpolation arises for continuous adjoints. Be aware that if you want consistent gradients, you'll want a discrete adjoint, which integrates backward in time using the adjoint RK method. The stages of the adjoint RK method coincide with the forward method, so no interpolation is needed. Hong (Cc'd) is developing an example that can integrate discrete adjoint equations and we intend to eventually incorporate that functionality into the library. -- *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
[petsc-users] When is the TSMonitor called?
Hi all I want to save a certain quantity at the very beginning of each time step, before any TSFunction of any kind is called. Can I do this within the TSMonitor? It would also be ok to save that quantity at the very end of the time step. Thanks 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
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
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
[petsc-users] Best way to save entire solution history of a TS
Hello I'm performing an adjoint sensitivity analysis in a transient problem and I need to save the entire solution history. I don't have many degrees of freedom or time steps so I should be able to avoid memory problems. I've thought of creating a std::vectorVec and save all the PETSc Vec that I get from TSGetSolution there. Would this be a problem if I'm running the simulation in parallel? Is there a better way to save the entire solution history in memory? Thanks 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
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
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
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
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
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. 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
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
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
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
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
[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
[petsc-users] Initial estimation on SNES and KSP
Dear all SNES uses internally a KSP to solve the linear system of equations right? Now the case that we had a linear system of equations that we are solving with SNES, how could we set the initial estimation for the KSP? If we just included the option -ksp_initial_guess_nonzero, the KSP will grab the vector X we passed to the SNES? 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] Adaptive mesh refinement in Petsc
Thanks a lot for your responses. I will get started with MOOSE. On Thu, May 1, 2014 at 8:04 PM, Derek Gaston fried...@gmail.com wrote: Miguel, I'm the lead for the MOOSE Framework project Barry spoke of... we would love to help you get up and running with adaptive finite elements for solid mechanics with MOOSE. If you are doing fairly normal solid mechanics using small or large strain formulations with some plasticity... most of what you need is already there. You may need to plug in your particular material model but that's about it. Mesh adaptivity is built-in and should work out of the box. The major benefit of using MOOSE is that you can easily couple in other physics (like heat conduction, chemistry and more) and of course you have full access to all the power of PETSc. I recommend going through the Getting Started material on http://www.mooseframework.org to get set up... and go ahead and create yourself a new Application using these instructions: http://mooseframework.org/create-an-app/ . That Application will already have full access to our solid mechanics capabilities (as well as tons of other stuff like heat conduction, chemistry, etc.). After that - join up on the moose-users mailing list and you can get in touch with everyone else doing solid mechanics with MOOSE who can point you in the right direction depending on your particular application. Let me know if you have any questions... Derek On Thu, May 1, 2014 at 6:31 PM, Barry Smith bsm...@mcs.anl.gov wrote: You also could likely benefit from Moose http://www.mooseframework.orgit sits on top of libMesh which sits on top of PETSc and manages almost all of what you need for finite element analysis. Barry On May 1, 2014, at 7:19 PM, Matthew Knepley knep...@gmail.com wrote: On Thu, May 1, 2014 at 6:14 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hello everybody I want to implement an adaptive mesh refinement library in a code written in petsc. I have checked out some of the available libraries, but I want to work with the latest petsc-dev version and I am sure there will be many incompatibilities. So far I think I'll end up working with one of these libraries: SAMRAI, Chombo, libMesh and deal II. Before I start checking out each of them and learn how to use them I though I would ask you guys which one you would recommend. My code would be a finite element analysis in solid mechanics. I would like to take full advantage of petsc capabilities, but I would not mind start with some restrictions. I hope my question is not too broad. SAMRAI, Chombo, and Deal II are all structured adaptive refinement codes, whereas LibMesh is unstructured. If you want unstructured, there is really no other game in town. If you use deal II, I would suggest trying out p4est underneath which gives great scalability. My understanding is that Chombo is mostly used for finite volume and SAMRAI and deal II for finite element, but this could be out of date. Matt Take care 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 -- *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
[petsc-users] Adaptive mesh refinement in Petsc
Hello everybody I want to implement an adaptive mesh refinement library in a code written in petsc. I have checked out some of the available libraries, but I want to work with the latest petsc-dev version and I am sure there will be many incompatibilities. So far I think I'll end up working with one of these libraries: SAMRAI, Chombo, libMesh and deal II. Before I start checking out each of them and learn how to use them I though I would ask you guys which one you would recommend. My code would be a finite element analysis in solid mechanics. I would like to take full advantage of petsc capabilities, but I would not mind start with some restrictions. I hope my question is not too broad. Take care 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] Elasticity tensor in ex52
I understand. So if we had a linear elastic material, the weak form would
be something like
phi_{i,j} C_{ijkl} u_{k,l}
where the term C_{ijkl} u_{k,l} would correspond to the term f_1(U,
gradient_U) in equation (1) in your paper I mentioned above, and g_3 would
be C_{ijkl}. Therefore the indices {i,j} would be equivalent to the indices
ic, id you mentioned before and jc and jd would be {k,l}?
For g3[ic, id, jc, jd], transforming the four dimensional array to linear
memory would be like this:
g3[((ic*Ncomp+id)*dim+jc)*dim+jd] = 1.0;
where Ncomp and dim are equal to the problem's spatial dimension.
However, in the code, there are only two loops, to exploit the symmetry of
the fourth order identity tensor:
for (compI = 0; compI Ncomp; ++compI) {
for (d = 0; d dim; ++d) {
g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
}
}
Therefore the tensor entries that are set to one are:
g3[0,0,0,0]
g3[1,1,0,0]
g3[0,0,1,1]
g3[1,1,1,1]
This would be equivalent to the fourth order tensor \delta_{ij}
\delta_{kl}, but I think the one we need is \delta_{ik} \delta_{jl},
because it is the derivative of a matrix with respect itself (or the
derivative of a gradient with respect to itself). This is assuming the
indices of g3 correspond to what I said.
Thanks in advance.
Miguel
On Apr 19, 2014 6:19 PM, Matthew Knepley knep...@gmail.com wrote:
On Sat, Apr 19, 2014 at 5:25 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks for your response. Now I understand a bit better your
implementation, but I am still confused. Is the g3 function equivalent to
the f_{1,1} function in the matrix in equation (3) in this paper:
http://arxiv.org/pdf/1309.1204v2.pdf ? If so, why would it have terms
that depend on the trial fields? I think this is what is confusing me.
Yes, it is. It has no terms that depend on the trial fields. It is just
indexed by the number of components in that field. It is
a continuum thing, which has nothing to do with the discretization. That
is why it acts pointwise.
Matt
On Sat, Apr 19, 2014 at 11:35 AM, Matthew Knepley knep...@gmail.comwrote:
On Fri, Apr 18, 2014 at 1:23 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Hello everybody.
First, I am taking this example from the petsc-dev version, I am not
sure if I should have posted this in another mail-list, if so, my
apologies.
In this example, for the elasticity case, function g3 is built as:
void g3_elas(const PetscScalar u[], const PetscScalar gradU[], const
PetscScalar a[], const PetscScalar gradA[], const PetscReal x[],
PetscScalar g3[])
{
const PetscInt dim = spatialDim;
const PetscInt Ncomp = spatialDim;
PetscInt compI, d;
for (compI = 0; compI Ncomp; ++compI) {
for (d = 0; d dim; ++d) {
g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
}
}
}
Therefore, a fourth-order tensor is represented as a vector. I was
checking the indices for different iterator values, and they do not seem to
match the vectorization that I have in mind. For a two dimensional case,
the indices for which the value is set as 1 are:
compI = 0 , d = 0 - index = 0
compI = 0 , d = 1 - index = 3
compI = 1 , d = 0 - index = 12
compI = 1 , d = 1 - index = 15
The values for the first and last seem correct to me, but they other
two are confusing me. I see that this elasticity tensor (which is the
derivative of the gradient by itself in this case) would be a four by four
identity matrix in its matrix representation, so the indices in between
would be 5 and 10 instead of 3 and 12, if we put one column on top of each
other.
I have read this a few times, but I cannot understand that you are
asking. The simplest thing I can
respond is that we are indexing a row-major array, using the indices:
g3[ic, id, jc, jd]
where ic indexes the components of the trial field, id indexes the
derivative components,
jc indexes the basis field components, and jd its derivative components.
Matt
I guess my question is then, how did you vectorize the fourth order
tensor?
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
--
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] Elasticity tensor in ex52
Thanks for your response. Now I understand a bit better your
implementation, but I am still confused. Is the g3 function equivalent to
the f_{1,1} function in the matrix in equation (3) in this paper:
http://arxiv.org/pdf/1309.1204v2.pdf ? If so, why would it have terms that
depend on the trial fields? I think this is what is confusing me.
On Sat, Apr 19, 2014 at 11:35 AM, Matthew Knepley knep...@gmail.com wrote:
On Fri, Apr 18, 2014 at 1:23 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Hello everybody.
First, I am taking this example from the petsc-dev version, I am not sure
if I should have posted this in another mail-list, if so, my apologies.
In this example, for the elasticity case, function g3 is built as:
void g3_elas(const PetscScalar u[], const PetscScalar gradU[], const
PetscScalar a[], const PetscScalar gradA[], const PetscReal x[],
PetscScalar g3[])
{
const PetscInt dim = spatialDim;
const PetscInt Ncomp = spatialDim;
PetscInt compI, d;
for (compI = 0; compI Ncomp; ++compI) {
for (d = 0; d dim; ++d) {
g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
}
}
}
Therefore, a fourth-order tensor is represented as a vector. I was
checking the indices for different iterator values, and they do not seem to
match the vectorization that I have in mind. For a two dimensional case,
the indices for which the value is set as 1 are:
compI = 0 , d = 0 - index = 0
compI = 0 , d = 1 - index = 3
compI = 1 , d = 0 - index = 12
compI = 1 , d = 1 - index = 15
The values for the first and last seem correct to me, but they other two
are confusing me. I see that this elasticity tensor (which is the
derivative of the gradient by itself in this case) would be a four by four
identity matrix in its matrix representation, so the indices in between
would be 5 and 10 instead of 3 and 12, if we put one column on top of each
other.
I have read this a few times, but I cannot understand that you are asking.
The simplest thing I can
respond is that we are indexing a row-major array, using the indices:
g3[ic, id, jc, jd]
where ic indexes the components of the trial field, id indexes the
derivative components,
jc indexes the basis field components, and jd its derivative components.
Matt
I guess my question is then, how did you vectorize the fourth order
tensor?
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
--
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
