- Thank you received: 3
Sparse inverse issue in 1D
2 years 11 months ago #4137
by ddrake
Sparse inverse issue in 1D was created 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])
Thanks for taking a look!
Best,
Dow
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)
Best,
Dow
Attachments:
2 years 11 months ago #4138
by joachim
Replied by joachim on topic Sparse inverse issue in 1D
Hi Dow,
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):
best,
Joachim
Code:
Inverse(freedofs=...)
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.096 seconds