Hi Everyone,
I am stuck trying to verify my code for implementing Darcy’s flow though a unit square.
The 2 equations for this flow are: u_bar = -K * grad p and div (u_bar) = 0;
Given: K = max(0.10, exp(-pow(10y-1.0sin(10*x)-5.0, 2)) * I. Where I is an identity matrix, u_bar represents velocity, p represents pressure and the flow uses pressure boundary conditions with p = 1 at x = 0 and p = 0 at x = 1.
Now using the method of manufactured solutions, I selected the analytical expression for pressure to be p = 1 - x^2.
I have tried to solve it on FEniCS with no success. Please see the code below and let me know where I got it wrong and what I should add to make it run successfully:
from fenics import *
import numpy as np
from ufl import nabla_div
Creating mesh and defining function space
mesh = UnitSquareMesh(24, 24)
V = FunctionSpace(mesh, ‘P’, 1)
Defining Exact Solution for Pressure Distribution
p_e = Expression(‘1 - x[0]*x[0]’, degree=2)
Defining Dirichlet boundary
p_L = Constant(1.0)
def boundary_L(x, on_boundary):
return on_boundary and near(x[0], 0)
bc_L = DirichletBC(V, p_L, boundary_L)
p_R = Constant(0.0)
def boundary_R(x, on_boundary):
return on_boundary and near(x[0], 1)
bc_R = DirichletBC(V, p_R, boundary_R)
bcs = [bc_L, bc_R]
Defining variational problem
p = TrialFunction(V)
v = TestFunction(V)
d = 2
I = Identity(d)
M = Expression(‘fmax(0.10, exp(-pow(10x[1]-1.0sin(10*x[0])-5.0, 2)))’, degree=2)
K = M * I
a = dot(K * grad p , grad(v))dx
f1 = Expression(’(-2x[0], 0)’, degree=1)
f = nabla_div(-K * interpolate(f1, V))
L = inner(f, v)*dx
Computing solutions
p = Function(V)
solve(a == L, p, bcs)
p_e_f = interpolate(p_e, V)
diff_p = Function(V)
diff_p.vector()[:] = p.vector() - p_e_f.vector()
Saving Solution in VTK format
file = File(‘Pressure_Gradient.pvd’)
file << p
File(‘diff_p.pvd’) << diff_p
The error message is:
Shapes do not match: and .
Traceback (most recent call last):
File “Numerical_Exact.py”, line 43, in
L = inner(f, v)*dx
File “/usr/lib/python3/dist-packages/ufl/operators.py”, line 144, in inner
return Inner(a, b)
File “/usr/lib/python3/dist-packages/ufl/tensoralgebra.py”, line 161, in new
error(“Shapes do not match: %s and %s.” % (ufl_err_str(a), ufl_err_str(b)))
File “/usr/lib/python3/dist-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: NablaDiv id=140429426823800> and <Argument id=140429427152376.