Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Scipy LinearOperator with FreeDofs

More
1 year 10 months ago #4399 by pumpkinpieter
Hi there,

I'm working on optimizing an inverse of a large block sparse matrix : [ [A C], [B, D] ]   I've found that I can use the smaller blocks, and now I just need to be able to do an inverse of the matrix W = A - C D-1 B.  The matrix D is a mass matrix so I am able to invert it and use the freedofs from the space to give to the inverse.   The other matrices are also derived from Bilinear forms, in particular A is derived from the curl curl form plus a weighted mass form.  

To provide a method for inverting this matrix, I was thinking of using the Matrix Free method outlined in the section on interfacing with Numpy and Scipy in the documentation.  It requires only a method for applying the matrix W, which is no problem.  

The problem though is that I don't know how to get this process to recognize the freedofs.  In the documentation the example that is given is this:
Code:
import scipy import scipy.sparse.linalg tmp1 = f.vec.CreateVector() tmp2 = f.vec.CreateVector() def matvec(v): tmp1.FV().NumPy()[:] = v tmp2.data = a.mat * tmp1 return tmp2.FV().NumPy() A = scipy.sparse.linalg.LinearOperator( (a.mat.height,a.mat.width), matvec) u.vec.FV().NumPy()[:], succ = scipy.sparse.linalg.cg(A, f.vec.FV().NumPy())

It mentions in particular that the matrix a.mat is derived without Dirichlet dofs.

My question then is how can I incorporated Dirichlet dofs into the above example, and hence into my own process?  Is there a way to do this that someone knows?

Thanks,
Pieter Vandenberge
Portland State
 
More
1 year 10 months ago #4400 by joachim
Typically the freedofs enter into the preconditioner,
The range of the preconditioner defines the sub-space satisfying the (homogeneous) Dirichlet constraints.

The simplest possibility is to use a projector:

pre = Projector(mask=fes.FreeDofs(), range=True)

best,
Joachim
 
Time to create page: 0.180 seconds