[petsc-users] Problems imposing boundary conditions

[Maximilian Hartig]

Hello,

I am working on a transient structural FEM code with PETSc. I managed to create 
a slow but functioning program with the use of petscFE and a TS solver. The 
code runs fine until I try to restrict movement in all three spatial directions 
for one face. I then get the error which is attached below.
So apparently DMPlexMatSetClosure tries to write/read beyond what was priorly 
allocated. I do however not call MatSeqAIJSetPreallocation myself in the code. 
So I’m unsure where to start looking for the bug. In my understanding, PETSc 
should know from the DM how much space to allocate.
Could you kindly give me a hint?

Thanks,

Max

0 SNES Function norm 2.508668036663e-06 
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: New nonzero at (41754,5) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off 
this check
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.7.5-3223-g99077fc  GIT Date: 
2017-02-28 13:41:43 -0600
[0]PETSC ERROR: ./S3 on a arch-linux-gnu-intel named XXX by hartig Fri Mar  
3 17:55:57 2017
[0]PETSC ERROR: Configure options PETSC_ARCH=arch-linux-gnu-intel 
--with-cc=/opt/intel/compilers_and_libraries/linux/mpi/intel64/bin/mpiicc 
--with-cxx=/opt/intel/compilers_and_libraries/linux/mpi/intel64/bin/mpiicpc 
--with-fc=/opt/intel/compilers_and_libraries/linux/mpi/intel64/bin/mpiifort 
--download-ml
[0]PETSC ERROR: #1 MatSetValues_SeqAIJ() line 455 in 
/home/hartig/petsc/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #2 MatSetValues() line 1270 in 
/home/hartig/petsc/src/mat/interface/matrix.c
[0]PETSC ERROR: [0]ERROR in DMPlexMatSetClosure
[0]mat for sieve point 60
[0]mat row indices[0] = 41754
[0]mat row indices[1] = 41755
[0]mat row indices[2] = 41756
[0]mat row indices[3] = 41760
[0]mat row indices[4] = 41761
[0]mat row indices[5] = 41762
[0]mat row indices[6] = 41766
[0]mat row indices[7] = -41768
[0]mat row indices[8] = 41767
[0]mat row indices[9] = 41771
[0]mat row indices[10] = -41773
[0]mat row indices[11] = 41772
[0]mat row indices[12] = 41776
[0]mat row indices[13] = 41777
[0]mat row indices[14] = 41778
[0]mat row indices[15] = 41782
[0]mat row indices[16] = -41784
[0]mat row indices[17] = 41783
[0]mat row indices[18] = 261
[0]mat row indices[19] = -263
[0]mat row indices[20] = 262
[0]mat row indices[21] = 24318
[0]mat row indices[22] = 24319
[0]mat row indices[23] = 24320
[0]mat row indices[24] = -7
[0]mat row indices[25] = -8
[0]mat row indices[26] = 6
[0]mat row indices[27] = 1630
[0]mat row indices[28] = -1632
[0]mat row indices[29] = 1631
[0]mat row indices[30] = 41757
[0]mat row indices[31] = 41758
[0]mat row indices[32] = 41759
[0]mat row indices[33] = 41763
[0]mat row indices[34] = 41764
[0]mat row indices[35] = 41765
[0]mat row indices[36] = 41768
[0]mat row indices[37] = 41769
[0]mat row indices[38] = 41770
[0]mat row indices[39] = 41773
[0]mat row indices[40] = 41774
[0]mat row indices[41] = 41775
[0]mat row indices[42] = 41779
[0]mat row indices[43] = 41780
[0]mat row indices[44] = 41781
[0]mat row indices[45] = 41784
[0]mat row indices[46] = 41785
[0]mat row indices[47] = 41786
[0]mat row indices[48] = 263
[0]mat row indices[49] = 264
[0]mat row indices[50] = 265
[0]mat row indices[51] = 24321
[0]mat row indices[52] = 24322
[0]mat row indices[53] = 24323
[0]mat row indices[54] = 5
[0]mat row indices[55] = 6
[0]mat row indices[56] = 7
[0]mat row indices[57] = 1632
[0]mat row indices[58] = 1633
[0]mat row indices[59] = 1634
[0] 1.29801 0.0998428 -0.275225 1.18171 -0.0093323 0.055045 1.18146 0.00525527 
-0.11009 -0.588378 0.264666 -0.0275225 -2.39586 -0.210511 0.22018 -0.621071 
0.0500786 0.137613 -0.180869 -0.0974804 0.0344031 -0.0302673 -0.09 0. -0.145175 
-0.00383346 -0.00688063 0.300442 -0.00868618 -0.0275225 8.34577e-11 0. 0. 
4.17288e-11 0. 0. 4.17288e-11 0. 0. 4.17288e-11 0. 0. 4.17288e-11 0. 0. 
2.08644e-11 0. 0. -1.04322e-11 0. 0. -1.04322e-11 0. 0. -1.56483e-11 0. 0. 
-1.56483e-11 0. 0.
[0] 0.0998428 0.590663 -0.0320009 0.0270594 0.360043 0.0297282 -0.0965489 
0.270936 -0.0652389 0.32647 -0.171351 -0.00845384 -0.206902 -0.657189 0.0279137 
0.0500786 -0.197561 -0.0160508 -0.12748 -0.138996 0.0408591 -0.06 -0.105935 
0.0192308 0.00161757 -0.0361182 0.0042968 -0.0141372 0.0855084 -0.000284088 0. 
8.34577e-11 0. 0. 4.17288e-11 0. 0. 4.17288e-11 0. 0. 4.17288e-11 0. 0. 
4.17288e-11 0. 0. 2.08644e-11 0. 0. -1.04322e-11 0. 0. -1.04322e-11 0. 0. 
-1.56483e-11 0. 0. -1.56483e-11 0.
[0] -0.275225 -0.0320009 0.527521 -0.055045 0.0285918 0.234796 -0.165135 
-0.0658071 0.322754 0.0275225 -0.0207062 -0.114921 0.33027 0.0418706 -0.678455 
0.137613 -0.0160508 -0.235826 0.0344031 0.0312437 -0.0845583 0. 0.0288462 
-0.0302673 0.00688063 0.00443884 -0.0268103 -0.0412838 -0.000426133 0.0857668 
0. 0. 8.34577e-11 0. 0. 4.17288e-11 

Re: [petsc-users] Problems imposing boundary conditions

[Maximilian Hartig]

Yes Sander, your assessment is correct. I use DMPlex and specify the BC using 
DMLabel. I do however get this error also when running in serial.

Thanks,
Max

> On 3 Mar 2017, at 22:14, Matthew Knepley  wrote:
> 
> On Fri, Mar 3, 2017 at 12:56 PM, Sander Arens  <mailto:sander.ar...@ugent.be>> wrote:
> Max, 
> 
> I'm assuming you use DMPlex for your mesh? If so, did you only specify the 
> faces in the DMLabel (and not vertices or edges). Do you get this error only 
> in parallel? 
> 
> If so, I can confirm this bug. I submitted a pull request for this yesterday.
> 
> Yep, I saw Sander's pull request. I will get in merged in tomorrow when I get 
> home to Houston.
> 
>   Thanks,
> 
>  Matt
>  
> On 3 March 2017 at 18:43, Lukas van de Wiel  <mailto:lukas.drinkt.t...@gmail.com>> wrote:
> You have apparently preallocated the non-zeroes of you matrix, and the room 
> was insufficient to accommodate all your equations. 
>  
> What happened after you tried:
> 
> MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)
> 
> 
> Cheers
> Lukas
> 
> On Fri, Mar 3, 2017 at 6:37 PM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Hello,
> 
> I am working on a transient structural FEM code with PETSc. I managed to 
> create a slow but functioning program with the use of petscFE and a TS 
> solver. The code runs fine until I try to restrict movement in all three 
> spatial directions for one face. I then get the error which is attached below.
> So apparently DMPlexMatSetClosure tries to write/read beyond what was priorly 
> allocated. I do however not call MatSeqAIJSetPreallocation myself in the 
> code. So I’m unsure where to start looking for the bug. In my understanding, 
> PETSc should know from the DM how much space to allocate.
> Could you kindly give me a hint?
> 
> Thanks,
> 
> Max
> 
> 0 SNES Function norm 2.508668036663e-06
> [0]PETSC ERROR: - Error Message 
> --
> [0]PETSC ERROR: Argument out of range
> [0]PETSC ERROR: New nonzero at (41754,5) caused a malloc
> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off 
> this check
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html 
> <http://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting.
> [0]PETSC ERROR: Petsc Development GIT revision: v3.7.5-3223-g99077fc  GIT 
> Date: 2017-02-28 13:41:43 -0600
> [0]PETSC ERROR: ./S3 on a arch-linux-gnu-intel named XXX by hartig Fri 
> Mar  3 17:55:57 2017
> [0]PETSC ERROR: Configure options PETSC_ARCH=arch-linux-gnu-intel 
> --with-cc=/opt/intel/compilers_and_libraries/linux/mpi/intel64/bin/mpiicc 
> --with-cxx=/opt/intel/compilers_and_libraries/linux/mpi/intel64/bin/mpiicpc 
> --with-fc=/opt/intel/compilers_and_libraries/linux/mpi/intel64/bin/mpiifort 
> --download-ml
> [0]PETSC ERROR: #1 MatSetValues_SeqAIJ() line 455 in 
> /home/hartig/petsc/src/mat/impls/aij/seq/aij.c
> [0]PETSC ERROR: #2 MatSetValues() line 1270 in 
> /home/hartig/petsc/src/mat/interface/matrix.c
> [0]PETSC ERROR: [0]ERROR in DMPlexMatSetClosure
> [0]mat for sieve point 60
> [0]mat row indices[0] = 41754
> [0]mat row indices[1] = 41755
> [0]mat row indices[2] = 41756
> [0]mat row indices[3] = 41760
> [0]mat row indices[4] = 41761
> [0]mat row indices[5] = 41762
> [0]mat row indices[6] = 41766
> [0]mat row indices[7] = -41768
> [0]mat row indices[8] = 41767
> [0]mat row indices[9] = 41771
> [0]mat row indices[10] = -41773
> [0]mat row indices[11] = 41772
> [0]mat row indices[12] = 41776
> [0]mat row indices[13] = 41777
> [0]mat row indices[14] = 41778
> [0]mat row indices[15] = 41782
> [0]mat row indices[16] = -41784
> [0]mat row indices[17] = 41783
> [0]mat row indices[18] = 261
> [0]mat row indices[19] = -263
> [0]mat row indices[20] = 262
> [0]mat row indices[21] = 24318
> [0]mat row indices[22] = 24319
> [0]mat row indices[23] = 24320
> [0]mat row indices[24] = -7
> [0]mat row indices[25] = -8
> [0]mat row indices[26] = 6
> [0]mat row indices[27] = 1630
> [0]mat row indices[28] = -1632
> [0]mat row indices[29] = 1631
> [0]mat row indices[30] = 41757
> [0]mat row indices[31] = 41758
> [0]mat row indices[32] = 41759
> [0]mat row indices[33] = 41763
> [0]mat row indices[34] = 41764
> [0]mat row indices[35] = 41765
> [0]mat row indices[36] = 41768
> [0]mat row indices[37] = 41769
> [0]mat row indices[38] = 41770
> [0]mat row indices[39] = 41773
> [0]mat row indices[40] = 41774
> [0]mat row indices[41] = 41775
> [0]mat row indices[42] = 41779
> [0]mat row indices[43] = 41780
> [0]mat row indices[44]

Re: [petsc-users] Problems imposing boundary conditions

[Maximilian Hartig]

It seems you are correct. In theory, the problem should not be over 
constrained. It is 1/4 of a simple hollow cylinder geometry with rotational 
symmetry around the z-axis. I restrict movement completely on the upper and 
lower (z) end as well as movement in x- and y- direction respectively on the 
symmetry planes.
I am not completely sure what I am looking at with the output of 
-dm_petscsection_view. But these lines struck me as odd:

 
  (5167) dim  3 offset   0 constrained 0 1 1 2
  (5168) dim  3 offset   6 constrained 0 1 1 2
.
.
.
 (5262) dim  3 offset   0 constrained 0 0 1 2
 (5263) dim  3 offset   6 constrained 0 0 1 2


It seems that vertices that are part of the closures of both Face Sets get 
restricted twice in their respective degree of freedom.
This does however also happen when restricting movement in x- direction only 
for upper and lower faces. In that case without the solver producing an error:
  (20770) dim  3 offset  24 constrained 0 0
  (20771) dim  3 offset  30 constrained 0 0
  (20772) dim  3 offset  36 constrained 0 0
  (20773) dim  3 offset  42 constrained 0 0

Thanks,
Max

> On 6 Mar 2017, at 14:43, Matthew Knepley  wrote:
> 
> On Mon, Mar 6, 2017 at 8:38 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Of course, please find the source as well as the mesh attached below. I run 
> with:
> 
> -def_petscspace_order 2 -vel_petscspace_order 2 -snes_monitor 
> -snes_converged_reason -ksp_converged_reason -ksp_monitor _true_residual 
> -ksp_type fgmres -pc_type sor
> 
> This sounds like over-constraining a point to me. I will try and run it soon, 
> but I have a full schedule this week. The easiest
> way to see if this is happening should be to print out the Section that gets 
> made
> 
>   -dm_petscsection_view
> 
>   Thanks,
> 
>  Matt
>  
> Thanks,
> Max
> 
> 
> 
> 
>> On 4 Mar 2017, at 11:34, Sander Arens > <mailto:sander.ar...@ugent.be>> wrote:
>> 
>> Hmm, strange you also get the error in serial. Can you maybe send a minimal 
>> working which demonstrates the error?
>> 
>> Thanks,
>> Sander
>> 
>> On 3 March 2017 at 23:07, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> Yes Sander, your assessment is correct. I use DMPlex and specify the BC 
>> using DMLabel. I do however get this error also when running in serial.
>> 
>> Thanks,
>> Max
>> 
>>> On 3 Mar 2017, at 22:14, Matthew Knepley >> <mailto:knep...@gmail.com>> wrote:
>>> 
>>> On Fri, Mar 3, 2017 at 12:56 PM, Sander Arens >> <mailto:sander.ar...@ugent.be>> wrote:
>>> Max, 
>>> 
>>> I'm assuming you use DMPlex for your mesh? If so, did you only specify the 
>>> faces in the DMLabel (and not vertices or edges). Do you get this error 
>>> only in parallel? 
>>> 
>>> If so, I can confirm this bug. I submitted a pull request for this 
>>> yesterday.
>>> 
>>> Yep, I saw Sander's pull request. I will get in merged in tomorrow when I 
>>> get home to Houston.
>>> 
>>>   Thanks,
>>> 
>>>  Matt
>>>  
>>> On 3 March 2017 at 18:43, Lukas van de Wiel >> <mailto:lukas.drinkt.t...@gmail.com>> wrote:
>>> You have apparently preallocated the non-zeroes of you matrix, and the room 
>>> was insufficient to accommodate all your equations. 
>>>  
>>> What happened after you tried:
>>> 
>>> MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)
>>> 
>>> 
>>> Cheers
>>> Lukas
>>> 
>>> On Fri, Mar 3, 2017 at 6:37 PM, Maximilian Hartig >> <mailto:imilian.har...@gmail.com>> wrote:
>>> Hello,
>>> 
>>> I am working on a transient structural FEM code with PETSc. I managed to 
>>> create a slow but functioning program with the use of petscFE and a TS 
>>> solver. The code runs fine until I try to restrict movement in all three 
>>> spatial directions for one face. I then get the error which is attached 
>>> below.
>>> So apparently DMPlexMatSetClosure tries to write/read beyond what was 
>>> priorly allocated. I do however not call MatSeqAIJSetPreallocation myself 
>>> in the code. So I’m unsure where to start looking for the bug. In my 
>>> understanding, PETSc should know from the DM how much space to allocate.
>>> Could you kindly give me a hint?
>>> 
>>> Thanks,
>>> 
>>> Max
>>> 
>>> 0 SNES Function norm 2.508668036663e-06
>>> [0]PETSC ERROR: - Error Message 

Re: [petsc-users] Problems imposing boundary conditions

[Maximilian Hartig]

> On 7 Mar 2017, at 16:29, Matthew Knepley  wrote:
> 
> On Tue, Mar 7, 2017 at 3:28 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> It seems you are correct. In theory, the problem should not be over 
> constrained. It is 1/4 of a simple hollow cylinder geometry with rotational 
> symmetry around the z-axis. I restrict movement completely on the upper and 
> lower (z) end as well as movement in x- and y- direction respectively on the 
> symmetry planes.
> I am not completely sure what I am looking at with the output of 
> -dm_petscsection_view. But these lines struck me as odd:
> 
>  
>   (5167) dim  3 offset   0 constrained 0 1 1 2
>   (5168) dim  3 offset   6 constrained 0 1 1 2
> .
> .
> .
>  (5262) dim  3 offset   0 constrained 0 0 1 2
>  (5263) dim  3 offset   6 constrained 0 0 1 2
> 
> 
> It seems that vertices that are part of the closures of both Face Sets get 
> restricted twice in their respective degree of freedom.
> 
> Yes, that is exactly what happens.
>  
> This does however also happen when restricting movement in x- direction only 
> for upper and lower faces. In that case without the solver producing an error:
>   (20770) dim  3 offset  24 constrained 0 0
>   (20771) dim  3 offset  30 constrained 0 0
>   (20772) dim  3 offset  36 constrained 0 0
>   (20773) dim  3 offset  42 constrained 0 0
> 
> The fact that this does not SEGV is just luck.
> 
> Now, I did not put in any guard against this because I was not sure what 
> should happen. We could error if a local index is repeated, or
> we could ignore it. This seems unsafe if you try to constrain it with two 
> different values, but there is no way for me to tell if the values are
> compatible. Thus I just believed whatever the user told me.
> 
> What is the intention here? It would be straightforward to ignore duplicates 
> I guess.
Yes, ignoring duplicates would solve the problem then. I can think of no 
example where imposing two different Dirichlet BC on the same DOF of the same 
vertex would make sense (I might be wrong of course). That means the only issue 
is to determine wether the first or the second BC is the correct one to be 
imposed.
I don’t know how I could filter out the vertices in question from the Label. I 
use GMSH to construct my meshes and could create a label for the edges without 
too much effort. But I cannot see an easy way to exclude them when imposing the 
BC.
I tried to figure out where PETSC actually imposes the BC but got lost a bit in 
the source. Could you kindly point me towards the location?

Thanks,
Max


> 
>   Thanks,
> 
>  Matt
>  
> Thanks,
> Max
> 
>> On 6 Mar 2017, at 14:43, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Mon, Mar 6, 2017 at 8:38 AM, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> Of course, please find the source as well as the mesh attached below. I run 
>> with:
>> 
>> -def_petscspace_order 2 -vel_petscspace_order 2 -snes_monitor 
>> -snes_converged_reason -ksp_converged_reason -ksp_monitor _true_residual 
>> -ksp_type fgmres -pc_type sor
>> 
>> This sounds like over-constraining a point to me. I will try and run it 
>> soon, but I have a full schedule this week. The easiest
>> way to see if this is happening should be to print out the Section that gets 
>> made
>> 
>>   -dm_petscsection_view
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>> Thanks,
>> Max
>> 
>> 
>> 
>> 
>>> On 4 Mar 2017, at 11:34, Sander Arens >> <mailto:sander.ar...@ugent.be>> wrote:
>>> 
>>> Hmm, strange you also get the error in serial. Can you maybe send a minimal 
>>> working which demonstrates the error?
>>> 
>>> Thanks,
>>> Sander
>>> 
>>> On 3 March 2017 at 23:07, Maximilian Hartig >> <mailto:imilian.har...@gmail.com>> wrote:
>>> Yes Sander, your assessment is correct. I use DMPlex and specify the BC 
>>> using DMLabel. I do however get this error also when running in serial.
>>> 
>>> Thanks,
>>> Max
>>> 
>>>> On 3 Mar 2017, at 22:14, Matthew Knepley >>> <mailto:knep...@gmail.com>> wrote:
>>>> 
>>>> On Fri, Mar 3, 2017 at 12:56 PM, Sander Arens >>> <mailto:sander.ar...@ugent.be>> wrote:
>>>> Max, 
>>>> 
>>>> I'm assuming you use DMPlex for your mesh? If so, did you only specify the 
>>>> faces in the DMLabel (and not vertices or edges). Do you get this error 
>>>> only in parallel? 
>>>> 
>>>> 

Re: [petsc-users] Problems imposing boundary conditions

[Maximilian Hartig]

Ok thank you, so can I just do something along the lines of:
PetscSectionGetFieldConstraintDof(…)
to find the overconstrained vertices and then correct them manually with 
PetscSectionSetFieldConstraintDof()
PetscSectionSetFieldConstraintIndices()
?


Or will this mess up the Jacobian and the iFunction?

Thanks,
Max


> On 7 Mar 2017, at 18:21, Matthew Knepley  wrote:
> 
> On Tue, Mar 7, 2017 at 11:11 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> 
>> On 7 Mar 2017, at 16:29, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Tue, Mar 7, 2017 at 3:28 AM, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> It seems you are correct. In theory, the problem should not be over 
>> constrained. It is 1/4 of a simple hollow cylinder geometry with rotational 
>> symmetry around the z-axis. I restrict movement completely on the upper and 
>> lower (z) end as well as movement in x- and y- direction respectively on the 
>> symmetry planes.
>> I am not completely sure what I am looking at with the output of 
>> -dm_petscsection_view. But these lines struck me as odd:
>> 
>>  
>>   (5167) dim  3 offset   0 constrained 0 1 1 2
>>   (5168) dim  3 offset   6 constrained 0 1 1 2
>> .
>> .
>> .
>>  (5262) dim  3 offset   0 constrained 0 0 1 2
>>  (5263) dim  3 offset   6 constrained 0 0 1 2
>> 
>> 
>> It seems that vertices that are part of the closures of both Face Sets get 
>> restricted twice in their respective degree of freedom.
>> 
>> Yes, that is exactly what happens.
>>  
>> This does however also happen when restricting movement in x- direction only 
>> for upper and lower faces. In that case without the solver producing an 
>> error:
>>   (20770) dim  3 offset  24 constrained 0 0
>>   (20771) dim  3 offset  30 constrained 0 0
>>   (20772) dim  3 offset  36 constrained 0 0
>>   (20773) dim  3 offset  42 constrained 0 0
>> 
>> The fact that this does not SEGV is just luck.
>> 
>> Now, I did not put in any guard against this because I was not sure what 
>> should happen. We could error if a local index is repeated, or
>> we could ignore it. This seems unsafe if you try to constrain it with two 
>> different values, but there is no way for me to tell if the values are
>> compatible. Thus I just believed whatever the user told me.
>> 
>> What is the intention here? It would be straightforward to ignore duplicates 
>> I guess.
> 
> Yes, ignoring duplicates would solve the problem then. I can think of no 
> example where imposing two different Dirichlet BC on the same DOF of the same 
> vertex would make sense (I might be wrong of course). That means the only 
> issue is to determine wether the first or the second BC is the correct one to 
> be imposed.
> I don’t know how I could filter out the vertices in question from the Label. 
> I use GMSH to construct my meshes and could create a label for the edges 
> without too much effort. But I cannot see an easy way to exclude them when 
> imposing the BC.
> I tried to figure out where PETSC actually imposes the BC but got lost a bit 
> in the source. Could you kindly point me towards the location?
> 
> It is in stages.
> 
> 1) You make a structure with AddBoundary() that has a Label and function for 
> boundary values
> 
> 2) The PetscSection gets created with stores which points have constraints 
> and which components they affect
> 
> 3) When global Vecs are made, these constraints are left out
> 
> 4) When local Vecs are made, they are left in
> 
> 5) DMPlexInsertBoundaryValues() is called on local Vecs, and puts in the 
> values from your functions. This usually happens
> when you copy the solutions values from the global Vec to a local Vec to 
> being assembly.
> 
>   Thanks,
> 
>  Matt
>  
> Thanks,
> Max
> 
> 
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>> Thanks,
>> Max
>> 
>>> On 6 Mar 2017, at 14:43, Matthew Knepley >> <mailto:knep...@gmail.com>> wrote:
>>> 
>>> On Mon, Mar 6, 2017 at 8:38 AM, Maximilian Hartig >> <mailto:imilian.har...@gmail.com>> wrote:
>>> Of course, please find the source as well as the mesh attached below. I run 
>>> with:
>>> 
>>> -def_petscspace_order 2 -vel_petscspace_order 2 -snes_monitor 
>>> -snes_converged_reason -ksp_converged_reason -ksp_monitor _true_residual 
>>> -ksp_type fgmres -pc_type sor
>>> 
>>> This sounds like over-constraining a point to me. I will try and ru

[petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

Hello,

I’m trying to implement plasticity and have problems getting the Petsc SNES to 
converge. To check if my residual formulation is correct I tried running with 
-snes_fd for an easy example as the Petsc FAQ suggest. I cannot seem to get the 
solver to converge at any cost.
I already tried to impose bounds on the solution and moved to vinewtonrsls as a 
nonlinear solver. I checked and rechecked my residuals but I do not find an 
error there. I now have the suspicion that the -snes_fd option is not made for 
handling residuals who’s first derivatives are not continuous (e.g. have an 
“if” condition in them for the plasticity/ flow-condition). Can you confirm my 
suspicion? And is there another way to test my residual formulation separate 
from my hand-coded jacobian?


Thanks,
Max

Re: [petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

Hi Luv,

I’m modelling linear hardening(sigma_y = sigma_y0 + K_iso*epsilon_plast_eqiv) 
with isotropic  plasticity only. So I should not need to use an iterative 
method to find the point on the yield surface. I have three fields and 7 
unknowns in total:
Field 0:3 displacement components
Field 1:3 velocity components 
Field 2:1 equivalent plastic strain

It is the solver for these three fields that is not converging. I am using 
PetscFE. As residuals for the plastic case (sigma_vM > sigma_yield) I have:

Field 0 (displacement):
f0[i] = rho*u_t[u_Off[1]+i]
f1[i*dim+j] = sigma_tr[i*dim+j] - 2*mu*sqrt(3/2)*u_t[uOff[2]]*N[i*dim+j]

where sigma_tr is the trial stress, mu is the shear modulus and N is the unit 
deviator tensor normal to the yield surface.

Field 1 (velocity):
f0[i] = u[uOff[1]+i]-u_t[i]
f1[i*dim+j] = 0

Field 2 (effective plastic strain):
f0[0] = ||s_tr|| -2*mu*sqrt(3/2)*u_t[uOff[2]]-sqrt(2/3)*sigma_y
f1[i] = 0
where ||s_tr|| is the norm of the deviator stress tensor.

Field 0 residual is essentially newton’s second law of motion and Field 2 
residual should be the yield criterion. I might have just fundamentally 
misunderstood the equations of plasticity but I cannot seem to find my mistake.

Thanks,
Max


> On 30. Jun 2017, at 11:09, Luv Sharma  wrote:
> 
> Hi Max,
> 
> Is your field solver not converging or the material point solver ;)? 
> 
> Best regards,
> Luv
>> On 30 Jun 2017, at 10:45, Maximilian Hartig  wrote:
>> 
>> Hello,
>> 
>> I’m trying to implement plasticity and have problems getting the Petsc SNES 
>> to converge. To check if my residual formulation is correct I tried running 
>> with -snes_fd for an easy example as the Petsc FAQ suggest. I cannot seem to 
>> get the solver to converge at any cost.
>> I already tried to impose bounds on the solution and moved to vinewtonrsls 
>> as a nonlinear solver. I checked and rechecked my residuals but I do not 
>> find an error there. I now have the suspicion that the -snes_fd option is 
>> not made for handling residuals who’s first derivatives are not continuous 
>> (e.g. have an “if” condition in them for the plasticity/ flow-condition). 
>> Can you confirm my suspicion? And is there another way to test my residual 
>> formulation separate from my hand-coded jacobian?
>> 
>> 
>> Thanks,
>> Max
> 

Re: [petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

 iket127138 by
hartig Wed Jul  5 18:09:53 2017
[0]PETSC ERROR: Configure options PETSC_ARCH=arch-linux2-c-debug
--with-cc=mpigcc --with-cxx=mpicxx --with-fc=mpif90 --download-ml
--download-parmetis --download-metis --download-hypre --download-mumps
--download-scalapack
[0]PETSC ERROR: #1 TSStep() line 3818 in
/home/hartig/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #2 TSSolve() line 4054 in
/home/hartig/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #3 main() line 1829 in
/home/hartig/Nextcloud/Dissertation/Code/Debugging/miniFEM.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -def_petscspace_order 2
[0]PETSC ERROR: -eplast_petscspace_order 2
[0]PETSC ERROR: -ksp_converged_reason
[0]PETSC ERROR: -ksp_monitor_true_residual
[0]PETSC ERROR: -pc_type lu
[0]PETSC ERROR: -plast
[0]PETSC ERROR: -snes_converged_reason
[0]PETSC ERROR: -snes_fd
[0]PETSC ERROR: -snes_monitor
[0]PETSC ERROR: -snes_type newtonls
[0]PETSC ERROR: -ts_monitor
[0]PETSC ERROR: -vel_petscspace_order 2
[0]PETSC ERROR: End of Error Message ---send entire
error message to petsc-ma...@mcs.anl.gov--


Thanks,
Max


2017-06-30 17:14 GMT+02:00 Barry Smith :

>
>What is the output from
>
> -snes_monitor -ksp_monitor_true_residual -pc_type lu
> -ksp_converged_reason -snes_converged_reason
>
>
> > On Jun 30, 2017, at 3:45 AM, Maximilian Hartig 
> wrote:
> >
> > Hello,
> >
> > I’m trying to implement plasticity and have problems getting the Petsc
> SNES to converge. To check if my residual formulation is correct I tried
> running with -snes_fd for an easy example as the Petsc FAQ suggest. I
> cannot seem to get the solver to converge at any cost.
> > I already tried to impose bounds on the solution and moved to
> vinewtonrsls as a nonlinear solver. I checked and rechecked my residuals
> but I do not find an error there. I now have the suspicion that the
> -snes_fd option is not made for handling residuals who’s first derivatives
> are not continuous (e.g. have an “if” condition in them for the plasticity/
> flow-condition). Can you confirm my suspicion? And is there another way to
> test my residual formulation separate from my hand-coded jacobian?
> >
> >
> > Thanks,
> > Max
>
>

Re: [petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

I do not clearly understand the discrimination between local and global
plasticity. I do have areas where I expect the behaviour to be elastic and
other areas where I expect elasto-plastic behaviour.
Inertia effects are of importance and hence I need second order temporal
derivatives of my displacements. The only way I have found to implement
this in Petsc is to create a separate velocity field which I use to then
compute ü.
To account for plasticity, in my understanding I need to introduce at least
one additional history variable. In my case this is the effective plastic
strain e_p. I then solve the equation of motion (grad(sigma)-rho*ü+F=0) and
the consistency condition (sigma_eq - sigma_yield = 0) at the same time. Or
try to at least.

Thanks,
Max

2017-06-30 20:49 GMT+02:00 Luv Sharma :

> Hi Max,
>
> I do not understand the equations that you write very clearly.
>
> Are you looking to implement a “local” and “if” type of isotropic
> hardening plasticity? If that is the case, then in my understanding you
> need to solve only 1 field equation for the displacement components or for
> the strain components. You can look at the following code:
> https://github.com/tdegeus/GooseFFT/blob/master/small-
> strain/laminate/elasto-plasticity.py
>
> If you are looking for a PETSc based implementation for plasticity
> (isotropic/anisotropic) you can look at
> https://damask.mpie.de/
> I had presented a talk about the same at the PETSc User Meeting last year.
>
> As I understand it, additional field equations will only be necessary if
> the plasticity or elasticity were “nonlocal”. You may want to look at:
> On the role of moving elastic–plastic boundaries in strain gradient
> plasticity, R H J Peerlings
>
> Best regards,
> Luv
>
> On 30 Jun 2017, at 11:52, Maximilian Hartig 
> wrote:
>
> Hi Luv,
>
> I’m modelling linear hardening(sigma_y = sigma_y0 +
> K_iso*epsilon_plast_eqiv) with isotropic  plasticity only. So I should not
> need to use an iterative method to find the point on the yield surface. I
> have three fields and 7 unknowns in total:
> Field 0: 3 displacement components
> Field 1: 3 velocity components
> Field 2: 1 equivalent plastic strain
>
> It is the solver for these three fields that is not converging. I am using
> PetscFE. As residuals for the plastic case (sigma_vM > sigma_yield) I have:
>
> Field 0 (displacement):
> f0[i] = rho*u_t[u_Off[1]+i]
> f1[i*dim+j] = sigma_tr[i*dim+j] - 2*mu*sqrt(3/2)*u_t[uOff[2]]*N[i*dim+j]
>
> where sigma_tr is the trial stress, mu is the shear modulus and N is the
> unit deviator tensor normal to the yield surface.
>
> Field 1 (velocity):
> f0[i] = u[uOff[1]+i]-u_t[i]
> f1[i*dim+j] = 0
>
> Field 2 (effective plastic strain):
> f0[0] = ||s_tr|| -2*mu*sqrt(3/2)*u_t[uOff[2]]-sqrt(2/3)*sigma_y
> f1[i] = 0
> where ||s_tr|| is the norm of the deviator stress tensor.
>
> Field 0 residual is essentially newton’s second law of motion and Field 2
> residual should be the yield criterion. I might have just fundamentally
> misunderstood the equations of plasticity but I cannot seem to find my
> mistake.
>
> Thanks,
> Max
>
>
> On 30. Jun 2017, at 11:09, Luv Sharma  wrote:
>
> Hi Max,
>
> Is your field solver not converging or the material point solver ;)?
>
> Best regards,
> Luv
>
> On 30 Jun 2017, at 10:45, Maximilian Hartig 
> wrote:
>
> Hello,
>
> I’m trying to implement plasticity and have problems getting the Petsc
> SNES to converge. To check if my residual formulation is correct I tried
> running with -snes_fd for an easy example as the Petsc FAQ suggest. I
> cannot seem to get the solver to converge at any cost.
> I already tried to impose bounds on the solution and moved to vinewtonrsls
> as a nonlinear solver. I checked and rechecked my residuals but I do not
> find an error there. I now have the suspicion that the -snes_fd option is
> not made for handling residuals who’s first derivatives are not continuous
> (e.g. have an “if” condition in them for the plasticity/ flow-condition).
> Can you confirm my suspicion? And is there another way to test my residual
> formulation separate from my hand-coded jacobian?
>
>
> Thanks,
> Max
>
>
>
>
>

Re: [petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

If I set the function count higher it results in a diverged line search
error. And if I reduce the absolute residual tolerance it continues for one
or two more timesteps before it fails again.

Thanks,
Max

2017-07-05 20:24 GMT+02:00 Barry Smith :

>
> >39 SNES Function norm 3.758101652998e-08
> >   0 KSP preconditioned resid norm 1.931834799209e-02 true resid norm
> 3.758101652998e-08 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 5.958979968059e-16 true resid norm
> 2.588244438237e-18 ||r(i)||/||b|| 6.887105983875e-11
>
>The SNES function norm is pretty small here are you sure you need it
> smaller.
>
> > Nonlinear solve did not converge due to DIVERGED_FUNCTION_COUNT i
>
>You should just set the function count to say 1 billion so it doesn't
> stop the solution process.
>
>
>
>
> > On Jul 5, 2017, at 11:24 AM, Maximilian Hartig 
> wrote:
> >
> > The solver goes on for some timesteps, then it gets hung up in step 10.
> The output is as follows:
> >
> >   28 SNES Function norm 7.825649542035e-06
> >   0 KSP preconditioned resid norm 5.335581451751e+00 true resid norm
> 7.825649542035e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 3.875747966342e-13 true resid norm
> 3.082821236889e-15 ||r(i)||/||b|| 3.939380648635e-10
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >29 SNES Function norm 7.825641650401e-06
> >   0 KSP preconditioned resid norm 5.335568660325e+00 true resid norm
> 7.825641650401e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 1.169064646930e-13 true resid norm
> 2.949237553452e-15 ||r(i)||/||b|| 3.768684646199e-10
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >30 SNES Function norm 7.825633758775e-06
> >   0 KSP preconditioned resid norm 5.335575714181e+00 true resid norm
> 7.825633758775e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 3.626528075820e-13 true resid norm
> 7.454935467844e-16 ||r(i)||/||b|| 9.526302530432e-11
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >31 SNES Function norm 7.825625867144e-06
> >   0 KSP preconditioned resid norm 5.335574289757e+00 true resid norm
> 7.825625867144e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 1.220021201426e-13 true resid norm
> 6.140701297052e-16 ||r(i)||/||b|| 7.846913973787e-11
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >32 SNES Function norm 7.825617975525e-06
> >   0 KSP preconditioned resid norm 5.335556211708e+00 true resid norm
> 7.825617975525e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 3.238615954344e-14 true resid norm
> 1.639170816829e-15 ||r(i)||/||b|| 2.094621564655e-10
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >33 SNES Function norm 7.825610083913e-06
> >   0 KSP preconditioned resid norm 5.335545135100e+00 true resid norm
> 7.825610083913e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 8.328043430360e-14 true resid norm
> 7.902728010647e-16 ||r(i)||/||b|| 1.009854557780e-10
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >34 SNES Function norm 7.825609294752e-06
> >   0 KSP preconditioned resid norm 5.335519330054e+00 true resid norm
> 7.825609294752e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 1.176339150741e-13 true resid norm
> 2.137240765649e-15 ||r(i)||/||b|| 2.731085446703e-10
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >35 SNES Function norm 7.825608505592e-06
> >   0 KSP preconditioned resid norm 5.335552250936e+00 true resid norm
> 7.825608505592e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 2.388837462586e-13 true resid norm
> 3.899327292367e-16 ||r(i)||/||b|| 4.982778386601e-11
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >36 SNES Function norm 7.825607716432e-06
> >   0 KSP preconditioned resid norm 1.017554264638e-01 true resid norm
> 7.825607716432e-06 ||r(i)||/||b|| 1.e+00
> >   1 KSP preconditioned resid norm 2.758703070887e-14 true resid norm
> 1.990491888342e-16 ||r(i)||/||b|| 2.543562059932e-11
> > Linear solve converged due to CONVERGED_RTOL iterations 1
> >37 SNES Function norm 3.815837655573e-08
> >   0 KSP preconditioned resid norm 7.197730150406e+00 true resid norm
> 3.815837655573e-08 ||r(i)||/||b|| 1.e+00
> >   1 KSP 

Re: [petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

I have looked into several books for plasticity, but it cannot hurt to
check out another one. I will do so as soon as I can get a hold of it.
I do have a basic understanding on how to implement plasticity with the
return mapping algorithm. My main problem is that I am lost when trying to
implement it within Petsc. I just might have to create two separate snes
-solver contexts to run inside a single TS solver. In my understanding,
that is not how Petsc is designed to work though.
I cannot seem to find a clear formulation for dynamic plasticity either. I
looked into the Cristescu book "dynamic plasticity" but that was not very
helpful either. I have read through:
Rust, W. (2009). Nichtlineare Finite- Elemente-Berechnungen. Wiesbaden:
Springer Vieweg.
Kim, N.-H. (2015). Introduction to Nonlinear Finite Element Analysis. New
York, NY: Springer New York Heidelberg Dodrecht London.
http://doi.org/10.1007/978-1-4419-1746-1
Steinke, P. (2015). Finite-Elemente-Methode (5th ed.). Springer Vieweg.
and they all describe how to set up a stationary (or pseudo-timestepping)
solver for plasticity.
But firstly I'd like to avoid writing my own solver and secondly I need it
to be transient.

I have looked at the links you provided Luv, thanks. DAMASK is to complex a
subject for me to understand it quickly. I will need some time to
understand what it does and how to make use of it exactly.
Thanks,
Max

2017-07-05 20:49 GMT+02:00 Luv Sharma :

> Hi ,
>
> I agree with Sanjay and you can still look at the links that I mentioned;
> to get an understanding of different plasticity formulations.
> Honestly, I am also lost in plasticity ;)
>
> The first link is an FFT based spectral solver implementation of
> plasticity.
> I have a PETSc implementation (in python) of the same but without
> plasticity (only hyper-elasticity). I hope to make it open source soon.
> May be that will be helpful for you as well!
>
> Best regards,
> Luv
>
> On 5 Jul 2017, at 19:08, Sanjay Govindjee  wrote:
>
> Let me suggest that you grab a hold of Simo and Hughes, Computational
> Inelasticity, Springer-Verlag (1998).  It explains a lot about how to set
> up this problem -- in particular Chapter 1 gives a comprehensive
> one-dimensional tutorial on everything you need to know.
>
>
> On 7/5/17 9:39 AM, Maximilian Hartig wrote:
>
> I do not clearly understand the discrimination between local and global
> plasticity. I do have areas where I expect the behaviour to be elastic and
> other areas where I expect elasto-plastic behaviour.
> Inertia effects are of importance and hence I need second order temporal
> derivatives of my displacements. The only way I have found to implement
> this in Petsc is to create a separate velocity field which I use to then
> compute ü.
> To account for plasticity, in my understanding I need to introduce at
> least one additional history variable. In my case this is the effective
> plastic strain e_p. I then solve the equation of motion
> (grad(sigma)-rho*ü+F=0) and the consistency condition (sigma_eq -
> sigma_yield = 0) at the same time. Or try to at least.
>
> Thanks,
> Max
>
> 2017-06-30 20:49 GMT+02:00 Luv Sharma :
>
>> Hi Max,
>>
>> I do not understand the equations that you write very clearly.
>>
>> Are you looking to implement a “local” and “if” type of isotropic
>> hardening plasticity? If that is the case, then in my understanding you
>> need to solve only 1 field equation for the displacement components or for
>> the strain components. You can look at the following code:
>> https://github.com/tdegeus/GooseFFT/blob/master/small-strain
>> /laminate/elasto-plasticity.py
>>
>> If you are looking for a PETSc based implementation for plasticity
>> (isotropic/anisotropic) you can look at
>> https://damask.mpie.de/
>> I had presented a talk about the same at the PETSc User Meeting last year.
>>
>> As I understand it, additional field equations will only be necessary if
>> the plasticity or elasticity were “nonlocal”. You may want to look at:
>> On the role of moving elastic–plastic boundaries in strain gradient
>> plasticity, R H J Peerlings
>>
>> Best regards,
>> Luv
>>
>> On 30 Jun 2017, at 11:52, Maximilian Hartig 
>> wrote:
>>
>> Hi Luv,
>>
>> I’m modelling linear hardening(sigma_y = sigma_y0 +
>> K_iso*epsilon_plast_eqiv) with isotropic  plasticity only. So I should not
>> need to use an iterative method to find the point on the yield surface. I
>> have three fields and 7 unknowns in total:
>> Field 0: 3 displacement components
>> Field 1: 3 velocity components
>> Field 2: 1 equivalent plastic strain
>>
>> It is the solver for these three fields that is not 

Re: [petsc-users] -snes_fd for problems with residuals with non-continuous first derivative?

[Maximilian Hartig]

We do have one second order timestepper (Newmark) put in by Lisandro

Great, I was not aware of that. I'll try to make use of it as soon as I
resolve this issue.

I have a suspicion that my problem is the fact that I'm updating my plastic
strain during the nonlinear solver iterations for the stress. I tried to
create an auxiliary field for the plastic strain and then update it at the
end of each timestep with PEtscDSSetUpdate. But the update function does
not seem to get called. Do I have to update the auxiliary field in a
TSMonitor context? I couldn't find any examples of how to update an
Auxiliary field.

Thanks,
Max


2017-07-05 18:46 GMT+02:00 Matthew Knepley :

> On Wed, Jul 5, 2017 at 9:39 AM, Maximilian Hartig <
> imilian.har...@gmail.com> wrote:
>
>> I do not clearly understand the discrimination between local and global
>> plasticity. I do have areas where I expect the behaviour to be elastic and
>> other areas where I expect elasto-plastic behaviour.
>> Inertia effects are of importance and hence I need second order temporal
>> derivatives of my displacements. The only way I have found to implement
>> this in Petsc is to create a separate velocity field which I use to then
>> compute ü.
>>
>
> We do have one second order timestepper (Newmark) put in by Lisandro:
>
>   http://www.mcs.anl.gov/petsc/petsc-current/docs/
> manualpages/TS/TSALPHA2.html
>
>
>> To account for plasticity, in my understanding I need to introduce at
>> least one additional history variable. In my case this is the effective
>> plastic strain e_p. I then solve the equation of motion
>> (grad(sigma)-rho*ü+F=0) and the consistency condition (sigma_eq -
>> sigma_yield = 0) at the same time. Or try to at least.
>>
>
> Did you look at the DAMASK formulation that Luv suggested?
>
>   Thanks,
>
> Matt
>
>
>> Thanks,
>> Max
>>
>> 2017-06-30 20:49 GMT+02:00 Luv Sharma :
>>
>>> Hi Max,
>>>
>>> I do not understand the equations that you write very clearly.
>>>
>>> Are you looking to implement a “local” and “if” type of isotropic
>>> hardening plasticity? If that is the case, then in my understanding you
>>> need to solve only 1 field equation for the displacement components or for
>>> the strain components. You can look at the following code:
>>> https://github.com/tdegeus/GooseFFT/blob/master/small-strain
>>> /laminate/elasto-plasticity.py
>>>
>>> If you are looking for a PETSc based implementation for plasticity
>>> (isotropic/anisotropic) you can look at
>>> https://damask.mpie.de/
>>> I had presented a talk about the same at the PETSc User Meeting last
>>> year.
>>>
>>> As I understand it, additional field equations will only be necessary if
>>> the plasticity or elasticity were “nonlocal”. You may want to look at:
>>> On the role of moving elastic–plastic boundaries in strain gradient
>>> plasticity, R H J Peerlings
>>>
>>> Best regards,
>>> Luv
>>>
>>> On 30 Jun 2017, at 11:52, Maximilian Hartig 
>>> wrote:
>>>
>>> Hi Luv,
>>>
>>> I’m modelling linear hardening(sigma_y = sigma_y0 +
>>> K_iso*epsilon_plast_eqiv) with isotropic  plasticity only. So I should not
>>> need to use an iterative method to find the point on the yield surface. I
>>> have three fields and 7 unknowns in total:
>>> Field 0: 3 displacement components
>>> Field 1: 3 velocity components
>>> Field 2: 1 equivalent plastic strain
>>>
>>> It is the solver for these three fields that is not converging. I am
>>> using PetscFE. As residuals for the plastic case (sigma_vM > sigma_yield) I
>>> have:
>>>
>>> Field 0 (displacement):
>>> f0[i] = rho*u_t[u_Off[1]+i]
>>> f1[i*dim+j] = sigma_tr[i*dim+j] - 2*mu*sqrt(3/2)*u_t[uOff[2]]*N[i*dim+j]
>>>
>>> where sigma_tr is the trial stress, mu is the shear modulus and N is the
>>> unit deviator tensor normal to the yield surface.
>>>
>>> Field 1 (velocity):
>>> f0[i] = u[uOff[1]+i]-u_t[i]
>>> f1[i*dim+j] = 0
>>>
>>> Field 2 (effective plastic strain):
>>> f0[0] = ||s_tr|| -2*mu*sqrt(3/2)*u_t[uOff[2]]-sqrt(2/3)*sigma_y
>>> f1[i] = 0
>>> where ||s_tr|| is the norm of the deviator stress tensor.
>>>
>>> Field 0 residual is essentially newton’s second law of motion and Field
>>> 2 residual should be the yield criterion. I might have just fundamentally
>>> misunderstood the equations of plasticity but I 

[petsc-users] accessing fields from previous, converged solution.

[Maximilian Hartig]

Hello,

I have a problem with several fields that I solve with PetscFE and TS. I now 
need to access the solution from the previous timestep to compute the residual 
for the current timestep.
I tried a TSMonitor with the following code in it:

TSGetDM(ts,&dm);
DMClone(dm,&dm_aux);
DMGetDS(dm,&prob_aux);
DMSetDS(dm_aux,prob_aux);
DMCreateGlobalVector(dm_aux,&old_solution);
VecCopy(u,oldsolution);
PetscObjectCompose((PetscObject) dm, “A”, (PetscObject) old_solution);

VecDestroy(&old_solution);
DMDestroy(&dm_aux);

hoping that it would create an auxiliary field that I could access in the 
evaluation of the residual. It did that but messed with the discretisation of 
the initial problem in some way. So I figure that adding auxiliary fields to a 
dm after having fed it to a TS context is not something you should be doing.

Is there a way to access the fields of the solution for the previous timestep 
during the evaluation of the current residual?

Thanks,
Max

Re: [petsc-users] accessing fields from previous, converged solution.

[Maximilian Hartig]

Thanks for your help Matt.
> On 16. Aug 2017, at 16:17, Matthew Knepley  wrote:
> 
> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Hello,
> 
> I have a problem with several fields that I solve with PetscFE and TS. I now 
> need to access the solution from the previous timestep to compute the 
> residual for the current timestep.
> I tried a TSMonitor with the following code in it:
> 
> TSGetDM(ts,&dm);
> DMClone(dm,&dm_aux);
> DMGetDS(dm,&prob_aux);
> DMSetDS(dm_aux,prob_aux);
> DMCreateGlobalVector(dm_aux,&old_solution);
> VecCopy(u,oldsolution);
> PetscObjectCompose((PetscObject) dm, “A”, (PetscObject) old_solution);
> 
> VecDestroy(&old_solution);
> DMDestroy(&dm_aux);
> 
> hoping that it would create an auxiliary field that I could access in the 
> evaluation of the residual. It did that but messed with the discretisation of 
> the initial problem in some way. So I figure that adding auxiliary fields to 
> a dm after having fed it to a TS context is not something you should be doing.
> 
> Is there a way to access the fields of the solution for the previous timestep 
> during the evaluation of the current residual?
> 
> First, I can show you how to do what you are asking for. I think you can 
> simply
> 
> PetscObjectQuery((PetscObject) dm, "old_solution", (PetscObject *) 
> &old_solution);
> if (!old_solution) {
>   DMCreateGlobalVector(dm, &old_solution);
>   PetscObjectCompose((PetscObject) dm, "old_solution", old_solution);
> }
> VecCopy(u, oldsolution);

Unfortunately, this produces an error and tells me “An object cannot be 
composed with an object that was composed with it."

> 
> Second, I think a better way to do this than composition is to use
> 
>   DMGetNamedGlobalVector(dm, "old_solution", &old_solution);

Yes, this seems to work and I do not get an error while running the monitor. 
Also the discretisation seems to be fine. But the old solution now is not 
available as an auxiliary field.

> 
> Third, I can say that I am profoundly troubled by this. This interferes with 
> the operation of
> the time integrator, and I cannot see a reason for this. If you are keeping 
> track of the integral
> of some quantity, I would update that in TSPostStep() and request that 
> integral instead of the
> previous solution. We do this for some non-Newtonian rheologies.

I have an “if” condition in my residual which I need to evaluate to determine 
the correct formulation for my residuals. I use actual fields to keep track of 
the integrals. But when I use the actual field to evaluate the “if” condition, 
due to the implicit nature of the solver I jump “back and forth” over that 
condition. I need the “if” condition to be independent from the field for which 
it determines the residual.
The actual application is as follows:
 I have, amongst others, a displacement and a stress field.
I evaluate for my stress given a displacement increment delta u.
If the resulting stress is > yield stress, I need a plasticity formulation, if 
it is smaller, I can use elasticity.

Thanks,
Max

> 
>   Thanks,
> 
>  Matt
>  
> Thanks,
> Max
> 
> 
> 
> -- 
> 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
> 
> http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/~mk51/>

Re: [petsc-users] accessing fields from previous, converged solution.

[Maximilian Hartig]

> On 17. Aug 2017, at 13:51, Matthew Knepley  wrote:
> 
> On Thu, Aug 17, 2017 at 3:13 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Thanks for your help Matt.
>> On 16. Aug 2017, at 16:17, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> Hello,
>> 
>> I have a problem with several fields that I solve with PetscFE and TS. I now 
>> need to access the solution from the previous timestep to compute the 
>> residual for the current timestep.
>> I tried a TSMonitor with the following code in it:
>> 
>> TSGetDM(ts,&dm);
>> DMClone(dm,&dm_aux);
>> DMGetDS(dm,&prob_aux);
>> DMSetDS(dm_aux,prob_aux);
>> DMCreateGlobalVector(dm_aux,&old_solution);
>> VecCopy(u,oldsolution);
>> PetscObjectCompose((PetscObject) dm, “A”, (PetscObject) old_solution);
>> 
>> VecDestroy(&old_solution);
>> DMDestroy(&dm_aux);
>> 
>> hoping that it would create an auxiliary field that I could access in the 
>> evaluation of the residual. It did that but messed with the discretisation 
>> of the initial problem in some way. So I figure that adding auxiliary fields 
>> to a dm after having fed it to a TS context is not something you should be 
>> doing.
>> 
>> Is there a way to access the fields of the solution for the previous 
>> timestep during the evaluation of the current residual?
>> 
>> First, I can show you how to do what you are asking for. I think you can 
>> simply
>> 
>> PetscObjectQuery((PetscObject) dm, "old_solution", (PetscObject *) 
>> &old_solution);
>> if (!old_solution) {
>>   DMCreateGlobalVector(dm, &old_solution);
>>   PetscObjectCompose((PetscObject) dm, "old_solution", old_solution);
>> }
>> VecCopy(u, oldsolution);
> 
> Unfortunately, this produces an error and tells me “An object cannot be 
> composed with an object that was composed with it."
> 
>> 
>> Second, I think a better way to do this than composition is to use
>> 
>>   DMGetNamedGlobalVector(dm, "old_solution", &old_solution);
> 
> Yes, this seems to work and I do not get an error while running the monitor. 
> Also the discretisation seems to be fine. But the old solution now is not 
> available as an auxiliary field.
> 
>> 
>> Third, I can say that I am profoundly troubled by this. This interferes with 
>> the operation of
>> the time integrator, and I cannot see a reason for this. If you are keeping 
>> track of the integral
>> of some quantity, I would update that in TSPostStep() and request that 
>> integral instead of the
>> previous solution. We do this for some non-Newtonian rheologies.
> 
> I have an “if” condition in my residual which I need to evaluate to determine 
> the correct formulation for my residuals. I use actual fields to keep track 
> of the integrals. But when I use the actual field to evaluate the “if” 
> condition, due to the implicit nature of the solver I jump “back and forth” 
> over that condition. I need the “if” condition to be independent from the 
> field for which it determines the residual.
> The actual application is as follows:
>  I have, amongst others, a displacement and a stress field.
> I evaluate for my stress given a displacement increment delta u.
> If the resulting stress is > yield stress, I need a plasticity formulation, 
> if it is smaller, I can use elasticity.
> 
> Hmm, I think this is a formulation problem. You should be using the "correct" 
> formulation at the implicit step, since otherwise
> the residual will be inconsistent with what you tested when evaluated at the 
> new step. I have the same kind of elasto-plastic
> system in PyLith (http://www.geodynamics.org/cig/software/pylith/ 
> <http://www.geodynamics.org/cig/software/pylith/>) which we solve implicitly 
> without much problem.

Possible, but I rechecked the formulation and it looks fine to me. In 
literature, this kind of algorithm uses an elastic predictor step to determine 
whether to use the elastic or the plastic formulation. The elastic “trial” 
stress is compared with the yield stress from the past timestep. In my current 
setup, I do not compare the elastic “trial” stress with the past yield stress 
but the current stress (including plastic correction) with the current yield 
stress. I think this is what leads to my problem. When I cheat and just tell 
petsc to use the plastic formulation from the time on where I expect plastic 
behaviour, I get convergence and the results look reaso

Re: [petsc-users] accessing fields from previous, converged solution.

[Maximilian Hartig]

Sorry to bother you again, but I have checked and rechecked my formulation and 
it corresponds to what I have found in literature. The only difference is that 
my flow criteria is not evaluated using the trial stress and the current yield 
stress but rather the corrected stress and the updated yield stress. I am 
convinced this is what causes the problem.
Unfortunately I do not have any idea on how to go about implementing this. 
First I don’t know how to access the solution from the past iteration and 
secondly I don’t know how to do the elastic predictor step. I tried doing it 
with auxiliary fields that get updated every timestep but I gather that is not 
a good idea. Could you please point me in some direction?  

I appreciate you patience with my questioning.

Thanks,
Max

> On 22. Aug 2017, at 13:52, Maximilian Hartig  wrote:
> 
>> 
>> On 17. Aug 2017, at 13:51, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Thu, Aug 17, 2017 at 3:13 AM, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> Thanks for your help Matt.
>>> On 16. Aug 2017, at 16:17, Matthew Knepley >> <mailto:knep...@gmail.com>> wrote:
>>> 
>>> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig 
>>> mailto:imilian.har...@gmail.com>> wrote:
>>> Hello,
>>> 
>>> I have a problem with several fields that I solve with PetscFE and TS. I 
>>> now need to access the solution from the previous timestep to compute the 
>>> residual for the current timestep.
>>> I tried a TSMonitor with the following code in it:
>>> 
>>> TSGetDM(ts,&dm);
>>> DMClone(dm,&dm_aux);
>>> DMGetDS(dm,&prob_aux);
>>> DMSetDS(dm_aux,prob_aux);
>>> DMCreateGlobalVector(dm_aux,&old_solution);
>>> VecCopy(u,oldsolution);
>>> PetscObjectCompose((PetscObject) dm, “A”, (PetscObject) old_solution);
>>> 
>>> VecDestroy(&old_solution);
>>> DMDestroy(&dm_aux);
>>> 
>>> hoping that it would create an auxiliary field that I could access in the 
>>> evaluation of the residual. It did that but messed with the discretisation 
>>> of the initial problem in some way. So I figure that adding auxiliary 
>>> fields to a dm after having fed it to a TS context is not something you 
>>> should be doing.
>>> 
>>> Is there a way to access the fields of the solution for the previous 
>>> timestep during the evaluation of the current residual?
>>> 
>>> First, I can show you how to do what you are asking for. I think you can 
>>> simply
>>> 
>>> PetscObjectQuery((PetscObject) dm, "old_solution", (PetscObject *) 
>>> &old_solution);
>>> if (!old_solution) {
>>>   DMCreateGlobalVector(dm, &old_solution);
>>>   PetscObjectCompose((PetscObject) dm, "old_solution", old_solution);
>>> }
>>> VecCopy(u, oldsolution);
>> 
>> Unfortunately, this produces an error and tells me “An object cannot be 
>> composed with an object that was composed with it."
>> 
>>> 
>>> Second, I think a better way to do this than composition is to use
>>> 
>>>   DMGetNamedGlobalVector(dm, "old_solution", &old_solution);
>> 
>> Yes, this seems to work and I do not get an error while running the monitor. 
>> Also the discretisation seems to be fine. But the old solution now is not 
>> available as an auxiliary field.
>> 
>>> 
>>> Third, I can say that I am profoundly troubled by this. This interferes 
>>> with the operation of
>>> the time integrator, and I cannot see a reason for this. If you are keeping 
>>> track of the integral
>>> of some quantity, I would update that in TSPostStep() and request that 
>>> integral instead of the
>>> previous solution. We do this for some non-Newtonian rheologies.
>> 
>> I have an “if” condition in my residual which I need to evaluate to 
>> determine the correct formulation for my residuals. I use actual fields to 
>> keep track of the integrals. But when I use the actual field to evaluate the 
>> “if” condition, due to the implicit nature of the solver I jump “back and 
>> forth” over that condition. I need the “if” condition to be independent from 
>> the field for which it determines the residual.
>> The actual application is as follows:
>>  I have, amongst others, a displacement and a stress field.
>> I evaluate for my stress given a displacement increment delta u.
>> If the resulting stress is > yield stress, I need a 

Re: [petsc-users] accessing fields from previous, converged solution.

[Maximilian Hartig]

> On 28. Aug 2017, at 12:31, Matthew Knepley  wrote:
> 
> On Mon, Aug 28, 2017 at 5:15 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Sorry to bother you again, but I have checked and rechecked my formulation 
> and it corresponds to what I have found in literature. The only difference is 
> that my flow criteria is not evaluated using the trial stress and the current 
> yield stress but rather the corrected stress and the updated yield stress. I 
> am convinced this is what causes the problem.
> Unfortunately I do not have any idea on how to go about implementing this. 
> First I don’t know how to access the solution from the past iteration
> 
> I think you can do what you want using
> 
>   
> http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPostStep.html
>  
> <http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPostStep.html>
> 
> In your function, you would TSGetSolution() and then copy it into an 
> auxiliary vector you
> set. This avoids the problem of composition you had before.

ok, so now I have a tspoststep function in which I do:

TSGetDM(ts, &dm);
PetscObjcetQuery((PetscObject) dm, “A”, (PetscObject*) &copied_solution);
TSGetSolution(ts, &solution);
VecCopy(solution, copied_solution);

but Petsc complains that the vectors have different sizes. I guess this is due 
to the Dirichlet BC I imposed. I tried imposing the same boundary condition on 
the auxiliary field which made the vectors the same length. However it seemed 
to mess up the way auxiliary fields get accessed. Any idea on how I can get the 
“full” solution vector where the boundary values are filled in? I tried writing 
it to a hdf5- file and then calling VecRead on it but I cannot get the hdf5 
file format to work.
>  
> and secondly I don’t know how to do the elastic predictor step.
> 
> You can always create another SNES and do a subsolve, maybe in
> 
>   
> http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPreStep.html#TSSetPreStep
>  
> <http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPreStep.html#TSSetPreStep>
That seems very useful as well. Are there any examples for the sub solve?

Thanks,
Max
> 
> which you then copy into an auxiliary vector, just as above.
> 
> Tell me if I am missing something.
> 
>   Thanks,
> 
>  Matt
>  
> I tried doing it with auxiliary fields that get updated every timestep but I 
> gather that is not a good idea. Could you please point me in some direction?  
> 
> I appreciate you patience with my questioning.
> 
> Thanks,
> Max
> 
>> On 22. Aug 2017, at 13:52, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> 
>>> 
>>> On 17. Aug 2017, at 13:51, Matthew Knepley >> <mailto:knep...@gmail.com>> wrote:
>>> 
>>> On Thu, Aug 17, 2017 at 3:13 AM, Maximilian Hartig 
>>> mailto:imilian.har...@gmail.com>> wrote:
>>> Thanks for your help Matt.
>>>> On 16. Aug 2017, at 16:17, Matthew Knepley >>> <mailto:knep...@gmail.com>> wrote:
>>>> 
>>>> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig 
>>>> mailto:imilian.har...@gmail.com>> wrote:
>>>> Hello,
>>>> 
>>>> I have a problem with several fields that I solve with PetscFE and TS. I 
>>>> now need to access the solution from the previous timestep to compute the 
>>>> residual for the current timestep.
>>>> I tried a TSMonitor with the following code in it:
>>>> 
>>>> TSGetDM(ts,&dm);
>>>> DMClone(dm,&dm_aux);
>>>> DMGetDS(dm,&prob_aux);
>>>> DMSetDS(dm_aux,prob_aux);
>>>> DMCreateGlobalVector(dm_aux,&old_solution);
>>>> VecCopy(u,oldsolution);
>>>> PetscObjectCompose((PetscObject) dm, “A”, (PetscObject) old_solution);
>>>> 
>>>> VecDestroy(&old_solution);
>>>> DMDestroy(&dm_aux);
>>>> 
>>>> hoping that it would create an auxiliary field that I could access in the 
>>>> evaluation of the residual. It did that but messed with the discretisation 
>>>> of the initial problem in some way. So I figure that adding auxiliary 
>>>> fields to a dm after having fed it to a TS context is not something you 
>>>> should be doing.
>>>> 
>>>> Is there a way to access the fields of the solution for the previous 
>>>> timestep during the evaluation of the current residual?
>>>> 
>>>> First, I can show you how to do what you are asking for. I think you can 
>>>> simply
&g

[petsc-users] Gradient in DMProjectField

[Maximilian Hartig]

Hello,

in my attempt to use an incremental formulation for elastoplasticity I wanted 
to create an auxiliary field in which the stresses from the past tilmestep are 
stored. I managed to create an auxiliary DM on which I have a field for the 
displacement, one for the stress and one for the inner variables. I transfer 
the displacements to a global vector on the auxiliary field, then calculate the 
stresses and inner variables using the DMProjectField function. Transfer of the 
displacements works fine and they are correct. The problem is with the 
stresses. While on some gauss points I get the correct stress, others have 
random stresses. I construct the strain using the gradient of the displacements 
and I believe this is where the problem lies. However, I have not been able to 
pinpoint it.
When I just transfer the displacements to the auxiliary fields and then use a_x 
in the residual function instead, the stresses are correct. Is there a trick 
with using u_x in the DMProjectField?

Thanks,
Max

Re: [petsc-users] Gradient in DMProjectField

[Maximilian Hartig]

On 12. Sep 2017, at 21:26, Matthew Knepley <knep...@gmail.com> wrote:On Tue, Sep 12, 2017 at 12:52 PM, Maximilian Hartig <imilian.har...@gmail.com> wrote:Hello,

in my attempt to use an incremental formulation for elastoplasticity I wanted to create an auxiliary field in which the stresses from the past tilmestep are stored. I managed to create an auxiliary DM on which I have a field for the displacement, one for the stress and one for the inner variables. I transfer the displacements to a global vector on the auxiliary field, then calculate the stresses and inner variables using the DMProjectField function. Transfer of the displacements works fine and they are correct. The problem is with the stresses. While on some gauss points I get the correct stress, others have random stresses. I construct the strain using the gradient of the displacements and I believe this is where the problem lies. However, I have not been able to pinpoint it.
When I just transfer the displacements to the auxiliary fields and then use a_x in the residual function instead, the stresses are correct. Is there a trick with using u_x in the DMProjectField?Its possible that there is a bug. Just so I understand, you are saying that DMProjectField() has the wrong gradient of the input Vec?It would really be great if you could demonstrate it with a standalone example.That is correct, and I suspect it has something to do with the Dirichlet boundary conditions I impose.If I comment out the following line, projected result and expected result are the same: ierr = DMAddBoundary(dm, PETSC_TRUE, "fixed", "Face Sets",0,Ncomp,restrictall,(void (*)()) zero_vector, Nfid,fid,NULL);CHKERRQ(ierr);Please find attached the test example and the test mesh. I run it with: -disp_petscspace_order 2 -stress_petscspace_order 2Thanks,Max

projectfieldtest.c
Description: Binary data


testcube.msh
Description: Binary data
  Thanks,    Matt 
Thanks,
Max-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.-- Norbert Wienerhttp://www.caam.rice.edu/~mk51/

[petsc-users] Hints for using petscfe for plasticity -- how to update/access internal variables?

[Maximilian Hartig]

Hello,

I’m trying to implement plasticity using petscFE but I am quite stuck since a 
while. Here’s what I’m trying to do:

I have a TS which solves the following equation:
gradient(stress) +Forces = density*acceleration
where at the moment stress is a linear function of the strain and hence the 
gradient of the displacement. This works fine. Now I want to compare the stress 
to a reference value and if it lies above this yield stress, I have to 
reevaluate the stress at the respective location. Then I need to update the 
plastic strain / yield stress  at this location.
I tried doing that first by solving three fields at the same time: 
displacements, stresses and yield stress. This failed.
Then, I tried solving only for displacement increments, storing the 
displacements, stresses and yield stress from the past time step in an 
auxiliary field. The auxiliary fields are updated after each time step with a 
second SNES, using the displacement increments from the current, converged time 
step. This also failed.
In both cases the code had problems converging and when it did, I ended up with 
negative plastic strain. This is not possible and I don’t know how it happens 
because I explicitly only increment the plastic strain when the increment is 
positive.

I am sure there is an easy solution to how I can update the internal variables 
and determine the correct stress for the residual but I just cannot figure it 
out. I’d be thankful for any hints.

Thanks,
Max

Re: [petsc-users] Hints for using petscfe for plasticity -- how to update/access internal variables?

[Maximilian Hartig]

> On 20. Sep 2017, at 18:17, Matthew Knepley  wrote:
> 
> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Hello,
> 
> I’m trying to implement plasticity using petscFE but I am quite stuck since a 
> while. Here’s what I’m trying to do:
> 
> I have a TS which solves the following equation:
> gradient(stress) +Forces = density*acceleration
> where at the moment stress is a linear function of the strain and hence the 
> gradient of the displacement. This works fine. Now I want to compare the 
> stress to a reference value and if it lies above this yield stress, I have to 
> reevaluate the stress at the respective location. Then I need to update the 
> plastic strain / yield stress  at this location.
> I tried doing that first by solving three fields at the same time: 
> displacements, stresses and yield stress. This failed.
> Then, I tried solving only for displacement increments, storing the 
> displacements, stresses and yield stress from the past time step in an 
> auxiliary field. The auxiliary fields are updated after each time step with a 
> second SNES, using the displacement increments from the current, converged 
> time step. This also failed.
> In both cases the code had problems converging and when it did, I ended up 
> with negative plastic strain. This is not possible and I don’t know how it 
> happens because I explicitly only increment the plastic strain when the 
> increment is positive.
> 
> I am sure there is an easy solution to how I can update the internal 
> variables and determine the correct stress for the residual but I just cannot 
> figure it out. I’d be thankful for any hints.
> 
> It looks like there are two problems above:
> 
> 1) Convergence
> 
> For any convergence question, we at minimum need to see the output of
> 
>   -snes_view -snes_converged_reason -snes_monitor -ksp_monitor_true_residual 
> -snes_linesearch_monitor
> 
> However, this does not seem to be the main issue.
> 
> 2) Negative plastic strain

This is what I’m mainly concerned with.
> 
> If the system really converged (I cannot tell without other information), 
> then the system formulation is wrong. Of course, its
> really easy to check by just plugging your solution into the residual 
> function too. I do not understand your explanation above
> completely however. Do you solve for the plastic strain or the increment?

I am trying to find a formulation that works and I think there is a core 
concept I am just not “getting”. 
I want to solve for the displacements. 
This works fine in an elastic case. When plasticity is involved, I need to 
determine the actual stress for my residual evaluation and I have not found a 
way to do that.
All formulations for stress I found in literature use strain increments so I 
tried to just solve for increments each timestep and then add them together in 
tspoststep. But I still need to somehow evaluate the stress for my displacement 
increment residuals. So currently, I have auxiliary fields with the stress and 
the plastic strain. I evaluate the current trial stress by adding a stress 
increment assuming elastic behaviour. If the trial stress lies beyond the yield 
stress I calculate the corrected stress to evaluate my residual for the 
displacements. But now I somehow need to update my plastic strain and the 
stress in the auxiliary fields. So in tspoststep I created another SNES to now 
calculate the stress and plastic strain while the displacement is the auxiliary 
field. 

I’m sure there’s an elegant solution on how to update internal variables but I 
have not found it.

Thanks,
Max
> 
>   Thanks,
> 
>  Matt
>  
> Thanks,
> Max
> 
> 
> 
> -- 
> 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
> 
> http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/~mk51/>

Re: [petsc-users] Hints for using petscfe for plasticity -- how to update/access internal variables?

[Maximilian Hartig]

> On 20. Sep 2017, at 19:05, Matthew Knepley  wrote:
> 
> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
>> On 20. Sep 2017, at 18:17, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig 
>> mailto:imilian.har...@gmail.com>> wrote:
>> Hello,
>> 
>> I’m trying to implement plasticity using petscFE but I am quite stuck since 
>> a while. Here’s what I’m trying to do:
>> 
>> I have a TS which solves the following equation:
>> gradient(stress) +Forces = density*acceleration
>> where at the moment stress is a linear function of the strain and hence the 
>> gradient of the displacement. This works fine. Now I want to compare the 
>> stress to a reference value and if it lies above this yield stress, I have 
>> to reevaluate the stress at the respective location. Then I need to update 
>> the plastic strain / yield stress  at this location.
>> I tried doing that first by solving three fields at the same time: 
>> displacements, stresses and yield stress. This failed.
>> Then, I tried solving only for displacement increments, storing the 
>> displacements, stresses and yield stress from the past time step in an 
>> auxiliary field. The auxiliary fields are updated after each time step with 
>> a second SNES, using the displacement increments from the current, converged 
>> time step. This also failed.
>> In both cases the code had problems converging and when it did, I ended up 
>> with negative plastic strain. This is not possible and I don’t know how it 
>> happens because I explicitly only increment the plastic strain when the 
>> increment is positive.
>> 
>> I am sure there is an easy solution to how I can update the internal 
>> variables and determine the correct stress for the residual but I just 
>> cannot figure it out. I’d be thankful for any hints.
>> 
>> It looks like there are two problems above:
>> 
>> 1) Convergence
>> 
>> For any convergence question, we at minimum need to see the output of
>> 
>>   -snes_view -snes_converged_reason -snes_monitor -ksp_monitor_true_residual 
>> -snes_linesearch_monitor
>> 
>> However, this does not seem to be the main issue.
>> 
>> 2) Negative plastic strain
> 
> This is what I’m mainly concerned with.
>> 
>> If the system really converged (I cannot tell without other information), 
>> then the system formulation is wrong. Of course, its
>> really easy to check by just plugging your solution into the residual 
>> function too. I do not understand your explanation above
>> completely however. Do you solve for the plastic strain or the increment?
> 
> I am trying to find a formulation that works and I think there is a core 
> concept I am just not “getting”. 
> I want to solve for the displacements. 
> This works fine in an elastic case. When plasticity is involved, I need to 
> determine the actual stress for my residual evaluation and I have not found a 
> way to do that.
> All formulations for stress I found in literature use strain increments so I 
> tried to just solve for increments each timestep and then add them together 
> in tspoststep. But I still need to somehow evaluate the stress for my 
> displacement increment residuals. So currently, I have auxiliary fields with 
> the stress and the plastic strain.
> 
> First question: Don't you get stress by just applying a local operator, 
> rather than a solve?
That depends on the type of plasticity. For a linear hardening formulation it 
is correct that I could just apply a local operator. I’d be happy with that for 
now. But I’d still need to save stress state and plastic strain to determine 
whether or not I’m dealing with a plasticity. I don’t know how to do that 
inside the residual evaluation. Plus DMProjectField seems to have problems 
evaluating the gradient when boundary conditions are imposed.

Thanks,
Max
> 
>   Thanks,
> 
>  Matt
>  
> I evaluate the current trial stress by adding a stress increment assuming 
> elastic behaviour. If the trial stress lies beyond the yield stress I 
> calculate the corrected stress to evaluate my residual for the displacements. 
> But now I somehow need to update my plastic strain and the stress in the 
> auxiliary fields. So in tspoststep I created another SNES to now calculate 
> the stress and plastic strain while the displacement is the auxiliary field. 
> 
> I’m sure there’s an elegant solution on how to update internal variables but 
> I have not found it.
> 
> Thanks,
> Max
>> 
>>   Thanks,
>> 
>

Re: [petsc-users] Hints for using petscfe for plasticity -- how to update/access internal variables?

[Maximilian Hartig]

> On 20. Sep 2017, at 22:51, Matthew Knepley  wrote:
> 
> On Wed, Sep 20, 2017 at 3:46 PM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> 
>> On 20. Sep 2017, at 19:05, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig 
>> mailto:imilian.har...@gmail.com>> wrote:
>>> On 20. Sep 2017, at 18:17, Matthew Knepley >> <mailto:knep...@gmail.com>> wrote:
>>> 
>>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig 
>>> mailto:imilian.har...@gmail.com>> wrote:
>>> Hello,
>>> 
>>> I’m trying to implement plasticity using petscFE but I am quite stuck since 
>>> a while. Here’s what I’m trying to do:
>>> 
>>> I have a TS which solves the following equation:
>>> gradient(stress) +Forces = density*acceleration
>>> where at the moment stress is a linear function of the strain and hence the 
>>> gradient of the displacement. This works fine. Now I want to compare the 
>>> stress to a reference value and if it lies above this yield stress, I have 
>>> to reevaluate the stress at the respective location. Then I need to update 
>>> the plastic strain / yield stress  at this location.
>>> I tried doing that first by solving three fields at the same time: 
>>> displacements, stresses and yield stress. This failed.
>>> Then, I tried solving only for displacement increments, storing the 
>>> displacements, stresses and yield stress from the past time step in an 
>>> auxiliary field. The auxiliary fields are updated after each time step with 
>>> a second SNES, using the displacement increments from the current, 
>>> converged time step. This also failed.
>>> In both cases the code had problems converging and when it did, I ended up 
>>> with negative plastic strain. This is not possible and I don’t know how it 
>>> happens because I explicitly only increment the plastic strain when the 
>>> increment is positive.
>>> 
>>> I am sure there is an easy solution to how I can update the internal 
>>> variables and determine the correct stress for the residual but I just 
>>> cannot figure it out. I’d be thankful for any hints.
>>> 
>>> It looks like there are two problems above:
>>> 
>>> 1) Convergence
>>> 
>>> For any convergence question, we at minimum need to see the output of
>>> 
>>>   -snes_view -snes_converged_reason -snes_monitor 
>>> -ksp_monitor_true_residual -snes_linesearch_monitor
>>> 
>>> However, this does not seem to be the main issue.
>>> 
>>> 2) Negative plastic strain
>> 
>> This is what I’m mainly concerned with.
>>> 
>>> If the system really converged (I cannot tell without other information), 
>>> then the system formulation is wrong. Of course, its
>>> really easy to check by just plugging your solution into the residual 
>>> function too. I do not understand your explanation above
>>> completely however. Do you solve for the plastic strain or the increment?
>> 
>> I am trying to find a formulation that works and I think there is a core 
>> concept I am just not “getting”. 
>> I want to solve for the displacements. 
>> This works fine in an elastic case. When plasticity is involved, I need to 
>> determine the actual stress for my residual evaluation and I have not found 
>> a way to do that.
>> All formulations for stress I found in literature use strain increments so I 
>> tried to just solve for increments each timestep and then add them together 
>> in tspoststep. But I still need to somehow evaluate the stress for my 
>> displacement increment residuals. So currently, I have auxiliary fields with 
>> the stress and the plastic strain.
>> 
>> First question: Don't you get stress by just applying a local operator, 
>> rather than a solve?
> That depends on the type of plasticity.
> 
> What type of plasticity is not local?
I did not express myself correctly. The plasticity is local, but for nonlinear 
hardening I would have to iteratively solve to get the correct stress and 
plastic strain from the displacement increments.
>  
> For a linear hardening formulation it is correct that I could just apply a 
> local operator. I’d be happy with that for now. But I’d still need to save 
> stress state and plastic strain to determine whether or not I’m dealing with 
> a plasticity. I don’t know how to do that inside the residual evaluation.
> 
> I do not know what you mean b

Re: [petsc-users] Hints for using petscfe for plasticity -- how to update/access internal variables?

[Maximilian Hartig]

> On 21. Sep 2017, at 16:02, Mark Adams  wrote:
> 
>>> 
> 
> That depends on the type of plasticity. For a linear hardening formulation it 
> is correct that I could just apply a local operator. I’d be happy with that 
> for now.
> 
> I would just do this for now and get it working.


I agree, and the local version would be enough for me to continue for now. The 
issue was that I could not get DMProjectField to work to calculate my stress so 
I decided to just do a solve and set the residual of the fields to 
f0[]=u[]-calculatedstress[]
and
f0[] = u[uOff[1]]-calculatedplasticstrain[]
respectively.

My calculatedplasticstrain was equal to the one from the last step in case of 
elastic behaviour. In case of plastic behaviour it was the one from the last 
step plus a positive increment. Yet I ended up with negative plastic strains.

I wouldn’t need to do a global solve if I managed to just project the correct 
stresses and strains in the auxiliary field. I don’t even want to do it. I just 
couldn’t get DMProjectFields to give me the correct gradients and tried this as 
a workaround.
>  
> But I’d still need to save stress state and plastic strain to determine 
> whether or not I’m dealing with a plasticity. I don’t know how to do that 
> inside the residual evaluation. Plus DMProjectField seems to have problems 
> evaluating the gradient when boundary conditions are imposed.
> 
> After you get the local version working try your nonlocal thing without BCs. 
> You may find bugs while you do this that are causing your problem, and if not 
> we can think more about it.

I don’t understand what you mean by “without BC”. If I don’t impose boundary 
conditions there is nothing to solve. The DMProjectField thing does work 
without boundary conditions if you mean that. It is only when I impose 
Dirichlet BC on the dm that I get the wrong gradients. 
This is also true if I use a second, unconstrained dm_unconstrained to project 
the fields but the original vector comes from a constrained dm.
I also tried creating a second vector on dm_unconstrained and using 
DMGlobalToLocal to transfer from the original vector to the second, 
unconstrained one and then this one as input for DMProjectField. It does not 
prevent the gradient issue from happening. 

I use petsc from bitbucket origin/master


Thanks,
Max

 
>  
> 
> Thanks,
> Max
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>> I evaluate the current trial stress by adding a stress increment assuming 
>> elastic behaviour. If the trial stress lies beyond the yield stress I 
>> calculate the corrected stress to evaluate my residual for the 
>> displacements. But now I somehow need to update my plastic strain and the 
>> stress in the auxiliary fields. So in tspoststep I created another SNES to 
>> now calculate the stress and plastic strain while the displacement is the 
>> auxiliary field. 
>> 
>> I’m sure there’s an elegant solution on how to update internal variables but 
>> I have not found it.
>> 
>> Thanks,
>> Max
>>> 
>>>   Thanks,
>>> 
>>>  Matt
>>>  
>>> Thanks,
>>> Max
>>> 
>>> 
>>> 
>>> -- 
>>> 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
>>> 
>>> http://www.caam.rice.edu/~mk51/ 
>> 
>> 
>> 
>> -- 
>> 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
>> 
>> http://www.caam.rice.edu/~mk51/ 
> 

Re: [petsc-users] Hints for using petscfe for plasticity -- how to update/access internal variables?

[Maximilian Hartig]

> On 21. Sep 2017, at 17:33, Matthew Knepley  wrote:
> 
> On Thu, Sep 21, 2017 at 11:28 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
>> On 20. Sep 2017, at 22:51, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Wed, Sep 20, 2017 at 3:46 PM, Maximilian Hartig > <mailto:imilian.har...@gmail.com>> wrote:
>> 
>>> On 20. Sep 2017, at 19:05, Matthew Knepley >> <mailto:knep...@gmail.com>> wrote:
>>> 
>>> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig 
>>> mailto:imilian.har...@gmail.com>> wrote:
>>>> On 20. Sep 2017, at 18:17, Matthew Knepley >>> <mailto:knep...@gmail.com>> wrote:
>>>> 
>>>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig 
>>>> mailto:imilian.har...@gmail.com>> wrote:
>>>> Hello,
>>>> 
>>>> I’m trying to implement plasticity using petscFE but I am quite stuck 
>>>> since a while. Here’s what I’m trying to do:
>>>> 
>>>> I have a TS which solves the following equation:
>>>> gradient(stress) +Forces = density*acceleration
>>>> where at the moment stress is a linear function of the strain and hence 
>>>> the gradient of the displacement. This works fine. Now I want to compare 
>>>> the stress to a reference value and if it lies above this yield stress, I 
>>>> have to reevaluate the stress at the respective location. Then I need to 
>>>> update the plastic strain / yield stress  at this location.
>>>> I tried doing that first by solving three fields at the same time: 
>>>> displacements, stresses and yield stress. This failed.
>>>> Then, I tried solving only for displacement increments, storing the 
>>>> displacements, stresses and yield stress from the past time step in an 
>>>> auxiliary field. The auxiliary fields are updated after each time step 
>>>> with a second SNES, using the displacement increments from the current, 
>>>> converged time step. This also failed.
>>>> In both cases the code had problems converging and when it did, I ended up 
>>>> with negative plastic strain. This is not possible and I don’t know how it 
>>>> happens because I explicitly only increment the plastic strain when the 
>>>> increment is positive.
>>>> 
>>>> I am sure there is an easy solution to how I can update the internal 
>>>> variables and determine the correct stress for the residual but I just 
>>>> cannot figure it out. I’d be thankful for any hints.
>>>> 
>>>> It looks like there are two problems above:
>>>> 
>>>> 1) Convergence
>>>> 
>>>> For any convergence question, we at minimum need to see the output of
>>>> 
>>>>   -snes_view -snes_converged_reason -snes_monitor 
>>>> -ksp_monitor_true_residual -snes_linesearch_monitor
>>>> 
>>>> However, this does not seem to be the main issue.
>>>> 
>>>> 2) Negative plastic strain
>>> 
>>> This is what I’m mainly concerned with.
>>>> 
>>>> If the system really converged (I cannot tell without other information), 
>>>> then the system formulation is wrong. Of course, its
>>>> really easy to check by just plugging your solution into the residual 
>>>> function too. I do not understand your explanation above
>>>> completely however. Do you solve for the plastic strain or the increment?
>>> 
>>> I am trying to find a formulation that works and I think there is a core 
>>> concept I am just not “getting”. 
>>> I want to solve for the displacements. 
>>> This works fine in an elastic case. When plasticity is involved, I need to 
>>> determine the actual stress for my residual evaluation and I have not found 
>>> a way to do that.
>>> All formulations for stress I found in literature use strain increments so 
>>> I tried to just solve for increments each timestep and then add them 
>>> together in tspoststep. But I still need to somehow evaluate the stress for 
>>> my displacement increment residuals. So currently, I have auxiliary fields 
>>> with the stress and the plastic strain.
>>> 
>>> First question: Don't you get stress by just applying a local operator, 
>>> rather than a solve?
>> That depends on the type of plasticity.
>> 
>> What type of plasticity is not local?
> I did not express myself correctly. Th

[petsc-users] TS and petscFE

[Maximilian Hartig]

Hello all,

I would like to run a transient problem with PetscFE. Example ex11.c seems 
relevant since it uses the PestcFV context to create boundary conditions and 
RHS Functions for the TS.
Is there an easy way to do transient analysis with TS and petscFE or do I have 
to code my own time-stepping routine?

Thanks,
Max

Re: [petsc-users] TS and petscFE -unable to access u_t

[Maximilian Hartig]

Thank you,

I used these routines to setup a CN type TS. The problem I face now is that I 
seem to be unable to access the temporal derivatives u_t[..] for the definition 
of the residuals. I get a segmentation violation inside the  
PetscFEIntegrateResidual routine whenever I try to. I have attached the error 
message below.
Also, I am wondering whether it is possible to update the neumann boundary 
condition of the FE object for each timestep. This would be useful for coupling 
purposes.

Thank you,

Max


[0]PETSC ERROR: 

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably 
memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see 
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to 
find memory corruption errors
[0]PETSC ERROR: likely location of problem given in stack below
[0]PETSC ERROR: -  Stack Frames 

[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
[0]PETSC ERROR:   INSTEAD the line number of the start of the function
[0]PETSC ERROR:   is given.
[0]PETSC ERROR: [0] PetscFEIntegrateResidual_Basic line 3503 
/Users/maxhartig/PETSc/src/dm/dt/interface/dtfe.c
[0]PETSC ERROR: [0] PetscFEIntegrateResidual line 5753 
/Users/maxhartig/PETSc/src/dm/dt/interface/dtfe.c
[0]PETSC ERROR: [0] DMPlexComputeResidual_Internal line 1706 
/Users/maxhartig/PETSc/src/snes/utils/dmplexsnes.c
[0]PETSC ERROR: [0] DMPlexSNESComputeResidualFEM line 2152 
/Users/maxhartig/PETSc/src/snes/utils/dmplexsnes.c
[0]PETSC ERROR: [0] SNESComputeFunction_DMLocal line 65 
/Users/maxhartig/PETSc/src/snes/utils/dmlocalsnes.c
[0]PETSC ERROR: [0] SNES user function line 2144 
/Users/maxhartig/PETSc/src/snes/interface/snes.c
[0]PETSC ERROR: [0] SNESComputeFunction line 2129 
/Users/maxhartig/PETSc/src/snes/interface/snes.c
[0]PETSC ERROR: [0] SNESSolve_NEWTONLS line 150 
/Users/maxhartig/PETSc/src/snes/impls/ls/ls.c
[0]PETSC ERROR: [0] SNESSolve line 3961 
/Users/maxhartig/PETSc/src/snes/interface/snes.c
[0]PETSC ERROR: [0] TS_SNESSolve line 188 
/Users/maxhartig/PETSc/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: [0] TSStep_Theta line 206 
/Users/maxhartig/PETSc/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: [0] TSStep line 3700 
/Users/maxhartig/PETSc/src/ts/interface/ts.c
[0]PETSC ERROR: [0] TSSolve line 3921 
/Users/maxhartig/PETSc/src/ts/interface/ts.c
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.7.2, unknown 
[0]PETSC ERROR: Configure options --download-triangle
[0]PETSC ERROR: #1 User provided function() line 0 in  unknown file
application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=59

> On 03 Aug 2016, at 16:44, Matthew Knepley  wrote:
> 
> On Tue, Aug 2, 2016 at 8:22 AM, Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> Hello all,
> 
> I would like to run a transient problem with PetscFE. Example ex11.c seems 
> relevant since it uses the PestcFV context to create boundary conditions and 
> RHS Functions for the TS.
> Is there an easy way to do transient analysis with TS and petscFE or do I 
> have to code my own time-stepping routine?
> 
> You can use
> 
> ierr = DMTSSetBoundaryLocal(adaptedDM,  DMPlexTSComputeBoundary, 
> user);CHKERRQ(ierr);
> ierr = DMTSSetIFunctionLocal(adaptedDM, DMPlexTSComputeIFunctionFEM, 
> user);CHKERRQ(ierr);
> ierr = DMTSSetIJacobianLocal(adaptedDM, DMPlexTSComputeIJacobianFEM, 
> user);CHKERRQ(ierr);
> 
> I have been meaning to write a heat equation example, but I have not finished 
> yet,
> 
>   Thanks,
> 
>  Matt
>  
> Thanks,
> Max
> 
> 
> 
> -- 
> 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] TS and petscFE -unable to access u_t

[Maximilian Hartig]

Gladly. I also send the simple gmsh I am using.

cimply.c
Description: Binary data


testmesh_2D_box_quad.msh
Description: Binary data
Thank you,MaxOn 08 Aug 2016, at 17:23, Matthew Knepley <knep...@gmail.com> wrote:On Mon, Aug 8, 2016 at 10:15 AM, Maximilian Hartig <imilian.har...@gmail.com> wrote:Thank you,I used these routines to setup a CN type TS. The problem I face now is that I seem to be unable to access the temporal derivatives u_t[..] for the definition of the residuals. I get a segmentation violation inside the  PetscFEIntegrateResidual routine whenever I try to. I have attached the error message below.Also, I am wondering whether it is possible to update the neumann boundary condition of the FE object for each timestep. This would be useful for coupling purposes.Could you send you example so I can run it myself?  Thanks,     Matt Thank you,Max[0]PETSC ERROR: [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger[0]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors[0]PETSC ERROR: likely location of problem given in stack below[0]PETSC ERROR: -  Stack Frames [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,[0]PETSC ERROR:       INSTEAD the line number of the start of the function[0]PETSC ERROR:       is given.[0]PETSC ERROR: [0] PetscFEIntegrateResidual_Basic line 3503 /Users/maxhartig/PETSc/src/dm/dt/interface/dtfe.c[0]PETSC ERROR: [0] PetscFEIntegrateResidual line 5753 /Users/maxhartig/PETSc/src/dm/dt/interface/dtfe.c[0]PETSC ERROR: [0] DMPlexComputeResidual_Internal line 1706 /Users/maxhartig/PETSc/src/snes/utils/dmplexsnes.c[0]PETSC ERROR: [0] DMPlexSNESComputeResidualFEM line 2152 /Users/maxhartig/PETSc/src/snes/utils/dmplexsnes.c[0]PETSC ERROR: [0] SNESComputeFunction_DMLocal line 65 /Users/maxhartig/PETSc/src/snes/utils/dmlocalsnes.c[0]PETSC ERROR: [0] SNES user function line 2144 /Users/maxhartig/PETSc/src/snes/interface/snes.c[0]PETSC ERROR: [0] SNESComputeFunction line 2129 /Users/maxhartig/PETSc/src/snes/interface/snes.c[0]PETSC ERROR: [0] SNESSolve_NEWTONLS line 150 /Users/maxhartig/PETSc/src/snes/impls/ls/ls.c[0]PETSC ERROR: [0] SNESSolve line 3961 /Users/maxhartig/PETSc/src/snes/interface/snes.c[0]PETSC ERROR: [0] TS_SNESSolve line 188 /Users/maxhartig/PETSc/src/ts/impls/implicit/theta/theta.c[0]PETSC ERROR: [0] TSStep_Theta line 206 /Users/maxhartig/PETSc/src/ts/impls/implicit/theta/theta.c[0]PETSC ERROR: [0] TSStep line 3700 /Users/maxhartig/PETSc/src/ts/interface/ts.c[0]PETSC ERROR: [0] TSSolve line 3921 /Users/maxhartig/PETSc/src/ts/interface/ts.c[0]PETSC ERROR: - Error Message --[0]PETSC ERROR: Signal received[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.[0]PETSC ERROR: Petsc Release Version 3.7.2, unknown [0]PETSC ERROR: Configure options --download-triangle[0]PETSC ERROR: #1 User provided function() line 0 in  unknown fileapplication called MPI_Abort(MPI_COMM_WORLD, 59) - process 0[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=59On 03 Aug 2016, at 16:44, Matthew Knepley <knep...@gmail.com> wrote:On Tue, Aug 2, 2016 at 8:22 AM, Maximilian Hartig <imilian.har...@gmail.com> wrote:Hello all,I would like to run a transient problem with PetscFE. Example ex11.c seems relevant since it uses the PestcFV context to create boundary conditions and RHS Functions for the TS.Is there an easy way to do transient analysis with TS and petscFE or do I have to code my own time-stepping routine?You can use    ierr = DMTSSetBoundaryLocal(adaptedDM,  DMPlexTSComputeBoundary, user);CHKERRQ(ierr);    ierr = DMTSSetIFunctionLocal(adaptedDM, DMPlexTSComputeIFunctionFEM, user);CHKERRQ(ierr);    ierr = DMTSSetIJacobianLocal(adaptedDM, DMPlexTSComputeIJacobianFEM, user);CHKERRQ(ierr);I have been meaning to write a heat equation example, but I have not finished yet,  Thanks,     Matt Thanks,Max-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.-- Norbert Wiener-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.-- Norbert Wiener

Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Maximilian Hartig]

> On 24. Oct 2018, at 12:49, Matthew Knepley  wrote:
> 
> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell  <mailto:we...@gmx.li>> wrote:
> Hi Max,
> 
> (I'm cc'ing in the petsc-users mailing list which may have more advice, if 
> you are using PETSc you should definitely subscribe!
> 
> > On 24 Oct 2018, at 09:27, Maximilian Hartig  > <mailto:imilian.har...@gmail.com>> wrote:
> > 
> > Hello Lawrence,
> > 
> > sorry to message you out of the blue. My name is Max and I found your post 
> > on GitHub (https://github.com/firedrakeproject/firedrake/issues/1246 
> > <https://github.com/firedrakeproject/firedrake/issues/1246> ) on DMPlex 
> > being able to read periodic gmsh files. I am currently trying to do just 
> > that (creating a periodic DMPlex mesh with gmsh) in the context of my PhD 
> > work. So far I haven’ t found any documentation on the periodic BC’s with 
> > DMPlex and gmsh in the official petsc documentation. 
> > I was wondering whether you’d be so kind as to point me in a general 
> > direction concerning how to achieve this. You seem experienced in using 
> > petsc and I would greatly appreciate your help. 
> 
> 
> I think the answer is "it depends". If you're just using DMPlex directly and 
> all the of the functionality with PetscDS, then I /think/ that reading 
> periodic meshes via gmsh (assuming you're using the appropriate gmsh mesh 
> format [v2]) "just works".
> 
> There are two phases here: topological and geometric. DMPlex represents the 
> periodic topological entity directly. For example,  a circle is just a 
> segment with one end hooked to the other. Vertices are not duplicated, or 
> mapped to each other. This makes topology simple and easy to implement. 
> However, then geometry is more complicated. What Plex does is allow 
> coordinates to be represented by a discontinuous field taking values on 
> cells, in addition to vertices. In our circle example, each cells near the 
> cut will have 2 coordinates, one for each vertex, but they will not agree 
> across the cut. If you define a periodic domain, then Plex can construct this 
> coordinate field automatically using DMPlexLocalize(). These DG coordinates 
> are then used by the integration routines.

Ok, I think I understand the concept. DMPlex reads information about both 
topology and coordinates from the .msh file. Creating a periodic mesh in gmsh 
then should allow DMPlex to identify the periodic boundaries as the “cut” and 
build the mesh topology accordingly. Coordinate information is handled 
separately.
That means, as Lawrence suggested, building periodic meshes in gmsh and reading 
them in to petsc’s DMPlex should indeed “just work”.  (From the user 
perspective). The only extra step is to call DMLocalizeCoordinates() after 
DMPlexReadFromFile(). Sorry to reiterate, I am just trying to make sense of 
this. 
>  
> From my side, the issue is to do with mapping that coordinate field into one 
> that I understand (within Firedrake). You may not have this problem.
> 
> Firedrake uses its own coordinate mapping and integration routines, so they 
> must manage the second part independently. I hope to get change this slightly 
> soon by making the Firedrake representation a DMField, so that it looks the 
> same to Plex.
> 
>   Thanks,
> 
> Matt
>  
> Thanks,
> 
> Lawrence
> 
> 
> -- 
> 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
> 
> https://www.cse.buffalo.edu/~knepley/ <https://www.cse.buffalo.edu/~knepley/>
> To read periodic meshes from GMSH, you need to use the option 
> -dm_plex_gmsh_periodic and DMPlexCreateFromFile

Ahh, thanks. I was missing the option " -dm_plex_gmsh_periodic “. But using 
this option I now generate a segmentation fault error when calling VecView() on 
the solution vector with vtk and hdf5 viewers. Any suggestions?

 
> 
> See  src/dm/impls/plex/examples/tests/ex1.c. An example runs
> 
> $ ./ex1 -filename 
> ${PETSC_DIR}/share/petsc/datafiles/meshes/cube_periodic_bin.msh 
> -dm_plex_gmsh_periodic -dm_view ::ascii_info_detail -interpolate -test_shape
> 
> generating periodic meshes in gmsh may be tricky, Lisandro for sure may 
> advice.

Thanks,
Max

Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Maximilian Hartig]

>> Ahh, thanks. I was missing the option " -dm_plex_gmsh_periodic “. But using 
>> this option I now generate a segmentation fault error when calling VecView() 
>> on the solution vector with vtk and hdf5 viewers. Any suggestions?
>>  

> Small example? VTK is deprecated. HDF5 should work, although it will require 
> you to have proper coordinates I think. We have to
> think about what you mean. If its for a checkpoint, there is no problem, but 
> for viz, those programs do not understand periodicity. Thus I embed it in a 
> higher dimensional space.
> 
>Matt

Building the small example I realised that hdf5 wasn’t working altogether. I’m 
trying to fix this and see if I can run VecView with the periodic DMPlex then.
I planned using this for visualisation. 


> On 25. Oct 2018, at 15:36, Stefano Zampini  wrote:
> 
> Maybe this is a fix
> 
> diff --git a/src/dm/impls/plex/plexvtu.c b/src/dm/impls/plex/plexvtu.c
> index acdea12c2f..1a8bbada6a 100644
> --- a/src/dm/impls/plex/plexvtu.c
> +++ b/src/dm/impls/plex/plexvtu.c
> @@ -465,10 +465,11 @@ PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer 
> viewer)
>  if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
>PetscScalar *xpoint;
>  
> -  ierr = DMPlexPointLocalRead(dm,v,x,&xpoint);CHKERRQ(ierr);
> +  ierr = 
> DMPlexPointLocalRead(dm,closure[v],x,&xpoint);CHKERRQ(ierr);
>y[cnt + off++] = xpoint[i];
>  }
>}
> +  cnt += off;
>ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, 
> &closureSize, &closure);CHKERRQ(ierr);
>  }
>}
> 
> Max, does this fix your problem? If you confirm, I'll fix this in the maint 
> branch

It did indeed fix the issue! Thank you. I’ve had problems with hdf5 before and 
switched to vtk. It’s good news I can continue to use it.
> 
> If I run the below command line with the patch and with snes tutorials ex12 I 
> get the nice picture attached
> $ ./ex12 -quiet -run_type test -interpolate 1 -bc_type dirichlet 
> -petscspace_degree 1 -vec_view vtk:test.vtu:vtk_vtu -f 
> ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh 
> -dm_plex_gmsh_periodic -x_periodicity periodic -y_periodicity periodic 
> -dm_refine 4
> 
> Il giorno gio 25 ott 2018 alle ore 15:11 Stefano Zampini 
> mailto:stefano.zamp...@gmail.com>> ha scritto:
> Matt,
> 
> you can reproduce it via
> 
> $ valgrind ./ex12 -quiet -run_type test -interpolate 1 -bc_type dirichlet 
> -petscspace_degree 1 -vec_view vtk:test.vtu:vtk_vtu -f 
> ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh 
> -dm_plex_gmsh_periodic
> 
> Long time ago I added support for viewing meshes with periodic vertices in 
> the VTK_VTU viewer, but I did not fix the part that writes fields
> 
> 
> Il giorno mer 24 ott 2018 alle ore 21:04 Matthew Knepley  <mailto:knep...@gmail.com>> ha scritto:
> On Wed, Oct 24, 2018 at 11:36 AM Maximilian Hartig  <mailto:imilian.har...@gmail.com>> wrote:
> 
> 
>> On 24. Oct 2018, at 12:49, Matthew Knepley > <mailto:knep...@gmail.com>> wrote:
>> 
>> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell > <mailto:we...@gmx.li>> wrote:
>> Hi Max,
>> 
>> (I'm cc'ing in the petsc-users mailing list which may have more advice, if 
>> you are using PETSc you should definitely subscribe!
>> 
>> > On 24 Oct 2018, at 09:27, Maximilian Hartig > > <mailto:imilian.har...@gmail.com>> wrote:
>> > 
>> > Hello Lawrence,
>> > 
>> > sorry to message you out of the blue. My name is Max and I found your post 
>> > on GitHub (https://github.com/firedrakeproject/firedrake/issues/1246 
>> > <https://github.com/firedrakeproject/firedrake/issues/1246> ) on DMPlex 
>> > being able to read periodic gmsh files. I am currently trying to do just 
>> > that (creating a periodic DMPlex mesh with gmsh) in the context of my PhD 
>> > work. So far I haven’ t found any documentation on the periodic BC’s with 
>> > DMPlex and gmsh in the official petsc documentation. 
>> > I was wondering whether you’d be so kind as to point me in a general 
>> > direction concerning how to achieve this. You seem experienced in using 
>> > petsc and I would greatly appreciate your help. 
>> 
>> 
>> I think the answer is "it depends". If you're just using DMPlex directly and 
>> all the of the functionality with PetscDS, then I /think/ that reading 
>> periodic meshes via gmsh (assuming you