Hi, everyone
I am trying to resolve the Richard’s equation in steady-state(no time dependence) and equation takes the form:
which is a similar problem to the Poisson nonlinear problem explained in the FEniCS Tutorial, the difference here is that the nonlinear coefficient do not depend only on our target function ‘u’, it also depends on x that is a point in Omega=[0,1]x[0,1].
I tried to change the nonlinear coefficient in ft05_poisson_nonlinear.py example and I do not know why it does not work, here is what I have tried so far:
""" Fenics Code used to search the solution of Richard's Equation in a steady state (no time dependence)."""
# Packages used
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
# Defining Relativity hydraulic conductivity
# Defining k_s(x), positive function.
def k_s(x,y):
return x
# Haverkamp model
# Constants used
A=1
B=1
Gamma=1
def k(x,y,u):
return k_s(x,y)*(A/(A+((abs(u)/B)**Gamma)))
# Use SymPy to compute f from the manufactured solution u
import sympy as sym
x, y = sym.symbols('x[0], x[1]')
u = 0
f = 1
# f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)
# Create mesh and define function space
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression(u_code, degree=1)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression(f_code, degree=1)
F = k(x,y,u)*dot(grad(u), grad(v))*dx - f*v*dx
# Compute solution
solve(F==0,u,bc)
# Plot solution
plot(u, title='function with Haverkamp Model')
plt.show()
# Save solution to file in VTK format
vtkfile = File('poisson/final_work1.pvd')
vtkfile << u
When I set the k_s(x,y) function as constant the program works but in other case it shows me the next error:
Traceback (most recent call last):
File "final_work.py", line 50, in <module>
F = k(x,y,u)*dot(grad(u), grad(v))*dx - f*v*dx
File "final_work.py", line 21, in k
return k_s(x,y)*(A/(A+((abs(u)/B)**Gamma)))
TypeError: unsupported operand type(s) for *: 'Symbol' and 'Division'
I am working with Ubuntu 18.14.1 LTS
Thank you in advance for helping me to solve this problem.