Hi Dmitry,
Thanks for your input, very helpful. Regarding an example: Let me know if
the attached patch for systems_of_equations_ex6 is helpful for you. It's
the patch I sent through earlier + Paul's fix for setting the vector block
size. It seems to be working well: the number of iterations isn't quite
independent of mesh refinement, but it's not far off. I guess there are
some more options that can be tweaked to improve scaling/performance
further.
David
On Mon, Mar 2, 2015 at 10:14 PM, Dmitry Karpeyev <[email protected]>
wrote:
> GAMG needs the NEARNullSpace (not the best terminology), which it uses with
> smoothed aggregation to obtain a good coarse space.
> I think the best way to think about these is as the low energy modes of the
> operator. In practical terms, it is the NullSpace of the operator
> with the boundary conditions removed. Note that the same PETSc object
> (MatNullSpace) is used to encode both the NearNullSpace and the NullSpace
> of an operator (possibly adding to the confusion). There is support in
> libMesh to pass user-specified (Near)NullSpaces down to the PETSc solver,
> although it would be good to have an example.
>
> I used to get very good performance (compared to HYPRE's BoomerAMG and even
> ASM) for linear elasticity, but lately I haven't been able to reproduce
> it. In particular, one should expect essentially the same number of
> iterations with mesh refinement. The iterations do get to be more
> expensive, naturally, but with BoomerAMG they grow, as they do with ASM,
> unless overlap is aggressively increased. This might be small consolation
> for someone with a moderate problem size. We need to have better GAMG
> defaults that result in good performance out of the box.
>
> I will try to come up with a libMesh example and reproduce my results of a
> couple of years ago.
>
> Dmitry.
>
> On Mon, Mar 2, 2015 at 9:00 PM Cody Permann <[email protected]> wrote:
>
> > Missed the list on my reply.
> > ---------- Forwarded message ---------
> > From: Cody Permann <[email protected]>
> > Date: Mon, Mar 2, 2015 at 7:58 PM
> > Subject: Re: [Libmesh-users] PETSc's GAMG with systems_of_equations_ex6
> > To: Paul T. Bauman <[email protected]>, Dmitry Karpeev <
> [email protected]
> > >
> >
> >
> > Dmitry is helping us with Gamg issues. I'm CC'ing him on this thread. He
> > should have some insight into this problem for you.
> >
> > Cody
> > On Mon, Mar 2, 2015 at 7:55 PM Paul T. Bauman <[email protected]>
> wrote:
> >
> > > On Mon, Mar 2, 2015 at 9:45 PM, David Knezevic <
> > [email protected]
> > > >
> > > wrote:
> > >
> > > >
> > > > I gather that the main requirement is to specify rigid-body modes as
> > the
> > > > near-nullspace.
> > >
> > >
> > > Ha! I was, quite literally, just having a look at similar
> functionality
> > > for removing rigid rotations.
> > >
> > > I tried this in the attached diff, and I specified
> > > > "-pc_type gamg" (and I tried a few other options like -pc_gamg_type
> > agg),
> > > > but it doesn't seem to give very good convergence. I was wondering if
> > > > anyone has experience with this and could let me know if I'm missing
> > > > something?
> > > >
> > >
> > > I don't think NumericVector sets the block size. In the PetscVector
> clone
> > > for the coordinate vector, you need to call VecSetBlockSize and set it
> > to 3
> > > (for the 3D case).
> > >
> > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/
> > > VecSetBlockSize.html
> > >
> > > Petsc's null space function for rigid modes relies on the block size
> for
> > > figuring out the vectors.
> > >
> > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/
> > > MatNullSpaceCreateRigidBody.html
> > >
> > > So, in theory, you should just need to set the block size and it should
> > be
> > > much better.
> > >
> > > I'll be very curious how this turns out.
> > >
> > > Hope that helps,
> > >
> > > Paul
> > > ------------------------------------------------------------
> > > ------------------
> > > Dive into the World of Parallel Programming The Go Parallel Website,
> > > sponsored
> > > by Intel and developed in partnership with Slashdot Media, is your hub
> > for
> > > all
> > > things parallel software development, from weekly thought leadership
> > blogs
> > > to
> > > news, videos, case studies, tutorials and more. Take a look and join
> the
> > > conversation now. http://goparallel.sourceforge.net/
> > > _______________________________________________
> > > Libmesh-users mailing list
> > > [email protected]
> > > https://lists.sourceforge.net/lists/listinfo/libmesh-users
> > >
> > ------------------------------------------------------------
> > ------------------
> > Dive into the World of Parallel Programming The Go Parallel Website,
> > sponsored
> > by Intel and developed in partnership with Slashdot Media, is your hub
> for
> > all
> > things parallel software development, from weekly thought leadership
> blogs
> > to
> > news, videos, case studies, tutorials and more. Take a look and join the
> > conversation now. http://goparallel.sourceforge.net/
> > _______________________________________________
> > Libmesh-users mailing list
> > [email protected]
> > https://lists.sourceforge.net/lists/listinfo/libmesh-users
> >
>
> ------------------------------------------------------------------------------
> Dive into the World of Parallel Programming The Go Parallel Website,
> sponsored
> by Intel and developed in partnership with Slashdot Media, is your hub for
> all
> things parallel software development, from weekly thought leadership blogs
> to
> news, videos, case studies, tutorials and more. Take a look and join the
> conversation now. http://goparallel.sourceforge.net/
> _______________________________________________
> Libmesh-users mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
diff --git
a/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
b/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
index 404218f..50e787b 100644
---
a/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
+++
b/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
@@ -62,6 +62,10 @@
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
+#include "libmesh/petsc_matrix.h"
+#include "libmesh/petsc_vector.h"
+#include "libmesh/petsc_macro.h"
+
#define x_scaling 1.3
// boundary IDs
@@ -522,6 +526,8 @@ int main (int argc, char** argv)
unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
+ NumericVector<Number>& nodal_coords = system.add_vector("nodal_coords");
+
LinearElasticity le(equation_systems);
system.attach_assemble_object(le);
@@ -565,6 +571,53 @@ int main (int argc, char** argv)
// Print information about the system to the screen.
equation_systems.print_info();
+ // If we're using PETSc, then define a near-nullspace, which is required
+ // in order for us to get good performance with a multigrid solver.
+ {
+ MeshBase::node_iterator node_it = system.get_mesh().nodes_begin();
+ const MeshBase::node_iterator node_end = system.get_mesh().nodes_end();
+
+ for( ; node_it != node_end; node_it++)
+ {
+ Node* node = *node_it;
+
+ libmesh_assert( node->n_dofs(system.number()) > 0 );
+
+ unsigned int dof_index_u = node->dof_number(system.number(), 0, 0);
+ unsigned int dof_index_v = node->dof_number(system.number(), 1, 0);
+ unsigned int dof_index_w = node->dof_number(system.number(), 2, 0);
+
+ if(node->processor_id() == system.comm().rank())
+ {
+ nodal_coords.set(dof_index_u, (*node)(0));
+ nodal_coords.set(dof_index_v, (*node)(1));
+ nodal_coords.set(dof_index_w, (*node)(2));
+ }
+
+ }
+ nodal_coords.close();
+
+ MatNullSpace matnull = PETSC_NULL;
+
+ PetscVector<Number>& petsc_nodal_coords =
+ cast_ref< PetscVector<Number>& >(nodal_coords);
+ PetscMatrix<Number>& petsc_system_matrix =
+ cast_ref< PetscMatrix<Number>& >(*system.matrix);
+
+ // Set the block size based on the spatial dimension of the problem.
+ // This info is used by MatNullSpaceCreateRigidBody.
+ VecSetBlockSize(petsc_nodal_coords.vec(), dim);
+
+ MatNullSpaceCreateRigidBody(petsc_nodal_coords.vec(), &matnull);
+
+ PetscErrorCode ierr = 0;
+ ierr = MatSetNearNullSpace(petsc_system_matrix.mat(), matnull);
+ CHKERRABORT(init.comm().get(),ierr);
+
+ ierr = MatNullSpaceDestroy(&matnull);
+ CHKERRABORT(init.comm().get(), ierr);
+ }
+
// Solve the system
system.solve();
------------------------------------------------------------------------------
Dive into the World of Parallel Programming The Go Parallel Website, sponsored
by Intel and developed in partnership with Slashdot Media, is your hub for all
things parallel software development, from weekly thought leadership blogs to
news, videos, case studies, tutorials and more. Take a look and join the
conversation now. http://goparallel.sourceforge.net/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
