Hi all,
I wanted to see if I could get PETSc's GAMG functionality to work well with
systems_of_equations_ex6. There is a thread on the PETSc list here
<https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2013-August/018486.html>
about using GAMG for a "vanilla" linear elasticity problem, which is
exactly what systems_of_equations_ex6 provides. There was also some
discussion of GAMG on the libMesh list here
<http://sourceforge.net/p/libmesh/mailman/message/30819290/>.
I gather that the main requirement is to specify rigid-body modes as the
near-nullspace. 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?
Thanks,
David
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..aa766f4 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("near_nullspace");
+
LinearElasticity le(equation_systems);
system.attach_assemble_object(le);
@@ -565,6 +571,52 @@ 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();
+
+ DenseVector<Number> uvw;
+ DenseVector<Number> rotated_uvw;
+
+ 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);
+
+ 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
