Solving PDE with a parameter

Question

I am new in FEniCS and i am trying to solve a 2D time-independant Schrödinger with a parameter, the idea is to study the evolution in the energy (eigenvalues) with the increase of the parameter (in my case alpha). So, i made a loop before the definition of the variational form, created a file for print the eigenvalues and i thought that the PDE will be solve n-times, but after the first case, it stoped.

domain = Circle(Point(0,0),5.5)+Circle(Point(10,0),5.5)+Circle(Point(5,8.66),5.5)


mesh = generate_mesh(domain,8)

n = mesh.num_vertices()
d = mesh.geometry().dim()

# Create the triangulation
mesh_coordinates = mesh.coordinates().reshape((n, d))
triangles = np.asarray([cell.entities(0) for cell in cells(mesh)])
triangulation = tri.Triangulation(mesh_coordinates[:, 0],
                                  mesh_coordinates[:, 1],
                                  triangles)

# Plot the mesh
plt.figure()
plt.triplot(triangulation)
plt.savefig('mesh.png')


V = FunctionSpace(mesh, 'Lagrange', 1)

#build essential boundary conditions
def u0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V,Constant(0.0) , u0_boundary)

#define functions
u = TrialFunction(V)
v = TestFunction(V)

Pot = Expression('0.0')


alpha = 0.0

for alpha in [float(j) / 100 for j in range(0, 100, 1)]:
    print alpha
    terminocampoe = Expression('1*x[0]')
    #define problem
    a = (inner(grad(u), grad(v)) \
    + Pot*u*v + alpha*inner(terminocampoe,v))*dx
    m = u*v*dx
#
    A = PETScMatrix()
    M = PETScMatrix()
    _ = PETScVector()
    L = Constant(1)*v*dx

    assemble_system(a, L, bc, A_tensor=A, b_tensor=_)
    assemble_system(m, L,bc, A_tensor=M, b_tensor=_)
    assemble_system(m, L, A_tensor=M, b_tensor=_)
    #create eigensolver
    eigensolver = SLEPcEigenSolver(A,M)
    eigensolver.parameters['spectrum'] = 'smallest magnitude'
    eigensolver.parameters['tolerance'] = 1.e-15

    #solve for eigenvalues
    eigensolver.solve()
    p = eigensolver.get_number_converged()
    print "numero de eigenvalores", p
    u = Function(V)

    energiaArchivo = open('energia_campoe','w')
    energiaArchivo.write(str(alpha)+' ')
    energiaArchivo.close()
    for i in range(0,20):
        energiaArchivo = open('energia_campoe','a')
        #extract next eigenpair
        r, c, rx, cx = eigensolver.get_eigenpair(i)
        energiaArchivo.write(str(r)+' ')

and the error is:

Traceback (most recent call last):
  File "prueba_tresminf_campoe.py", line 72, in <module>
    assemble_system(a, L, bc, A_tensor=A, b_tensor=_)
  File "/home/adriana/.hashdist/bld/profile/mcfwu2rmfcy2/lib/python2.7/site-packages/dolfin/fem/assembling.py", line 254, in assemble_system
    assembler = cpp.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)
  File "/home/adriana/.hashdist/bld/profile/mcfwu2rmfcy2/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 2338, in __init__
    _fem.SystemAssembler_swiginit(self,_fem.new_SystemAssembler(*args))
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to assemble system.                                                                                                                                                      
*** Reason:  expected a bilinear form for a.                                                                                                                                                 
*** Where:   This error was encountered inside SystemAssembler.cpp.                                                                                                                          
*** Process: 0                                                                                                                                                                               
***                                                                                                                                                                                          
*** DOLFIN version: 1.6.0                                                                                                                                                                    
*** Git changeset:                                                                                                                                                                           
*** -------------------------------------------------------------------------

I don't understand what happen. Please help me.

Answer 1

Hi, if you examine $a$ you will see that alpha*inner(terminocampoe,v)*dx is not bilinear. Further, even if $a$ had all its terms defined in terms of TrialFunction and Testfunction the code will only work in the first iteration of the loop, for by u = Function(V) the $u$ seizes to be a TrialFunction and thus $a$ is no longer bilinear form. Note also that the form could be defined outside the loop and inside you would only change the value of alpha, alpha.assign(new_alpha_value) where $\alpha$ is defined as Constant. Finally, as it is now you might be extracting more eigenvalues then converged, consider range(min(20, p)).

Answer 2

Hi,

the fundamental problem is that you overwrite the trial function u, by the function u, since you specify somewhere the following line:

u = Function(V)

As an immediate result, there is no bilinear form anymore to assemble in the next cycle and your computation crashes.

Please consider the following simple poisson example in order to see what's going on and to learn how you can include a varying parameter in a simple and straightforward manner:

from dolfin import * 

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define trial and test function
u = TrialFunction(V)
v = TestFunction(V)

# The following loop reproduces the error
#for kappa in [0.8, 1.0, 1.2]:
    #print 'Kappa', kappa 
    #f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
    #g = Expression("sin(5*x[0])", degree=2)
    #a = inner(kappa*grad(u), grad(v))*dx
    #L = f*v*dx + g*v*ds

    ## Compute solution
    #u = Function(V)
    #solve(a == L, u, bc)

# This is a working solution, but not really smart either:
#for kappa in [0.8, 1.0, 1.2]:
    #print 'Kappa', kappa 
    #f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
    #g = Expression("sin(5*x[0])", degree=2)
    #a = inner(kappa*grad(u), grad(v))*dx
    #L = f*v*dx + g*v*ds

    ## Compute solution
    #uh = Function(V)
    #solve(a == L, uh, bc)

# This is THE way to test a varying parameter
uh = Function(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2) 
kappa = Constant(0.)
a = inner(kappa*grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

for kappav in [0.8, 1.0, 1.2]:
    print 'kappa', kappav
    # The following line does the trick: 
    kappa.assign(kappav) 
    # Compute solution
    solve(a == L, uh, bc)