Bilinear form into sparse matrix

Rauf · November 16, 2023, 2:41pm

Hi,
I hope all are doing well. I want to compute eigenvalue using Scipy solver. For that, I converted the bilinear form to a sparse matrix. However, I did not get correct values compared to other eigenvalue solvers (Arnoldi) and analytical values. You can see in following attachment. The code for that is:

from ngsolve import *
import numpy as np
from ngsolve.meshes import Make1DMesh
from ngsolve.webgui import Draw
import scipy.sparse as sp

L = 10.0
thick = 0.03
width = 0.01
E = 70e3
nu = 0
N0 = 1e-3

EI = Ewidththick**3/12
GS = E/2/(1+nu)thickwidth
kappa = 5.0/6.0

mapping = lambda x : (x*L)
mesh = Make1DMesh(100, mapping=mapping, periodic=False)

fes = H1(mesh, order=2, dirichlet=“left|right”)H1(mesh, order=2, dirichlet=“left”)
(w,beta),(dw,dbeta) = fes.TnT()
a = BilinearForm(fes, symmetric=False)
a += EI(grad(beta)grad(dbeta))dx + kappaGS((grad(w)-beta)*(grad(dw)-dbeta))*dx

m = BilinearForm(fes,symmetric=False)
m += N0*grad(w)*grad(dw)*dx

a.Assemble()
m.Assemble()

u = GridFunction(fes,multidim=3)

Arnoldi method

with TaskManager():
lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),list(u.vecs), shift=1)

For sparse matrix

rows,cols,vals = a.mat.COO()
A_csr = sp.csr_matrix((vals,(rows,cols)))

rows,cols,vals = m.mat.COO()
M_csr = sp.csr_matrix((vals,(rows,cols)))

Scipy Eigenvalue solver

Timoshenko_Beam.ipynb (4.1 KB)

n_eig = 4
from scipy.sparse.linalg import eigs
r_, w_ = eigs(A=A_csr, k=n_eig, M=M_csr, sigma=0.0)

christopher · November 18, 2023, 2:50pm

Hi,
you don’t handle the dirichlet dofs when parsing to scipy, have a look here how this can be done:

github.com

NGSolve/ngsolve/blob/master/python/directsolvers.py

from ngsolve import BaseMatrix, BitArray, BilinearForm, BaseVector

class SuperLU(BaseMatrix):
    # def __init__(self, a: BaseMatrix | BilinearForm, freedofs: BitArray = None):
    def __init__(self, a, freedofs: BitArray = None):        
        super().__init__()
        self.a = a
        self.freedofs = freedofs

    def Update(self):
        import scipy.sparse as sp
        import scipy.sparse.linalg as spla
        a = self.a if isinstance(self.a, BaseMatrix) else self.a.mat
        mat = sp.csr_matrix(a.CSR())
        if self.freedofs is not None:
            self.fd = list(self.freedofs)
            mat = mat[self.fd,:][:,self.fd]
        self.lu = spla.factorized(sp.csc_matrix(mat))

    def Mult(self, x: BaseVector, y: BaseVector):

This file has been truncated. show original

best, Christopher

Rauf · November 22, 2023, 1:50pm

Hi Christopher,
Thanks for your reply. Yes, you are right the Dirichlet dofs can not include when parsing to scipy. I tried by this given algorithm but I faced some error. I am trying to call this class but I faced the following error:

TypeError Traceback (most recent call last)
Cell In[7], line 1
----> 1 A_sp = SuperLU(a=ab, freedofs = fes.FreeDofs(), BitArray = None)

Cell In[6], line 4, in SuperLU.init(self, a, freedofs, BitArray)
3 def init(self, a, freedofs, BitArray = None):
----> 4 super().init()
5 self.a = a
6 self.freedofs = freedofs

TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. ngsolve.comp.BilinearForm(space: ngsolve.comp.FESpace, name: str = ‘biform_from_py’, **kwargs)
2. ngsolve.comp.BilinearForm(trialspace: ngsolve.comp.FESpace, testspace: ngsolve.comp.FESpace, name: str = ‘biform_from_py’, **kwargs)
3. ngsolve.comp.BilinearForm(arg0: ngsolve.comp.SumOfIntegrals, **kwargs)

Invoked with:

I am also trying by import the SuperLU from ngsolve directsolver but I faced again error which is given as:

No module named ‘ngsolve.directsolvers’

You can see in the following attachment. Can you tell me how can I add dirichlet dofs? It will be very helpful for me.

Thanks,
Best,
Rauf.
Timoshenko_Sparse.ipynb (6.1 KB)

christopher · November 24, 2023, 12:51pm

no i did not mean you should use this class, but look how dirichlet handling to scipy is done there

Rauf · November 28, 2023, 5:22pm

Hi Christopher,
Thanks for your reply. I successfully handled the Dirichlet to Scipy. I got the same results as other methods (i.e. Arnoldi) and analytical value. Now, I am trying to plot the eigenvector but I have an issue with the vector dimensions. The dimension of the eigenvector did not match the dimension of the grid function. Do you have any idea how can I plot the eigenvector? The file is attached below.
Timoshenko_Scipy.ipynb (4.1 KB)

Thanks,

Best,
Rauf

christopher · November 28, 2023, 8:15pm

I think you need to also use the freedofs array for assigning like in the superlu implementation:

            y.FV().NumPy()[self.fd] = self.lu(x.FV().NumPy()[self.fd])