libmesh-users — General information for libmesh users

[Libmesh-users] Creating a Petsc matrix and vector

[Prashant K. J. <pjh...@gm...>]

Hi All,

I was wondering if there is a way in LibMesh to create a petsc vector and
matrix given the number of DoFs and sparsity? I have a 1D-3D coupled
problem and to create a resulting block matrix, I am relying on GMM
library. It would be nice if I can create a block matrix (one block for 3D
field in LibMesh and another for 1D field on 1D network) in Petsc and use
its parallel computation feature. Any suggestions are welcome.

Thank you.

Regards,
Prashant

Re: [Libmesh-users] (no subject)

[Prashant K. J. <pjh...@gm...>]

Hi All,

I figured it out. I was not calling localize from all the processor and
that could have been the problem. Nevermind the question.

Regards,
Prashant

On Sat, Jul 11, 2020 at 9:02 PM Prashant K. Jha <pjh...@gm...> wrote:

> Correction: For approach 1 and 2, I get error
>
> Assertion `this->verify(r.size())' failed.
>
> On Sat, Jul 11, 2020 at 8:55 PM Prashant K. Jha <pjh...@gm...>
> wrote:
>
> > Hi All,
> >
> > I have been struggling with a very simple problem. I create a system as
> > follows:
> >
> > EquationSystems es(mesh);
> > auto &taf_sys = es.add_system<TransientLinearImplicitSystem>("TAF");
> > taf_sys.add_variable("taf", FIRST);
> >
> > After I solve taf_sys, I want to localize the solution of taf_sys system
> > (ie "taf" field) on all processors in the taf_localize vector. I have
> tried
> > many approaches such as
> >
> > 1.
> > std::vector<double> taf_localize;
> > taf_sys.solution->localize(taf_localize);
> >
> > 2.
> > std::vector<double> taf_localize;
> > taf_sys.current_local_solution->localize(taf_localize);
> >
> > 3.
> > static int init = -1;
> > static std::unique_ptr<NumericVector<Number>> taf_localize =
> >         NumericVector<Number>::build(comm());
> > if (init < 0) {
> >   taf_localize->init (taf_sys.solution->size(), false, SERIAL);
> >   init = 0;
> > }
> >
> > taf_sys.solution->localize(*taf_localize);
> >
> > They all fail with
> >
> > Assertion `(this->comm()).verify(std::to_string(642))' failed
> >
> > Can anyone help with this? Thank you.
> >
> > Regards,
> > Prashant
> >
> >
>
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

Re: [Libmesh-users] (no subject)

[Prashant K. J. <pjh...@gm...>]

Correction: For approach 1 and 2, I get error

Assertion `this->verify(r.size())' failed.

On Sat, Jul 11, 2020 at 8:55 PM Prashant K. Jha <pjh...@gm...> wrote:

> Hi All,
>
> I have been struggling with a very simple problem. I create a system as
> follows:
>
> EquationSystems es(mesh);
> auto &taf_sys = es.add_system<TransientLinearImplicitSystem>("TAF");
> taf_sys.add_variable("taf", FIRST);
>
> After I solve taf_sys, I want to localize the solution of taf_sys system
> (ie "taf" field) on all processors in the taf_localize vector. I have tried
> many approaches such as
>
> 1.
> std::vector<double> taf_localize;
> taf_sys.solution->localize(taf_localize);
>
> 2.
> std::vector<double> taf_localize;
> taf_sys.current_local_solution->localize(taf_localize);
>
> 3.
> static int init = -1;
> static std::unique_ptr<NumericVector<Number>> taf_localize =
>         NumericVector<Number>::build(comm());
> if (init < 0) {
>   taf_localize->init (taf_sys.solution->size(), false, SERIAL);
>   init = 0;
> }
>
> taf_sys.solution->localize(*taf_localize);
>
> They all fail with
>
> Assertion `(this->comm()).verify(std::to_string(642))' failed
>
> Can anyone help with this? Thank you.
>
> Regards,
> Prashant
>
>

[Libmesh-users] (no subject)

[Prashant K. J. <pjh...@gm...>]

Hi All,

I have been struggling with a very simple problem. I create a system as
follows:

EquationSystems es(mesh);
auto &taf_sys = es.add_system<TransientLinearImplicitSystem>("TAF");
taf_sys.add_variable("taf", FIRST);

After I solve taf_sys, I want to localize the solution of taf_sys system
(ie "taf" field) on all processors in the taf_localize vector. I have tried
many approaches such as

1.
std::vector<double> taf_localize;
taf_sys.solution->localize(taf_localize);

2.
std::vector<double> taf_localize;
taf_sys.current_local_solution->localize(taf_localize);

3.
static int init = -1;
static std::unique_ptr<NumericVector<Number>> taf_localize =
        NumericVector<Number>::build(comm());
if (init < 0) {
  taf_localize->init (taf_sys.solution->size(), false, SERIAL);
  init = 0;
}

taf_sys.solution->localize(*taf_localize);

They all fail with

Assertion `(this->comm()).verify(std::to_string(642))' failed

Can anyone help with this? Thank you.

Regards,
Prashant

Re: [Libmesh-users] Poisson w/ Penalty Dirichlet Example

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

On Thu, Jul 9, 2020 at 11:06 AM Caleb M Phillips <cal...@ut...>
wrote:

> John,
>
> I wasn't too sure how to accomplish that, but here is the main file. This
> is a module in a larger code, but the function I'm looking at is
> "assemble_pressure".
> https://github.com/CalebPhillips5/ABM_Ang/blob/master/pressure_ABM . In
> this instance I'm just trying to set the top boundary to 1 and the bottom
> boundary to 0 with 2nd order elements.
>

Looking quickly at the assembly code, it seems like it should be fine,
although I could be missing something subtle.

Are you hitting the "build_cube()" code path for the Mesh? If so, you will
need to change the element type to TRI6 in order to use SECOND-order FEs.
Similarly for any meshes you read in... one way to achieve this is by
generating/reading in a mesh with linear elements and then calling
all_second_order() to convert it.

If you aren't already, I'd also strongly suggest compiling and running in
dbg mode, as I believe an error like this should be caught by an assert in
dbg mode.

--
John

Re: [Libmesh-users] Poisson w/ Penalty Dirichlet Example

[Caleb M P. <cal...@ut...>]

John,

I wasn't too sure how to accomplish that, but here is the main file. This
is a module in a larger code, but the function I'm looking at is
"assemble_pressure".
https://github.com/CalebPhillips5/ABM_Ang/blob/master/pressure_ABM . In
this instance I'm just trying to set the top boundary to 1 and the bottom
boundary to 0 with 2nd order elements.

Cheers,
Caleb

On Wed, Jul 8, 2020 at 5:10 PM John Peterson <jwp...@gm...> wrote:

>
>
> On Wed, Jul 8, 2020 at 4:48 PM Caleb M Phillips <cal...@ut...>
> wrote:
>
>> Hey all,
>>
>> Quick question that is completely throwing me for a loop. I'm working on a
>> poisson equation for pressure and wanting to specify dirichlet boundary
>> conditions. I used the penalty method (per libmesh documentation), and it
>> works beautifully, aka matches the real solution. However, when I change
>> to
>> 2nd order polynomials, everything falls apart. The penalty method does
>> implement the boundary conditions, but the solution is complete garbage.
>>
>> Any explanation for changing (.., FIRST) to (.., SECOND) and the solution
>> (the code runs fine) but the solution being garbage?
>>
>
>
> I just ran introduction_ex3, which uses penalty BCs and second-order
> elements, and wrote the result to an Exodus file, and it looks OK to me.
>
> So my guess is it's a bug of some sort in your code. Is it possible for
> you to share it somehow?
>
> --
> John
>
>

Re: [Libmesh-users] Poisson w/ Penalty Dirichlet Example

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

On Wed, Jul 8, 2020 at 4:48 PM Caleb M Phillips <cal...@ut...>
wrote:

> Hey all,
>
> Quick question that is completely throwing me for a loop. I'm working on a
> poisson equation for pressure and wanting to specify dirichlet boundary
> conditions. I used the penalty method (per libmesh documentation), and it
> works beautifully, aka matches the real solution. However, when I change to
> 2nd order polynomials, everything falls apart. The penalty method does
> implement the boundary conditions, but the solution is complete garbage.
>
> Any explanation for changing (.., FIRST) to (.., SECOND) and the solution
> (the code runs fine) but the solution being garbage?
>


I just ran introduction_ex3, which uses penalty BCs and second-order
elements, and wrote the result to an Exodus file, and it looks OK to me.

So my guess is it's a bug of some sort in your code. Is it possible for you
to share it somehow?

-- 
John

[Libmesh-users] Poisson w/ Penalty Dirichlet Example

[Caleb M P. <cal...@ut...>]

Hey all,

Quick question that is completely throwing me for a loop. I'm working on a
poisson equation for pressure and wanting to specify dirichlet boundary
conditions. I used the penalty method (per libmesh documentation), and it
works beautifully, aka matches the real solution. However, when I change to
2nd order polynomials, everything falls apart. The penalty method does
implement the boundary conditions, but the solution is complete garbage.

Any explanation for changing (.., FIRST) to (.., SECOND) and the solution
(the code runs fine) but the solution being garbage?

Cheers,
Caleb Phillips

Re: [Libmesh-users] Create an unstructured mesh of map of country/region

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

On Sat, 6 Jun 2020, Prashant K. Jha wrote:

> Thank you so much for your time in replying to my email. I have been able to create a Delaunay triangulation of the map of the Texas state. How I
> did it is as follows:
> 
> 1. Get the shapefile for the map from http://gis-txdot.opendata.arcgis.com/datasets/texas-state-boundary
> 2. Load this into QGIS application which is free and runs across different platforms.
> 3. Use tool 'Vector->Geometry Tools->Simplify' to get coarse lines. You can control how coarse you want the outer lines.
> 4. Use tool 'Vector->Geometry Tools->Extract Vertices' to create a layer with vertices of the simplified map you obtained in step 3.
> 5. Save the layer corresponding to vertices by 'Layer->Save as'. In the save as option, you can deselect most of the features. Look for option
> 'GEOMETRY' and select 'As_XY' to store the (x,y) coordinates of the vertices. This link was
> helpful: https://gis.stackexchange.com/questions/8844/getting-list-of-coordinates-for-points-in-layer-using-qgis
> 6. Finally, I created the .geo file by reading the .csv file in step 5. The file generated in step 5 has the first and last vertices equal so
> while creating the .geo file I excluded the last vertex. I then ran gmsh on .geo file. 

This is pretty great; thanks for the update!

> I think this is not very efficient but for now it works.

Sounds about as efficient as you could hope for for a one-off use
case.  If you were trying to make meshes for dozens of other shapes
too then at that point I'd look into scripting it.
---
Roy

Re: [Libmesh-users] Create an unstructured mesh of map of country/region

[Prashant K. J. <pjh...@gm...>]

Hello Vikram and Roy,

Thank you so much for your time in replying to my email. I have been able
to create a Delaunay triangulation of the map of the Texas state. How I did
it is as follows:

1. Get the shapefile for the map from
http://gis-txdot.opendata.arcgis.com/datasets/texas-state-boundary
2. Load this into QGIS application which is free and runs across different
platforms.
3. Use tool 'Vector->Geometry Tools->Simplify' to get coarse lines. You can
control how coarse you want the outer lines.
4. Use tool 'Vector->Geometry Tools->Extract Vertices' to create a layer
with vertices of the simplified map you obtained in step 3.
5. Save the layer corresponding to vertices by 'Layer->Save as'. In the
save as option, you can deselect most of the features. Look for option
'GEOMETRY' and select 'As_XY' to store the (x,y) coordinates of the
vertices. This link was helpful:
https://gis.stackexchange.com/questions/8844/getting-list-of-coordinates-for-points-in-layer-using-qgis
6. Finally, I created the .geo file by reading the .csv file in step 5. The
file generated in step 5 has the first and last vertices equal so while
creating the .geo file I excluded the last vertex. I then ran gmsh on .geo
file.

I think this is not very efficient but for now it works. Alternatively, I
could directly extract the vertices from the shapefile by using some python
libraries such as basemap or some other libraries. I am guessing that
python libraries may also have features similar to step 3 using which you
can coarsen the boundary lines.

Regards,
Prashant

On Sat, May 16, 2020 at 12:39 AM Vikram Garg <vik...@gm...>
wrote:

>
> Vikram Garg
>
> vikramvgarg.github.io/
>
>
> On Fri, May 15, 2020 at 11:02 PM Roy Stogner <roy...@od...>
> wrote:
>
>>
>> On Fri, 15 May 2020, Vikram Garg wrote:
>>
>> > Once you have the boundaries defined, TRIANGLE, which uses Delaunay
>> > triangulation should be able to do this.
>>
>> This is true, but defining the boundaries in code is certainly a pain.
>> IIRC you wrote your own vectorizer as an undergrad, and it took us a
>> week to figure out that even though the meshes looked okay to the
>> naked eye, sharp jumps from pixel row to pixel row were behind the
>> crazy pressure spikes we were seeing in airfoil simulations.
>>
>>
> Yes, I was wary of that, but then we had sharp fractal like features due
> to the ice crystals we were trying to represent. We can have similar
> features on a coastal
> boundary representation, but I depending on the eventual purpose, these
> can be smoothed considerably.
>
>
>> Is there any "standard" 2D format for vector shape data?  My wife
>> uses (or when necessary converts everything to) SVG for laser cutter
>> and plotter use, but SVG has a zillion features beyond defining loops.
>> I'm not sure if there's a specific subset of SVG that's well-defined
>> enough to write a reader for.  I guess it's just XML under the hood.
>> Maybe the thing to do would be to load a vector Texas outline into
>> Inkscape, save it, and write a quick-and-dirty parser for whatever's
>> there?
>>
>
> I am not aware of a standard format. The inkscape idea could work, I
> myself used Matlab's imread, not sure if it supports svg currently.
>
>
>> ---
>> Roy
>>
>

Re: [Libmesh-users] Reduced Basis: Petsc Matrix Ignore Zero Entries

[Nikhil V. <nik...@gm...>]

Hello,

Thanks for your response.

@David
Reducing the number of Aq matrices is not really ideal for my application.
I would like to try and see if I can reduce the memory requirement in some
other way.

@Paul
Indeed PetscMatrix is used for Reduced Basis also. I went through the
PetscMatrix class documentation, and tried to insert the flag
MAT_IGNORE_ZERO_ENTRIES in the MatSetOption function. It turns out that
that flag does not work if the matrix is pre-allocated.

Since pre-allocation of the matrix is necessary for good performance, I do
not want to disable it. Is it possible to construct a sparsity pattern only
for a particular subdomain of the mesh and then attach it to a sparse
matrix? If this works, I think it will solve my problem.

Best regards,
Nikhil

On Thu, Apr 23, 2020 at 3:22 PM Paul T. Bauman <ptb...@gm...> wrote:

> Hi Nikhil,
>
> Typically, you would grab the "raw" underlying Mat object and do the PETSc
> call yourself. This would be accessible from the libMesh::PetscMatrix
> object. Typically the ImplicitSystem has the system matrix, but I'm not
> familiar with where these matrices would be cached on the RB side of
> things. I would use the Oxygen documentation as a starting point:
> libmesh.github.io (although it seems GitHub is misbehaving ATM). I hope
> that helps.
>
> Best,
>
> Paul
>
> On Thu, Apr 23, 2020 at 12:06 AM Nikhil Vaidya <nik...@gm...>
> wrote:
>
>> Hi,
>>
>> I am using reduced basis for my work. The number of Aq matrices in my
>> problem is very large (~250). Because of this, the memory requirements of
>> my program are ~200 times the mesh-file size. I did some memory usage
>> study
>> and found out that the allocate_data_structures() function in
>> RBConstruction is responsible for most of the memory usage.
>>
>> I think this might be due to the large number of Petsc sparse matrices
>> being constructed. Each Aq matrix is known to have non-zero entries
>> corresponding to only a small sub-domain of the geometry. All these
>> sub-domains are separate blocks in the mesh. I get the feeling that for
>> each sparse matrix memory is being allocated even for the zero entries. I
>> saw that Petsc has an option MAT_IGNORE_ZERO_ENTRIES. How does one pass
>> this in libMesh? Please note that I am not using libMesh directly. My code
>> is a MOOSE app, and I am using the RB functionality of the underlying
>> libmesh code.
>>
>> Best regards,
>> Nikhil
>>
>> _______________________________________________
>> Libmesh-users mailing list
>> Lib...@li...
>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>
>

Re: [Libmesh-users] A question for the ordering of unknowns

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

On Mon, 25 May 2020, Li Luo wrote:

> For a multicomponent system, such as Navier-Stokes equations, I found that
> the ordering of unknowns in the system depends on the order of polynomials
> used for the variables. Take a two-component system with variables u and p
> for example:
>   -If using first order for both variables, the ordering is point by
> point, that is u_1, p_1, u_2, p_2, ...., u_n, p_n.
>   -If using second order for u and first order for p, the ordering is
> field by field, i.e., all the unknowns of u are ordered first then followed
> by all the unknowns of p.
>
> For the second case, is there any way to change to a point by point
> ordering? Though in such case there is no unknown of p on the mid-edge
> point.

Try setting --node-major-dofs on the command line.  It should give you
a point by point ordering, though as you say there's no way in that
case for the number of variables to be *identical* on each point.
---
Roy

[Libmesh-users] A question for the ordering of unknowns

[Li L. <li...@ka...>]

Dear developers,

I have a question for the ordering of unknowns.

For a multicomponent system, such as Navier-Stokes equations, I found that
the ordering of unknowns in the system depends on the order of polynomials
used for the variables. Take a two-component system with variables u and p
for example:
   -If using first order for both variables, the ordering is point by
point, that is u_1, p_1, u_2, p_2, ...., u_n, p_n.
   -If using second order for u and first order for p, the ordering is
field by field, i.e., all the unknowns of u are ordered first then followed
by all the unknowns of p.

For the second case, is there any way to change to a point by point
ordering? Though in such case there is no unknown of p on the mid-edge
point.

Regards,
Li

-- 

This message and its contents, including attachments are intended solely 
for the original recipient. If you are not the intended recipient or have 
received this message in error, please notify me immediately and delete 
this message from your computer system. Any unauthorized use or 
distribution is prohibited. Please consider the environment before printing 
this email.

Re: [Libmesh-users] augment sparsity

[Gary Hu <hug...@gm...>]

Thanks John,

I modified the operator() and it works now:

void
AugmentSparsityToDense::operator()(const MeshBase::const_element_iterator &
range_begin,
                                   const MeshBase::const_element_iterator &
range_end,
                                   processor_id_type p,
                                   map_type & coupled_elements)
{
  const CouplingMatrix * const null_mat = nullptr;
  for (const auto & elem : as_range(range_begin, range_end))
  {
    coupled_elements.emplace(elem, null_mat);
    for (const auto & elem_remote : _mesh.active_element_ptr_range())
      coupled_elements.emplace(elem_remote, null_mat);
  }
}

I originally thought the range_begin and range_end are for all active
elements, but it turns out that (at least in my use case for debugging) the
range_begin and range_end is giving me one element at a time.
Therefore I have to explicitly loop over all the active elements to attach
the coupling matrix. I think the documentation on the operator() needs to
be improved unless I missed something there.

Now it works and it is indeed much faster compared to when I had "*total
number of mallocs used during MatSetValues calls != 0*".

Regarding the dense matrix -- I am happy to submit a PR but I need time to
go through the source code (today is my third day using libMesh).
I am not trying to implement a fast multipole method. I am just assemblying
a Galerkin projection of a Fredholm Integral Equation for spatial
covariance.
I only need to factorize the covariance matrix once and for all, so I can
live with the N^2 assembly time.
I am aware of some of the multigrid methods in literature, but I don't
really want to deal with the accuracy issues that those methods bring to me.

Thanks,
Gary

On Mon, May 18, 2020 at 9:11 AM John Peterson <jwp...@gm...> wrote:

>
>
> On Sat, May 16, 2020 at 5:12 PM Gary Hu <hug...@gm...> wrote:
>
>> Hello all,
>>
>> I am new to libMesh. I am trying to solve a generalized eigen-value
>> problem
>> for an integral equation. I tried example "Solving a generalized Eigen
>> Problem" and it works fine.
>>
>> Next, I modified the assembly routine to have a skeleton like this:
>>
>> for (const auto & elem : mesh.active_local_element_ptr_range())
>>     for (const auto & elem_remote : mesh.active_element_ptr_range())
>>
>> So basically for each local element, I am looping over all the active
>> elements (local and nonlocal) to assemble the integral equation. I am
>> pretty confident that I used fe->reinit and other housekeeping stuff
>> properly.
>>
>
>
> It sounds like you are trying to build a dense matrix... in such a case
> all the extra "bookkeeping" indices etc. that are required for sparse
> matrices is not required, so this approach won't scale well. I think you
> may want to look into building a dense matrix for your method, e.g.
> MatCreateDense in PETSc [0]. We don't have C++ interfaces for MatDense in
> libMesh yet since our main focus is FEA but it's something I think we'd
> definitely like to have a PR adding. Regarding your original question, I
> don't see anything obviously wrong with the GhostingFunctor approach you
> are taking, so it may just be a coding bug that's causing it to not work.
> I'd start by adding print statements to check whether it's actually being
> called and what elements are being added to coupled_elements. I don't know
> how much time you want to spend debugging this approach though, because
> even if it works I fear it is going to be incredibly slow, not only due to
> the storing-dense-matrix-in-sparse-format issue, but also because of the
> N^2 assembly loop you've proposed above... perhaps there is an existing way
> to speed up whatever you are doing (fast multipole method?) which you
> should look into first.
>
> [0]:
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateDense.html
>
> --
> John
>
>
>
>

Re: [Libmesh-users] augment sparsity

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

On Sat, May 16, 2020 at 5:12 PM Gary Hu <hug...@gm...> wrote:

> Hello all,
>
> I am new to libMesh. I am trying to solve a generalized eigen-value problem
> for an integral equation. I tried example "Solving a generalized Eigen
> Problem" and it works fine.
>
> Next, I modified the assembly routine to have a skeleton like this:
>
> for (const auto & elem : mesh.active_local_element_ptr_range())
>     for (const auto & elem_remote : mesh.active_element_ptr_range())
>
> So basically for each local element, I am looping over all the active
> elements (local and nonlocal) to assemble the integral equation. I am
> pretty confident that I used fe->reinit and other housekeeping stuff
> properly.
>


It sounds like you are trying to build a dense matrix... in such a case all
the extra "bookkeeping" indices etc. that are required for sparse matrices
is not required, so this approach won't scale well. I think you may want to
look into building a dense matrix for your method, e.g. MatCreateDense in
PETSc [0]. We don't have C++ interfaces for MatDense in libMesh yet since
our main focus is FEA but it's something I think we'd definitely like to
have a PR adding. Regarding your original question, I don't see anything
obviously wrong with the GhostingFunctor approach you are taking, so it may
just be a coding bug that's causing it to not work. I'd start by adding
print statements to check whether it's actually being called and what
elements are being added to coupled_elements. I don't know how much time
you want to spend debugging this approach though, because even if it works
I fear it is going to be incredibly slow, not only due to the
storing-dense-matrix-in-sparse-format issue, but also because of the N^2
assembly loop you've proposed above... perhaps there is an existing way to
speed up whatever you are doing (fast multipole method?) which you should
look into first.

[0]:
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateDense.html

--
John

[Libmesh-users] augment sparsity

[Gary Hu <hug...@gm...>]

Hello all,

I am new to libMesh. I am trying to solve a generalized eigen-value problem
for an integral equation. I tried example "Solving a generalized Eigen
Problem" and it works fine.

Next, I modified the assembly routine to have a skeleton like this:

for (const auto & elem : mesh.active_local_element_ptr_range())
    for (const auto & elem_remote : mesh.active_element_ptr_range())

So basically for each local element, I am looping over all the active
elements (local and nonlocal) to assemble the integral equation. I am
pretty confident that I used fe->reinit and other housekeeping stuff
properly.

However, as expected, I encountered the following error from PETSc:

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: New nonzero at (0,4) caused a malloc

because I didn't augment the sparsity pattern of my matrix.

Next, I added the following lines to suppress the PETSc error:

PetscMatrix<Number> * matrix_A_petsc = dynamic_cast<PetscMatrix<Number>
*>(eigen_system.matrix_A.get());
MatSetOption(matrix_A_petsc->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR,
PETSC_FALSE);

Now my calculation could run to completion, and both the assembled matrix
and the solution are correct.
However, this is undesirable because my "*total number of mallocs used
during MatSetValues calls != 0*" and therefore the assembly process is
unbelievably slow.
I think I have to augment the sparsity correctly to speed up the assembly.

Next, I followed miscellaneous ex9 to augment the sparsity pattern using
the coupling functor. I created a class called AugmentSparsityToDense,
deriving from GhostingFunctor, and its operator() looks like:

void
AugmentSparsityToDense::operator()(const MeshBase::const_element_iterator &
range_begin,
                                   const MeshBase::const_element_iterator &
range_end,
                                   processor_id_type p,
                                   map_type & coupled_elements)
{
  const CouplingMatrix * const null_mat = nullptr;
  for (const auto & elem : as_range(range_begin, range_end))
  {
    coupled_elements.emplace(elem, null_mat);
    for (const auto & elem_remote : as_range(range_begin, range_end))
      coupled_elements.emplace(elem_remote, null_mat);
  }
}

and I attached it to the DofMap using the following lines:

std::shared_ptr<AugmentSparsityToDense> augment_sparsity(new
AugmentSparsityToDense());
eigen_system.get_dof_map().add_coupling_functor(augment_sparsity);

but then I still get the PETSc error:

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: New nonzero at (0,4) caused a malloc

and when I examine the n_nz and n_oz from DofMap, the sparsity is indeed
not augmented. Can anyone give me some suggestions?

Thanks,
Gary

Re: [Libmesh-users] Create an unstructured mesh of map of country/region

[Vikram G. <vik...@gm...>]

Vikram Garg

vikramvgarg.github.io/


On Fri, May 15, 2020 at 11:02 PM Roy Stogner <roy...@od...>
wrote:

>
> On Fri, 15 May 2020, Vikram Garg wrote:
>
> > Once you have the boundaries defined, TRIANGLE, which uses Delaunay
> > triangulation should be able to do this.
>
> This is true, but defining the boundaries in code is certainly a pain.
> IIRC you wrote your own vectorizer as an undergrad, and it took us a
> week to figure out that even though the meshes looked okay to the
> naked eye, sharp jumps from pixel row to pixel row were behind the
> crazy pressure spikes we were seeing in airfoil simulations.
>
>
Yes, I was wary of that, but then we had sharp fractal like features due to
the ice crystals we were trying to represent. We can have similar features
on a coastal
boundary representation, but I depending on the eventual purpose, these can
be smoothed considerably.


> Is there any "standard" 2D format for vector shape data?  My wife
> uses (or when necessary converts everything to) SVG for laser cutter
> and plotter use, but SVG has a zillion features beyond defining loops.
> I'm not sure if there's a specific subset of SVG that's well-defined
> enough to write a reader for.  I guess it's just XML under the hood.
> Maybe the thing to do would be to load a vector Texas outline into
> Inkscape, save it, and write a quick-and-dirty parser for whatever's
> there?
>

I am not aware of a standard format. The inkscape idea could work, I myself
used Matlab's imread, not sure if it supports svg currently.


> ---
> Roy
>

Re: [Libmesh-users] Create an unstructured mesh of map of country/region

[Vikram G. <vik...@gm...>]

Once you have the boundaries defined, TRIANGLE, which uses Delaunay
triangulation should be able to do this.

Vikram Garg

vikramvgarg.github.io/


On Fri, May 15, 2020 at 5:18 PM Prashant K. Jha <pjh...@gm...> wrote:

> Hi All,
>
> This may be an unrelated question to LibMesh. I am looking for ways to
> create an unstructured mesh of the map of Texas state. If anyone is aware
> of how to do this or know specific libraries that do this, please let me
> know. Any help would be appreciated.
>
> Regards,
> Prashant
>
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

[Libmesh-users] Create an unstructured mesh of map of country/region

[Prashant K. J. <pjh...@gm...>]

Hi All,

This may be an unrelated question to LibMesh. I am looking for ways to
create an unstructured mesh of the map of Texas state. If anyone is aware
of how to do this or know specific libraries that do this, please let me
know. Any help would be appreciated.

Regards,
Prashant

Re: [Libmesh-users] make check reduced_basis_ex7 flags

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

On Tue, May 12, 2020 at 10:10 AM <ed...@op...> wrote:

> On 2020-05-12 14:18, John Peterson wrote:
> > The options still exist for people using older versions of PETSc, so
> > yes,
> > it makes sense to still test for them.
>
> Of course! Thanks! By the way, I am maintaining an AUR package
> (aur.archlinux.org) of libmesh right now.
>

OK, I don't know much about Arch but please let us know if we can help by
having some packaging metadata files in the repository or something along
those lines.

-- 
John

Re: [Libmesh-users] make check reduced_basis_ex7 flags

[<ed...@op...>]

On 2020-05-12 14:18, John Peterson wrote:
> The options still exist for people using older versions of PETSc, so 
> yes,
> it makes sense to still test for them.

Of course! Thanks! By the way, I am maintaining an AUR package 
(aur.archlinux.org) of libmesh right now.

-------------------------------------------------
This free account was provided by VFEmail.net - report spam to ab...@vf...
 
ONLY AT VFEmail! - Use our Metadata Mitigator to keep your email out of the NSA's hands!
$24.95 ONETIME Lifetime accounts with Privacy Features!  
15GB disk! No bandwidth quotas!
Commercial and Bulk Mail Options!  

Re: [Libmesh-users] make check reduced_basis_ex7 flags

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

On Mon, May 11, 2020 at 11:50 PM <ed...@op...> wrote:

>
> That is very ingenious. Noob question: Is there still a need to test for
> options that don't exist? I know there must be a reason behind it :P .
>

The options still exist for people using older versions of PETSc, so yes,
it makes sense to still test for them.

-- 
John

Re: [Libmesh-users] make check reduced_basis_ex7 flags

[<ed...@op...>]

On 2020-05-11 11:06, David Knezevic wrote:
> I think PETSc changed the option name from 
> "-pc_factor_mat_solver_package"
> to "-pc_factor_mat_solver_type", which I guess is the issue here.

Yes, thank you.

> 
> Also, you can use "-ksp_view" on the command line to get a printout of 
> the
> solver options that were used by PETSc, so you can see if it actually 
> used
> MUMPS or superlu.

Yes, thank you again.

> 
> Best,
> David


On 2020-05-11 13:47, John Peterson wrote:
> On Sun, May 10, 2020 at 10:08 PM <ed...@op...> wrote:
> 
> This message is from libMesh (BTW what is libMesh++?):

libMesh++ is a hole in my brain from Elmer's documentation, sorry (not 
blaming it on them).

>> 
> 
> And these messages are from PETSc:
> 
>> 
> 
> I think this is all working as intended, since if you look at run.sh 
> for
> reduced_basis_ex7, you can see that we run all four combinations of
> "type"/"package" and "mumps"/superlu" from the script, relying on the 
> logic
> within the example itself to skip those cases which don't apply.

That is very ingenious. Noob question: Is there still a need to test for 
options that don't exist? I know there must be a reason behind it :P .

> The examples must be run from the build directory, so
> for example to test mumps with PETSC >= 3.9 (assuming 
> $LIBMESH_SRC/build is
> the location of your build directory):
> 
> cd $LIBMESH_SRC/build/examples/reduced_basis/reduced_basis_ex7
> ./example-opt -online_mode 0 -ksp_type preonly -pc_type lu
> -pc_factor_mat_solver_type mumps

Thank you very much for the very clear instructions.

> 
> Finally, there is one last check you can do, which is to look at 
> (again,
> from the build directory) the include/libmesh_config.h file and inspect 
> the
> values of the
> 
> LIBMESH_PETSC_HAVE_MUMPS
> LIBMESH_PETSC_HAVE_SUPERLU_DIST
> 
> preprocessor defines.

I really got it with this. It seems that SuperLU is not being used, but 
MUMPS is. I will check my PETSc compilation before asking further 
questions.

All of this was very useful. Thanks, dear David and John.

-------------------------------------------------
This free account was provided by VFEmail.net - report spam to ab...@vf...
 
ONLY AT VFEmail! - Use our Metadata Mitigator to keep your email out of the NSA's hands!
$24.95 ONETIME Lifetime accounts with Privacy Features!  
15GB disk! No bandwidth quotas!
Commercial and Bulk Mail Options!  

Re: [Libmesh-users] make check reduced_basis_ex7 flags

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

On Sun, May 10, 2020 at 10:08 PM <ed...@op...> wrote:

> Hi.
>
> I want to know if libMesh++ is actually using mumps and superlu, and if
> it is not how to fix it. I just recompiled libMesh++, and I am getting
> the following message with the reduced_basis_ex7 (this is actually an
> improvement, because it used to crash at that point). Thanks!
>
> #+begin_EXAMPLE
>
>    ... snip 8<  ...
>
>    ,***************************************************************
>    ,* Running Example reduced_basis_ex7:
>    ,*   ./example-opt -online_mode 1 -pc_factor_mat_solver_package
> superlu
>    ,***************************************************************
>
>
This message is from libMesh (BTW what is libMesh++?):


>    LibMesh was configured with PETSc >= 3.9, but command-line options use
> deprecated syntax. Skipping now.
>

And these messages are from PETSc:

   WARNING! There are options you set that were not used!
>    WARNING! could be spelling mistake, etc!
>    There are 2 unused database options. They are:
>    Option left: name:-online_mode value: 1
>    Option left: name:-pc_factor_mat_solver_package value: superlu
>

I think this is all working as intended, since if you look at run.sh for
reduced_basis_ex7, you can see that we run all four combinations of
"type"/"package" and "mumps"/superlu" from the script, relying on the logic
within the example itself to skip those cases which don't apply. If you
want to test whether it's working you can run the example manually with the
different cases and inspect the output of each one rather than running the
full "make check". The examples must be run from the build directory, so
for example to test mumps with PETSC >= 3.9 (assuming $LIBMESH_SRC/build is
the location of your build directory):

cd $LIBMESH_SRC/build/examples/reduced_basis/reduced_basis_ex7
./example-opt -online_mode 0 -ksp_type preonly -pc_type lu
-pc_factor_mat_solver_type mumps

Finally, there is one last check you can do, which is to look at (again,
from the build directory) the include/libmesh_config.h file and inspect the
values of the

LIBMESH_PETSC_HAVE_MUMPS
LIBMESH_PETSC_HAVE_SUPERLU_DIST

preprocessor defines.



>    ,***************************************************************
>    ,* Done Running Example reduced_basis_ex7:
>    ,*   ./example-opt -online_mode 1 -pc_factor_mat_solver_package
> superlu
>    ,***************************************************************
> #+end_EXAMPLE
>
> I get a similar message above with mumps. I am attaching my config.log
> as a tar ball
>

This mailing list does not allow attachments.

-- 
John

Re: [Libmesh-users] make check reduced_basis_ex7 flags

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

I think PETSc changed the option name from "-pc_factor_mat_solver_package"
to "-pc_factor_mat_solver_type", which I guess is the issue here.

Also, you can use "-ksp_view" on the command line to get a printout of the
solver options that were used by PETSc, so you can see if it actually used
MUMPS or superlu.

Best,
David


On Sun, May 10, 2020 at 11:08 PM <ed...@op...> wrote:

> Hi.
>
> I want to know if libMesh++ is actually using mumps and superlu, and if
> it is not how to fix it. I just recompiled libMesh++, and I am getting
> the following message with the reduced_basis_ex7 (this is actually an
> improvement, because it used to crash at that point). Thanks!
>
> #+begin_EXAMPLE
>
>    ... snip 8<  ...
>
>    ,***************************************************************
>    ,* Running Example reduced_basis_ex7:
>    ,*   ./example-opt -online_mode 1 -pc_factor_mat_solver_package
> superlu
>    ,***************************************************************
>
>    LibMesh was configured with PETSc >= 3.9, but command-line options use
> deprecated syntax. Skipping now.
>    WARNING! There are options you set that were not used!
>    WARNING! could be spelling mistake, etc!
>    There are 2 unused database options. They are:
>    Option left: name:-online_mode value: 1
>    Option left: name:-pc_factor_mat_solver_package value: superlu
>
>    ,***************************************************************
>    ,* Done Running Example reduced_basis_ex7:
>    ,*   ./example-opt -online_mode 1 -pc_factor_mat_solver_package
> superlu
>    ,***************************************************************
> #+end_EXAMPLE
>
> I get a similar message above with mumps. I am attaching my config.log
> as a tar ball
>
> -------------------------------------------------
> This free account was provided by VFEmail.net - report spam to
> ab...@vf...
>
> ONLY AT VFEmail! - Use our Metadata Mitigator to keep your email out of
> the NSA's hands!
> $24.95 ONETIME Lifetime accounts with Privacy Features!
> 15GB disk! No bandwidth quotas!
> Commercial and Bulk Mail Options!
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>