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.

Sparse inverse issue in 1D

More
2 years 2 months ago #4137 by ddrake
Hi, I'm having trouble inverting the sparse matrix using the script below. For any order p, I think the number of freedofs of the P space should equal the number of dofs of the A space, namely N*(p+1), where N is the number of 1D elements, resulting in a square matrix which can be shown to be invertible.

In the case p=0, both umfpack and pardiso give a warning 'sizes don't match', but the inverse is computed anyway yielding the correct solution. When p>0, both umfpack and pardiso throw an exception 'Symbolic factorization failed' when the inverted matrix is multiplied by a vector of length N*(p+1). It behaves as if the non-freedof of the P space due to the Dirichlet boundary condition is not being accounted for. This issue does not occur if I use a product space, e.g. AP = FESpace([A, P])

Code:
from ngsolve import x, dx, grad import ngsolve as ng from ngsolve.meshes import Make1DMesh N = 10 p = 1 mesh = Make1DMesh(N) A = ng.L2(mesh, order=p, all_dofs_together=True) P = ng.H1(mesh, order=p+1, dirichlet='right') theta = A.TestFunction() phi = P.TrialFunction() a2p = ng.BilinearForm(trialspace=P, testspace=A) a2p += grad(phi) * theta * dx a2p.Assemble() a2pinv = a2p.mat.Inverse(freedofs=P.FreeDofs(), inverse='umfpack') m = ng.BilinearForm(A) u, v = A.TnT() m += u * v * dx m.Assemble() gfalpha = ng.GridFunction(A) gfphi = ng.GridFunction(P) gfalpha.Set(x) tvec = gfalpha.vec.CreateVector() tvec.data = m.mat * gfalpha.vec gfphi.vec.data = a2pinv * tvec print("solution:", gfphi.vec)
Thanks for taking a look!

Best,
Dow
Attachments:
More
2 years 2 months ago #4138 by joachim
Replied by joachim on topic Sparse inverse issue in 1D
Hi Dow,
Code:
Inverse(freedofs=...)
does not distinguish between rows and cols, so this does not work.

You can eliminate the right-most dof completely from the system, by setting it as unused-dof, and then compress the space to kick-out the unused dofs completely (P.ndof changes):
You are not authorised to see this attachment.
You are not authorised to see this attachment.
Code:
P = ng.H1(mesh, order=p+1, dirichlet='right') P.couplingtype[N] = COUPLING_TYPE.UNUSED_DOF P = ng.Compress(P) print (P.ndof, A.ndof)

best,
Joachim
Attachments:
The following user(s) said Thank You: ddrake
Time to create page: 0.116 seconds