libmesh-users — General information for libmesh users

Re: [Libmesh-users] problems with running libmesh in parallel

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

Thank you, Roy.
This is exactly what had happened: the test runs till the end with MPICH 
and hangs here with openMPI.
I'm going to inform our admins about this issue.
Michael.


On 09/11/2017 09:49 AM, Roy Stogner wrote:
>
> On Sun, 10 Sep 2017, Michael Povolotskyi wrote:
>
>> I checked that there is no infinite loop.
>> Both ranks pass  this->send (dest_processor_id, sendvec, type1, req, 
>> send_tag);
>> and hang on  this->receive (source_processor_id, recv, type2, recv_tag);
>
>> Both processes are sending zero elements, is this correct?
>
> Not necessarily, but it definitely *could* be correct; there are a
> number of places in the code where one processor has nothing to send
> *except* for the fact that it has nothing to send.
>
>> Could you, please, suggest a simple MPI test to mimic this situation? 
>> Then I can check the MPI implementation.
>
> I haven't checked if this even compiles, but I *believe* what we're
> doing is:
>
> // On 2 procs:
> const unsigned int other = (rank+1)%2;
>
> MPI_Request req;
> MPI_ISend(NULL, 0, MPI_INT, other, 0, comm, &req);
>
> MPI_Status stat;
> MPI_Probe (other, MPI_ANY_TAG, comm, &stat);
> MPI_Status stat2 = stat;
> MPI_Recv(NULL, 0, MPI_INT, other, MPI_ANY_TAG, comm, &stat2);
>
> MPI_Status stat3;
> MPI_Wait(&req, &stat3);
>
> That should do it, I think; if the same behavior happens then you'll
> hang in MPI_Recv.
> ---
> Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

On Sun, 10 Sep 2017, Michael Povolotskyi wrote:

> I checked that there is no infinite loop.
> Both ranks pass  this->send (dest_processor_id, sendvec, type1, req, 
> send_tag);
> and hang on  this->receive (source_processor_id, recv, type2, recv_tag);

> Both processes are sending zero elements, is this correct?

Not necessarily, but it definitely *could* be correct; there are a
number of places in the code where one processor has nothing to send
*except* for the fact that it has nothing to send.

> Could you, please, suggest a simple MPI test to mimic this situation? Then I 
> can check the MPI implementation.

I haven't checked if this even compiles, but I *believe* what we're
doing is:

// On 2 procs:
const unsigned int other = (rank+1)%2;

MPI_Request req;
MPI_ISend(NULL, 0, MPI_INT, other, 0, comm, &req);

MPI_Status stat;
MPI_Probe (other, MPI_ANY_TAG, comm, &stat);
MPI_Status stat2 = stat;
MPI_Recv(NULL, 0, MPI_INT, other, MPI_ANY_TAG, comm, &stat2);

MPI_Status stat3;
MPI_Wait(&req, &stat3);

That should do it, I think; if the same behavior happens then you'll
hang in MPI_Recv.
---
Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

Hello Roy,
I checked that there is no infinite loop.
Both ranks pass  this->send (dest_processor_id, sendvec, type1, req, 
send_tag);
and hang on  this->receive (source_processor_id, recv, type2, recv_tag);

Both processes are sending zero elements, is this correct?

Could you, please, suggest a simple MPI test to mimic this situation? 
Then I can check the MPI implementation.
Thank you,
Michael.


On 09/10/2017 09:07 PM, Roy Stogner wrote:
>
> On Sat, 9 Sep 2017, Michael Povolotskyi wrote:
>
>> (gdb) p  recv_tag
>> $2 = (const libMesh::Parallel::MessageTag &) @0x2ac9d6d3d980: {
>>   static invalid_tag = -2147483648, _tagvalue = -1, _comm = 0x0}
>> (gdb) p  recv_tag
>> $2 = (const libMesh::Parallel::MessageTag &) @0x2b9331482980: {
>>  static invalid_tag = -2147483648, _tagvalue = -1, _comm = 0x0}
>>
>> So, does it look like a logic error or like a problem with openMPI 
>> installation?
>
> In openmpi MPI_ANY_TAG==-1, so they're both set to the right thing,
> not waiting for the a message that will never come.
>
> Could you, after attaching the debugger, try stepping forward in both
> processes, and make sure they're both hanging, not infinite looping?
> ---
> Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

On Sat, 9 Sep 2017, Michael Povolotskyi wrote:

> (gdb) p  recv_tag
> $2 = (const libMesh::Parallel::MessageTag &) @0x2ac9d6d3d980: {
>   static invalid_tag = -2147483648, _tagvalue = -1, _comm = 0x0}
> (gdb) p  recv_tag
> $2 = (const libMesh::Parallel::MessageTag &) @0x2b9331482980: {
>  static invalid_tag = -2147483648, _tagvalue = -1, _comm = 0x0}
>
> So, does it look like a logic error or like a problem with openMPI 
> installation?

In openmpi MPI_ANY_TAG==-1, so they're both set to the right thing,
not waiting for the a message that will never come.

Could you, after attaching the debugger, try stepping forward in both
processes, and make sure they're both hanging, not infinite looping?
---
Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

Hello Roy, following your request, I have checked recv_tag on both MPI 
processes.
Here is what I have got:

process #1 (I do not know the MPI rank, so I just call it as "1")

#9  0x00002ac9d4852a19 in 
libMesh::Parallel::Communicator::send_receive<unsigned int, unsigned 
int> (this=0x7ffc0e5e18b8, dest_processor_id=1, sendvec=...,
     type1=..., source_processor_id=1, recv=..., type2=..., send_tag=...,
     recv_tag=...) at ./include/libmesh/parallel_implementation.h:2788
2788      this->receive (source_processor_id, recv, type2, recv_tag);
(gdb)  p source_processor_id
$1 = 1
(gdb) p  recv_tag
$2 = (const libMesh::Parallel::MessageTag &) @0x2ac9d6d3d980: {
   static invalid_tag = -2147483648, _tagvalue = -1, _comm = 0x0}

process #2

#9  0x00002b932ef97a19 in 
libMesh::Parallel::Communicator::send_receive<unsigned int, unsigned 
int> (this=0x7fff571c8ae8, dest_processor_id=0, sendvec=...,
     type1=..., source_processor_id=0, recv=..., type2=..., send_tag=...,
     recv_tag=...) at ./include/libmesh/parallel_implementation.h:2788
2788      this->receive (source_processor_id, recv, type2, recv_tag);
(gdb) p source_processor_id
$1 = 0
(gdb) p  recv_tag
$2 = (const libMesh::Parallel::MessageTag &) @0x2b9331482980: {
   static invalid_tag = -2147483648, _tagvalue = -1, _comm = 0x0}
(gdb)

So, does it look like a logic error or like a problem with openMPI 
installation?
Thank you,
Michael.

On 09/08/2017 06:20 PM, Roy Stogner wrote:
>
> On Fri, 8 Sep 2017, Michael Povolotskyi wrote:
>
>> still let me share with you the stack trace.
>> It works for me now with mpich, but if I can help to improve libmesh 
>> portability I would be glad to do it.
>
> I originally thought this looked a little like a local linkage mixup,
> in which case there'd be no point in working on it... but I certainly
> don't see anything obviously mixed up in those library names.
>
>> Here is the trace from the 1st MPI process.
>> On the 2nd it looks the same.
>
> They're both stuck in the receive() of a send_receive()??
>
>> If you want me to print something in the debugger, I'll be happy to 
>> do it.
>
> It would be good to verify that they're blocking, not infinite
> looping.  If they're blocking, I suppose you could verify that the
> recv_tag is MPI_ANY_TAG.  If a receive with MPI_ANY_TAG really isn't
> receiving a just-sent message, then I'd just as soon chalk the problem
> up to "screwy local MPI build" and forget about it, so long as MPICH2
> is still working for you.  I wasted too much time recently trying to
> figure out how to workaround a bug in an old MKL version before I
> became sane and just switched to a newer MKL.
>
> However, if there's an infinite loop going on, then we might be able
> to find something strange going on in the libMesh logic, and if you
> could reproduce it with the git master then I'd love more help
> tracking it down.
> ---
> Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

On Fri, 8 Sep 2017, Michael Povolotskyi wrote:

> still let me share with you the stack trace.
> It works for me now with mpich, but if I can help to improve libmesh 
> portability I would be glad to do it.

I originally thought this looked a little like a local linkage mixup,
in which case there'd be no point in working on it... but I certainly
don't see anything obviously mixed up in those library names.

> Here is the trace from the 1st MPI process.
> On the 2nd it looks the same.

They're both stuck in the receive() of a send_receive()??

> If you want me to print something in the debugger, I'll be happy to do it.

It would be good to verify that they're blocking, not infinite
looping.  If they're blocking, I suppose you could verify that the
recv_tag is MPI_ANY_TAG.  If a receive with MPI_ANY_TAG really isn't
receiving a just-sent message, then I'd just as soon chalk the problem
up to "screwy local MPI build" and forget about it, so long as MPICH2
is still working for you.  I wasted too much time recently trying to
figure out how to workaround a bug in an old MKL version before I
became sane and just switched to a newer MKL.

However, if there's an infinite loop going on, then we might be able
to find something strange going on in the libMesh logic, and if you
could reproduce it with the git master then I'd love more help
tracking it down.
---
Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

Dear Roy,
still let me share with you the stack trace.
It works for me now with mpich, but if I can help to improve libmesh 
portability I would be glad to do it.
Here is the trace from the 1st MPI process.
On the 2nd it looks the same.
(gdb) c
Continuing.
^C
Program received signal SIGINT, Interrupt.
mlx5_stall_poll_cq (ibcq=0x2476cb0, ne=8, wc=0x7ffd610c0370) at src/cq.c:341
341     src/cq.c: No such file or directory.
         in src/cq.c
(gdb) up
#1  mlx5_poll_cq (ibcq=0x2476cb0, ne=8, wc=0x7ffd610c0370) at src/cq.c:515
515     in src/cq.c
(gdb)
#2  0x00002af7b5a2ffc4 in ibv_poll_cq (ep=<value optimized out>)
     at /usr/include/infiniband/verbs.h:1245
1245            return cq->context->ops.poll_cq(cq, num_entries, wc);
(gdb)
#3  fi_ibv_rdm_tagged_poll_recv (ep=<value optimized out>)
     at prov/verbs/src/ep_rdm/verbs_tagged_ep_rdm.c:608
608     prov/verbs/src/ep_rdm/verbs_tagged_ep_rdm.c: No such file or 
directory.
         in prov/verbs/src/ep_rdm/verbs_tagged_ep_rdm.c
(gdb)
#4  0x00002af7b5a29d65 in fi_ibv_rdm_tagged_cq_readfrom (
     cq=<value optimized out>, buf=<value optimized out>,
     count=<value optimized out>, src_addr=<value optimized out>)
     at prov/verbs/src/ep_rdm/verbs_cq_ep_rdm.c:82
82      prov/verbs/src/ep_rdm/verbs_cq_ep_rdm.c: No such file or directory.
         in prov/verbs/src/ep_rdm/verbs_cq_ep_rdm.c
(gdb)
#5  0x00002af7b5a29e19 in fi_ibv_rdm_tagged_cq_read (cq=<value optimized 
out>,
     buf=<value optimized out>, count=<value optimized out>)
     at prov/verbs/src/ep_rdm/verbs_cq_ep_rdm.c:98
98      in prov/verbs/src/ep_rdm/verbs_cq_ep_rdm.c
(gdb)
#6  0x00002af7b57c83c7 in ompi_mtl_ofi_progress_no_inline ()
    from /apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/openmpi/mca_mtl_ofi.so
(gdb)
#7  0x00002af7afab5f2a in opal_progress ()
    from /apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libopen-pal.so.13
(gdb)
#8  0x00002af7b4b5f055 in mca_pml_cm_recv ()
    from /apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/openmpi/mca_pml_cm.so
(gdb)
#9  0x00002af7a962022c in PMPI_Recv ()
    from /apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libmpi.so.12
(gdb)

#10 0x00002af7aab4b297 in libMesh::Parallel::Status 
libMesh::Parallel::Communicator::receive<unsigned int>(unsigned int, 
std::vector<unsigned int, std::allocator<unsigned int> >&, 
libMesh::Parallel::DataTypenst&, libMesh::Parallel::MessageTag const&) 
const () at ./include/libmesh/parallel_implementation.h:2649
2649      libmesh_call_mpi
(gdb)
#11 0x00002af7aab45a19 in void 
libMesh::Parallel::Communicator::send_receive<unsigned int, unsigned 
int>(unsigned int, std::vector<unsigned int, std::allocator<unsigned 
int> > const&, libMesh::Parallel::Datae const&, unsigned int, 
std::vector<unsigned int, std::allocator<unsigned int> >&, 
libMesh::Parallel::DataType const&, libMesh::Parallel::MessageTag 
const&, libMesh::Parallel::MessageTag const&) const () at 
nclude/libmesh/parallel_implementation.h:2788
2788      this->receive (source_processor_id, recv, type2, recv_tag);
(gdb)
#12 0x00002af7aab40643 in void 
libMesh::Parallel::Communicator::send_receive<unsigned int>(unsigned 
int, std::vector<unsigned int, std::allocator<unsigned int> > const&, 
unsigned int, std::vector<unsigned instd::allocator<unsigned int> >&, 
libMesh::Parallel::MessageTag const&, libMesh::
2861      this->send_receive (dest_processor_id, sendvec,
(gdb)
#13 0x00002af7aaf28774 in unsigned int 
libMesh::DistributedMesh::renumber_dof_ob
     at src/mesh/distributed_mesh.C:1090
1090              this->comm().send_receive(procup, requested_ids[procup],
(gdb)
#14 0x00002af7aaf2337b in 
libMesh::DistributedMesh::renumber_nodes_and_elements(
1267      _n_elem = this->renumber_dof_objects (this->_elements);
(gdb)
#15 0x00002af7aaf84a22 in libMesh::MeshBase::prepare_for_use(bool, bool) ()
     at src/mesh/mesh_base.C:209
209         this->renumber_nodes_and_elements();
(gdb)
#16 0x00002af7aafe36b8 in 
libMesh::MeshTools::Generation::build_cube(libMesh::Un) ()
     at src/mesh/mesh_generation.C:1426
1426      mesh.prepare_for_use (/*skip_renumber =*/ false);
(gdb)
#17 0x0000000000417a99 in main () at create_mesh.C:19
19          MeshTools::Generation::build_cube (mesh, 5, 5, 5);
(gdb)


If you want me to print something in the debugger, I'll be happy to do it.
Michael.


On 09/06/2017 11:32 AM, Roy Stogner wrote:
>
> On Wed, 6 Sep 2017, Michael Povolotskyi wrote:
>
>> I found that if I rebuilt everything with MPICH, instead of using 
>> installed openmpi, then everything works perfectly. Is libmesh 
>> supposed to work with openmpi? If yes, and I can recompile it again 
>> and produce the stack trace.
>
> libMesh does work with openmpi; I'm using MPICH2 right now but I was
> using OpenMPI up until a month or two ago, and the problem that made
> me switch was troublesome MPI_Abort behavior, not any bugs in regular
> operation.
>
> However, libMesh does *not* work with mixed MPI installs.  Running an
> openmpi-compiled libMesh with an mpich-based mpiexec, for instance,
> will typically run N 1-processor jobs rather than 1 N-processor job.
> If you somehow managed to *link* to multiple MPI versions (seems very
> unlikely, but perhaps via another dependency which was built against a
> different version?) then all bets are off.
>
> Either way, no need for a stack trace if you've found a workaround.
> ---
> Roy

Re: [Libmesh-users] question about Point locator

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

On 09/06/2017 02:41 PM, John Peterson wrote:
>
>
> On Tue, Sep 5, 2017 at 10:51 PM, Michael Povolotskyi 
> <mpo...@pu... <mailto:mpo...@pu...>> wrote:
>
>     Dear Libmesh community,
>
>     I have a question about PointLocator.
>
>     Does it search among all active elements or among local active
>     elements?
>
>
>
> The PointLocatorBase::build() function takes a PointLocatorType enum 
> which can be one of {TREE, TREE_ELEMENTS, TREE_LOCAL_ELEMENTS}.
>
> The PointLocator in Meshbase is built with the TREE_ELEMENTS enumeration:
>
> _point_locator.reset (PointLocatorBase::build(TREE_ELEMENTS, 
> *this).release());
>
> so that will look search in all elements in the mesh by default.
>
> -- 
> John

Thank you, it is clear now.

Michael.

Re: [Libmesh-users] question about Point locator

[John P. <jwp...@gm...>]

On Tue, Sep 5, 2017 at 10:51 PM, Michael Povolotskyi <mpo...@pu...>
wrote:

> Dear Libmesh community,
>
> I have a question about PointLocator.
>
> Does it search among all active elements or among local active elements?
>


The PointLocatorBase::build() function takes a PointLocatorType enum which
can be one of {TREE, TREE_ELEMENTS, TREE_LOCAL_ELEMENTS}.

The PointLocator in Meshbase is built with the TREE_ELEMENTS enumeration:

_point_locator.reset (PointLocatorBase::build(TREE_ELEMENTS,
*this).release());

so that will look search in all elements in the mesh by default.

-- 
John

Re: [Libmesh-users] problems with running libmesh in parallel

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

On Wed, 6 Sep 2017, Michael Povolotskyi wrote:

> I found that if I rebuilt everything with MPICH, instead of using installed 
> openmpi, then everything works perfectly. Is libmesh supposed to work with 
> openmpi? If yes, and I can recompile it again and produce the stack trace.

libMesh does work with openmpi; I'm using MPICH2 right now but I was
using OpenMPI up until a month or two ago, and the problem that made
me switch was troublesome MPI_Abort behavior, not any bugs in regular
operation.

However, libMesh does *not* work with mixed MPI installs.  Running an
openmpi-compiled libMesh with an mpich-based mpiexec, for instance,
will typically run N 1-processor jobs rather than 1 N-processor job.
If you somehow managed to *link* to multiple MPI versions (seems very
unlikely, but perhaps via another dependency which was built against a
different version?) then all bets are off.

Either way, no need for a stack trace if you've found a workaround.
---
Roy

Re: [Libmesh-users] problems with running libmesh in parallel

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

On 09/06/2017 10:57 AM, Roy Stogner wrote:
>
> On Mon, 4 Sep 2017, Michael Povolotskyi wrote:
>
>> using namespace libMesh;
>>
>> int main (int argc, char ** argv)
>> {
>>   {
>>     LibMeshInit init (argc, argv);
>>
>>     Mesh mesh(init.comm());
>>     MeshTools::Generation::build_cube (mesh, 5, 5, 5);
>>
>>    mesh.print_info();
>>
>>    TecplotIO mesh_output(mesh);
>>
>>     mesh_output.write("mesh.dat");
>>   }
>>
>>  return 0;
>> }
>>
>> The code is frozen at some MPI related calls inside 
>> mesh.prepare_for_use(). The libmesh was configured as:
>>
>> ./configure PETSC_DIR=/home/mpovolot/photonics/petsc/petsc-3.7.6 
>> MPIHOME=/apps/conte/openmpi/1.10.1/gcc-5.2.0 
>> PETSC_ARCH=arch-linux2-cxx-opt F77=mpif77 CC=mpicc GCC=mpicc 
>> CXX=mpiCC --disable-vtk 
>> --with-vtk-include=/depot/itap/tsaiwei/apps/vtk/6.3.0_test/include/vtk-6.3/ 
>> --with-vtk-lib=/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib 
>> --enable-tetgen --enable-triangle --enable-tecplot --disable-nemesis 
>> --disable-strict-lgpl --enable-parmesh --enable-amr 
>> --enable-shared=yes --disable-glibcxx-debugging 
>> --prefix=/home/mpovolot/photonics/libmesh/libmesh-1.2.0/built
>>
>> Could you, please, tell me what is wrong with my code?
>
> I don't see a thing wrong with it.  Are you running in dbg mode, and
> if not, would you?  Could you give us a stack trace pointing to where
> each of two processes ends up frozen?
>
> Thanks,
> ---
> Roy

Thank you,

I found that if I rebuilt everything with MPICH, instead of using 
installed openmpi, then everything works perfectly. Is libmesh supposed 
to work with openmpi? If yes, and I can recompile it again and produce 
the stack trace.

Michael.

Re: [Libmesh-users] problems with running libmesh in parallel

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

On Mon, 4 Sep 2017, Michael Povolotskyi wrote:

> using namespace libMesh;
>
> int main (int argc, char ** argv)
> {
>   {
>     LibMeshInit init (argc, argv);
>
>     Mesh mesh(init.comm());
>     MeshTools::Generation::build_cube (mesh, 5, 5, 5);
>
>    mesh.print_info();
>
>    TecplotIO mesh_output(mesh);
>
>     mesh_output.write("mesh.dat");
>   }
>
>  return 0;
> }
>
> The code is frozen at some MPI related calls inside mesh.prepare_for_use(). 
> The libmesh was configured as:
>
> ./configure PETSC_DIR=/home/mpovolot/photonics/petsc/petsc-3.7.6 
> MPIHOME=/apps/conte/openmpi/1.10.1/gcc-5.2.0 PETSC_ARCH=arch-linux2-cxx-opt 
> F77=mpif77 CC=mpicc GCC=mpicc CXX=mpiCC --disable-vtk 
> --with-vtk-include=/depot/itap/tsaiwei/apps/vtk/6.3.0_test/include/vtk-6.3/ 
> --with-vtk-lib=/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib --enable-tetgen 
> --enable-triangle --enable-tecplot --disable-nemesis --disable-strict-lgpl 
> --enable-parmesh --enable-amr --enable-shared=yes --disable-glibcxx-debugging 
> --prefix=/home/mpovolot/photonics/libmesh/libmesh-1.2.0/built
>
> Could you, please, tell me what is wrong with my code?

I don't see a thing wrong with it.  Are you running in dbg mode, and
if not, would you?  Could you give us a stack trace pointing to where
each of two processes ends up frozen?

Thanks,
---
Roy

[Libmesh-users] question about Point locator

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

Dear Libmesh community,

I have a question about PointLocator.

Does it search among all active elements or among local active elements?

Thank you,

Michael.

Re: [Libmesh-users] question about transient linear system

[John P. <jwp...@gm...>]

On Tue, Sep 5, 2017 at 9:36 AM, Michael Povolotskyi <mpo...@pu...>
wrote:

> Thank you.
> I found that before calling system.update() I have to call
> system.solution->close()
> Is this technically correct?
>

Yes, setting values manually in the vector leaves it in the "unassembled"
state, so before any communication can occur (in system.update()) the
solution vector must be closed.

-- 
John

Re: [Libmesh-users] question about transient linear system

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

Thank you.
I found that before calling system.update() I have to call 
system.solution->close()
Is this technically correct?
Michael.



On 09/05/2017 11:17 AM, John Peterson wrote:
>
>
> On Mon, Sep 4, 2017 at 12:46 PM, Michael Povolotskyi 
> <mpo...@pu... <mailto:mpo...@pu...>> wrote:
>
>     Hello,
>
>     I'm going to solve transient linear system.
>
>     I need to set initial conditions.
>
>     I'm going to code a special initialization function for this.
>
>     In this function, do I need to set system.solution or
>     system.current_local_solution?
>
>
> system.solution and then call system.update() to set values in 
> current_local_solution. As far as I know, one never sets values in 
> current_local_solution directly...
>
> Even better would be to call system.project_solution() with an 
> appropriate function pointer object. This will take care of both for 
> you. See, for example, the init_cd() function in 
> transient/transient_ex1/transient_ex1.C.
>
> -- 
> John

Re: [Libmesh-users] question about transient linear system

[John P. <jwp...@gm...>]

On Mon, Sep 4, 2017 at 12:46 PM, Michael Povolotskyi <mpo...@pu...>
wrote:

> Hello,
>
> I'm going to solve transient linear system.
>
> I need to set initial conditions.
>
> I'm going to code a special initialization function for this.
>
> In this function, do I need to set system.solution or
> system.current_local_solution?
>

system.solution and then call system.update() to set values in
current_local_solution. As far as I know, one never sets values in
current_local_solution directly...

Even better would be to call system.project_solution() with an appropriate
function pointer object. This will take care of both for you. See, for
example, the init_cd() function in transient/transient_ex1/transient_ex1.C.

-- 
John

[Libmesh-users] problems with running libmesh in parallel

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

Dear Libmesh developers and users,

I have a simple code which runs perfectly on a single MPI process, but 
fails on more than one MPI process.

The code is as follows:

using namespace libMesh;

int main (int argc, char ** argv)
{
   {
     LibMeshInit init (argc, argv);

     Mesh mesh(init.comm());
     MeshTools::Generation::build_cube (mesh, 5, 5, 5);

     mesh.print_info();

     TecplotIO mesh_output(mesh);

     mesh_output.write("mesh.dat");
   }

   return 0;
}

The code is frozen at some MPI related calls inside 
mesh.prepare_for_use(). The libmesh was configured as:

./configure PETSC_DIR=/home/mpovolot/photonics/petsc/petsc-3.7.6 
MPIHOME=/apps/conte/openmpi/1.10.1/gcc-5.2.0 
PETSC_ARCH=arch-linux2-cxx-opt F77=mpif77 CC=mpicc GCC=mpicc CXX=mpiCC 
--disable-vtk 
--with-vtk-include=/depot/itap/tsaiwei/apps/vtk/6.3.0_test/include/vtk-6.3/ 
--with-vtk-lib=/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib 
--enable-tetgen --enable-triangle --enable-tecplot --disable-nemesis 
--disable-strict-lgpl --enable-parmesh --enable-amr --enable-shared=yes 
--disable-glibcxx-debugging 
--prefix=/home/mpovolot/photonics/libmesh/libmesh-1.2.0/built

Could you, please, tell me what is wrong with my code?

Thank you,

Michael.

[Libmesh-users] question about transient linear system

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

Hello,

I'm going to solve transient linear system.

I need to set initial conditions.

I'm going to code a special initialization function for this.

In this function, do I need to set system.solution or 
system.current_local_solution?

Thank you,

Michael.

Re: [Libmesh-users] problem with vtk linking

[John P. <jwp...@gm...>]

On Wed, Aug 30, 2017 at 5:02 PM, Michael Povolotskyi <mpo...@pu...>
wrote:

> I can try that, thank you.
>
> Still, why some libraries are found and some are not? Should this be
> covered when libmesh_opt.so is built?
>

This depends on the way VTK was built on your system. When I build VTK, I
use the following Cmake flags to "compile in" the installation paths to the
libraries themselves. This works on Mac, on Linux your mileage may vary...

      -DCMAKE_INSTALL_RPATH:STRING="/path/to/vtk/lib"
      -DCMAKE_INSTALL_RPATH_USE_LINK_PATH:BOOL=ON
      -DCMAKE_INSTALL_NAME_DIR:STRING="/path/to/vtk/lib"

Additionally, the use of system modules (e.g. module load vtk) is
recommended to control the contents of environment variables like
LD_LIBRARY_PATH in a transparent way.

-- 
John

Re: [Libmesh-users] problem with vtk linking

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

I can try that, thank you.

Still, why some libraries are found and some are not? Should this be 
covered when libmesh_opt.so is built?

Since I use vtk via libmesh only, ideally I do not need to know where 
the vtk libraries are located.

Michael.


On 8/30/2017 6:58 PM, John Peterson wrote:
>
>
> On Wed, Aug 30, 2017 at 4:47 PM, Michael Povolotskyi 
> <mpo...@pu... <mailto:mpo...@pu...>> wrote:
>
>     Dear Libmesh team,
>
>     I have built libmesh with vtk support.
>
>     The build runs okay, but the library can see path to some of the
>     necessary vtk libraries only, and some libraries are not found.
>
>     Is it possible to fix this somehow? Thank you!
>
>
> Have you tried adding /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/ to 
> your LD_LIBRARY_PATH? On Linux that is usually the fix when dynamic 
> libraries are not found.
>
>
> I'm not exactly sure what's going on, but the libraries that are 
> found, i.e.
>
>     libvtkIOCore-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOCore-6.3.so.1 
> (0x00002b750b77b000)
>     libvtkCommonCore-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkCommonCore-6.3.so.1 
> (0x00002b750b9e6000)
>     libvtkCommonDataModel-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkCommonDataModel-6.3.so.1 
> (0x00002b750c007000)
>     libvtkFiltersCore-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkFiltersCore-6.3.so.1 
> (0x00002b750c513000)
>     libvtkIOXML-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOXML-6.3.so.1 
> (0x00002b750cb81000)
>     libvtkImagingCore-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkImagingCore-6.3.so.1 
> (0x00002b750ce4a000)
>     libvtkIOImage-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOImage-6.3.so.1 
> (0x00002b750d2aa000)
>     libvtkImagingMath-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkImagingMath-6.3.so.1 
> (0x00002b750d600000)
>     libvtkIOParallelXML-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOParallelXML-6.3.so.1 
> (0x00002b750d83d000)
>     libvtkParallelMPI-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkParallelMPI-6.3.so.1 
> (0x00002b750da63000)
>     libvtkParallelCore-6.3.so.1 => 
> /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkParallelCore-6.3.so.1 
> (0x00002b750dc84000)
>
> look like the list of VTK libraries libmesh needs to link against (see 
> vtk.m4). The ones that aren't found, e.g
>
> libvtkCommonExecutionModel-6.3.so.1 => not found
> libvtkCommonMisc-6.3.so.1 => not found
>
> are perhaps getting pulled in by those?
>
> -- 
> John

Re: [Libmesh-users] problem with vtk linking

[John P. <jwp...@gm...>]

On Wed, Aug 30, 2017 at 4:47 PM, Michael Povolotskyi <mpo...@pu...>
wrote:

> Dear Libmesh team,
>
> I have built libmesh with vtk support.
>
> The build runs okay, but the library can see path to some of the necessary
> vtk libraries only, and some libraries are not found.
>
> Is it possible to fix this somehow? Thank you!
>

Have you tried adding /depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/ to your
LD_LIBRARY_PATH? On Linux that is usually the fix when dynamic libraries
are not found.


I'm not exactly sure what's going on, but the libraries that are found, i.e.

    libvtkIOCore-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOCore-6.3.so.1
(0x00002b750b77b000)
    libvtkCommonCore-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkCommonCore-6.3.so.1
(0x00002b750b9e6000)
    libvtkCommonDataModel-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkCommonDataModel-6.3.so.1
(0x00002b750c007000)
    libvtkFiltersCore-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkFiltersCore-6.3.so.1
(0x00002b750c513000)
    libvtkIOXML-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOXML-6.3.so.1
(0x00002b750cb81000)
    libvtkImagingCore-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkImagingCore-6.3.so.1
(0x00002b750ce4a000)
    libvtkIOImage-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOImage-6.3.so.1
(0x00002b750d2aa000)
    libvtkImagingMath-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkImagingMath-6.3.so.1
(0x00002b750d600000)
    libvtkIOParallelXML-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOParallelXML-6.3.so.1
(0x00002b750d83d000)
    libvtkParallelMPI-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkParallelMPI-6.3.so.1
(0x00002b750da63000)
    libvtkParallelCore-6.3.so.1 =>
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkParallelCore-6.3.so.1
(0x00002b750dc84000)

look like the list of VTK libraries libmesh needs to link against (see
vtk.m4). The ones that aren't found, e.g

libvtkCommonExecutionModel-6.3.so.1 => not found
libvtkCommonMisc-6.3.so.1 => not found

are perhaps getting pulled in by those?

-- 
John

[Libmesh-users] problem with vtk linking

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

Dear Libmesh team,

I have built libmesh with vtk support.

The build runs okay, but the library can see path to some of the 
necessary vtk libraries only, and some libraries are not found.

Is it possible to fix this somehow? Thank you!

Michael.

p.s.

Here is the ldd output and the config.log is attached.

ldd /home/mpovolot/photonics/libmesh/libmesh-1.2.0/built/lib/libmesh_opt.so
     linux-vdso.so.1 =>  (0x00007ffccefa6000)
     libnetcdf.so.7 => 
/home/mpovolot/photonics/libmesh/libmesh-1.2.0/built/lib/libnetcdf.so.7 
(0x00002b750af11000)
     libcurl.so.4 => /usr/lib64/libcurl.so.4 (0x00002b750b25b000)
     libglpk.so.0 => /usr/lib64/libglpk.so.0 (0x00002b750b4b0000)
     libvtkIOCore-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOCore-6.3.so.1 
(0x00002b750b77b000)
     libvtkCommonCore-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkCommonCore-6.3.so.1 
(0x00002b750b9e6000)
     libvtkCommonDataModel-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkCommonDataModel-6.3.so.1 
(0x00002b750c007000)
     libvtkFiltersCore-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkFiltersCore-6.3.so.1 
(0x00002b750c513000)
     libvtkIOXML-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOXML-6.3.so.1 
(0x00002b750cb81000)
     libvtkImagingCore-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkImagingCore-6.3.so.1 
(0x00002b750ce4a000)
     libvtkIOImage-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOImage-6.3.so.1 
(0x00002b750d2aa000)
     libvtkImagingMath-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkImagingMath-6.3.so.1 
(0x00002b750d600000)
     libvtkIOParallelXML-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkIOParallelXML-6.3.so.1 
(0x00002b750d83d000)
     libvtkParallelMPI-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkParallelMPI-6.3.so.1 
(0x00002b750da63000)
     libvtkParallelCore-6.3.so.1 => 
/depot/itap/tsaiwei/apps/vtk/6.3.0_test/lib/libvtkParallelCore-6.3.so.1 
(0x00002b750dc84000)
     libz.so.1 => /lib64/libz.so.1 (0x00002b750ded7000)
     libmkl_scalapack_lp64.so => 
/apps/rhel6/intel/compilers_and_libraries_2016.1.150/linux/mkl/lib/intel64/libmkl_scalapack_lp64.so 
(0x00002b750e0ee000)
     libmkl_blacs_intelmpi_lp64.so => 
/apps/rhel6/intel/compilers_and_libraries_2016.1.150/linux/mkl/lib/intel64/libmkl_blacs_intelmpi_lp64.so 
(0x00002b750e9da000)
     libmkl_intel_lp64.so => 
/apps/rhel6/intel/compilers_and_libraries_2016.1.150/linux/mkl/lib/intel64/libmkl_intel_lp64.so 
(0x00002b750ec16000)
     libmkl_gnu_thread.so => 
/apps/rhel6/intel/compilers_and_libraries_2016.1.150/linux/mkl/lib/intel64/libmkl_gnu_thread.so 
(0x00002b750f555000)
     libmkl_core.so => 
/apps/rhel6/intel/compilers_and_libraries_2016.1.150/linux/mkl/lib/intel64/libmkl_core.so 
(0x00002b7510488000)
     libX11.so.6 => /usr/lib64/libX11.so.6 (0x00002b7511dbf000)
     libhwloc.so.5 => /usr/lib64/libhwloc.so.5 (0x00002b75120fd000)
     libssl.so.10 => /usr/lib64/libssl.so.10 (0x00002b7512325000)
     libcrypto.so.10 => /usr/lib64/libcrypto.so.10 (0x00002b7512591000)
     libmpi_usempif08.so.11 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libmpi_usempif08.so.11 
(0x00002b7512976000)
     libmpi_usempi_ignore_tkr.so.6 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libmpi_usempi_ignore_tkr.so.6 
(0x00002b7512ba5000)
     libmpi_mpifh.so.12 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libmpi_mpifh.so.12 
(0x00002b7512dab000)
     libgfortran.so.3 => /apps/rhel6/gcc/5.2.0/lib64/libgfortran.so.3 
(0x00002b7513003000)
     libquadmath.so.0 => /apps/rhel6/gcc/5.2.0/lib64/libquadmath.so.0 
(0x00002b7513323000)
     libgomp.so.1 => /apps/rhel6/gcc/5.2.0/lib64/libgomp.so.1 
(0x00002b7513561000)
     libmpi_cxx.so.1 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libmpi_cxx.so.1 
(0x00002b7513783000)
     libmpi.so.12 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libmpi.so.12 (0x00002b751399d000)
     libosmcomp.so.3 => /usr/lib64/libosmcomp.so.3 (0x00002b7513c7b000)
     libibverbs.so.1 => /usr/lib64/libibverbs.so.1 (0x00002b7513e89000)
     libopen-rte.so.12 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libopen-rte.so.12 
(0x00002b7514097000)
     libopen-pal.so.13 => 
/apps/conte/openmpi/1.10.1/gcc-5.2.0/lib/libopen-pal.so.13 
(0x00002b7514312000)
     libnuma.so.1 => /usr/lib64/libnuma.so.1 (0x00002b75145ef000)
     libpciaccess.so.0 => /usr/lib64/libpciaccess.so.0 (0x00002b75147fa000)
     libdl.so.2 => /lib64/libdl.so.2 (0x00002b7514a03000)
     librt.so.1 => /lib64/librt.so.1 (0x00002b7514c08000)
     libutil.so.1 => /lib64/libutil.so.1 (0x00002b7514e10000)
     libstdc++.so.6 => /apps/rhel6/gcc/5.2.0/lib64/libstdc++.so.6 
(0x00002b7515013000)
     libm.so.6 => /lib64/libm.so.6 (0x00002b75153a1000)
     libpthread.so.0 => /lib64/libpthread.so.0 (0x00002b7515625000)
     libc.so.6 => /lib64/libc.so.6 (0x00002b7515842000)
     libgcc_s.so.1 => /apps/rhel6/gcc/5.2.0/lib64/libgcc_s.so.1 
(0x00002b7515bd7000)
     libidn.so.11 => /lib64/libidn.so.11 (0x00002b7515ded000)
     libldap-2.4.so.2 => /lib64/libldap-2.4.so.2 (0x00002b751601f000)
     libgssapi_krb5.so.2 => /lib64/libgssapi_krb5.so.2 (0x00002b7516270000)
     libkrb5.so.3 => /lib64/libkrb5.so.3 (0x00002b75164b4000)
     libk5crypto.so.3 => /lib64/libk5crypto.so.3 (0x00002b751679b000)
     libcom_err.so.2 => /lib64/libcom_err.so.2 (0x00002b75169c8000)
     libssl3.so => /usr/lib64/libssl3.so (0x00002b7516bcc000)
     libsmime3.so => /usr/lib64/libsmime3.so (0x00002b7516e0f000)
     libnss3.so => /usr/lib64/libnss3.so (0x00002b751703c000)
     libnssutil3.so => /usr/lib64/libnssutil3.so (0x00002b751737c000)
     libplds4.so => /lib64/libplds4.so (0x00002b75175a8000)
     libplc4.so => /lib64/libplc4.so (0x00002b75177ad000)
     libnspr4.so => /lib64/libnspr4.so (0x00002b75179b2000)
     libssh2.so.1 => /usr/lib64/libssh2.so.1 (0x00002b7517bf0000)
     /lib64/ld-linux-x86-64.so.2 (0x000000321f400000)
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkzlib-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkIOGeometry-6.3.so.1 => not found
     libvtkIOXMLParser-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkjpeg-6.3.so.1 => not found
     libvtkpng-6.3.so.1 => not found
     libvtktiff-6.3.so.1 => not found
     libvtkmetaio-6.3.so.1 => not found
     libvtkDICOMParser-6.3.so.1 => not found
     libvtkzlib-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkIOGeometry-6.3.so.1 => not found
     libvtkIOXMLParser-6.3.so.1 => not found
     libvtkIOLegacy-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkIOLegacy-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libvtkIOLegacy-6.3.so.1 => not found
     libvtksys-6.3.so.1 => not found
     libvtkCommonExecutionModel-6.3.so.1 => not found
     libvtkCommonMisc-6.3.so.1 => not found
     libvtkCommonSystem-6.3.so.1 => not found
     libvtkCommonTransforms-6.3.so.1 => not found
     libvtkCommonMath-6.3.so.1 => not found
     libxcb.so.1 => /usr/lib64/libxcb.so.1 (0x00002b7517e31000)
     libpci.so.3 => /lib64/libpci.so.3 (0x00002b7518050000)
     libxml2.so.2 => /usr/lib64/libxml2.so.2 (0x00002b751825c000)
     libibumad.so.3 => /usr/lib64/libibumad.so.3 (0x00002b75185b0000)
     liblber-2.4.so.2 => /lib64/liblber-2.4.so.2 (0x00002b75187b8000)
     libresolv.so.2 => /lib64/libresolv.so.2 (0x00002b75189c8000)
     libsasl2.so.2 => /usr/lib64/libsasl2.so.2 (0x00002b7518be2000)
     libkrb5support.so.0 => /lib64/libkrb5support.so.0 (0x00002b7518dfc000)
     libkeyutils.so.1 => /lib64/libkeyutils.so.1 (0x00002b7519008000)
     libXau.so.6 => /usr/lib64/libXau.so.6 (0x00002b751920b000)
     libcrypt.so.1 => /lib64/libcrypt.so.1 (0x00002b751940f000)
     libselinux.so.1 => /lib64/libselinux.so.1 (0x00002b7519646000)
     libfreebl3.so => /lib64/libfreebl3.so (0x00002b7519866000)

Re: [Libmesh-users] Trouble with boundary integration

[Renato P. <re...@gm...>]

Hi Mike,

Thanks once again. I will do the calculations.

The problem is not in the normals, at least not only.
I am pretty sure the gradient at the quadrature point is wrong.
It must be around 1E+5. Must not cancel out.
The values in the dofs are reasonable.
Shape function issue. Must be.

I dont believe in any non-physical stuff in this case. It is such a simple
case!

Will let you know about my findings!
Thanks

Renato


On Sat, Aug 26, 2017 at 10:43 AM, Mike Marchywka <mar...@ho...>
wrote:

>
> It should not be too hard to just write the value u = s1*u1+s2*u2 ... and
> differentiate s_i to
> get the grad at that q point and see if the code matches the paper and
> look for any obvious missing factors or terms.
> I'm just doing QUAD4 - linear shape functions and rectangles with a lot of
> my own code-
> but any polynomial not much harder. It looks like you have a sum of values
> that nearly cancel out.
> I know I did some interpolation and the linear worked better than
> quadratic due to nonphysical
> excess curvature ( often leading to nonphysical negative values between
> nodes, In some cases I thought interpolating
> the logs would work better but never got back to it) .  But with an
> analytical results again comparisons should be
> informative.
> Also I guess with an integral over a closed surface you probably
> need to check the normal sign but that should be obvious just plotting the
> flux at each node
> if you expect it to be uniform rather than oscillating :)
>
> ________________________________________
> From: Renato Poli <re...@gm...>
> Sent: Friday, August 25, 2017 7:46 PM
> To: Mike Marchywka
> Cc: John Peterson; libmesh-users
> Subject: Re: [Libmesh-users] Trouble with boundary integration
>
> On Fri, Aug 25, 2017 at 8:08 PM, Mike Marchywka <mar...@ho...
> <mailto:mar...@ho...>> wrote:
> >Any hint on how to easily reproduce the shape function calculations by
> hand (that is, manually producing the dphi_face values)?
> Very carefully :) I ended up writing some polynomial classes and doing the
> integrals and derivatives for QUAD4 that way. With my
> rational number class it was easy to get exact results and make generated
> code for the simple things like Laplacian of rectangle.
> That seems a good idea. I am using TRI6.
> The results I am getting are so off the expected that I won't matter
> getting approximations for now as reference.
> This is the reason why I think I got a conceptual issue. It is really far
> from target.
>
> I guess you could probably do the same with MatLab or Mathematica or maybe
> R.  I saw the point by point quadrature
> thing in the examples but AFAICT in many cases you can do it exactly and
> just skip a lot of that. The Laplacian
> is just an integral of a second derivative which would be zero for linear
> shape functions except for integration by
> parts which then gives you an easily computable integral. I'm using all
> analytical integrals now- just putting numbers
> into the system matrix and rhs with QUAD4 elements- but have a possible
> model problem that probably needs point by point  integration :)
>
> But do the derivatives look right? The normals? The solution? Does it work
> without AMR ? p? h? etc...
> The normals look right. I have simplified the mesh to have right angles
> and thus simple normals (1,0,0) or (-1,0,0) etc.
> I am not using AMR at all. Mesh is dynamically generated with a couple of
> parameters (maximum triangle size and internal angles, for example).
>
> What does work, if anything?
> Everything works in a wrong way. :-)
>
> Look this example:
>
> +++ Triangle coordinates (each dof)
> [0]    (x,y,z)=(      10,        5,        0)
> [1]    (x,y,z)=(      10,       10,        0)
> [2]    (x,y,z)=( 6.73058,  6.91831,        0)
> [3]    (x,y,z)=(      10,      7.5,        0)
> [4]    (x,y,z)=( 8.36529,  8.45915,        0)
> [5]    (x,y,z)=( 8.36529,  5.95915,        0)
>
> See the highlighted system solutions for dof[1] and dof[2].
> It shows a delta of 0.1206E7 in (dx=10-6.73=3.27).
> Something like grad_x=0.0369E7.
> I dont get why the shape function and the summation sort of kills this
> gradient...
> Still struggling here...
>
> +++ intermediate calculations
>      qface_point = (x,y,z)=(      10,      7.5,        0)
>         == i=0 ===
>         dof_indices[0]: 12
>         dphi_face[0][0]: (x,y,z)=(0.188516,     -0.2,       -0)
>         _sys_solution[0->12]: 3e+07
>         _sys_solution[0->12] * dphi_face[0]: (x,y,z)=(5.65548e+06,
>  -6e+06,       -0)
>         == i=1 ===
>         dof_indices[1]: 13
>         dphi_face[1][0]: (x,y,z)=(0.117348,      0.2,        0)
>         _sys_solution[1->13]: 3e+07
> /// DOF[1] => 3E7 @ (10,5) /// <===
>         _sys_solution[1->13] * dphi_face[1]: (x,y,z)=(3.52045e+06,
> 6e+06,        0)
>         == i=2 ===
>         dof_indices[2]: 14
>         dphi_face[2][0]: (x,y,z)=(0.305864,       -0,       -0)
>         _sys_solution[2->14]: 3.12065e+07
> /// DOF[2] => 3.12E7 @ (6.73,  6.91) /// <===
>         _sys_solution[2->14] * dphi_face[2]: (x,y,z)=(9.54495e+06,
>  -0,       -0)
>         == i=3 ===
>         dof_indices[3]: 15
>         dphi_face[3][0]: (x,y,z)=(0.611729,        0,        0)
>         _sys_solution[3->15]: 3e+07
>         _sys_solution[3->15] * dphi_face[3]: (x,y,z)=(1.83519e+07,
> 0,        0)
>         == i=4 ===
>         dof_indices[4]: 16
>         dphi_face[4][0]: (x,y,z)=(-0.611729,        0,        0)
>         _sys_solution[4->16]: 3.02508e+07
>         _sys_solution[4->16] * dphi_face[4]: (x,y,z)=(-1.85053e+07,
> 0,        0)
>         == i=5 ===
>         dof_indices[5]: 17
>         dphi_face[5][0]: (x,y,z)=(-0.611729,        0,        0)
>         _sys_solution[5->17]: 3.03524e+07
>         _sys_solution[5->17] * dphi_face[5]: (x,y,z)=(-1.85674e+07,
> 0,        0)
>         =======
>    JxW_face[0]: 5
>    grad_at_qn=        (x,y,z)=(-6.29947e-06, -7.7514e-06,        0)
>    normals[0]= (x,y,z)=(       1,        0,       -0)
>    apply_normal=      -3.14973e-05
>
>
>
>
> ________________________________________
> From: Renato Poli <re...@gm...<mailto:re...@gm...>>
> Sent: Friday, August 25, 2017 6:12 PM
> To: Mike Marchywka
> Cc: John Peterson; libmesh-users
> Subject: Re: [Libmesh-users] Trouble with boundary integration
>
> Hi Mike,
>
> Thanks for the reply.
> Any hint on how to easily reproduce the shape function calculations by
> hand (that is, manually producing the dphi_face values)?
>
> On Fri, Aug 25, 2017 at 6:17 PM, Mike Marchywka <mar...@ho...
> <mailto:mar...@ho...><mailto:mar...@ho...<mailto:m
> arc...@ho...>>> wrote:
> I've been trying to learn this myself and encounter a variety of issues
> like this. If you have
> an exact solution, why not just print all the pieces- values and normals
> etc-  and see how they differ
> point by point instead of
> trying to do an inverse problem of guessing what is wrong with an integral
> ?
>
> Well - that is what I am trying to do.
> I isolated a single element and, as the problem is symmetrical, it should
> produce a non-zero result.
> There is not a lot of freedom in the calculation.
> It should be a simple multiplication of a bunch of numbers...
>
> Varying the mesh size seems to be useful sometimes.  How do the element
> sizes
> compare to the radius of curvature ? Can you test the code on a simpler
> geometry such
> as between two planes?
>
> Agreed. I spent two days on that, varying all sort of stuff.
> I tried from very small elements, to huge ones.
> I think my problem might be related to the order of the shape function or
> the gauss quadrature.
> Really not sure.
> Any hint on how to easily reproduce the shape function calculations by
> hand (that is, manually producing the dphi_face values)?
>
> AFAICT in the general
> case doing integrals of normal components over surfaces is not always easy
> with p and h refinements etc. Although maybe someone can clarify a bit.
>
> I am using CGAL to generate and refine the mesh. It is quite flexible.
> I am really considering that I have a conceptual issue.
> It seems there is something really close to me that is fooling me ...
>
> I will keep trying.
>
> Please let me know if anyone has any suggestion ...
>
> Thanks,
> Renato
>
>
> Thanks.
>
> ________________________________________
> From: Renato Poli <re...@gm...<mailto:re...@gm...><mailto:
> re...@gm...<mailto:re...@gm...>>>
> Sent: Friday, August 25, 2017 3:02 PM
> To: John Peterson
> Cc: libmesh-users
> Subject: Re: [Libmesh-users] Trouble with boundary integration
>
> Hi John,
>
> I am still struggling and running out of ideas here.
> I wrote down some debugging - perhaps you can find something weird.
> The calculations seem correct, at least to what I could understand so far.
>
> But things are zero'ing out.
>
> If, in the mesh creation, I declare the variable as "FIRST" order, things
> do not zero out.
> The integration of the gradient only zeros out if I use:
>         sys.add_variable("u", SECOND);
>
> For face integration, I am using first order gauss quadrature:
>         QGauss qface(dim-1, FIRST);
>
> Any hint? Any debugging suggestion?
>
> Rgds,
> Renato
>
> # Debug[1]: Elem: 4
> # Debug[1]:     (x,y,z)=(       5,       10,        0)
> # Debug[1]:     (x,y,z)=(       0,       10,        0)
> # Debug[1]:     (x,y,z)=( 3.08169,  6.73058,        0)
> # Debug[1]:     (x,y,z)=(     2.5,       10,        0)
> # Debug[1]:     (x,y,z)=( 1.54085,  8.36529,        0)
> # Debug[1]:     (x,y,z)=( 4.04085,  8.36529,        0)
> # Debug[1]: num_dofs=6
> # Debug[1]:     Side: 0
> # Debug[1]:     Boundary: 1
> # Debug[1]:          qface_point = (x,y,z)=(     2.5,       10,        0)
> # Debug[1]:             == i=0 ===
> # Debug[1]:             dof_indices[0]: 23
> # Debug[1]:             dphi_face[0][0]: (x,y,z)=(     0.2, 0.188516,
> -0)
> # Debug[1]:             _sys_solution[0->23]: 3e+07
> # Debug[1]:             _sys_solution[0->23] * dphi_face[0]: (x,y,z)=(
> 6e+06, 5.65548e+06,       -0)
> # Debug[1]:             == i=1 ===
> # Debug[1]:             dof_indices[1]: 24
> # Debug[1]:             dphi_face[1][0]: (x,y,z)=(    -0.2,
> 0.117348,        0)
> # Debug[1]:             _sys_solution[1->24]: 3e+07
> # Debug[1]:             _sys_solution[1->24] * dphi_face[1]: (x,y,z)=(
> -6e+06, 3.52045e+06,        0)
> # Debug[1]:             == i=2 ===
> # Debug[1]:             dof_indices[2]: 25
> # Debug[1]:             dphi_face[2][0]: (x,y,z)=(       0, 0.305864,
> -0)
> # Debug[1]:             _sys_solution[2->25]: 3.12065e+07
> # Debug[1]:             _sys_solution[2->25] * dphi_face[2]:
> (x,y,z)=(       0, 9.54495e+06,       -0)
> # Debug[1]:             == i=3 ===
> # Debug[1]:             dof_indices[3]: 26
> # Debug[1]:             dphi_face[3][0]: (x,y,z)=(      -0,
> 0.611729,        0)
> # Debug[1]:             _sys_solution[3->26]: 3e+07
> # Debug[1]:             _sys_solution[3->26] * dphi_face[3]: (x,y,z)=(
> -0, 1.83519e+07,        0)
> # Debug[1]:             == i=4 ===
> # Debug[1]:             dof_indices[4]: 27
> # Debug[1]:             dphi_face[4][0]: (x,y,z)=(       0,
> -0.611729,        0)
> # Debug[1]:             _sys_solution[4->27]: 3.02508e+07
> # Debug[1]:             _sys_solution[4->27] * dphi_face[4]:
> (x,y,z)=(       0, -1.85053e+07,        0)
> # Debug[1]:             == i=5 ===
> # Debug[1]:             dof_indices[5]: 28
> # Debug[1]:             dphi_face[5][0]: (x,y,z)=(       0,
> -0.611729,        0)
> # Debug[1]:             _sys_solution[5->28]: 3.03524e+07
> # Debug[1]:             _sys_solution[5->28] * dphi_face[5]:
> (x,y,z)=(       0, -1.85674e+07,        0)
> # Debug[1]:             =======
> # Debug[1]:        JxW_face[0]: 5
> # Debug[1]:        grad_at_qn=        (x,y,z)=(8.95187e-06,
> -6.02752e-06,        0)
> # Debug[1]:        normals[0]= (x,y,z)=(       0,        1,       -0)
> # Debug[1]:        JxW_face[0] * grad_at_qn * normals[0] =
> -3.01376e-05
>
>
> On Fri, Aug 25, 2017 at 2:25 PM, Renato Poli <re...@gm...<mailto:
> re...@gm...><mailto:re...@gm...<mailto:re...@gm...>>>
> wrote:
>
> > Hi John,
> >
> > _sys_solution is the system->solution vector.
> >
> > It is fed by:
> >         _system.solution->localize_to_one( _sys_solution );
> >
> >
> > On Fri, Aug 25, 2017 at 1:48 PM, John Peterson <jwp...@gm...
> <mailto:jwp...@gm...><mailto:jwp...@gm...<mailto:jwpe
> te...@gm...>>>
> > wrote:
> >
> >>
> >>
> >> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli <re...@gm...
> <mailto:re...@gm...><mailto:re...@gm...<mailto:rebpoli@
> gmail.com>>> wrote:
> >>
> >>> Hi,
> >>>
> >>> I need some help here ...
> >>>
> >>> I am having trouble to understand what is going on with a boundary
> >>> integral.
> >>> I have solved a Poisson problem for pressure in a 2D mesh, as a model
> >>> for a
> >>> incompressible water flow in porous media.
> >>> The mesh (Tri6 elements) has a near-circular polygon in the middle
> >>> representing an injection well.
> >>>
> >>> RHS of the equations are zero except in boundaries dofs. Dirichlet
> >>> boundary
> >>> conditions are imposed using a penalty both on the well and on the
> >>> external
> >>> boundaries to represent fixed pressure. (the solution looks ok - seems
> to
> >>> work)
> >>>
> >>> To validate the solution, I am trying to estimate the water flow both
> in
> >>> the well and in the external boundaries (they must be identical and
> match
> >>> the analytical solution). I therefore calculate the normal pressure
> >>> gradient and integrate in the whole boundary path.
> >>>
> >>> When I use first order variable, the results are "almost correct". But
> >>> when
> >>> I use second order , it is absolutely off.
> >>>
> >>> I am using the code below to do the integration. I tried many
> >>> alternatives,
> >>> no success.
> >>>
> >>> Any idea of what I am doing wrong?
> >>>
> >>
> >> The code seems OK to me, although you should not need to call
> >> fe_face->reinit() inside the loop over boundary ids.
> >>
> >> What is _sys_solution?
> >>
> >> --
> >> John
> >>
> >
> >
> ------------------------------------------------------------
> ------------------
> 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:Libmesh
> -u...@li...><mailto:Lib...@li...
> <mailto:Lib...@li...>>
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
>
>

Re: [Libmesh-users] Trouble with boundary integration

[Mike M. <mar...@ho...>]

It should not be too hard to just write the value u = s1*u1+s2*u2 ... and differentiate s_i to
get the grad at that q point and see if the code matches the paper and look for any obvious missing factors or terms.  
I'm just doing QUAD4 - linear shape functions and rectangles with a lot of my own code- 
but any polynomial not much harder. It looks like you have a sum of values that nearly cancel out. 
I know I did some interpolation and the linear worked better than quadratic due to nonphysical 
excess curvature ( often leading to nonphysical negative values between nodes, In some cases I thought interpolating
the logs would work better but never got back to it) .  But with an analytical results again comparisons should be
informative. 
Also I guess with an integral over a closed surface you probably
need to check the normal sign but that should be obvious just plotting the flux at each node
if you expect it to be uniform rather than oscillating :) 

________________________________________
From: Renato Poli <re...@gm...>
Sent: Friday, August 25, 2017 7:46 PM
To: Mike Marchywka
Cc: John Peterson; libmesh-users
Subject: Re: [Libmesh-users] Trouble with boundary integration

On Fri, Aug 25, 2017 at 8:08 PM, Mike Marchywka <mar...@ho...<mailto:mar...@ho...>> wrote:
>Any hint on how to easily reproduce the shape function calculations by hand (that is, manually producing the dphi_face values)?
Very carefully :) I ended up writing some polynomial classes and doing the integrals and derivatives for QUAD4 that way. With my
rational number class it was easy to get exact results and make generated code for the simple things like Laplacian of rectangle.
That seems a good idea. I am using TRI6.
The results I am getting are so off the expected that I won't matter getting approximations for now as reference.
This is the reason why I think I got a conceptual issue. It is really far from target.

I guess you could probably do the same with MatLab or Mathematica or maybe R.  I saw the point by point quadrature
thing in the examples but AFAICT in many cases you can do it exactly and just skip a lot of that. The Laplacian
is just an integral of a second derivative which would be zero for linear shape functions except for integration by
parts which then gives you an easily computable integral. I'm using all analytical integrals now- just putting numbers
into the system matrix and rhs with QUAD4 elements- but have a possible model problem that probably needs point by point  integration :)

But do the derivatives look right? The normals? The solution? Does it work without AMR ? p? h? etc...
The normals look right. I have simplified the mesh to have right angles and thus simple normals (1,0,0) or (-1,0,0) etc.
I am not using AMR at all. Mesh is dynamically generated with a couple of parameters (maximum triangle size and internal angles, for example).

What does work, if anything?
Everything works in a wrong way. :-)

Look this example:

+++ Triangle coordinates (each dof)
[0]    (x,y,z)=(      10,        5,        0)
[1]    (x,y,z)=(      10,       10,        0)
[2]    (x,y,z)=( 6.73058,  6.91831,        0)
[3]    (x,y,z)=(      10,      7.5,        0)
[4]    (x,y,z)=( 8.36529,  8.45915,        0)
[5]    (x,y,z)=( 8.36529,  5.95915,        0)

See the highlighted system solutions for dof[1] and dof[2].
It shows a delta of 0.1206E7 in (dx=10-6.73=3.27).
Something like grad_x=0.0369E7.
I dont get why the shape function and the summation sort of kills this gradient...
Still struggling here...

+++ intermediate calculations
     qface_point = (x,y,z)=(      10,      7.5,        0)
        == i=0 ===
        dof_indices[0]: 12
        dphi_face[0][0]: (x,y,z)=(0.188516,     -0.2,       -0)
        _sys_solution[0->12]: 3e+07
        _sys_solution[0->12] * dphi_face[0]: (x,y,z)=(5.65548e+06,   -6e+06,       -0)
        == i=1 ===
        dof_indices[1]: 13
        dphi_face[1][0]: (x,y,z)=(0.117348,      0.2,        0)
        _sys_solution[1->13]: 3e+07
/// DOF[1] => 3E7 @ (10,5) /// <===
        _sys_solution[1->13] * dphi_face[1]: (x,y,z)=(3.52045e+06,    6e+06,        0)
        == i=2 ===
        dof_indices[2]: 14
        dphi_face[2][0]: (x,y,z)=(0.305864,       -0,       -0)
        _sys_solution[2->14]: 3.12065e+07
/// DOF[2] => 3.12E7 @ (6.73,  6.91) /// <===
        _sys_solution[2->14] * dphi_face[2]: (x,y,z)=(9.54495e+06,       -0,       -0)
        == i=3 ===
        dof_indices[3]: 15
        dphi_face[3][0]: (x,y,z)=(0.611729,        0,        0)
        _sys_solution[3->15]: 3e+07
        _sys_solution[3->15] * dphi_face[3]: (x,y,z)=(1.83519e+07,        0,        0)
        == i=4 ===
        dof_indices[4]: 16
        dphi_face[4][0]: (x,y,z)=(-0.611729,        0,        0)
        _sys_solution[4->16]: 3.02508e+07
        _sys_solution[4->16] * dphi_face[4]: (x,y,z)=(-1.85053e+07,        0,        0)
        == i=5 ===
        dof_indices[5]: 17
        dphi_face[5][0]: (x,y,z)=(-0.611729,        0,        0)
        _sys_solution[5->17]: 3.03524e+07
        _sys_solution[5->17] * dphi_face[5]: (x,y,z)=(-1.85674e+07,        0,        0)
        =======
   JxW_face[0]: 5
   grad_at_qn=        (x,y,z)=(-6.29947e-06, -7.7514e-06,        0)
   normals[0]= (x,y,z)=(       1,        0,       -0)
   apply_normal=      -3.14973e-05




________________________________________
From: Renato Poli <re...@gm...<mailto:re...@gm...>>
Sent: Friday, August 25, 2017 6:12 PM
To: Mike Marchywka
Cc: John Peterson; libmesh-users
Subject: Re: [Libmesh-users] Trouble with boundary integration

Hi Mike,

Thanks for the reply.
Any hint on how to easily reproduce the shape function calculations by hand (that is, manually producing the dphi_face values)?

On Fri, Aug 25, 2017 at 6:17 PM, Mike Marchywka <mar...@ho...<mailto:mar...@ho...><mailto:mar...@ho...<mailto:mar...@ho...>>> wrote:
I've been trying to learn this myself and encounter a variety of issues like this. If you have
an exact solution, why not just print all the pieces- values and normals etc-  and see how they differ
point by point instead of
trying to do an inverse problem of guessing what is wrong with an integral ?

Well - that is what I am trying to do.
I isolated a single element and, as the problem is symmetrical, it should produce a non-zero result.
There is not a lot of freedom in the calculation.
It should be a simple multiplication of a bunch of numbers...

Varying the mesh size seems to be useful sometimes.  How do the element sizes
compare to the radius of curvature ? Can you test the code on a simpler geometry such
as between two planes?

Agreed. I spent two days on that, varying all sort of stuff.
I tried from very small elements, to huge ones.
I think my problem might be related to the order of the shape function or the gauss quadrature.
Really not sure.
Any hint on how to easily reproduce the shape function calculations by hand (that is, manually producing the dphi_face values)?

AFAICT in the general
case doing integrals of normal components over surfaces is not always easy
with p and h refinements etc. Although maybe someone can clarify a bit.

I am using CGAL to generate and refine the mesh. It is quite flexible.
I am really considering that I have a conceptual issue.
It seems there is something really close to me that is fooling me ...

I will keep trying.

Please let me know if anyone has any suggestion ...

Thanks,
Renato


Thanks.

________________________________________
From: Renato Poli <re...@gm...<mailto:re...@gm...><mailto:re...@gm...<mailto:re...@gm...>>>
Sent: Friday, August 25, 2017 3:02 PM
To: John Peterson
Cc: libmesh-users
Subject: Re: [Libmesh-users] Trouble with boundary integration

Hi John,

I am still struggling and running out of ideas here.
I wrote down some debugging - perhaps you can find something weird.
The calculations seem correct, at least to what I could understand so far.

But things are zero'ing out.

If, in the mesh creation, I declare the variable as "FIRST" order, things
do not zero out.
The integration of the gradient only zeros out if I use:
        sys.add_variable("u", SECOND);

For face integration, I am using first order gauss quadrature:
        QGauss qface(dim-1, FIRST);

Any hint? Any debugging suggestion?

Rgds,
Renato

# Debug[1]: Elem: 4
# Debug[1]:     (x,y,z)=(       5,       10,        0)
# Debug[1]:     (x,y,z)=(       0,       10,        0)
# Debug[1]:     (x,y,z)=( 3.08169,  6.73058,        0)
# Debug[1]:     (x,y,z)=(     2.5,       10,        0)
# Debug[1]:     (x,y,z)=( 1.54085,  8.36529,        0)
# Debug[1]:     (x,y,z)=( 4.04085,  8.36529,        0)
# Debug[1]: num_dofs=6
# Debug[1]:     Side: 0
# Debug[1]:     Boundary: 1
# Debug[1]:          qface_point = (x,y,z)=(     2.5,       10,        0)
# Debug[1]:             == i=0 ===
# Debug[1]:             dof_indices[0]: 23
# Debug[1]:             dphi_face[0][0]: (x,y,z)=(     0.2, 0.188516,
-0)
# Debug[1]:             _sys_solution[0->23]: 3e+07
# Debug[1]:             _sys_solution[0->23] * dphi_face[0]: (x,y,z)=(
6e+06, 5.65548e+06,       -0)
# Debug[1]:             == i=1 ===
# Debug[1]:             dof_indices[1]: 24
# Debug[1]:             dphi_face[1][0]: (x,y,z)=(    -0.2,
0.117348,        0)
# Debug[1]:             _sys_solution[1->24]: 3e+07
# Debug[1]:             _sys_solution[1->24] * dphi_face[1]: (x,y,z)=(
-6e+06, 3.52045e+06,        0)
# Debug[1]:             == i=2 ===
# Debug[1]:             dof_indices[2]: 25
# Debug[1]:             dphi_face[2][0]: (x,y,z)=(       0, 0.305864,
-0)
# Debug[1]:             _sys_solution[2->25]: 3.12065e+07
# Debug[1]:             _sys_solution[2->25] * dphi_face[2]:
(x,y,z)=(       0, 9.54495e+06,       -0)
# Debug[1]:             == i=3 ===
# Debug[1]:             dof_indices[3]: 26
# Debug[1]:             dphi_face[3][0]: (x,y,z)=(      -0,
0.611729,        0)
# Debug[1]:             _sys_solution[3->26]: 3e+07
# Debug[1]:             _sys_solution[3->26] * dphi_face[3]: (x,y,z)=(
-0, 1.83519e+07,        0)
# Debug[1]:             == i=4 ===
# Debug[1]:             dof_indices[4]: 27
# Debug[1]:             dphi_face[4][0]: (x,y,z)=(       0,
-0.611729,        0)
# Debug[1]:             _sys_solution[4->27]: 3.02508e+07
# Debug[1]:             _sys_solution[4->27] * dphi_face[4]:
(x,y,z)=(       0, -1.85053e+07,        0)
# Debug[1]:             == i=5 ===
# Debug[1]:             dof_indices[5]: 28
# Debug[1]:             dphi_face[5][0]: (x,y,z)=(       0,
-0.611729,        0)
# Debug[1]:             _sys_solution[5->28]: 3.03524e+07
# Debug[1]:             _sys_solution[5->28] * dphi_face[5]:
(x,y,z)=(       0, -1.85674e+07,        0)
# Debug[1]:             =======
# Debug[1]:        JxW_face[0]: 5
# Debug[1]:        grad_at_qn=        (x,y,z)=(8.95187e-06,
-6.02752e-06,        0)
# Debug[1]:        normals[0]= (x,y,z)=(       0,        1,       -0)
# Debug[1]:        JxW_face[0] * grad_at_qn * normals[0] =      -3.01376e-05


On Fri, Aug 25, 2017 at 2:25 PM, Renato Poli <re...@gm...<mailto:re...@gm...><mailto:re...@gm...<mailto:re...@gm...>>> wrote:

> Hi John,
>
> _sys_solution is the system->solution vector.
>
> It is fed by:
>         _system.solution->localize_to_one( _sys_solution );
>
>
> On Fri, Aug 25, 2017 at 1:48 PM, John Peterson <jwp...@gm...<mailto:jwp...@gm...><mailto:jwp...@gm...<mailto:jwp...@gm...>>>
> wrote:
>
>>
>>
>> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli <re...@gm...<mailto:re...@gm...><mailto:re...@gm...<mailto:re...@gm...>>> wrote:
>>
>>> Hi,
>>>
>>> I need some help here ...
>>>
>>> I am having trouble to understand what is going on with a boundary
>>> integral.
>>> I have solved a Poisson problem for pressure in a 2D mesh, as a model
>>> for a
>>> incompressible water flow in porous media.
>>> The mesh (Tri6 elements) has a near-circular polygon in the middle
>>> representing an injection well.
>>>
>>> RHS of the equations are zero except in boundaries dofs. Dirichlet
>>> boundary
>>> conditions are imposed using a penalty both on the well and on the
>>> external
>>> boundaries to represent fixed pressure. (the solution looks ok - seems to
>>> work)
>>>
>>> To validate the solution, I am trying to estimate the water flow both in
>>> the well and in the external boundaries (they must be identical and match
>>> the analytical solution). I therefore calculate the normal pressure
>>> gradient and integrate in the whole boundary path.
>>>
>>> When I use first order variable, the results are "almost correct". But
>>> when
>>> I use second order , it is absolutely off.
>>>
>>> I am using the code below to do the integration. I tried many
>>> alternatives,
>>> no success.
>>>
>>> Any idea of what I am doing wrong?
>>>
>>
>> The code seems OK to me, although you should not need to call
>> fe_face->reinit() inside the loop over boundary ids.
>>
>> What is _sys_solution?
>>
>> --
>> John
>>
>
>
------------------------------------------------------------------------------
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...><mailto:Lib...@li...<mailto:Lib...@li...>>
https://lists.sourceforge.net/lists/listinfo/libmesh-users

Re: [Libmesh-users] Trouble with boundary integration

[Renato P. <re...@gm...>]

On Fri, Aug 25, 2017 at 8:08 PM, Mike Marchywka <mar...@ho...>
wrote:

> >Any hint on how to easily reproduce the shape function calculations by
> hand (that is, manually producing the dphi_face values)?
> Very carefully :) I ended up writing some polynomial classes and doing the
> integrals and derivatives for QUAD4 that way. With my
> rational number class it was easy to get exact results and make generated
> code for the simple things like Laplacian of rectangle.
>
That seems a good idea. I am using TRI6.
The results I am getting are so off the expected that I won't matter
getting approximations for now as reference.
This is the reason why I think I got a conceptual issue. It is really far
from target.


> I guess you could probably do the same with MatLab or Mathematica or maybe
> R.  I saw the point by point quadrature
> thing in the examples but AFAICT in many cases you can do it exactly and
> just skip a lot of that. The Laplacian
> is just an integral of a second derivative which would be zero for linear
> shape functions except for integration by
> parts which then gives you an easily computable integral. I'm using all
> analytical integrals now- just putting numbers
> into the system matrix and rhs with QUAD4 elements- but have a possible
> model problem that probably needs point by point  integration :)
>
> But do the derivatives look right? The normals? The solution? Does it work
> without AMR ? p? h? etc...
>
The normals look right. I have simplified the mesh to have right angles and
thus simple normals (1,0,0) or (-1,0,0) etc.
I am not using AMR at all. Mesh is dynamically generated with a couple of
parameters (maximum triangle size and internal angles, for example).


> What does work, if anything?
>
Everything works in a wrong way. :-)

Look this example:

+++ Triangle coordinates (each dof)
[0]    (x,y,z)=(      10,        5,        0)
[1]    (x,y,z)=(      10,       10,        0)
[2]    (x,y,z)=( 6.73058,  6.91831,        0)
[3]    (x,y,z)=(      10,      7.5,        0)
[4]    (x,y,z)=( 8.36529,  8.45915,        0)
[5]    (x,y,z)=( 8.36529,  5.95915,        0)

See the highlighted system solutions for dof[1] and dof[2].
It shows a delta of 0.1206E7 in (dx=10-6.73=3.27).
Something like grad_x=0.0369E7.
I dont get why the shape function and the summation sort of kills this
gradient...
Still struggling here...

+++ intermediate calculations
     qface_point = (x,y,z)=(      10,      7.5,        0)
        == i=0 ===
        dof_indices[0]: 12
        dphi_face[0][0]: (x,y,z)=(0.188516,     -0.2,       -0)
        _sys_solution[0->12]: 3e+07
        _sys_solution[0->12] * dphi_face[0]: (x,y,z)=(5.65548e+06,
-6e+06,       -0)
        == i=1 ===
        dof_indices[1]: 13
        dphi_face[1][0]: (x,y,z)=(0.117348,      0.2,        0)
        _sys_solution[1->13]: 3e+07
/// DOF[1] => 3E7 @ (10,5) /// <===
        _sys_solution[1->13] * dphi_face[1]: (x,y,z)=(3.52045e+06,
6e+06,        0)
        == i=2 ===
        dof_indices[2]: 14
        dphi_face[2][0]: (x,y,z)=(0.305864,       -0,       -0)
        _sys_solution[2->14]: 3.12065e+07
/// DOF[2] => 3.12E7 @ (6.73,  6.91) /// <===
        _sys_solution[2->14] * dphi_face[2]: (x,y,z)=(9.54495e+06,
-0,       -0)
        == i=3 ===
        dof_indices[3]: 15
        dphi_face[3][0]: (x,y,z)=(0.611729,        0,        0)
        _sys_solution[3->15]: 3e+07
        _sys_solution[3->15] * dphi_face[3]: (x,y,z)=(1.83519e+07,
0,        0)
        == i=4 ===
        dof_indices[4]: 16
        dphi_face[4][0]: (x,y,z)=(-0.611729,        0,        0)
        _sys_solution[4->16]: 3.02508e+07
        _sys_solution[4->16] * dphi_face[4]: (x,y,z)=(-1.85053e+07,
0,        0)
        == i=5 ===
        dof_indices[5]: 17
        dphi_face[5][0]: (x,y,z)=(-0.611729,        0,        0)
        _sys_solution[5->17]: 3.03524e+07
        _sys_solution[5->17] * dphi_face[5]: (x,y,z)=(-1.85674e+07,
0,        0)
        =======
   JxW_face[0]: 5
   grad_at_qn=        (x,y,z)=(-6.29947e-06, -7.7514e-06,        0)
   normals[0]= (x,y,z)=(       1,        0,       -0)
   apply_normal=      -3.14973e-05



>
>
> ________________________________________
> From: Renato Poli <re...@gm...>
> Sent: Friday, August 25, 2017 6:12 PM
> To: Mike Marchywka
> Cc: John Peterson; libmesh-users
> Subject: Re: [Libmesh-users] Trouble with boundary integration
>
> Hi Mike,
>
> Thanks for the reply.
> Any hint on how to easily reproduce the shape function calculations by
> hand (that is, manually producing the dphi_face values)?
>
> On Fri, Aug 25, 2017 at 6:17 PM, Mike Marchywka <mar...@ho...
> <mailto:mar...@ho...>> wrote:
> I've been trying to learn this myself and encounter a variety of issues
> like this. If you have
> an exact solution, why not just print all the pieces- values and normals
> etc-  and see how they differ
> point by point instead of
> trying to do an inverse problem of guessing what is wrong with an integral
> ?
>
> Well - that is what I am trying to do.
> I isolated a single element and, as the problem is symmetrical, it should
> produce a non-zero result.
> There is not a lot of freedom in the calculation.
> It should be a simple multiplication of a bunch of numbers...
>
> Varying the mesh size seems to be useful sometimes.  How do the element
> sizes
> compare to the radius of curvature ? Can you test the code on a simpler
> geometry such
> as between two planes?
>
> Agreed. I spent two days on that, varying all sort of stuff.
> I tried from very small elements, to huge ones.
> I think my problem might be related to the order of the shape function or
> the gauss quadrature.
> Really not sure.
> Any hint on how to easily reproduce the shape function calculations by
> hand (that is, manually producing the dphi_face values)?
>
> AFAICT in the general
> case doing integrals of normal components over surfaces is not always easy
> with p and h refinements etc. Although maybe someone can clarify a bit.
>
> I am using CGAL to generate and refine the mesh. It is quite flexible.
> I am really considering that I have a conceptual issue.
> It seems there is something really close to me that is fooling me ...
>
> I will keep trying.
>
> Please let me know if anyone has any suggestion ...
>
> Thanks,
> Renato
>
>
> Thanks.
>
> ________________________________________
> From: Renato Poli <re...@gm...<mailto:re...@gm...>>
> Sent: Friday, August 25, 2017 3:02 PM
> To: John Peterson
> Cc: libmesh-users
> Subject: Re: [Libmesh-users] Trouble with boundary integration
>
> Hi John,
>
> I am still struggling and running out of ideas here.
> I wrote down some debugging - perhaps you can find something weird.
> The calculations seem correct, at least to what I could understand so far.
>
> But things are zero'ing out.
>
> If, in the mesh creation, I declare the variable as "FIRST" order, things
> do not zero out.
> The integration of the gradient only zeros out if I use:
>         sys.add_variable("u", SECOND);
>
> For face integration, I am using first order gauss quadrature:
>         QGauss qface(dim-1, FIRST);
>
> Any hint? Any debugging suggestion?
>
> Rgds,
> Renato
>
> # Debug[1]: Elem: 4
> # Debug[1]:     (x,y,z)=(       5,       10,        0)
> # Debug[1]:     (x,y,z)=(       0,       10,        0)
> # Debug[1]:     (x,y,z)=( 3.08169,  6.73058,        0)
> # Debug[1]:     (x,y,z)=(     2.5,       10,        0)
> # Debug[1]:     (x,y,z)=( 1.54085,  8.36529,        0)
> # Debug[1]:     (x,y,z)=( 4.04085,  8.36529,        0)
> # Debug[1]: num_dofs=6
> # Debug[1]:     Side: 0
> # Debug[1]:     Boundary: 1
> # Debug[1]:          qface_point = (x,y,z)=(     2.5,       10,        0)
> # Debug[1]:             == i=0 ===
> # Debug[1]:             dof_indices[0]: 23
> # Debug[1]:             dphi_face[0][0]: (x,y,z)=(     0.2, 0.188516,
> -0)
> # Debug[1]:             _sys_solution[0->23]: 3e+07
> # Debug[1]:             _sys_solution[0->23] * dphi_face[0]: (x,y,z)=(
> 6e+06, 5.65548e+06,       -0)
> # Debug[1]:             == i=1 ===
> # Debug[1]:             dof_indices[1]: 24
> # Debug[1]:             dphi_face[1][0]: (x,y,z)=(    -0.2,
> 0.117348,        0)
> # Debug[1]:             _sys_solution[1->24]: 3e+07
> # Debug[1]:             _sys_solution[1->24] * dphi_face[1]: (x,y,z)=(
> -6e+06, 3.52045e+06,        0)
> # Debug[1]:             == i=2 ===
> # Debug[1]:             dof_indices[2]: 25
> # Debug[1]:             dphi_face[2][0]: (x,y,z)=(       0, 0.305864,
> -0)
> # Debug[1]:             _sys_solution[2->25]: 3.12065e+07
> # Debug[1]:             _sys_solution[2->25] * dphi_face[2]:
> (x,y,z)=(       0, 9.54495e+06,       -0)
> # Debug[1]:             == i=3 ===
> # Debug[1]:             dof_indices[3]: 26
> # Debug[1]:             dphi_face[3][0]: (x,y,z)=(      -0,
> 0.611729,        0)
> # Debug[1]:             _sys_solution[3->26]: 3e+07
> # Debug[1]:             _sys_solution[3->26] * dphi_face[3]: (x,y,z)=(
> -0, 1.83519e+07,        0)
> # Debug[1]:             == i=4 ===
> # Debug[1]:             dof_indices[4]: 27
> # Debug[1]:             dphi_face[4][0]: (x,y,z)=(       0,
> -0.611729,        0)
> # Debug[1]:             _sys_solution[4->27]: 3.02508e+07
> # Debug[1]:             _sys_solution[4->27] * dphi_face[4]:
> (x,y,z)=(       0, -1.85053e+07,        0)
> # Debug[1]:             == i=5 ===
> # Debug[1]:             dof_indices[5]: 28
> # Debug[1]:             dphi_face[5][0]: (x,y,z)=(       0,
> -0.611729,        0)
> # Debug[1]:             _sys_solution[5->28]: 3.03524e+07
> # Debug[1]:             _sys_solution[5->28] * dphi_face[5]:
> (x,y,z)=(       0, -1.85674e+07,        0)
> # Debug[1]:             =======
> # Debug[1]:        JxW_face[0]: 5
> # Debug[1]:        grad_at_qn=        (x,y,z)=(8.95187e-06,
> -6.02752e-06,        0)
> # Debug[1]:        normals[0]= (x,y,z)=(       0,        1,       -0)
> # Debug[1]:        JxW_face[0] * grad_at_qn * normals[0] =
> -3.01376e-05
>
>
> On Fri, Aug 25, 2017 at 2:25 PM, Renato Poli <re...@gm...<mailto:
> re...@gm...>> wrote:
>
> > Hi John,
> >
> > _sys_solution is the system->solution vector.
> >
> > It is fed by:
> >         _system.solution->localize_to_one( _sys_solution );
> >
> >
> > On Fri, Aug 25, 2017 at 1:48 PM, John Peterson <jwp...@gm...
> <mailto:jwp...@gm...>>
> > wrote:
> >
> >>
> >>
> >> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli <re...@gm...
> <mailto:re...@gm...>> wrote:
> >>
> >>> Hi,
> >>>
> >>> I need some help here ...
> >>>
> >>> I am having trouble to understand what is going on with a boundary
> >>> integral.
> >>> I have solved a Poisson problem for pressure in a 2D mesh, as a model
> >>> for a
> >>> incompressible water flow in porous media.
> >>> The mesh (Tri6 elements) has a near-circular polygon in the middle
> >>> representing an injection well.
> >>>
> >>> RHS of the equations are zero except in boundaries dofs. Dirichlet
> >>> boundary
> >>> conditions are imposed using a penalty both on the well and on the
> >>> external
> >>> boundaries to represent fixed pressure. (the solution looks ok - seems
> to
> >>> work)
> >>>
> >>> To validate the solution, I am trying to estimate the water flow both
> in
> >>> the well and in the external boundaries (they must be identical and
> match
> >>> the analytical solution). I therefore calculate the normal pressure
> >>> gradient and integrate in the whole boundary path.
> >>>
> >>> When I use first order variable, the results are "almost correct". But
> >>> when
> >>> I use second order , it is absolutely off.
> >>>
> >>> I am using the code below to do the integration. I tried many
> >>> alternatives,
> >>> no success.
> >>>
> >>> Any idea of what I am doing wrong?
> >>>
> >>
> >> The code seems OK to me, although you should not need to call
> >> fe_face->reinit() inside the loop over boundary ids.
> >>
> >> What is _sys_solution?
> >>
> >> --
> >> John
> >>
> >
> >
> ------------------------------------------------------------
> ------------------
> 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:Libmesh
> -u...@li...>
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
>