[petsc-users] ML and linear elasticity

[Mark F. Adams]

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>