### Form Issues - DG and Linear Solver

361
views
0
8 months ago by
I get the following error:
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-101-13f384f73a12> in <module>()
1 u = Function(V)
2 ffc_options = {"optimize": True, "quadrature_degree": 8}
----> 3 problem = LinearVariationalProblem(lhs_F,rhs_F,u,[],form_compiler_parameters=ffc_options)
4 solver  = LinearVariationalSolver(problem)
5 solver.solve()

/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc in __init__(self, a, L, u, bcs, form_compiler_parameters)
83
84         # Initialize C++ base class
---> 85         cpp.LinearVariationalProblem.__init__(self, a, L, u, bcs)
86
87

/usr/lib/python2.7/dist-packages/dolfin/cpp/fem.pyc in __init__(self, a, L, u, bcs)
2097
2098         """
-> 2099         _fem.LinearVariationalProblem_swiginit(self, _fem.new_LinearVariationalProblem(a, L, u, bcs))
2100
2101     def bilinear_form(self):

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
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

​

when I run the following code:

import numpy as np
from fenics import *

parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['optimize'] = True

L = 1.0
start = 0.0
h = 0.5
num_elements = int(np.ceil((L-start)/h))
mesh = IntervalMesh(num_elements,start,L)
V = FunctionSpace(mesh,"DG",1)

v = TestFunction(V)
u = TrialFunction(V)

u_left = 1.0
u_right = 0.0
bcs = Expression("(u_left-u_right)/(L-start)*(x[0]-start)+u_right",degree=1,
u_left=u_left, u_right=u_right, start=start, L=L)

K = Constant(1.0)
u_true = Expression("(1-x[0])*exp(-x[0]*x[0])", degree=1)
eps = 1.0 # Penalizes jumps in the solution.
# Penalizes jumps in the solution and the test function.
sig0 = Constant(1.0)
J0 = sig0/h*jump(v)*jump(u)*dS
# Penalizes jumps in the solution's and the test function's derivative.
sig1 = Constant(0.0)
J1 = sig1/h*jump(v.dx(0))*jump(u.dx(0))*dS
# The final weak form
f = Expression("exp(-x[0]*x[0])*(-4*x[0]*x[0]*x[0] + 4*x[0]*x[0]+ 6*x[0] - 2)", degree=1)
lhs_F = K*u.dx(0)*v.dx(0)*dx - avg(K*u.dx(0))*jump(v)*dS + eps*avg(K*v.dx(0))*dS + J0 + J1
rhs_F = f*v*dx - eps*K*bcs*v.dx(0)*ds

u = Function(V)
ffc_options = {"optimize": True, "quadrature_degree": 8}
problem = LinearVariationalProblem(lhs_F,rhs_F,u,[],form_compiler_parameters=ffc_options)
solver  = LinearVariationalSolver(problem)
solver.solve()
plot(u,interactive=True)​

Any thoughts why?

Community: FEniCS Project
1
Dear Saro,
Just a small question regarding your form.
What does it mean u.dx(0) as you wrote on the first term of the lhs?
I do not have plenty of experience, but I would probably check for the rank of J0, J1 and each term of the lhs. One of them could be missing the trial function (maybe "eps*avg(K*v.dx(0))*dS" ?).
Sorry if this do not make sense.
All the best, Murilo.
written 8 months ago by Murilo Moreira
1
I think your comment makes sense. The term eps*avg(K*v.dx(0))*dS is missing a dependence on u. In this form it belongs on the right hand side.
written 8 months ago by klunkean

3
8 months ago by
@Murilo Moreira and klunkean: you both are right. There should have been a dependence on u in that term. The program now runs.