Getting error to solve Darcy equation

1

Hi,

I am trying to solve Darcy equation but I am getting this error. Can anyone help me to understand what the problem is?

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

from dolfin import *
from fenics import *
k1 = 'exp(x[1])'
k = Expression(k1, degree = 2)
mu = 0.001

inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

# Create mesh and define function spaces
domain = Rectangle(Point(0, 0), Point(1.0, 1.0))	
mesh = generate_mesh(domain, 50)
basis_order = 2 

Element1 =  VectorElement("CG", mesh.ufl_cell(), basis_order)
Element2 =  FiniteElement("CG", mesh.ufl_cell(), 1)
W_elem = MixedElement([Element1, Element2])
W = FunctionSpace(mesh, W_elem)
V = FunctionSpace(mesh, 'P', 1)  #1st order FEM

(v,q) = TestFunctions(W)
w = Function(W)
(u, p) = split(w)

bcu_walls = DirichletBC(W.sub(0).sub(1),Constant(0.0), walls) #free slip(partial(u)/partial(y) =0 $ v =0)
bcp_outflow = DirichletBC(W.sub(1), Constant(0.0), outflow) # outlet pressure
bcp_inflow = DirichletBC(W.sub(1), Constant(1.0), inflow) # inlet pressure
bc = [bcu_walls,bcp_inflow, bcp_outflow]

f = Constant(0.0)

F1 = dot(v, (k/mu)*u)*dx - div(v)*p*dx + q*div(u)*dx 

F2 = q*f*dx 

print ("Solve primal problem")
problem_primal = solve ( F1 == F2, w, bc)

(u_h, p_h) = problem_primal.solve().split()
file_pvd = File("velocity.pvd")
ufile_pvd << u_h
pfile_pvd = File("pressure.pvd")
pfile_pvd << p_h
1 Like
2

As you are solving a linear variational problem,
w should be a TrialFunction, i.e.

(u, p) = TrialFunctions(W)
4

dokken,
Thanks for the reply.

Right! it solved the error. But I am getting the wrong results. Any idea where I am making mistake??

5

I want to solve this problem:

This program solves steady-state Darcy's flow on a
square plate with spatially varying permeability.
        (K/mu)*u + grad(p) = 0
                  div(u) = 0
Which, in weak form reads:
 (v, (K/mu)*u) - (div(v), p) = 0           for all v
               (q, div(u)) = 0             for all q
7

Any comments on this???