Jed,
I've added the appended code to get the null space from the Mat. I did not see:
ierr = MatNullSpaceDestroy(&matnull);CHKERRQ(ierr);
in your ML code ... shouldn't that be here?
Mark
PetscFunctionBegin;
ierr = MatGetNearNullSpace( a_A, &mnull ); CHKERRQ(ierr);
if( !mnull ) {
ierr = PCSetCoordinates_AGG( pc, -1, PETSC_NULL ); CHKERRQ(ierr);
}
else {
PetscReal *nullvec;
PetscBool has_const;
PetscInt i,j,mlocal,nvec,bs;
const Vec *vecs; const PetscScalar *v;
ierr = MatGetLocalSize(a_A,&mlocal,PETSC_NULL);CHKERRQ(ierr);
ierr = MatNullSpaceGetVecs(mnull,&has_const,&nvec,&vecs);CHKERRQ(ierr);
ierr = MatGetBlockSize( a_A, &bs ); CHKERRQ( ierr );
ierr = PetscMalloc((nvec+!!has_const)*mlocal*sizeof
*nullvec,&nullvec);CHKERRQ(ierr);
if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
for (i=0; i<nvec; i++) {
ierr = VecGetArrayRead(vecs[i],&v);CHKERRQ(ierr);
for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] =
PetscRealPart(v[j]);
ierr = VecRestoreArrayRead(vecs[i],&v);CHKERRQ(ierr);
}
pc_gamg->data = nullvec;
pc_gamg->data_cell_cols = (nvec+!!has_const);
pc_gamg->data_cell_rows = bs;
}
PetscFunctionReturn(0);
}
On Apr 3, 2012, at 9:21 AM, Jed Brown wrote:
> On Tue, Apr 3, 2012 at 01:02, Karin&NiKo <niko.karin at gmail.com> wrote:
> I am not sure if the definition of the near null space of the operator
> suffice in order to get an optimal preconditioner for linear elasticity.
> Considering ML, one must also set the "num_PDEs" attribute, mustn't he? But
> it seems to me that this attribute cannot be set in the PETSc interface.
> Am I wrong?
>
> It is set using the "block size" of the matrix (MatSetBlockSize()). But that
> is only part of the game. For robustness, you should also provide null space
> information. Here is an example of using a coordinate Vec to create a
> MatNullSpace defining the rigid body modes and passing it along (from
> src/ksp/ksp/examples/tutorials/ex49.c).
>
> ierr = DMDAGetCoordinates(elas_da,&vel_coords);CHKERRQ(ierr);
> ierr = MatNullSpaceCreateRigidBody(vel_coords,&matnull);CHKERRQ(ierr);
> ierr = MatSetNearNullSpace(A,matnull);CHKERRQ(ierr);
> ierr = MatNullSpaceDestroy(&matnull);CHKERRQ(ierr);
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120403/6ff7e507/attachment.htm>
