libmesh-users — General information for libmesh users

Re: [Libmesh-users] 答复: 答复: Not decreasing error bound

[David K. <dav...@ak...>]

On Sun, Jan 28, 2018 at 11:24 PM, <cau...@gm...> wrote:

> OK, now that the choice of solver does not pose a problem, I consider
> using SCM to calculate the lower bound of the coercivity constant (αLB)
> instead of setting it 1 since αLB may vary a lot with the parameters and
> the selection of basis will be affected if we simply set αLB equal to 1.
>
>
>
> Hence, I refer to reduced_basis_ex2 where the SCM is implemented.
> However, I meet the error below when running the example.
>
>
>
> “Assembling affine operator 1 of 3
> Assembling affine operator 2 of 3
> Assembling affine operator 3 of 3
> Assembling affine vector 1 of 1
> Assembling output vector, (1,1) of (4,1)
> Assembling output vector, (2,1) of (4,1)
> Assembling output vector, (3,1) of (4,1)
> Assembling output vector, (4,1) of (4,1)
>
> B_min(0) = 1.20285e-33
> Eigen solver for computing B_max did not converge”
>
>
>
> Where does the problem come from and how to address it? I would appreciate
> your suggestion.
>


Are you using "-eps_type lapack"? If not, I would try that.

David








*发件人**:* David Knezevic [mailto:dav...@ak...]
*发送时间:* 2018年1月29日 10:45
*收件人:* Gauvain Wu <cau...@gm...>; libmesh-users <
lib...@li...>
*主题:* Re: 答复: [Libmesh-users] Not decreasing error bound



(Note: Please reply-all so that the answers are recorded on the
libMesh-users list for future reference.)





On Sun, Jan 28, 2018 at 9:33 PM, <cau...@gm...> wrote:

Hi David,



I solve the system with the command below:



-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps



It means that the system is solved by KSP with LU as preconditioning.



OK, that's good. I guess before you were using an iterative solver which
just wasn't converging fully. Using LU is a good approach here. I would
recommend MUMPS with LU for the reduced basis training unless your problem
is too large, in which case you have to switch to an iterative solver.





Should I precise the asymmetric solver in my command? How to write it
correctly? Thanks for your help.



No, what you have done is fine, since you asked for LU which is appropriate
for non-symmetric matrices. (In other problems if your system is symmetric
you can use MUMPS's cholesky solver instead, via "-pc_type cholesky".)



David







*发件人**:* David Knezevic [mailto:dav...@ak...]
*发送时间:* 2018年1月29日 10:05
*收件人:* Gauvain Wu <cau...@gm...>
*抄送:* libmesh-users <lib...@li...>
*主题:* Re: [Libmesh-users] Not decreasing error bound



Hello,



The convergence behavior that you describe is typical of reduced basis
convergence: It will plateau after an error reduction of about six orders
of magnitude or so. So it sounds like the convergence is working fine in
the sense that you got a reduction from 1.3569e7 to 41. When you get the
message "Exiting greedy because the same parameters were selected twice"
that is another indication that the greedy algorithm has plateaued.



I do not know why the RB solution and FE solution did not match well at the
end, though --- that of course indicates that something is wrong. One
thought; Did you make sure to use an asymmetric solver, since
thermoelasticity is not symmetric?



David





On Sat, Jan 27, 2018 at 4:13 AM, <cau...@gm...> wrote:

Hi all,



I made a thermoelasticity model based on the cantilever example,
reduced_basis_ex5, by adding a new temperature variable. At the beginning of
the basis training procedure, the maximum error bound drops sharply from
1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5. After
that, although the basis dimension keeps growing, the error bound stops
decreasing and stays at a certain number. The relative training tolerance is
set at 1.e-7 and the mesh is a T-shaped pipe.



---- Basis dimension: 5 ----

Performing RB solves on training set

Maximum error bound is 2.42578



Performing truth solve at parameter:

h: 1.055972e+01

h_Tinf: 2.472563e+02

heat_flux: 4.261782e+01



---- Basis dimension: 6 ----

Performing RB solves on training set

Maximum error bound is 2.43818



Performing truth solve at parameter:

h: 1.151397e+01

h_Tinf: 2.473108e+02

heat_flux: 4.481571e+01



---- Basis dimension: 7 ----

Performing RB solves on training set

Maximum error bound is 2.44673



Exiting greedy because the same parameters were selected twice



The RB result obtained from this basis differs a lot from the FEM result. I
searched archives of the mailing list and found that this phenomenon might
result from an overly low training tolerance. However, the initial error
bound being nearly e+07, if I select a less strict tolerance, I will end up
having an unsatisfying error and probably a worse result. Could you please
suggest me some advice? I would be grateful for your response.



Best regards,

Gauvain

------------------------------------------------------------
------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Lib...@li...
https://lists.sourceforge.net/lists/listinfo/libmesh-users

[Libmesh-users] 答复: 答复: Not decreasing error bound

[<cau...@gm...>]

OK, now that the choice of solver does not pose a problem, I consider using SCM to calculate the lower bound of the coercivity constant (αLB) instead of setting it 1 since αLB may vary a lot with the parameters and the selection of basis will be affected if we simply set αLB equal to 1. 

 

Hence, I refer to reduced_basis_ex2 where the SCM is implemented. However, I meet the error below when running the example. 

 

“Assembling affine operator 1 of 3
Assembling affine operator 2 of 3
Assembling affine operator 3 of 3
Assembling affine vector 1 of 1
Assembling output vector, (1,1) of (4,1)
Assembling output vector, (2,1) of (4,1)
Assembling output vector, (3,1) of (4,1)
Assembling output vector, (4,1) of (4,1)

B_min(0) = 1.20285e-33
Eigen solver for computing B_max did not converge”

 

Where does the problem come from and how to address it? I would appreciate your suggestion.

 

Gauvain

 

发件人: David Knezevic [mailto:dav...@ak...] 
发送时间: 2018年1月29日 10:45
收件人: Gauvain Wu <cau...@gm...>; libmesh-users <lib...@li...>
主题: Re: 答复: [Libmesh-users] Not decreasing error bound

 

(Note: Please reply-all so that the answers are recorded on the libMesh-users list for future reference.)

 

 

On Sun, Jan 28, 2018 at 9:33 PM, <cau...@gm... <mailto:cau...@gm...> > wrote:

Hi David,

 

I solve the system with the command below:

 

-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps

 

It means that the system is solved by KSP with LU as preconditioning.

 

OK, that's good. I guess before you were using an iterative solver which just wasn't converging fully. Using LU is a good approach here. I would recommend MUMPS with LU for the reduced basis training unless your problem is too large, in which case you have to switch to an iterative solver.

 

 

Should I precise the asymmetric solver in my command? How to write it correctly? Thanks for your help.

 

No, what you have done is fine, since you asked for LU which is appropriate for non-symmetric matrices. (In other problems if your system is symmetric you can use MUMPS's cholesky solver instead, via "-pc_type cholesky".)

 

David

 

 

 

发件人: David Knezevic [mailto:dav...@ak... <mailto:dav...@ak...> ] 
发送时间: 2018年1月29日 10:05
收件人: Gauvain Wu <cau...@gm... <mailto:cau...@gm...> >
抄送: libmesh-users <lib...@li... <mailto:lib...@li...> >
主题: Re: [Libmesh-users] Not decreasing error bound

 

Hello,

 

The convergence behavior that you describe is typical of reduced basis convergence: It will plateau after an error reduction of about six orders of magnitude or so. So it sounds like the convergence is working fine in the sense that you got a reduction from 1.3569e7 to 41. When you get the message "Exiting greedy because the same parameters were selected twice" that is another indication that the greedy algorithm has plateaued.

 

I do not know why the RB solution and FE solution did not match well at the end, though --- that of course indicates that something is wrong. One thought; Did you make sure to use an asymmetric solver, since thermoelasticity is not symmetric?

 

David

 

 

On Sat, Jan 27, 2018 at 4:13 AM, <cau...@gm... <mailto:cau...@gm...> > wrote:

Hi all,



I made a thermoelasticity model based on the cantilever example,
reduced_basis_ex5, by adding a new temperature variable. At the beginning of
the basis training procedure, the maximum error bound drops sharply from
1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5. After
that, although the basis dimension keeps growing, the error bound stops
decreasing and stays at a certain number. The relative training tolerance is
set at 1.e-7 and the mesh is a T-shaped pipe.



---- Basis dimension: 5 ----

Performing RB solves on training set

Maximum error bound is 2.42578



Performing truth solve at parameter:

h: 1.055972e+01

h_Tinf: 2.472563e+02

heat_flux: 4.261782e+01



---- Basis dimension: 6 ----

Performing RB solves on training set

Maximum error bound is 2.43818



Performing truth solve at parameter:

h: 1.151397e+01

h_Tinf: 2.473108e+02

heat_flux: 4.481571e+01



---- Basis dimension: 7 ----

Performing RB solves on training set

Maximum error bound is 2.44673



Exiting greedy because the same parameters were selected twice



The RB result obtained from this basis differs a lot from the FEM result. I
searched archives of the mailing list and found that this phenomenon might
result from an overly low training tolerance. However, the initial error
bound being nearly e+07, if I select a less strict tolerance, I will end up
having an unsatisfying error and probably a worse result. Could you please
suggest me some advice? I would be grateful for your response.



Best regards,

Gauvain

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Lib...@li... <mailto:Lib...@li...> 
https://lists.sourceforge.net/lists/listinfo/libmesh-users

 

 

Re: [Libmesh-users] 答复: Not decreasing error bound

[David K. <dav...@ak...>]

On Sun, Jan 28, 2018 at 9:57 PM, <cau...@gm...> wrote:

> Sorry for habitually clicking “reply” instead of “reply all” and I resend
> it now.
>
>
>
> Hi David,
>
>
>
> I solve the system with the command below:
>
>
>
> -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps
>
>
>
> It means that the system is solved by KSP with LU as preconditioning.
> Should I precise the asymmetric solver in my command? How to write it
> correctly? Thanks for your help.
>

No, what you did is fine since LU is appropriate for non-symmetric
problems. (See my previous email for a bit more detail.)

David


>
>
>
> *发件人**:* David Knezevic [mailto:dav...@ak...]
> *发送时间:* 2018年1月29日 10:05
> *收件人:* Gauvain Wu <cau...@gm...>
> *抄送:* libmesh-users <lib...@li...>
> *主题:* Re: [Libmesh-users] Not decreasing error bound
>
>
>
> Hello,
>
>
>
> The convergence behavior that you describe is typical of reduced basis
> convergence: It will plateau after an error reduction of about six orders
> of magnitude or so. So it sounds like the convergence is working fine in
> the sense that you got a reduction from 1.3569e7 to 41. When you get the
> message "Exiting greedy because the same parameters were selected twice"
> that is another indication that the greedy algorithm has plateaued.
>
>
>
> I do not know why the RB solution and FE solution did not match well at
> the end, though --- that of course indicates that something is wrong. One
> thought; Did you make sure to use an asymmetric solver, since
> thermoelasticity is not symmetric?
>
>
>
> David
>
>
>
>
>
> On Sat, Jan 27, 2018 at 4:13 AM, <cau...@gm...> wrote:
>
> Hi all,
>
>
>
> I made a thermoelasticity model based on the cantilever example,
> reduced_basis_ex5, by adding a new temperature variable. At the beginning
> of
> the basis training procedure, the maximum error bound drops sharply from
> 1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5.
> After
> that, although the basis dimension keeps growing, the error bound stops
> decreasing and stays at a certain number. The relative training tolerance
> is
> set at 1.e-7 and the mesh is a T-shaped pipe.
>
>
>
> ---- Basis dimension: 5 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 2.42578
>
>
>
> Performing truth solve at parameter:
>
> h: 1.055972e+01
>
> h_Tinf: 2.472563e+02
>
> heat_flux: 4.261782e+01
>
>
>
> ---- Basis dimension: 6 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 2.43818
>
>
>
> Performing truth solve at parameter:
>
> h: 1.151397e+01
>
> h_Tinf: 2.473108e+02
>
> heat_flux: 4.481571e+01
>
>
>
> ---- Basis dimension: 7 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 2.44673
>
>
>
> Exiting greedy because the same parameters were selected twice
>
>
>
> The RB result obtained from this basis differs a lot from the FEM result. I
> searched archives of the mailing list and found that this phenomenon might
> result from an overly low training tolerance. However, the initial error
> bound being nearly e+07, if I select a less strict tolerance, I will end up
> having an unsatisfying error and probably a worse result. Could you please
> suggest me some advice? I would be grateful for your response.
>
>
>
> Best regards,
>
> Gauvain
>
> ------------------------------------------------------------
> ------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
>
>
>
>

[Libmesh-users] 答复: Not decreasing error bound

[<cau...@gm...>]

Sorry for habitually clicking “reply” instead of “reply all” and I resend it now. 

 

Hi David,

 

I solve the system with the command below:

 

-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps

 

It means that the system is solved by KSP with LU as preconditioning. Should I precise the asymmetric solver in my command? How to write it correctly? Thanks for your help. 

 

Regards,

Gauvain

 

 

发件人: David Knezevic [mailto:dav...@ak...] 
发送时间: 2018年1月29日 10:05
收件人: Gauvain Wu <cau...@gm...>
抄送: libmesh-users <lib...@li...>
主题: Re: [Libmesh-users] Not decreasing error bound

 

Hello,

 

The convergence behavior that you describe is typical of reduced basis convergence: It will plateau after an error reduction of about six orders of magnitude or so. So it sounds like the convergence is working fine in the sense that you got a reduction from 1.3569e7 to 41. When you get the message "Exiting greedy because the same parameters were selected twice" that is another indication that the greedy algorithm has plateaued.

 

I do not know why the RB solution and FE solution did not match well at the end, though --- that of course indicates that something is wrong. One thought; Did you make sure to use an asymmetric solver, since thermoelasticity is not symmetric?

 

David

 

 

On Sat, Jan 27, 2018 at 4:13 AM, <cau...@gm... <mailto:cau...@gm...> > wrote:

Hi all,



I made a thermoelasticity model based on the cantilever example,
reduced_basis_ex5, by adding a new temperature variable. At the beginning of
the basis training procedure, the maximum error bound drops sharply from
1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5. After
that, although the basis dimension keeps growing, the error bound stops
decreasing and stays at a certain number. The relative training tolerance is
set at 1.e-7 and the mesh is a T-shaped pipe.



---- Basis dimension: 5 ----

Performing RB solves on training set

Maximum error bound is 2.42578



Performing truth solve at parameter:

h: 1.055972e+01

h_Tinf: 2.472563e+02

heat_flux: 4.261782e+01



---- Basis dimension: 6 ----

Performing RB solves on training set

Maximum error bound is 2.43818



Performing truth solve at parameter:

h: 1.151397e+01

h_Tinf: 2.473108e+02

heat_flux: 4.481571e+01



---- Basis dimension: 7 ----

Performing RB solves on training set

Maximum error bound is 2.44673



Exiting greedy because the same parameters were selected twice



The RB result obtained from this basis differs a lot from the FEM result. I
searched archives of the mailing list and found that this phenomenon might
result from an overly low training tolerance. However, the initial error
bound being nearly e+07, if I select a less strict tolerance, I will end up
having an unsatisfying error and probably a worse result. Could you please
suggest me some advice? I would be grateful for your response.



Best regards,

Gauvain

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Lib...@li... <mailto:Lib...@li...> 
https://lists.sourceforge.net/lists/listinfo/libmesh-users

 

Re: [Libmesh-users] 答复: Not decreasing error bound

[David K. <dav...@ak...>]

(Note: Please reply-all so that the answers are recorded on the
libMesh-users list for future reference.)


On Sun, Jan 28, 2018 at 9:33 PM, <cau...@gm...> wrote:

> Hi David,
>
>
>
> I solve the system with the command below:
>
>
>
> -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps
>
>
>
> It means that the system is solved by KSP with LU as preconditioning.
>

OK, that's good. I guess before you were using an iterative solver which
just wasn't converging fully. Using LU is a good approach here. I would
recommend MUMPS with LU for the reduced basis training unless your problem
is too large, in which case you have to switch to an iterative solver.



> Should I precise the asymmetric solver in my command? How to write it
> correctly? Thanks for your help.
>

No, what you have done is fine, since you asked for LU which is appropriate
for non-symmetric matrices. (In other problems if your system is symmetric
you can use MUMPS's cholesky solver instead, via "-pc_type cholesky".)

David



*发件人**:* David Knezevic [mailto:dav...@ak...]
*发送时间:* 2018年1月29日 10:05
*收件人:* Gauvain Wu <cau...@gm...>
*抄送:* libmesh-users <lib...@li...>
*主题:* Re: [Libmesh-users] Not decreasing error bound



Hello,



The convergence behavior that you describe is typical of reduced basis
convergence: It will plateau after an error reduction of about six orders
of magnitude or so. So it sounds like the convergence is working fine in
the sense that you got a reduction from 1.3569e7 to 41. When you get the
message "Exiting greedy because the same parameters were selected twice"
that is another indication that the greedy algorithm has plateaued.



I do not know why the RB solution and FE solution did not match well at the
end, though --- that of course indicates that something is wrong. One
thought; Did you make sure to use an asymmetric solver, since
thermoelasticity is not symmetric?



David





On Sat, Jan 27, 2018 at 4:13 AM, <cau...@gm...> wrote:

Hi all,



I made a thermoelasticity model based on the cantilever example,
reduced_basis_ex5, by adding a new temperature variable. At the beginning of
the basis training procedure, the maximum error bound drops sharply from
1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5. After
that, although the basis dimension keeps growing, the error bound stops
decreasing and stays at a certain number. The relative training tolerance is
set at 1.e-7 and the mesh is a T-shaped pipe.



---- Basis dimension: 5 ----

Performing RB solves on training set

Maximum error bound is 2.42578



Performing truth solve at parameter:

h: 1.055972e+01

h_Tinf: 2.472563e+02

heat_flux: 4.261782e+01



---- Basis dimension: 6 ----

Performing RB solves on training set

Maximum error bound is 2.43818



Performing truth solve at parameter:

h: 1.151397e+01

h_Tinf: 2.473108e+02

heat_flux: 4.481571e+01



---- Basis dimension: 7 ----

Performing RB solves on training set

Maximum error bound is 2.44673



Exiting greedy because the same parameters were selected twice



The RB result obtained from this basis differs a lot from the FEM result. I
searched archives of the mailing list and found that this phenomenon might
result from an overly low training tolerance. However, the initial error
bound being nearly e+07, if I select a less strict tolerance, I will end up
having an unsatisfying error and probably a worse result. Could you please
suggest me some advice? I would be grateful for your response.



Best regards,

Gauvain

------------------------------------------------------------
------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Lib...@li...
https://lists.sourceforge.net/lists/listinfo/libmesh-users

Re: [Libmesh-users] Not decreasing error bound

[David K. <dav...@ak...>]

Hello,

The convergence behavior that you describe is typical of reduced basis
convergence: It will plateau after an error reduction of about six orders
of magnitude or so. So it sounds like the convergence is working fine in
the sense that you got a reduction from 1.3569e7 to 41. When you get the
message "Exiting greedy because the same parameters were selected twice"
that is another indication that the greedy algorithm has plateaued.

I do not know why the RB solution and FE solution did not match well at the
end, though --- that of course indicates that something is wrong. One
thought; Did you make sure to use an asymmetric solver, since
thermoelasticity is not symmetric?

David


On Sat, Jan 27, 2018 at 4:13 AM, <cau...@gm...> wrote:

> Hi all,
>
>
>
> I made a thermoelasticity model based on the cantilever example,
> reduced_basis_ex5, by adding a new temperature variable. At the beginning
> of
> the basis training procedure, the maximum error bound drops sharply from
> 1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5.
> After
> that, although the basis dimension keeps growing, the error bound stops
> decreasing and stays at a certain number. The relative training tolerance
> is
> set at 1.e-7 and the mesh is a T-shaped pipe.
>
>
>
> ---- Basis dimension: 5 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 2.42578
>
>
>
> Performing truth solve at parameter:
>
> h: 1.055972e+01
>
> h_Tinf: 2.472563e+02
>
> heat_flux: 4.261782e+01
>
>
>
> ---- Basis dimension: 6 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 2.43818
>
>
>
> Performing truth solve at parameter:
>
> h: 1.151397e+01
>
> h_Tinf: 2.473108e+02
>
> heat_flux: 4.481571e+01
>
>
>
> ---- Basis dimension: 7 ----
>
> Performing RB solves on training set
>
> Maximum error bound is 2.44673
>
>
>
> Exiting greedy because the same parameters were selected twice
>
>
>
> The RB result obtained from this basis differs a lot from the FEM result. I
> searched archives of the mailing list and found that this phenomenon might
> result from an overly low training tolerance. However, the initial error
> bound being nearly e+07, if I select a less strict tolerance, I will end up
> having an unsatisfying error and probably a worse result. Could you please
> suggest me some advice? I would be grateful for your response.
>
>
>
> Best regards,
>
> Gauvain
>
> ------------------------------------------------------------
> ------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

[Libmesh-users] Not decreasing error bound

[<cau...@gm...>]

Hi all,

 

I made a thermoelasticity model based on the cantilever example,
reduced_basis_ex5, by adding a new temperature variable. At the beginning of
the basis training procedure, the maximum error bound drops sharply from
1.35694e+07 to 41 as the dimension of the basis increases from 0 to 5. After
that, although the basis dimension keeps growing, the error bound stops
decreasing and stays at a certain number. The relative training tolerance is
set at 1.e-7 and the mesh is a T-shaped pipe. 

 

---- Basis dimension: 5 ----

Performing RB solves on training set

Maximum error bound is 2.42578

 

Performing truth solve at parameter:

h: 1.055972e+01

h_Tinf: 2.472563e+02

heat_flux: 4.261782e+01

 

---- Basis dimension: 6 ----

Performing RB solves on training set

Maximum error bound is 2.43818

 

Performing truth solve at parameter:

h: 1.151397e+01

h_Tinf: 2.473108e+02

heat_flux: 4.481571e+01

 

---- Basis dimension: 7 ----

Performing RB solves on training set

Maximum error bound is 2.44673

 

Exiting greedy because the same parameters were selected twice

 

The RB result obtained from this basis differs a lot from the FEM result. I
searched archives of the mailing list and found that this phenomenon might
result from an overly low training tolerance. However, the initial error
bound being nearly e+07, if I select a less strict tolerance, I will end up
having an unsatisfying error and probably a worse result. Could you please
suggest me some advice? I would be grateful for your response. 

 

Best regards,

Gauvain

Re: [Libmesh-users] Requirements for adding dof_constraint_rows in parallel

[David K. <dav...@ak...>]

On Fri, Jan 26, 2018 at 11:06 AM, Roy Stogner <roy...@ic...>
wrote:

>
> On Thu, 25 Jan 2018, David Knezevic wrote:
>
> Thanks for that info. I'm only using ReplicatedMesh for now. Regarding
>> this:
>>
>>       IIRC: on ReplicatedMesh we try to avoid communication as much
>> as possible, so you have to add user constraints identically on
>>       every processor.  There might be ways to get around that (add
>> local constraints only, then re-add after any repartitioning?) but
>>       I definitely wouldn't suggest trying to depend on them.
>>
>> hmm, the reason I asked about this was because I had a case that hit
>> a "cyclic loop" assertion in parallel when I added the constraints
>> identically on every processor,
>>
>> and then when I changed it so that I only added the constraint on
>> the processor that owned the constrained dof, it worked fine. I
>> assumed that I must have been violating an assumption by adding the
>> constraints identically on all processors, but your comment would
>> indicate that isn't the case...
>>
>
> Looking back at dof_map_constraints.C, it looks like my comment was
> wrong - I did parallelize the constraint calculations, even on
> ReplicatedMesh.  We're looping only over local elements when we
> compute hanging node and Dirichlet constraints.  That means we *must*
> communicate (e.g. to handle cases where a local constraint only gets
> computed on a ghost element) later.



OK, thanks for confirming that. Yes, I see that we only loop over local
elements in DofMap::create_dof_constraints.



> I'd like to understand this better. Can you point me to the relevant
>> code in the library? I've been looking through
>> allgather_recursive_constraints and similar, and I figured they send
>> and receive constraint rows as needed, which is why I assumed we
>> didn't need to add constraint rows identically on all processors.
>>
>
> The first section of allgather_recursive_constraints() is what's
> relevant here - it does the initial push of constraints that we've
> calculated to their owners.
>
> And based on that code, it looks like the rule is "we don't need to
> add constraint rows on all processors, but if we *do* then we need to
> add them identically".
>

OK, that makes sense. Yeah, I see where we have the comment "If we don't
already have a constraint for this dof, add the one we were sent"



> It might be worth adding some error checking there to make sure that
> happens.  Right now if a processor has a constrained row, and gets
> pushed a constraint for that row, it assumes the remote constraint is
> equivalent to the local and it ignores it.  We should probably
> libmesh_assert equivalency (within TOLERANCE*TOLERANCE ought to be
> safe; these are nondimensional numbers) instead.
>

Agreed, that seems like a good idea. I guess that check could be added near
the comment that I pasted above.


Thanks for your help! That clarifies things a lot!

Regards,
David

Re: [Libmesh-users] Requirements for adding dof_constraint_rows in parallel

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, David Knezevic wrote:

> Thanks for that info. I'm only using ReplicatedMesh for now. Regarding this:
>
>       IIRC: on ReplicatedMesh we try to avoid communication as much as possible, so you have to add user constraints identically on
>       every processor.  There might be ways to get around that (add local constraints only, then re-add after any repartitioning?) but
>       I definitely wouldn't suggest trying to depend on them.
> 
> hmm, the reason I asked about this was because I had a case that hit
> a "cyclic loop" assertion in parallel when I added the constraints
> identically on every processor,
>
> and then when I changed it so that I only added the constraint on
> the processor that owned the constrained dof, it worked fine. I
> assumed that I must have been violating an assumption by adding the
> constraints identically on all processors, but your comment would
> indicate that isn't the case...

Looking back at dof_map_constraints.C, it looks like my comment was
wrong - I did parallelize the constraint calculations, even on
ReplicatedMesh.  We're looping only over local elements when we
compute hanging node and Dirichlet constraints.  That means we *must*
communicate (e.g. to handle cases where a local constraint only gets
computed on a ghost element) later.

> I'd like to understand this better. Can you point me to the relevant
> code in the library? I've been looking through
> allgather_recursive_constraints and similar, and I figured they send
> and receive constraint rows as needed, which is why I assumed we
> didn't need to add constraint rows identically on all processors.

The first section of allgather_recursive_constraints() is what's
relevant here - it does the initial push of constraints that we've
calculated to their owners.

And based on that code, it looks like the rule is "we don't need to
add constraint rows on all processors, but if we *do* then we need to
add them identically".

It might be worth adding some error checking there to make sure that
happens.  Right now if a processor has a constrained row, and gets
pushed a constraint for that row, it assumes the remote constraint is
equivalent to the local and it ignores it.  We should probably
libmesh_assert equivalency (within TOLERANCE*TOLERANCE ought to be
safe; these are nondimensional numbers) instead.
---
Roy

Re: [Libmesh-users] Requirements for adding dof_constraint_rows in parallel

[David K. <dav...@ak...>]

Hi Roy,

Thanks for that info. I'm only using ReplicatedMesh for now. Regarding this:

IIRC: on ReplicatedMesh we try to avoid communication as much as possible,
> so you have to add user constraints identically on every processor.  There
> might be ways to get around that (add local constraints only, then re-add
> after any repartitioning?) but I definitely wouldn't suggest trying to
> depend on them.


hmm, the reason I asked about this was because I had a case that hit a
"cyclic loop" assertion in parallel when I added the constraints
identically on every processor, and then when I changed it so that I only
added the constraint on the processor that owned the constrained dof, it
worked fine. I assumed that I must have been violating an assumption by
adding the constraints identically on all processors, but your comment
would indicate that isn't the case...

I'd like to understand this better. Can you point me to the relevant code
in the library? I've been looking through allgather_recursive_constraints
and similar, and I figured they send and receive constraint rows as needed,
which is why I assumed we didn't need to add constraint rows identically on
all processors.

David



On Thu, Jan 25, 2018 at 6:06 PM, Roy Stogner <roy...@ic...>
wrote:

>
> On Thu, 25 Jan 2018, David Knezevic wrote:
>
> Hi all (though Roy I think you're the expert on this one ;) ),
>>
>
> Well, I was the least-amateur when I refactored the code, but an
> *expert* would have also properly documented it, so apparently none of
> those were available.
>
> I know that in DofMap::process_constraints libMesh handles communicating
>> the constraint rows to the processes that need to know about them, but I
>> just want to make sure that I satisfy any requirements that we have on how
>> user-defined constraint rows should be set up in parallel initially.
>>
>
> Good question!
>
> IIRC: on ReplicatedMesh we try to avoid communication as much as
> possible, so you have to add user constraints identically on every
> processor.  There might be ways to get around that (add local
> constraints only, then re-add after any repartitioning?) but I
> definitely wouldn't suggest trying to depend on them.
>
> In particular: are we required to set up the constraint row only on
>> the process that owns the constrained dof index, or can the
>> constraint row be on any process (or all processes)?
>>
>
> On DistributedMesh we set up library-added constraint rows on the
> process that owns the constrained index, and I *think* we have the
> same restriction on user constraints (to avoid the need for an initial
> usually-redundant communication sweep sending out ghost dof
> constraints).
> ---
> Roy
>

Re: [Libmesh-users] Requirements for adding dof_constraint_rows in parallel

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, David Knezevic wrote:

> Hi all (though Roy I think you're the expert on this one ;) ),

Well, I was the least-amateur when I refactored the code, but an
*expert* would have also properly documented it, so apparently none of
those were available.

> I know that in DofMap::process_constraints libMesh handles communicating
> the constraint rows to the processes that need to know about them, but I
> just want to make sure that I satisfy any requirements that we have on how
> user-defined constraint rows should be set up in parallel initially.

Good question!

IIRC: on ReplicatedMesh we try to avoid communication as much as
possible, so you have to add user constraints identically on every
processor.  There might be ways to get around that (add local
constraints only, then re-add after any repartitioning?) but I
definitely wouldn't suggest trying to depend on them.

> In particular: are we required to set up the constraint row only on
> the process that owns the constrained dof index, or can the
> constraint row be on any process (or all processes)?

On DistributedMesh we set up library-added constraint rows on the
process that owns the constrained index, and I *think* we have the
same restriction on user constraints (to avoid the need for an initial
usually-redundant communication sweep sending out ghost dof
constraints).
---
Roy

Re: [Libmesh-users] question on libMeshInit

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, Michael Povolotskyi wrote:

> can libMeshInit accept NULL as constructor arguments?

For argv?  I sure wouldn't bet on it.  We pass that to MPI, Slepc,
PETSc, VTK, GetPot...

I wouldn't even bet on all those codes supporting argc < 1 or an empty
string in argv[0]; IIRC the C/C++ standards require main to be handed
*some* program name even if it hasn't been invoked with any arguments.
---
Roy

[Libmesh-users] question on libMeshInit

[Michael P. <mpo...@pu...>]

Hello,

can libMeshInit accept NULL as constructor arguments?

Michael.

Re: [Libmesh-users] Defining hanging nodes in the XDA format

[Salazar De T. M. <sal...@ll...>]

Thanks, the only alternative I am left with is therefore having a flattened mesh and defining my own hanging node constraints with System::attach_constraint_object(). The only issue is that I need to be careful with the refinement.
-- 
 

On 1/25/18, 12:48 PM, "Roy Stogner" <roy...@ic...> wrote:

    
    On Thu, 25 Jan 2018, Salazar De Troya, Miguel wrote:
    
    > An easier question to make would be: Can I have a parent element
    > with just two children?
    
    Ah!  I'm not sure if this was an easier question to make, but it's a
    *much* easier question to answer: no.
    
    The last time I tried counting, we had literally more than a hundred
    places where libMesh core library code made the core assumption that
    any hierarchical refinement was done isotropically.  IIRC there are
    even many places in the code that make the assumption that the
    refinement pattern within each parent element of a hierarchic mesh
    matches one of the few refinement patterns which libMesh creates
    itself.
    
    It's something we'd all love to change; anisotropic AMR is nice.  But
    I don't see it changing any time soon.
    ---
    Roy
    

[Libmesh-users] Requirements for adding dof_constraint_rows in parallel

[David K. <dav...@ak...>]

Hi all (though Roy I think you're the expert on this one ;) ),

I know that in DofMap::process_constraints libMesh handles communicating
the constraint rows to the processes that need to know about them, but I
just want to make sure that I satisfy any requirements that we have on how
user-defined constraint rows should be set up in parallel initially.

In particular: are we required to set up the constraint row only on the
process that owns the constrained dof index, or can the constraint row be
on any process (or all processes)?

Thanks,
David

Re: [Libmesh-users] Defining hanging nodes in the XDA format

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, Salazar De Troya, Miguel wrote:

> An easier question to make would be: Can I have a parent element
> with just two children?

Ah!  I'm not sure if this was an easier question to make, but it's a
*much* easier question to answer: no.

The last time I tried counting, we had literally more than a hundred
places where libMesh core library code made the core assumption that
any hierarchical refinement was done isotropically.  IIRC there are
even many places in the code that make the assumption that the
refinement pattern within each parent element of a hierarchic mesh
matches one of the few refinement patterns which libMesh creates
itself.

It's something we'd all love to change; anisotropic AMR is nice.  But
I don't see it changing any time soon.
---
Roy

Re: [Libmesh-users] Defining hanging nodes in the XDA format

[Salazar De T. M. <sal...@ll...>]

An easier question to make would be: Can I have a parent element with just two children? For instance, the following mesh has two elements and two children for the second element. I am not able to read the mesh attached. The error is:

Program received signal SIGSEGV, Segmentation fault.
0x00002aaaadab14a4 in libMesh::Elem::parent (this=0x0) at ./include/libmesh/elem.h:2357
2357	  return _elemlinks[0];
Missing separate debuginfos, use: debuginfo-install blas-3.4.2-8.el7.x86_64 cyrus-sasl-lib-2.1.26-21.el7.x86_64 glpk-4.52.1-2.el7.x86_64 hdf5-1.8.12-8.el7.x86_64 keyutils-libs-1.5.8-3.el7.x86_64 krb5-libs-1.15.1-8.2chaos.ch6.x86_64 lapack-3.4.2-8.el7.x86_64 libX11-1.6.5-1.el7.x86_64 libXau-1.0.8-2.1.el7.x86_64 libcom_err-1.42.9-10.el7.x86_64 libcurl-7.29.0-42.el7_4.1.x86_64 libibverbs-13-7.el7.x86_64 libidn-1.28-4.el7.x86_64 libnl3-3.2.28-4.el7.x86_64 libpsm2-10.2.63-2.el7.x86_64 libselinux-2.5-11.el7.x86_64 libssh2-1.4.3-10.el7_2.1.x86_64 libuuid-2.23.2-43.el7_4.2.x86_64 libxcb-1.12-1.el7.x86_64 nspr-4.13.1-1.0.el7_3.x86_64 nss-3.28.4-15.el7_4.x86_64 nss-softokn-freebl-3.28.3-8.el7_4.x86_64 nss-util-3.28.4-3.el7.x86_64 numactl-libs-2.0.9-6.el7_2.x86_64 openldap-2.4.44-5.el7.x86_64 openssl-libs-1.0.2k-8.el7.x86_64
(gdb) bt
#0  0x00002aaaadab14a4 in libMesh::Elem::parent (this=0x0) at ./include/libmesh/elem.h:2357
#1  0x00002aaaae45ae8a in libMesh::MeshTools::libmesh_assert_valid_refinement_tree (mesh=...) at src/mesh/mesh_tools.C:1682
#2  0x00002aaaae1dd7a5 in libMesh::DistributedMesh::delete_remote_elements (this=0x7fffffffbcd8) at src/mesh/distributed_mesh.C:1363
#3  0x00002aaaae27e2c2 in libMesh::MeshBase::prepare_for_use (this=0x7fffffffbcd8, skip_renumber_nodes_and_elements=false, skip_find_neighbors=false) at src/mesh/mesh_base.C:269
#4  0x00002aaaae4e5705 in libMesh::UnstructuredMesh::read (this=0x7fffffffbcd8, name=..., skip_renumber_nodes_and_elements=false, skip_find_neighbors=false) at src/mesh/unstructured_mesh.C:625
#5  0x0000000000424716 in main (argc=1, argv=0x7fffffffc368) at /g/g92/miguel/code/topsm/test/mwe_restart/mwe_restart.C:49

libMesh-1.3.0
4	 # number of elements
10	 # number of nodes
.	 # boundary condition specification file
.	 # subdomain id specification file
n/a	 # processor id specification file
n/a	 # p-level specification file
8	 # type size
0	 # uid size
0	 # pid size
8	 # sid size
0	 # p-level size
8	 # eid size
8	 # side size
8	 # bid size
0	 # subdomain id to name map
2	 # n_elem at level 0, [ type sid (n0 ... nN-1) ]
5 0 0 1 2 3
5 0 1 7 9 2
2	 # n_elem at level 1, [ type sid (n0 ... nN-1) ]
5 1 0 4 7 8 5
5 1 0 5 8 9 6
5.0000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
1.0000000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
1.0000000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
5.0000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
1.2500000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
1.3750000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
1.3750000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
1.3750000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
0	 # presence of unique ids
0	 # sideset id to name map
0	 # number of side boundary conditions
0	 # nodeset id to name map
0	 # number of nodesets
0	 # sideset id to name map
0	 # number of edge boundary conditions
0	 # sideset id to name map
0	 # number of shellface boundary conditions

-- 
 

On 1/24/18, 4:15 PM, "Salazar De Troya, Miguel" <sal...@ll...> wrote:

    Attaching a picture. I drew the eight elements I removed. This way you can see how I generated the initial mesh (with two levels of refinement in the outermost elements and one in the interior).  Without flattening the mesh, I could not remove those children I drew.
    
    Miguel
    --
    
    From: "Paul T. Bauman" <ptb...@gm...>
    Date: Wednesday, January 24, 2018 at 3:45 PM
    To: "Salazar De Troya, Miguel" <sal...@ll...>
    Cc: "lib...@li..." <lib...@li...>
    Subject: Re: [Libmesh-users] Defining hanging nodes in the XDA format
    
    I'm not sure I'm following (a picture is worth a thousand words here), but it sounds like you're trying to create a child of a quad by only splitting in one direction and not both?
    
    On Tue, Jan 23, 2018 at 7:22 PM, Salazar De Troya, Miguel <sal...@ll...<mailto:sal...@ll...>> wrote:
    Hello,
    
    I have the mesh attached below. It is the result of having three quadrangular elements in line, refine the left and right ones, then on their leftmost and rightmost children, do another refinement. After this, I flattened out the mesh and removed the leftmost and rightmost layer of elements. The issue is that I’ve got hanging nodes that libmesh does not recognize. Is there any way to find them or tell libMesh which ones they are?
    
    Is there an alternative way to generate a similar mesh? I was not able to remove the leftmost and rightmost children without flattening the mesh first. It seems as if each parent element must have four elements.
    
    Thanks
    Miguel
    
    libMesh-1.3.0
    13           # number of elements
    26           # number of nodes
    n/a         # boundary condition specification file
    .               # subdomain id specification file
    n/a         # processor id specification file
    n/a         # p-level specification file
    8              # type size
    0              # uid size
    0              # pid size
    8              # sid size
    0              # p-level size
    0              # eid size
    0              # side size
    0              # bid size
    0              # subdomain id to name map
    13           # n_elem at level 0, [ type sid (n0 ... nN-1) ]
    5 0 0 1 2 3
    5 0 4 0 5 6
    5 0 6 5 3 7
    5 0 1 8 9 10
    5 0 10 9 11 2
    5 0 12 4 13 14
    5 0 14 13 6 15
    5 0 15 6 16 17
    5 0 17 16 7 18
    5 0 8 19 20 21
    5 0 21 20 22 9
    5 0 9 22 23 24
    5 0 24 23 25 11
    5.0000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
    1.0000000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
    1.0000000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
    5.0000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
    2.5000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
    5.0000000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
    2.5000000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
    2.5000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
    1.2500000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
    1.0000000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
    2.5000000000000000e+01 1.2500000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+01 1.2500000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
    2.5000000000000000e+01 3.7500000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+01 3.7500000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
    1.3750000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
    1.3750000000000000e+02 1.2500000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+02 1.2500000000000000e+01 0.0000000000000000e+00
    1.3750000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
    1.3750000000000000e+02 3.7500000000000000e+01 0.0000000000000000e+00
    1.2500000000000000e+02 3.7500000000000000e+01 0.0000000000000000e+00
    1.3750000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
    0              # presence of unique ids
    0              # sideset id to name map
    0              # number of side boundary conditions
    0              # nodeset id to name map
    0              # number of nodesets
    0              # sideset id to name map
    0              # number of edge boundary conditions
    0              # sideset id to name map
    0              # number of shellface boundary conditions
    
    
    --
    
    ------------------------------------------------------------------------------
    Check out the vibrant tech community on one of the world's most
    engaging tech sites, Slashdot.org! http://sdm.link/slashdot
    _______________________________________________
    Libmesh-users mailing list
    Lib...@li...<mailto:Lib...@li...>
    https://lists.sourceforge.net/lists/listinfo/libmesh-users
    
    ------------------------------------------------------------------------------
    Check out the vibrant tech community on one of the world's most
    engaging tech sites, Slashdot.org! http://sdm.link/slashdot
    _______________________________________________
    Libmesh-users mailing list
    Lib...@li...
    https://lists.sourceforge.net/lists/listinfo/libmesh-users
    

Re: [Libmesh-users] question about the numeric vector

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, Michael Povolotskyi wrote:

> I added a plain System class, and it worked.

Thanks!  That was fast!  (possibly causality-violatingly so?  I could
swear I saw this message header in my inbox *before* I hit send on my
own previous message...)
---
Roy

> Michael.
>
>
> On 01/25/2018 12:17 PM, Roy Stogner wrote:
>>
>>  On Thu, 25 Jan 2018, Michael Povolotskyi wrote:
>>
>>>  I will add another system and will not solve it.
>>
>>  Make sure it's ExplicitSystem (no, wait, maybe just plain System?
>>  ExplicitSystem does still allocate an rhs, not just a solution).
>>
>>  ImplicitSystem will try to allocate a Jacobian for you whether you
>>  plan to eventually assemble and solve with it or not, and that can get
>>  expensive.
>>  ---
>>  Roy
>> 
>>
>>>  On 01/25/2018 08:25 AM, Roy Stogner wrote:
>>>>
>>>>   On Wed, 24 Jan 2018, Paul T. Bauman wrote:
>>>>
>>>>>   Sorry for slow response. I have nothing better to offer other than 3
>>>>>   vectors, but maybe others will.
>>>>
>>>>   You can add another ExplicitSystem to the mesh, giving it a
>>>>   vector-valued variable corresponding to your solution system's
>>>>   scalar-valued variable.
>>>>
>>>>   This is more flexible (if you have multiple variables in your solution
>>>>   system but only want to add a vector-valued postprocessed field for
>>>>   one of them, for instance) but probably a bit harder to work with.
>>>>   ---
>>>>   Roy
>>>>
>>>>>   On Mon, Jan 22, 2018 at 12:40 PM, Michael Povolotskyi
>>>>>   <mpo...@pu...>
>>>>>   wrote:
>>>>>
>>>>>>   Dear Libmesh developers and users,
>>>>>>
>>>>>>   I need to compute a field on mesh based on a solution.
>>>>>>
>>>>>>   To do so I usually add a vector to a system and fill it with data.
>>>>>>
>>>>>>   This works perfectly for a scalar field.
>>>>>>
>>>>>>   What would you recommend to do for a vector field.
>>>>>>
>>>>>>   Of course, I can add 3 vectors, one for each component of the field.
>>>>>>  Is
>>>>>>   there a better way?
>>>>>>
>>>>>>   Thank you,
>>>>>>
>>>>>>   Michael.
>>>>>> 
>>>>>>
>>>>>>   ------------------------------------------------------------
>>>>>>   ------------------
>>>>>>   Check out the vibrant tech community on one of the world's most
>>>>>>   engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>>   _______________________________________________
>>>>>>   Libmesh-users mailing list
>>>>>>   Lib...@li...
>>>>>>   https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>>>
>>>>>   ------------------------------------------------------------------------------ 
>>>>>
>>>>>   Check out the vibrant tech community on one of the world's most
>>>>>   engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>   _______________________________________________
>>>>>   Libmesh-users mailing list
>>>>>   Lib...@li...
>>>>>   https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>> 
>>>>> 
>>> 
>>> 
>>> 
>
>
>

Re: [Libmesh-users] question about the numeric vector

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, Michael Povolotskyi wrote:

> Can I add just a plain System to the EquationSystem?

That's actually a good question, and if you'd try it and report back
I'd appreciate it.

The answer used to be "no", which is why everybody habitually
recommends adding ExplicitSystem for storing auxilliary data.

But at some point a few years back IIRC I added the ability for us to
instantiate a bare System.  I've only used that in various
postprocessing tools, though.  It should be the "proper" way to create
an auxilliary data system now, but I don't know if anybody's ever
tried it!

So: give it a shot?  Then if it works let us know, and if it doesn't
then file a ticket and use ExplicitSystem instead?

Thanks,
---
Roy

> On 01/25/2018 12:17 PM, Roy Stogner wrote:
>>
>>  On Thu, 25 Jan 2018, Michael Povolotskyi wrote:
>>
>>>  I will add another system and will not solve it.
>>
>>  Make sure it's ExplicitSystem (no, wait, maybe just plain System?
>>  ExplicitSystem does still allocate an rhs, not just a solution).
>>
>>  ImplicitSystem will try to allocate a Jacobian for you whether you
>>  plan to eventually assemble and solve with it or not, and that can get
>>  expensive.
>>  ---
>>  Roy
>> 
>>
>>>  On 01/25/2018 08:25 AM, Roy Stogner wrote:
>>>>
>>>>   On Wed, 24 Jan 2018, Paul T. Bauman wrote:
>>>>
>>>>>   Sorry for slow response. I have nothing better to offer other than 3
>>>>>   vectors, but maybe others will.
>>>>
>>>>   You can add another ExplicitSystem to the mesh, giving it a
>>>>   vector-valued variable corresponding to your solution system's
>>>>   scalar-valued variable.
>>>>
>>>>   This is more flexible (if you have multiple variables in your solution
>>>>   system but only want to add a vector-valued postprocessed field for
>>>>   one of them, for instance) but probably a bit harder to work with.
>>>>   ---
>>>>   Roy
>>>>
>>>>>   On Mon, Jan 22, 2018 at 12:40 PM, Michael Povolotskyi
>>>>>   <mpo...@pu...>
>>>>>   wrote:
>>>>>
>>>>>>   Dear Libmesh developers and users,
>>>>>>
>>>>>>   I need to compute a field on mesh based on a solution.
>>>>>>
>>>>>>   To do so I usually add a vector to a system and fill it with data.
>>>>>>
>>>>>>   This works perfectly for a scalar field.
>>>>>>
>>>>>>   What would you recommend to do for a vector field.
>>>>>>
>>>>>>   Of course, I can add 3 vectors, one for each component of the field.
>>>>>>  Is
>>>>>>   there a better way?
>>>>>>
>>>>>>   Thank you,
>>>>>>
>>>>>>   Michael.
>>>>>> 
>>>>>>
>>>>>>   ------------------------------------------------------------
>>>>>>   ------------------
>>>>>>   Check out the vibrant tech community on one of the world's most
>>>>>>   engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>>   _______________________________________________
>>>>>>   Libmesh-users mailing list
>>>>>>   Lib...@li...
>>>>>>   https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>>>
>>>>>   ------------------------------------------------------------------------------ 
>>>>>
>>>>>   Check out the vibrant tech community on one of the world's most
>>>>>   engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>   _______________________________________________
>>>>>   Libmesh-users mailing list
>>>>>   Lib...@li...
>>>>>   https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>> 
>>>>> 
>>> 
>>> 
>>> 
>
>
>

Re: [Libmesh-users] question about the numeric vector

[Michael P. <mpo...@pu...>]

I added a plain System class, and it worked.

Michael.


On 01/25/2018 12:17 PM, Roy Stogner wrote:
>
> On Thu, 25 Jan 2018, Michael Povolotskyi wrote:
>
>> I will add another system and will not solve it.
>
> Make sure it's ExplicitSystem (no, wait, maybe just plain System?
> ExplicitSystem does still allocate an rhs, not just a solution).
>
> ImplicitSystem will try to allocate a Jacobian for you whether you
> plan to eventually assemble and solve with it or not, and that can get
> expensive.
> ---
> Roy
>
>
>> On 01/25/2018 08:25 AM, Roy Stogner wrote:
>>>
>>>  On Wed, 24 Jan 2018, Paul T. Bauman wrote:
>>>
>>>>  Sorry for slow response. I have nothing better to offer other than 3
>>>>  vectors, but maybe others will.
>>>
>>>  You can add another ExplicitSystem to the mesh, giving it a
>>>  vector-valued variable corresponding to your solution system's
>>>  scalar-valued variable.
>>>
>>>  This is more flexible (if you have multiple variables in your solution
>>>  system but only want to add a vector-valued postprocessed field for
>>>  one of them, for instance) but probably a bit harder to work with.
>>>  ---
>>>  Roy
>>>
>>>>  On Mon, Jan 22, 2018 at 12:40 PM, Michael Povolotskyi
>>>>  <mpo...@pu...>
>>>>  wrote:
>>>>
>>>>>  Dear Libmesh developers and users,
>>>>>
>>>>>  I need to compute a field on mesh based on a solution.
>>>>>
>>>>>  To do so I usually add a vector to a system and fill it with data.
>>>>>
>>>>>  This works perfectly for a scalar field.
>>>>>
>>>>>  What would you recommend to do for a vector field.
>>>>>
>>>>>  Of course, I can add 3 vectors, one for each component of the 
>>>>> field. Is
>>>>>  there a better way?
>>>>>
>>>>>  Thank you,
>>>>>
>>>>>  Michael.
>>>>>
>>>>>
>>>>>  ------------------------------------------------------------
>>>>>  ------------------
>>>>>  Check out the vibrant tech community on one of the world's most
>>>>>  engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>  _______________________________________________
>>>>>  Libmesh-users mailing list
>>>>>  Lib...@li...
>>>>>  https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>>
>>>>  ------------------------------------------------------------------------------ 
>>>>
>>>>  Check out the vibrant tech community on one of the world's most
>>>>  engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>  _______________________________________________
>>>>  Libmesh-users mailing list
>>>>  Lib...@li...
>>>>  https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>
>>>>
>>
>>
>>

Re: [Libmesh-users] question about the numeric vector

[Roy S. <roy...@ic...>]

On Thu, 25 Jan 2018, Michael Povolotskyi wrote:

> I will add another system and will not solve it.

Make sure it's ExplicitSystem (no, wait, maybe just plain System?
ExplicitSystem does still allocate an rhs, not just a solution).

ImplicitSystem will try to allocate a Jacobian for you whether you
plan to eventually assemble and solve with it or not, and that can get
expensive.
---
Roy


> On 01/25/2018 08:25 AM, Roy Stogner wrote:
>>
>>  On Wed, 24 Jan 2018, Paul T. Bauman wrote:
>>
>>>  Sorry for slow response. I have nothing better to offer other than 3
>>>  vectors, but maybe others will.
>>
>>  You can add another ExplicitSystem to the mesh, giving it a
>>  vector-valued variable corresponding to your solution system's
>>  scalar-valued variable.
>>
>>  This is more flexible (if you have multiple variables in your solution
>>  system but only want to add a vector-valued postprocessed field for
>>  one of them, for instance) but probably a bit harder to work with.
>>  ---
>>  Roy
>>
>>>  On Mon, Jan 22, 2018 at 12:40 PM, Michael Povolotskyi
>>>  <mpo...@pu...>
>>>  wrote:
>>>
>>>>  Dear Libmesh developers and users,
>>>>
>>>>  I need to compute a field on mesh based on a solution.
>>>>
>>>>  To do so I usually add a vector to a system and fill it with data.
>>>>
>>>>  This works perfectly for a scalar field.
>>>>
>>>>  What would you recommend to do for a vector field.
>>>>
>>>>  Of course, I can add 3 vectors, one for each component of the field. Is
>>>>  there a better way?
>>>>
>>>>  Thank you,
>>>>
>>>>  Michael.
>>>> 
>>>>
>>>>  ------------------------------------------------------------
>>>>  ------------------
>>>>  Check out the vibrant tech community on one of the world's most
>>>>  engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>  _______________________________________________
>>>>  Libmesh-users mailing list
>>>>  Lib...@li...
>>>>  https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>
>>>  ------------------------------------------------------------------------------ 
>>>
>>>  Check out the vibrant tech community on one of the world's most
>>>  engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>  _______________________________________________
>>>  Libmesh-users mailing list
>>>  Lib...@li...
>>>  https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>> 
>>> 
>
>
>

Re: [Libmesh-users] Defining hanging nodes in the XDA format

[Roy S. <roy...@ic...>]

Could you Cc: the picture to me?  Our list server generally strips
those, and I wasn't able to understand what you were trying to do
based on the text alone.

Thanks,
---
Roy

On Thu, 25 Jan 2018, Salazar De Troya, Miguel wrote:

> Attaching a picture. I drew the eight elements I removed. This way you can see how I generated the initial mesh (with two levels of refinement in the outermost elements and one in the interior).  Without flattening the mesh, I could not remove those children I drew.
>
> Miguel
> --
>
> From: "Paul T. Bauman" <ptb...@gm...>
> Date: Wednesday, January 24, 2018 at 3:45 PM
> To: "Salazar De Troya, Miguel" <sal...@ll...>
> Cc: "lib...@li..." <lib...@li...>
> Subject: Re: [Libmesh-users] Defining hanging nodes in the XDA format
>
> I'm not sure I'm following (a picture is worth a thousand words here), but it sounds like you're trying to create a child of a quad by only splitting in one direction and not both?
>
> On Tue, Jan 23, 2018 at 7:22 PM, Salazar De Troya, Miguel <sal...@ll...<mailto:sal...@ll...>> wrote:
> Hello,
>
> I have the mesh attached below. It is the result of having three quadrangular elements in line, refine the left and right ones, then on their leftmost and rightmost children, do another refinement. After this, I flattened out the mesh and removed the leftmost and rightmost layer of elements. The issue is that I’ve got hanging nodes that libmesh does not recognize. Is there any way to find them or tell libMesh which ones they are?
>
> Is there an alternative way to generate a similar mesh? I was not able to remove the leftmost and rightmost children without flattening the mesh first. It seems as if each parent element must have four elements.
>
> Thanks
> Miguel
>
> libMesh-1.3.0
> 13           # number of elements
> 26           # number of nodes
> n/a         # boundary condition specification file
> .               # subdomain id specification file
> n/a         # processor id specification file
> n/a         # p-level specification file
> 8              # type size
> 0              # uid size
> 0              # pid size
> 8              # sid size
> 0              # p-level size
> 0              # eid size
> 0              # side size
> 0              # bid size
> 0              # subdomain id to name map
> 13           # n_elem at level 0, [ type sid (n0 ... nN-1) ]
> 5 0 0 1 2 3
> 5 0 4 0 5 6
> 5 0 6 5 3 7
> 5 0 1 8 9 10
> 5 0 10 9 11 2
> 5 0 12 4 13 14
> 5 0 14 13 6 15
> 5 0 15 6 16 17
> 5 0 17 16 7 18
> 5 0 8 19 20 21
> 5 0 21 20 22 9
> 5 0 9 22 23 24
> 5 0 24 23 25 11
> 5.0000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
> 1.0000000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
> 1.0000000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
> 5.0000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
> 2.5000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
> 5.0000000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
> 2.5000000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
> 2.5000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
> 1.2500000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
> 1.0000000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
> 2.5000000000000000e+01 1.2500000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+01 1.2500000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
> 2.5000000000000000e+01 3.7500000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+01 3.7500000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
> 1.3750000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
> 1.3750000000000000e+02 1.2500000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+02 1.2500000000000000e+01 0.0000000000000000e+00
> 1.3750000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
> 1.3750000000000000e+02 3.7500000000000000e+01 0.0000000000000000e+00
> 1.2500000000000000e+02 3.7500000000000000e+01 0.0000000000000000e+00
> 1.3750000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
> 0              # presence of unique ids
> 0              # sideset id to name map
> 0              # number of side boundary conditions
> 0              # nodeset id to name map
> 0              # number of nodesets
> 0              # sideset id to name map
> 0              # number of edge boundary conditions
> 0              # sideset id to name map
> 0              # number of shellface boundary conditions
>
>
> --
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...<mailto:Lib...@li...>
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

Re: [Libmesh-users] question about the numeric vector

[Roy S. <roy...@ic...>]

On Wed, 24 Jan 2018, Paul T. Bauman wrote:

> Sorry for slow response. I have nothing better to offer other than 3
> vectors, but maybe others will.

You can add another ExplicitSystem to the mesh, giving it a
vector-valued variable corresponding to your solution system's
scalar-valued variable.

This is more flexible (if you have multiple variables in your solution
system but only want to add a vector-valued postprocessed field for
one of them, for instance) but probably a bit harder to work with.
---
Roy

> On Mon, Jan 22, 2018 at 12:40 PM, Michael Povolotskyi <mpo...@pu...>
> wrote:
>
>> Dear Libmesh developers and users,
>>
>> I need to compute a field on mesh based on a solution.
>>
>> To do so I usually add a vector to a system and fill it with data.
>>
>> This works perfectly for a scalar field.
>>
>> What would you recommend to do for a vector field.
>>
>> Of course, I can add 3 vectors, one for each component of the field. Is
>> there a better way?
>>
>> Thank you,
>>
>> Michael.
>>
>>
>> ------------------------------------------------------------
>> ------------------
>> Check out the vibrant tech community on one of the world's most
>> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>> _______________________________________________
>> Libmesh-users mailing list
>> Lib...@li...
>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
>

Re: [Libmesh-users] Best mesh format for writing and reading in parallel

[Roy S. <roy...@ic...>]

On Wed, 24 Jan 2018, Paul T. Bauman wrote:

> Others can and should chime in, but I believe CheckpointIO will allow this.
> My only uncertainty is on whether is supports AMR meshes or not, but I
> think it does.

It does.

> Beware we're not guaranteeing backward compatibility with it yet, so
> please watch for potential changes in the future if you do use it.

And please let us know about any potential bugs in the present!  The
combination with AMR in particular doesn't have as much test coverage
as it should for such a corner-case-laden feature.

>> Basically, the question in the title. XDA/XDR do not work for distributed
>> meshes in parallel.

Well, technically they "work", just not efficiently enough to be
useful for frequent I/O in problems where you really need a
distributed mesh.
---
Roy

Re: [Libmesh-users] Defining hanging nodes in the XDA format

[Salazar De T. M. <sal...@ll...>]

Attaching a picture. I drew the eight elements I removed. This way you can see how I generated the initial mesh (with two levels of refinement in the outermost elements and one in the interior).  Without flattening the mesh, I could not remove those children I drew.

Miguel
--

From: "Paul T. Bauman" <ptb...@gm...>
Date: Wednesday, January 24, 2018 at 3:45 PM
To: "Salazar De Troya, Miguel" <sal...@ll...>
Cc: "lib...@li..." <lib...@li...>
Subject: Re: [Libmesh-users] Defining hanging nodes in the XDA format

I'm not sure I'm following (a picture is worth a thousand words here), but it sounds like you're trying to create a child of a quad by only splitting in one direction and not both?

On Tue, Jan 23, 2018 at 7:22 PM, Salazar De Troya, Miguel <sal...@ll...<mailto:sal...@ll...>> wrote:
Hello,

I have the mesh attached below. It is the result of having three quadrangular elements in line, refine the left and right ones, then on their leftmost and rightmost children, do another refinement. After this, I flattened out the mesh and removed the leftmost and rightmost layer of elements. The issue is that I’ve got hanging nodes that libmesh does not recognize. Is there any way to find them or tell libMesh which ones they are?

Is there an alternative way to generate a similar mesh? I was not able to remove the leftmost and rightmost children without flattening the mesh first. It seems as if each parent element must have four elements.

Thanks
Miguel

libMesh-1.3.0
13           # number of elements
26           # number of nodes
n/a         # boundary condition specification file
.               # subdomain id specification file
n/a         # processor id specification file
n/a         # p-level specification file
8              # type size
0              # uid size
0              # pid size
8              # sid size
0              # p-level size
0              # eid size
0              # side size
0              # bid size
0              # subdomain id to name map
13           # n_elem at level 0, [ type sid (n0 ... nN-1) ]
5 0 0 1 2 3
5 0 4 0 5 6
5 0 6 5 3 7
5 0 1 8 9 10
5 0 10 9 11 2
5 0 12 4 13 14
5 0 14 13 6 15
5 0 15 6 16 17
5 0 17 16 7 18
5 0 8 19 20 21
5 0 21 20 22 9
5 0 9 22 23 24
5 0 24 23 25 11
5.0000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
1.0000000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
1.0000000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
5.0000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
2.5000000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
5.0000000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
2.5000000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
2.5000000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
1.2500000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
1.0000000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+01 0.0000000000000000e+00 0.0000000000000000e+00
2.5000000000000000e+01 1.2500000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+01 1.2500000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+01 2.5000000000000000e+01 0.0000000000000000e+00
2.5000000000000000e+01 3.7500000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+01 3.7500000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+01 5.0000000000000000e+01 0.0000000000000000e+00
1.3750000000000000e+02 0.0000000000000000e+00 0.0000000000000000e+00
1.3750000000000000e+02 1.2500000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+02 1.2500000000000000e+01 0.0000000000000000e+00
1.3750000000000000e+02 2.5000000000000000e+01 0.0000000000000000e+00
1.3750000000000000e+02 3.7500000000000000e+01 0.0000000000000000e+00
1.2500000000000000e+02 3.7500000000000000e+01 0.0000000000000000e+00
1.3750000000000000e+02 5.0000000000000000e+01 0.0000000000000000e+00
0              # presence of unique ids
0              # sideset id to name map
0              # number of side boundary conditions
0              # nodeset id to name map
0              # number of nodesets
0              # sideset id to name map
0              # number of edge boundary conditions
0              # sideset id to name map
0              # number of shellface boundary conditions


--

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Lib...@li...<mailto:Lib...@li...>
https://lists.sourceforge.net/lists/listinfo/libmesh-users