dolfin team mailing list archive

[Raphael Kruse <question141904@xxxxxxxxxxxxxxxxxxxxx>]

Question #141904 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/141904

    Status: Answered => Open

Raphael Kruse is still having a problem:
Hi,

I wrote the following code to solve the problem. I would be very
grateful if someone with experience could go over it. Very likely it is
possible to make the code more efficient.


"""
This scripts builds a FEM-approximation of a Carleman operator Q, which is
given by an integral kernel q(x,y). Then it solves the corresponding
eigenproblem of finding all eigenpairs (lambda,u) such that Qu = lambda*u.

Q is given by [Qu](y) = \int q(x,y) u(x) dx

Author: Raphael Kruse (Univerity Bielefeld), 26 Jan 2011
"""

from dolfin import *
import numpy as np
from numpy import linalg as LA

# parameters for the kernel q
sigma = 1       # standard deviation
gamma = 0.2     # correlation length

dof = 200    # degrees of freedom
tol = 1E-14  # tolerance for eigenvalues

# Defining the integral kernel
q_1d = Expression('sigma*sigma*exp(-pow(x[0]-y,2)/(gamma*gamma))')
q_1d.sigma = sigma
q_1d.gamma = gamma

# Test for PETSc and SLEPc
if not has_la_backend("PETSc"):
    print "DOLFIN has not been configured with PETSc. Exiting."
    exit()

if not has_slepc():
    print "DOLFIN has not been configured with SLEPc. Exiting."
    exit()

# Define mesh, function space
mesh = UnitInterval(dof-1)
coor = mesh.coordinates()

V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)

# building the mass matrix M as PETScMatrix

print "Building the mass matrix M...",
M = PETScMatrix()
a = u*v*dx
assemble(a, tensor=M)

print "Done"

# Building the Qu matrix for u running through the testfunctions
print "Building the Q matrix...",
Qu = np.zeros( (dof,dof) )  # holds all values of convolutions of q with
                            # testfunctions
for i in range(dof):
    q_1d.y = coor[i][0]     # y runs through the coordinates of the mesh
    L = q_1d*v*dx
    b = assemble(L, mesh=mesh)
    Qu[i,:] = b.array()

Qu = Qu.transpose()

# Computing the values of \int Qu(y) v(y) dy for each u,v

Q = np.zeros( (dof, dof) ) # the matrix Q
Fct_Qu = Function(V)    
for k in range(dof):
    quk = np.array(Qu[k][:])  # runs through the values of u
    Fct_Qu.vector()[:] = quk
    L = Fct_Qu*v*dx
    b = assemble(L, mesh=mesh)
    Q[:,k] = b.array()

print "Done"

# transforming Q into PETScMatrix
print "Transforming Q into PETScMatrix format...",
Qh = PETScMatrix(M)

for i in range(dof):
    for k in range(dof):
        Qh[i,k] = Q[i,k]

print "Done"
#print Qh.array()

# Create eigensolver
eigensolver = SLEPcEigenSolver()
eigensolver.parameters["solver"] = "lapack"

# Compute all eigenvalues of Q x = \lambda M x

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.